# Theory of photonic crystal polaritons in periodically patterned multilayer
waveguides

Simone Zanotti  Dipartimento di Fisica, Università di Pavia, via Bassi 6,
I-27100 Pavia, Italy  Hai Son Nguyen  Univ Lyon, Ecole Centrale de Lyon, CNRS,
INSA Lyon, Université Claude Bernard Lyon 1, CPE Lyon, CNRS, INL, UMR5270,
Ecully 69130, France  Institut Universitaire de France (IUF), Paris, France
Momchil Minkov  Flexcompute Inc., 130 Trapelo Rd., Belmont, Massachusetts
02478, USA  Lucio C. Andreani  Dipartimento di Fisica, Università di Pavia,
via Bassi 6, I-27100 Pavia, Italy  Dario Gerace  Dipartimento di Fisica,
Università di Pavia, via Bassi 6, I-27100 Pavia, Italy

###### Abstract

We present a formalism for studying the radiation-matter interaction in
multilayered dielectric structures with active semiconductor quantum wells
patterned with an in-plane periodic lattice. The theory is based on the
diagonalization of the generalized Hopfield matrix, and it includes loss
channels in a non-Hermitian formulation. Hybrid elementary excitations named
photonic crystal polaritons arise in these systems, whose detailed dispersion
and loss characteristics are shown to depend on material composition as well
as on symmetry properties of the lattice. We show the generality of the
approach by calculating polariton dispersions in very diverse material
platforms, such as multilayered perovskite-based lattices or inorganic
semiconductor heterostructures. As an application of the method, we show how
to engineer lossless polariton modes through excitonic coupling to bound
states in the continuum at either zero or finite in-plane wavevector, and
discuss their topological properties. Detailed comparison with a semiclassical
approach based on the scattering matrix method is provided, which allows to
interpret the optical spectra in terms of polarization-dependent excitation of
the different polariton branches. This work introduces an efficient and
invaluably versatile numerical approach to engineer photonic crystal
polaritons, with potential applications ranging from low-threshold lasers to
symmetry-protected propagating modes of hybrid radiation-matter states.

##  I Introduction

Non-perturbative radiation-matter coupling between dipole-active material
excitations and eigenmodes of the electromagnetic field leads to the concept
of polaritons, i.e., hybrid excitations of mixed nature involving at least two
fields with different characteristics [1]. Microcavity exciton-polaritons, in
particular, arising from the strong coupling between bound electron-hole pairs
in semiconductor quantum wells (QW excitons) and planar microcavity photons,
have been shown to behave as low-mass elementary excitations with bosonic-like
statistics [2], thus allowing to unravel phenomena that are typically
displayed in weakly interacting bosonic gases [3], such as Bose-Einstein
condensation [4] and superfluidity [5], with potential applications in all-
optical devices such as low-threshold lasing [6], all-optical switching [7]
and routing [8, 9], and more recently quantum technologies [10]. The concept
of exciton-polaritons has then been extended to excitons coupled to
periodically patterned metasurfaces with sub-wavelength pitch, also named
photonic crystal polaritons [11]. These excitations have been shown to arise
in very diverse materials and platforms, ranging from inorganic semiconductor
heterostructures [12, 13, 14], to two-dimensional materials (TDM) deposited on
a patterned metasurface [15, 16, 17, 18] or layered perovskites within a
passive backbone [19, 20, 21]. As it will be shown in the following, in all
these different scenarios the physics of photonic crystal polaritons can
essentially be described within a common theoretical background.

While polaritons can be understood as semiclassical excitations of the
electromagnetic field in a polarizable medium [22], a full quantum-mechanical
theory of polaritons was originally proposed by Hopfield [23] and Agranovich
[24] in two seminal papers. In particular, the Hopfield formulation is based
on a transformation of operators and can be extended to include imaginary
parts of the oscillator energies, thus making the theory a non-hermitian one.
While a generalized formulation of the Hopfield method has been reported in
the literature [25, 26], a Hopfield-based quantum theory accounting for an
arbitrary multilayered structures with an in-plane pattern is still lacking.
In this work, we extend a previously formulated theory of radiation-matter
interaction in fully etched thin photonic crystal slabs containing a single
active layer [11] to describe an arbitrarily etched multilayered structure. In
particular, we consider a guided mode expansion (GME) basis for the
approximate photonic modes that are the closed system solution of a
multilayered planar waveguide with partially etched core [27, 28], in which
losses are calculated within a perturbative approach by coupling to the
continuum of radiative modes in the light cone [27] to obtain the eigenvalues
imaginary parts. Then we formulate a generalized Hopflied Hamiltonian
accounting for radiation-matter coupling to an arbitrary number of active
layers characterized by an excitonic response. We generalize previous
treatments by including the polarization response of the excitonic transition.

As applications of this non-Hermitian theory, we show that the scaling of the
polariton Rabi splitting deviates from the expected $\sqrt{K}$ dependence as a
function of the number of active layers ($K$), due to the guided mode profile
that weights the light-matter coupling on the waveguide thickness. Then, we
show that the polariton dispersion in square lattices of perovskite
metasurfaces can be obtained, and it compares favorably well with either
rigorous coupled wave analysis simulations. Finally, we further demonstrate
the generality and versatility of the approach by analysing the effect of
strong light-matter coupling between QW excitons and bound-states in the
continuum (BIC) supported in partially etched triangular lattice photonic
crystal slab. We perform a systematic analysis of the topological properties
of such polariton BICs, relating in particular the topological charge to the
polarization singularity in reciprocal space. These studies acquire particular
significance in light of recent experiments unravelling the topological nature
of photonic crystal polaritons [29, 20, 21, 14]. Finally, we clarify the role
of non-radiative versus radiative loss contributions in the BIC optical
response.

##  II Theory

The class of systems we consider in this work are schematically represented in
Fig. 1, where the active layers (i.e., QWs containing bound electron-hole
pairs) are coupled to a periodic photonic landscape. In the following, we
summarize the basic steps to formulate a Hopfield matrix description of this
problem.

Figure 1: (Colour online) (a) Photonic crystal slab with quantum wells
embedded from the centre of the slab in $z=0$ to the top. (b) Brillouin zone
with high-symmetry points for a square lattice.

###  II.1 Photonic modes: the Guided-Mode Expansion

We start with the formalism to numerically calculate the non-Hermitian
photonic modes by the guided-mode expansion (GME) method [27][28]. This method
is particularly efficient in calculating the (approximated) photonic
dispersion (real part of the eigenvalues) of layered structures, with an
arbitrary in-plane periodicity, as well as their imaginary part (related to
out-of-plane scattering losses) within perturbation theory. Here we give a
brief summary of the method for completeness, and to fix the notation for the
following.

The GME method assumes the presence of linear, isotropic, non-dispersive,
lossless, and non magnetic materials in absence of free charges. Under such
hypotheses, Maxwell equations can be written as an eigenvalue problem for the
harmonic magnetic field,
$\mathbf{H}(\mathbf{r},t)=\mathbf{H}(\mathbf{r})e^{-i\omega t}$:

$\hat{\Theta}\mathbf{H}(\mathbf{r})=\frac{\omega^{2}}{c^{2}}\mathbf{H}(\mathbf{r}),\quad\hat{\Theta}\doteq\nabla\times\frac{1}{\varepsilon(\mathbf{r})}\times,$
(1)

where $\varepsilon(\mathbf{r})$ is the space-dependent relative permittivity.
The electric field $\mathbf{E}(\mathbf{r},t)=\mathbf{E}(\mathbf{r})e^{-i\omega
t}$ can be recovered from the magnetic one:

$\mathbf{E}(\mathbf{r})=\frac{i}{\omega\varepsilon_{0}\varepsilon(\mathbf{r})}\nabla\times\mathbf{H}(\mathbf{r}).$
(2)

In the GME method, the magnetic field is expanded in a basis of orthonormal
modes:

$\mathbf{H}(\mathbf{r})=\sum_{\mu}c_{\mu}\mathbf{H}_{\mu}(\mathbf{r}),\quad\int_{\text{unit
cell}}\mathbf{H}_{\mu}(\mathbf{r})\cdot\mathbf{H}_{\nu}(\mathbf{r})d\mathbf{r}=\delta_{\mu,\nu},$
(3)

leading to the linear eigenvalue problem:

$\displaystyle\sum_{\nu}\mathcal{H}_{\mu\nu}c_{\nu}=\frac{\omega^{2}}{c^{2}}c_{\mu},$
$\displaystyle\mathcal{H}_{\mu,\nu}=\int\frac{1}{\varepsilon(\mathbf{r})}\left(\nabla\times\mathbf{H}_{\nu}^{*}(\mathbf{r})\right)\cdot\left(\nabla\times\mathbf{H}_{\mu}(\mathbf{r})\right)d\mathbf{r}.$
(4)

In order to evaluate the matrix elements in (II.1), the magnetic field is
expanded for each in-plane wavevector $\mathbf{k}$ in the suitable basis of
the guided modes:

$\mathbf{H}_{\mathbf{k}}(\mathbf{r})=\sum_{\mathbf{G},\alpha}c(\mathbf{k}+\mathbf{G},\alpha)\mathbf{H}^{\text{guided}}_{\mathbf{k}+\mathbf{G},\alpha}(\mathbf{r}),$
(5)

where $\mathbf{G}$ is a reciprocal lattice vector, $\alpha$ is index of the
guided mode which can be explicitly given as:

$\mathbf{H}^{\text{guided}}_{\mathbf{k}+\mathbf{G},\alpha}(\mathbf{r})=\mathbf{h}_{\mathbf{k}+\mathbf{G},\alpha}(z)e^{i(\mathbf{k}+\mathbf{G})\cdot\rho},$
(6)

$\rho$ being the in-plane position vector $(x,y)$, and
$\mathbf{h}_{\mathbf{k}+\mathbf{G},\alpha}(z)$ the profile of the guided mode.

This procedure leads to purely real frequencies of the photonic modes,
$\omega_{\mathbf{k}n}$. The radiative losses of leaky modes lying above the
light line are calculated with the photonic analogue of the Fermi golden rule,
yielding the (wave vector-dependent) imaginary part of each eigenvalue as:

$\Im\left(\frac{\omega^{2}_{\mathbf{k}}}{c^{2}}\right)=-\sum_{\mathbf{G},p}\>\sum_{o=l,u}|\mathcal{H}_{\mathbf{k}r}|^{2}\rho_{0}\left(\mathbf{k}+\mathbf{G},\omega_{1}\right),$
(7)

where $\rho_{0}$ is the 1D density of photonic states,
$\mathcal{H}_{\mathbf{k}r}$ is the coupling matrix element defined as:

$\mathcal{H}_{\mathbf{k}r}=\int\frac{1}{\varepsilon(\mathbf{r})}\left(\nabla\times\mathbf{H}^{*}_{\mathbf{k}}(\mathbf{r})\right)\cdot\left(\nabla\times\mathbf{H}^{\text{rad}}_{\mathbf{k}+\mathbf{G},p,o}(\mathbf{r})\right)d\mathbf{r},$
(8)

$p$ represents the polarisation of the radiative modes (TE,TM), while the
summation over $o$ takes into account the leakage through the upper ($u$) and
lower ($l$) cladding.

Starting from the second-quantised formalism we can write the photonic
Hamiltonian as:

$H_{\text{ph}}=\sum_{\mathbf{k},n}\hbar\omega_{\mathbf{k}n}\left(\hat{a}_{\mathbf{k}n}^{\dagger}\hat{a}_{\mathbf{k}n}+\frac{1}{2}\right),$
(9)

where $\omega_{\mathbf{k}n}$ is the (real) frequency of the n-th band at given
wave vector $\mathbf{k}$ calculated with GME, and $\hat{a}_{\mathbf{k}n}$ is
the annihilation operator which follows bosonic commutation rules:

$\left[\hat{a}_{\mathbf{k}n},\hat{a}^{\dagger}_{\mathbf{k}^{\prime}n^{\prime}}\right]=\delta_{\mathbf{k},\mathbf{k}^{\prime}}\delta_{n,n^{\prime}},\quad\left[\hat{a}_{\mathbf{k}n},\hat{a}_{\mathbf{k}^{\prime}n^{\prime}}\right]=\left[\hat{a}^{\dagger}_{\mathbf{k}n},\hat{a}^{\dagger}_{\mathbf{k}^{\prime}n^{\prime}}\right]=0$
(10)

When generalizing to the non-Hermitian Hopfield matrix, the photonic
eigenvalues will be taken as the complex eigenmodes
$\tilde{\omega}_{\mathbf{k}n}=\omega_{\mathbf{k}n}-i\gamma_{\mathbf{k}n}$ with
$\gamma_{\mathbf{k}n}=-\Im\left[\sqrt{\Re\left(\omega^{2}_{\mathbf{k}n}\right)+i\Im\left(\omega^{2}_{\mathbf{k}n}\right)}\right]>0$,
in which $\Im\left(\omega^{2}_{\mathbf{k}n}\right)<0$ is obtained from Eq.
(7), see Appendix.

###  II.2 Exciton envelope function in a periodic potential

Excitons are collective excitations formed by an electron-hole bound state. We
will treat excitons as a non-interacting two-dimensional gas of bosonic quasi-
particles. This assumption is justified in the low-density regime [3], where
the excitonic density, $n_{\text{ex}}$, is much smaller than the saturation
density, $n_{\text{sat}}\approx 1/(2\pi a_{\text{2D}})$. Excitons in two-
dimensionally confined semiconductors and insulators can be treated within the
envelope function approximation, yielding the two-body function of electron
and hole [30].

In the present case of multilayered systems, we start from the factorisation
of the excitonic envelope function written as:

$F(\mathbf{r}_{e},\mathbf{r}_{h})=F_{\mathbf{k}}(\mathbf{R}_{||})f(\rho,z_{e},z_{h}),$
(11)

where $\mathbf{R}_{||}$ and $\rho$ are the in-plane exciton center of mass and
bound electron-hole pair relative coordinates, respectively. The in-plane
dynamics of the excitons can be calculated by solving the Schrödinger equation
for the envelope function:

$\left[-\frac{\hbar^{2}\nabla^{2}}{2M}+V\left(\mathbf{R}_{||}\right)\right]F_{\mathbf{k}}\left(\mathbf{R}_{||}\right)=\left(E_{\mathbf{k}}-E_{\text{ex}}\right)F_{\mathbf{k}}\left(\mathbf{R}_{||}\right),$
(12)

with $E_{\text{ex}}$ being the bare exciton energy of the active material and
$M=m^{*}_{e}+m^{*}_{h}$ the exciton mass; $V\left(\mathbf{R}_{||}\right)$ is
the in-plane potential felt by excitons, see Sec. III for a number of
examples. Similarly to the plane wave expansion of the magnetic field in the
GME method, the envelope wave function can be expanded on the same set of
plane wave modes, as:

$F_{\mathbf{k}}\left(\mathbf{R}_{||}\right)=\sum_{\mathbf{G}}F\left(\mathbf{k}+\mathbf{G}\right)e^{i\left(\mathbf{k}+\mathbf{G}\right)\cdot\mathbf{R}_{||}},$
(13)

leading to the Schrödinger-type wave equation in reciprocal space, which is
numerically easy to solve:

$\displaystyle\sum_{\mathbf{G}^{\prime}}\left[\frac{\hbar^{2}|\mathbf{k}+\mathbf{G}|^{2}}{2M}\delta_{\mathbf{G},\mathbf{G}^{\prime}}+V\left(\mathbf{G}-\mathbf{G}^{\prime}\right)\right]F\left(\mathbf{k}+\mathbf{G}^{\prime}\right)$
$\displaystyle=\left(E_{\mathbf{k}}-E_{\text{ex}}\right)F\left(\mathbf{k}+\mathbf{G}\right),$
(14)

where $V\left(\mathbf{G}-\mathbf{G}^{\prime}\right)$ is the Fourier component
of the in-plane potential. Following the diagonalization of the wave equation,
we can write a quantised excitonic Hamiltonian by introducing the annihilation
operator $\hat{b}_{\mathbf{k},\nu}$ for excitonic quasi-particles:

$H_{\text{ex}}=\sum_{\mathbf{k},\nu}E_{\mathbf{k}\nu}\hat{b}^{\dagger}_{\mathbf{k}\nu}\hat{b}_{\mathbf{k}\nu}.$
(15)

In view of formulating the generalized Hopfield matrix, the intrinsic
excitonic losses (either radiative or non-radiative) can be added a posteriori
through an imaginary term of the energies,
$\gamma_{\text{ex}}=\Gamma_{\text{ex}}/2$, where $\Gamma_{\text{ex}}$ is
intended as the resonance linewidth (i.e., the full width at half maximum),
i.e., by defining
$\tilde{E}_{\mathbf{k}\nu}=E_{\mathbf{k}\nu}-i\gamma_{\text{ex}}$ (see
Appendix).

###  II.3 Generalized Hopfield theory of exciton-photon coupling

Now that we have described the purely photonic and excitonic pictures in
second quantization, see Eqs. (9) and (15), we can take into account their
mutual interaction. Radiation-matter interaction can be described by the
minimal coupling Hamiltonian written in terms of the vector potential,
$\mathbf{A}$, and electron position ($\mathbf{r}$) and momentum ($\mathbf{p}$)
within a generalised Coulomb gauge, where
$\nabla\cdot(\varepsilon\mathbf{A})=0$:

$\displaystyle\begin{split}H_{I}=-\frac{e}{4m_{0}}\sqrt{\frac{1}{\pi\varepsilon_{0}c^{2}}}\sum_{l=1}^{L}\\{\mathbf{A}(\mathbf{r}_{l})\cdot\mathbf{p}_{l}+\mathbf{p}_{l}\cdot\mathbf{A}(\mathbf{r}_{l})\\}\\\
+\frac{e^{2}}{8\pi
m_{0}\varepsilon_{0}c^{2}}\sum_{l=1}^{L}|\mathbf{A}(\mathbf{r}_{l})|^{2}=H_{I}^{(1)}+H_{I}^{(2)}.\end{split}$
(16)

The two parts of the Hamiltonian (16) can be rewritten in terms of the
photonic and excitonic destruction (creation) operators, $\hat{a}$
($\hat{a}^{\dagger}$) and $\hat{b}$ ($\hat{b}^{\dagger}$), respectively [11]:

$\displaystyle
H_{I}^{(1)}=i\mkern-18.0mu\;\sum_{\mathbf{k},n,\nu,\sigma,j}\\!\\!\\!\\!C_{\mathbf{k}n\nu\sigma
j}(\hat{a}_{\mathbf{k}n}+\hat{a}_{-\mathbf{k}n}^{\dagger})(\hat{b}_{-\mathbf{k}\nu\sigma
j}-\hat{b}_{\mathbf{k}\nu\sigma j}^{\dagger})$ (17) $\displaystyle
H_{I}^{(2)}=\mkern-18.0mu\\!\\!\sum_{\mathbf{k},nn^{\prime},\nu,\sigma,j}\\!\\!\\!\\!\\!D_{\mathbf{k}nn^{\prime}\nu\sigma
j}(\hat{a}_{-\mathbf{k}n}+\hat{a}_{\mathbf{k}n}^{\dagger})(\hat{a}_{\mathbf{k}n^{\prime}}+\hat{a}_{-\mathbf{k}n^{\prime}}^{\dagger}),$
(18)

where $n\in\left[1,N\right]$ is the photonic band index ,
$\nu\in\left[1,M\right]$ denotes the excitonic band index. We also notice that
the present formulation generalizes the previous one [11], since we now also
consider indices $j\in\left[1,K\right]$ corresponding to the active layer
within the vertical heterostructure, and $\sigma=1,2,3$ corresponding to the
exciton polarization (along $x$, $y$, or $z$).

The coupling terms can be formally calculated as
$D_{\mathbf{k}nn^{\prime}\nu\sigma j}=C^{*}_{\mathbf{k}n\nu\sigma
j}C_{\mathbf{k}n^{\prime}\nu\sigma j}/E_{\mathbf{k},\nu,\sigma,j}$, where
$C_{\mathbf{k}n\nu\sigma j}$ is obtained as [11]:

$C_{\mathbf{k}n\nu\sigma
j}\simeq-i\left(\frac{\hbar^{2}e^{2}}{4m_{0}\varepsilon_{0}}\frac{f_{\sigma}}{S}\right)^{\frac{1}{2}}\sum_{\mathbf{G}}\hat{\mathbf{e}}_{\sigma}\cdot\mathbf{E}_{\mathbf{k}+\mathbf{G},n}(z_{j})F^{*}_{\mathbf{k}+\mathbf{G},\nu,\sigma,j}\,.$
(19)

in which $\mathbf{E}_{\mathbf{k}+\mathbf{G},n}(z_{j})$
($F^{*}_{\mathbf{k}+\mathbf{G},\nu,\sigma,j}$) is the Fourier component of the
electric field (excitonic envelope function) in the $j$-th active layer
centered in the $z_{j}$ position, while $f_{\sigma}/{S}$ is the oscillator
strength per unit surface, in general depending on the polarization. The
photon and exciton Fourier components can be retrieved, as previously shown in
Secs. II.1 and II.2.

The total Hamiltonian,
$H=H_{\text{ph}}+H_{\text{ex}}+H_{I}^{(1)}+H_{I}^{(2)}$, can be exactly
diagonalised by following a generalised Hopfield method [23, 11], i.e., by
introducing the exciton-polariton destruction bosonic operators
$p_{\mathbf{k}}$:

$\displaystyle\begin{split}p_{\mathbf{k}}&=\\!\sum_{n}w_{\mathbf{k}n}\hat{a}_{\mathbf{k}n}+\sum_{\nu\sigma
j}x_{\mathbf{k}\nu\sigma j}\hat{b}_{\mathbf{k}\nu\sigma j}\\\
&+\sum_{n}y_{\mathbf{k}n}\hat{a}^{\dagger}_{-\mathbf{k}n}+\sum_{\nu\sigma
j}z_{\mathbf{k}\nu\sigma j}\hat{b}^{\dagger}_{-\mathbf{k}\nu\sigma
j},\end{split}$ (20)

which diagonalise the new Hamiltonian with the condition
$\left[p_{\mathbf{k}},H\right]=\hbar\Omega_{\textbf{k}}p_{\mathbf{k}}$,
leading to the construction of a non-Hermitian matrix (the generalized
Hopfield operator), as described in detail in Sec. V.

As it is well known, the diagonalisation of the Hopfield matrix yields the
correct polariton eigenvalues, but not the polariton eigenstates expressed on
the exciton/photon basis [31, 25]. However, in the rotating wave approximation
(RWA, which is hereby fulfilled), the polaritonic modes can be generally
written as superposition of a photonic $\ket{\psi}_{\text{ph}}$ and excitonic
$\ket{\phi}_{\text{ex}}$ components:

$\ket{\Psi}_{\text{pol}}=\alpha\ket{\psi}_{\text{ph}}+\beta\ket{\phi}_{\text{ex}}.$
(21)

From the last equation, we can easily retrieve the excitonic ($|\beta|^{2}$)
and photonic ($|\alpha|^{2}$) fractions of the polariton modes, respectively:

$\displaystyle|\beta|^{2}=\left|\langle\phi\middle|\Psi\rangle\right|^{2},$
(22)
$\displaystyle|\alpha|^{2}=\left|\langle\psi\middle|\Psi\rangle\right|^{2},$
(23)

with the constraint $|\beta|^{2}+|\alpha|^{2}=1$ satisfied under RWA
conditions. A direct calculation of these quantities can then be obtained from
the Hopfield eigenvectors, as detailed in the Appendix, and can be very useful
for a wide range of applications, as it will be shown in the following. For
example, polariton-polariton scattering processes, which are responsible for
the polariton condensation, strongly depend on the excitonic fraction
$|\alpha|^{2}$, while for long-range transport experiments a lower excitonic
fraction is expected to yield lower intrinsic losses [3].

##  III Results

We now report a few examples to show the usefulness of this formulation of
radiation-matter interaction in patterned multilayers, which can be applied to
a variety of different systems, allowing to directly calculate the complex
eigenmodes (i.e., real and imaginary parts) of the hybrid system. We also
benchmark the results of this theoretical model through direct comparison with
the optical response calculated via rigorous coupled wave analysis (RCWA). For
the latter, we use a customized code based on the freely available S4 software
[32] that extends previous RCWA formulations [33, 34].

###  III.1  Vacuum Rabi splitting and number of QWs

Figure 2: (Colour online) (a) Polaritonic dispersion obtained by the strong
coupling between the photonic and excitonic modes for a single quantum well.
(b) Log-log plot of the Rabi splitting of the polaritonic dispersion for (a)
an increasing number of QWs. Filled circles are relative to QWs placed at the
center of the slab (i.e., in $z=0$), empty circles refer to QWs uniformly
positioned for increasing $z$ positions, starting from the slab center ($z=0$)
to the top ($z=d/2$). The slope of the simulated splittings for all QWs
assumed at the centre of the slab is simulated to be $0.500$, as expected.

As a first case study, we analyse a multi-quantum wells layer embedded in a
photonic crystal (PhC) slab. A schematic representation of this inorganic
semiconductor heterostructure is illustrated in Fig. 1. In the following
results, the dielectric slab has been assumed to be made of a material with
refractive index $n=3.43$, thickness $d=$658\text{\,}\mathrm{nm}$$, and air
holes etched in a square array (see Fig. 1). We hereby assume a lattice
constant $a=$219\text{\,}\mathrm{nm}$$, and holes radius
$r=$75\text{\,}\mathrm{nm}$$. The QWs thickness is assumed to be negligible (a
good approximation given that III-V semiconductor QWs are usually in the 10-20
nm range thick, i.e. much smaller than their typical exciton transition
wavelength), their excitonic resonance is set to
$1.485\text{\,}\mathrm{e}\mathrm{V}$ with oscillator strength per unit area
$f/S=$$4.2\text{\times}{10}^{16}\text{\,}{\mathrm{m}}^{-2}$. Similar
parameters were also considered in [11], but including a single QW. The
excitonic modes undergo strong radiation-atter coupling with the truly guided
photonic mode, ideally lossless since its dispersion falls below the light
line. This result is shown in Fig. 2(a), where N=5 photonic and M=5 excitonic
modes were kept in the basis to build the Hopfield matrix. In addition, only
in-plane polarized excitons are kept in the basis ($\sigma=1,2$), since
typical III-V semiconductor QWs possess heavy-hole fundamental excitonic
transitions. In the plot, only the polaritonic branches are evidenced, while
the exciton center of mass modes have been neglected for easier visualization.
We notice that excitonic sub-bands are almost dispersionless and quantized
within a $\mathrm{meV}$ window[11]. The Rabi splitting, i.e. the energy
separation between the two polariton branches at the resonance condition, can
be easily inferred from this plot.

Here we simulate the same structure for a varying number of QWs, $K\in[1,10]$,
in two different configurations. In a first, unrealistic situation, all the
QWs are placed at the centre of the slab, i.e. at position $z_{\text{QW}}=0$.
In such a case, all QWs interact with the same photonic mode through the same
coupling terms (19). Hence, the Rabi splitting is expected to have a
$\sqrt{K}$ dependence. The second, realistic configuration is obtained by
displacing the QWs in the upper half of the slab, as depicted in Fig. 1. In
this case a deviation from the $\sqrt{K}$ behaviour is expected. This is
indeed confirmed by our numerical solutions. In Fig. 2(b) we plot the Rabi
splitting with increasing number of QWs for both configurations, the slope in
the log-log plot reveals the predicted $\sqrt{K}$ behavior when all QWs lie at
the centre of the slab. Thus, the correct formalization of this problem allows
to correctly capture the deviation from the simplistic $\sqrt{K}$ scaling for
multi-QW structures, due to the weighted coupling of the guided mode along the
slab thickness: QWs that are next to the slab center feel an intense electric
field, whose amplitude decreases on moving towards the upper and lower edges
such that QWs positioned far from the center are less coupled.

