[b,c]Shoji Hashimoto

# Composition of the inclusive semi-leptonic decay of $B$ meson

Gabriela Bailas

###### Abstract

Utilizing the approach recently proposed for the inclusive semi-leptonic decay
rate on the lattice, we compute the differential decay rate of a $B_{s}$ meson
for various kinematical channels. The results are compared with the
contributions from the ground states ($D$ and $D^{*}$) as well as from the
orbitally excited states ($D^{**}$’s). The computation so far is carried out
with an unphysically light bottom quark and strange spectator quark.

KEK-CP-0386

##  1 Introduction

Semi-leptonic decays of $B$ meson have been measured in various specific final
states, such as $D^{(*)}\ell\nu$ or $D^{**}\ell\nu$. They can be used to
determine the Cabibbo-Kobayashi-Maskawa (CKM) matrix element $|V_{cb}|$.
(Here, we focus on the $b\to c$ decays, but the formulation and calculation
method are applicable also for the $b\to u$ channels.) On the other hand, the
experimentalists can also perform the so-called inclusive analysis, i.e. all
possible final states including a charm quark are counted. Theoretically, such
experimental results may be compared with the OPE analysis [1, 2] to determine
$|V_{cb}|$.

More recently, one of the authors proposed a method to compute the inclusive
decay rate using lattice QCD [3]. It utilizes a method to implicitly sum over
all possible states to appear in two operator insertions on the lattice [4,
5]. To be explicit, we consider the forward-Compton amplitude of the form

$\langle
B(\bm{0})|\tilde{J}_{\mu}^{\dagger}(-\bm{q};t)\;\;\Big{|}\;\;\tilde{J}_{\nu}(\bm{q};0)|B(\bm{0})\rangle,$
(1)

where two flavor-changing currents $J_{\mu}$ are inserted with a specific
spatial momentum $\bm{q}$. The vertical line in the middle represents all
possible states of the specified quantum number that can contribute. The key
idea is to use the time separation $t$ to control the weight among the states
having different energies. The method allows a fully non-perturbative
computation of the inclusive decay rate, so that one can verify the OPE-based
calculation that has been used so far.

There is a long-standing tension between the $|V_{cb}|$ determinations from
the exclusive and inclusive channels. The aim of our study is to understand
the cause of the problem by providing a framework of theoretical calculation
that can be applied for both analyses. We may also identify the contribution
from the excited state $D$ mesons and $D\pi$ continuum states from the lattice
data, which would enable another test of the calculation by comparing with the
corresponding experimental data.

This contribution mainly describes the lattice computatin, and the comparison
with the OPE calculation for the same set of parameters is separately
presented in [6].

##  2 Inclusive decay rate: outline of the formalism

Here we outline the method proposed in [3].

The differential decay rate of the $B$ meson can be decomposed as
$d\Gamma\sim|V_{cb}|^{2}l^{\mu\nu}W_{\mu\nu}dq_{0}d\bm{q}$, where the leptonic
tensor $l^{\mu\nu}$ is determined by the kinematics of the final state leptons
$\ell$ and $\nu$. The momentum transfer to the lepton pair is
$q=(q_{0},\bm{q})$. The hadronic tensor $W_{\mu\nu}$, on the other hand, has a
complicated structure

$W_{\mu\nu}(q_{0},\bm{q})\sim\sum_{X}(2\pi)^{3}\delta^{(4)}(p_{B}-q-p_{X})\frac{1}{2m_{B}}\langle
B(\bm{p}_{B})|J_{\mu}^{\dagger}(0)|X\rangle\langle
X|J_{\nu}(0)|B(\bm{p}_{B})\rangle,$ (2)

where the flavor-changing current $J_{\mu}$ induces the decay $b\to c$, and
the sum of the state $X$ runs over all possible final states including a charm
quark. When the initial $B$ meson is at rest, $\bm{p}_{B}=\bm{0}$, the final
hadronic state $X$ has a momentum $\bm{p}_{X}=-\bm{q}$ given by the momentum
transfer $q$. Since it depends on the structure of the initial $B$ meson
state, the hadronic tensor $W_{\mu\nu}$ is also called the structure function.

We notice that the sum over the states $X$ may be considered as an integral
over its energy $p_{X}^{0}$. Because of the $\delta$-function, it is given by
$p_{X}^{0}\equiv\omega=m_{B}-q_{0}$. Thus, the structure function picks up the
state of energy $\omega$, and

$W_{\mu\nu}(\omega,\bm{q})\sim\langle
B(\bm{0})|\tilde{J}_{\mu}^{\dagger}(-\bm{q})\delta(\omega-\hat{H})\tilde{J}_{\nu}(\bm{q})|B(\bm{0})\rangle.$
(3)

Here $\tilde{J}_{\mu}(\bm{q})$ is a Fourier-transform of the current. We
specify the energy $\omega$ by inserting the $\delta$-function
$\delta(\omega-\hat{H})$ between the currents. ($\hat{H}$ is the Hamiltonian
of QCD.)

The total decay rate can then be written as

$\Gamma\propto|V_{cb}|^{2}\int_{0}^{\bm{q}^{2}_{\mathrm{max}}}d\bm{q}^{2}\int_{\sqrt{m_{D}^{2}+\bm{q}^{2}}}^{m_{B}-\sqrt{\bm{q}^{2}}}d\omega\,K(\omega;\bm{q}^{2})\langle
B(\bm{0})|\tilde{J}_{\mu}^{\dagger}(-\bm{q})\delta(\omega-\hat{H})\tilde{J}_{\nu}(\bm{q})|B(\bm{0})\rangle,$
(4)

where $K(\omega;\bm{q}^{2})$ is a kinematical factor originating from the
leptonic tensor.

The problem is then how to compute the integral over the final-state energy
$\omega$ with the weight of $K(\omega;\bm{q})$. The upper limit of the
$\omega$-integral is also imposed by the kinematics; we may include them in
$K(\omega,\bm{q}^{2})$ using the Heaviside function as
$\theta(m_{B}-\sqrt{\bm{q}^{2}}-\omega)$ while extending the upper limit to
infinity. Then, we can rewrite the $\omega$-integral using

$\int_{0}^{\infty}\\!d\omega\,K(\omega;\bm{q}^{2})\langle
B(\bm{0})|\tilde{J}_{\mu}^{\dagger}(-\bm{q})\delta(\omega-\hat{H})\tilde{J}_{\nu}(\bm{q})|B(\bm{0})\rangle=\langle
B(\bm{0})|\tilde{J}_{\mu}^{\dagger}(-\bm{q})K(\hat{H};\bm{q}^{2})\tilde{J}_{\nu}(\bm{q})|B(\bm{0})\rangle.$
(5)

Here, the kinematical factor is promoted to an operator by replacing $\omega$
by $\hat{H}$. On the other hand, what one can calculate on the lattice is the
Compton amplitude (1), which is also written in the form

$\langle
B(\bm{0})|\tilde{J}_{\mu}^{\dagger}(-\bm{q};t)\tilde{J}_{\nu}(\bm{q};0)|B(\bm{0})\rangle=\langle
B(\bm{0})|\tilde{J}_{\mu}^{\dagger}(-\bm{q};t)e^{-\hat{H}t}\tilde{J}_{\nu}(\bm{q};0)|B(\bm{0})\rangle,$
(6)

because the time separation between the currents may be described by the
transfer matrix $\exp(-\hat{H}t)$.

Given the similarity between (5) and (6), one notices that (5) can be
evaluated if the kernel operator $K(\hat{H})$ can be approximated in the form
$K(\hat{H})\simeq
k_{0}+k_{1}e^{-\hat{H}}+k_{2}e^{-2\hat{H}}+\cdots+k_{N}e^{-N\hat{H}}$ with
some coefficients $k_{j}$. Such approximation can be constructed using the
Chebyshev polynomials:

$K(\omega)\simeq\frac{c_{0}}{2}+\sum_{j=1}^{N}c_{j}^{*}T_{j}^{*}(e^{-\omega}),$
(7)

where $T_{j}^{*}(x)$’s are shifted Chebyshev polynomials defined as
$T_{j}^{*}(x)=T_{j}(2x-1)$ from the standard Chebyshev polynomial $T_{j}(z)$.
Thus, the shifted Chebyshev polynomials are defined for $0\leq x\leq 1$ or
$\infty\geq\omega\geq 1$ since $x=e^{-\omega}$. The coefficients $c_{j}^{*}$
can be easily computed for arbitrary kernel operator $K(\omega)$. Since the
Chebyshev polynomials $T_{j}^{*}(e^{-\omega})$ are constructed from
polynomials $(e^{-\omega})^{m}$ with positive integer $m$’s, they are related
to the transfer matrix $e^{-\hat{H}t}$ appearing in (6). The forward Compton
amplitude (6) can therefore be used to approximate the target quantity (5).

The kernel function $K(\omega;\bm{q}^{2})$ has the following structure

$K(\omega;\bm{q}^{2})\sim e^{2\omega
t_{0}}(m_{B}-\omega)^{l}\theta(m_{B}-\sqrt{\bm{q}^{2}}-\omega).$ (8)

Here, $(m_{B}-\omega)^{l}$ ($l=$ 0, 1 or 2) originates from the leptonic
tensor; the factor $e^{2\omega t_{0}}$ is introduced to avoid any divergence
due to a contact term between the two currents by normalizing the matrix
element by the value at a small time separation $t_{0}$. The Heaviside
$\theta$-function $\theta(m_{B}-\sqrt{\bm{q}^{2}}-\omega)$ implements the
upper limit of the $\omega$ integral (see (4)).

The Chebyshev approximation of the kernel function is harder when the function
contains a discontinuous (or even rapid) change such as that given by the
$\theta$-function, and we introduce a smearing to modify the $\theta$-function
to a smooth function with a smearing width $\sigma$. (To be explicit we use
the sigmoid function, but the details are not important.) We need to take the
limit $\sigma\to 0$ to obtain the final result.

Figure 1: Kernel function for $l=0$ with smearing with the width $\sigma=0.05$
(left panel) 0.1 (middle) and 0.2 (right). Together with the true function
(thich curve), Chebyshev approximations of order $N$ = 5, 10 and 20 are shown
(dashed and dotted curves). The threshold is at $\omega\simeq$ 1.0. (All in
the lattice unit.)

The kernel function $K(\omega)$ is plotted for $l=0$ with the smearing width
$\sigma$ = 0.05, 0.1 and 0.2 in Fig. 1. Around the threshold of the
$\theta$-function, the shape is smeared so that the function is smoothed. The
Chebyshev approximations are also shown in the plots for the order of
polynomial $N$ = 5, 10 and 20. The approximation is rather precise when the
width is large ($\sigma$ = 0.2 on the right panel) even with the limited order
of polynomials, while it gets harder for small $\sigma$. Although it would
depend on the precision one wants to achieve, it seems that at least $N$ = 20
is necessary to achieve a sensible approximation for $\sigma$ = 0.05.

##  3 Compton amplitude

We compute the Compton amplitude of the form (6) on the lattice. It is
obtained from a four-point function with two flavor-changing $V-A$ currents
inserted with time separation $t$ between interpolating operators to create or
annihilate the $B$ meson. To compute the amplitude with all possible
combinations of spin orientations $\mu$ and $\nu$ we have to repeat the
computation of sequential sources many times, as well as for the choices of
final momenta. (The initial $B$ meson is set at rest.)

The computation is done on a lattice of $48^{3}\times 96$ at $1/a\simeq$ 3.6
GeV generated including up, down and strange sea quarks described by the
Mobius domain-wall fermions. The ensemble is the same as that used in [3] and
is a part of the large set of ensembles generated and used in [7]. The valence
quarks are also domain-wall fermions; the charm quark mass is tuned to its
physical value while the bottom quark mass is taken (unphysically) light and
about 2.7 GeV. The spectator quark is strange in this study, so that the
initial state is (unphysically light) $B_{s}$ meson.

The statistics is 100 gauge configurations and the measurement is repeated
four times on each configuration with different source time slices.

Figure 2: Compton amplitude as a function of the time separation. The final
state has a vanishing spatial momentum. Different current insertions are shown
together: vector ($V$) and axial-vector ($A$) currents in the temporal ($0$)
and spatial ($k=1$) directions.

Fig. 2 shows the Compton amplitude (6) as a function of the time separation
$t\equiv t_{2}-t_{1}$ for the case of zero spatial momentum in the final
state. The current insertions are either $V(t_{2})V(t_{1})$ or
$A(t_{2})A(t_{1})$ combinations of vector ($V$) or axial-vector ($A$)
currents. Non-zero values are obtained for $(\mu,\nu)=(0,0)$ (temporal) or
$(\mu,\nu)=(k,k)$ (spatial) orientations of the currents. Among them, the
$V_{0}V_{0}$ and $A_{k}A_{k}$ combinations show the largest signal due to the
$S$-wave ground-state contributions from the pseudo-scalar $D_{s}$ ($0^{-}$)
or vector $D_{s}^{*}$ ($1^{-}$) meson. If we look into the details, the
amplitude of $V_{0}V_{0}$ is slightly larger than $A_{k}A_{k}$ because the
pseudo-scalar state is lighter.

