# Soliton formation in an exciton-polariton condensate at a bound state in the continuum

I. Septembre,¹ I. Foudjo,¹ V. Develay,² T. Guillet,² S. Bouchoule,³

J. Zúñiga-Pérez,⁴,⁵ D. D. Solnyshkov,¹,⁶ and G. Malpuech¹

¹Université Clermont Auvergne, Clermont Auvergne INP,
CNRS, Institut Pascal, F-63000 Clermont-Ferrand, France

²Laboratoire Charles Coulomb, Université de Montpellier, CNRS, 34095 Montpellier, France

³Centre de Nanosciences et de Nanotechnologies, CNRS, Université Paris-Saclay, France

⁴UCA, CRHEA-CNRS, Rue Bernard Gregory, 06560 Valbonne, France

⁵MajuLab, International Research Laboratory IRL 3654, CNRS,
Université Côte d'Azur, Sorbonne Université, National University of Singapore,
Nanyang Technological University, Singapore, Singapore

⁶Institut Universitaire de France (IUF), 75231 Paris, France

Bound states in the continuum (BIC) are of special interest in photonics due to their theoretically infinite radiative lifetime. Here, we design a structure composed of a GaN layer with guided exciton-polaritons and a TiO₂ 1D photonic crystal slab. The photonic BIC hosted by the photonic crystal slab couples with the excitons of GaN to form a polaritonic BIC with a negative mass. This allows condensation to be reached with a low threshold in a structure suitable for electrical injection, paving the way for room-temperature polariton microdevices. We study in detail how the repulsive interaction between exciton-polaritons affects the condensate distribution in reciprocal space and, consequently, the condensate's overlap with the BIC resonance and, therefore, the condensate lifetime. We study an intrinsic contribution related to the formation of a bright soliton and the extrinsic contribution related to the interaction with an excitonic reservoir induced by spatially focused non-resonant pumping. We then study the peculiar dynamics of the condensation process in a BIC state for interacting particles using Boltzmann equations and hybrid Boltzmann-Gross Pitaevskii equations. We find optimal conditions allowing one to benefit from the long lifetime of the BIC for polariton condensation in a real structure.

## I. INTRODUCTION

A large part of Physics is built upon the conservation of quantities, such as energy or momentum, related to symmetries¹. This is the case of the Hamiltonian description of classical mechanics and then of quantum mechanics. However, the irreversible decay of a particle or a mode because of its coupling to a continuum is ubiquitous and, therefore, essential to describe various phenomena. Techniques have been developed, for example in quantum optics, which separate the reversible Hamiltonian dynamics and the dissipative dynamics related to the irreversible decay towards a continuum². A somewhat simpler description of multimode non-conservative systems came with the development of non-Hermitian theories, where modes have complex energies with imaginary parts describing their decay³. The coupling between two such modes hybridizes them but also redistributes⁴ their imaginary parts Γ₁,₂. In the case of a Hermitian coupling, the maximal redistribution occurs at exceptional points (EPs), where the mode decay becomes the mean value of that of the uncoupled modes (Γ₁ + Γ₂)/2. It can be zero only in systems with equal gain and loss (PT-symmetric) for the original modes. The study of non-Hermitian systems and in particular of EPs became extremely popular in the last decade⁵, partly because of the possibility offered to control the decay of the modes. This is essential for many applications including lasing⁶, and can even lead to fundamentally new concepts, such as the

non-Hermitian skin effect⁷,⁸.

Extra opportunities in this quest of controlling lifetimes came from the so-called dissipative coupling introduced in the 80s⁹. The idea is that two modes are coupled to the same continuum and because of this coupling a transfer of particles from one mode to another can take place with a rate computed as $i\sqrt{\Gamma_1\Gamma_2} 10^{10}$. This non-Hermitian term offers a new resonance condition, where one eigenmode of the coupled system becomes entirely non-decaying. This type of mode is called bound state in the continuum (BIC)¹¹. It was described long ago in the context of quantum mechanics¹², but BICs in photonic systems appear very appealing and commonly present¹³. BICs in photonic crystal slabs or metasurfaces are modes that lie above the light cone but are not coupled to the far field. One should note that other decay mechanisms, such as disorder-induced elastic scattering towards other modes, are not suppressed for BICs, which means that, in reality, their lifetime is finite. In the 2000s, when photonic crystal slabs (PCS) became widely studied, a very important goal was to improve the lifetime of the modes to realize lasers, but the BIC concept had not yet been established. Nevertheless, experimental and theoretical works found that modes with a long lifetime could be realized near the Γ point of photonic crystals dispersion¹⁴,¹⁵. Here, the BICs form when two modes of the same polarization and parity couple through the continuum.

Recently, BICs have been implemented in strongly
coupled exciton-photon systems where exciton-polaritons (polaritons) form<sup>16,17</sup>. As expected, non-resonantly pumped polaritons were found to accumulate in the state with the longest lifetime, namely the BIC-state, forming a Bose-Einstein Condensate (BEC)<sup>18,19</sup>. A principal advantage of a polariton condensate versus a photon laser is that the transparency condition required for lasing is not necessary for polariton lasing<sup>20,21</sup>. The only lasing condition for polaritons (which is similar for photon lasers) is that scattering towards a state overcomes its losses, which are expected to be essentially radiative losses for polaritons and that vanish for BICs. Another key difference between polaritons and bare photons is that polaritons are interacting particles. They interact with themselves and with reservoir excitons created by an external pumping. Depending on the sign of the effective mass these interactions, which are typically repulsive, can either tend to delocalize particles<sup>22</sup> in real space or localize them, forming bright solitons<sup>23,24</sup> or simply reservoir-induced localized states<sup>25</sup>. In systems with both positive and negative mass states, the localization near the exciton reservoir induced by the pump on negative mass states tends to favor condensation in these states<sup>23,25,26</sup>. An important consequence of the formation of either soliton-like BEC or of BEC bound to a spatially-localized exciton reservoir is that it broadens the condensate distribution in k-space. In a system containing a BIC, this k-space broadening can considerably reduce the average particle lifetime in the BEC. Finally, BIC-based polariton BEC has been demonstrated in GaAs-based structures at low temperature<sup>18,19</sup>. The potential use of polariton-based devices requires room temperature operation, which is typically implemented either using large-bandgap semiconductors (ZnO<sup>27-29</sup>, GaN<sup>30-32</sup>), organic materials<sup>33,34</sup> or perovskites<sup>35-37</sup>. In that perspective, GaN-based devices remain the best potential choice since this material offers the best stability as an inorganic semiconductor, and also the possibility to fabricate electrically-injected devices.

In this work, we first propose a design of a realistic BIC-based polariton laser in a GaN/TiO₂ polaritonic crystal slab, using finite element methods to solve Maxwell's equations. We then focus on polariton-specific features, particularly the role played by interactions. In the conservative limit we determine the shape of the polariton BEC with and without the reservoir potential. This allows us to analytically compute the average BEC particle lifetime in both cases. We then use the computed lifetimes to model quantitatively the polariton condensation using Boltzmann equations. We find that the advantage of having a BIC is strongly reduced, if not cancelled, when the interaction with a reservoir reduces too strongly the condensate extension in real space, explaining recent experimental results<sup>19</sup>. We find that the formation of a soliton provides anyway an intrinsic limitation to the lifetime enhancement provided by the BIC, although the threshold can still be reduced by more than one order of magnitude. At large pumping density, the