###  III.2 Square lattice of perovskite pillars

Figure 3: (Colour online) Sketch of perovskite metasurface considered in the
main text, and experimentally studied in [20].

As a second example, we apply this theoretical framework to a 2D layered
perovskite, so-called PEPI, supporting a strong excitonic resonance and used
as an active material within a periodically textured backbone with square
lattice [20]. In this case, the metasurface is composed of a SiO2
($n_{\ch{SiO2}}=1.48$) core with etched air holes, in which the PEPI can be
infiltrated. Here we consider an etching depth
$t_{\text{etch}}=$150\text{\,}\mathrm{nm}$$. A residual PEPI layer of
thickness $t_{\text{PEPI}}=$30\text{\,}\mathrm{nm}$$ lies on top of the etched
region, as sketched in Fig 3. Finally, the device is encapsulated with a layer
of a PMMA polymer of thickness $t_{\text{PMMA}}=$200\text{\,}\mathrm{nm}$$
($n_{\text{PMMA}}=1.49$). The optical response of this active metasurface can
be calculated through RCWA, in which the excitonic response of the PEPI
material can be modeled with an isotropic Lortentzian dielectric function:

$\varepsilon_{\text{PEPI}}\left(E\right)=n^{2}+\frac{A_{\text{ex}}}{E_{\text{ex}}^{2}-E^{2}-i\Gamma_{\text{ex}}E},$
(24)

in which $n=2.41$ is the background refractive index,
$E_{\text{ex}}=$2.394\text{\,}\mathrm{eV}$$ the exciton energy,
$\Gamma_{\text{ex}}=$30\text{\,}\mathrm{meV}$$ its linewidth, and
$A_{\text{ex}}=$0.85\text{\,}{\mathrm{eV}}^{2}$$ is commonly referred to as
the “oscillator strength”. Strictly speaking, however, the oscillator strength
in either Lorentz model or quantum theory is a dimensionless quantity. For a
layered system with excitonic response, the oscillator strength per unit
surface can be related to $A_{\text{ex}}$ with the following relation [22]:

$\frac{f}{S}=\frac{\varepsilon_{0}m_{0}A_{\text{ex}}}{\hbar^{2}}\,L\,,$ (25)

where $L$ is the thickness of the active material considered.

Figure 4: (Colour online) (a) Polaritonic dispersion calculated for the
perovskite metasurface with the generalised Hopfield method. The dashed line
represents the bare excitonic energy. (b) Reflectivity spectrum calculated for
the same structure with the RCWA method.

For the polariton dispersion calculated through diagonalization of the
Hopfield matrix, we modeled the active material as a pile up of $M=15$
effective QW layers, each assumed with thickness
$L=\left(t_{\text{etch}}+t_{\text{PEPI}}\right)/15=$12\text{\,}\mathrm{nm}$$
leading to $f/S=$7.4\text{\times}{10}^{18}\text{\,}{\mathrm{m}}^{-2}$$ for
each effective active layer. Notice that in actual perovskite-based structures
the number of QW layers is much larger, which we could model exactly, in
principle; however, we realize that no appreciable change would be obtained.
Depending on the $z$-position of the effective active layer, excitons will
experience a different 2D potential profile in the $xy$-plane, which affects
the excitonic quantized energies, $E_{\mathbf{k}}$, due to different potential
terms $V\left(\mathbf{G}-\mathbf{G}^{\prime}\right)$ in Eq. (II.2). In
addition, we assume $\gamma_{\text{ex}}=\Gamma_{\text{ex}}/2=15$ meV as the
exciton imaginary part in the Hopfield matrix. For the photonic calculations
only odd modes with respect to the plane of incidence were taken into account,
namely, transverse electric (TE, i.e., s-polarised) modes. The resulting
polaritonic dispersion is shown in Fig. 4(a). In particular, the extrapolated
excitonic fraction at the bottom of the uppermost band is $0.35$, as it was
previously reported [20], which confirms the accuracy of the method.

For direct comparison and benchmarking of our model, we also report the
optical reflectivity spectra calculated through RCWA in Fig. 4(b). Guided
resonances appear as narrow Fano-like features on a reflection background,
which can be clearly evidenced in the color scale plot. The angular dispersion
of these features (which is reported as an energy-wave vector dependence) can
be put in one-to-one correspondence with the calculated polariton modes in
Fig. 4(a). In order to compare the two simulations, we slightly shifted the
photonic bands calculated by GME method and then used in the Hopfield matrix,
by assuming a PEPI material index $n=2.46$, as compared to the $n=2.41$ of the
Lorentzian response in Eq. (24) that we employed in the RCWA simulations.
Apart from the approximation of a layered sequence of 15 QWs used to mimic the
PEPI (a uniform and isotropic layer is used, instead, in the RCWA simulation),
no other adjustable parameter was employed and the two simulations are
perfectly consistent with each other. We highlight the overall agreement over
a large spectral window, which justifies a posteriori our assumption of using
15 layers, and it shows the great potential power of the Hopfield method. In
fact, numerical diagonalization of the Hopfield matrix (with $N=8$ and $M=10$)
only takes a few minutes to recover the results in Fig. 4(a) on a normal
Desktop processor. Correspondingly, RCWA simulations are known to be
computationally more demanding, and typically require about one order of
magnitude longer simulation times on comparable hardware and for comparable
parameters (number of plane waves in the Fourier basis).

Figure 5: (Colour online) (a) Cross-sectional view ($xz$ plane) of the
partially etched dielectric slab with two quantum well layers in the
unpatterned region. (b) Brillouin zone for a triangular lattice, in which
high-symmetry points have been highlighted. (c) Purely photonic band
dispersion around the $\Gamma$ point of the partially etched structure with
triangular lattice of holes, where M$/3$ (K$/3$) represents one third of the
$\Gamma$$\to$M ($\Gamma$$\to$K) high-symmetry line. A bound-state in the
continuum is found at the $\Gamma$ point of the first band.

###  III.3 Coupling to bound-states in the continuum

Figure 6: (Colour online) Left panels refer to the first photonic band of the
standard triangular lattice shown in (a), right panels refer to the same band
for the lattice with broken symmetry shown in (b). (c)-(d) Energy dispersion
around $\Gamma$ represented as color scale plot as a function of
($k_{x},k_{y}$). (e)-(f) Quality factor calculated for the bands shown in
(c)-(d); the BIC of the triangular lattice in (e) splits into two BICs in (f),
since the triangular symmetry is broken. (g)-(h) Polarisation angle of
electric field and associated topological charges. Black lines represent the
polarisation vector.

Bound-states in the continuum (BIC)s are peculiar solutions of a wave equation
that remain indefinitely confined in a resonant system despite lying in a
continuum spectrum, where modes are supposed to radiate. This paradigm has
been recently highlighted in Photonics as well [35, 36], in which non-trivial
topological properties were evidenced[37, 38]. Polariton BICs have recently
received significant attention in different physical realizations[15, 29, 20,
21, 14], also in connection to their topological properties. Recently, it has
been demonstrated that photonic BICs can be tailored by exploiting symmetry
breaking of the underlying structure [39]. Here we used this degree of freedom
to flexibly shift the BIC condition along high symmetry lines in reciprocal
space, thus enabling the strong coupling between an off-$\Gamma$ photonic BIC
and the excitonic resonance. To the best of our knowledge, this is the first
application of such concepts to polariton physics.

Figure 7: (Colour online) Lower polariton mode calculated along two symmetry
lines. Left panels refer to the standard triangular lattice of Fig. 6(a) whose
Brillouin zone is shown in (a), right panels refer to the lattice with $C_{6}$
symmetry broken (while the $C_{2}$ symmetry is preserved), i.e. corresponding
to the photonic bands of Fig. 6(b), whose Brillouin zone is shown in (b). The
vertical lines in (e)-(h) denote the wave vector for which the excitonic
fraction is $|\alpha|^{2}=0.5$. (c)-(d) Lower polariton branches obtained from
excitonic and photonic dispersions along the symmetry lines indicated in
(a)-(b). (e)-(f) Q factor of the polaritonic dispersions shown in (c)-(d); the
Q factors at excitonic fraction $|\alpha|^{2}=0.5$ are explicitly reported in
the figures. (g)-(h) Excitonic fraction corresponding to the lower polariton
bands in (c)-(d).

The starting structure is a dielectric slab partially etched with two active
QW layers embedded, as pictured in Fig. 5(a). We considered air holes etched
in a triangular lattice. The main parameters are the thicknesses of the upper
slab, QWs, barrier and lower film, respectively, which we fix as:
$t_{\text{sl}}=$132.5\text{\,}\mathrm{nm}$$,
$t_{\text{qw}}=$20\text{\,}\mathrm{nm}$$,
$t_{\text{ba}}=$10\text{\,}\mathrm{nm}$$,
$t_{\text{lw}}=$22.5\text{\,}\mathrm{nm}$$. The etching depth is assumed
$t_{\text{et}}=$110\text{\,}\mathrm{nm}$$. The lattice constant is set to
$a=$305\text{\,}\mathrm{nm}$$ and the hole radius to
$r=$106\text{\,}\mathrm{nm}$$. The refractive index of the barriers is
$n_{\text{ba}}=3.3$, while for all the other layers the assumed refractive
index is $n_{\text{qw}}=3.56$. We notice that these are typical values for
material platforms based on Al_0.4Ga_0.6As barriers and GaAs QWs in the
considered energy range. The exciton energy for such QWs is typically centered
at $E=$1.525\text{\,}\mathrm{eV}$$.