The other channels $A_{0}A_{0}$ and $V_{k}V_{k}$ show substantially ($\times
50$) smaller contributions. They correspond to the scalar ($0^{+}$) and axial-
vector ($1^{+}$) states, respectively. In the quark model, they are $P$-wave
states and their coupling to the initial $B_{s}$ meson is weak. According to
the analysis based on the heavy quark effective theory (HQET) [8], the $B\to
D^{**}\ell\nu$ ($D^{**}$ denotes the $P$-wave states generically) form factor
is suppressed as $\bar{\Lambda}/2m_{c}$ with $\bar{\Lambda}$ a typical QCD
scale $\sim$ 300 MeV compared to the $B\to D^{(*)}\ell\nu$ form factors, which
is $O(1)$, so that the Compton amplitude is suppressed by
$(\bar{\Lambda}/2m_{c})^{2}\sim$ 0.01.

These opposite parity channels can also couple to $DK$ continuum, which is
heavier than the above mentioned $P$-wave mesons. Interestingly, the data show
the contributions of such excited states at small time separations
$t_{2}-t_{1}\lesssim$ 5.

##  4 Inclusive decay rate

According to (4), we obtain the total decay rate from the Compton amplitudes
as described in the previous section. The kinematical factor
$K(\omega;\bm{q}^{2})$ is promoted to the integral kernel for the
$\omega$-integral. The spatial momentum integration over $\bm{q}^{2}$ is yet
to be performed.

Figure 3:  Differential decay rate divided by $|\bm{q}|$ as a function of
$\bm{q}^{2}$. The results are shown for four different channels:
$VV_{\parallel}$, $VV_{\perp}$, $AA_{\parallel}$ and $AA_{\perp}$. (For more
details, see the text.) The bands are contributions of ground state, $D$ or
$D^{*}$ meson, as estimated using the corresponding form factors.

The results for the differential decay rate (divided by $|\bm{q}|$) are
plotted in Fig. 3 as a function of $\bm{q}^{2}$. Four contributions are
plotted separately: $VV$ and $AA$ denote the two currents inserted, which are
decomposed from the $V-A$ weak currents. The other possible combinations, $VA$
and $AV$, do not contribute for the differential decay rate after integrating
over the lepton energy $E_{\ell}$. The current orientations, $\parallel$ and
$\perp$, distinguish the currents in the direction parallel or perpendicular
to the momentum $\bm{q}$. For instance, when the momentum $\bm{q}$ is in the
$z$-direction, $(0,0,1)$, the $VV\\!\perp$ represent the contribution from
$V^{1}V^{1}$ and $V^{2}V^{2}$, while $VV\\!\parallel$ is that of $V^{3}V^{3}$,
$V^{3}V^{0}$, $V^{0}V^{3}$ and $V^{0}V^{0}$.

Also shown by the bands are the corresponding contributions from the ground
states, i.e. from $B_{s}\to D_{s}\ell\nu$ for $VV\\!\parallel$ and from
$B_{s}\to D_{s}^{*}\ell\nu$ for others. We have separately computed the form
factors of these decay modes, $h_{+}(w)$, $h_{-}(w)$ and $h_{A_{1}}(w)$,
$h_{A_{2}}(w)$, $h_{A_{3}}(w)$, $h_{V}(w)$. (The computation is done along
with [9], and the analysis method is essentially the same.) Using these
results, the differential decay rate for each channel can be constructed.

As one can see, the inclusive decay rate is dominated by the ground state
contributions for each channel. This is not unreasonable because the initial
state is lighter than the physical $B_{s}$ meson. Kinematically there is not
so much room to produce extra particles other than $D_{s}^{(*)}$. Also, the
heavy quark symmetry implies that the wave function of light degrees of
freedom are very similar between the initial $B_{s}$ and final $D_{s}^{(*)}$
states, so that the overlap between these states are enhanced when the initial
and final quark masses are similar, especially near the zero-recoil
($\bm{q}\sim 0$) limit.

As we mentioned in the previous section, the contribution from the excited
states including those from the $P$-wave meson is visible in the Compton
amplitudes. It is however not significant for the differential decay rate,
since their contributions are roughly 50 times smaller.

##  5 Summary

Framework to compute inclusive decay rate on the lattice is now available. The
essential step is the reconstruction of the energy integral from Euclidean
lattice correlators corresponding to the Compton amplitudes. In our case, it
is implemented using the Chebyshev approximation. Other methods, such as the
Backus-Gilbert method can also do the job [10, 11]. The method to compute the
energy integral can be also applied for calculations of (not so) deep
inelastic scattering cross section, $\sigma(\ell N\to\ell^{\prime}X)$ [12] as
well as those of spectral sum of hadronic vacuum polarization [13].

The Compton amplitudes contain the contributions from ground states and
excited states. With the lattice parameters taken in this work the ground
state contributions dominate the differential cross section, although the
excited state contributions such as those from the $P$-wave states are visible
in the Compton amplitude.

The present results can already been used to compare with the analytic OPE
calculations. Even though the $b$-quark mass is not tuned to the physical
value, the comparison would still yield useful tests of the theoretical
approaches [6].

## Acknowledgement

This work is a part of the efforts to make a comparison between lattice and
OPE calculations of inclusive semileptonic $B$ decays, being done in
collaboration with Paolo Gambino and Sandro Machler. We thank the members of
the JLQCD collaboration for helpful discussions and providing the
computational framework and lattice data. Takashi Kaneko, in particular,
provided the unpublished data of the $B_{s}\to D_{s}^{(*)}$ form factors
computed on the same ensemble. This work is supported in part by JSPS KAKENHI
grant number 18H03710 and by the PostK and Fugaku supercomputer project
through the Joint Institute for Computational Fundamental Science (JICFuS).

## References

  * [1] B. Blok, L. Koyrakh, M. A. Shifman and A. I. Vainshtein, Phys. Rev. D 49, 3356 (1994) [erratum: Phys. Rev. D 50, 3572 (1994)] doi:10.1103/PhysRevD.50.3572 [arXiv:hep-ph/9307247 [hep-ph]].
  * [2] A. V. Manohar and M. B. Wise, Phys. Rev. D 49, 1310-1329 (1994) doi:10.1103/PhysRevD.49.1310 [arXiv:hep-ph/9308246 [hep-ph]].
  * [3] P. Gambino and S. Hashimoto, Phys. Rev. Lett. 125, no.3, 032001 (2020) doi:10.1103/PhysRevLett.125.032001 [arXiv:2005.13730 [hep-lat]].
  * [4] S. Hashimoto, PTEP 2017, no.5, 053B03 (2017) doi:10.1093/ptep/ptx052 [arXiv:1703.01881 [hep-lat]].
  * [5] G. Bailas, S. Hashimoto and T. Ishikawa, PTEP 2020, no.4, 043B07 (2020) doi:10.1093/ptep/ptaa044 [arXiv:2001.11779 [hep-lat]].
  * [6] P. Gambino, S. Hashimoto and S. Mächler, [arXiv:2111.02833 [hep-ph]]; PoS(LATTICE2021)512.
  * [7] K. Nakayama, B. Fahy and S. Hashimoto, Phys. Rev. D 94, no.5, 054507 (2016) doi:10.1103/PhysRevD.94.054507 [arXiv:1606.01002 [hep-lat]].
  * [8] A. K. Leibovich, Z. Ligeti, I. W. Stewart and M. B. Wise, Phys. Rev. D 57, 308-330 (1998) doi:10.1103/PhysRevD.57.308 [arXiv:hep-ph/9705467 [hep-ph]].
  * [9] T. Kaneko et al., PoS(LATTICE2021)561.
  * [10] M. T. Hansen, H. B. Meyer and D. Robaina, Phys. Rev. D 96, no.9, 094513 (2017) doi:10.1103/PhysRevD.96.094513 [arXiv:1704.08993 [hep-lat]].
  * [11] M. Hansen, A. Lupo and N. Tantalo, Phys. Rev. D 99, no.9, 094508 (2019) doi:10.1103/PhysRevD.99.094508 [arXiv:1903.06476 [hep-lat]].
  * [12] H. Fukaya, S. Hashimoto, T. Kaneko and H. Ohki, Phys. Rev. D 102, no.11, 114516 (2020) doi:10.1103/PhysRevD.102.114516 [arXiv:2010.01253 [hep-lat]].
  * [13] T. Ishikawa and S. Hashimoto, Phys. Rev. D 104, no.7, 074521 (2021) doi:10.1103/PhysRevD.104.074521 [arXiv:2103.06539 [hep-lat]].


# Dirac exciton-polariton condensates in photonic crystal gratings

Helgi Sigurðsson<EMAIL_ADDRESS>Institute of Experimental
Physics, Faculty of Physics, University of Warsaw, ul. Pasteura 5, PL-02-093
Warsaw, Poland  Science Institute, University of Iceland, Dunhagi 3, IS-107
Reykjavik, Iceland  Hai Chau Nguyen  Naturwissenschaftlich-Technische
Fakultät, Universität Siegen, Walter-Flex-Straße 3, 57068 Siegen, Germany  Hai
Son Nguyen  Univ Lyon, Ecole Centrale de Lyon, INSA Lyon, Université Claude
Bernard Lyon 1, CPE Lyon, CNRS, INL, UMR5270, Ecully 69130, France  Institut
Universitaire de France (IUF), 75231 Paris, France

(August 28, 2024)

###### Abstract

Bound states in the continuum have recently been utilized in photonic crystal
gratings to achieve strong coupling and ultralow power-driven condensation of
bosonic exciton-polariton quasiparticles with atypical Dirac-like features in
their dispersion relation. Here, we develop the single- and many-body theory
of these new effective relativistic exciton-polaritons modes and describe
their mean field condensation dynamics facilitated by the interplay between
protection from the radiative continuum and negative-mass pump induced optical
trapping. Our theory accounts for many tunable grating parameters giving full
control over the diffractive coupling properties between guided polaritons and
the radiative continuum previously unexplored in the context of driven
condensation. In particular, we discover stable cyclical condensate solutions
mimicking a driven-dissipative analog of the zitterbewegung effect
characterized by coherent superposition of both ballistic (rapid phase front)
and trapped (slow phase front) polariton waves. Finally, important
distinctions are drawn between the concepts of near field and far field in the
photonic grating, clarifying recent experimental observations on the emission
characteristics of these long lived nonlinear Dirac polaritons.

##  I Introduction

Quantum fluids with both unusual dispersive properties and strong non-
Hermitian effects form a new exciting testbed to investigate many-body
physics. Macroscopic quantum fluids of exciton-polaritons [1] (hereafter,
polaritons) are particularly suited for this task given their optical
malleability [2, 3, 4, 5] and readout, high interaction strengths, and the
wide breadth of materials permitting polariton condensation at elevated
temperatures [6, 7]. These light-mass bosonic quasiparticles form in the
strong coupling regime between confined photonic modes and exciton resonances
in semiconductor microcavities [1]. In particular, there has been some
excitement in simulating relativistic phenomena in artificial Dirac materials
exploiting the polariton spin in both real and synthetic magnetic fields [8,
9, 10, 11, 12, 13, 14] and patterned photonic structures [15, 16, 17, 18, 19,
20, 21, 22]. Such materials, hosting associated Dirac cones, offer valuable
insight to a plethora of exotic phenomena such as quantum Hall physics [9],
nontrivial topological phases [18], Weyl semimetals [23], and relativistic
trapping [24, 25] while supplemented with strong polariton nonlinearities.
Moreover, alternative neighbouring platforms for exploration into light-matter
Dirac physics involve phonon-polaritons [26, 27] and plasmon-polaritons [28].

Recently, a photonic crystal platform was realized to explore Dirac physics
using exciton-polaritons. It consists of a subwavelength grated GaAs-based
semiconductor waveguide embedded with multiple quantum wells hosting Wannier-
Mott excitons (see Fig. 1). Both strong light-matter coupling and ultralow
threshold polariton Bose-Einstein condensation into the waveguide’s associated
bound-states-in-the-continuum (BIC) was demonstrated in the same study [29].
Strong coupling is facilitated by the photonic structure’s protection from the
continuum [30, 31] which allows photons to survive long enough to form
polariton states [32, 33]. The initial experiment [29] was soon followed with
fascinating results on the behaviour of the fluid’s elementary excitations
[34] and demonstration of macroscopic hybridization between coupled BIC
condensates [35]. Besides III-V semiconductor photonic crystals [36], other
platforms able to show BIC-facilitated strong exciton-photon coupling consist
of through dry transfer deposited transition metal dichalcogenide monolayers
such as MoSe2 [37, 38] or WS2 [39, 40, 41]; or spin coated hybrid organic-
inorganic perovskites [42, 43, 44, 45].