formation of a soliton provides saturation of the condensate density growth rate versus pumping, which scales in $P^{1/3}$ instead of being linear. Finally, we model the polariton condensation process, solving numerically hybrid Boltzmann-Gross-Pitaevskii equations<sup>38</sup>, which allows us to take into account self-consistently the energy relaxation and the inter-particle interactions.

## II. DESIGN OF THE STRUCTURE

The structure we study is sketched in Figure 1(a). It is composed of a TiO₂ slab deposited on top of a GaN waveguide. The GaN slab is isolated from the bulk GaN substrate by a cladding layer, that we choose to be $Al_xGa_{1-x}N$, with $x \approx 0.2$, as often done in experiments<sup>32</sup>. The TiO₂ layer is patterned with a 1D step pattern, forming a photonic crystal slab. Its lattice constant is denoted $a_0$. The width of the steps is $a_1$, while the space between them is $a_2$, with $a_1 + a_2 = a_0$.

We simulate the structure in COMSOL and find the dispersion relation of the photonic modes $E_p(k)$, together with their radiative losses $\Gamma_r(k)$. The structural parameters we use are $a_0 = 135$ nm, $a_1 = 80$ nm, $a_2 = 55$ nm, $h_0 = 50$ nm, and $h_{GaN} = 100$ nm. The dispersion relation is plotted in Figure 1(b) (dashed lines). We see that there are two photonic branches, both behaving quadratically for low wave vectors. We call the negative-mass photonic branch the "BIC" branch because it contains the BIC (at $k=0$), while the branch at higher energy is called the "lossy" branch because the states are lossy at all wave vectors. A dashed arrow indicates the photonic BIC (state) at $k=0$. We can approximate the photonic dispersion relation with a quadratic dispersion of the form:

$$E_{\text{p}}^{\text{BIC/lossy}}(k) = E_{0,p}^{\text{BIC/lossy}} + \frac{\hbar^2 k^2}{2m_{\text{p}}^{\text{BIC/lossy}}}, \quad (1)$$

where $\hbar$ is the reduced Planck's constant and $m_p$ the photonic effective mass. $E_{0,p}$ is the energy of the branches at $k=0$. The lossy branch has a positive mass, while the BIC branch has a negative one in our case. This can be inverted by tuning the parameters, particularly the filling factor $a_1/a_0$.

In Figure 1(c), we show the radiative losses (imaginary part of the energy, $\hbar\Gamma$) of the photonic modes. We can see that while the radiative losses of the BIC branch are zero for $k=0$, the losses of the lossy branch are always nonzero. Similarly to the real part of the energy, the imaginary part can be approximated by a square dependence on the wave vector:

$$\Gamma_{\text{r}}^{\text{BIC/lossy}}(k) = \Gamma_{0}^{\text{BIC/lossy}} + \gamma_{\text{r}}^{\text{BIC/lossy}} k^2, \quad (2)$$

where $\Gamma_{0r}$ are the radiative losses at $k=0$. For the BIC branch, this term is zero because the losses vanish at the BIC at $k=0$:

$$\Gamma_{\text{r}}^{\text{BIC}}(k) = \gamma_{\text{r}} k^2. \quad (3)$$
We can calculate the lifetime of the modes from the radiative losses as:

$$ \tau_r(k) = \frac{1}{2\Gamma_r(k)}, \qquad (4) $$

and we see that the radiative lifetime of the BIC is theoretically infinite, as can be seen in Figure 1(d). However, several limitations prevent real BICs from having an infinite lifetime. The first limitation results from the finite size of the structure. Indeed, the size of the full structure L limits the minimal wave vector that can be achieved $k_{min} \approx \frac{2\pi}{L}$. The radiative lifetime of the BIC branch is infinite only at $k = 0$, so since this value is not achievable, the radiative lifetime of an ideal photonic BIC is limited to:

$$ \tau_{\max}^{\text{BIC}} = \frac{L^2}{8\pi^2\gamma_r} \qquad (5) $$

This is represented in Figure 1(e) as a green line: the larger the structure, the longer the radiative lifetime, but a real (finite) structure will never exhibit an infinite lifetime. Moreover, non-radiative losses can also play an important role, as we explore in the following.

We designed the structure so that the energies of interest are close to the excitonic resonance of GaN. The photonic modes will thus couple to the excitonic modes and form exciton-polaritons, which are part-light, part-matter quasiparticles. To find the polaritonic dispersion relation out of the photonic dispersion relation containing the lossy branch and the BIC branch, we consider the strong coupling between excitons and photons through the effective strong coupling Hamiltonian:

$$ M_{\text{SC}} = \begin{pmatrix} E_X(k) & \rho \hbar \Omega_R / 2 & 0 & 0 \\ \rho \hbar \Omega_R / 2 & E_p^{\text{lossy}}(k) & 0 & 0 \\ 0 & 0 & E_X(k) & \rho' \hbar \Omega_R / 2 \\ 0 & 0 & \rho' \hbar \Omega_R / 2 & E_p^{\text{BIC}}(k) \end{pmatrix}, \qquad (6) $$

where $E_X(k) \approx E_X = 3480$ meV at 5 K and 3520 meV at room temperature$^{39,40}$ is the effective excitonic resonance of GaN, which groups both A and B exciton contributions. The two different photonic modes have distinct spatial distributions and are coupled to distinct spatial distributions of excitons, which explains the structure of the Hamiltonian. $E_p(k)$ are the photonic dispersions found in COMSOL, and we take the Rabi splitting $\hbar\Omega_R = 100$ meV for GaN$^{39}$. $\rho$ and $\rho'$ are the fractions of the lossy/BIC modes, respectively, that are located in the GaN layer. Indeed, since the waveguide is composed not only of GaN, but also of a TiO$_2$ photonic crystal slab, which does not contain excitons, we take it into account in the overlap integral between the photon and excitons wave function, which sets the Rabi splitting value, as we did previously$^{41}$. In our design, the value of $\rho$ is typically around 60%, while the value of $\rho'$ is close to 75%. The polaritonic dispersion arising from the strong coupling between excitons and photons is determined by finding

the eigenvalues of the matrix $M_{\text{SC}}$. The polaritonic dispersions read:

$$ E_{\pm}(k) = \frac{E_X + E_p(k)}{2} \pm \frac{1}{2} \sqrt{(\rho \hbar \Omega_R)^2 + (E_X - E_p(k))^2}, \qquad (7) $$

where the + sign corresponds to the upper polariton branches (UPB) with energies above the exciton resonance, and the - sign corresponds to the lower polariton branches (LPB) whose energies lie below the exciton resonance.

We plot the resulting polaritonic dispersion in solid lines in Figure 1(b). The exciton energy is represented as a dotted black line, which is horizontal because the mass of the excitons is orders of magnitude larger than the effective mass of photons. One obtains four polaritonic branches from the two excitonic and two photonic modes. The UPB and LPB corresponding to the BIC (blue) and lossy (red) branches are shown in the figure. The main branch of interest for us is the BIC LPB, represented as a solid blue line with an energy below the exciton energy. The BIC is located at $k = 0$ and indicated by a black arrow. We see that at low wave vectors, the dispersion relation of this branch can be approximated by a quadratic dependence on the wave vector:

$$ E(k) \approx E_0^{\text{BIC}} + \frac{\hbar^2 k^2}{2m}, \qquad (8) $$

where $m = -(1.5 \pm 0.1) \times 10^{-6} m_0$ is the fitted polariton mass, $m_0$ being the free electron mass. This effective mass is negative. $E_0^{\text{BIC}}$ is the energy of the BIC LPB branch at $k = 0$.

In the strong coupling regime, the imaginary part of the polariton energy $\Gamma$ can be expressed as$^{21}$:

$$ \Gamma = C_X(k)\Gamma_X + (1 - C_X(k))\Gamma_p(k), \qquad (9) $$

where $\Gamma_{p,X}$ are the losses of the photonic/excitonic modes respectively, and $C_X$ is the exciton fraction. This fraction can be modulated by engineering the dispersion relation; more precisely, it depends on the detuning between the exciton energy and the photonic energy at $k = 0$. Thus, by engineering the structure (period, etc), one can increase the excitonic part of the polaritonic modes and consequently increase the interactions, which have been experimentally measured in GaN just recently$^{42}$.

The lifetime of the polaritons, as that of the photonic modes, can be fitted by a quadratic dependence on the wave vector for small $k$. The approximate expression is given by:

$$ \hbar\Gamma_{\text{pol}}(k) \approx \hbar\Gamma_0 + \Lambda\frac{\hbar^2 k^2}{2m}, \qquad (10) $$

where $\Gamma_0$ are the constant losses and $\Lambda = -0.27 \pm 0.03$ is obtained from a fitting. $\Lambda$ is negative here because $m$ is negative while the losses increase with the wave vector for the BIC branch.

Ultimately, the non-radiative lifetime introduced by the excitons will limit the lifetime of a polaritonic BIC.
<figure><img src="image_1.png" /><figcaption>Figure 1: (a) Sketch of the photonic crystal slab structure studied. A two-layer waveguide (TiO₂ in yellow and GaN in pink) is separated from the substrate by a cladding layer (AlGaN in purple). The TiO₂ layer is etched with a step-like profile (step width a₁). (b) Photonic dispersion relation (dashed lines) for the branch containing a BIC (BIC branch, in blue) and the branch containing no BIC (lossy branch, in red). The corresponding polaritonic dispersion relations are plotted in solid lines (red/blue for the polariton dispersion relation calculated from the photonic lossy/BIC branch, respectively). There are two polariton branches (Upper/Lower Polariton Branch) for each photonic branch. The exciton effective energy is indicated as a dotted line. We extract m from the BIC LPB. (c) Imaginary part of the energy of the photonic branches (dashed lines) and lower polaritonic branches (solid lines) versus wave vector. The red/blue curve corresponds to the lossy/BIC branch. We extract Λ from the BIC LPB. (d) Photonic/polaritonic (dashed/solid line) lifetime (log scale) against wave vector. (e) Maximum lifetime achievable depending on the structure length L. The blue/pink lines indicate the radiative/non-radiative lifetimes. In (b-d), the polaritonic/photonic BIC at k = 0 is indicated with a solid/dashed black arrow.</figcaption></figure>

According to Ref.⁴³, the order of magnitude of the non-radiative excitonic lifetime is 50 ps at room temperature in GaN. This means that the longest lifetime possible for a polaritonic BIC in this structure is limited, even in an infinite structure. Figure 1(e) shows that the lifetime of the BIC is limited by the size of the structure when the structure is small but finally, for large structures, the non-radiative lifetime will limit the polaritonic lifetime.

Still, depending on the parameters, very high lifetimes can be observed in realistically large structures. For a structure of 100 µm length, the lifetime is of the order of 20 ps (corresponding to a quality factor of Q ≈ 20000 in a λ/2 microcavity). This is a very high value at room temperature, even for guided modes⁴⁴ (and here, the mode is above the light cone). For larger structures, the lifetime is ultimately limited by the non-radiative lifetime, not by the structure's size.

### III. INFLUENCE OF THE INTERACTIONS ON THE LIFETIME

In this section, we use the k-dependent losses of the BIC to compute the losses of the states hosted by the structure. Indeed, as shown in the previous section, the decay rates (and lifetimes) of modes around the BIC are strongly k-dependent. An ideal non-interacting homogeneous condensate would occupy a single state k = 0 and benefit from the corresponding extremely long lifetime. However, different mechanisms can lead to a localization of the condensate in real space and thus to a broadening in k-space, leading to an increase in the overall condensate decay rate (decrease of lifetime). This contribution can be computed as:

$$ \langle \Gamma_{\text{pol},r} \rangle = \frac{1}{N} \int \Gamma_{\text{pol},r}(k) |\psi(k)|^2 dk. \quad (11) $$

where $\psi(k)$ is the k-component of the condensate wave function normalised to the number N of condensed particles and $\hbar\Gamma_{\text{pol},r}(k) = \Lambda\hbar^2k^2/2m$ is the wavevector-dependent part of the polariton decay rate (10).

There are two localization mechanisms in real space for negative-mass particles with repulsive interactions. The first one is an intrinsic mechanism related to the formation of a bright soliton, whose width is related to N. The second mechanism takes place if the pump is non-uniform (localized) in real space. The generated excitonic reservoir repulsively interacts with the condensate, playing the role of an attractive potential. Generally, both mechanisms are present²³. We first treat this problem analytically, using the conservative Gross-Pitaevskii (GP) equation:

$$ i\hbar \frac{\partial \psi}{\partial t} = -\frac{\hbar^2}{2m} \frac{\partial^2 \psi}{\partial x^2} + \alpha |\psi|^2 \psi + U(x) \psi \quad (12) $$

where $\alpha$ is the polariton-polariton interaction constant and $U(x)$ an external potential induced by the exciton reservoir.

So far, the losses have been neglected. We start with the intrinsic mechanism of bright soliton formation (without any external potential $U(x) = 0$). The bright soliton solution of the GP equation without external potential reads:

$$ \psi_s = \frac{\psi_s^{(0)}}{\cosh(x/\sqrt{2\xi})}, \quad (13) $$
where $\xi$ is the healing length defined as:

$$ \xi = \frac{\hbar}{\sqrt{2|m\alpha|}\psi_s^{(0)}} \quad (14) $$

and $\psi_s^{(0)}$ is the value of the wave function at $x = 0$, defined as:

$$ \psi_s^{(0)} = \frac{\sqrt{|m\alpha|N}}{2\hbar}, \quad (15) $$

which ensures the normalization of the wave function to $N$. Finally, the healing length can be rewritten explicitly as a function of $N$ as:

$$ \xi = \frac{\sqrt{2\hbar^2}}{|m\alpha|N}. \quad (16) $$

This shape allows us to compute analytically the decay rate in the presence of a BIC using formula (11) as:

$$ \langle \hbar \Gamma_r \rangle = \Lambda \frac{\alpha^2 m}{12 \pi^2 \hbar^2} N^2 \quad (17) $$

This non-linear loss rate will saturate the growth of $N$ after the condensation threshold. Qualitatively, the accumulation of particles in the soliton renders it larger in $k$-space, which increases the decay rate, which then limits the growth of $N$.

We now consider the localization by an external potential induced by the finite size of a Gaussian reservoir created by a non-resonant pumping:

$$ U(x) = U_0 e^{-\frac{x^2}{2\sigma^2}}, \quad (18) $$

where $U_0$ is the height of the potential "trap" formed by the pump and $\sigma$ is its width in real space. The full width at half maximum (FWHM) of the pump is $2\sqrt{2\ln(2)}\sigma$. We approximate the potential by a parabolic profile found from the Gaussian profile by a series expansion:

$$ U \approx U_0 \left(1 - \frac{x^2}{\sigma^2}\right). \quad (19) $$

The solutions to the Schrödinger equation for negative mass particles in this harmonic potential correspond to the solutions for positive mass particles but with a potential $-U$. The solutions are the well-known states of the harmonic oscillator, which are the products of Hermite polynomials and of a Gaussian function. This harmonic oscillator description is valid only for particle energies within the gap between the BIC and lossy branches. This limits the number of modes confined by the potential and, in practice, we will only consider the two first modes $\tilde{\psi}_{0,1}$ defined as:

$$ \tilde{\psi}_0 = \frac{1}{\sqrt{w\sqrt{\pi}}} e^{-\frac{x^2}{2w^2}}, \quad \tilde{\psi}_1 = \sqrt{2}\frac{x}{w}\tilde{\psi}_0. \quad (20) $$

where $w$ is the width of the Gaussian wave function. The $p$-state is expected to host a condensate when the energy

of the $s$-state becomes larger than the gap and is therefore not confined anymore$^{23,25,45,46}$.

Because of the presence of both the interactions and the potential, the Gross-Pitaevskii equation is not exactly solvable. We use the variational approach to find the solution analytically, using the width $w$ of the Gaussian wave function as a variational parameter. We start with a trial wave function of a Gaussian shape, given by $\psi_0 = \sqrt{N}\tilde{\psi}_0$ where $N$ is the number of particles (note that $\int |\psi_0|^2 dx = N$ while $\int |\tilde{\psi}_0|^2 dx = 1$). We then find the width $w$ that minimizes the total energy of the system $E_{\text{tot}}$, which is the sum of the kinetic, potential, and interaction energy:

$$ E_{\text{tot}} = E_{\text{kin}} + E_{\text{pot}} + E_{\text{int}}. \quad (21) $$

The minimization condition allows us to write an equation for $w$ that we solve:

$$ \frac{\partial E_{\text{tot}}}{\partial w} = 0. \quad (22) $$

The different energy terms are calculated from the trial wave function as follows$^{22}$:

$$ E_{\text{kin}} = \int -\psi_0^* \frac{\hbar^2}{2m} \frac{\partial^2 \psi_0}{\partial x^2} dx, \quad (23) $$

$$ E_{\text{pot}} = \int |\psi_0|^2 U dx, \quad (24) $$

$$ E_{\text{int}} = \int \frac{\alpha}{2} |\psi_0|^4 dx. \quad (25) $$

In the end, we find the total energy, written here per particle:

$$ E_{\text{tot}}/N = \frac{\hbar^2}{4mw^2} + \frac{w^2U_0}{\sigma^2} + \frac{\alpha N}{2\sqrt{2\pi w}}. \quad (26) $$

In the absence of interactions ($\alpha = 0$), the energy is minimized for a width$^{47,48}$:

$$ w_{|\alpha|=0} \equiv w_0 = \sqrt{\frac{\hbar\sigma}{\sqrt{|m|U_0}}}. \quad (27) $$

The trial function is an exact solution in this case. The decay rate can be found analytically. It reads:

$$ \langle \Gamma_{r,0} \rangle = \frac{\pi \Lambda \hbar}{2mw^2}. \quad (28) $$

A similar analysis can be made for state 1, and the decay rate reads:

$$ \langle \Gamma_{r,1} \rangle = 3 \langle \Gamma_{r,0} \rangle, \quad (29) $$

and more generally, the losses of the state $\mathcal{M}$ is given by:

$$ \langle \Gamma_{r,\mathcal{M}} \rangle = \left( \mathcal{M} + \frac{1}{2} \right) \frac{\pi \Lambda \hbar}{mw^2}. \quad (30) $$
When interactions are present, Eq. (22) becomes:

$$ \frac{U_0}{\sigma^2} w^4 - \frac{\alpha N}{\sqrt{2\pi}} w - \frac{\hbar^2}{2m} = 0. \quad (31) $$

This is a 4th-order polynomial equation. It is perfectly solvable; however, the expressions of the solutions are quite cumbersome. Among the four solutions only one is real positive, which is the only one we keep since w is a length. To understand the behaviour of the system, we consider analytically the two limits $N \to 0$ and $N \to \infty$, which give two different expressions for w by expansion in Laurent series. We find that:

$$ w_{|N \to 0} \equiv w_-(N) = w_0 - \frac{\alpha \sigma \sqrt{|m|}}{4\hbar\sqrt{2\pi U_0}} N, \quad (32) $$

and

$$ w_{|N \to \infty} \equiv w_+(N) = \frac{\sqrt{2\pi}\hbar^2}{|m\alpha|N}. \quad (33) $$

The two limits coincide for a value of N that we denote:

$$ N_{\text{crit}} = \frac{2\hbar}{\alpha} \sqrt{\frac{2\pi\hbar}{|m|\sigma} \sqrt{\frac{U_0}{|m|}}}, \quad (34) $$

obtained by solving the equation $w_− = w_+$ . This is the critical value for which interactions within the condensate become the dominant mechanism controlling its width. We can therefore write an approximate expression for the width $w(N)$ using the two asymptotic expressions:

$$ w \approx w_{-} \mathcal{U}(-(N - N_{\text{crit}})) + w_{+} \mathcal{U}(N - N_{\text{crit}}), \quad (35) $$

where $\mathcal{U}$ denotes the Heaviside step function. Figure. 2(a) shows the width of the 0 mode versus N, which is here considered as an external parameter. At low N, the mode width is $w_0$ governed by the reservoir potential. The width then decreases and becomes dominated by inter- action within the condensate for $N > N_{\text{crit}}$. At very large N, the width becomes the same as the one of the bare bright soliton $\xi$. The full solution of (31) is shown as a thick black line, while the two asymptotic solutions combined by (35) are shown by gold and orange lines.

We now consider the second confined state $\psi_1$, for which we perform the same type of analysis. We use the spectroscopic notation for the states, as it is often done since the polaritonic trap can be considered an ar- tificial atom. The states 0 and 1 therefore correspond to s and p, respectively. We remind that state 1 can host the condensate when the energy of state 0 overcomes the band gap between the BIC and lossy branches. We plot in Figure. 2(b-e) the real and k-space profiles of the dif- ferent modes for low N (b,d) and large N (c,e). For a low number of particles, the confined states are mostly deter- mined by the width of the confining potential, whereas a free soliton is much broader in real space. On the contrary, for large N, the interactions make all modes