Figure 8: (Colour online) Real part of polariton eigenmodes (a) calculated by
diagonalization of the Hopfield matrix for the lattice of Fig. 7(b), and
corresponding RCWA simulation (b) of angle-resolved reflectivity spectra,
clearly showing the off-$\Gamma$ BIC.

The purely photonic dispersion calculated by GME along high-symmetry lines
around the $\Gamma$ point is shown in Fig. 5(c). For the present example, we
considered the first band at lowest energy, displaying a BIC in $\Gamma$. The
latter is clearly identified in the GME calculation, since its imaginary part
calculated by perturbation theory through Eq. 7 is identically zero. This band
is a good candidate since the at-$\Gamma$ BIC has a topological charge $q=-2$,
as defined with an integral on a closed loop $C$ in the
$\mathbf{k}_{||}=\left(k_{x},k_{y}\right)$ space around the BIC[39]:

$q=\frac{1}{2\pi}\int_{C}\nabla_{\mathbf{k}_{||}}\phi\left(\mathbf{k}_{||}\right)\cdot
d\mathbf{k}_{||}\,,$ (26)

where $\phi\left(\mathbf{k}_{||}\right)$ is the polarisation angle of the
electric field in the slab. This BIC can be split into two off-$\Gamma$ BICs,
each with a topological charge $q=-1$, by just breaking the $C_{6}$ symmetry
of the triangular lattice, i.e., its invariance under a $\pi/3$ rotation [39].
This opens up the possibility to achieve strong radiation-matter coupling
between the excitonic resonance and a photonic state with vanishing losses at
finite wave vector. As a result, a polaritonic mode with extremely high
quality (Q) factor, and consequently extremely low losses (i.e., long
lifetime), will form at the same wave vector. Notice that the Q factor is
calculated for each GME mode as
$Q_{\mathbf{k}n}=\omega_{\mathbf{k}n}/2\gamma_{\mathbf{k}n}$; hence, the Q
factor diverges at the exact BIC condition, but in the color scale plot this
is cut to an upper bound of $10^{6}$ to be displayed in log-scale.

To break the $C_{6}$ symmetry it is sufficient to slightly shift the holes
within the elementary unit cell, e.g., by modifying the angle $\theta$ defined
in Fig. 6(a). Here, the holes have been shifted in such a way to form an angle
$\theta=$62.45\text{\,}\mathrm{\SIUnitSymbolDegree}$$, as shown schematically
in Fig. 6(b). Under such conditions, the at-$\Gamma$ BIC with energy
$E\approx$1.56\text{\,}\mathrm{eV}$$, characterised by a diverging Q factor,
splits into two BICs at $k_{y}a/2\pi\approx\pm 0.073$ and energy
$E\approx$1.52\text{\,}\mathrm{eV}$$, as shown in Figs. 6(e)-(f).

Figure 9: (Colour online) Linewidth associated to the lowest energy
polaritonic band (where the polariton BIC is highlighted in Fig. 8b). In From
the RCWA simulated absorption spectra, the linewidth is extracted from
Lorentzian fits of the resonant features; in the generalised Hofpield method
the linewidth is explicitly obtained after diagonalization by taking twice the
imaginary part of the corresponding complex energy eigenvalue.

These two off-$\Gamma$ BICs can undergo strong coupling with the exciton
resonance at $E=$1.525\text{\,}\mathrm{eV}$$, thus forming off-$\Gamma$
polariton BICs. This is explicitly shown shown in Fig. 7, where we show the
results of the Hopfield matrix diagonalization for the vertical
heterostructure of Fig. 5(a) and a triangular lattice with $a=302$ nm. The
left panels refer to the lattice with
$\theta=$60\text{\,}\mathrm{\SIUnitSymbolDegree}$$, while the right ones
correspond to $\theta=$62.45\text{\,}\mathrm{\SIUnitSymbolDegree}$$.

Figures 7(c)-(d) show the lower polariton branches around $\Gamma$ and along
the symmetry directions indicated in panels (a)-(b), derived from the strong
coupling of the photonic and excitonic resonances for both configurations.
Despite the two polaritonic dispersions look very similar, they strongly
differ in their loss properties. In fact, the polariton branch at excitonic
fraction $|\alpha|^{2}=0.5$ exhibits a Q factor that is almost seven times
greater when the exciton is coupled to the BIC, as evidenced in Figs.
7(e)-(f). This implies a polariton lifetime that is almost seven times longer
for the polariton-BIC as compared to the lossy polariton mode at the same wave
vector. The off-$\Gamma$ polariton BIC is even more clearly identified by a
strongly peaked Q-factor in panel (f), whose finite value is bound by
excitonic losses in the corresponding Hopfield eigenvalue even if the purely
photonic Q-factor diverges. The corresponding excitonic fraction is plotted in
Figs. 7(g)-(h). These results are particularly relevant, since they shown that
photonic crystal polariton losses can be fully engineered within the Brillouin
zone through slight distortions of the lattice.
Finally, the topological charge conservation is demonstrated from the
polarisation maps reported in Figs. 6(g)-(h). While these results have been
obtained for the purely photonic modes, their topological properties can be
transferred one-to-one to the polariton BICs, as already shown [20].

As a final result, we show the direct comparison of Hopfield diagonalization
and RCWA simulation of the optical response of such a multilayered structure.
This is reported in Fig. 8, where the polariton dispersion (i.e., real part of
the complex energy eigenvalues) obtained from the diagonalization of the
Hopfield matrix in panel (a) are shown to be in one-to-one correspondence with
the reflectivity features in the RCWA simulated spectra in panel (b), for the
very same structure parameters of Fig. 7 corresponding to the distorted
triangular lattice. The off-$\Gamma$ polariton BIC is clearly seen in Fig.
8(b) from the vanishing reflectivity signal on top of the background. The
overall agreement is remarkable, highlighting once more the effectiveness of
the generalized Hopfield method. We just notice that a slightly different
refractive index for the QW material was assumed in the GME calculation of
purely photonic modes ($n_{\text{qw}}=3.56$), as compared to the RCWA
simulation ($n_{\text{qw}}=3.53$), to slightly adjust a few meV blue shift of
the whole spectrum that is typically obtained in GME as compared to RCWA.
The comparison between the two approaches is further benchmarked by explicitly
showing the polariton mode linewidth around the BIC condition, in Fig. 9.
Despite a small discrepancy due to the BIC condition occurring at a slightly
different point in reciprocal space, the overall agreement is remarkable. We
also notice that the Lorentzian fitting of the RCWA absorption spectra (not
shown) allows to evidence a clear saturation of losses et the off-$\Gamma$
BIC. While the fitting is not possible at the exact BIC condition (due to the
vanishing resonant response, which is due to the radiative losses becoming
identically zero), the residual loss can be evidently attributed to the non-
radiative component of excitonic losses. In the diagonalized Hopfield matrix,
the latter is automatically contained into the imaginary part of the
corresponding eigenvalue, which is due to the excitonic fraction of the lower
polariton BIC mode even when the photonic losses become identically zero. This
figure clarifies the connection between the complex mode eigenvalues obtained
within a coupled oscillators model, and the optical spectra derived from
either RCWA or experimental data.

##  IV Conclusions

We have reported a quantum theory of photonic crystal polaritons, based on a
generalization of the Hopfield matrix that takes into account the in-plane
periodic modulation of a multilayered structure containing active layers with
2D excitonic response. We have shown that this approach is extremely general
and flexible, allowing to simultaneously capture the essential physics of very
diverse material platforms and heterostructures, from inorganic multi-quantum
wells with partial etching, to infiltrated metasurfaces used as a backbone for
active media, and deposited thin layers, e.g., two-dimensional materials, on
top of patterned gratings. We have shown a few examples of application of the
method, all connected to the polariton physics and strong radiation-matter
coupling between excitonic and photonic modes, such as the correct derivation
of the Rabi splitting scaling as a function of the number of active layers
distributed along the photonic structure thickness, the dispersion of
perovskite based polaritons, and the non-trivial topological properties of
polariton bound states in the continuum. Interestingly, our results can be put
in one-to-one correspondence with quantitative information derived from the
optical spectra. We have shown that rigorous coupled wave simulations of
optical reflectivity (or absorption) match remarkably well the complex
eigenvalues obtained from diagonalizing the Hopfield matrix, by comparing
either the real part to the energy of the resonant features, or the imaginary
part to half their linewidth. This numerically fast and efficient approach
promises then to serve as a useful tool to analyse the physics of polariton
eigenmodes, to be used as an alternative or in combination with either
numerical results obtained with classical approaches, or experimental spectra
obtained by optical excitation of such heterostructures. It can be the basis
for describing nonlinear polariton interactions in photonic crystal
structures. Furthermore, the very same concepts could be extended to deal with
strongly coupled electromagnetic excitations in other contexts, from phonon to
intersubband polaritons.

###### Acknowledgements.

This work was partly supported by the Italian Ministry of Research (MUR)
through PRIN project “INPhoPOL”. H.S.N is supported by the French National
Research Agency (ANR) under the project POPEYE (ANR-17-CE24-0020) and the
IDEXLYON from Université de Lyon, Scientific Breakthrough project TORE within
the Programme Investissements d’Avenir (ANR-19-IDEX-0005). Very useful
discussions with Davide Nigro and Luca Zagaglia are gratefully acknowledged.
Some of the colormaps used in the manuscript are designed by B. R. Millis, and
they can be found at the link: https://github.com/BlakeRMills/MetBrewer.