Inspired by these rapid developments in BIC facilitated condensation, and the
surging interest to simulate nonlinear and non-Hermitian relativistic physics
in an optically addressable setting [46], we develop the two-band theory of
BIC Dirac exciton-polaritons in photonic crystal gratings. We then propose a
many-body description in the mean field picture, allowing us to construct a
simplified generalized two-band Gross-Pitaevskii model describing the
interplay between BIC facilitated condensation and negative-mass optical
trapping of Dirac polaritons. Our theory is in excellent agreement with recent
experiments on polariton crystal gratings under nonresonant excitation [29]
and goes beyond by uncovering fascinating cyclical condensate dynamics through
continuous adjustment of experimentally accessible parameters.

In particular, our mean field simulations reveal that the condensate negative-
mass population “drops” iteratively to higher order optical trap modes as a
function of pump power density, underlining complex interplay between pump
induced polariton energy shifts and gain. Around the drop point the condensate
can converge into a limit cycle solution, characterized by a coherent
superposition of two distinct trap levels. More interestingly, through careful
tuning of the photonic grating the BIC can be gradually moved from the lower
negative-mass branch to the upper positive-mass branch, until the condensate
suddenly stabilizes into a zitterbewegung-like [47, 14] limit cycle, forming a
strange mixture of confined low momentum negative-mass polaritons and
ballistic high momentum positive-mass polaritons. Lastly, we elucidate on how
the polariton field within the photonic crystal grating relates to the emitted
far field measured in experiment, in sharp contrast to typical condensation
experiments in planar microcavities [48].

##  II Massive Dirac polaritons

###  II.1 Photonic modes in optical grating

We consider a wave guide stack consisting of multiple periodic layers of GaAs
quantum wells separated by AlGaAs barriers along the normal $z$-direction [see
Fig. 1(a)]. On its upper surface, the wave guide is patterned with a one-
dimensional subwavelength grating with a period $a$ along the $x$-direction
and filling factor $FF$, forming a one-dimensional (1D) photonic crystal slab.
This kind of a sample has been demonstrated in [29, 35], however our developed
theory is also suitable for photonic crystals deposited with 2D materials [49,
39, 37, 38, 40, 41] or perovskites [42, 43, 44, 45].

From here on we will consider only propagation of electromagnetic modes along
the $x$-direction. We concentrate also on the TE mode, where the polarization
of the electromagnetic field is assumed to be along the $y$-axis. In the
absence of the grating, the guided electromagnetic modes are plane waves
$\lvert e^{iqx}\rangle$ with momentum $q$ and frequency $\omega_{q}$. As is
well known [50], the grating with period $a$ can be modeled by a periodic
potential acting on the photonic modes,

$V_{\mathrm{ph}}(x,z)=u(x)w(z)$ (1)

where $u(x)$ is a periodic function of period $a$, and $w(z)$ is a square-step
function that is equal to $1$ for $z$ within the patterned layers and vanishes
otherwise.

The periodic potential couples photonic modes with wavevectors different by an
integer number of the primitive reciprocal lattice number $K=2\pi/a$, a
condition known as Bragg reflection which folds photonic bands over the first
Brillouin zone. Here the period $a$ is chosen such that the exciton energy is
around the frequencies $\omega_{\pm K}$. Therefore, the relevant photonic
modes for exciton-photon coupling correspond to those with wavevectors $q=\pm
K+k$ with $k\ll K$. In this range of frequencies, the periodic potential
couples nearly degenerate guided modes $\lvert e^{i(K+k)x}\rangle$ and $\lvert
e^{i(-K+k)x}\rangle$, giving rise to the Dirac Hamiltonian describing the
dynamics of the electromagnetic waves. Remarkably, the periodic potential not
only couples the two modes $\lvert e^{\pm iKx}\rangle$ together, but also
couples them with lossy modes residing at normal incidence $q=0$. This renders
the two guided modes $\lvert e^{\pm iKx}\rangle$ eventually lossy; see Fig.
1(b). In total, the dynamics of the electromagnetic waves in the photonic
crystal grating can be described by a lossy Dirac Hamiltonian,

Figure 1: (a) Schematic of the 1D subwavelength grated waveguide with period
$a$ and filling fraction $FF$. The different colored layers denote e.g. GaAs
quantum wells separated by AlGaAs barriers. (b) Sketch of the two guided modes
$\lvert e^{\pm iKx}\rangle$ coupled with each other and coupled to the Fabry-
Pérot lossy modes through the grating. (c,d) Dashed black lines show the real
and imaginary parts of the dispersion relation of uncoupled counter-
propagating photons. Blue and red solid curves correspond to finite
diffractive coupling leading to gap opening and formation of symmetric and
antisymmetric standing wave modes, respectively. The BIC can be seen from the
absence of losses in the antisymmetric mode at $k=0$. We have used:
$\varphi=0$, $\hbar v=56.6$ meV $\mu$m, $\hbar\kappa=1.75$ meV, and
$\hbar\gamma=0.18$ meV corresponding to physical values obtained from rigorous
coupled-wave analysis (RCWA) (see the section C.1 of the Appendix).

$H_{\text{ph}}=\begin{pmatrix}vk&\kappa\\\
\kappa&-vk\end{pmatrix}-i\gamma\begin{pmatrix}1&e^{i\varphi}\\\
e^{-i\varphi}&1\end{pmatrix},$ (2)

with

$\kappa=|U_{2}|,\quad\varphi=2\operatorname{arg}(U_{1})-\operatorname{arg}(U_{2}),\quad\gamma\propto|U_{1}|^{2},$
(3)

where $U_{p}$ is $p$th Fourier coefficient of $u(x)$, i.e., $U_{p}=\langle
1\rvert u(x)\lvert e^{-ipKx}\rangle$. We have also linearised the dispersal
relation around $K$ so that $\omega_{k\pm K}=\omega_{K}\pm vk$ and ignore the
constant term. A more detailed derivation of this effective low-momentum
photonic Hamiltonian is given in Appendix A.

Equation (2) is a non-Hermitian Dirac Hamiltonian of which the coupling
between counter-propagating guided photon modes depends on the waveguide
diffraction mechanism of strength $\kappa$ and the loss exchange mechanism of
strength $\gamma e^{i\varphi}$ via the radiative continuum [32, 51]. Notice
that we adopt here the convention of using $e^{-i\omega t}$ to describe the
temporal oscillation, hence losses are given by the negative imaginary
component of the dispersion relation.

Following Eq. (3), the parameters $\kappa$, $\gamma$ and $\varphi$ are
dictated by the two first Fourier components $U_{1}$ and $U_{2}$ of the
periodic modulation $u(x)$. The diffractive coupling $\kappa>\gamma$ is the
main parameter responsible for the bandgap opening at $k=0$ [see solid curves
in Fig. 1(c)]. Its value can be engineered by tuning the filling fraction of
the grating [32]; see Appendix C.3 for numerical values of $\kappa$ in
realistic sample designs.

The dispersion relation of the Hamiltonian (2), corresponding to the new
symmetric ($+$) and antisymmetric ($-$) standing-wave eigenmodes, can be found
as

$\omega_{\pm}=-i\gamma\pm\sqrt{(vk)^{2}+\kappa^{2}-\gamma^{2}-2i\gamma\kappa\cos{(\varphi)}}.$
(4)

In the absence of losses one has recover the Dirac dispersal relation with the
effective Dirac mass $\kappa$ of the photons. The real and imaginary parts of
the energies are plotted in Fig. 1(c,d), respectively, for $\varphi=0$,
showing clearly a zero-loss BIC at $k=0$ in the antisymmetric branch. This can
be understood by explicitly compute the imaginary part of the frequency,

$\displaystyle\text{Im}{(\omega_{\pm})}$
$\displaystyle=-\gamma\mp\sin{\left[\frac{1}{2}\tan^{-1}{\left(\frac{2\kappa\gamma\cos{(\varphi)}}{(vk)^{2}+\kappa^{2}-\gamma^{2}}\right)}\right]}$
$\displaystyle\times\sqrt[4]{\left[(vk)^{2}+\kappa^{2}-\gamma^{2}\right]^{2}+4\kappa^{2}\gamma^{2}\cos^{2}{(\varphi)}}.$
(5)

For small loss, $\gamma\ll\kappa$, and $k\to 0$, we have

$\lim_{k\to 0}\text{Im}{(\omega_{\pm})}\approx-\gamma(1\pm\cos{\varphi}).$ (6)

If the grating design possesses mirror symmetry, $x\rightarrow-x$, we have
$u(x)=u(-x)$ and all Fourier coefficients of $u(x)$ are real. Hence, $\varphi$
can only take values that are integer multiple of $\pi$. Notably, a $\pi$-jump
of $\varphi$ can be obtained by sweeping the filling fraction through a band-
inversion point; see Appendix C.3. From (6), in the case of $\varphi=0$, a BIC
mode of infinite lifetime appears in the center of the lower $\omega_{-}$
branch while a lossy $2\gamma$ mode appears in the upper $\omega_{+}$ branch
[see Fig. 1(d)].

$\lim_{k\to
0}\text{Im}{(\omega_{\pm})}\Big{|}_{\varphi=0}=\begin{cases}-2\gamma\\\
0\end{cases}.$ (7)

Conversely, when $\varphi=\pi$ the BIC switches branches. In both case, as
long as $\varphi$ is an integer multiple of $\pi$ due to the presence of the
mirror symmetry $x\rightarrow-x$, a formation of BIC at $\Gamma$ point is
guaranteed. This BIC is therefore of a symmetry-protected nature.

Breaking the mirror symmetry will relax the aforementioned constraint on
$\varphi$. As a result, none of the branches exhibit a BIC. If the symmetry
breaking is only a small perturbation, the symmetry-protected BIC becomes a
quasi-BIC with finite but extremely long lifetime [52, 53, 54, 55, 56]. We
note that while a similar form of (2) has been previously reported to describe
the formation of symmetry-protected BICs [32, 29], our work provide the first
effective Hamiltonian for the general case where the in-plane mirror-symmetry
can be broken. Importantly, for realistic grating structures with lateral
symmetry breaking design, a fine tuning of $\varphi$ from 0 to $\pi/2$ can be
achieved (see Appendix C.4).

To start with, we will explore the case of grating with lateral mirror
symmetry that has $\varphi=0$. At a later stage we will relax this constraint
and explore Dirac-polariton BIC condensation when $0<\varphi<\pi$.

###  II.2 Coupling between photonic BIC and excitons

We now consider the strong light-matter coupling regime between the photons
and quantum well excitons leading to new hybrid modes known as exciton-
polaritons [1]. Our goal is to design a simple mean field model describing the
dynamics of a driven Bose-Einstein condensate of Dirac-polaritons. We start in
the single particle limit where a standard coupled oscillator model can be
written to describe the mixing of standing-wave photons and excitons with a
light-matter coupling parameter $\Omega$ (also known as the exciton-photon
Rabi frequency) [57, 32],

$H_{\text{pl}}=\begin{pmatrix}\omega_{+}&\Omega&0&0\\\
\Omega&\omega_{X}&0&0\\\ 0&0&\omega_{-}&\Omega\\\ 0&0&\Omega&\omega_{X}\\\
\end{pmatrix}.$ (8)

Here, $\omega_{X}=\omega^{(0)}_{X}-i\gamma_{\text{nr}}$ where
$\omega^{(0)}_{X}$ and $\gamma_{\text{nr}}^{-1}$ denote the detuning of the
excitons from the photon branches at $k=0$ and their nonradiative lifetime,
respectively. Here we have assumed that the mass of excitons is practically
infinite compared to the confined photons. We note that there is no direct
coupling from the excitons to the radiative continuum, only to the localized
waveguided modes. The eigenmodes of the above matrix are referred to as upper
$|U,\pm\rangle$ and lower $|L,\pm\rangle$ symmetric-antisymmetric polaritons
with a dispersion relation,

$\displaystyle\begin{split}\omega_{U,\pm}&=\frac{\omega_{\pm}+\omega_{X}}{2}+\frac{1}{2}\sqrt{(\omega_{\pm}-\omega_{X})^{2}+4\Omega^{2}},\\\
\omega_{L,\pm}&=\frac{\omega_{\pm}+\omega_{X}}{2}-\frac{1}{2}\sqrt{(\omega_{\pm}-\omega_{X})^{2}+4\Omega^{2}},\end{split}$
(9)

which is plotted in Fig. 2(a) and 2(b). Notice that due to the exciton losses,
the BIC now becomes a quasi-BIC with finite losses.

Figure 2: Solid lines show the (a) real and (b) imaginary energies of the
polariton dispersion from Eq. (9) with the exciton line and the lower
symmetric $\omega_{L,+}$ and antisymmetric $\omega_{L,-}$ branches marked.
(c,d) Show a zoom of the real and imaginary energies belonging to the lower
symmetric and antisymmetric polariton states. The dashed black lines show the
approximation obtained using (10). Here we set: $\varphi=0$,
$\Omega/\kappa=1.8$, $\hbar\omega_{X}^{(0)}/\kappa=5.5$, and
$\gamma_{\text{nr}}=\gamma$ corresponding to physical values obtained from
RCWA (see the section C.1 of the Appendix).