<figure><img src="image_2.png" /><figcaption>Figure 2: (a) Width of the Gaussian wave function w (black thick line) and healing length of the soliton $\xi$ (blue line) as a function of the number of particles N. $w_+$ (respectively gold/orange lines) are the approximations of w for large/small N, and $w_0$ is the width in the non-interacting case. (b-c) Profile of the wave function in real space for $N = N_{\text{crit}}/10$ (b) and $N = 10N_{\text{crit}}$ (c). The confining potential (light blue) has a predominant effect at small N but a small effect at large N. (d-e) Profile of the wave function in reciprocal space for $N = N_{\text{crit}}/10$ (d) and $N = 10N_{\text{crit}}$ (e). The $k^2$ losses close to the BIC (green line) give very low losses at low N but high losses at large N. In (b-d), the dark blue curve is the soliton wave function and the orange/red curves are the wave functions of the states 0 and 1, respectively.</figcaption></figure>

narrower so that the confining potential (which does not change between (b) and (c)) has only a minor influence on the width of the confined states, which is close to that of the free soliton.

In k-space, the trend is inverse: the narrow states in real space are broad in k-space. The broadening in k- space for all states increases linearly with the growth of
$N$ (Eqs. (16),(33)). Another important aspect is that the p-state shows a node at $k = 0$ and shows, in all cases, a smaller overlap with the BIC resonance compared with the soliton and the s-states. On the other hand, for large $N$, all states are broad in reciprocal space. It means that most of their probability density is located outside the BIC point and, therefore, these states almost do not benefit from the extremely long BIC lifetime.

These features are summarized in Figure. 3(a), which shows the radiative losses of the soliton, the s-, and the p-states versus $N$. We observe the $N^2$ dependence of $\Gamma$ for the soliton case. We observe at $N = N_{crit}$ the transition for s- and p-states from a regime where $\Gamma$ is fixed by the pump size to a soliton-like regime dominated by the interactions within the condensate. We observe the constant ratio (equal to 3, see Eq. (29)) for the decay rates between the s- and p-states. The transition occurs at $N = N_{crit}$. All decay rates $\Gamma$ depend on the size of the confining potential $\sigma$. This is shown in Fig 3b, which also includes a comparison with the band that does not contain a BIC (lossy branch) but instead is characterized by a non-radiative decay rate $\langle \Gamma_{other} \rangle$ (for instance, due to the excitonic resonance when considering exciton-polaritons) and a radiative decay rate $\langle \Gamma_r \rangle$ which does not depend on $k$ (for instance, the losses of a usual mode above the light cone), so that:

$$ \langle \Gamma \rangle = \langle \Gamma_r \rangle + \langle \Gamma_{other} \rangle \quad (36) $$

We observe that for a small reservoir size, the condensate is considerably wider in k-space with respect to the range where BIC modifies the decay rate. This situation corresponds approximately to the observation of Ref. 45, where the pump size is approximately $\sigma \approx 5$ µm. On the contrary, when the pump size increases, the condensate wave function expands in real space and localizes in k-space, overlapping better and better with the BIC resonance. In that case, a band with a BIC shows a much smaller decay rate. This is the regime that was demonstrated in Ref. 18, where the pumping spot has a width $\sigma \approx 80$ µm. To summarize, using a small pumping spot prevents benefiting from the high lifetime of a BIC.

## IV. INFLUENCE OF THE INTERACTIONS ON THE CONDENSATION THRESHOLD

In the previous section, we found analytical formulas for the losses of the soliton and s- and p-states. We show that the losses depend on the number of particles $N$. In this section, we use the previously found density-dependent losses to solve the semi-classical Boltzmann equations 49, which describe the evolution of the population of the condensate and reservoir states with time. In the previous sections, the system was assumed to be 1D and uniform in the transverse direction. Here, we need to take into account a normalization surface in order to operate with the populations and scattering rates. Experimentally, operating with a finite-size pump simply

<figure><img src="image_3.png" /><figcaption>Figure 3: (a) Losses of the soliton (blue line) and s- and p-states with respect to the number of particles $N$ (in log-log scale). (b) Losses of the s-state (orange lines) and p-state (red lines) for a band with a BIC (solid lines) and without a BIC (dashed lines). The left (right) grey line corresponds to the situations of Ref. 45 (Ref. 18), and the black arrows represent the difference between losses of the s state for a non-BIC state and the p state for a BIC state.</figcaption></figure>

leads to additional non-radiative losses due to the repulsion of the particles from the excitation spot due to the interactions. The simplified 2-state model used in this section does not allow us to take into account the details of the relaxation processes that depend on the dispersion relation. Its goal is to focus on the effects of the density-dependent lifetime in the condensate, which is due to the soliton formation. The two coupled equations read:

$$ \dot{N}_c = W_{LO}N_x(1+N_c) - N_c\Gamma_c(N_c), \quad (37) $$

$$ \dot{N}_x = -W_{LO}N_x(1+N_c) - N_x\Gamma_x + P, \quad (38) $$

where $N_{c,x}$ are, respectively, the numbers of particles in the condensate/reservoir, $\Gamma_{c,x}$ are the condensate/reservoir losses, and $P$ denotes the pumping strength. $W_{LO}$ describes the relaxation of excitons through LO phonons, which we assume to be the dominant relaxation mechanism. Indeed, the design of the structure was optimized in order to have a resonant enhancement for the LO-phonon assisted relaxation: $E_0^{BIC} = E_x + \frac{3}{2}k_B T - E_{LO}$ at $T = 300$ K. $W_{LO}$ is estimated from Ref. 49 to $W_{LO} = 6 \times 10^9$ s$^{-1}$ for the surface we consider. We neglect the dependence of the scattering rate on the spatial shape of the condensate. We solve the Boltzmann equations over time until a stationary configuration is obtained. These simulations are repeated for different values of pumping $P$. This allows us to obtain the stationary population of the condensate for each pumping.

We consider three different cases, which we plot in Figure 4. In the first case, we consider a configuration where the condensate is not a BIC. The losses for this state (both radiative and non-radiative) are assumed to be k-independent $\Gamma_c = \langle \Gamma_{nr} \rangle + \langle \Gamma_{r,0} \rangle$. When increasing the pumping strength, we see that the blue curve exhibits a threshold around $P \approx 10^{13}$ s$^{-1}$, but the slopes before and after the threshold are the same, $N_c \sim P$. For large $N_c$, the stationary solution can be found analytically as:

$$ N_c \approx \frac{\hbar P}{\Gamma_c} - \frac{\Gamma_x}{W_{LO}}. \quad (39) $$
The two other curves show the solutions when the system hosts a BIC. The golden curve represents the solution without confining potential, where the condensate wavefunction is the one of a soliton. We see that using a BIC instead of an ordinary lossy state reduces the threshold required for condensation by approximately two orders of magnitude (the threshold is estimated to be reached for $N_c \approx 3$). Before and just after the threshold, the condensate population evolves linearly versus $P$, as in the case without the BIC. However, far above threshold, when $N_c \approx 10^4$, the dependence changes due to the formation of the soliton, whose losses scale as $N_c^2$, as shown in Eq. (17). In that regime, keeping only leading terms in $N_c$ one can find the asymptotic dependence of $N_c$ versus $P$ as

