Dirac exciton–polariton condensates in photonic crystal gratings

2.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 one-dimensional (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], [35], [36], [37], [38], [46] or perovskites [39], [40], [41], [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 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 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 ± i K x together, but also couples them with lossy modes residing at normal incidence q = 0. This renders the two guided modes e ± i K x 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 p K x . 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^−iωt 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 ₁ and U ₂ 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

If the grating design possesses mirror symmetry, x → −x, we have u(x) = u(−x) and all Fourier coefficients of u(x) are real. Hence, φ can only take values that are integer multiple of π. Notably, a π-jump of φ can be obtained by sweeping the filling fraction through a band-inversion point; see Methods. From (6), in the case of φ = 0, a BIC mode of infinite lifetime appears in the center of the lower ω ₋ branch while a lossy 2γ mode appears in the upper ω ₊ branch [see Figure 1(d)].

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 symmetry-protected 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], [50], [51], [52], [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 mirror-symmetry 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 < φ < π.

2.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 Ω (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 = ω X ( 0 ) − i γ nr where ω X ( 0 ) and γ 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. 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.

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, ℏ ω X ( 0 ) / κ = 5.5 , and γ _nr = γ corresponding to physical values obtained from RCWA (see Methods).

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 i K x and backward-propagating e − i K x photons to excitons, leading to lower-forward L , e iKx and lower-backward L , e − i K x 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,

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 v ̃ 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 − i K x polaritons, we can describe the system with the following massive non-Hermitian Dirac operator,

with new symmetric and antisymmetric lower polariton eigenstates | L , ± 〉 = A ± L , e iKx + B ± L , e − i K x with eigenenergies ω _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 κ ≪ Ω ≪ ω X ( 0 ) then Eq. (10) remains accurate and implies that waveguided TE polaritons behave approximately as Dirac particles with renormalized velocity v ̃ 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, ℏ ω X ( 0 ) / κ = 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.

3.1 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 ₀. 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].

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 ₀ = {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 Figure 2 are 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).

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], [59], [60] where stronger pumping results in a condensate dropping into lower order trap modes until it reaches the ground state.