We will assume that the upper polariton branches $\omega_{U,\pm}$ are far away
in energy and weakly populated and thus only focus on the lower polariton
branches $\omega_{L,\pm}$ around $k=0$ where condensation preferentially takes
place [29, 58, 35]. The lower branches can be approximated by considering
first the coupling of forward- $\lvert e^{iKx}\rangle$ and backward-
propagating $\lvert e^{-iKx}\rangle$ photons to excitons, leading to lower-
forward $\lvert L,e^{iKx}\rangle$ and lower-backward $\lvert
L,e^{-iKx}\rangle$ propagating polaritons. This is followed up by the photonic
diffractive coupling mechanism evaluated at small momenta. The lower polariton
dispersion can then be written,

$\omega_{L,\pm}\approx\tilde{\omega}_{L}\pm\sqrt{(\tilde{v}k)^{2}+[\kappa^{2}-\gamma^{2}-2i\gamma\kappa\cos{(\varphi)}]|C_{0}|^{4}},$
(10)

where the first term in (10) corresponds simply to an overall complex energy
shift due to the light-matter coupling which is written,

$\tilde{\omega}_{L}=\frac{\omega_{X}}{2}\left(1-\sqrt{1+\dfrac{4\Omega^{2}}{\omega_{X}^{2}}}\right)-\frac{i\gamma}{2}\left(1+\frac{1}{\sqrt{1+\dfrac{4\Omega^{2}}{\omega_{X}^{2}}}}\right).$
(11)

The renormalized light-matter velocity $\tilde{v}$ and the photon Hopfield
coefficient of forward and backward propagating lower polaritons around $k=0$
are given by,

$\displaystyle\tilde{v}$
$\displaystyle=\frac{v}{2}\left(1+\frac{\omega_{X}}{\sqrt{\omega_{X}^{2}+4\Omega^{2}}}\right),$
(12) $\displaystyle|C_{0}|^{2}$
$\displaystyle=\frac{4\Omega^{2}}{4\Omega^{2}+\left|\omega_{X}+i\gamma-\sqrt{(\omega_{X}+i\gamma)^{2}+4\Omega^{2}}\right|^{2}}.$
(13)

The form of Eq. (10) implies that, in the truncated basis of lower forward
$\lvert L,e^{iKx}\rangle$ and backward $\lvert L,e^{-iKx}\rangle$ polaritons,
we can describe the system with the following massive non-Hermitian Dirac
operator,

$H_{D}=\tilde{\omega}_{L}\mathds{1}_{2\times
2}+\begin{pmatrix}\tilde{v}k&(\kappa-i\gamma e^{i\varphi})|C_{0}|^{2}\\\
(\kappa-i\gamma e^{-i\varphi})|C_{0}|^{2}&-\tilde{v}k\end{pmatrix},$ (14)

with new symmetric and antisymmetric lower polariton eigenstates
$|L,\pm\rangle=A_{\pm}\lvert L,e^{iKx}\rangle+B_{\pm}\lvert L,e^{-iKx}\rangle$
with eigenenergies $\omega_{L,\pm}$. One limitation of Eq. (10) is that it
neglects the dependence of the photonic Hopfield coefficient (13) on both
momentum and the original diffractive coupling between the counterpropagating
photons. However, if $\kappa\ll\Omega\ll\omega_{X}^{(0)}$ then Eq. (10)
remains accurate and implies that waveguided TE polaritons behave
approximately as Dirac particles with renormalized velocity $\tilde{v}$ and
gap opening. In Fig. 2(c) and 2(d) we compare our approximated dispersion (10)
(dashed lines) with the exact lower polariton dispersion relation (9) (solid
lines) for values extracted from RCWA simulations, $\Omega/\kappa=1.8$,
$\hbar\omega_{X}^{(0)}/\kappa=5.5$, and $\gamma_{\text{nr}}=\gamma$, and
observe very good agreement. We note that the parameters used in our study
accurately represent a real example of a photonic grating analyzed using RCWA
in Appendix C. The above underpins the feasibility in creating photonic
samples that permit study of Dirac polariton quasi-BIC physics.

##  III Mean-field formalism

Figure 3: Normalized total condensate density in (a-d) real space and (e-h)
energy resolved momentum space for increasing pump power
$P_{0}=\\{1.5,2,2.4,3\\}P_{\text{th}}$. The overlaid curves in panel (c) are
taken at different time steps to show the nonstationary oscillations. Same
parameters as in Fig. 2 are used here with a pump spot size of
$\text{FWHM}=20$ $\mu$m (full width at half maximum). Note that we have
shifted the zero energy point into the center of the bandgap. The labels
$|0\rangle$ and $|1\rangle$ denote the condensate fraction occupying the trap
ground state and first excited state, respectively. White solid curves denote
the single-particle guided polariton dispersions in the photonic crystal from
Eq. (10).

To create a macroscopic coherent quantum state of the polaritons by means of
Bose–Einstein condensation, the system is excited by an external nonresonant
laser [29, 58]. This creates hot free charge carriers which relax in energy to
form a reservoir of excitons at the so-called “bottleneck region” in the lower
polariton dispersion relation [59] denoted by the density parameter $n_{R}$.
When then density of the pumped reservoir is sufficiently high, the polariton
occupation number accumulated at a particular level (such as the quasi-BIC)
can exceed unity and stimulated scattering of polaritons into this state
starts. This signifies the spontaneous breaking of symmetry and non-
equilibrium Bose-Einstein condensation into a single quantum state [1], marked
by a threshold power $P_{\text{th}}$. Following the mean-field theory [59, 60]
and our approximate single particle Dirac operator (14), the condensate can be
described by a two-component macroscopic spinor wave function, or an order
parameter, $\Psi(x,t)=(\psi_{+1},\psi_{-1})^{\text{T}}$ where $\pm 1$ denotes
the lower forward- and backward propagating polaritons, respectively. Note
that subscripts $(\pm)$ without "1" denote the symmetric-antisymmetric basis
which should not be confused here.

The generalized Gross-Pitaevskii equation for the condensate in real space can
be found to be

$\displaystyle i\frac{\partial\Psi}{\partial t}$
$\displaystyle=\bigg{[}H_{D}+g\Psi^{\dagger}\Psi$
$\displaystyle+g_{R}\left(n_{R}+\frac{\eta
P(x)}{\Gamma_{R}}\right)+i\frac{Rn_{R}}{2}\bigg{]}\Psi,$ (15)
$\displaystyle\frac{\partial n_{R}}{\partial t}$
$\displaystyle=P(x)-(\Gamma_{R}+R\Psi^{\dagger}\Psi)n_{R}.$ (16)

The first term corresponds to the single particle dynamics (14) in real space
obtained by substituting $k\to-i\partial_{x}$. In the second term the
polaritons are assumed to interact via short-range interaction described
phenomenologically by the repulsive non-linear term $g\Psi^{\dagger}\Psi$.
Moreover, polaritons also interact repulsively with strength $g_{R}>g$ against
any background excitons whose density can be divided in to the bottleneck part
$n_{R}$ and a static inactive dark exciton background parametrized by the
dimensionless number $\eta$. The last term describes the stimulated scattering
of reservoir excitons into the condensate at a rate $R$. The term
$P(x)=P_{0}e^{-x^{2}/2w^{2}}$ describes the continuous wave nonresonant
Gaussian pump of waist $w$, and $\Gamma_{R}$ is the average reservoir exciton
redistribution and decay rate.

For large waist $w$, the pumping can be considered to be uniform,
$P(x)=P_{0}$. In this case, the condensation threshold is given by
$P_{\text{th}}=-2\Gamma_{R}\text{max}\\{\text{Im}{[\omega_{L,\pm}]}\\}/R$
where $\text{Im}{[\omega_{L,\pm}]}<0$ and the maximum is taken over $\pm$. In
this uniform case the maximum (i.e., minimum losses) will always correspond to
the quasi-BIC. For a finite size pump the threshold actually increases due to
finite gain region effects and must be determined numerically. For details on
numerical modeling of the mean field equations please see Appendix D.

Below threshold one has $\Psi^{\dagger}\Psi=0$ in the long time limit and the
reservoir converges to the steady state $n_{R}=P(x)/\Gamma_{R}$. We can then
define a pump induced potential term acting on the Dirac polaritons,

$V(x)=g_{R}(1+\eta)\frac{P(x)}{\Gamma_{R}}>0.$ (17)

The above expression gives a good estimate for the optical trap felt by the
condensate when pumped only weakly above threshold. For the positive branch
polaritons, $V(x)$ acts as a repulsive gain region which, if tightly focused
into a small enough spot, results in so-called ballistic condensation [3, 5].
For the negative branch polaritons however, it acts as an attractive gain
region, pulling in generated polaritons and trapping them efficiently with a
much lower threshold [29, 35, 60].

###  III.1 Mode dropping

Figure 3 shows the total condensate density $\rho=\Psi^{\dagger}\Psi$ in
(upper row) real space and (lower row) energy-resolved momentum space for
increasing pump power density $P_{0}$. At low powers above threshold [Fig.
3(a,e)], the condensate first occupies the ground state in the effective
optical trap since it is closest to the quasi-BIC and has the lowest particle
losses [60].

Increasing the power we observe monotonic blueshift of the condensate level
[compare Fig. 3(a) and 3(b)]. By increasing the power, more trap states become
available for the condensate to populate and we can locate stable cyclical
solutions (i.e., limit cycles) in which the condensate becomes nonstationary
and coherently divided between two neighbouring trap modes [61] in the same
branch. Such cyclical solutions [see Fig. 3(c)] usually appear through Hopf
bifurcations when one fixed point attractor deteriorates and another takes
over as parameters of the system are tuned. When the power is further
increased, the blueshift is so strong that the fundamental trap mode is swept
into the upper positive-mass band with increased losses. Consequently, the
condensate abandons the fundamental mode and shifts its population into the
neighbouring higher-order mode at at lower energies [see Fig. 3(d)]. In this
sense, the condensate "drops" from one trap mode to the next, as predicted by
Nigro et al. [60]. Increasing the power further, we observe periodically the
same mode-dropping behaviour as subsequent higher order trap modes form in
vicinity of the quasi-BIC and blueshift up into the lossy positive-mass band.

This power driven change in the condensate structure is in agreement with
recent experimental observations [29]. Energetically, this behaviour is in
sharp contrast to optically trapped polariton condensates in planar cavities
[62, 61, 5] where stronger pumping results in a condensate dropping into lower
order trap modes until it reaches the ground state.

###  III.2 Negative-positive mass superposition

Next, we characterize the interplay between the quasi-BIC state and the
negative-mass trapping mechanism coming from the localized pumping area [see
Eq. (17)]. As mentioned around Eq. (2) the photonic grating introduces a
complex coupling parameter between the counterpropagating photons (which
carries into the polariton modes). Up until now, we have taken $\varphi=0$
which leads to a quasi-BIC in the lower (symmetric) branch in Fig. 2(c) and
2(d). If $\varphi=\pi$ then the quasi-BIC would instead form in the upper
branch [32]. We perform a scan across both $\varphi$ and the pump power
$P_{0}$ and investigate the difference between symmetric and antisymmetric
condensate occupation

$\Delta_{\rho}=\rho_{+}-\rho_{-}.$ (18)

where

$\rho_{\pm}=|\langle\Psi_{\pm}|\Psi\rangle|^{2}.$ (19)

and $|\Psi_{\pm}\rangle$ are the symmetric and antisymmetric eigenstates of
the potential-free Dirac operator (14). In this sense, the quantity
$\Delta_{\rho}$ is similar to the projection of a two level quantum system
onto the axis connecting the north and south antipodal points of the Bloch
sphere. If $\Delta_{\rho}<0$ then most of the condensate forms in the lower
branch, whereas if $\Delta_{\rho}>0$ then the condensate forms in the upper
branch.

Figure 4: Average condensate population difference
$\langle\Delta_{\rho}\rangle$ between the symmetric and antisymmetric
projections (branches) as a function of pump power, and the dissipative
coupling parameter $\varphi$, and two different sizes of the pump spot. The
colorscale is saturated around $\pm 1$ to more clearly show the white-dark
regions. The red contours show total ($\rho=\Psi^{\dagger}\Psi$) isodensity
curves. Other parameters are the same as in Fig. 3.

The results are shown in Fig. 4 for two different sizes of pump spots (FWHM
$=50$ and $200$ $\mu$m) where each pixel corresponds to the spatiotemporal
average of $\langle\Delta_{\rho}\rangle$ over the entire simulation grid and
integration time. The results show that for $\varphi\sim 0$, when the quasi-
BIC is in the lower antisymmetric branch, we always have condensation in the
same branch (dark region) as seen in experiments [29]. Interestingly, for
smaller pump spots in Fig. 4(a) condensation still takes place in the lower
antisymmetric branch even when the quasi-BIC has moved to the upper symmetric
branch as can be seen from the weakly dark region at high powers around
$\varphi\sim\pi$. This implies that the optical trapping mechanism can be more
efficient in reducing transverse ($x$-direction) losses than the quasi-BIC in
reducing out-of-plane ($z$-direction) losses. Nevertheless, the presence of
the BIC has a dramatic effect on the condensation threshold curve
approximately given by far-left red contour. Indeed, when the BIC is in the
negative mass branch the optical trapping and protection from the continuum
complement each other to lower the power needed to reach bosonic stimulation.