$$N_c \approx \left(\frac{\hbar P}{\beta}\right)^{1/3}, \qquad (40)$$

which corresponds to the $1/3$ slope we observe in the numerical results at large $N_c$.

Finally, the green curve is obtained when considering a BIC with both confinement by an exciton reservoir and self-interactions. In this case, the losses are given by Equation (29), and the width is taken from the solution of Equation (31) (the analytical exact solution, not the approximate solutions). We take the losses of state 1 (the $p$-state) in order to reproduce the situation studied experimentally in Ref. $^{45}$. The solution for state 0 ($s$-state) is not plotted because it is very similar to the one we found considering only the soliton without confinement. Indeed, for large $N_c$, the losses of state 0 are close to the soliton losses, and for small $N_c$, the solutions of Boltzmann equations are less influenced by $N$-dependent losses. We see that the green curve is an intermediate case between the soliton located at a BIC and a usual configuration with a BIC: the threshold is higher than with the BIC and without confinement, but the interactions limit the growth of $N_c$ at large values with the same power law. Modifying the width of the confinement potential allows changing the threshold value: the narrower the potential, the higher the threshold.

To conclude, in this section we have shown that the BIC can be beneficial to reduce the threshold. However, the interactions deeply modify the behavior at large pumping strengths/number of particles because the dependence between the two parameters becomes sublinear. Finally, the size of the reservoir can reduce the effects of the BIC.

## V. NUMERICAL SIMULATIONS OF POLARITON CONDENSATION WITH GROSS-PITAEVSKII EQUATION

In this section, we solve numerically the Gross-Pitaevskii equation to verify qualitatively the validity of

the hypotheses used and results obtained in the previous sections. The studies performed in the previous sections

<figure><img src="image_4.png" /><figcaption>Figure 4: Number of particles in the condensate $N_c$ versus pumping the strength $P$. The blue line corresponds to the situation without a BIC, the gold line to the situation with a BIC and a free soliton (no confinement), and the green line the situation with a BIC and confined states (the $p$-state in this case). Short vertical dashed lines indicate the condensation thresholds in each situation. Short vertical dashed-dotted lines indicate the slope breaks (1 to $1/3$) due to the presence of the BIC.</figcaption></figure>

were based on the hypothesis that the shape of the condensate wave function is a solution to the conservative GP equation. The role of the pumping and lifetime is to set the number of particles in the condensate, which is a wave function parameter. We did not consider that the $k$-dependent lifetime itself could deeply affect the shape of the condensate wave function. In order to take into account on an equal footing interaction and lifetime effects, we solve numerically the hybrid Boltzmann Gross-Pitaevskii equation, including the scattering from a localized excitonic reservoir and $k$-dependent lifetime, similar to what is usually done to describe energy relaxation$^{38,50}$. The goal here is to check the validity of the analytical predictions for the shape of the condensate and simulate the emission pattern observed in experiments, and not to obtain a quantitative estimate of the threshold, which is why we eliminate the equation for the reservoir. The equation reads:
$$i\hbar \frac{\partial}{\partial t} \psi(x,t) = \left( (i\Lambda - 1) \frac{\hbar^2}{2m} \frac{\partial^2}{\partial x^2} + \alpha (|\psi(x,t)|^2 + 2n_x(x)) + i\hbar(W_{LO}n_x(x)e^{-N_x/N_{sat}} - \Gamma_0) \right) \psi(x,t) + \chi(x,t). \quad (41)$$

where $m = -(1.5 \pm 0.1) 10^{-6} m_0$, $\Lambda = -(0.27 \pm 0.03)$ as obtained in the first section. $n_x(x)$ is the reservoir density, which we assume to keep the spatial profile of the pumping laser and to be constant in time. The reservoir depletion is taken into account via the saturated gain term. We add a noise term $\chi(x,t)$ to simulate the spontaneous scattering into the polaritonic modes from the reservoir. The filling of the condensate from an excitonic reservoir generated by a non-resonant optical pumping is assumed to be assisted by the LO phonons through the term containing the LO phonon scattering matrix element $W_{LO}$. The saturation is controlled by the total number of particles in the reservoir $N_x = \int n(x) dx$ and its saturation value $N_{sat}$, which in the Boltzmann equations is controlled by the pumping $P$ via the second equation (38) that we eliminate adiabatically here. The interaction constant $\alpha$ describes polariton-polariton interaction. A twice larger value for the reservoir excitons is taken in the hypothesis of zero detuning operation. Here, the normalization area for both $\alpha$ and $W_{LO}$ is taken equal to the spatial grid size $2 \ \mu m$ times the transverse pump size $20 \ \mu m$. The main difference between our design and the one presented in Ref.⁵¹ is that our BIC is a polaritonic BIC and not just a photonic BIC, which enlarges the potential applications of our proposal compared to this previous realization. The goal of these simulations is to demonstrate that the interactions lead to the formation of a soliton, that its reciprocal-space image corresponds to the one observed experimentally, and that its shape is very close to the one considered analytically.

We show the results of numerical simulations based on Eq. (41) in Figure 5 using a homogeneous reservoir profile. The time-dependent simulation is run until a stationary solution is achieved. Above threshold (which in our simulations occurs well below the Mott density), a condensate forms in the ground state of the system, but the interactions lead to its localization in real space and spreading in reciprocal space, as expected. Panel (a) shows the condensate distribution in reciprocal space $|\psi(k)|^2$. We can see that the population is maximal at the top of the dispersion (the dispersion is indicated by a cyan dashed line), which is precisely the position of the BIC. Indeed, the BIC has fewer losses than any other mode, so condensation occurs there.

In most of the systems exhibiting polariton condensation, a higher population corresponds to a higher intensity of emitted light (but even for usual polaritons, one has to take into account the photonic fraction to determine the polariton population from the emission intensity⁵²). In the case of BICs, the population is very

<figure><img src="image_5.png" /><figcaption>Figure 5: Polariton Bose-Einstein condensate in a BIC. (a,b) Numerical results of the simulation of the Gross-Pitaevskii equation (41) for the BIC LPB (containing a BIC at $k = 0$). (a) Population in the reciprocal space. Note the unique intensity peak at the maximum of the dispersion (at the BIC). (b) Numerically calculated emission in the reciprocal space. Note the two peaks surrounding $k = 0$. In both panels, the dispersion is indicated as a cyan dashed line, and a white arrow shows the position of the BIC (state) in reciprocal space. (c) Numerically obtained distribution of particles in real space (blue points) after the stationary regime is reached and fitted by a bright soliton wave function (red line). Note the excellent agreement.</figcaption></figure>

different from the emission. Indeed, the light emitted by the structure is proportional to the losses: the light that is collected is the light that escaped the waveguide and impinges onto the detector. The emission $\Gamma_r(k)|\psi(k)|^2$ is plotted in panel (b). We can see two peaks close to the maximum of the dispersion, corresponding to the experimental observation of the polariton condensation at the BICs¹⁸. The population is huge in the BIC, but the emission is very low precisely because the lifetime is very high, larger than modes close in reciprocal space, which still have a very high population. In the end, we only see two peaks of emission close to the BIC, but not the BIC itself.

