Dirac exciton–polariton condensates in photonic crystal gratings

Abstract Bound states in the continuum have recently been utilized in photonic crystal gratings to achieve strong coupling and ultralow threshold condensation of 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 polaritonic modes and describe their mean-field condensation dynamics facilitated by the interplay between protection from the radiative continuum and negative-mass optical trapping. Our theory accounts for tunable grating parameters giving full control over the diffractive coupling properties between guided polaritons and the radiative continuum, unexplored for polariton condensates. In particular, we discover stable cyclical condensate solutions mimicking a driven-dissipative analog of the zitterbewegung effect characterized by coherent superposition of ballistic and trapped polariton waves. We clarify important distinctions between the polariton nearfield and farfield explaining recent experiments on the emission characteristics of these long lived nonlinear Dirac polaritons.


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] and readout, high interaction strengths, and the wide breadth of materials permitting polariton condensation at elevated temperatures [4], [5].These lightmass 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 [6]- [11] and patterned photonic structures [12]- [18].Such materials, hosting associated Dirac cones, offer valuable insight to a plethora of exotic phenomena such as quantum Hall physics [7], nontrivial topological phases [14], Weyl semimetals [19], and relativistic trapping [20], [21] while supplemented with strong polariton nonlinearities.Moreover, alternative neighbouring platforms for exploration into light-matter Dirac physics involve phonon-polaritons [22], [23] and plasmon-polaritons [24].
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 Figure 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 [25].A parallel study demonstrated the impressive tunability such grated waveguides offered over the polariton properties [26].In these systems, strong coupling is facilitated by the photonic structure's protection from the continuum [27], [28] which allows photons to survive long enough to form polariton states [29], [30].The initial experiment [25] was soon followed with fascinating results on the behaviour of the fluid's elementary excitations [31] and demonstration of macroscopic hybridization between coupled BIC condensates [32].Besides III-V semiconductor photonic crystals [33], other platforms able to show BIC-facilitated strong exciton-photon coupling consist of dry transfer deposited transition metal dichalcogenide monolayers such as MoSe 2 [34], [35] or WS 2 [36]- [38]; or spin coated hybrid organicinorganic perovskites [39]- [42].
Inspired by these rapid developments in BIC facilitated polariton condensation, and the surging interest to simulate nonlinear and non-Hermitian relativistic physics in an optically addressable setting [43], 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 negativemass optical trapping of Dirac polaritons.Our theory is in excellent agreement with recent experiments on polariton crystal gratings under nonresonant excitation [25], [32], [44] and uncovers 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 [11], [45] limit cycle, forming a strange mixture of confined low momentum negativemass 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 [1].