To explore the reduction of the optical trapping effect on lower branch
polaritons we repeat the calculation using a much larger pump spot of FWHM
$=200$ $\mu$m in Fig. 4(b). The much wider pump spot imposes a weaker
confinement compared to the narrow spot. Indeed, now around $\varphi\sim\pi$
we see that condensation starts taking place in the upper branch (white
region), following the quasi-BIC. This means that condensation can be
optically adjusted between the lower and the upper branch by simply changing
the size of the pump spot in a given sample.

Figure 5: A stable superposition of trapped and ballistic polaritons with a
single pump spot. (a,b) Density and (c,d) phase of the symmetric and
antisymmetric polaritons. The transparency of the phase map is proportional to
the condensate density. Red trembling trajectory shows the center-of-mass of
the condensate. (e,f) Corresponding condensate densities in Fourier space
showing that each component belongs to different branches. Other parameters
are the same as in Fig. 3.

At the interface of such qualitatively different condensate solutions [i.e.,
dark and bright regions in Fig. 4(b)] more exotic patterns might appear. We
propose that by tuning the size of the trap, the grating pitch ($\varphi$),
and pump power one can achieve simultaneous condensation in the upper and
lower branches. This corresponds to a stable coherent mixture of positive and
negative mass condensate polaritons (i.e., ballistic and trapped polaritons).
Interestingly, such a solution bears similarities to the famous Zitterbewegung
effect, the trembling motion of relativistic particles, but here in a driven-
dissipative setting [47]. Interference between the upper and lower branches
causes the center-of-mass of the condensate, $\langle
x\rangle=\int\Psi^{\dagger}\hat{x}\Psi dx/\int\Psi^{\dagger}\Psi dx$, to
jitter in time (see red oscillating curve). In Fig. 5(a,b) and (c,d) we show
the density and phase of symmetric and antisymmetric condensate polaritons
belonging to such a solution, obtained at the location of the red diamond
marker in Fig. 4(b). This solution is also a limit cycle, a stable coherent
superposition of ballistic upper branch and trapped lower branch polaritons,
with THz Rabi oscillations, $\propto\cos^{2}{(\pi\Delta Et/h)}$ (not to be
confused with the light-matter Rabi frequency $\Omega$). Here $h$ is the
Planck’s constant and $\Delta E$ is the energy splitting between the ballistic
and the trapped condensate levels [see Fig. 5(e,f)]. The ballistic nature of
the upper branch fluid manifests in its more delocalized nature and more
rapidly varying phase profile along the $x$-direction.

##  IV Nearfield and Farfield pattern

In the previous sections, we demonstrated that the dynamics of polariton
condensation is determined by the wavefunction and population of polaritons,
derived from the generalized Gross-Pitaevskii equation. To probe polaritons in
practical realizations, most experimental works rely on detecting their
photonic component using far-field setups in either real or momentum space
[29, 35]. We remind the reader that the polariton state is explicitly related
to the emitted photon state vector $\Psi\propto\Phi$ through the photonic
Hopfield coefficient [1]. For polaritons in microcavities, the distinctions
between nearfield and farfield, both directly determined by the polariton
population, were often overlooked. Furthermore, recent studies show that
nearfield setups can probe the local wavefunction of polaritons when they are
not embedded in thick vertical microcavities [63]. As a result, it is crucial
to bridge the polariton wavefunction with the pattern of the electric field in
both nearfield and farfield scenarios. Interestingly, for a BIC, the
relationship between nearfield and farfield turns out to be radically
different.

We note that there are certain inconsistent use of the terms ‘nearfield’ and
‘farfield’ in the literature of polaritonics. In this work, we adopt the
conventional definition presented in Ref. [64]: ‘nearfield’ refers to the
light field that is confined within structures and can only be probed using
evanescent techniques such as scanning near-field optical microscopy (SNOM),
while ‘farfield’ pertains to the light field that propagates through space and
can be probed using conventional imaging techniques. This is different from
corresponding notions that have been used in several polaritonic experiments
[65, 66, 67, 68], where ‘nearfield’ actually refers to the usual farfield
measurements in real space, and ‘farfield’ refers to the usual farfield
measurements in momentum space.

In fact, the nearfield and farfield pattern of BIC can be deduced rather
straightforwardly from the effective Dirac Hamiltonian for the photonic
component (2). Indeed, we start with rewriting the Hamiltonian (2) in real-
space by substituting $k\to-i\partial_{x}$ and denoting the two-component
photon state vector as $\Phi=(\phi_{+1},\phi_{-1})^{\text{T}}$ where
$\phi_{\pm 1}$ are the coefficients of the forward and backward propagating
photons $\lvert e^{i(k\pm K)x}\rangle$,

$H_{\text{ph}}=\begin{pmatrix}-vi\partial_{x}&\kappa\\\
\kappa&+vi\partial_{x}\end{pmatrix}-i\gamma\begin{pmatrix}1&e^{i\varphi}\\\
e^{-i\varphi}&1\end{pmatrix}.$ (20)

Then using the dynamical equation
$\partial_{t}\Phi(x,t)=-iH_{\text{ph}}\Phi(x,t)$ and its conjugation, we
arrive directly at the continuity equation for the photon field,

$\partial_{t}(\Phi^{\dagger}\Phi)=-\partial_{x}(\Phi^{\dagger}\sigma_{z}\Phi)-\gamma\Phi^{\dagger}\begin{pmatrix}1&e^{i\varphi}\\\
e^{-i\varphi}&1\end{pmatrix}\Phi.$ (21)

where $\sigma_{z}$ is the third Pauli matrix. From the continuity equation,
one infers, as usual, that $\Phi^{\dagger}\Phi$ describe the nearfield
intensity, and $\Phi^{\dagger}\sigma_{z}\Phi$ describes the photon density
current. The last term,

$\Phi^{\dagger}\begin{pmatrix}1&e^{i\varphi}\\\
e^{-i\varphi}&1\end{pmatrix}\Phi=\gamma|\phi_{+1}+e^{i\varphi}\phi_{-1}|^{2},$
(22)

is then identified as losses by means of radiation into the farfield. In this
way, the expression for nearfield and farfield intensities have been obtained
only by formally investigating the structure of the effective Dirac equation.
Their difference can be appreciated from the additional interference term
between the forward and backward propagating polaritons when (22) is expanded.
We show in Appendix B how they can also be understood from the microscopic
theory. Figure 6 shows the difference between the nearfield (red curves) and
farfield (blue curves) intensities of the emitted light from the condensate in
both real space and momentum space corresponding to our results in Fig. 3.
Notably, in [29] only the farfield was measured showing emission profiles
which agree very well with our results.

Figure 6: Comparison of the nearfield (35) and farfield (38) intensity
(individually normalized) of the emitted cavity light coming from the
luminescent polariton condensate. (a,b) show the trap ground state and first
excited state condensate emission corresponding to Figs. 3(a) and 3(d). (c,d)
Show the corresponding momentum space emission belonging to Figs. 3(e) and
3(h).

##  V Conclusions

We have introduced the concept of Dirac polaritons in photonic crystal
gratings—one dimensional photonic crystal slabs—containing excitonic
resonances and symmetry protected photonic modes or bound states in the
continuum. We developed the single particle theory of these effectively
relativistic bosonic elementary excitations of light and matter, followed by
intuitive extension to the many-body picture through the mean-field formalism.
We propose a generalized Gross-Pitaevskii model to describe BIC-facilitated
condensation of polaritons into pump-induced optical traps. Our findings are
in excellent agreement with recent experimental observations on multi-quantum-
well structures [29] and are applicable to other forms of optically active
materials such as transition metal dichalcogenide monolayers including MoSe2
[37, 38] or WS2 [39, 40, 41]; or hybrid organic-inorganic perovskites [42, 43,
44, 45].

Our theory is fully generalized towards photonic gratings with broken lateral
symmetry which manifests in tunable diffractive coupling mechanism between
guided photons and the continuum, allowing us to continuously tune the BIC
from the lower energy Dirac branch to the upper. This gives powerful control
over the polariton condensation threshold, and final state stimulation. We
have also clarified on the distinction between farfield and nearfield emission
patterns which becomes more important when dealing with subwavelength-grated
photonic structures, and therefore will be important for future works on
polaritonic crystals.

In particular, in mean-field simulations, we have identified peculiar
zitterbewegung-like solutions in the driven Dirac polariton condensate which
manifests in spontaneous formation of coherent superposition of upper-branch
(positive) and lower-branch (negative) polaritons. The implications of such a
hybrid quantum fluid are the two very different coupling mechanisms when
neighbouring condensates are added. One one hand, the high energy component of
the condensate will interact ballistically with its neighbours, on the other,
the low energy component will interact evanescently with its neighbours. Such
system could offer novel patterns of synchronicity between multi-component
nonlinear oscillators with contrasting coupling mechanisms with competing
coherence scales and time-scales of domain wall formation.

##  VI Acknowledgements

H.S. acknowledges the project No. 2022/45/P/ST3/00467 co-funded by the Polish
National Science Centre and the European Union Framework Programme for
Research and Innovation Horizon 2020 under the Marie Skłodowska-Curie grant
agreement No. 945339. H.C.N. acknowledges the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation, project numbers 447948357 and 440958198),
the Sino-German Center for Research Promotion (Project M-0294), the German
Ministry of Education and Research (Project QuKuK, BMBF Grant No. 16KIS1618K)
and the ERC (Consolidator Grant 683107/TempoQ).

##  Appendix A Derivation of the non-hermitian Hamiltonian

Here we derive the effective photonic Dirac Hamiltonian (2) from the
microscopic consideration. To develop a perturbative theory for the guided
photons, we first remark that in the absence of the modulated refraction index
in the $x$-direction, the eigenmodes of the lowest band of the system with
frequency $\omega_{q}$ are the plane waves $\lvert e^{iqx}\rangle$. For now,
we ignore the confinement of the wave function in $z$-direction, which is
assumed to only weakly dependent on $q$. We also ignore the free evolution in
the $y$-direction and assume that the system is time-reversal symmetric
$\omega_{q}=\omega_{-q}$ without any loss of generality. Note that modes below
the light cone are lossless, while the one above the light cone are lossy; see
Fig. 1(a).

The periodic modulation—with period $a$—of the refraction index in the $x$
direction introduces an periodic potential $u(x)=u(x+a)$ acting on the
photons. As an expression of the Bragg reflection, the potential then couples
modes $\lvert e^{iqx}\rangle$ with that of $\lvert e^{i(q+nK)x}\rangle$ where
$K=2\pi/a$ is the primitive reciprocal lattice vector and $n$ is an integer.
The relevant matrix elements are written $U_{n}=\langle e^{i(q+nK)x}\rvert
u(x)\lvert e^{iqx}\rangle$, which are the $n$-th Fourier coefficients of the
potential. In particular, for any integer $n$, the potential couples
degenerate modes of the same frequency, $q=nK/2$ and $q=-nK/2$.

We are interested in the system excited at frequencies corresponding to $q$
around $\pm K$. The relevant wavevectors are therefore of $\pm(K+k)$ with
$k\ll K$. The two modes $\lvert e^{i(k+K)x}\rangle$ and $\lvert
e^{i(k-K)x}\rangle$ are coupled by the potential $u(x)$ with matrix elements
$U_{2}=\langle e^{i(k+K)x}\rvert u(x)\lvert e^{+i(k-K)x}\rangle=\langle
1\rvert u(x)\lvert e^{-i2Kx}\rangle$ and $U_{2}^{\ast}=\langle 1\rvert
u(x)\lvert e^{+i2Kx}\rangle$, which are simply the second Fourier coefficients
of $u(x)$.

Being guided modes located below the light-line, both modes $\lvert e^{i(k\pm
K)x}\rangle$ are technically lossless. They become, however, lossy through
coupling to lossy modes at low momentum, $\lvert e^{ikx}\rangle$. Notice that
these lossy modes at $\lvert e^{ikx}\rangle$ are distributed on the Fabry-
Pérot modes of the stack; see Fig. 1(b). Therefore, in order to describe
matrix elements of these scattering processes, we need to include the
confinement wave function in the $z$-direction. Including these confinement
factors, the full wavefunctions of the two modes $\lvert e^{i(k\pm
K)x}\rangle$ are $\lvert\chi_{K}(z)e^{i(k\pm K)x}\rangle$. We ignore the
dependence of the confinement wavefunction $\chi$ on $k$. Further, we have
used $\chi_{+K}(z)=\chi_{-K}(z)$ because the non-perturbed structure is
symmetric under $x$-reflection. The lossy modes at $\lvert e^{ikx}\rangle$ are
modeled by $\lvert\chi^{(0)}(z)e^{ikx}\rangle$,
$\lvert\chi^{(1)}(z)e^{ikx}\rangle$, …. The matrix element scattering
$\lvert\chi_{K}(z)e^{i(k\pm K)x}\rangle$ into the $p$-th lossy mode
$\lvert\chi^{(p)}(z)e^{ikx}\rangle$ is given by
$\langle\chi^{(p)}(z)|\chi_{K}(z)\rangle\langle 1\rvert u(x)\lvert e^{i(k\pm
K)x}\rangle=c^{p}U_{\pm 1}$. Notice that the first Fourier coefficients of
$u(x)$ are related by $U_{-1}=U_{1}^{\ast}$. For simplicity, in the following
we consider the coupling to only the $0$-th lossy mode; the analysis can be
extended to coupling to many lossy Fabry-Pérot modes in a straightforward way.