Our numerical simulations with the Gross-Pitaevskii equation confirm the formation of bright solitons due to the interactions. We plot in Figure 5(c) the distribution of particles in real space calculated numerically after condensation and after the system has stabilized (the condensation occurs after 20 ps and the picture is taken after 15 ns). We see that there is one peak (blue numerical
points) that is well-fitted by the wave function of a bright soliton (red line), as given by Eq. (13). The shape corresponds very well to the profile found in the conservative approximation in the previous sections, which confirms the validity of our analytical considerations. Note that the formation of solitons and their interplay with BICs in a different (resonant) pumping configuration have been theoretically discussed in Refs.⁵³,⁵⁴. However, we observe numerically an additional effect. When the number of particles exceeds a critical value, the system becomes unstable and does not converge to a stationary solution. This may lead to the formation of several solitons that exhibit periodic oscillations, without merging into a single soliton. According to our analysis, this is not due to the Kibble-Zurek mechanism⁵⁵,⁵⁶, which has already been proposed to explain different effects in polaritonic structures⁵⁷⁻⁵⁹. The effect is rather due to an analogue of modulational instability, but in a driven-dissipative system. A detailed study of this phenomenon exceeds the scope of the present article and requires further investigations.

## VI. CONCLUSIONS

To conclude, we design a structure for polariton Bose-Einstein condensation (polariton lasing) based on a 1D BIC in a TiO₂-GaN waveguide. We predict a possible realization of room temperature polariton Bose-Einstein

condensation in the BIC state. This could be implemented using state-of-the-art experimental techniques. It would allow obtaining an electrically-injected room-temperature polariton laser with a low lasing threshold. We also investigate how to benefit from the high lifetime of the BIC. It requires working with the largest structure possible to access the smallest wave vector and the highest lifetimes, and the pumping spot needs to be sufficiently large as well, to prevent confinement effects from reducing the lifetime of the BIC. Soliton formation due to the interactions can increase the losses above threshold. An interesting perspective is to consider condensation of interacting particles in 2D-BIC.

## Acknowledgments

We thank H. S. Nguyen and X. Letartre for participating in the early development of this work, P. Kokhanchik for contributing to the supervision of I. Fudjo, and E. Cambril and J.-Y. Duboz for useful discussions. This research was supported by the ANR Labex GaNext (ANR-11-LABX-0014), the ANR program "Investissements d'Avenir" through the IDEX-ISITE initiative 16-IDEX-0001 (CAP 20-25), the ANR project "NEWAVE" (ANR-21-CE24-0019) and the European Union's Horizon 2020 program, through a FET Open research and innovation action under the grant agreement No. 964770 (TopoLight).

1. E. Noether, Nachrichten von der Gesellschaft der Wissenschaften zu Göttingen, Mathematisch-Physikalische Klasse **1918**, 235 (1918), URL http://eudml.org/doc/59024.

2. G. Lindblad, Communications in Mathematical Physics **48**, 119 (1976).

3. M. Król, I. Septembre, P. Oliwa, M. Kędziora, K. Łempicka-Mirek, M. Muszyński, R. Mazur, P. Morawiak, W. Piecek, P. Kula, et al., Nature Communications **13**, 5340 (2022).

4. C. M. Bender and S. Boettcher, Phys. Rev. Lett. **80**, 5243 (1998), URL https://link.aps.org/doi/10.1103/PhysRevLett.80.5243.

5. E. J. Bergholtz, J. C. Budich, and F. K. Kunst, Reviews of Modern Physics **93**, 015005 (2021).

6. B. Peng, Ş. K. Özdemir, M. Liertzer, W. Chen, J. Kramer, H. Yılmaz, J. Wiersig, S. Rotter, and L. Yang, Proceedings of the National Academy of Sciences **113**, 6845 (2016).

7. S. Yao and Z. Wang, Physical Review Letters **121**, 086803 (2018).

8. P. Kokhanchik, D. D. Solnyshkov, and G. Malpuech, Physical Review B **108**, L041403 (2023).

9. H. Friedrich and D. Wintgen, Physical Review A **32**, 3231 (1985).

10. A. Z. Devdariani, V. N. Ostrovskii, and Y. N. Sebyakin, Sov. Phys. JETP **44**, 477 (1976).

11. C. W. Hsu, B. Zhen, A. D. Stone, J. D. Joannopoulos, and M. Soljačić, Nature Reviews Materials **1**, 1 (2016).

12. J. von Neuman and E. Wigner, Physikalische Zeitschrift **30**, 467 (1929).

13. S. I. Azzam and A. V. Kildishev, Advanced Optical Materials **9**, 2001469 (2021).

14. H. Matsubara, S. Yoshimoto, H. Saito, Y. Jianglin, Y. Tanaka, and S. Noda, Science **319**, 445 (2008).

15. D. D. Solnyshkov, T. Weiss, G. Malpuech, and N. A. Gippius, Applied Physics Letters **99**, 111110 (2011).

16. A. Grudinina, M. Efthymiou-Tsironi, V. Ardizzone, F. Riminucci, M. D. Giorgi, D. Trypogeorgos, K. Baldwin, L. Pfeiffer, D. Ballarini, D. Sanvitto, et al., Nature Communications **14**, 3464 (2023).

17. A. Gianfrate, H. Sigurdsson, V. Ardizzone, H. C. Nguyen, F. Riminucci, M. Efthymiou-Tsironi, K. W. Baldwin, L. N. Pfeiffer, D. Trypogeorgos, M. De Giorgi, et al., Nature Physics p. 1 (2024).

18. V. Ardizzone, F. Riminucci, S. Zanotti, A. Gianfrate, M. Efthymiou-Tsironi, D. Suàrez-Forero, F. Todisco, M. De Giorgi, D. Trypogeorgos, G. Gigli, et al., Nature **605**, 447 (2022).

19. F. Riminucci, V. Ardizzone, L. Francaviglia, M. Lorenzon, C. Stavrakas, S. Dhuey, A. Schwartzberg, S. Zanotti, D. Gerace, K. Baldwin, et al., Physical Review Applied **18**, 024039 (2022).

20. A. Imamoglu, R. J. Ram, S. Pau, Y. Yamamoto, et al., Physical Review A **53**, 4250 (1996).

21. A. Kavokin, J. J. Baumberg, G. Malpuech, and F. P. Laussy, *Microcavities* (Oxford university press, 2017).
22 L. Pitaevskii and S. Stringari, *Bose-Einstein Condensation* (Oxford Science Publications - International Series of Monographs on Physics 116, 2003).

23 D. Tanese, H. Flayac, D. D. Solnyshkov, A. Amo, A. Lemaître, E. Galopin, R. Braive, P. Senellart, I. Sagnes, G. Malpuech, et al., *Nature Communications* **4**, 1749 (2013).

24 N. Pernet, P. St-Jean, D. D. Solnyshkov, G. Malpuech, N. Carlon, Zambon, Q. Fontaine, B. Real, O. Jamadi, A. Lemaître, M. Morassi, et al., *Nature Physics* **18**, 678 (2022).