Photonic modes in optical grating
We consider a waveguide stack consisting of multiple periodic layers of GaAs quantum wells separated by AlGaAs barriers along the normal z-direction [see Figure 1(a)].On its upper surface, the waveguide is patterned with a one-dimensional subwavelength grating with a period a along the x-direction and filling factor FF, forming a onedimensional (1D) photonic crystal slab.This kind of a sample has been demonstrated in [25], [26], [32], however our developed theory is also suitable for photonic crystals deposited with 2D materials [34]- [38], [46] or perovskites [39]- [42].
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 | | e iqx ⟩ with momentum q and frequency  q .As is well known [47], the grating with period a can be modeled by a periodic potential acting on the photonic modes, where u(x) is a periodic function of period a, and (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 integer number of the primitive reciprocal lattice number K = 2π∕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  ±K .Therefore, the relevant photonic modes for exciton-photon coupling correspond to those with wavevectors q = ±K + k with k ≪ K.In this range of frequencies, the periodic potential couples nearly degenerate guided modes | | e i(K+k)x ⟩ and | | e i(−K+k)x ⟩ , giving rise to an effective Dirac Hamiltonian describing the dynamics of the electromagnetic waves.Remarkably, the periodic potential not only couples the two modes | | e ±iKx ⟩ together, but also couples them with lossy modes residing at normal incidence q = 0.This renders the two guided modes | | e ±iKx ⟩ eventually lossy; see Figure 1(b).In total, the dynamics of the electromagnetic waves in the photonic crystal grating can be described by a lossy Dirac Hamiltonian, with where U p is pth Fourier coefficient of u(x), i.e., U p = ⟨1|u(x) | | e −i pKx ⟩ .We have also linearised the dispersal relation around K so that  k±K =  K ± vk and ignore the constant term.A more detailed derivation of this effective low-momentum photonic Hamiltonian is given in the Methods.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  and the loss exchange mechanism of strength e i via the radiative continuum [29], [48].Notice that we adopt here the convention of using e −it to describe the temporal oscillation, hence losses are given by the negative imaginary component of the dispersion relation.Following Eq. ( 3), the parameters ,  and  are dictated by the two first Fourier components U 1 and U 2 of the periodic modulation u(x).The diffractive coupling  >  is the main parameter responsible for the bandgap opening at k = 0 [see solid curves in Figure 1(c)].Its value can be engineered by tuning the filling fraction of the grating [29]; see Methods for numerical values of  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 In the absence of losses one recovers the Dirac dispersal relation with the effective photon Dirac mass = .The real and imaginary parts of the photon frequencies are plotted in Figure 1(c and d) for  = 0 showing a zero-loss BIC at k = 0 in the antisymmetric branch.This can be understood from, For small loss,  ≪ , and k → 0, we have lim k→0 Im( ± ) ≈ −(1 ± cos ). ( Conversely, when  = π the BIC switches branches.In both case, as long as  is an integer multiple of π due to the presence of the mirror symmetry x → −x, a formation of BIC at Γ point is guaranteed.This BIC is therefore of a symmetryprotected nature [27]. Breaking the mirror symmetry will relax the aforementioned constraint on .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 [49]- [53].We note that a similar form of (2) has been previously reported to describe the formation of symmetry-protected BICs [25], [29], however our work provides the first effective Hamiltonian for the general case where the in-plane mirrorsymmetry can be broken.Importantly, for realistic grating structures with lateral symmetry breaking design, a fine tuning of  from 0 to π∕2 can be achieved (see Methods).
To start with, we will explore the case of a grating with lateral mirror symmetry that has  = 0.At a later stage we will relax this constraint and explore Dirac-polariton BIC condensation when 0 <  < π.

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 Diracpolaritons.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 Ω (also known as the exciton-photon Rabi frequency) [29], [54].In our photonic grating, the lossy mode and BIC mode exhibit opposite parities of mirror symmetry x → −x, with one being symmetric and the other being antisymmetric.Their near-field distributions are periodic functions with a period of the grating (a ≈ 300 nm), with the anti-nodes of the lossy mode field being the nodes of the BIC field, and vice versa.Therefore, excitons, tightly localized within the Bohr radius of about 10 nm, can only couple efficiently to either the lossy mode or the BIC mode due to the spatial overlap requirement.We then follow the approximation proposed by Lu et al. [29] to divide excitons into two groups: one that only couples to the lossy mode and another that only couples to the BIC mode with the same coupling strength.Therefore, the strong coupling picture is dictated by a 4 × 4 Hamiltonian, given by [29]: Here,  X =  (0) X − i nr where  (0) X and  −1 nr 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.We highlight that this simple model successfully matches the numerical simulation results of polaritonic branches (see Section 6.4).
In a more precise consideration, one should include a third group of excitons that are not localized at the antinodes of both photonic modes, thus do not undergo the strong coupling regime.This corresponds to the observation of uncoupled excitons shown in the simulations of Section 6.4.Since these uncoupled excitons do not affect the strong coupling physics, we do not include them in our effective theory.
The eigenmodes of the Hamiltonian (8) are referred to as upper |U, ±⟩ and lower |L, ±⟩ symmetric-antisymmetric polaritons with a dispersion relation, which is plotted in Figure 2(a) and (b).Notice that due to the exciton losses, the BIC now becomes a quasi-BIC with finite losses.
We will assume that the upper polariton branches  U,± are far away in energy and weakly populated and thus only focus on the lower polariton branches  L,± around k = 0 where condensation preferentially takes place [25], [32], [55].The lower branches can be approximated by considering first the coupling of forward-| | e iKx ⟩ and backwardpropagating | | e −iKx ⟩ photons to excitons, leading to lower-  propagating polaritons.This is followed up by the photonic diffractive coupling mechanism evaluated at small momenta.In other words, polariton modes are first calculated in the grating free waveguide then we perturbatively introduce the grating back (diffractive coupling).The lower polariton dispersion can then be written, (10) where the first term in (10) corresponds simply to an overall complex energy shift due to the light-matter coupling which is written, The renormalized light-matter velocity ṽ and the photon Hopfield coefficient of forward and backward propagating lower polaritons around k = 0 are given by, The form of Eq. ( 10) implies that, in the truncated basis of lower forward | | L, e iKx ⟩ and backward | | L, e −iKx ⟩ polaritons, we can describe the system with the following massive non-Hermitian Dirac operator, ) with new symmetric and antisymmetric lower polariton energies  L,± .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  ≪ Ω ≪  (0) X then Eq. ( 10) remains accurate and implies that waveguided TE polaritons behave approximately as Dirac particles with renormalized velocity ṽ and gap opening.
In Figure 2(c) and (d) we compare our approximated dispersion (10) (dashed lines) with the lower polariton dispersion relation coming from (9) (solid lines) for values extracted from RCWA simulations, Ω∕ = 1.8, ℏ (0)  X ∕ = 5.5, and  nr = , 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 (see Methods).The above underpins the feasibility in creating photonic samples that permit study of Dirac polariton quasi-BIC physics.

Mean-field formalism
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 [25], [55].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 [56] 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 nonequilibrium Bose-Einstein condensation into a single quantum state [1], marked by a threshold power Pth.Following the mean-field theory [56], [57] and our approximate single particle Dirac operator (14), the condensate can be described by a two-component macroscopic wave function, or an order parameter, Ψ(x, t) = ( +1 ,  −1 ) T where ±1 denotes the lower forward-and backward propagating polaritons, respectively.Note that subscripts (±) 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 The first term corresponds to the single particle dynamics (14) in real space obtained by substituting k → −i x .In the second term the polaritons are assumed to interact via short-range interaction described phenomenologically by the repulsive non-linear term gΨ † Ψ.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 .
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 ∕2 2 describes the continuous wave nonresonant Gaussian pump of waist , and Γ R is the average reservoir exciton redistribution and decay rate.
For large waist , the pumping can be considered to be uniform, P(x) = P 0 .In this case, the condensation threshold is given by P th = −2Γ R max{Im[ L,± ]}∕R where Im[ L,± ] < 0 and the maximum is taken over ±.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 Methods.
Below threshold one has Ψ † Ψ = 0 in the long time limit and the reservoir converges to the steady state n R = P(x)∕Γ R .We can then define a pump induced potential term acting on the Dirac polaritons, The above expression gives a good estimate for the optical potential 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], [58].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 [25], [32], [57].