Restricted to the space spanned by three modes
$\lvert+1\rangle=\lvert\chi_{K}(z)e^{i(k+K)x}\rangle$,
$\lvert-1\rangle=\lvert\chi_{K}({z})e^{i(k-K)x}\rangle$ and $\lvert
0\rangle=\lvert\chi^{0}(z)e^{ikx}\rangle$, a general wavefunction can be
written as

$|\tilde{\Phi}\rangle=\tilde{\phi}^{k}_{+1}\lvert+1\rangle+\tilde{\phi}^{k}_{-1}\lvert-1\rangle+\tilde{\phi}^{k}_{0}\lvert
0\rangle.$ (23)

The evolution of the coefficients
$(\tilde{\phi}^{k}_{+1},\tilde{\phi}^{k}_{-1},\tilde{\phi}^{k}_{0})$ follows
the Schrödinger-like equation

$\frac{d}{dt}\begin{pmatrix}\tilde{\phi}^{k}_{+1}\\\ \tilde{\phi}^{k}_{-1}\\\
\tilde{\phi}^{k}_{0}\end{pmatrix}=-i\begin{pmatrix}\omega_{K}+vk&U_{2}&c_{0}U_{1}\\\
U_{2}^{\ast}&\omega_{K}-vk&c_{0}U_{1}^{\ast}\\\
c_{0}^{\ast}U_{1}^{\ast}&c_{0}^{\ast}U_{1}&\omega_{0}-i\gamma_{0}\end{pmatrix}\begin{pmatrix}\tilde{\phi}^{k}_{+1}\\\
\tilde{\phi}^{k}_{-1}\\\ \tilde{\phi}^{k}_{0}\end{pmatrix},$ (24)

where $\gamma_{0}$ is the decay rate of the low-momentum mode
$\lvert\chi^{(0)}(z)e^{ikx}\rangle$, which is assumed to vary negligibly for
small $k$. We have also linearized the dispersal relation so that
$\omega_{K+k}\approx\omega_{K}+vk$ and $\omega_{-K+k}\approx\omega_{K}-vk$,
with $v$ being the light velocity near $K$. By adjusting a global phase, one
can also assume that $\omega_{K}=0$ and $\omega_{0}$ can be then replaced by
the difference in frequency $\Delta=\omega_{0}-\omega_{K}$.

Assuming that the decay rate of the lossy mode $\gamma_{0}$ is much faster
than the matrix element $U_{2}$, one can adiabatically eliminate
$\tilde{\phi}^{k}_{0}$. This is done by solving $\tilde{\phi}^{k}_{0}$ in
terms of $\tilde{\phi}^{k}_{+1}$ and $\tilde{\phi}^{k}_{-1}$ as

$\tilde{\phi}^{k}_{0}=-ic_{0}^{\ast}\int_{0}^{t}d\tau
e^{-(\gamma_{0}+i\Delta)\tau}\begin{pmatrix}U_{1}^{*}&U_{1}\end{pmatrix}\begin{pmatrix}\tilde{\phi}_{+1}^{k}(t-\tau)\\\
\tilde{\phi}_{-1}^{k}(t-\tau)\end{pmatrix}.$ (25)

As the decaying of the lossy mode $\gamma_{0}$ is fast in comparison to the
dynamics of the confined mode, one can make the Markovian approximation
$\tilde{\phi}_{+1}^{k}(t-\tau)\approx\tilde{\phi}_{+1}^{k}(t)$,
$\tilde{\phi}_{-1}^{k}(t)=\tilde{\phi}_{+1}^{k}(t-\tau)$ and
$\int_{0}^{t}d\tau e^{-(\gamma_{0}+i\Delta)\tau}\approx\int_{0}^{\infty}d\tau
e^{-(\gamma_{0}+i\Delta)\tau}=1/(\gamma_{0}+i\Delta)\approx 1/\gamma_{0}$. We
have we also assumed $\gamma_{0}\gg\Delta$ in the last approximation. In the
end, we then obtain

$\tilde{\phi}^{k}_{0}=\frac{-ic_{0}^{\ast}}{\gamma_{0}}\begin{pmatrix}U_{1}^{*}&U_{1}\end{pmatrix}\begin{pmatrix}\tilde{\phi}_{+1}^{k}\\\
\tilde{\phi}_{-1}^{k}\end{pmatrix}.$ (26)

Inserting (26) into (24), we obtain the evolution equation for
$\tilde{\phi}^{k}_{+1}$ and $\tilde{\phi}^{k}_{-1}$ in the form

$\frac{d}{dt}\begin{pmatrix}\tilde{\phi}_{+1}^{k}\\\ \tilde{\phi}_{-1}^{k}\\\
\end{pmatrix}=-iH\begin{pmatrix}\tilde{\phi}_{+1}^{k}\\\
\tilde{\phi}_{-1}^{k}\\\ \end{pmatrix}$ (27)

with the non-hermitian Hamiltonian in momentum space

$H=\begin{pmatrix}+ivk&U_{2}\\\ U_{2}^{\ast}&-ivk\\\
\end{pmatrix}-i\gamma\begin{pmatrix}1&e^{+i2\varphi_{1}}\\\
e^{-i2\varphi_{1}}&1\end{pmatrix}$ (28)

where $\varphi_{1}$ is defined by $U_{1}=|U_{1}|e^{i\varphi_{1}}$ and
$\gamma=|c_{0}|^{2}\frac{|U_{1}|^{2}}{\gamma_{0}}$. It is interesting to
notice that the indirect loss rate $\gamma$ for the two modes
$\lvert\chi_{K}(z)e^{\pm iKx}\rangle$ is inversely proportional to the lossy
rate $\gamma_{0}$ of the mode $\lvert\chi^{(0)}(z)e^{ikx}\rangle$, indicating
an analogy of the Zeno effect in quantum system [69]. Indeed, a strong
"measurement regime" that corresponds to a very leaky channel,
$\gamma_{0}\gg|U_{1}|$, will "freeze" the population of
$\lvert\chi_{K}(z)e^{\pm iKx}\rangle$ because $\gamma\approx 0$. A similar
setup, combining lossless waveguides and a lossy one, has been recently
proposed to demonstrate the optical Zeno effect [70].

The Fourier coeffiecient $U_{2}$ in (28) is generally a complex number. We
denote $U_{2}=|U_{2}|e^{i\varphi_{2}}$ and eliminate the phase $\varphi_{2}$
by the following unitary transformation,

$\begin{pmatrix}\tilde{\phi}_{+}^{k}\\\
\tilde{\phi}_{-}^{k}\end{pmatrix}=\begin{pmatrix}e^{+i\varphi_{2}/2}&0\\\
0&e^{-i\varphi_{2}/2}\end{pmatrix}\begin{pmatrix}{\phi}_{+}^{k}\\\
{\phi}_{-}^{k}\end{pmatrix}.$ (29)

One then obtains the non-hermitian Hamiltonian for $\phi_{+1}^{k}$ and
$\phi_{-1}^{k}$ as

$H=\begin{pmatrix}+ivk&|U_{2}|\\\ |U_{2}|&-ivk\\\
\end{pmatrix}-i\gamma\begin{pmatrix}1&e^{+i\varphi}\\\
e^{-i\varphi}&1\end{pmatrix}$ (30)

with $\varphi=2\varphi_{1}-\varphi_{2}$. This is the effective Hamiltonian (2)
introduced in the main text with $\kappa=|U_{2}|$.

##  Appendix B Derivation of the farfield and nearfield intensity from the
microscopic description

We start with remarking again that the loss due to radiation into the farfield
of the spinor polaritons $\Psi$ inherits directly from the loss of the
photonic component $\Phi$. Therefore the farfield pattern of the polariton can
be understood directly from its photonic components.

Generally, the effective wave function $\Phi(x)$ is a superposition of
different wavevectors $k$,

$\Phi(x)=\sum_{k}e^{ikx}\begin{pmatrix}\phi_{+1}^{k}\\\
\phi_{-1}^{k}\end{pmatrix}.$ (31)

At the microscopic level, this is a plane wave of

$\lvert\Phi\rangle=\sum_{k}\phi^{k}_{+1}e^{i\varphi_{2}/2}\lvert
e^{i(k+K)x}\rangle+\phi^{k}_{-1}e^{-i\varphi_{2}/2}\lvert e^{i(k-K)x}\rangle,$
(32)

if we ignore the confining mode function in the $z$-direction and the lossy
mode as comparison to (23). The latter is only relevant to the lossy dynamics.
Explicitly in terms of the electric field, one has

$\vec{E}^{\text{near}}\propto\left(\phi_{+1}^{k}e^{i\varphi_{2}/2}e^{iKx}+\phi_{-1}^{k}e^{-i\varphi_{2}/2}e^{-iKx}\right)e^{ikx}\vec{u}_{y},$
(33)

where $\vec{u}_{y}$ denotes the polarization direction. Or using the real
space representation $\Phi(x)$, we can write

$\vec{{E}}^{\text{near}}\propto\left[\phi_{+1}(x)e^{i\varphi_{2}/2}e^{iKx}+\phi_{-1}(x)e^{-i\varphi_{2}/2}e^{-iKx}\right]\vec{u}_{y}.$
(34)

Recall that we are working in the regime where $\phi^{k}_{+1}$ and
$\phi^{k}_{-1}$ are only significant at wavelengths much smaller than
$K=\pi/a$. Equivalently, $\phi_{+1}(x)$ and $\phi_{-1}(x)$, which are referred
to as _envelope functions_ , vary much slower than $e^{\pm iKx}$, which are
referred to as _core functions_. Averaging out the fast fluctuation at the
wavevector $K$, the nearfield intensity can then be obtained as:

$I^{\text{near}}(x)=|\phi_{+1}(x)|^{2}+|\phi_{-1}(x)|^{2}.$ (35)

This agrees with the formal derivation from the effective photonic Dirac
equation (21).

As for the farfield, we notice that the farfield is the observation of the
lossy mode $\phi^{k}_{0}$, which is given by (26) for a single $k$. With a
superposition of different wavevectors (31), one has

$\phi_{0}(x)\propto\sum_{k}[U_{1}^{\ast}e^{+i\varphi_{2}/2}\phi^{k}_{+1}+U_{1}e^{-i\varphi_{2}/2}\phi^{k}_{-1}].$
(36)

Explicitly in terms of the electric field, this corresponds to

$\vec{E}^{\text{far}}(x)\propto\sum_{k}e^{ikx}(\phi^{k}_{+1}+e^{i\varphi}\phi^{k}_{-1})\vec{u}_{y},$
(37)

where we have again used $\varphi=2\varphi_{1}-\varphi_{2}$. The farfield
intensity is then obtained as

$I^{\text{far}}(x)\propto|\phi_{+1}(x)+e^{i\varphi}\phi_{-1}(x)|^{2},$ (38)

which coincides with (22) under the approximation that $\Psi$ is replaced by
its photonic components.

##  Appendix C Design and numerical simulations for pratical realizations

We propose realistic a stack of dielectric layers and quantum wells, based on
the experimental works from [29, 35]. The sample stack is composed of a
waveguide core made of 12 GaAs quantum wells (QWs) 20 nm thick and 13
Al0.4Ga0.6As barriers 20 nm thick grown on an Al0.8Ga0.2As 500 nm thick
cladding layer. The cladding and the GaAs substrated are seperated by a 50 nm-
AlAs layer. The whole stack is capped by a 10 nm GaAs layer.

For numerical simulations of photonic modes in the gratings, the excitonic
resonances in the GaAs QWs are removed in the dielectric function. The
photonic modes are calculated by numerical simulations based on Rigorous
Coupled-Wave Analysis (RCWA) method with the S4 package provided by the Fan
Group at the Stanford Electrical Engineering Department [71]. The refractive
index of the QWs, barriers and the cladding are:
$n_{\text{QWs}}=3.547+0.0001i$, $n_{\text{barrier}}=3.3$,
$n_{\text{cladding}}=3.063$. The imaginary part in the refractive index of the
QWs are simply added to probe the photonic modes in absorption simulations. We
only calculate TE (transverse electric) photonic modes since the TM
(transverse magnetic) photonic modes are inefficiently coupled to in-plane
excitonic dipoles of the QWs.

We also employ RCWA method for the numerical simulations of polariton modes.
To do so, a Lorentz oscillator at 1527.4 meV with 0.001 eV2 oscillator
strength and 0.35 meV linewidth is added in the dielectric function of the
QWs. This excitonic resonance corresponds to the heavy-hole excitons of the
QWs.

