# Theory of exciton-polariton condensation in gap-confined eigenmodes

Davide Nigro  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

Exciton-polaritons are bosonic-like elementary excitations in semiconductors,
which have been recently shown to display large occupancy of topologically
protected polariton bound states in the continuum in suitably engineered
photonic lattices [Nature 605, 447 (2022)], compatible with the definition of
polariton condensation. However, a full theoretical description of such
condensation mechanism that is based on a non equilibrium Gross-Pitaevskii
formulation is still missing. Given that the latter is well known to account
for polariton condensation in conventional semiconductor microcavities, here
we report on its multi-mode generalization, showing that it allows to fully
interpret the recent experimental findings in patterned photonic lattices,
including emission characteristics and condensation thresholds. Beyond that,
it is shown that the polariton condensation in these systems is actually the
result of an interplay between negative mass confinement of polariton
eigenstates (e.g., due to the photonic gap originated from the periodic
pattern in plane) and polariton losses. We are then able to show that
polariton condensation can also occur in gap-confined bright modes, i.e.,
coupling of QW excitons to a dark photonic mode is not necessarily required to
achieve a macroscopic occupation with low population threshold.

##  I Introduction

Exciton-polaritons arising in low-dimensional nanostructures, such as coupled
quantum well (QW) excitons and confined photonic modes, have been shown to
behave as a weakly interacting Bose gas [1]. Condensation phenomena in these
systems have to be inevitably described in terms of a balance between driving
and losses, due to photon leakage or exciton recombination. In this respect, a
suitable generalization of the standard nonlinear Schrödinger equation, known
as the non-equilibrium Gross-Pitaevskii equation (NEGPE), is able to largely
account for the observed phenomenology including external driving and
intrinsic losses, either under resonant or non-resonant pumping of the
exciton-polariton field [2]. In particular, spectacular phenomena such as
Bose-Einstein condensation [3, 4], superfluidity [5], formation of quantized
vortices [6] have been observed in planar semiconductor microcavities, in
which the photon field is confined in the QW plane between two high-
reflectivity Bragg mirrors [7]. Besides showing interest for fundamental
physics studies, exciton-polaritons have also been shown to potentially allow
for useful applications due to their superior nonlinear properties, such as
low-threshold lasers [8], all-optical switching [9], sensing [10], etc.

Figure 1: Schematic representation of the multi-QW heterostructure with
surface periodic patterning that creates a photonic gap for in-plane
propagating modes. Elementary excitations can be created by a laser spot
$P(x)$, which locally changes the refractive index in the layers underneath
($n_{1}$ and $n_{2}$), thus inducing an effective confining potential for the
negative mass polaritons.

Recently, polariton condensation has been observed also in multi-layered QW
heterostructures with periodic surface patterning, such as the one sketched in
Fig. 1. In fact, it has been shown that top/bottom mirrors are not necessary
to achieve the polariton condensation threshold, since out-of plane radiation
losses associated to the photonic component of the polariton field can be
fully engineered by the photonic lattice [11], even a one-dimensional one
[12]. In particular, radiative losses of the photonic component can be fully
suppressed by symmetry along certain directions, leading to the currently
widespread concept of photonic bound states in the continuum (BIC) [13, 14].
When coupled to QW excitons, these modes lead to the formation of polariton
eigenmodes that are dark for emission along specific spatial directions and at
fixed energy [15], inheriting the topologically protected nature of the purely
photonic BIC [16, 17, 18]. Polariton condensation has been observed in these
systems by efficient relaxation into modes arising from the saddle-like
dispersion of the lower polariton branch, i.e., not an absolute minimum of the
dispersion relation [19, 20]. In these works, the condensation mechanism has
been qualitatively interpreted as due to the creation of an effective
potential well for negative mass polaritons, induced by the exciting laser
spot over a finite width along the periodicity lattice, as well as to the fact
that such confined states acquire a longer lifetime than other modes due to
coupling to purely photonic BICs. However, no quantitative theory has been
provided to account for the condensation in such BIC-like modes, to date. In
particular, a rigorous discussion about the role of the negative mass branch
in connection with the photonic-induced gap obtained as a consequence of the
periodic refractive index modulation is still missing. This manuscript aims at
filling this theoretical gap, by developing a model of out-of-equilibrium
condensation in a multi-band exciton-photon coupled system with symmetry-
dependent losses. In particular, we hereby unravel the critical role played by
the negative mass and the energy gap between dark and bright modes in reaching
the condensation onset.

Here is a short outline of the results presented in the paper. In Sec. II we
describe the dynamical equations used to picture the onset of polariton
condensation in the continuous-wave (CW) regime of a model Hamiltonian with a
given number of photonic branches coupled to QW exciton modes. Then, in Sec.
III we address in detail the behavior of our effective theory in a specific
case, where the relevant physics can be traced back to the presence of two
counter-propagating photonic modes. In order to understand the polariton
condensation in this theoretical framework, we first characterize the
polariton dispersion arising from the Hamiltonian diagonalization, by deriving
polariton bands in Sec. III.1, and the effects of an external space-dependent
potential $V(x)$ coupled to the particle density in Sec. III.2. In Sec. IV we
show the results concerning the polariton condensation in discrete energy
eigenvalues lying within the energy gap, which is ultimately due to the
periodic photonic modulation. These eigenvalues correspond to spatially
confined eigenmodes below the potential $V(x)$. In particular, after showing
the results of our numerical simulations describing the behavior of the
condensate population in the steady-state, we provide an analysis of the
condensate loss rate (Sec. IV.1), as well as the spectral density describing
the energy-momentum behavior of the light emitted by the condensate (Sec.
IV.2). The behavior of such quantities in all the cases considered in the
present work evidently supports the conjecture that polaritons are condensing
into the discrete gap-confined eigenmodes previously characterized in Sec.
III.2. Depending on the dark or bright nature of the negative mass branch out
of which they originate, these modes display very different direct and
reciprocal space profiles. In light of these results, in Sec. V we finally
discuss the relevance of our findings and future perspectives.

##  II Theoretical model

In this section we introduce a general dynamical model that can be used to
describe and characterize the onset of polariton condensation and its relation
to the eigenstates properties in heterostructures like the one sketched in
Fig. 1. We will assume CW pumping in this work, although the model is general.
The system is assumed to be uniform in the direction transverse to $x$,
consisting of several layers with different refractive indices, in the figure
denoted as $n_{1}$ and $n_{2}$. For instance [19, 20], such a stack of
different index material might correspond to a set of GaAs/AlxGa1-xAs layers,
which give rise to the formation of quantum well (QW) excitons in the lower
band-gap material [21, 22]. Due to the top surface patterning with spatial
periodicity $a$, the free-propagating electromagnetic modes get folded,
leading to the formation of gapped Bloch resonances [11, 12]. When coupled to
QW excitons, these photonic modes are shown to give rise to the concept of
photonic crystal polaritons [23], which have been evidenced in different
systems and material platforms [24, 25, 26, 27], all of them falling within
the same theoretical treatment [16]. In particular, it has been evidenced that
a surface patterning is sufficient to induce gapped polariton branches with
positive and negative effective masses around normal incidence and energies
below the bare exciton resonance [19, 20]. It has been observed that polariton
condensation occurs within a set of discretized quantum states appearing
within such energy gaps. Their emission properties are shown to depend both on
the spatial profile of the input light-source [19, 28], $P(x)$ in Fig. 1, and
on the intrinsic emission features of the lower-lying polariton branch
delimiting the band-gap. Similarly to the description of polariton
condensation in the planar microcavity case, here it is reasonable to assume
that particle scattering towards the condensed eigenstate can be described by
means of an exciton reservoir non-linearly coupled to polariton states [29,
30, 31, 32]. However, in the present case it is also crucial to account for
the mechanism leading to the formation of such discrete states within the
polariton gap, as well as the peculiar emission properties shown by either the
bare polariton bands or the polariton condensate modes. Moreover, the presence
of different photonic modes may lead to some of them being in weak coupling
with the exciton, but still participating in the relaxation dynamics and thus
worth being taken into account. In this respect, as it is shown in the
following, our approach seems to provide a quite complete picture allowing to
account for the observed phenomenology, as well as predicting possible new
outcomes, when including all these crucial aspects.

In our model, similarly to the standard approach adopted to describe
condensation in microcavity systems [29, 33, 34], we assume $P(x)$ to be
incoherently coupled to a reservoir population, $n(x,\,t)$, describing high-
energy exciton states. By stimulated scattering, such a reservoir feeds lower-
energy eigenmodes, thus replenishing the exciton-polariton population in these
states. However, in order to capture the complex band-gap physics, emission
properties as well as the formation of discrete levels, instead of using a
single component wavefunction $\psi(x,\,t)$ (usually describing the lower
branch polariton field), here we assume that the reservoir density couples to
a set of relevant bare photonic and excitonic modes of the structure. In other
words, the system dynamics is hereby described by means of a multi-component
vector, $\overrightarrow{\psi}(x,\,t)$ formally reading

$\overrightarrow{\psi}(x,\,t)=\left(A_{1},A_{2},\cdots,A_{L},X_{1},X_{2},\cdots,X_{L}\right)^{T},$
(1)

in which $A_{l}\equiv A_{l}(x,\,t)$ and $X_{l}\equiv X_{l}(x,t)$
($l\in[1,2,\,\cdots,\,L]$) denote the wavefunction of the $l$-th photonic and
$l$-th excitonic mode respectively. In particular, $2L$ is the total number of
polariton branches participating in the system dynamical evolution.
By setting $n\equiv n(x,\,t)$,
$\overrightarrow{\psi}\equiv\overrightarrow{\psi}(x,\,t)$, we assume the
reservoir-polariton dynamics to be described by the following set of coupled
differential equations