Mode dropping
Figure 3 shows the total condensate density  = Ψ † Ψ in (upper row) real space and (lower row) energy-resolved momentum space for increasing pump power density P 0 .At low powers above threshold [Figure 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 [57].
Increasing the power we observe monotonic blueshift of the condensate level [compare Figure 3(a) and (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 [59] in the same branch.Such cyclical solutions [see Figure 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 lower energies [see Figure 3(d)].In this sense, the condensate "drops" from one trap mode to the next, as predicted by Nigro et al. [57].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 [25].Energetically, this behaviour is in sharp contrast to optically trapped polariton condensates in planar cavities [58]- [60] where stronger pumping results in a condensate dropping into lower order trap modes until it reaches the ground state.

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  = 0 which leads to a quasi-BIC in the lower (symmetric) branch in Figure 2(c) and (d).If  = π then the quasi-BIC would instead form in the upper branch [29].
We perform a scan across both  and the pump power P 0 and investigate the difference between symmetric and antisymmetric condensate occupation where and |Ψ ± ⟩ are the symmetric and antisymmetric eigenstates of the potential-free Dirac operator (14).In this sense, the quantity Δ  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 Δ  < 0 then most of the condensate forms in the lower branch, whereas if Δ  > 0 then the condensate forms in the upper branch.
The results are shown in Figure 4 for two different sizes of pump spots (FWHM = 50 and 200 μm) where each pixel corresponds to the spatiotemporal average of ⟨Δ  ⟩ over the entire simulation grid and integration time.The results show that for  ∼ 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 [25].Interestingly, for smaller pump spots in Figure 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  ∼ π.This implies that the optical trapping mechanism can be more efficient in reducing transverse (xdirection) 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 μm in Figure 4(b).
The much wider pump spot imposes a weaker confinement compared to the narrow spot.Indeed, now around  ∼ π 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.
At the interface of such qualitatively different condensate solutions [i.e., dark and bright regions in Figure 4  more exotic patterns might appear.We propose that by tuning the size of the trap, the grating pitch (), 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 [45].Interference between the upper and lower branches causes the center-of-mass of the condensate, ⟨x⟩ = ∫ Ψ † xΨdx∕ ∫ Ψ † Ψdx, to jitter in time (see red oscillating curve).In Figure 5(a-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 Figure 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, ∝ cos 2 (πΔEt∕h) (not to be confused with the light-matter Rabi frequency Ω).Here h is the Planck's constant and ΔE is the energy splitting between the ballistic The ballistic nature of the upper branch fluid manifests in its more delocalized nature and more rapidly varying phase profile along the x-direction.

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 [25], [32].We remind the reader that the polariton state is explicitly related to the emitted photon state vector Ψ ∝ Φ 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 [61].
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. [62]: '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 [63]- [66], 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 → −i x and denoting the two-component photon state vector as Φ = ( +1 ,  −1 ) T where  ±1 are the coefficients of the forward and backward propagating photons | | e i(k±K)x ⟩ , ) Then using the dynamical equation  t Φ(x, t) = −iH ph Φ(x, t) and its conjugation, we arrive directly at the continuity equation for the photon field, where  z is the third Pauli matrix.From the continuity equation, one infers, as usual, that Φ † Φ describe the nearfield intensity, and Φ †  z Φ describes the photon density current.The last term, Φ † ( 1 e i e −i 1 ) is then identified as losses by means of radiation into the farfield.In this way, the expression for nearfield and farfield intensities has 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 the Methods section 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 Figure 3. Notably, in [25] only the farfield was measured showing emission profiles which agree very well with our results.

Conclusions and discussion
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 manybody picture through the mean-field formalism.We propose a generalized Gross-Pitaevskii model to describe BICfacilitated condensation of polaritons into pump-induced optical traps.Our findings are in excellent agreement with recent experimental observations on multi-quantum-well structures [25] and are applicable to other forms of optically active materials such as transition metal dichalcogenide monolayers including MoSe 2 [34], [35] or WS 2 [36]- [38]; or hybrid organic-inorganic perovskites [39]- [42].
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 mass) and lower-branch (negative mass) polaritons.The implications of such a hybrid quantum fluid are the two very different coupling mechanisms when neighbouring condensates are added.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 multicomponent nonlinear oscillators with contrasting coupling mechanisms with competing coherence scales and timescales of domain wall formation.

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  q are the plane waves | | e iqx ⟩ . 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  q =  −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 Figure 1(a).
The periodic modulation -with period a -of the refraction index in the x direction introduces a periodic potential u(x) = u(x + a) acting on the photons.As an expression of the Bragg reflection, the potential then couples modes where K = 2π∕a is the primitive reciprocal lattice vector and n is an integer.The relevant matrix elements are written , which are the nth 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 ±K.The relevant wavevectors are therefore of ⟩ are coupled by the potential u(x) with matrix elements For simplicity, in the following we consider the coupling to only the 0th 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 , a general wavefunction can be written as The evolution of the coefficients where  0 is the decay rate of the low-momentum mode | |  (0) (z)e ikx ⟩ , which is assumed to vary negligibly for small k.We have also linearized the dispersal relation so that  K+k ≈  K + vk and  −K+k ≈  K − vk, with  being the light velocity near K.By adjusting a global phase, one can also assume that  K = 0 and  0 can be then replaced by the difference in frequency Δ =  0 −  K .
Assuming that the decay rate of the lossy mode  0 is much faster than the matrix element U 2 , one can adiabatically eliminate φk 0 .This is done by solving φk 0 in terms of φk +1 and φk −1 as As the decaying of the lossy mode  0 is fast in comparison to the dynamics of the confined mode, one can make the Markovian approximation φk We have we also assumed  0 ≫ Δ in the last approximation.In the end, we then obtain ) Inserting ( 26) into (24), we obtain the evolution equation for φk

+1
and φk with the non-Hermitian Hamiltonian in momentum space where  1 is defined by . It is interesting to notice that the indirect loss rate  for the two modes | |  K (z)e ±iKx ⟩ is inversely proportional to the lossy rate  0 of the mode | |  (0) (z)e ikx ⟩ , indicating an analogy of the Zeno effect in quantum system [67].Indeed, a strong "measurement regime" that corresponds to a very leaky channel,  0 ≫ |U 1 |, will "freeze" the population of | |  K (z)e ±iKx ⟩ because  ≈ 0. A similar setup, combining lossless waveguides and a lossy one, has been recently proposed to demonstrate the optical Zeno effect [68].The Fourier coefficient U 2 in (28) is generally a complex number.
We denote U 2 = |U 2 |e i 2 and eliminate the phase  2 by the following unitary transformation, ( φk One then obtains the non-Hermitian Hamiltonian for  k +1 and  k with  = 2 1 −  2 .This is the effective Hamiltonian (2) introduced in the main text with  = |U 2 |.

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 Ψ inherits directly from the loss of the photonic component Φ. Therefore the farfield pattern of the polariton can be understood directly from its photonic components.
Generally, the effective wave function Φ(x) is a superposition of different wavevectors k, At the microscopic level, this is a plane wave of 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 where ⃗ u y denotes the polarization direction.Or using the real space representation Φ(x), we can write Recall that we are working in the regime where  k +1 and  k −1 are only significant at wavelengths k ≪ K. Equivalently,  +1 (x) and  −1 (x), which are referred to as spatial envelope functions, vary much slower than e ±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: 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  k 0 , which is given by ( 26) for a single k.With a superposition of different wavevectors (31), one has Explicitly in terms of the electric field, this corresponds to where we have again used  = 2 1 −  2 .The farfield intensity is then obtained as which coincides with (22).

Design and numerical simulations for pratical realizations
We propose realistic a stack of dielectric layers and quantum wells, based on the experimental works from [25], [32].The sample stack is composed of a waveguide core made of 12 GaAs quantum wells (QWs) 20 nm thick and 13 Al 0.4 Ga 0.6 As barriers 20 nm thick grown on an Al 0.8 Ga 0.2 As 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 S 4 package provided by the Fan Group at the Stanford Electrical Engineering Department [69].The refractive index of the QWs, barriers and the cladding are: n QWs = 3.547 + 0.0001i, n barrier = 3.3, n cladding = 3.063.The imaginary part in the refractive index of the QWs is 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 eV 2 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.

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 Figure 7.It shows that both simulated photonic and polaritonic modes are perfectly followed by our analytical model with: ℏ = ∕2 can be freely varied, with all other parameters unchanged.This is clearly evidenced in Figure 8 that reports the RCWA results of the energy of |±⟩ symmetric and antisymmetric photonic modes at k = 0 when scanning the period a from 241 nm to 244.5 nm.
6.4.2Tuning the diffractive coupling  and implementing -jump to : As previously reported [25], [29], the value of  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 [29].Indeed, simulated photonic dispersions with FF = 0.52, 0. imaginary-part [see Figure 9(d)] of the photonic modes.These results, presented in Figure 9(e) and (f), show that  is continuously tuned between 0 and 5 meV; and  undergoes a π jump when  = 0 (FF = 0.62).

Tuning continuously the phase 𝝋:
To obtain  that is not a multiple of π, we break the in-plan mirror symmetry −x → x of the grating.This is achieved by employing double-period the design [70]: each unitcell now consists of two sub-cell of period (1 + )a and (1 − )a with  being the symmetry-breaking coefficient [see Figure 10(a)].Consequently, the phase-shift  can be continuously tuned by changing .
To illustrate this effect, we keep a = 250 nm, FF = 0.45, and monitoring the modification of the photonic modes when scanning  from 0 to 0.25.As expected, a non-zero  turn a BIC into quasi-BIC [see Figure 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  and  as the function of  are extracted from the real-part [see Figure 10(c)] and imaginary-part [see Figure 10(d)] of the photonic modes.These results, presented in Figure 10(e) and (f), show that while  only undergoes small variation,  is continuously tuned between 0 and π∕2.Interestingly, for |L, −⟩ is still almost lossless, thus being quasi-BIC, until  exceeds 0.15 [see Figure 10

Numerical mean field methods
Equations ( 15) and ( 16) form a coupled nonlinear systems of equation to be solved for n R and Ψ. 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 Δx is chosen small enough to encompass the necessary features in momentum space.In our case, Δx = 4 μm was sufficient.A variable timestep is set by the MATLAB ® ode45 solver to reach the desired accuracy.

Figure 1 :
Figure 1: Scheme of the system and guided photonic modes.(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 | | e ±iKx ⟩ coupled with each other and coupled to the Fabry-Pérot lossy modes through the grating.(c and 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:  = 0, ℏv = 56.6 meV μm, ℏ = 1.75 meV, and ℏ = 0.18 meV corresponding to physical values obtained from rigorous coupled-wave analysis (see Methods).

Figure 2 :
Figure 2: Dirac exciton-polariton dispersion relation.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  L,+ and antisymmetric  L,− branches marked.(c and 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:  = 0, Ω∕ = 1.8, ℏ (0) X ∕ = 5.5, and  nr =  corresponding to physical values obtained from RCWA (see Methods).

Figure 3 :
Figure 3: Mode dropping dynamics in simulated Dirac polariton BIC condensates.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 th .The overlaid curves in panel (c) are taken at different time steps to show the nonstationary oscillations.Same parameters as in Figure2are used here with a pump spot size of FWHM = 20 μm (full width at half maximum).Note that we have shifted the zero energy point into the center of the bandgap.The labels |0⟩ and |1⟩ 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).

Figure 4 :
Figure 4: Condensation phase map as a function of BIC location and power.Average condensate population difference ⟨Δ  ⟩ between the symmetric and antisymmetric projections (branches) as a function of pump power, and the dissipative coupling parameter , and two different sizes of the pump spot.The colorscale is saturated around ±1 to more clearly show the white-dark regions.The red contours show total ( = Ψ † Ψ) isodensity curves.Other parameters are the same as in Figure 3.

Figure 5 :
Figure 5: A stable superposition of trapped and ballistic polaritons with a single pump spot.(a and b) Density and (c and 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 and f) Corresponding condensate densities in Fourier space showing that each component belongs to different branches.Other parameters are the same as in Figure 3.

Figure 6 :
Figure 6: Comparison of the nearfield and farfield polariton photoluminescence intensities.Curves in all panels are individually normalized.(a and b) show the trap ground state and first excited state condensate emission corresponding to Figure 3(a and d).(c and d) show the corresponding momentum space emission belonging to Figure 3(e and h).
which are simply the second Fourier coefficients of u(x).Being guided modes located below the light-line, both modes | | e i(k±K)x ⟩ are technically lossless.They become, however, lossy through coupling to lossy modes at low momentum, | | e ikx ⟩ .Notice that these lossy modes at | | e ikx ⟩ are distributed on the Fabry-Pérot modes of the stack; see Figure 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 | | e i(k±K)x ⟩ are | |  K (z)e i(k±K)x ⟩ .We ignore the dependence of the confinement wavefunction  on k.Further, we have used  +K (z) =  −K (z) because the non-perturbed structure is symmetric under x-reflection.The lossy modes at | | e ikx ⟩ are modeled by | |  (0) (z)e ikx ⟩ , | |  (1) (z)e ikx ⟩ , . . . .The matrix element scattering |

Figure 7 :
Figure 7: RCWA simulations versus effective Dirac-photon theory.(a) Sketch of the grating design.(b and 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 and 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 and 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 ℏ = 0.19 meV, ℏ = 1.75 meV, ℏv = 56.61meV μm, ℏΩ = 3.2 meV, ℏ (0) X = 2.41 meV, ℏ nr = 0.18 meV and  = 0.

Figure 8 :
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 :
Figure 9: Continuously tuning  and implementing π-jump to  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 and d) Real and imaginary part for the energy of the photonic modes when scanning FF from 0.4 to 0.86.(e and f)  and FF when scanning FF from 0.4 to 0.86.

62 and 0 Figure 10 :
Figure 10: Continuously tuning  via the symmetry-breaking coefficient .(a) Sketch of the grating design with double-period symmetry-breaking.(b) Photonic dispersion, numerically calculated by RCWA, for  = 0 and  = 0.1.(c and d) Real and imaginary part for the energy of the photonic modes when scanning  from 0 to 0.3.(e and f)  and FF when scanning FF from 0 to 0.3.