##  V Appendix

The interaction Hamiltonian can be written in terms of photonic $\hat{a}$ and
excitonic $\hat{b}$ operators as:

$\displaystyle\begin{split}H_{\mathbf{k}}=\underbrace{\sum_{n}\hbar\omega_{\mathbf{k}n}\hat{a}_{\mathbf{k}n}^{\dagger}\hat{a}_{\mathbf{k}n}}_{H_{\text{ph}}}+\underbrace{\sum_{\nu,\sigma,j}E_{\mathbf{k}\nu\sigma
j}\hat{b}_{\mathbf{k}\nu\sigma j}^{\dagger}\hat{b}_{\mathbf{k}\nu\sigma
j}}_{H_{\text{ex}}}+\underbrace{i\sum_{n,\nu,\sigma,j}C_{\mathbf{k}n\nu\sigma
j}(\hat{a}_{\mathbf{k}n}+\hat{a}_{-\mathbf{k}n}^{\dagger})(\hat{b}_{-\mathbf{k}\nu\sigma
j}-\hat{b}_{\mathbf{k}\nu\sigma j}^{\dagger})}_{H_{\text{I}}^{(1)}}\\\
+\underbrace{\sum_{\mathbf{k},nn^{\prime},\nu,\sigma,j}D_{\mathbf{k}nn^{\prime}\nu\sigma
j}(\hat{a}_{-\mathbf{k}n}+\hat{a}_{\mathbf{k}n}^{\dagger})(\hat{a}_{\mathbf{k}n^{\prime}}+\hat{a}_{-\mathbf{k}n^{\prime}}^{\dagger})}_{H_{\text{I}}^{(2)}}.\end{split}$
(27)

where $n\in\left[1,N\right]$ is the photonic band index ,
$\nu\in\left[1,M\right]$ is the excitonic band index and
$j\in\left[1,K\right]$ is the index of the active layer. Introducing the
polaritonic operators defined as:

$p_{\mathbf{k}}=\sum_{n}w_{\mathbf{k}n}\hat{a}_{\mathbf{k}n}+\sum_{\nu\sigma
j}x_{\mathbf{k}\nu\sigma j}\hat{b}_{\mathbf{k}\nu\sigma
j}+\sum_{n}y_{\mathbf{k}n}\hat{a}^{\dagger}_{-\mathbf{k}n}+\sum_{\nu\sigma
j}z_{\mathbf{k}\nu\sigma j}\hat{b}^{\dagger}_{-\mathbf{k}\nu\sigma j}.$ (28)

The Hamiltonian can be diagonalised in terms of the polaritonic operators
$p_{\mathbf{k}}$ by imposing
$\left[p_{\mathbf{k}},H\right]=\hbar\Omega_{\textbf{k}}p_{\mathbf{k}}$ where
$\hbar\Omega_{\textbf{k}}$ are the polaritonic energies. It is convenient to
evaluate the commutator with each Hamiltonian block separately:

$\left[p_{\mathbf{k}},H_{\text{ph}}\right]=\sum_{n}\hbar\omega_{\mathbf{k}n}(w_{\mathbf{k}n}\hat{a}_{\mathbf{k}n}-y_{\mathbf{k}n}\hat{a}^{\dagger}_{-\mathbf{k}n}),\quad\left[p_{\mathbf{k}},H_{\text{ex}}\right]=\sum_{\nu,\sigma,j}E_{\mathbf{k}\nu\sigma
j}(x_{\mathbf{k}\nu\sigma j}\hat{b}_{\mathbf{k}\nu\sigma
j}-z_{\mathbf{k}\nu\sigma j}\hat{b}^{\dagger}_{-\mathbf{k}\nu\sigma j}),$ (29)
$\displaystyle\begin{split}\left[p_{\mathbf{k}},H_{I}^{(1)}\right]=i\sum_{n,\nu,\sigma,j}C_{-\mathbf{k}n\nu\sigma
j}\left\\{(w_{\mathbf{k}n}-y_{\mathbf{k}n})(\hat{b}_{\mathbf{k}\nu\sigma
j}-\hat{b}^{\dagger}_{-\mathbf{k}\nu\sigma
j})\right\\}-i\sum_{n,\nu,\sigma,j}C_{\mathbf{k}n\nu\sigma
j}\left\\{(x_{\mathbf{k}\nu\sigma j}+z_{\mathbf{k}\nu\sigma
j})(\hat{a}_{\mathbf{k}n}+\hat{a}^{\dagger}_{-\mathbf{k}n}),\right\\}\end{split}$
(30)
$\displaystyle\begin{split}\left[p_{\mathbf{k}},H_{I}^{(2)}\right]=\sum_{\nu,\sigma,j,n,n^{\prime}}\Bigl{\\{}(w_{\mathbf{k}n}-y_{\mathbf{k}n})\bigl{(}D_{-\mathbf{k}nn^{\prime}\nu\sigma
j}\hat{a}_{\mathbf{k}n^{\prime}}+D_{\mathbf{k}nn^{\prime}\nu\sigma
j}\hat{a}_{\mathbf{k}n^{\prime}}+D_{\mathbf{k}nn^{\prime}\nu\sigma
j}\hat{a}^{\dagger}_{-\mathbf{k}n^{\prime}}+D_{-\mathbf{k}nn^{\prime}\nu\sigma
j}\hat{a}^{\dagger}_{-\mathbf{k}n^{\prime}}\bigr{)}\Bigr{\\}}.\end{split}$
(31)

In the absence of losses, the Hamiltonian is Hermitian, which implies
$C_{\mathbf{k}}=C_{-\mathbf{k}}^{*}$ and $D_{\mathbf{k}}=D_{-\mathbf{k}}^{*}$,
moreover $D_{\mathbf{k}mn\nu\sigma j}=D_{\mathbf{k}nm\nu\sigma j}$. Now, in
order to simplify the linear problem
$\left[p_{\mathbf{k}},H\right]=\hbar\Omega_{\textbf{k}}p_{\mathbf{k}}$ we
calculate the following commutators:

$\displaystyle\left[\hat{a}^{\dagger}_{\mathbf{k}n},\left[p_{\mathbf{k}},H\right]\right]=\hbar\Omega_{\textbf{k}}\left[\hat{a}^{\dagger}_{\mathbf{k}n},p_{\mathbf{k}}\right],$
(32)
$\displaystyle\left[\hat{a}_{\mathbf{k}n},\left[p_{\mathbf{k}},H\right]\right]=\hbar\Omega_{\textbf{k}}\left[\hat{a}_{\mathbf{k}n},p_{\mathbf{k}}\right],$
(33) $\displaystyle\left[\hat{b}^{\dagger}_{\mathbf{k}\nu\sigma
j},\left[p_{\mathbf{k}},H\right]\right]=\hbar\Omega_{\textbf{k}}\left[\hat{b}^{\dagger}_{\mathbf{k}\nu\sigma
j},p_{\mathbf{k}}\right],$ (34)
$\displaystyle\left[\hat{b}_{\mathbf{k}\nu\sigma
j},\left[p_{\mathbf{k}},H\right]\right]=\hbar\Omega_{\textbf{k}}\left[\hat{b}_{\mathbf{k}\nu\sigma
j},p_{\mathbf{k}}\right],$ (35)

which lead to the following set of equations where $\Re$ denotes the real
part:

$\displaystyle\begin{split}w_{\mathbf{k}n}\hbar\omega_{\mathbf{k}}-i\sum_{\nu,\sigma,j}C_{\mathbf{k}\nu\sigma
nj}x_{\mathbf{k}\nu\sigma j}-i\sum_{\nu,\sigma,j}C_{\mathbf{k}\nu\sigma
nj}z_{\mathbf{k}\nu\sigma
j}+\sum_{n^{\prime},\nu,\sigma,j}2\Re\left[D_{\mathbf{k}nn^{\prime}\nu\sigma
j}\right]w_{\mathbf{k}n^{\prime}}-\sum_{n^{\prime},\nu,\sigma,j}2\Re\left[D_{\mathbf{k}nn^{\prime}\nu\sigma
j}\right]y_{\mathbf{k}n^{\prime}}=\hbar\Omega_{\mathbf{k}}w_{\mathbf{k}n},\end{split}$
(36) $\displaystyle\begin{split}x_{\mathbf{k}\nu\sigma
j}E_{\mathbf{k}\nu\sigma j}+i\sum_{n}C^{*}_{\mathbf{k}n\nu\sigma
j}w_{\mathbf{k}n\nu\sigma j}-i\sum_{n}C^{*}_{\mathbf{k}n\nu\sigma
j}y_{\mathbf{k}n\nu\sigma j}=\hbar\Omega_{\mathbf{k}}x_{\mathbf{k}\nu\sigma
j},\end{split}$ (37)
$\displaystyle\begin{split}-y_{\mathbf{k}n}\hbar\omega_{\mathbf{k}}-i\sum_{\nu,\sigma,j}C_{\mathbf{k}\nu\sigma
nj}x_{\mathbf{k}\nu\sigma j}-i\sum_{\nu,\sigma,j}C_{\mathbf{k}\nu\sigma
nj}z_{\mathbf{k}\nu\sigma
j}+\sum_{n^{\prime},\nu,\sigma,j}2\Re\left[D_{\mathbf{k}nn^{\prime}\nu\sigma
j}\right]w_{\mathbf{k}n^{\prime}}-\sum_{n^{\prime},\nu,\sigma,j}2\Re\left[D_{\mathbf{k}nn^{\prime}\nu\sigma
j}\right]y_{\mathbf{k}n^{\prime}}=\hbar\Omega_{\mathbf{k}}y_{\mathbf{k}n},\end{split}$
(38) $\displaystyle\begin{split}-z_{\mathbf{k}\nu\sigma
j}E_{\mathbf{k}\nu\sigma j}-i\sum_{n}C^{*}_{\mathbf{k}n\nu\sigma
j}w_{\mathbf{k}n\nu\sigma j}+i\sum_{n}C^{*}_{\mathbf{k}n\nu\sigma
j}y_{\mathbf{k}n\nu\sigma j}=\hbar\Omega_{\mathbf{k}}z_{\mathbf{k}\nu\sigma
j}.\end{split}$ (39)