Figure 7: RCWA simulations versus effective Dirac-photon theory. (a) Sketch of
the grating design. (b,e) Photonic and polaritonic dispersions, numerically
calculated by RCWA, for $a=243$ nm, $FF=0.34$, and 110 nm of etching. The
solid blue lines correspond to theoretical dispersion from the effective
theory. (c,d) Real and imaginary part for the energy of the photonic modes.
Symbols are numerical results extracted from RCWA simulation in (b). Solid
black line is the theoretical prediction using (2). (f,g) Real and imaginary
part for the energy of the two lower polaritonic modes. Symbols are numerical
results extracted from RCWA simulation in (e). Solid black line is the
theoretical prediction using (14). The parameters for effective theory are
$\hbar\gamma=0.19$ meV, $\hbar\kappa=1.75$ meV, $\hbar v=56.61$ meV $\mu$m,
$\hbar\Omega=3.2$ meV, $\hbar\omega^{(0)_{X}}=2.41$ meV,
$\hbar\gamma_{nr}=0.18$ meV and $\varphi=0$. Figure 8: Varying the exciton-
photon energy detuning. (a) Energy at $k=0$ of the photonic modes, given by
RCWA simulation results as function of the period $a$. (b) The dependence of
the exciton-photon detuning as function of $a$. The RCWA are performed for
grating of filling fraction $FF=0.37$ and the etching depth of 110 nm. Figure
9: Continuously tuning $\kappa$ and implementing $\pi$-jump to $\varphi$ via
the filling fraction $FF$. (a) Sketch of the grating design. (b) Photonic
dispersion, numerically calculated by RCWA, for $FF=0.52$, $FF=0.62$ and
$FF=0.8$. (c,d) Real and imaginary part for the energy of the photonic modes
when scanning $FF$ from 0.4 to 0.86. (e,f) $\kappa$ and $FF$ when scanning
$FF$ from 0.4 to 0.86. Figure 10: Continuously tuning $\varphi$ via the
symmetry-breaking coefficient $\alpha$. (a) Sketch of the grating design with
double-period symmetry-breaking. b) Photonic dispersion, numerically
calculated by RCWA, for $\alpha=0$ and $\alpha=0.1$. (c,d) Real and imaginary
part for the energy of the photonic modes when scanning $\alpha$ from 0 to
0.3. (e,f) $\kappa$ and $FF$ when scanning $FF$ from 0 to 0.3.

###  C.1 Numerical RCWA simulation versus developed effective theory

To validate the effective Dirac theory for the guided photonic modes (2) and
the polaritonic modes (14), we perform RCWA simulations of the grating of
period $a=243$ nm, filling fraction $FF=0.37$ and 110 nm of etching depth. The
numerical results of photonic and polaritonic modes, together with the
fittings using the effective Hamiltonians are presented in Fig. 7. It shows
that both simulated photonic and polaritonic modes are perfectly followed by
our analytical model with: $\hbar\gamma=0.19$ meV, $\hbar\kappa=1.75$ meV,
$\hbar v=56.61$ meV $\mu$m, $\hbar\Omega=3.2$ meV,
$\hbar\omega^{(0)_{X}}=2.41$ meV, $\hbar\gamma_{nr}=0.18$ meV and $\varphi=0$.
Importantly, the vanishing of the intensity in the lower antisymmetric
polariton branch $\lvert L,-\rangle$ at $k=0$ in Fig. 7(c) confirms the
infinite radiative lifetime of these polaritons that are inherited from the
photonic BIC. Therefore the finite linewidth of $\lvert L,-\rangle$ at $k=0$
[see Fig. 7(g)] is purely nonradiative and is inherited by the excitonic
component. These results are in good agreement with the experimental
observations in [43].

###  C.2 Varying the exciton-photon energy detuning

By scanning the period $a$ of the same design (i.e. $FF$, etching depth), the
energy detuning $\hbar\omega^{(0)}_{X}$ between the exciton energy $E_{X}$ and
the mid gap of photonic modes
$E_{0}=\text{Re}{\left(\hbar\omega_{+}+\hbar\omega_{-}\right)/2}$—different
from the guided-lossy photonic detuning $\Delta$ in Appendix A—can be freely
varied, with all other parameters unchanged. This is clearly evidenced in Fig.
8 that reports the RCWA results of the energy of $\lvert\pm\rangle$ symmetric
and antisymmetric photonic modes at $k=0$ when scanning the period $a$ from
241 nm to 244.5 nm.

###  C.3 Tuning the diffractive coupling $\kappa$ and implementing $\pi$-jump
to $\varphi$

As previously reported [32, 29], the value of $\kappa$ can be continuously
tuned by modifying the filling fraction $FF$. To illustrate this effect in the
case of our design, we keep the period $a=250$ nm fixed and monitoring the
modification of the photonic modes when scanning $FF$ from 0.4 to 0.86. This
scanning induces a band-inversion to the photonic modes [32]. Indeed,
simulated photonic dispersions with $FF=$0.52, 0.62 and 0.8 are shown in Fig.
9(b). These results represent respectively three case: i) photonic BIC in the
lower band, corresponding to non-zero $\kappa$ and $\varphi=0$; ii) gap
closing and formation of Exceptional Points [32], corresponding to $\kappa=0$;
iii) photonic BIC in the upper band, corresponding to non-zero $\kappa$ and
$\varphi=\pi$. The values of $\kappa$ and $\varphi$ as the function of $FF$
are extracted from the real-part [see Fig. 9(c)] and imaginary-part [see
Fig.9(d)] of the photonic modes. These results, presented in Figs. 9(e) and
9(f), show that $\kappa$ is continuously tuned between 0 and 5 meV; and
$\varphi$ undergoes a $\pi$ jump when $\kappa=0$ ($FF$=0.62).

###  C.4 Tuning continuously the phase $\varphi$

To obtain $\varphi$ that is not a multiple of $\pi$, we break the in-plan
mirror symmetry $-x\rightarrow x$ of the grating. This is achieved by
employing double-period the design [72]: each unitcell now consists of two
sub-cell of period $(1+\alpha)a$ and $(1-\alpha)a$ with $\alpha$ being the
symmetry-breaking coefficient [see Fig. 10(a)]. Consequently, the phase-shift
$\varphi$ can be continuously tuned by changing $\alpha$.

To illustrate this effect, we keep $a=250$ nm, $FF=0.45$, and monitoring the
modification of the photonic modes when scanning $\alpha$ from 0 to 0.25. As
expected, a non-zero $\alpha$ turn a BIC into quasi-BIC [see Fig. 10(b)]: the
farfield at $k=0$ of the quasi-BIC is not-vanished and the photonic band
exhibits non-zero linewidth. The values of $\kappa$ and $\varphi$ as the
function of $\alpha$ are extracted from the real-part [see Fig. 10(c)] and
imaginary-part [see Fig. 10(d)] of the photonic modes. These results,
presented in Figs. 10(e) and 10(f), show that while $\kappa$ only undergoes
small variation, $\varphi$ is continuously tuned between 0 and $\pi/2$.
Interestingly, for $\lvert L,-\rangle$ is still almost lossless, thus being
quasi-BIC, until $\alpha$ exceeds 0.15 [see Fig. 10(d)], that corresponds to
$\varphi<0.15\pi$ [see Fig. 10(f)].

##  Appendix D Numerical mean field methods

Equations (15) and (16) form a coupled nonlinear systems of equation to be
solved for $n_{R}$ and $\Psi$. Numerical simulations are performed using fast
Fourier transform spectral methods in space and an explicit Runge-Kutta
4th-5th order formula, the Dormand-Prince pair, in time. We apply damped
boundary conditions in real and reciprocal space to avoid periodic boundary
effects coming from the spectral methods, and always start from random white
noise initial conditions. The gridpoint spacing $\Delta x$ is chosen small
enough to encompass the necessary features in momentum space. In our case,
$\Delta x=4$ $\mu$m was sufficient. A variable timestep is set by the MATLAB®
ode45 solver to reach the desired accuracy.