25 M. Milićević, O. Bleu, D. D. Solnyshkov, I. Sagnes, A. Lemaître, L. Le Gratiet, A. Harouri, J. Bloch, G. Malpuech, and A. Amo, *SciPost Physics* **5**, 064 (2018).

26 T. Jacqmin, I. Carusotto, I. Sagnes, M. Abbarchi, D. D. Solnyshkov, G. Malpuech, E. Galopin, A. Lemaître, J. Bloch, and A. Amo, *Physical Review Letters* **112**, 116402 (2014).

27 M. Zamfirescu, A. Kavokin, B. Gil, G. Malpuech, and M. Kaliteevski, *Phys. Rev. B* **65**, 161205 (2002), URL https://link.aps.org/doi/10.1103/PhysRevB.65.161205.

28 F. Li, L. Orosz, O. Kamoun, S. Bouchoule, C. Brimont, P. Disseix, T. Guillet, X. Lafosse, M. Leroux, J. Leymarie, et al., *Phys. Rev. Lett.* **110**, 196406 (2013), URL http://link.aps.org/doi/10.1103/PhysRevLett.110.196406.

29 O. Jamadi, F. Reveret, P. Disseix, F. Medard, J. Leymarie, A. Moreau, D. Solnyshkov, C. Deparis, M. Leroux, E. Cambril, et al., *Light: Science & Applications* **7**, 82 (2018).

30 G. Malpuech, A. Di Carlo, A. Kavokin, J. Baumberg, M. Zamfirescu, and P. Lugli, *Appl. Phys. Lett.* **81**, 412 (2002).

31 S. Christopoulos, G. B. H. von Högersthal, A. J. D. Grundy, P. G. Lagoudakis, A. V. Kavokin, J. J. Baumberg, G. Christmann, R. Butté, E. Feltin, J.-F. Carlin, et al., *Phys. Rev. Lett.* **98**, 126405 (2007), URL http://link.aps.org/abstract/PRL/v98/e126405.

32 H. Souissi, M. Gromovyi, T. Gueye, C. Brimont, L. Doyennette, D. Solnyshkov, G. Malpuech, E. Cambril, S. Bouchoule, B. Alloing, et al., *Phys. Rev. Applied* **18**, 044029 (2022), URL https://link.aps.org/doi/10.1103/PhysRevApplied.18.044029.

33 S. Kéna-Cohen and S. Forrest, *Nature Photonics* **4**, 371 (2010).

34 J. D. Plumhof, T. Stöferle, L. Mai, U. Scherf, and R. F. Mahrt, *Nature Materials* **13**, 247 (2014).

35 R. Su, C. Diederichs, J. Wang, T. C. Liew, J. Zhao, S. Liu, W. Xu, Z. Chen, and Q. Xiong, *Nano letters* **17**, 3982 (2017).

36 R. Su, S. Ghosh, J. Wang, S. Liu, C. Diederichs, T. C. Liew, and Q. Xiong, *Nature Physics* **16**, 301 (2020).

37 K. Peng, R. Tao, L. Haeberlé, Q. Li, D. Jin, G. R. Fleming, S. Kéna-Cohen, X. Zhang, and W. Bao, *Nature Communications* **13**, 7388 (2022).

38 D. D. Solnyshkov, H. Terças, K. Dini, and G. Malpuech, *Physical Review A* **89**, 033626 (2014).

39 C. Brimont, L. Doyennette, G. Kreyder, F. Réveret, P. Disseix, F. Médard, J. Leymarie, E. Cambril, S. Bouchoule, M. Gromovyi, et al., *Physical Review Applied* **14**, 054060 (2020).

40 L. Mallet-Dida, P. Disseix, F. Réveret, F. Médard, B. Alloing, J. Zúñiga-Pérez, and J. Leymarie, *New Journal of Physics* **24**, 033031 (2022).

41 I. Septembre, C. Leblanc, L. Hermet, H. S. Nguyen, X. Letartre, D. D. Solnyshkov, and G. Malpuech, *Physical Review B* **107**, 155304 (2023).

42 H. Souissi, M. Gromovyi, I. Septembre, V. Develay, C. Brimont, L. Doyennette, E. Cambril, S. Bouchoule, B. Alloing, E. Frayssinet, et al., arXiv preprint arXiv:2310.18661 (2023).

43 O. Jamadi, F. Réveret, E. Mallet, P. Disseix, F. Médard, M. Mihailovic, D. D. Solnyshkov, G. Malpuech, J. Leymarie, X. Lafosse, et al., *Physical Review B* **93**, 115205 (2016).

44 J. Ciers, D. D. Solnyshkov, G. Callsen, Y. Kuang, J.-F. Carlin, G. Malpuech, R. Butté, and N. Grandjean, *Physical Review B* **102**, 155304 (2020).

45 F. Riminucci, A. Gianfrate, D. Nigro, V. Ardizzone, S. Dhuey, L. Francaviglia, K. Baldwin, L. N. Pfeiffer, D. Ballarini, D. Trypogeorgos, et al., *Physical Review Letters* **131**, 246901 (2023).

46 D. Nigro and D. Gerace, *Physical Review B* **108**, 085305 (2023).

47 L. D. Landau and E. M. Lifshitz, *Quantum mechanics: non-relativistic theory*, vol. 3 (Elsevier, 2013).

48 D. D. Solnyshkov, I. Septembre, K. Ndiaye, and G. Malpuech, *Physical Review B* **108**, 144306 (2023).

49 A. Kavokin and G. Malpuech, *Cavity polaritons* (Elsevier, 2003).

50 E. Wertz, A. Amo, D. D. Solnyshkov, L. Ferrier, T. C. H. Liew, D. Sanvitto, P. Senellart, I. Sagnes, A. Lemaître, A. V. Kavokin, et al., *Physical Review Letters* **109**, 216404 (2012).

51 X. Wu, J. Song, S. Zhang, W. Du, Y. Wang, Z. Zhu, Y. Liang, Q. Zhang, Q. Xiong, and X. Liu, *arXiv preprint arXiv:2303.09923* (2023).

52 J. Kasprzak, M. Richard, S. Kundermann, A. Baas, P. Jeambrun, J. M. J. Keeling, F. Marchetti, M. Szymańska, R. André, J. Staehli, et al., *Nature* **443**, 409 (2006).

53 D. Dolinina and A. Yulin, *Physical Review E* **103**, 052207 (2021).

54 D. Dolinina and A. Yulin, *Physical Review E* **104**, 054214 (2021).

55 T. W. B. Kibble, *Journal of Physics A: Mathematical and General* **9**, 1387 (1976).

56 W. H. Zurek, *Nature* **317**, 505 (1985).

57 D. D. Solnyshkov, A. V. Nalitov, and G. Malpuech, *Physical Review Letters* **116**, 046402 (2016).

58 D. D. Solnyshkov, L. Bessonart, A. Nalitov, and G. Malpuech, *Physical Review B* **104**, 035423 (2021).

59 D. D. Solnyshkov, C. Leblanc, I. Septembre, and G. Malpuech, *Physical Review Letters* **129**, 066802 (2022).