This set of equations can be written as a linear eigenvalue problem with a
generalised Hopfield matrix:

$M=\begin{pmatrix}\bm{\omega}+2\Re\mathbf{D}&-i\mathbf{C}_{1}&-i\mathbf{C}_{2}&\dots&-i\mathbf{C}_{K}&-2\Re\mathbf{D}&-i\mathbf{C}_{1}&-i\mathbf{C}_{2}&\dots&-i\mathbf{C}_{K}\\\
i\mathbf{C}_{1}^{\dagger}&\mathbf{E}_{1}&\mathbf{0}&\dots&\mathbf{0}&-i\mathbf{C}_{1}^{\dagger}&\mathbf{0}&\mathbf{0}&\dots&\mathbf{0}\\\
i\mathbf{C}_{2}^{\dagger}&\mathbf{0}&\mathbf{E}_{2}&\dots&\mathbf{0}&-i\mathbf{C}_{2}^{\dagger}&\mathbf{\mathbf{0}}&\mathbf{0}&\dots&\mathbf{0}\\\
\vdots&\vdots&\vdots&\ddots&\mathbf{0}&\vdots&\vdots&\vdots&\ddots&\mathbf{0}\\\
i\mathbf{C}_{K}^{\dagger}&\mathbf{0}&\mathbf{0}&\dots&\mathbf{E}_{K}&-i\mathbf{C}_{K}^{\dagger}&\mathbf{0}&\mathbf{0}&\dots&\mathbf{0}\\\
2\Re\mathbf{D}&-i\mathbf{C}_{1}&-i\mathbf{C}_{2}&\dots&-i\mathbf{C}_{K}&-\bm{\omega}-2\Re\mathbf{D}&-i\mathbf{C}_{1}&-i\mathbf{C}_{2}&\dots&-i\mathbf{C}_{K}\\\
-i\mathbf{C}_{1}^{\dagger}&\mathbf{0}&\mathbf{0}&\dots&\mathbf{0}&+i\mathbf{C}_{1}^{\dagger}&-\mathbf{E}_{1}&\mathbf{0}&\dots&\mathbf{0}\\\
-i\mathbf{C}_{2}^{\dagger}&\mathbf{0}&\mathbf{0}&\dots&\mathbf{0}&+i\mathbf{C}_{2}^{\dagger}&\mathbf{0}&-\mathbf{E}_{2}&\dots&\mathbf{0}\\\
\vdots&\vdots&\vdots&\ddots&\mathbf{0}&\vdots&\vdots&\vdots&\ddots&\mathbf{0}\\\
-i\mathbf{C}_{K}^{\dagger}&\mathbf{0}&\mathbf{0}&\dots&\mathbf{0}&+i\mathbf{C}_{K}^{\dagger}&\mathbf{0}&\mathbf{0}&\dots&-\mathbf{E}_{K}\\\
\end{pmatrix},$ (40)

whose eigenvalues correspond to the polaritonic energies. Note that that j-th
quantum well is coupled to the photonic modes by the $\mathbf{C}_{j}$ block.

Losses can be included a posteriori by assuming complex photon and exciton
energies. Specifically, the blocks composing the generalised Hopfield matrix
(40) including losses are: the $N\times N$ matrix (photon losses included as
imaginary parts of the eigenvalues, as discussed in Sec. II.1)

$\bm{\omega}+2\Re\mathbf{D}=[\\![\hbar\tilde{\omega}_{n}\delta_{n,n^{\prime}}+\sum_{\nu,\sigma,j}2\Re\left[D_{nn^{\prime}\nu\sigma
j}\right]]\\!],$ (41)

where $\tilde{\omega}_{n}=\omega_{n}-i\gamma_{n}$, the $3M\times 3M$ matrices
(exciton losses included as imaginary parts of the exciton energies)

$\mathbf{E}_{j}=[\\![\tilde{E}_{\nu\sigma}^{j}\delta_{(\nu\sigma),(\nu^{\prime}\sigma^{\prime})}]\\!],$
(42)

where $\tilde{E}_{\nu\sigma}^{j}=E_{\nu\sigma}^{j}-i\gamma_{\mathrm{ex}}$, the
$N\times 3M$ matrices

$\mathbf{C}_{j}=[\\![C_{n,\nu,\sigma}^{j}]\\!],$ (43)

and the zero-filled $3M\times 3M$ matrix denoted as $\mathbf{0}$. The notation
of the Hopfield matrix (40) can be further simplified by introducing
$\mathbf{\Gamma}$ and $\mathbf{\Xi}$ defined as:

$\mathbf{\Gamma}=\begin{bmatrix}\mathbf{C}_{1}&\mathbf{C}_{2}&\dots&\mathbf{C}_{K}\\\
\end{bmatrix},\quad\mathbf{\Xi}=\begin{bmatrix}\mathbf{E}_{1}&\mathbf{0}&\dots&\mathbf{0}\\\
\mathbf{0}&\mathbf{E}_{2}&\dots&\mathbf{0}\\\
\vdots&\vdots&\ddots&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}&\dots&\mathbf{E}_{K}\\\ \end{bmatrix},$ (44)

resulting in the following simple form:

$M=\begin{pmatrix}\bm{\omega}+2\Re\mathbf{D}&-i\mathbf{\Gamma}&-2\Re\mathbf{D}&-i\mathbf{\Gamma}\\\
i\mathbf{\Gamma}^{\dagger}&\mathbf{\Xi}&-i\mathbf{\Gamma}^{\dagger}&\mathbf{\tilde{0}}\\\
2\Re\mathbf{D}&-i\mathbf{\Gamma}&-\bm{\omega}-2\Re\mathbf{D}&-i\mathbf{\Gamma}\\\
-i\mathbf{\Gamma}^{\dagger}&\mathbf{\tilde{0}}&i\mathbf{\Gamma}^{\dagger}&-\mathbf{\Xi}\\\
\end{pmatrix},$ (45)

where $\mathbf{\tilde{0}}$ is the zero-filled $3M\cdot K\times 3M\cdot K$
matrix.

We notice that the excitonic ($(|\beta|^{2})$) and photonic ($(|\alpha|^{2})$)
single polariton fractions can be calculated from the normalized eigenvectors
of the Hopfield matrix, expressed as follows:

$\ket{v}=\begin{pmatrix}v_{1}\\\ \vdots\\\ v_{N}\\\ \vdots\\\ v_{N+3M\cdot
K}\\\ \vdots\\\ v_{2N+3M\cdot K}\\\ \vdots\\\ v_{2N+6M\cdot
K}\end{pmatrix}\,,$ (46)

from which we define:

$|\alpha|^{2}={\sum}_{i=1}^{N}|v_{i}|^{2}+{\sum}_{i=N+3M\cdot K+1}^{2N+3M\cdot
K}|v_{i}|^{2}$ (47) $|\beta|^{2}={\sum}_{i=N+1}^{N+3M\cdot
K}|v_{i}|^{2}+{\sum}_{i=2N+3M\cdot K+1}^{2N+6M\cdot K}|v_{i}|^{2}\,.$ (48)

The latter accurately reproduce the single exciton and single photon fractions
of the polariton eigenstates only when RWA can be assumed [31, 25].