## References

  * Carusotto and Ciuti [2013] I. Carusotto and C. Ciuti, Quantum fluids of light, Rev. Mod. Phys. 85, 299 (2013).
  * Tosi  _et al._ [2012] G. Tosi, G. Christmann, N. G. Berloff, P. Tsotsis, T. Gao, Z. Hatzopoulos, P. G. Savvidis, and J. J. Baumberg, Sculpting oscillators with light within a nonlinear quantum fluid, Nature Physics 8, 190 (2012).
  * Pickup  _et al._ [2020] L. Pickup, H. Sigurdsson, J. Ruostekoski, and P. G. Lagoudakis, Synthetic band-structure engineering in polariton crystals with non-hermitian topological phases, Nature Communications 11, 4431 (2020).
  * Pieczarka  _et al._ [2021] M. Pieczarka, E. Estrecho, S. Ghosh, M. Wurdack, M. Steger, D. W. Snoke, K. West, L. N. Pfeiffer, T. C. H. Liew, A. G. Truscott, and E. A. Ostrovskaya, Topological phase transition in an all-optical exciton-polariton lattice, Optica 8, 1084 (2021).
  * Alyatkin  _et al._ [2021] S. Alyatkin, H. Sigurdsson, A. Askitopoulos, J. D. Töpfer, and P. G. Lagoudakis, Quantum fluids of light in all-optical scatterer lattices, Nature Communications 12, 5571 (2021).
  * Sanvitto and Kéna-Cohen [2016] D. Sanvitto and S. Kéna-Cohen, The road towards polaritonic devices, Nature Materials 15, 1061 (2016).
  * Jiang  _et al._ [2022] Z. Jiang, A. Ren, Y. Yan, J. Yao, and Y. S. Zhao, Exciton-polaritons and their Bose–Einstein condensates in organic semiconductor microcavities, Advanced Materials 34, 2106095 (2022).
  * Terças  _et al._ [2014] H. Terças, H. Flayac, D. D. Solnyshkov, and G. Malpuech, Non-abelian gauge fields in photonic cavities and photonic superfluids, Phys. Rev. Lett. 112, 066402 (2014).
  * Gianfrate  _et al._ [2020] A. Gianfrate, O. Bleu, L. Dominici, V. Ardizzone, M. De Giorgi, D. Ballarini, G. Lerario, K. W. West, L. N. Pfeiffer, D. D. Solnyshkov, D. Sanvitto, and G. Malpuech, Measurement of the quantum geometric tensor and of the anomalous hall drift, Nature 578, 381 (2020).
  * Su  _et al._ [2021] R. Su, E. Estrecho, D. Biegańska, Y. Huang, M. Wurdack, M. Pieczarka, A. G. Truscott, T. C. H. Liew, E. A. Ostrovskaya, and Q. Xiong, Direct measurement of a non-hermitian topological invariant in a hybrid light-matter system, Science Advances 7, eabj8905 (2021).
  * Polimeno  _et al._ [2021] L. Polimeno, A. Fieramosca, G. Lerario, L. D. Marco, M. D. Giorgi, D. Ballarini, L. Dominici, V. Ardizzone, M. Pugliese, C. T. Prontera, V. Maiorano, G. Gigli, C. Leblanc, G. Malpuech, D. D. Solnyshkov, and D. Sanvitto, Experimental investigation of a non-abelian gauge field in 2D perovskite photonic platform, Optica 8, 1442 (2021).
  * Kokhanchik  _et al._ [2021] P. Kokhanchik, H. Sigurdsson, B. Piętka, J. Szczytko, and P. G. Lagoudakis, Photonic berry curvature in double liquid crystal microcavities with broken inversion symmetry, Phys. Rev. B 103, L081406 (2021).
  * Łempicka Mirek  _et al._ [2022] K. Łempicka Mirek, M. Król, H. Sigurdsson, A. Wincukiewicz, P. Morawiak, R. Mazur, M. Muszyński, W. Piecek, P. Kula, T. Stefaniuk, M. Kamińska, L. D. Marco, P. G. Lagoudakis, D. Ballarini, D. Sanvitto, J. Szczytko, and B. Piętka, Electrically tunable berry curvature and strong light-matter coupling in liquid crystal microcavities with 2D perovskite, Science Advances 8, eabq7533 (2022).
  * Lovett  _et al._ [2023] S. Lovett, P. M. Walker, A. Osipov, A. Yulin, P. U. Naik, C. E. Whittaker, I. A. Shelykh, M. S. Skolnick, and D. N. Krizhanovskii, Observation of Zitterbewegung in photonic microcavities, Light: Science & Applications 12, 126 (2023).
  * Jacqmin  _et al._ [2014] T. Jacqmin, I. Carusotto, I. Sagnes, M. Abbarchi, D. D. Solnyshkov, G. Malpuech, E. Galopin, A. Lemaître, J. Bloch, and A. Amo, Direct observation of Dirac cones and a flatband in a honeycomb lattice for polaritons, Phys. Rev. Lett. 112, 116402 (2014).
  * Nalitov  _et al._ [2015] A. V. Nalitov, G. Malpuech, H. Terças, and D. D. Solnyshkov, Spin-orbit coupling and the optical spin hall effect in photonic graphene, Phys. Rev. Lett. 114, 026803 (2015).
  * Yi and Karzig [2016] K. Yi and T. Karzig, Topological polaritons from photonic Dirac cones coupled to excitons in a magnetic field, Phys. Rev. B 93, 104303 (2016).
  * Klembt  _et al._ [2018] S. Klembt, T. H. Harder, O. A. Egorov, K. Winkler, R. Ge, M. A. Bandres, M. Emmerling, L. Worschech, T. C. H. Liew, M. Segev, C. Schneider, and S. Höfling, Exciton-polariton topological insulator, Nature 562, 552 (2018).
  * Milićević  _et al._ [2019] M. Milićević, G. Montambaux, T. Ozawa, O. Jamadi, B. Real, I. Sagnes, A. Lemaître, L. Le Gratiet, A. Harouri, J. Bloch, and A. Amo, Type-III and tilted Dirac cones emerging from flat bands in photonic orbital graphene, Phys. Rev. X 9, 031010 (2019).
  * Liu  _et al._ [2020] W. Liu, Z. Ji, Y. Wang, G. Modi, M. Hwang, B. Zheng, V. J. Sorger, A. Pan, and R. Agarwal, Generation of helical topological exciton-polaritons, Science 370, 600 (2020).
  * Li  _et al._ [2021] M. Li, I. Sinev, F. Benimetskiy, T. Ivanova, E. Khestanova, S. Kiriushechkina, A. Vakulenko, S. Guddala, M. Skolnick, V. M. Menon, D. Krizhanovskii, A. Alù, A. Samusev, and A. B. Khanikaev, Experimental observation of topological Z2 exciton-polaritons in transition metal dichalcogenide monolayers, Nature Communications 12, 4425 (2021).
  * Wang  _et al._ [2022] J. Wang, Y. Peng, H. Xu, J. Feng, Y. Huang, J. Wu, T. C. H. Liew, and Q. Xiong, Controllable vortex lasing arrays in a geometrically frustrated exciton–polariton lattice at room temperature, National Science Review 10, nwac096 (2022).
  * Lu  _et al._ [2015] L. Lu, Z. Wang, D. Ye, L. Ran, L. Fu, J. D. Joannopoulos, and M. Soljačić, Experimental observation of Weyl points, Science 349, 622 (2015).
  * Lee  _et al._ [2021] K. Y. Lee, K. W. Yoo, Y. Choi, G. Kim, S. Cheon, J. W. Yoon, and S. H. Song, Topological guided-mode resonances at non-hermitian nanophotonic interfaces, Nanophotonics 10, 1853 (2021).
  * Chen  _et al._ [2023] K. Chen, F. Komissarenko, D. Smirnova, A. Vakulenko, S. Kiriushechkina, I. Volkovskaya, S. Guddala, V. Menon, A. Alù, and A. B. Khanikaev, Photonic Dirac cavities with spatially varying mass term, Science Advances 9, eabq4243 (2023).
  * Schmidt  _et al._ [2015] M. Schmidt, V. Peano, and F. Marquardt, Optomechanical Dirac physics, New Journal of Physics 17, 023025 (2015).
  * Guddala  _et al._ [2021] S. Guddala, F. Komissarenko, S. Kiriushechkina, A. Vakulenko, M. Li, V. M. Menon, A. Alù, and A. B. Khanikaev, Topological phonon-polariton funneling in midinfrared metasurfaces, Science 374, 225 (2021).
  * In  _et al._ [2022] C. In, U. J. Kim, and H. Choi, Two-dimensional Dirac plasmon-polaritons in graphene, 3D topological insulator and hybrid systems, Light: Science & Applications 11, 313 (2022).
  * Ardizzone  _et al._ [2022] V. Ardizzone, F. Riminucci, S. Zanotti, A. Gianfrate, M. Efthymiou-Tsironi, D. G. Suàrez-Forero, F. Todisco, M. De 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).
  * 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).
  * Hwang  _et al._ [2022] M.-S. Hwang, K.-Y. Jeong, J.-P. So, K.-H. Kim, and H.-G. Park, Nanophotonic nonlinear and laser devices exploiting bound states in the continuum, Communications Physics 5, 106 (2022).
  * 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, Photon. Res. 8, A91 (2020).
  * 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).
  * Grudinina  _et al._ [2023] A. Grudinina, M. Efthymiou-Tsironi, V. Ardizzone, F. Riminucci, M. D. Giorgi, D. Trypogeorgos, K. Baldwin, L. Pfeiffer, D. Ballarini, D. Sanvitto, and N. Voronova, Collective excitations of a bound-in-the-continuum condensate, Nature Communications 14, 3464 (2023).
  * 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] .
  * 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 (2009).
  * Kravtsov  _et al._ [2020] V. Kravtsov, E. Khestanova, F. A. Benimetskiy, T. Ivanova, A. K. Samusev, I. S. Sinev, D. Pidgayko, A. M. Mozharov, I. S. Mukhin, M. S. Lozhkin, Y. V. Kapitonov, A. S. Brichkin, V. D. Kulakovskii, I. A. Shelykh, A. I. Tartakovskii, P. M. Walker, M. S. Skolnick, D. N. Krizhanovskii, and I. V. Iorsh, Nonlinear polaritons in a monolayer semiconductor coupled to optical bound states in the continuum, Light: Science & Applications 9, 56 (2020).
  * Koksal  _et al._ [2021] O. Koksal, M. Jung, C. Manolatou, A. N. Vamivakas, G. Shvets, and F. Rana, Structure and dispersion of exciton-trion-polaritons in two-dimensional materials: Experiments and theory, Phys. Rev. Res. 3, 033064 (2021).
  * 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).
  * Maggiolini  _et al._ [2023] E. Maggiolini, L. Polimeno, F. Todisco, A. Di Renzo, B. Han, M. De Giorgi, V. Ardizzone, C. Schneider, R. Mastria, A. Cannavale, M. Pugliese, L. De Marco, A. Rizzo, V. Maiorano, G. Gigli, D. Gerace, D. Sanvitto, and D. Ballarini, Strongly enhanced light–matter coupling of monolayer WS2 from a bound state in the continuum, Nature Materials 22, 964 (2023).
  * Weber  _et al._ [2023] T. Weber, L. Kühner, L. Sortino, A. Ben Mhenni, N. P. Wilson, J. Kühne, J. J. Finley, S. A. Maier, and A. Tittl, Intrinsic strong light-matter coupling with self-hybridized bound states in the continuum in van der waals metasurfaces, Nature Materials 22, 970 (2023).
  * 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 20, 2113 (2020).
  * 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).
  * Wang  _et al._ [2023a] Y. Wang, J. Tian, M. Klein, G. Adamo, S. T. Ha, and C. Soci, Directional Emission from Electrically Injected Exciton–Polaritons in Perovskite Metasurfaces, Nano Letters 23, 4431 (2023a).
  * Leykam  _et al._ [2017] D. Leykam, K. Y. Bliokh, C. Huang, Y. D. Chong, and F. Nori, Edge modes, degeneracies, and topological numbers in non-hermitian systems, Phys. Rev. Lett. 118, 040401 (2017).
  * Sedov  _et al._ [2018] E. S. Sedov, Y. G. Rubo, and A. V. Kavokin, $\mathit{Zitterbewegung}$ of exciton-polaritons, Phys. Rev. B 97, 245312 (2018).
  * 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).
  * Low  _et al._ [2017] T. Low, A. Chaves, J. D. Caldwell, A. Kumar, N. X. Fang, P. Avouris, T. F. Heinz, F. Guinea, L. Martin-Moreno, and F. Koppens, Polaritons in layered two-dimensional materials, Nature Materials 16, 182 (2017).
  * Joannopoulos  _et al._ [2011] J. Joannopoulos, S. Johnson, J. Winn, and R. Meade,  _Photonic Crystals: Molding the Flow of Light - Second Edition_ (Princeton University Press, 2011).
  * [51] Introducing a transverse momentum $k_{y}$ will induce a polarization mismatch between the two counterpropagating modes. This adds a supplemental coefficient of $\vec{u}_{f}\cdot\vec{u}_{b}$ to the radiative coupling strength, with $\vec{u}_{f,b}$ being the polarization vectors of forward and backward propagating guided modes [29, Mermet-Lyaudoz2023].
  * Liu  _et al._ [2019] Z. Liu, Y. Xu, Y. Lin, J. Xiang, T. Feng, Q. Cao, J. Li, S. Lan, and J. Liu, High-$q$ quasibound states in the continuum for nonlinear metasurfaces, Phys. Rev. Lett. 123, 253901 (2019).
  * Sun  _et al._ [2022] K. Sun, H. Jiang, D. A. Bykov, V. Van, U. Levy, Y. Cai, and Z. Han, 1D quasi-bound states in the continuum with large operation bandwidth in the $\omega\sim k$ space for nonlinear optical applications, Photon. Res. 10, 1575 (2022).
  * Wang  _et al._ [2023b] K. Wang, T. Gu, D. A. Bykov, X. Zhang, and L. Qian, Tunable nanolaser based on quasi-bic in a slanted resonant waveguide grating, Opt. Lett. 48, 4121 (2023b).
  * Głowadzka  _et al._ [2021] W. Głowadzka, M. Wasiak, and T. Czyszanowski, True- and quasi-bound states in the continuum in one-dimensional gratings with broken up-down mirror symmetry, Nanophotonics 10, 3979 (2021).
  * Liu  _et al._ [2023a] J. Liu, C. Chen, X. Li, J. Li, D. Dong, Y. Liu, and Y. Fu, Tunable dual quasi-bound states in continuum and electromagnetically induced transparency enabled by the broken material symmetry in all-dielectric compound gratings, Opt. Express 31, 4347 (2023a).
  * 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).
  * 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).
  * 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).
  * Nigro and Gerace [2023] D. Nigro and D. Gerace, Theory of exciton-polariton condensation in gap-confined eigenmodes, Phys. Rev. B 108, 085305 (2023).
  * Töpfer  _et al._ [2020] J. D. Töpfer, H. Sigurdsson, S. Alyatkin, and P. G. Lagoudakis, Lotka-volterra population dynamics in coherent and tunable oscillators of trapped polariton condensates, Phys. Rev. B 102, 195428 (2020).
  * Sun  _et al._ [2018] Y. Sun, Y. Yoon, S. Khan, L. Ge, M. Steger, L. N. Pfeiffer, K. West, H. E. Türeci, D. W. Snoke, and K. A. Nelson, Stable switching among high-order modes in polariton condensates, Phys. Rev. B 97, 045303 (2018).
  * Mrejen  _et al._ [2019] M. Mrejen, L. Yadgarov, A. Levanon, and H. Suchowski, Transient exciton-polariton dynamics in WSe2 by ultrafast near-field imaging, Science Advances 5, eaat9618 (2019).
  * Tománek and Grmela [2008] P. Tománek and L. Grmela, Optics of nano-objects, in  _Eighth International Conference on Correlation Optics_, Vol. 7008, edited by M. Kujawinska and O. V. Angelsky, International Society for Optics and Photonics (SPIE, 2008) p. 70081F.
  * 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 5, 805 (2009).
  * Antón  _et al._ [2014] C. Antón, D. Solnyshkov, G. Tosi, M. D. Martín, Z. Hatzopoulos, G. Deligeorgis, P. G. Savvidis, G. Malpuech, and L. Viña, Ignition and formation dynamics of a polariton condensate on a semiconductor microcavity pillar, Phys. Rev. B 90, 155311 (2014).
  * Estrecho  _et al._ [2019] E. Estrecho, T. Gao, N. Bobrovska, D. Comber-Todd, M. D. Fraser, M. Steger, K. West, L. N. Pfeiffer, J. Levinsen, M. M. Parish, T. C. H. Liew, M. Matuszewski, D. W. Snoke, A. G. Truscott, and E. A. Ostrovskaya, Direct measurement of polariton-polariton interaction strength in the thomas-fermi regime of exciton-polariton condensation, Phys. Rev. B 100, 035306 (2019).
  * Pieczarka  _et al._ [2020] M. Pieczarka, D. Poletti, C. Schneider, S. Höfling, E. A. Ostrovskaya, G. Sęk, and M. Syperek, Observation of gain-pinned dissipative solitons in a microcavity laser, APL Photonics 5, 086103 (2020).
  * Misra and Sudarshan [2008] B. Misra and E. C. G. Sudarshan, The Zeno’s paradox in quantum theory, Journal of Mathematical Physics 18, 756 (2008).
  * Liu  _et al._ [2023b] Q. Liu, W. Liu, K. Ziegler, and F. Chen, Engineering of zeno dynamics in integrated photonics, Phys. Rev. Lett. 130, 103801 (2023b).
  * Liu and Fan [2012] V. Liu and S. Fan, S4 : A free electromagnetic solver for layered periodic structures, Computer Physics Communications 183, 2233 (2012).
  * Nguyen  _et al._ [2018] H. S. Nguyen, F. Dubois, T. Deschamps, S. Cueff, A. Pardon, J.-L. Leclercq, C. Seassal, X. Letartre, and P. Viktorovitch, Symmetry breaking in photonic crystals: On-demand dispersion from flatband to Dirac cones, Phys. Rev. Lett. 120, 066102 (2018).