3.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 φ = 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 ₀ 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 (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.

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.

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(b)] 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 ̂ Ψ d x / ∫ Ψ † Ψ d x , 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²(πΔ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 and the trapped condensate levels [see Figure 5(e and 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.

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.

6.1 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 i q x . 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 e i q x with that of e i ( q + n K ) x where K = 2π/a is the primitive reciprocal lattice vector and n is an integer. The relevant matrix elements are written U n = e i ( q + n K ) x u ( x ) e iqx , 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 ±(K + k) with k ≪ K. The two modes e i ( k + K ) x and e i ( k − K ) x are coupled by the potential u(x) with matrix elements U 2 = e i ( k + K ) x u ( x ) e + i ( k − K ) x = 1 u ( x ) e − i 2 K x and U 2 * = 1 u ( x ) e + i 2 K x , 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 χ K ( z ) e i ( k ± K ) x into the pth lossy mode χ ( p ) ( z ) e ikx is given by ⟨ χ ( p ) ( z ) | χ K ( z ) ⟩ 1 u ( x ) e i ( k ± K ) x = c p U ± 1 . Notice that the first Fourier coefficients of u(x) are related by U − 1 = U 1 * . 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 + 1 = χ K ( z ) e i ( k + K ) x , − 1 = χ K ( z ) e i ( k − K ) x and 0 = χ 0 ( z ) e ikx , a general wavefunction can be written as

The evolution of the coefficients ϕ ̃ + 1 k , ϕ ̃ − 1 k , ϕ ̃ 0 k follows the Schrödinger-like equation

where γ ₀ 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 v being the light velocity near K. By adjusting a global phase, one can also assume that ω _K = 0 and ω ₀ can be then replaced by the difference in frequency Δ = ω ₀ − ω _K .

Assuming that the decay rate of the lossy mode γ ₀ is much faster than the matrix element U ₂, one can adiabatically eliminate ϕ ̃ 0 k . This is done by solving ϕ ̃ 0 k in terms of ϕ ̃ + 1 k and ϕ ̃ − 1 k as

As the decaying of the lossy mode γ ₀ is fast in comparison to the dynamics of the confined mode, one can make the Markovian approximation ϕ ̃ + 1 k ( t − τ ) ≈ ϕ ̃ + 1 k ( t ) , ϕ ̃ − 1 k ( t ) = ϕ ̃ + 1 k ( t − τ ) and ∫ 0 t d τ e − ( γ 0 + i Δ ) τ ≈ ∫ 0 ∞ d τ e − ( γ 0 + i Δ ) τ = 1 / ( γ 0 + i Δ ) ≈ 1 / γ 0 . We have we also assumed γ ₀ ≫Δ in the last approximation. In the end, we then obtain

Inserting (26) into (24), we obtain the evolution equation for ϕ ̃ + 1 k and ϕ ̃ − 1 k in the form

with the non-Hermitian Hamiltonian in momentum space

where φ ₁ is defined by U 1 = | U 1 | e i φ 1 and γ = | c 0 | 2 | U 1 | 2 γ 0 . It is interesting to notice that the indirect loss rate γ for the two modes χ K ( z ) e ± i K x is inversely proportional to the lossy rate γ ₀ 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, γ ₀ ≫|U ₁|, will “freeze” the population of χ K ( z ) e ± i K x 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 ₂ in (28) is generally a complex number. We denote U 2 = | U 2 | e i φ 2 and eliminate the phase φ ₂ by the following unitary transformation,

One then obtains the non-Hermitian Hamiltonian for ϕ + 1 k and ϕ − 1 k as

with φ = 2φ ₁ − φ ₂. This is the effective Hamiltonian (2) introduced in the main text with κ = |U ₂|.

6.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 ϕ + 1 k and ϕ − 1 k are only significant at wavelengths k ≪ K. Equivalently, ϕ ₊₁(x) and ϕ ₋₁(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 ϕ 0 k , 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φ ₁ − φ ₂. The farfield intensity is then obtained as

which coincides with (22).

6.3 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.4Ga_0.6As barriers 20 nm thick grown on an Al_0.8Ga_0.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 S⁴ 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² 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.

6.4 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: ℏγ = 0.19 meV, ℏκ = 1.75 meV, ℏv = 56.61 meV μm, ℏΩ = 3.2 meV, ℏ ω ( 0 ) X = 2.41  meV, ℏγ _nr = 0.18 meV and φ = 0. Importantly, the vanishing of the intensity in the lower antisymmetric polariton branch L , − at k = 0 in Figure 7(c) confirms the infinite radiative lifetime of these polaritons that are inherited from the photonic BIC. Therefore the finite linewidth of L , − at k = 0 [see Figure 7(g)] is purely nonradiative and is inherited by the excitonic component. These results are in good agreement with the experimental observations in [40].

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.61 meV μm, ℏΩ = 3.2 meV, ℏ ω ( 0 ) X = 2.41  meV, ℏγ _nr = 0.18 meV and φ = 0.

6.4.1 Varying the exciton–photon energy detuning

By scanning the period a of the same design (i.e. FF, etching depth), the energy detuning ℏ ω X ( 0 ) between the exciton energy E _X and the mid gap of photonic modes E 0 = Re ℏ ω + + ℏ ω − / 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.

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.

6.4.2 Tuning 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.62 and 0.8 are shown in Figure 9(b). These results represent respectively three case: (i) photonic BIC in the lower band, corresponding to non-zero κ and φ = 0; (ii) gap closing and formation of Exceptional Points [29], corresponding to κ = 0; (iii) photonic BIC in the upper band, corresponding to non-zero κ and φ = π. The values of κ and φ as the function of FF are extracted from the real-part [see Figure 9(c)] and 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).

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.

6.4.3 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 α.

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.

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(d)], that corresponds to φ < 0.15π [see Figure 10(f)].

6.5 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.