## References

  * [1] D. N. Basov, A. Asenjo-Garcia, P. J. Schuck, X. Zhu, and A. Rubio, “Polariton panorama,” Nanophotonics, vol. 10, no. 1, pp. 549–577, 2021.
  * [2] D. Sanvitto and S. Kéna-Cohen, “The road towards polaritonic devices,” Nature Materials, vol. 15, pp. 1061–1073, jul 2016.
  * [3] I. Carusotto and C. Ciuti, “Quantum fluids of light,” Rev. Mod. Phys., vol. 85, pp. 299–366, Feb 2013.
  * [4] J. Kasprzak, M. Richard, S. Kundermann, A. Baas, P. Jeambrun, J. M. J. Keeling, F. M. Marchetti, M. H. Szymańska, R. André, J. L. Staehli, V. Savona, P. B. Littlewood, B. Deveaud, and L. S. Dang, “Bose–Einstein condensation of exciton polaritons,” Nature, vol. 443, pp. 409–414, sep 2006.
  * [5] A. Amo, J. Lefrère, S. Pigeon, C. Adrados, C. Ciuti, I. Carusotto, R. Houdré, E. Giacobino, and A. Bramati, “Superfluidity of polaritons in semiconductor microcavities,” Nature Physics, vol. 5, pp. 805 – 810, 2009.
  * [6] S. Azzini, D. Gerace, M. Galli, I. Sagnes, R. Braive, A. Lemaître, J. Bloch, and D. Bajoni, “Ultra-low threshold polariton lasing in photonic crystal cavities,” Applied Physics Letters, vol. 99, no. 11, p. 111106, 2011.
  * [7] D. Ballarini, M. De Giorgi, E. Cancellieri, R. Houdré, E. Giacobino, R. Cingolani, A. Bramati, G. Gigli, and D. Sanvitto, “All-optical polariton transistor,” Nature Communications, vol. 4, p. 1778, jun 2013.
  * [8] H. S. Nguyen, D. Vishnevsky, C. Sturm, D. Tanese, D. Solnyshkov, E. Galopin, A. Lemaître, I. Sagnes, A. Amo, G. Malpuech, and J. Bloch, “Realization of a double-barrier resonant tunneling diode for cavity polaritons,” Phys. Rev. Lett., vol. 110, p. 236601, Jun 2013.
  * [9] F. Marsault, H. S. Nguyen, D. Tanese, A. Lemaître, E. Galopin, I. Sagnes, A. Amo, and J. Bloch, “Realization of an all optical exciton-polariton router,” Applied Physics Letters, vol. 107, p. 201115, nov 2015.
  * [10] D. Nigro, V. D’Ambrosio, D. Sanvitto, and D. Gerace, “Integrated quantum polariton interferometry,” Communications Physics, vol. 5, no. 1, p. 34, 2022.
  * [11] D. Gerace and L. C. Andreani, “Quantum theory of exciton-photon coupling in photonic crystal slabs with embedded quantum wells,” Phys. Rev. B, vol. 75, p. 235325, Jun 2007.
  * [12] D. Bajoni, D. Gerace, M. Galli, J. Bloch, R. Braive, I. Sagnes, A. Miard, A. Lemaître, M. Patrini, and L. C. Andreani, “Exciton polaritons in two-dimensional photonic crystals,” Phys. Rev. B, vol. 80, p. 201308(R), Nov 2009.
  * [13] C. E. Whittaker, T. Isoniemi, S. Lovett, P. M. Walker, S. Kolodny, V. Kozin, I. V. Iorsh, I. Farrer, D. A. Ritchie, M. S. Skolnick, and D. N. Krizhanovskii, “Exciton–polaritons in gaas-based slab waveguide photonic crystals,” Applied Physics Letters, vol. 119, no. 18, p. 181101, 2021.
  * [14] V. Ardizzone, F. Riminucci, S. Zanotti, A. Gianfrate, M. Efthymiou-Tsironi, D. G. Suàrez-Forero, F. Todisco, M. D. Giorgi, D. Trypogeorgos, G. Gigli, K. Baldwin, L. Pfeiffer, D. Ballarini, H. S. Nguyen, D. Gerace, and D. Sanvitto, “Polariton bose–einstein condensate from a bound state in the continuum,” Nature, vol. 605, pp. 447–452, May 2022.
  * [15] K. L. Koshelev, S. K. Sychev, Z. F. Sadrieva, A. A. Bogdanov, and I. V. Iorsh, “Strong coupling between excitons in transition metal dichalcogenides and optical bound states in the continuum,” Phys. Rev. B, vol. 98, p. 161113(R), Oct 2018.
  * [16] L. Zhang, R. Gogna, W. Burg, E. Tutuc, and H. Deng, “Photonic-crystal exciton-polaritons in monolayer semiconductors,” Nature Communications, vol. 9, p. 713, Feb. 2018.
  * [17] R. Gogna, L. Zhang, Z. Wang, and H. Deng, “Photonic crystals for controlling strong coupling in van der waals materials,” Opt. Express, vol. 27, pp. 22700–22707, Aug 2019.
  * [18] Y. Chen, S. Miao, T. Wang, D. Zhong, A. Saxena, C. Chow, J. Whitehead, D. Gerace, X. Xu, S.-F. Shi, and A. Majumdar, “Metasurface integrated monolayer exciton polariton,” Nano Letters, vol. 20, no. 7, pp. 5292–5300, 2020.  PMID: 32519865.
  * [19] T. Fujita, Y. Sato, T. Kuitani, and T. Ishihara, “Tunable polariton absorption of distributed feedback microcavities at room temperature,” Phys. Rev. B, vol. 57, pp. 12428–12434, May 1998.
  * [20] N. H. M. Dang, S. Zanotti, E. Drouard, C. Chevalier, G. Trippé-Allard, M. Amara, E. Deleporte, V. Ardizzone, D. Sanvitto, L. C. Andreani, C. Seassal, D. Gerace, and H. S. Nguyen, “Realization of polaritonic topological charge at room temperature using polariton bound states in the continuum from perovskite metasurface,” Advanced Optical Materials, vol. 10, no. 6, p. 2102386, 2022.
  * [21] S. Kim, B. H. Woo, S.-C. An, Y. Lim, I. C. Seo, D.-S. Kim, S. Yoo, Q.-H. Park, and Y. C. Jun, “Topological control of 2d perovskite emission in the strong coupling regime,” Nano Letters, vol. 21, pp. 10076–10085, 12 2021.
  * [22] L. C. Andreani, in “Strong light-matter coupling: from atoms to solid state systems,” (World Scientific, Singapore, 2013), pp. 37-82.
  * [23] J. J. Hopfield, “Theory of the contribution of excitons to the complex dielectric constant of crystals,” Phys. Rev., vol. 112, pp. 1555–1567, Dec 1958.
  * [24] V. M. Agranovich, “Dispersion of electromagnetic waves in crystals,” Sov. Phys. JETP, vol. 37, p. 307, February 1960.
  * [25] C. Ciuti, G. Bastard, and I. Carusotto, “Quantum vacuum properties of the intersubband cavity polariton field,” Phys. Rev. B, vol. 72, p. 115303, Sep 2005.
  * [26] D. Hagenmüller, S. De Liberato, and C. Ciuti, “Ultrastrong coupling between a cavity resonator and the cyclotron transition of a two-dimensional electron gas in the case of an integer filling factor,” Phys. Rev. B, vol. 81, p. 235303, Jun 2010.
  * [27] L. C. Andreani and D. Gerace, “Photonic-crystal slabs with a triangular lattice of triangular holes investigated using a guided-mode expansion method,” Phys. Rev. B, vol. 73, p. 235114, Jun 2006.
  * [28] M. Minkov, I. A. D. Williamson, L. C. Andreani, D. Gerace, B. Lou, A. Y. Song, T. W. Hughes, and S. Fan, “Inverse design of photonic crystals through automatic differentiation,” ACS Photonics, vol. 7, no. 7, pp. 1729–1741, 2020.
  * [29] V. Kravtsov, E. Khestanova, F. A. Benimetskiy, T. Ivanova, A. K. Samusev, I. S. Sinev, D. Pidgayko, A. M. Mozharov, I. S. Mukhin, M. S. Lozhkin, Y. V. Kapitonov, A. S. Brichkin, V. D. Kulakovskii, I. A. Shelykh, A. I. Tartakovskii, P. M. Walker, M. S. Skolnick, D. N. Krizhanovskii, and I. V. Iorsh, “Nonlinear polaritons in a monolayer semiconductor coupled to optical bound states in the continuum,” Light: Science & Applications, vol. 9, p. 56, dec 2020.
  * [30] G. Panzarini and L. C. Andreani, “Quantum theory of exciton polaritons in cylindrical semiconductor microcavities,” Phys. Rev. B, vol. 60, pp. 16799–16806, Dec 1999.
  * [31] A. Quattropani, L. C. Andreani, and F. Bassani, “Quantum theory of polaritons with spatial dispersion: Exact solutions,” Il Nuovo Cimento D, vol. 7, no. 1, pp. 55–69, 1986.
  * [32] V. Liu and S. Fan, “S4 : A free electromagnetic solver for layered periodic structures,” Computer Physics Communications, vol. 183, no. 10, pp. 2233 – 2244, 2012.
  * [33] D. M. Whittaker and I. S. Culshaw, “Scattering-matrix treatment of patterned multilayer photonic structures,” Phys. Rev. B, vol. 60, pp. 2610–2618, Jul 1999.
  * [34] M. Liscidini, D. Gerace, L. C. Andreani, and J. E. Sipe, “Scattering-matrix analysis of periodically patterned multilayers with asymmetric unit cells and birefringent media,” Phys. Rev. B, vol. 77, p. 035324, Jan 2008.
  * [35] C. W. Hsu, B. Zhen, A. D. Stone, J. D. Joannopoulos, and M. Soljačić, “Bound states in the continuum,” Nature Reviews Materials, vol. 1, no. 18, p. 16048, 2016.
  * [36] S. I. Azzam and A. V. Kildishev, “Photonic bound states in the continuum: From basics to applications,” Advanced Optical Materials, vol. 9, no. 1, p. 2001469, 2021.
  * [37] B. Zhen, C. W. Hsu, L. Lu, A. D. Stone, and M. Soljačić, “Topological Nature of Optical Bound States in the Continuum,” Physical Review Letters, vol. 113, no. 257401, 2014.
  * [38] Y. Zhang, A. Chen, W. Liu, C. W. Hsu, B. Wang, F. Guan, X. Liu, L. Shi, L. Lu, and J. Zi, “Observation of polarization vortices in momentum space,” Phys. Rev. Lett., vol. 120, p. 186103, May 2018.
  * [39] T. Yoda and M. Notomi, “Generation and annihilation of topologically protected bound states in the continuum and circularly polarized states by symmetry breaking,” Phys. Rev. Lett., vol. 125, p. 053902, Jul 2020.