$\displaystyle[left={\,\empheqlbrace}$
$\displaystyle\frac{d}{dt}n=P(x)-\frac{1}{\tau_{R}}\,n-gn\langle\overrightarrow{\psi},\,G\overrightarrow{\psi}\rangle$
(2a)
$\displaystyle\frac{d}{dt}\overrightarrow{\psi}=\tilde{H}(x)\overrightarrow{\psi}+\frac{g}{2}n\,G\,\overrightarrow{\psi}$
(2b)

in which $\langle\vec{w},\,\vec{z}\rangle=\sum_{j}w_{j}^{*}z_{j}$ denotes the
hermitian scalar product between the two vectors $\vec{w}$ and $\vec{z}$, with
$w_{j}^{*}$ denoting the complex-conjugated of the component $w_{j}$.

The different terms in Eq. (2a) account for three main physical contributions:
(i) the injection of population due to the CW driving through $P(x)$; (ii) the
presence of a loss term controlled by the intrinsic finite-lifetime $\tau_{R}$
of the reservoir; (iii) a non-linear decay term proportional to the effective
coupling $g>0$, that describes the stimulated population scattering from the
reservoir to photon and exciton states. The $2L\times 2L$ matrix $G$
represents an operator satisfying the following constraints:

$\displaystyle[left={\,\empheqlbrace}]$ $\displaystyle G^{\dagger}=G$
(hermiticity) (2ca) $\displaystyle G\geq 0$ (positive semi-definite) (2cb)
$\displaystyle\sum_{j}G_{jj}=1$ (unit trace) (2cc)

Equations (2ca) and (2cb) ensure that the reservoir density evolves as a real
quantity, with last term acting as an effective loss (i.e., gain on the
$\vec{\psi}$ field, as clear from Eq. (2b)). The normalization constraint in
Eq. (2cc) guarantees for a particle-conserving dynamics. Indeed, if on the one
hand $g$ controls the global scattering rate from the reservoir, on the other
the model accounts for several decay channels, each controlled by an element
of the matrix $G$. Since a particle scattered away from the reservoir must
create a single particle into the light-matter part of the model, the
$\\{G_{jj}\\}$ must add up to 1.

The coupled exciton-photon dynamics is encoded into Eq. (2b). The term
proportional to the (non-Hermitian) $2L\times 2L$ matrix operator
$\tilde{H}(x)$, given by

$\tilde{H}(x)=-i\frac{E(x)}{\hbar}-\Gamma(x)-i\frac{W(x)}{\hbar}$ (2d)

accounts for the peculiar complex energy-momentum dispersion laws of the
relevant bare modes and their mutual linear interaction, i.e.,
$-iE(x)/\hbar-\Gamma(x)$, as well as the presence of an external space-
dependent potential, $W(x)$, such that $W(x)=W^{\dagger}(x)$. In particular,
for the same reason leading to Eq. (2c), in order to interpret $\Gamma(x)$ as
a proper decay operator, one must require $\Gamma(x)$ to satisfy the following
constraints:

$\displaystyle[left={\,\empheqlbrace}]$
$\displaystyle\Gamma^{\dagger}(x)=\Gamma(x)$ (hermiticity) (2ea)
$\displaystyle\Gamma(x)\geq 0$ (positive semi-definite) (2eb)

Finally, the last term in Eq. (2b) accounts for the particle transfer from the
reservoir to photonic and excitonic modes. In particular, the structure of the
nonlinear terms proportional to $g$ in Eq. (2), and the constraints imposed on
$G$, lead to the expected $g$-independent behavior of the time derivative of
the total particle density, $N(x,\,t)\equiv
n(x,\,t)+\langle\overrightarrow{\psi}(x,\,t),\,\overrightarrow{\psi}(x,\,t)\rangle$.
By combining the two expressions in Eq. (2), the following rate equation for
$N(x,t)$ is obtained

$\begin{split}&\frac{d}{dt}N(x,\,t)=\frac{d}{dt}n(x,\,t)+\frac{d}{dt}\langle\overrightarrow{\psi}(x,\,t),\,\overrightarrow{\psi}(x,\,t)\rangle=\\\
&=P(x)-\frac{1}{\tau_{R}}\,n-\langle\tilde{H}(x)\overrightarrow{\psi},\,\overrightarrow{\psi}\rangle+\langle\overrightarrow{\psi},\,\tilde{H}(x)\overrightarrow{\psi}\rangle.\\\
\end{split}$ (2f)

In the long-time limit, since the CW source does not depend on time (i.e., the
only driven Fourier component of the reservoir density is the $\omega=0$ one),
it is reasonable to assume that the dynamical system does not produce
oscillating solutions (limit cycles). As a consequence, the steady-state (ss)
configuration, i.e., the solutions satisfying
$\frac{d}{dt}n(x,\,t)=\frac{d}{dt}\overrightarrow{\psi}(x,\,t)=0$, correspond
to a fixed point of the dynamical system in Eq. (2). In other words, by
recalling that due to Eq. (2cb) $\langle\vec{\psi},G\vec{\psi}\rangle\geq
0\,\forall\,\vec{\psi}$, we have the steady state relations

$\displaystyle[left={\,\empheqlbrace}]$ $\displaystyle
n_{ss}(x)=\frac{P(x)}{\frac{1}{\tau_{R}}+g\langle\overrightarrow{\psi}_{ss}(x),\,G\overrightarrow{\psi}_{ss}(x)\rangle}$
(2ga)
$\displaystyle\tilde{H}(x)\overrightarrow{\psi}_{ss}(x)=-\frac{g}{2}\frac{P(x)G\,\overrightarrow{\psi}_{ss}(x)}{\frac{1}{\tau_{R}}+g\langle\overrightarrow{\psi}_{ss}(x),\,G\overrightarrow{\psi}_{ss}(x)\rangle}\,$
(2gb)

in which $n_{ss}$ and $\overrightarrow{\psi}_{ss}$ denote the steady-state
reservoir density and exciton-photon configuration, respectively. As it is
easy to verify by direct inspection, the system of coupled equations (2g)
always admits the trivial configuration corresponding to an empty exciton-
photon subsystem and a reservoir density proportional to the CW profile as a
fixed point solution, i.e.,

$\overrightarrow{\psi}_{ss}=0,\quad n_{ss}(x)=\tau_{R}P(x).$ (2h)

Depending on the pump strength $P_{0}$ (see Fig. 1), we notice that such a
solution might become linearly unstable, i.e., a small perturbation
$\delta\overrightarrow{\psi}(x)$ placed on top of
$\overrightarrow{\psi}_{ss}=0$ gives rise to a macroscopic increasing of the
condensate density, i.e.,
$\langle\delta\overrightarrow{\psi}(x),\delta\overrightarrow{\psi}(x)\rangle$,
pushing the system away from the trivial configuration in Eq. (2h). By
linearizing the dynamical system in Eq. (2), we obtain the following
expression for the exciton-photon density evolution:

$\begin{split}\frac{d}{dt}\langle\delta\overrightarrow{\psi}(x),\delta\overrightarrow{\psi}(x)\rangle&=-2\langle\delta\overrightarrow{\psi}(x),\Gamma(x)\delta\overrightarrow{\psi}(x)\rangle\\\
&+g\tau_{R}P(x)\langle\delta\overrightarrow{\psi}(x),G\delta\overrightarrow{\psi}(x)\rangle,\end{split}$
(2i)

which implies that the empty solution becomes locally unstable whenever

$\frac{g\tau_{R}}{2}P(x)G\geq\Gamma(x),$ (2j)

that is when the particle gain induced by the CW source overcomes particle
losses.

We conclude this section with a final remark. Even though we have reported a
general formalism to account for $L$ photonic modes coupled to $L$ exciton
states, in what follows we will restrict our attention to the $L=2$ case.
Indeed, if on the one hand all the expressions previously reported hold true
in the general case, on the one other hand such case study is of particular
interest in light of the recent experiments [19, 20, 28], where it has been
observed that the relevant dynamics leading to the onset of condensation only
involves a pair of polariton bands below the exciton energy.

##  III Two photonic modes: spectral analysis

As discussed in Sec. II, in order to describe the onset of condensation we
first address the spectral properties of the system. In particular, since we
can safely assume the QW exciton energy-momentum dispersion to be flat, the
only non-trivial part concerns the representation of the counter-propagating
gapped Bloch resonances, in which the gap is induced by the periodic
dielectric modulation, e.g., as sketched of Fig. 1. Even though a complete
characterization of such modes in a one-dimensional (1D) periodic dielectric
structure can be performed through a Maxwell solver based, e.g., on a guided-
mode expansion [12, 11], here we rather follow the simplified approach
discussed in [35], where it was shown that close to normal incidence (i.e.,
around $k=0$), the dispersion of gapped photonic branches can be reproduced
and approximated by means of two counter-propagating modes with linear
dispersion that are diffractively coupled by an off-diagonal term $U$. The
latter can be either positive or negative, which is actually related to the
composition of the unit cell of the periodic 1D lattice (see, e.g., Ref.
[19]). In fact, we assume it as an effective model parameter, but we keep in
mind that it can be fully engineered through the characteristics of the
periodic pattern. We then consider the following matrices to build Eq. 2d

$E(x)=\hbar\left(\begin{array}[]{cccc}\omega_{A}-iv_{g}\partial_{x}&U&\Omega_{R}&0\\\
U&\omega_{A}+iv_{g}\partial_{x}&0&\Omega_{R}\\\ \Omega_{R}&0&\omega_{X}&0\\\
0&\Omega_{R}&0&\omega_{X}\\\ \end{array}\right),$ (2k)
$\Gamma(x)=\left(\begin{array}[]{cccc}\gamma_{A}&\tilde{\gamma}_{A}&0&0\\\
\tilde{\gamma}_{A}&\gamma_{A}&0&0\\\ 0&0&\gamma_{X}&0\\\ 0&0&0&\gamma_{X}\\\
\end{array}\right)\equiv\Gamma_{0},$ (2p)

and

$W(x)=V(x)\mathbb{1}_{4},$ (2u)

with $\mathbb{1}_{4}$ being the $4\times 4$ identity operator. Notice that we
hereby assume that the potential $V$ is the same for either photonic or
excitonic components, for simplicity of computations and without loss of
generality. In fact, while in practical cases one might differentiate between
photonic and excitonic contributions, the qualitative behavior of results we
are going to discuss do not depend from the choice made here. In this
scenario, the wavefunction $\overrightarrow{\psi}(x,\,t)$ reduces to the
following four-components vector:

$\overrightarrow{\psi}(x,\,t)=\left(A_{+}(x,t),A_{-}(x,t),X_{+}(x,t),X_{-}(x,t)\right)^{T}$
(2v)

This model describes two photonic modes separately coupled to two degenerate
excitonic states, in the presence of an external potential $V(x)$ that is
directly related to the presence of the CW source; $A_{+}$ and $A_{-}$
describe counter-propagating bands with linear dispersion whose slope is
directly defined by their group-velocity, $v_{g}$. These modes cross at $k=0$
with energy $\hbar\omega_{A}$, and they have an intrinsic loss rate given by
$\gamma_{A}$. The two modes $X_{+}$ and $X_{-}$ are then associated to a flat
QW exciton resonance at energy $\hbar\omega_{X}$, whose bare loss rate is
$\gamma_{X}$. The photonic modes are assumed to be linearly coupled via a term
proportional to $\hbar U-i\hbar\tilde{\gamma}_{A}$, which converts right ($+$)
propagating photons into left ($-$) propagating ones, and viceversa. In
particular, in order to ensure that $\Gamma_{0}$ describes a proper decay
term, one must require that $\gamma_{A}\geq|\tilde{\gamma}_{A}|$ (see Eq.
(2e)). Each photonic component is assumed to be independently coupled to an
exciton mode, with a coupling rate set by the Rabi frequency
$\Omega_{R}/(2\pi)$.

As it is shown in the following, such a minimal model is already able to
account for the onset of polariton condensation into discrete levels appearing
inside the energy-gap between polariton bands with opposite effective-mass. In
close agreement with experimental results reported in the literature, the
properties of such levels is shown to depend on both the bare polariton
dispersion and the pump characteristics. We first pay attention to the
spectral properties of the non-Hermitian operator $\tilde{H}(x)$. For the sake
of clarity, we proceed in two steps. In Sec. III.1, we first characterise the
polariton dispersion in the absence of the external potential $V(x)$. The role
of such space-dependent term in the appearence of a set of discrete levels is
addressed in Sec. III.2.

###  III.1 Eigenmodes: complex Polariton dispersion for $W(x)=0$

The complex eigenmodes of the exciton-photon coupled system can be obtained as
real and imaginary polariton bands by diagonalization of the operator

$\tilde{H}_{0}(x)=-i\frac{E(x)}{\hbar}-\Gamma_{0},$ (2w)

where $E(x)$ and $\Gamma_{0}$ are given in Eq.s (2k) and (2p), respectively.
To this purpose, let us consider a generic four-component wavefunction
$\overrightarrow{\psi}(x,t)$, as the one in Eq. 2v, and rewrite it in terms of
its Fourier components with respect to the real space coordinate $x$. By
applying the operator (2w) formally gives

$\tilde{H}_{0}(x)\overrightarrow{\psi}(x,\,t)=\int\frac{\mbox{d}k}{2\pi}e^{ikx}\tilde{H}_{0}(k)\overrightarrow{\psi}(k,\,t),$
(2x)

with $\tilde{H}_{0}(k)\equiv-i\frac{E_{0}(k)}{\hbar}-\Gamma_{0}$, where

$E_{0}(k)=\hbar\left(\begin{array}[]{cccc}\omega_{A}+v_{g}\,k&U&\Omega_{R}&0\\\
U&\omega_{A}-v_{g}\,k&0&\Omega_{R}\\\ \Omega_{R}&0&\omega_{X}&0\\\
0&\Omega_{R}&0&\omega_{X}\\\ \end{array}\right).$ (2y)

Since $\tilde{H}_{0}(k)$ satisfies the conditions
$\tilde{H}^{\dagger}_{0}(k)\tilde{H}_{0}(k)=\tilde{H}_{0}(k)\tilde{H}^{\dagger}_{0}(k)$,
by spectral theorem it can be diagonalized by means of a unitary operator,
$U_{k}$. In the present case the diagonalization can be performed
analytically. The resulting polariton eigenmodes (complex eigenvalues) are
explicitly reported in App. A.

Figure 2: Energy-momentum dispersion of the four polariton complex eigenmodes
obtained from the model in Eq. (2w). First, the real part of the polariton
dispersion ($Real[\lambda]$) is shown for (a) $\hbar U=4.45$ meV and (b)
$\hbar U=-4.45$ meV, respectively. The negative of the corresponding imaginary
part ($-Imag[\lambda]$) is shown for (c) the positive and (d) negative $U$
cases of (a) and (b), respectively. In all the panels, solid lines correspond
to the analytical expressions reported in Eqs. (2bf) and (2bg). Points labeled
with “num. diag.” correspond to the results of the numerical diagonalization.
Dashed lines in panels (a), (b),(c) and (d) describe the properties of the
bare photonic dispersions, i.e., modes $A_{+}(k)$ and $A_{-}(k)$. Dot-dashed
lines describe the two degenerate exciton bands. The other relevant parameters
are set as follows: $\hbar\omega_{X}=1527.5$ meV; $\hbar\gamma_{X}=0.2$ meV;
$\hbar\omega_{A}=1530$ meV; $\hbar\gamma_{A}=\hbar\tilde{\gamma}_{A}=0.6$ meV;
$v_{g}=110$ $\mu$m/ps and $\hbar=0.66$ meV$\cdot$ps.

Numerical diagonalization leads to the results shown in Fig. 2, where perfect
matching is obtained for the eigenvalues of $i\hbar\tilde{H}_{0}(k)$ obtained
numerically and the analytical expressions $\lambda_{\alpha,\,\beta}(k)$ given
in Eqs. (2bf) and (2bg). In particular, the parameters used in the present
case are compatible with typical values reported for inorganic semiconductor
samples [20]. As expected, at large $|k|$ the polariton eigenvalues
asymptotically approach the bare photonic and excitonic ones, both in terms of
real and imaginary parts. Also, and at difference with microcavity polariton
dispersions [1], this polariton band strucutre displays a number of very
peculiar features close to exciton-photon resonance. Independently of the sign
of $U$, the bands $\lambda_{+,-}$ and $\lambda_{+,+}$ are characterized by an
effective negative-mass $m^{*}$ (where we define
$m^{*}=\partial_{k}^{2}\lambda(k)|_{k=0}<0$), while $\lambda_{-,+}$ and
$\lambda_{-,-}$ are both characterized by a positive one, as evident from Fig.
2(a,b). In both cases, bands having opposite effective mass are separated by
an energy gap. In addition, in the vicinity of $k=0$ the bands are
characterized by a $k$-dependent imaginary part (i.e., loss rate). As shown in
Fig. 2(c), for $U>0$ intrinsic losses of negative-mass polaritons are smaller
than the exciton one, while positive-mass polariton loss rates are clearly
larger than $\hbar\gamma_{X}$. As shown in Fig. 2(d), the scenario is reversed
for negative values of $U$, where we see that positive-mass polaritons are
characterized by small losses, when compared to the exciton states. Such a
behavior is compatible with the one reported in Ref. [35], where a band
inversion between a bright and a dark mode is observed upon inverting the sign
of the diffractive coupling rate, $U$.

###  III.2 Gap-confined polariton states: $W(x)\neq 0$, $P(x)=0$

In this section we show numerical results obtained when an external potential
is added to the model addressed in the previous section. In particular, here
we focus on the energy region below the exciton resonance. Indeed, even though
our model prescribes the existence of polariton bands above the exciton energy
(upper polariton branches), we focus on experiments performed by analyzing
dynamics of states below $\hbar\omega_{X}$, which are the lowest energy
excitations and the ones towards which relaxation naturally occurs. The
external potential we used is given by a standard Gaussian function centered
at $x=0$

$V(x)=V_{0}\exp\\{-x^{2}/(2\sigma^{2})\\}\,,$ (2ad)

where $V_{0}\geq 0$ and $\sigma$ denote the height and the standard deviation
of the potential. The presence of such type of repulsive barrier can be
interpreted as an extra term accounting for local changes in the refractive
index for the photonic component, as well as the effects of a local blueshift
of the excitonic states. In experiments, both such effects are related to the
presence of the external input source $P(x)$. Here, we consider the effects of
such a barrier in order to provide a clearer interpretation of the results
shown in the next section, with the aim of characterizing the structure and
physical properties of eigenmodes.

Due to the local nature of $V(x)$, states with different $k$ get mixed and the
band structure gets modified. In particular, our analysis shows that discrete
levels appear within the energy gap between $\lambda_{+,-}(0)$ and
$\lambda_{-,-}(0)$. The properties of such states have been determined by
solving numerically the following eigenvalue problem:

$\tilde{H}(x)\vec{\psi}_{n}(x)=-i\frac{\,\lambda_{n}}{\hbar}\vec{\psi}_{n}(x),\quad\lambda_{n}\equiv
E_{n}-i\hbar\gamma_{n}$ (2ae)

where $E_{n}\equiv Re[\lambda_{n}]$ and $\hbar\gamma_{n}\equiv-
Im[\lambda_{n}]$ denote the real and imaginary parts of energy eigenvalues
corresponding to the eigenvectors $\vec{\psi}_{n}(x)$, respectively.

Figure 3: Spatial behavior of real ($Re[\phi]$) and imaginary ($Im[\phi]$)
parts of the 4 components of $\vec{\psi}_{0}(x)$, for a Gaussian potential
with $V_{0}=5$ meV and $\sigma=10$ $\mu$m. Panels (a) and (c) describe data
for $\hbar U=4.45$ meV, while panels (b) and (d) correspond to $\hbar U=-4.45$
meV. The star points labeled as “theory” are a plot of Eq. (2af) obtained from
the numerical solution for $A_{0,\,+}(x)$ (a similar agreement is observed for
the other component, not shown). Numerical results are obtained for
$\hbar\gamma_{A}=\hbar\tilde{\gamma}_{A}=\hbar\gamma_{X}=0.1$ meV. The other
relevant parameters are set as in Fig. 2.  Figure 4: Spatial behavior of real
($Re[\phi]$) and imaginary ($Im[\phi]$) parts of the 4 components of
$\vec{\psi}_{1}(x)$, calculated for the very same parameters assumed for the
results of Fig. 3. Again, panels (a) and (c) describe data for $\hbar U=4.45$
meV, while panels (b) and (d) correspond to $\hbar U=-4.45$ meV. The star
points labeled as “theory” are a plot of Eq. (2af) obtained from the numerical
solution for $A_{1,\,+}(x)$, with a similar agreement obtained for the other
component (not shown).

As an example, we report in Figs. 3 and 4 the numerical results of Eq. (2ae)
for the first two eigenmode components of $\vec{\psi}_{n}(x)$ (i.e., for $n=0$
and $n=1$, respectively). As a benchmark for the numerical solution, we notice
that

$X_{n,\,\pm}(x)=\frac{-i\Omega_{R}\,A_{n,\,\pm}(x)}{-i\frac{\lambda_{n}}{\hbar}+i\left(\omega_{X}+i\gamma_{X}+\frac{V(x)}{\hbar}\right)}\,,$
(2af)

as it is derived by direct inspection of the eigenvalue problem in Eq. (2ae).
Hence, we also plot in Figs. 3 and 4 the behavior of the $X_{n,\,+}(x)$
component obtained by plugging the numerical solution corresponding to
$A_{n,\,+}(x)$ into Eq. (2af) (markers labeled as “theory”), which perfectly
match the numerically computed solutions for the same component.

In addition, our numerical analysis also suggests that the photonic components
$A_{n,\,+}(x)$ and $A_{n,\,-}(x)$ are connected by the inversion operation,
that is $x\to-x$ (see App. B for details). In particular, a different symmetry
is obtained for positive and negative values of $U$. For $U>0$ data obtained
for $\overrightarrow{\psi}_{n}$ suggest that the photonic components are
related by the following expressions

$A_{n,\,+}(x)=(-1)^{n+1}A_{n,\,-}(-x),$ (2ag)

while for $U<0$ the relation becomes

$A_{n,\,+}(x)=(-1)^{n}A_{n,\,-}(-x)\,.$ (2ah)

Interestingly, given the expression reported in Eqs. (2af), (2ag) and (2ah),
it is easy to verify that

$\int
dx\langle\overrightarrow{\psi}_{n}(x),\overrightarrow{\psi}_{n+(2l+1)}(x)\rangle=0$
(2ai)

for $n$ and $l$ integers, i.e., any two eigenstates
$\overrightarrow{\psi}_{n}$ and $\overrightarrow{\psi}_{m}$ are orthogonal
whenever $n$ and $m$ have an opposite parity.

We proceed by showing the behavior of the eigenenergies corresponding to
$\vec{\psi}_{n}(x)$ within the lower polariton gap as a function of the
parameter $V_{0}$ and for different values of $\sigma$. The numerical results
are shown in Fig. 5. Numerical data reported in Fig. 5(a) describe the
behavior for $\hbar U=4.45$ meV, while in Fig. 5(b) are reported data
corresponing to $\hbar U=-4.45$ meV. Similarly to what observed for
$Real[\lambda]$, the energy of such discrete levels does not depend on the
sign of $U$. So for the sake of simplicity, let us pay attention to Fig. 5(a).
In analogy to the behavior expected when considering a confining potential
$\tilde{V}(x)=-V(x)$ perturbing the motion of a positive-effective mass
($m^{*}$) quantum particle, that is

$-\frac{\hbar^{2}}{2m^{*}}\frac{d^{2}}{dx^{2}}\phi_{n}(x)+\tilde{V}(x)\phi_{n}(x)=E_{n}\phi_{n}(x),\quad
m^{*}>0$ (2aj)

the number of bound-states as well as the distance with respect to the minimum
of the parabolic band $E(k)={\hbar^{2}k^{2}}/{(2m^{*})}$, only depend on the
properties of the potential well $\tilde{V}(x)$, namely its width and depth.
Obviously, the wider the well, the larger the number of bound states supported
by the well. A similar behavior is observed in the present case, where
negative effective mass excitations are trapped within a potential barrier.
For $\sigma=5$ $\mu$m, the potential supports a single discrete state (i.e.,
with $n=0$). For $\sigma=10$ $\mu$m, in the range of $V_{0}$ values
considered, the system supports two discrete levels, the second entering in
the gap in the vicinity of $V_{0}\approx 4$ meV. When the width of the
repulsive barrier is increased further, also the number of discrete states
increases. This behavior is confirmed by the results for $\sigma=15$ $\mu$m,
where the second and third discrete levels supported by the repulsive
potential appear in the vicinity of $V_{0}\approx 2$ meV and $V_{0}\approx
5.5$ meV. The first main difference with respect to the solutions of Eq. (2aj)
concerns the range of energies spanned by such bound states while sweeping
$V_{0}$. In the standard positive mass case, the deeper the potential, the
smaller the quantization energy. In the present case, as it is obtained in
Fig. 5, the eigenenergies of such states are always bounded between a lower
and an upper value. Such limits correspond to the maximum of the negative-mass
band $\lambda_{+,\,-}$ and to the minimum of the positive-mass band
$\lambda_{-,\,-}$, respectively. In particular, when a state reaches and
crosses the minimum of the upper band at $\lambda_{-,\,-}$, it can no longer
be normalized and it represents an extended solution across the whole real
space domain.

Figure 5: Real part of the discrete levels eigenenergies, $E_{n}$ (in meV)
created by the repulsive potential within the polariton gap below the exciton
resonance, as a function of the potential height $V_{0}$ and for
$\sigma=5,\,10,\,15$ $\mu$m (see legend), respectively. Results are shown for
(a) $U>0$ and (b) $U<0$, respectively. The relevant model parameters are set
as in Fig. 2. In both panels, the dot-dashed horizontal line below 1520 meV,
and the dashed horizontal line above 1523 meV correspond to the band extrema
of the eigenmodes $\lambda_{+,\,-}(k)$ (maximum at $k=0$) and
$\lambda_{-,\,-}(k)$ (minimum), respectively (see also Fig. 2(a)).  Figure 6:
Imaginary parts of discrete eigenmodes within the energy gap,
$\hbar\gamma_{n}$ (in meV) plotted as a function of $V_{0}$ for $\sigma=10$
$\mu$m. Results for $\hbar U=4.45$ meV are shown for (a)
$\hbar\gamma_{A}=\hbar\tilde{\gamma}_{A}=0.1$ meV, (c) $0.15$ meV, and (e)
$0.3$ meV, respectively. Results for $\hbar U=-4.45$ meV are shown for (b)
$\hbar\gamma_{A}=\hbar\tilde{\gamma}_{A}=0.1$ meV, (d) $0.15$ meV, and (f)
$0.3$ meV, respectively. In these simulations we have assumed
$\hbar\gamma_{X}=0.1$ meV for the exciton modes. All the other parameters are
set as in Fig. 2. In all the panels: the dot-dashed horizontal line
corresponds to the imaginary part of the $\lambda_{+,-}(k=0)$ eigenmode, the
dashed horizontal line corresponds to the one of $\lambda_{-,-}(k=0)$, the
solid horizontal line is placed at the value $\hbar\gamma_{A}$.

While the resonant energy (real part of the eigenvalues) within the energy gap
is not affected by the sign of $U$, the losses (i.e., the imaginary part of
the eigenvalues) of such discrete levels strongly depend on the nature of the
negative mass branch giving out of which they arise on increasing $V_{0}$. In
fact, numerical results for either $U>0$ or $U<0$ are reported in Fig. 6, and
for different values of $\hbar\gamma_{A}$. In particular, results are seen to
be quite different depending on the sign of $U$, displaying low qualitative
dependence on the value of $\gamma_{A}$. We thus focus on commenting the
differences between Figs. 6(a) and 6(b). When $V_{0}$ is slightly larger than
zero, the first mode entering in the gap ($n=0$) is characterized by an
imaginary part $\hbar\gamma_{n}$ that is close to the one corresponding to the
polariton band $\lambda_{+,-}(k=0)$ (dot-dashed horizontal line in both
panels). On the increasing of $V_{0}$, deviations from this value are
observed. The latter depends on the sign of $U$ (see also Fig. 2). Hence, when
$U>0$ an increasing of the repulsive barrier height yields an increased
$\gamma_{n}$, as seen in Fig. 6(a). Conversely, when the system is
characterized by a negative diffractive coupling $U$, $\gamma_{n}$ decreases
on increasing $V_{0}$, as shown in Fig. 6(b). In particular, similarly to what
already shown for the real part of the eigenvalues, also the imaginary parts
of the discrete gap-confined levels are bound to vary between those of the
negative- and positive-mass bands at $k=0$, i.e., $\lambda_{+,-}(k=0)$ and
$\lambda_{-,-}(k=0)$, respectively. In summary: the discrete levels created by
the repulsive potential $V(x)$ are confined within the polariton energy-gap,
and their complex eigenvalues are such that real and imaginary parts vary
continuously between those of $\lambda_{+,-}(k=0)$ and $\lambda_{-,-}(k=0)$.

##  IV Two-mode case: polariton condensation in the CW regime

As it is known from previous studies, polariton condensation is identified
from accumulation of excitations in one specific state, usually the lowest
energy one in planar microcavities, and this phenomenology is well captured by
a one-dimensional single-mode GPE formulation. In this Section we show
numerical results concerning the onset of polariton condensation in the
steady-state regime, when the coupled multi-band exciton-photon system is
connected to a particle reservoir driven by a spatially-dependent CW input
source. The expression of such driving term is given as

$P(x)=P_{0}\exp\\{-x^{2}/(2\sigma^{2})\\}\,,$ (2ak)

in which $P_{0}$ is a time-independent pump-rate per $\mu$m, i.e., $[P_{0}]=$
($\mu$m$\cdot$ps)-1. In particular, following the discussion from the previous
Section, it is hereby assumed that the pumping rate is linearly related to the
barrier height, $V_{0}$ by the following relation:

$P_{0}=\frac{\eta V_{0}}{\hbar}\,,$ (2al)

in which $\eta$ is a parameter having the dimension $[\eta]=\mu$m-1.
Concerning the reservoir-polariton interactions, we assume the two excitonic
modes to be similarly coupled to the reservoir density, and we made the same
assumption for the photonic components for simplicity of computation and
without loss of generality. In this scenario, the operator $G$ is diagonal
and, due to the constraints reported in Eqs. (2cb) and (2cc), it has the
following structure:

$G=\left(\begin{array}[]{cccc}(1-\alpha)/2&0&0&0\\\ 0&(1-\alpha)/2&0&0\\\
0&0&\alpha/2&0\\\ 0&0&0&\alpha/2\\\ \end{array}\right)$ (2am)

with $\alpha$ being an effective parameter such that $0\leq\alpha\leq 1$.

Since we are interested in the stationary solutions of Eq. (2) at large time
$t$, and since the model always admits as a possible solution
$\vec{\psi}(x)=0$, we consider an initial configuration characterized by a
small population. By doing so, we are able to understand whether a tiny
perturbation added to the “empty” solution actually leads to the appearance of
a stable, macroscopically populated condensate. Details about the system
initialization are given in App. C. Such initial configuration is then evolved
in time by means of a standard numerical integration algorithm (i.e., explicit
Runge-Kutta 4(5) method), and the convergence towards a possibly nonempty
condensate configuration is monitored by looking at the time behavior of the
total population, that is

$\begin{split}&N_{\psi}(t)=\int
dx\langle\vec{\psi}(x,\,t),\vec{\psi}(x,\,t)\rangle=\\\ &=\int
dx\left[\,|A_{+}|^{2}+|A_{-}|^{2}+|X_{+}|^{2}+|X_{-}|^{2}\right],\end{split}$
(2ar)

where $A_{\sigma}\equiv A_{\sigma}(x,\,t)$ and $X_{\sigma}\equiv
X_{\sigma}(x,\,t)$.

Figure 7: Behavior of the condensate occupation $N_{\psi}(t)$ (Eq. (2ar)) as a
function of time (in ps) for different values of $V_{0}$ (top legend), for (a)
$\hbar U=4.45$ meV and (b) $\hbar U=-4.45$ meV. Results are reported for
$\hbar\gamma_{A}=\hbar\tilde{\gamma}_{A}=0.1$ meV, $\hbar\gamma_{X}=0.1$ meV,
$\alpha=0.01$ (Eq. (2am)), $\eta=2$ $\mu$m-1 and $g=0.1$ $\mu$m/ps. All the
other relevant parameters are set as in Fig. 2.

An example of time evolution is reported in Fig. 7, where the behavior of
$N_{\psi}(t)$ as a function of time is again compared for the cases $U>0$ and
$U<0$, assuming the same CW driving protocol and the same values for $V_{0}$.
Even if the initial population injected in the light-matter subsystem is small
(see App. C), independently of the sign of $U$, for sufficiently large values
of $V_{0}$ the system dynamics tends towards a stable steady-state
configuration displaying a macroscopic occupation, i.e.,

$\lim_{t\to\infty}N_{\psi}(t)=N_{\psi,\,ss}\gg 0.$ (2as)

Interestingly, according to the results shown in Fig. 7, the convergence
towards a nonempty condensate seems to depend on the sign of $U$. Indeed, if
on the one hand for $V_{0}>1$ meV in both cases the system approaches a
configuration with $N_{\psi,\,ss}>0$, on the other hand at exactly $V_{0}=1$
meV polariton condensation is only triggered for $U>0$. In particular, since
$V_{0}$ is in one-to-one correspondence with $P_{0}$ (Eq. (2al)), such
behavior suggests that in systems with $U>0$ condensation would occur for
lower values of $P_{0}$, i.e., they would display a lower condensation
threshold. An interesting prediction to be verified experimentally.

Figure 8: Steady-state occupation of the exciton-photon coupled system,
$N_{\psi\,ss}$, as a function of $V_{0}$ (in meV), for different values of
$\hbar\gamma_{A}$ (top legend). Triangle markers correspond to data obtained
for $\hbar U=4.45$ meV, while square points are for $\hbar U=-4.45$ meV. All
the other relevant parameters are set as in Fig. 7.

In order to understand whether this is indeed the case, we performed a
characterization of the dependence of the steady-state occupation
$N_{\psi,\,ss}$ on $V_{0}$. Results obtained for different values of the
photonic losses are shown in Fig. 8. As in the previous cases, we compare the
trend in data obtained for $U>0$ (triangles) and $U<0$ (squares),
respectively. For what concerns the condensation threshold, data corresponding
to $U>0$ do not display a significant dependence on $\hbar\gamma_{A}$. In
particular, for all the cases considered the critical value of $V_{0}$ leading
to condensation is slightly below 1 meV. On the contrary, at first glance,
numerical results for $U<0$ display a considerably different behavior. In
particular, the condensation threshold in this case shows more pronounced
dependence on the photonic losses, going from $V_{0}\simeq 1$ meV to $\simeq
3$ meV while increasing $\hbar\gamma_{A}$ from $0.1$ to $0.3$ meV.

In addition, data obtained for values of $U$ having a different sign show a
quantitatively different dependence on $V_{0}$. If on the one hand they
clearly show that condensation is occurring in both cases, on the other hand
the value of $U$ does not only affect the condensation threshold. In the next
two sections we show that all these peculiar features of the model can be
interpreted by looking at the properties of the states previously
characterized in Sec. III.2. In particular, we first address the behavior of
the losses of the condensed states. Then, we analyse the spectral properties
of the radiation emitted from the condensate, and we show that it is peaked at
energies corresponding to the discrete levels considered in the previous
section.

###  IV.1 Analysis of the condensate losses

In this section, by analyzing the condensate losses, we show that the steady-
state configuration is given by one or few of the discrete eigenmodes
characterized in the previous section. However, in order to keep the
discussion as clear as possible and to show the main idea, we first assume
that the configuration approached by the system during relaxation towards its
lower energy eigenstates corresponds to a single normalized eigenfunction
$\vec{\psi}_{n}(x)$ of a gap-confined state, as obtained from the operator
$\tilde{H}(x)$ characterized in Sec. III.2. In this case, the steady state of
the system is given by

$\lim_{t\to+\infty}\vec{\psi}(x,\,t)=\sqrt{N_{\psi,ss}}\vec{\psi}_{n}(x)\,.$
(2at)

Then, let us consider the spatial integral of Eq. 2f, with the steady-state
condition imposing that the derivative of the total particle density should
become zero. Eq. 2f then reduces to

$\begin{split}0&=\int P(x)\,dx-\frac{1}{\tau_{R}}\int
n_{ss}(x)\,dx-2\gamma_{n}N_{\psi,\,ss}=\\\ &=\int
P(x)\,dx-\frac{1}{\tau_{R}}N_{R,\,ss}-2\gamma_{n}N_{\psi,\,ss}\,,\end{split}$
(2au)

in which $N_{R,\,ss}$ denotes the total number of particles accumulated in the
reservoir. Since $N_{\psi,\,ss}>0$, we can rearrange the terms in Eq. (2au)
and get an expression for the imaginary part of the n-th mode steady state,
$\gamma_{n}$, which is

$\gamma_{n}=\frac{\int P(x)\,dx-\gamma_{R}\,N_{R,\,ss}}{2\,N_{\psi,\,ss}}\,,$
(2av)

in which $\gamma_{R}=1/\tau_{R}$. By plugging the expressions reported in Eqs.
(2ak) and (2al) into Eq. (2av), we finally obtain an analytic expression for
the condensate losses as

$\hbar\gamma_{n}=\frac{\eta
V_{0}\sqrt{2\pi}\sigma-\hbar\gamma_{R}\,N_{R,\,ss}}{2\,N_{\psi,\,ss}}\,.$
(2aw)

In general, depending on the values of $V_{0}$ and $\sigma$ in the expression
for the repulsive potential, the stationary configuration might correspond to
a superposition of more than one gap-confined eigenstate. For instance, as
shown in Fig. 5, when $\sigma=10$ $\mu$m and $V_{0}\in[0,\,7]$ meV, two
discrete levels appear within the polariton gap, whose eigenfunctions are
$\psi_{0}(x)$ and $\psi_{1}(x)$. Then, let us suppose that

$\vec{\psi}_{ss}(x)=\alpha_{0}\vec{\psi}_{0}(x)+\alpha_{1}\vec{\psi}_{1}(x)\,.$
(2ax)

Since the two eigenstates are orthogonal by symmetry constraints, see Eq.
(2ai), by following the same steps leading to Eq. (2aw) we can write the
overall steady state condensate loss in terms of the weighted imaginary parts
of the two discrete modes that are present in the gap as

$\hbar\bar{\gamma}=\frac{\hbar\gamma_{0}|\alpha_{0}|^{2}+\hbar\gamma_{1}|\alpha_{1}|^{2}}{N_{\psi,\,ss}}=\frac{\eta
V_{0}\sqrt{2\pi}\sigma-\hbar\gamma_{R}\,N_{R,\,ss}}{2\,N_{\psi,\,ss}}\,,$
(2ay)

in which $N_{\psi,\,ss}=|\alpha_{0}|^{2}+|\alpha_{1}|^{2}$.

Figure 9: Comparison between the imaginary parts of the discrete eigenmodes
induced by the potential $V(x)$ ($\hbar\gamma_{n}$) and the expected
condensate losses obtained from Eq. (2ay), (star-shaped markers), plotted as a
function of $V_{0}$. Results are shown for $U>0$ and (a)
$\hbar\gamma_{A}=\hbar\tilde{\gamma}_{A}=0.1$ meV, (c) 0.15 meV, and (e) 0.3
meV, respectively (solid horizontal line). Similarly, results are shown for
$U<0$ and (b) $\hbar\gamma_{A}=\hbar\tilde{\gamma}_{A}=0.1$ meV, (d) 0.15 meV,
and (f) 0.3 meV, respectively (solid horizontal line). The other parameters
are set as in Fig. 6. The dashed and dot-dashed horizontal lines correspond to
the imaginary parts of the two gap-delimiting eigenmodes, $\lambda_{-,-}(k=0)$
and $\lambda_{+,-}(k=0)$.

A direct comparison between Eq. (2ay) and the numerical results previously
reported in Fig. 6 is given in Fig. 9. For $U<0$, as shown in Figs. 9(b,d,f)
for different values of $\gamma_{A}=\tilde{\gamma}_{A}$, the steady-state
condensate loss obtained from of Eq. (2ay) (star-shaped markers) nicely follow
the numerical solutions for the imaginary parts of the first discrete level
appearing within the polariton band gap on increasing $V_{0}$, for all the
values of $V_{0}$ such that $N_{\psi,\,ss}>0$ (see, e.g., Fig. 8). In
particular, in this case the imaginary part of eigenmodes decreases when
increasing $V_{0}$, and we clearly notice that our protocol leads to the
stabilization of a configuration compatible with the macroscopic occupation of
the first mode confined in the gap. A similar behavior is displayed for $U>0$,
results shown in Figs. 9(a,c,e). However, by further increasing the pumping
rate (i.e., $V_{0}$ in our model) the star-shaped markers move along a path
going smoothly from the branch corresponding to the first gap-confined
eigenmode imaginary part to the one corresponding to the second discrete
eigenmode. Such behavior is somehow compatible with the fact that the first
discrete state of the potential is moving towards the upper polariton branch,
$\lambda_{-,-}(k)$, and it is getting out from the energy gap (see also Fig.
5), such that also the second mode within the gap becomes populated, as it
will be detailed in the following Section.

###  IV.2 Radiation emission from the condensate

In this section we show numerical results describing the energy-momentum
spectral density of the the emitted radiation from the polariton system. In
the hypothesis that the field emitted by the structure is proportional to the
field inside the system (as it is done in a standard input-output scenario),
the emitted radiation has an overall space-time profile $A(x,\,t)$ given by

$A(x,\,t)\propto A_{+}(x,\,t)+A_{-}(x,\,t)\,.$ (2az)

Then the energy-momentum characteristics of the emitted radiation are
essentially encoded into the Fourier transform of $A(x,\,t)$

$A(k,\,E)=\int\frac{dk}{2\pi}\int\frac{dE}{2\pi\hbar}e^{-ikx}e^{iEt/\hbar}A(x,\,t).$
(2ba)

The corresponding spectral density is given as

$I(k,\,E)=|A(k,\,E)|^{2}\,.$ (2bb)

In order to display results related to the dependence of the spectral density
on $V_{0}$, we plot in Fig. 10 the normalized quantity defined as

$\tilde{I}(k,\,E)=I(k,\,E)/\mbox{max}[I(k,\,E)]\,.$ (2bc)

Figure 10: Direct comparison between the energy of the discrete levels created
by the potential $V(x)$ within the polariton gap (real part of eigenvalues,
coinciding for either positive or negative $U$, as shown in panels a, d, g)
calculated for $\sigma=10$ $\mu$m and
$\hbar\gamma_{A}=\hbar\tilde{\gamma}_{A}=0.3$ meV. The corresponding
(normalized) spectral density obtained from Eq. 10 is shown for $\hbar U=4.45$
meV and corresponding to the three values of $V_{0}$ highlighted with the
vertical red lines in panels (a,d,g), respectively, in panels (b) $V_{0}=0$,
(e) $V_{0}=3.0$ meV, and (h) $V_{0}=6.25$ meV. Similarly, results are shown
for $\hbar U=-4.45$ meV in panels (c), (f), and (i) in correspondence of the
same values of $V_{0}$. All the other relevant parameters are set as in Fig.
7.

We report results for the spectral density corresponding to positive $U$ and
different values of the pumping strength, quantified as $V_{0}=0$ in Fig.
10(b), $V_{0}=3.0$ meV in Fig. 10(e), and $V_{0}=6.25$ meV in Fig. 10(h). The
same quantity is also reported for the corresponding values of $V_{0}$, but
assuming $U<0$ in Figs 10(c,f,i). From these results, the close connection
between the discrete gap-confined eigenstates characterized in the previous
Sections to the position of the peaks of the spectral density should be quite
evident.

In all the cases considered, the system was initially prepared by injecting a
small population within in the field $\vec{\psi}$, below the exciton states
(see Eq. (2bj)). Then, for each value of the parameter $V_{0}$ and for the two
values $\hbar U=\pm 4.45$ meV, time evolution is solved until $t=100$ ps.
Finally, the spectral density distribution is obtained by taking the Fourier
transform along the entire space-time evolution by following the prescription
reported in Eqs. (2az), (2ba), and (2bb). For $P_{0}=V_{0}=0$, the system
evolves under the action of $\tilde{H}_{0}(x)$ (Eq. (2w)). Since no driving
source as well as repulsive barrier are included into the dynamics, the
spectral density displays a behavior compatible with the band structure
characterized in Sec. III.1 for both $U>0$ [Fig. 10(b)] and $U<0$ [Fig.
10(c)]. In particular, we notice the absence of emission at $k\simeq 0$ in the
lower polariton branch evidenced in Fig. 10(b) as well as in the upper
polariton branch in Fig. 10(c), which marks the behavior of a dark mode and it
is compatible with recent experimental results for the emission below
threshold [19, 20, 28].

For $V_{0}=3$ meV, the system is expected to have a single discrete level
within the gap. For the model employed so far, this gap-confined eigenmode is
occurs at a resonant energy $E_{n=0}\simeq 1521$ meV, both for $U>0$ and
$U<0$, as shown in Figs. 10(d,e,f). In fact, in this case the spectral density
is essentially nonzero only in the neighborhood of this eigenmode energy, and
the spectral density distribution displays a behavior as a function of the
wave vector that is, again, fully compatible with experimental results [19,
20, 28].

Finally, by further increasing the pump-rate, a new state enters in the gap,
as shown in Fig. 10(g) in which the vertical line at $V_{0}=6.25$ meV crosses
both branches describing the energy dependence of the first and the second
gap-confined modes supported by the structure. In agreement with the
interpretation provided in the previous Section for $U>0$, when
$\hbar\bar{\gamma}$ assumes an intermediate value between the first and second
mode within the gap [Fig. 9(e)], the spectral density is peaked in
correspondence of the energy of both states, as shown in Fig. 10(h). This
behavior is observe also experiments (see, e.g., Fig.2 of “Extended data
figures and table” in Ref. [19]). On the other hand, for $U<0$ we do see that
the emission profile is peaked at the energy of the first discrete state
within the gap, while very little population is coming from the second-order
mode. This result is compatible with the data displayed in panel Fig. 9(f).

Figure 11: Behavior of the spectral density $\tilde{I}_{n}(k)$ as a function
of $k$ (in $\mu$m-1), computed from the numerical solution for
$\vec{\psi}_{0}(x)$ and (a) $\hbar U=4.45$ meV and (b) $\hbar U=-4.45$ meV, as
well as from $\vec{\psi}_{1}(x)$ for (c) $\hbar U=4.45$ meV and (d) $\hbar
U=-4.45$ meV. Data reported in all panels have been obtained for a repulsive
barrier $V(x)$ with $V_{0}=6.25$ meV and $\sigma=10$ $\mu$m. All the other
parameters are set as in Fig. 10.

For completeness and clarity, we hereby report also the normalized spectral
density associated to the bound-state solutions of the effective model,
$\tilde{H}(x)$, for the same repulsive barrier used to obtain the out-of-
equilibrium results reported in Figs. 10(h) and (i), which is shown in Fig.
11. In particular, in each panel of Fig. 11 we report the behavior of

$\tilde{I}_{n}(k)=\frac{I_{n}(k)}{\mathrm{max}\\{I_{n}(k)\\}},\quad
I_{n}(k)=\left|A_{n,\,+}(k)+A_{n,\,-}(k)\right|^{2},$ (2bd)

as a function of $k$. By direct comparison of Figs. 11(a) and (c) to the
numerical results shown in panel Fig. 10(h), we notice that there is a very
good agreement between the spectral properties displayed by the output field
for the two bound states, $\vec{\psi}_{0}(x)$ and $\vec{\psi}_{1}(x)$, and
those displayed by the radiation emitted from the condensate. These results
are also in remarkable good agreement with experimental outcomes obtained on
similar devices [19, 20, 28]. A similarly good agreement is displayed by the
data obtained for $\hbar U=-4.45$ meV, and the corresponding ones in Fig.
10(i). We notice that the plot in Fig. 11(d) corresponds to the second mode
($n=1$) supported by the potential, which is not displaying macroscopic
occupation in this case.

##  V Summary and conclusion

We have reported an extensive theoretical analysis of the polariton
condensation mechanism occurring in periodically patterned QW multilayers
under incoherent driving. From a theoretical point of view, this work has
required the generalization of the non-equilibrium Gross-Pitaevskii
formulation under incoherent reservoir driving, originally proposed in Ref.
[29], in order to account for multiple photonic bands interacting with the
same QW exciton degenerate modes, which is an original contribution of the
present work.

After analyzing the complex eigenmodes of the non-Hermitian Hamiltonian model
including exciton-photon coupling and losses, we have numerically solved the
driven-dissipative dynamics in the presence of a continuous wave Gaussian pump
spot. We have modelled the effects of such a driving field as an effective
potential barrier for the exciton-photon field, showing that it naturally
induces the confinement of negative mass polaritons arising from the lower
branch and becoming trapped due to the energy gap between the two polariton
bands at energies below the exciton resonance.

In addition, we have studied the emission characteristics of these polaritonic
branches, which are associated to the sign of the diffractive coupling term
introduced to mimic the effect of the periodic pattern on photonic eigenmodes.
In particular, the positive sign of such diffractive coupling term is
associated with the lower, negative mass branch being dark at normal incidence
(BIC condition) with minimal losses, while if it is negative this same branch
becomes bright at normal incidence and its imaginary part is much larger.
Analysing the behavior of the model as a function of the pump intensity
(effectively described by the height of this potential barrier), we have found
that condensation actually occurs in such gap-confined eigenmodes, which
become spatially quantized as a function of the depth of the confining
potential. We then make the relevant conclusion that condensation in these
systems is indeed the result of such gap-confinement of negative mass
polaritons, irrespective of their dark or bright emission at normal incidence.
In particular, these results are in remarkable agreement with recent
experimental findings corresponding to samples in which the diffractive
coupling term is positive [19, 20]. Hence, we predict that a similar
phenomenology is going to occur even in the opposite case of a negative
diffractive coupling term, which motivates new experiments.

As a continuation of the present work, we envision applications of this
theoretical framework to more complex pumping configurations, such as, e.g.,
periodic repetitions of a finite number of pumping spots, as recently reported
in experiments [28]. In addition, a careful benchmarking of this effective
model with experimental data might be a further test of the usefulness of the
proposed approach. In particular, since experiments are often performed under
pulsed excitations, an analysis of the solutions of the present model for
pulsed driving, which is beyond the scope of the present manuscript, might be
one of the possible next steps to undertake as a future project. We also
notice that the model employed is a one-dimensional one, which essentially
captures the relaxation mechanism in an effective way. On the other hand, the
detailed dynamics along the transverse direction is completely neglected,
which might still play a role in experiments. Hence, a possible further
extension of the model might include the description of the modes dispersion
along the transverse spatial direction as well.

###### Acknowledgements.

We acknowledge useful discussions with L. C. Andreani, V. Ardizzone, A.
Gianfrate, E. Maggiolini, H. C. Nguyen, H. S. Nguyen, F. Riminucci, D.
Sanvitto, H. Sigurdsson, S. Zanotti. This research was financially supported
from the Italian Ministry of University and Scientific Research (MUR) through
PRIN-2017 project “Interacting Photons in Polariton Circuits” (INPhoPOL),
Grant No. 2017P9FJBS_001.

##  Appendix A Analytic expression of complex polariton dispeersion

As mentioned in Sec. III.1, the bands associated to the Hamiltonian model
$i\hbar\tilde{H}_{0}(k)$ (Eq. (2w)) can be determined analytically by looking
for the four complex roots of the characteristic polynomial equation
$P(\lambda)=\mbox{det}\left[E_{0}(k)-i\Gamma_{0}-\lambda\mathbb{1}_{4}\right]$,
in which “det” denotes the matrix determinant. In what follows, we assume that
$\sqrt{z}$ gives the square root a complex number $z$ with non-negative
imaginary part, and we define

$E_{X}=\hbar(\omega_{X}-i\gamma_{X}),\,E_{A}=\hbar(\omega_{A}-i\gamma_{A}),\,\tilde{U}=U-i\tilde{\gamma}_{A}\,.$
(2be)

Hence, the four eigenvalues of $i\hbar\tilde{H}_{0}(k)$ at each wave vector
$k$, i.e., the polariton bands, are given by the following analytic
expressions

$\begin{split}&\lambda_{+,\,\pm}(k)=\frac{E_{X}+E_{A}-\sqrt{(\hbar\tilde{U})^{2}+(\hbar
v_{g}k)^{2}}}{2}\\\
&\pm\frac{1}{2}\sqrt{\left(E_{X}-E_{A}+\sqrt{(\hbar\tilde{U})^{2}+(\hbar
v_{g}k)^{2}}\right)^{2}+4(\hbar\Omega_{R})^{2}},\end{split}$ (2bf)

and

$\begin{split}&\lambda_{-,\,\pm}(k)=\frac{E_{X}+E_{A}+\sqrt{(\hbar\tilde{U})^{2}+(\hbar
v_{g}k)^{2}}}{2}\\\
&\pm\frac{1}{2}\sqrt{\left(E_{X}-E_{A}-\sqrt{(\hbar\tilde{U})^{2}+(\hbar
v_{g}k)^{2}}\right)^{2}+4(\hbar\Omega_{R})^{2}}\,,\end{split}$ (2bg)

respectively, corresponding to a total of four complex branches. In Eqs. (2bf)
and (2bg), given the solution $\lambda_{\alpha,\,\beta}(k)$, the subscript
$\alpha$ controls the sign in front of the $\sqrt{(\hbar U)^{2}+(\hbar
v_{g}k)^{2}}$ term, while the subscript $\beta$ controls the sign in front of
the square root containing the Rabi term, $4(\hbar\Omega_{R})^{2}$.

Figure 12: Spatially dependent real part of the first four eigenfunctions
obtained from the model $\tilde{H}(x)$ for $V_{0}=2.5$ meV and $\sigma=35.0$
$\mu$m, assuming $\hbar U=+4.45$ meV and all the other relevant model
parameters set as in Fig. 3. Figure 13: Spatially dependent real part of the
first four eigenfunctions obtained from the model $\tilde{H}(x)$ for
$V_{0}=2.5$ meV and $\sigma=35.0$ $\mu$m, assuming $\hbar U=-4.45$ meV and all
the other relevant model parameters set as in Fig. 3.

##  Appendix B Symmetries and orthogonality of gap-confined states

We hereby report additional results concerning the inversion symmetry of the
left- and right-moving photonic components, i.e., Eqs. (2ag) and (2ah)
reported in the main text. In particular, here we report the real part of the
photonic components of the first four eigenstates supported by a repulsive
potential with $V_{0}=2.5$ meV and $\sigma=35$ $\mu$m, i.e.,
$\\{\vec{\psi}_{0},\,\vec{\psi}_{1},\,\vec{\psi}_{2},\,\vec{\psi}_{3}\\}$. We
first show the numerical results obtained for $\hbar U=+4.45$ meV in Fig. 12,
while results obtained for $\hbar U=-4.45\,meV$ are reported in Fig. 13.
Similar behaviors are obtained by considering the imaginary parts of the same
components (not shown). As it is possible to see by direct comparison between
the star-shaped markers and the dashed lines, for both positive and negative
values of $U$ the numerical solutions are fully consistent with the two
equations reported in Sec. III.2. In addition, due to the connection between
the excitonic and photonic components reported in Eq. (2af), since
$V(x)=V(-x)$, the exact same symmetry relations hold true also for
$X_{\pm,n}(x)$.

Figure 14: Absolute value of the overlap integral, $\mathcal{M}_{n,\,m}$,
calculated for the first four eigenfunctions supported by $\tilde{H}(x)$ when
$V_{0}=2.5$ meV and $\sigma=35.0$ $\mu$m, plotted in log scale (color bar on
the side). Data have been obtained for $\hbar U=4.45$ meV. The other relevant
parameters are set as in Fig. 3.

This is particularly relevant when considering the overlap integral between
different eigenfunctions, i.e., formally defined as

$\mathcal{M}_{n,\,m}\equiv\int\langle\vec{\psi}_{n}(x),\,\vec{\psi}_{m}(x)\rangle
dx\,.$ (2bh)

Indeed, under the action of the inversion operation $x\to-x$, one obtains that
the overlap integral transforms as follows

$\mathcal{M}_{n,\,m}\to(-1)^{m+n}\mathcal{M}_{n,\,m},$ (2bi)

which implies that $\mathcal{M}_{n,\,m}$ is identically zero whenever $n$ and
$m$ have opposite parity. It would be tempting to conclude that, as in the
Hermitian case, $\mathcal{M}_{n,\,m}=\delta_{n,m}$, which means that
eigenfunctions corresponding to different eigenvalues are orthogonal. However,
this does not seem be the case, as we have verified numerically. These results
are summarized in Fig. 14, where we show the overlap integral (2bh) between
the pairs of 4 eigenmodes shown in Fig. 12 is reported in log scale (see color
bar). While $\mathcal{M}_{n,\,m}$ is essentially zero, up to numerical
deviations, when $n$ and $m$ have opposite parity (in agreement with the
discussion above), when $n$ and $m$ are different but have the same parity the
overlap integral $\mathcal{M}_{n,\,m}\neq 0$. So in general we cannot conclude
that eigenfunctions corresponding to different eigenvalues are always
orthogonal.

##  Appendix C System initialization

In this section we provide some details concerning the system initialization
used to derive all the results shown in Sec. IV. In fact, in the time-
evolution simulations the reservoir has been considered as initially empty,
i.e., $n(x,t=0)=0$. For what concerns the exciton-photon subsystem, we
consider as initial configuration a state where a tiny population is injected
into the two bands $\lambda_{+,\,-}(k)$ and $\lambda_{-,\,-}(k)$. To do this,
we took advantage of the explicit expression of polariton eigenstates in
momentum space derived for $W(x)=0$. By denoting the eigenvector associated to
$\lambda_{\alpha,\,\beta}(k)$ as $\vec{w}_{\alpha,\,\beta}(k)$, a generic
system configuration where only the bands below the exciton resonance are
populated can be expressed as:

$\vec{\psi}(x)=\int\frac{dk}{2\pi}e^{ikx}[c_{+,\,-}(k)\vec{w}_{+,\,-}(k)+c_{-,\,-}(k)\vec{w}_{-,\,-}(k)],$
(2bj)

where $c_{\alpha,\,\beta}(k)$ are complex coefficients. In the present
analysis, we used a set of $\\{c_{\alpha,\,\beta}(k)\\}$ following a Gaussian
distribution in $k$ space, namely

$c_{\alpha,\,\beta}(k)\propto\exp(-(k-k_{0})^{2}/(2\Delta k)^{2}),$ (2bk)

with average $k_{0}$ and standard deviation $\Delta k$ drawn from two
different uniform random distributions. Such sets of coefficients are then
normalized in order to have an initial occupation corresponding to one
particle, that is

$\begin{split}&N_{\psi}(t=0)=\int
dx\langle\vec{\psi}(x),\,\vec{\psi}(x)\rangle=\\\
&=\int\frac{dk}{2\pi}\left[|\tilde{c}_{+,\,-}(k)|^{2}+|\tilde{c}_{-,\,-}(k)|^{2}\right]\,=1.\end{split}$
(2bl)

As already mentioned in the main text, the latter configuration is then
evolved forward in time, and polariton condensation is monitored by looking at
the time-behavior of the condensate occupation $N_{\psi}(t)$ for $t>0$.

## References

  * Sanvitto and Kéna-Cohen [2016] D. Sanvitto and S. Kéna-Cohen, The road towards polaritonic devices, Nature Materials 15, 1061 (2016).
  * Carusotto and Ciuti [2013] I. Carusotto and C. Ciuti, Quantum fluids of light, Rev. Mod. Phys. 85, 299 (2013).
  * Kasprzak  _et al._ [2006] 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 443, 409 (2006).
  * Balili  _et al._ [2007] R. Balili, V. Hartwell, D. Snoke, L. Pfeiffer, and K. West, Bose-Einstein condensation of microcavity polaritons in a trap, Science 316, 1007 (2007).
  * Amo  _et al._ [2009] 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 10.1038/nphys1364 (2009).
  * Lagoudakis  _et al._ [2008] K. G. Lagoudakis, M. Wouters, M. Richard, A. Baas, I. Carusotto, R. André, L. S. Dang, and B. Deveaud-Plédran, Quantized vortices in an exciton–polariton condensate, Nature Physics 4, 706 (2008).
  * Kavokin  _et al._ [2008] A. Kavokin, J. Baumberg, G. Malpuech, and F. Laussy, Microcavities, Microcavities , 1 (2008).
  * Azzini  _et al._ [2011] 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 10.1063/1.3638469 (2011).
  * Ballarini  _et al._ [2013] 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 4, 1778 (2013).
  * Paschos  _et al._ [2018] G. G. Paschos, T. C. H. Liew, Z. Hatzopoulos, A. V. Kavokin, P. G. Savvidis, and G. Deligeorgis, An exciton-polariton bolometer for terahertz radiation detection, Scientific Reports 8, 10092 (2018).
  * Andreani and Gerace [2006] 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 73, 235114 (2006).
  * Gerace and Andreani [2004] D. Gerace and L. C. Andreani, Gap maps and intrinsic diffraction losses in one-dimensional photonic crystal slabs, Phys. Rev. E 69, 056603 (2004).
  * Hsu  _et al._ [2016] C. W. Hsu, B. Zhen, A. D. Stone, J. D. Joannopoulos, and M. Soljačić, Bound states in the continuum, Nature Reviews Materials 1, 10.1038/natrevmats.2016.48 (2016).
  * Azzam and Kildishev [2021] S. I. Azzam and A. V. Kildishev, Photonic bound states in the continuum: From basics to applications, Advanced Optical Materials 9, 2001469 (2021).
  * Koshelev  _et al._ [2018] 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 98, 161113(R) (2018).
  * Zanotti  _et al._ [2022] S. Zanotti, H. S. Nguyen, M. Minkov, L. C. Andreani, and D. Gerace, Theory of photonic crystal polaritons in periodically patterned multilayer waveguides, Phys. Rev. B 106, 115424 (2022).
  * Dang  _et al._ [2022] 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 10, 2102386 (2022).
  * Kim  _et al._ [2021] 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 21, 10076 (2021).
  * Ardizzone  _et al._ [2022] 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 605, 447 (2022).
  * Riminucci  _et al._ [2022] F. Riminucci, V. Ardizzone, L. Francaviglia, M. Lorenzon, C. Stavrakas, S. Dhuey, A. Schwartzberg, S. Zanotti, D. Gerace, K. Baldwin, L. N. Pfeiffer, G. Gigli, D. F. Ogletree, A. Weber-Bargioni, S. Cabrini, and D. Sanvitto, Nanostructured $\mathrm{Ga}\mathrm{As}$/($\mathrm{Al},\mathrm{Ga}$)$\mathrm{As}$ waveguide for low-density polariton condensation from a bound state in the continuum, Phys. Rev. Appl. 18, 024039 (2022).
  * Andreani and Pasquarello [1990] L. C. Andreani and A. Pasquarello, Accurate theory of excitons in GaAs-Ga1-xAlxAs quantum wells, Physical Review B 42, 8928 (1990).
  * Savona  _et al._ [1995] V. Savona, L. C. Andreani, P. Schwendimann, and A. Quattropani, Quantum well excitons in semiconductor microcavities: Unified treatment of weak and strong coupling regimes, Solid State Communications 93, 733 (1995).
  * Gerace and Andreani [2007] D. Gerace and L. C. Andreani, Quantum theory of exciton-photon coupling in photonic crystal slabs with embedded quantum wells, Phys. Rev. B 75, 235325 (2007).
  * Bajoni  _et al._ [2009] 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 80, 201308(R) (2009).
  * Zhang  _et al._ [2018] L. Zhang, R. Gogna, W. Burg, E. Tutuc, and H. Deng, Photonic-crystal exciton-polaritons in monolayer semiconductors, Nature Communications 9, 713 (2018).
  * Dang  _et al._ [2020] N. H. M. Dang, D. Gerace, E. Drouard, G. Trippé-Allard, F. Lédée, R. Mazurczyk, E. Deleporte, C. Seassal, and H. S. Nguyen, Tailoring dispersion of room-temperature exciton-polaritons with perovskite-based subwavelength metasurfaces,  _Nano Letters_ , Nano Letters 20, 2113 (2020).
  * Chen  _et al._ [2020] 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_ , Nano Letters 20, 5292 (2020).
  * Gianfrate  _et al._ [2023] A. Gianfrate, H. Sigurdsson, V. Ardizzone, H. C. Nguyen, F. Riminucci, M. Efthymiou-Tsironi, K. W. Baldwin, L. N. Pfeiffer, D. Trypogeorgos, M. D. Giorgi, D. Ballarini, H. S. Nguyen, and D. Sanvitto, Optically reconfigurable molecules of topological bound states in the continuum (2023), arXiv:2301.08477 [physics.optics] .
  * Wouters and Carusotto [2007] M. Wouters and I. Carusotto, Excitations in a nonequilibrium bose-einstein condensate of exciton polaritons, Phys. Rev. Lett. 99, 140402 (2007).
  * Baboux  _et al._ [2018] F. Baboux, D. D. Bernardis, V. Goblot, V. N. Gladilin, C. Gomez, E. Galopin, L. L. Gratiet, A. Lemaître, I. Sagnes, I. Carusotto, M. Wouters, A. Amo, and J. Bloch, Unstable and stable regimes of polariton condensation, Optica 5, 1163 (2018).
  * Fontaine  _et al._ [2022] Q. Fontaine, D. Squizzato, F. Baboux, I. Amelio, A. Lemaître, M. Morassi, I. Sagnes, L. Le Gratiet, A. Harouri, M. Wouters, I. Carusotto, A. Amo, M. Richard, A. Minguzzi, L. Canet, S. Ravets, and J. Bloch, Kardar-parisi - zhang universality in a one-dimensional polariton condensate, Nature 608, 687 (2022).
  * Amelio  _et al._ [2023] I. Amelio, A. Chiocchetta, and I. Carusotto, Kardar-parisi-zhang universality in the linewidth of non-equilibrium 1d quasi-condensates, (2023), arXiv:2303.03275 [cond-mat.stat-mech] .
  * Wouters  _et al._ [2008] M. Wouters, I. Carusotto, and C. Ciuti, Spatial and spectral shape of inhomogeneous nonequilibrium exciton-polariton condensates, Phys. Rev. B 77, 115340 (2008).
  * Lagoudakis  _et al._ [2011] K. G. Lagoudakis, F. Manni, B. Pietka, M. Wouters, T. C. H. Liew, V. Savona, A. V. Kavokin, R. André, and B. Deveaud-Plédran, Probing the dynamics of spontaneous quantum vortices in polariton superfluids, Phys. Rev. Lett. 106, 115301 (2011).
  * Lu  _et al._ [2020] L. Lu, Q. Le-Van, L. Ferrier, E. Drouard, C. Seassal, and H. S. Nguyen, Engineering a light–matter strong coupling regime in perovskite-based plasmonic metasurface: quasi-bound state in the continuum and exceptional points, Photonics Research 8, A91 (2020).