Generalized non-Hermitian Hamiltonian for guided resonances in photonic crystal slabs

2.1 Eigenmodes of PhC slabs

PhC slabs are finite structures where light is confined vertically (z-direction) within a slab that consists of multilayers of high refractive index. The slab is periodically corrugated in the x-y plane, and the vertical confinement is achieved by the refractive index contrast between the slab material and its surrounding lower refractive index environment (e.g., air and substrate). The in-plane periodicity gives rise to photonic band structures, while the vertical confinement allows the PhC slab to support fully guided modes that are confined by total internal reflection, useful for integrated photonic on-chip applications, and leaky modes (guided resonances) that interact with the radiation field in free space and are critical for applications such as light-emitting devices, sensors, or detectors. Guided resonances are specific to PhC slabs and introduce non-Hermitian physics. In this work, we focus solely on PhC slabs and sometimes refer to them as PhCs for brevity.

The behavior of electromagnetic waves in PhC slabs is governed by Maxwell’s equations. Here, for the sake of simplicity, we assume the materials are isotropic, lossless, non-magnetic, and non-dispersive, simplifying the permeability μ = 1 and the permittivity to spatially varying constants ϵ(r). For time-harmonic fields E(r, t) = E(r)e^−iωt , the master equation for the electric field is given by:

This eigenvalue equation describes the frequency ω of electromagnetic waves as a function of the spatial variation of ϵ(r), laying the foundation for the photonic band structure in PhCs. It is important to note that this equation is not a standard eigenvalue equation but a generalized one, requiring specific formulations for the inner product and orthogonality, which involve weighting by the permittivity ϵ(r). Here the permittivity exhibits in-plane periodicity, satisfying ϵ(r) = ϵ(r + R) with R being a lattice vector in the xy plane. Since the periodicity exists in the in-plane directions (x, y) but not along z, the eigen Bloch modes are given by:

where k _∥ = (k _x , k _y ) is the in-plane wavevector, r _∥ = (x, y) are the in-plane spatial coordinates, and u k ∥ ( r ) is a periodic function with the same periodicity as ϵ(r). The eigenvalue problem yields the photonic band structure, ω(k _∥), which maps the allowed frequencies for each wavevector k _∥ within the first Brillouin zone.

The emergence of complex eigenvalues for guided resonances arises from their interaction with the radiative continuum. This coupling to the radiative continuum is a hallmark of non-Hermitian physics in PhC slabs. Mathematically, the eigenvalue problem for guided resonances can be written as H ̂ E = ω E , where the Hamiltonian H ̂ is non-Hermitian due to the inclusion of radiative losses. A non-Hermitian Hamiltonian is characterized by H ̂ ≠ H ̂ † , meaning it is not equal to its adjoint. This property introduces complex eigenvalues, where the imaginary part typically represents energy gain or loss in the system. In PhC slabs, the non-Hermitian nature of H ̂ stems from the open-system interaction between guided resonances and the radiative continuum, which allows energy leakage. This is a direct departure from Hermitian physics, where systems are typically closed and do not interact with an external environment. The non-Hermitian framework enables the study of unique phenomena such as EPs and BICs, which are absent in purely Hermitian systems.

2.2 Effective Hamiltonian from Δϵ perturbation

Before discussing the Hamiltonian of generalized guided modes expansion, we first introduce the common formalism of the coupled mode theory via Δϵ perturbation. The unperturbed system, characterized by ϵ ₀(r), has known eigenvalues ω _n and eigenmodes E _n (r): ∇ × ∇ × E n = ω n 2 c 2 ϵ 0 ( r ) E n . Introducing Δϵ(r) perturbs these eigenmodes and eigenvalues, allowing us to write: E(r) = ∑ _n c _n E _n (r), where c _n are expansion coefficients that account for the perturbation. Substituting E(r) into the perturbed master equation ∇ × ∇ × E = ω 2 c 2 ϵ 0 ( r ) + Δ ϵ ( r ) E and projecting onto the unperturbed basis modes {E _n (r)} in the approximation ω ≈ ω ₀ to achieve an eigenvalue problem of ω instead of ω ², we obtain:

where C = {c _n } is the vector of expansion coefficients. Here, the effective Hamiltonian H _eff includes the perturbative contributions from Δϵ(r):

with:

2.3 Field expansion in the generalized guided mode expansion

In the generalized guided mode expansion, the unperturbed system corresponds to the unpatterned slab with a permittivity profile ϵ ₀(z), representing the zeroth Fourier component of ϵ(r). The perturbation arises from the periodic modulation in the x-y plane, expressed as:

where G denotes the reciprocal lattice vectors, and r _∥ = (x, y) represents the in-plane coordinates.

The eigenmodes of the unpatterned slab ϵ ₀(z) include the guided modes of the planar waveguide and the radiative Fabry-Pérot (FP) modes. In this work, we focus on the low-energy regime, where the wave vectors G n + k ∥ primarily describe guided modes with evanescent out-of-plane components, while the wave vectors k _∥ correspond to modes inside the light cone – i.e., radiative FP modes. The total electric field in the PhC slab can thus be expanded as:

where E n = u n ( z ) e i ( G n + k ∥ ) ⋅ r ∥ are the guided modes confined within the slab, having corresponding eigenvalues ω _n . And E m , FP = u m , rad ( z ) e i k ∥ ⋅ r ∥ are the radiative FP modes propagating outside the slab, having corresponding eigenvalues ω ̃ m , FP = ω m , FP − i γ m , FP . We note that in this study, we restrict ourselves to the low-energy regime, where only the zeroth-order FP modes contribute to radiation. However, the formalism can be naturally extended to higher-energy regimes by including FP modes associated with higher-order diffraction channels. We also focus exclusively on transverse electric (TE) guided modes. The same formalism can be applied to transverse magnetic (TM) guided modes by formulating the perturbation theory in terms of the magnetic field H, rather than the electric field E as done here. We also neglect TE–TM coupling, which is justified in our system due to the strong spectral separation between TE and TM modes in high-contrast dielectric slabs.

It is important to note that despite the radiative Fabry–Pérot modes being quasinormal modes – characterized by complex eigenfrequencies due to their leaky nature, they can form a complete and orthogonal basis under specific conditions. The completeness and orthogonality of quasinormal modes in leaky optical cavities were rigorously established by Leung, Liu, and Young in 1994 [34]. Their work demonstrated that, despite the non-Hermitian nature of such systems, a discrete set of quasinormal modes suffices to accurately describe the field distribution within the cavity. Therefore, the guided modes (E _n ) and radiative Fabry–Pérot modes (E _m,FP) together provide a complete and orthogonal basis that ensures the expansion from Eq. (7) is both accurate and comprehensive.

The perturbation Δϵ(r) introduces coupling between these guided and radiative modes. The coupling matrix elements are:

These coupling terms correspond to three distinctive mechanisms:

1. Guided-to-Guided Coupling: The interaction between two guided modes is mediated by the Fourier component Δ ϵ G n − G m ( z ) , which matches their momentum difference. This term is crucial for photonic band structure modifications, such as band gaps.

2. Guided-to-FP Coupling: The coupling of a guided mode to a radiative (Fabry–Pérot) mode depends on the Fourier component Δ ϵ − G n ( z ) . This term introduces radiative losses to guided modes, converting them into quasi-normal modes.

3. FP-to-Guided Coupling: Radiative modes contribute to guided modes through the Δ ϵ G m ( z ) term. This term describes how energy can leak from FP modes back into guided modes. This is the reversed mechanism of the guided-to-FP coupling, evidenced by Δ H nm guided − FP = Δ H nm FP − guided * .

Note that as we restrict ourselves to the low-energy regime, where only the zeroth-order FP is involved, there is no FP-to-FP coupling.

2.4 Effective Hamiltonian for guided modes

The total Hamiltonian for the coupled system of guided and radiative modes can be written as:

which satisfies the eigenvalue problem:

where A = {a _n } and B = {b _m } are the expansion coefficients of the total electric field E in Eq. (7). From the second row of Eq. (10), we obtain a relation between the coefficients of the radiative modes and those of the guided modes B = ω − H FP − 1 H FP − guided A , where ω − H _FP is a diagonal matrix with entries Δ _n + iγ _n,rad, and Δ _n = ω − ω _n,rad denotes the detuning from the radiative FP mode frequencies.

Substituting this relation into the first row of Eq. (10) yields an effective operator acting on the guided-mode amplitudes:

where Σ(ω) is the self-energy term, given by Σ ( ω ) = H guided − FP ( ω − H FP ) − 1 H FP − guided . This effective operator satisfies

but the equation remains nonlinear due to the explicit ω-dependence of the self-energy term Σ(ω). Physically, this term encapsulates the modification of the guided-mode dynamics via coupling to the radiative Fabry–Pérot modes. The imaginary part of Σ(ω) captures the irreversible coupling to radiation and gives rise to the non-Hermitian character of the guided resonances. Its real part produces only a small dispersive frequency shift, which can be fully absorbed into the guided-mode Hamiltonian through the renormalization H _guided → H _guided + Re Σ, thereby shifting only the diagonal terms of H _guided without altering modal profiles, radiation rates, or topological properties. This procedure is standard in open electromagnetic and quantum systems – directly analogous to Lamb-shift corrections in QED, real-part self-energies in Green’s-function theory, resonance-frequency shifts in temporal coupled-mode theory, and on-site energy shifts in tight-binding models.

In the limit where the FP modes form a broadband and rapidly decaying reservoir (e.g., for thin or weakly confined slabs), the imaginary part of Σ dominates over the real part. However, even when this strict quasi-continuum condition is not satisfied, in practice, the diagonal parameters of H _guided in numerical examples are extracted from full-wave simulations, so any dispersive contribution from Re Σ is naturally absorbed into these extracted coefficients. This makes the approximation robust. In the concrete examples presented in Sec.3, this renormalization simply corresponds to a small adjustment of the guided-mode resonance frequencies ω _Γ and ω _K . Retaining only the imaginary part of the self-energy therefore yields a compact closed-form non-Hermitian effective Hamiltonian acting on the guided-mode subspace:

This is the non-Hermitian effective Hamiltonian we set out to derive. It captures both the dispersive properties of the guided modes and their radiative losses through coupling to the continuum. Importantly, the approximation eliminates the ω-dependence in the self-energy, so that H ̂ now defines a standard linear eigenvalue problem.

Moreover, within the quasi-continuum regime of FP modes, the radiation amplitudes B – corresponding to the far-field leakage – can be directly computed from the eigenvector A of H ̂ using:

2.5 Compact expression and physical meaning of the non-Hermitian Hamiltonian

The expressions of the coupling terms of H ̂ , given in Eq. (8), are simplified by summing out the polarization cross products (see Appendix A). We now discuss the compact expression and physical meaning of each of these terms.

1. The first term H _guided represents the unperturbed guided modes and the diffractive couplings between them. The diagonal elements correspond to the frequencies of guided modes in the unpatterned slab:

   (15) H nn guided = ω n .

   The off-diagonal elements describe direct coupling between guided modes via diffractive mechanism:

   (16) H nm guided = p n ⋅ p m U nm ,

   where p _n , p _m are the polarization vector of the guided mode E _n , E _m , respectively; and the coupling strength U _nm is given by:

   (17) U nm = ω 0 2 ∫ u n * ⋅ ϵ G n − G m ⋅ u m d z .

2. The second term H ̂ rad represents the radiative losses and radiative coupling of the guided modes. The diagonal elements introduce imaginary components in the eigenfrequencies of the guided modes due to coupling with radiative modes:

   (18) H nn rad = − i ∑ l γ n ( l ) .

   Here γ n ( l ) represents the radiative losses of the guided mode n via the radiative Fabry–Pérot mode l, and the expression of γ n ( l ) is given by:

   (19) γ n ( l ) = ω 0 2 ∫ u n * ⋅ ϵ − G n ⋅ u l , rad d z 2 4 c 4 γ l , rad .

   The off-diagonal elements describe indirect coupling between guided modes mediated by radiative modes:

   (20) H nm rad = − i p n ⋅ p m ∑ l γ n ( l ) γ m ( l ) e i ( ϕ n ( l ) − ϕ m ( l ) ) ,

   where coupling phase ϕ n ( l ) is given by:

   (21) ϕ n ( l ) = arg ∫ u n * ⋅ ϵ − G n ⋅ u l , rad d z .

It is worth emphasizing that the coupling terms expressed in Eqs.(17) and (20), derived here from first-principles Maxwell equations via the guided-mode expansion, correspond directly to the so-called diffractive coupling and radiative coupling often introduced phenomenologically in the literature on 1D photonic gratings [21], [35] and 2D PhC slabs [11], [36]. In particular, Eq.(20), which captures the off-diagonal elements of the radiative loss operator, is in excellent agreement with the inter-mode coupling terms found in the radiative Hamiltonian of temporal coupled-mode theory applied to non-orthogonal resonators [20]. This highlights the consistency between the full-wave modal expansion approach and reduced-order models, and provides a rigorous microscopic foundation for the coupling coefficients often assumed in heuristic or fitted models.

Furthermore, we highlight that the radiative losses of the guided mode n via the radiative channel l is governed by the overlap integral ∫ u n * ⋅ ϵ − G n ⋅ u l , rad d z , as shown in Eq. (19). It can be accidentally suppressed if the overlap integral is zero, leading to the suppression of a radiative channel for a given guided mode. In particular, for some particular design, it is possible to render null the radiative losses of every guided modes in the spectral window of interest. In such a configuration, the eigenmode is evidently lossless. This corresponds to an accidental BIC configuration, well documented in the literature, whose emergence is extremely sensitive to the slab thickness, effective index, and requires a vertical symmetry design [37], [38], [39].

2.6 From eigenmodes to nearfield pattern and farfield pattern of Bloch modes

The eigenmodes of the system are obtained by diagonalizing the effective non-Hermitian Hamiltonian H ̂ through the characteristic equation det(H − ω) = 0. Due to the non-Hermitian nature of H ̂ , its eigenvalues Ω _n are generally complex, and the corresponding eigenvectors Ω n describe the resonant modes of the structure. If the imaginary part of an eigenvalue vanishes, i.e. Im(Ω _n ) = 0, the corresponding eigenmode Ω n is lossless. Such states are known as BICs. Conversely, if Im(Ω _n ) < 0, the mode exhibits radiative loss and is referred to as a leaky mode. The eigenvector A n = ( a n 1 , a n 2 , a n 3 , . . ) associated with Ω n defines the near-field spatial distribution of the electric field:

A similar superposition principle applies to the magnetic field. Due to the TE nature of the guided modes, the magnetic field is predominantly polarized along the z-direction. The corresponding near-field profile of the magnetic field can thus be expressed as a scalar field:

The far-field radiation pattern is governed by the FP components in the field expansion of Eq. (7). For a given eigenmode Ω n , the FP coefficients B n = ( b n 1 , b n 2 , b n 3 , … ) can be obtained from the eigenvector A _n using Eq. (14). The resulting far-field electric field is given by:

where α m = ∑ l γ m ( l ) , e i ϕ m ( l ) accounts for the radiation amplitude and phase into each far-field channel l. We note a direct correspondence between the near-field and far-field expressions in Eqs. (22) and (24): the far-field radiation is derived from the near-field by replacing the mode profiles u _m (z), e^{iG m⋅r∥} with radiating plane wave components e i k z z , p m . This procedure effectively represents the folding of guided Bloch modes into the first Brillouin zone and their subsequent coupling to the radiation continuum.

The polarization texture – including polarization orientation, ellipticity, and topological charge of polarization singularities – can be readily computed from the far-field electric field (see Appendix E).

3.1 Hexagonal lattices with C ₆ symmetry

3.1.1 System description

We now apply our general non-Hermitian Hamiltonian to described guided resonances in the vicinity of the Γ and K point of a PhC slab with hexagonal lattices. The case of grating and square lattices is presented in the Appendix F. The PhC slab of consideration is of a hexagonal lattice with C ₆ symmetry, for example, a triangular lattice with a single circular hole in the unit cell (see Figure 2a), or a honeycomb lattice with two identical circular holes (see Figure 2b). The unit vectors are given by a 1 = Λ 3 2 , 1 2 and a 2 = Λ 3 2 , − 1 2 with Λ = a for the triangular lattice, and Λ = 3 a for the honeycomb lattice.

Figure 2:

Geometry of hexagonal lattices. a) Triangular lattice. b) Honeycomb lattice. c) First Brillouin zone in momentum space.

In the momentum space, the high symmetry points are K, K′, and M at the edge of the Brillouin zone, and Γ at the centre of the Brillouin zone (see Figure 2c). The corresponding unit vectors in momentum space are given by: b 1 = b 1 2 , 3 2 , and b 2 = b 1 2 , − 3 2 , with b = 4 π 3 a for the triangular lattice, and b = 4 π 3 a for the honeycomb lattice. We note that |b ₁ | = |b ₂ | = b is the distance between the Γ point of the first Brillouin zone (BZ) to the six closest Γ points of the neighbor (i.e. second) BZs (see Figure 3a).

Figure 3:

Guided mode basis. a) Eigenmodes operating at the Γ point and above the light cone are described by six guided modes Γ n of wave vector Γ _n , with n = 1 → 6. b) Eigenmodes operating at the K point and above the light cone are described by three guided modes K n of wave vector K _n , with n = 1 → 3. The green region indicates the first Brillouin zone, the blue regions indicate the second Brillouin zone, and the yellow region indicates the third Brillouin zone.

3.2 Guided mode basis

In this work, we only focus on photonic modes above the line cone. In this region, the lowest photonic bands at the Γ point is described by the basis consisting of guided modes operating at the six Γ points of the second BZs(see Figure 3a). These correspond to six guided modes Γ n , with n = 1 → 6, that both originate from the fundamental TE-guided mode of the slab ϵ ₀(z) and have wavevector Γ _n (see Figure 3a).

On the other hand, the lowest photonic bands at the K point is described by the basis consisting of guided modes operating in three K points, two of them are in the second BZ and the third one is in the third BZ (see Figure 3b). These correspond to three guided modes K n with n = 1 → 3, that both originate from the fundamental TE-guided mode of the slab ϵ ₀(z) and have wavevector K _n (see Figure 3b). In a general way, Γ _n = n ₁ b ₁ + n ₂ b ₂ and K _n = n ₁ b ₁ + n ₂ b ₂ + K with K = b 1 2 , 3 4 is the position of the K point in the first BZ. The couples (n ₁, n ₂) are given by:

for Γ 1 → Γ 6 and K 1 → K 3 , respectively. We note that Γ n = b and K n = 2 3 b . Explicit expressions of Γ _n and K _n are provided in the Appendix B.

The six guided modes Γ n have the same energy ω _Γ at the Γ point, and the three guided modes K n have the same energy ω _K at the K point. Adding an in-plane momentum k = (k cos φ, k sin φ) with |k|≪ b, the propagation vector of Γ n and K n becomes Γ _n + k and K _n + k, respectively. Thus the energy of Γ n and K n at k are given by ω Γ n ( k ) ≈ ω Γ + v Γ Γ n ⋅ k Γ n and ω K n ( k ) ≈ ω K + v K K n ⋅ k K n where v _Γ and v _K are the group velocities of the guided mode at Γ and K, respectively. Explicit expressions of ω Γ n ( k ) and ω K n ( k ) are provided in the Appendix C. Due to the TE nature, our guided modes exhibit in-plane polarization that are perpendicular to their wavevector. Thus p Γ n ⋅ Γ n + k = p K n ⋅ K n + k = 0 . In the limit of k ≪ b, the previous condition is approximately reduced to p Γ n ⋅ Γ n = p K n ⋅ K n = 0 . Thus one may easily obtain p Γ n and p K n (see Appendix D).

Finally, since Γ n belongs to the same planar waveguide mode, they exhibit the same vertical confinement profile u _Γ(z). The same, all three K n exhibit the same vertical confinement profile u _K (z). We note that each guided mode is associated with a pair (n ₁, n ₂) that corresponds to a Bloch vector G _n = n ₁ b ₁ + n ₂ b ₂ . Therefore, the electric field of our two basis are given by:

Here the momentum k _∥ within the first BZ is given by k for Γ n and K + k for K n .

3.3 Effective Hamiltonian

We now apply the general expressions derived in Section 2.5 to the two mode bases Γ _n and K _n . The diagonal elements of the guided-mode Hamiltonian H _guided are simply the dispersion relations of the corresponding modes:

The off-diagonal elements, corresponding to diffractive coupling between guided modes, take the form:

where the coupling coefficients are defined as:

Here, ϵ n 1 − m 1 , n 2 − m 2 denotes the Fourier component ϵ G n − G m of the in-plane permittivity modulation Δϵ(r _∥). Thanks to the C ₆ symmetry of the lattice and using the specific mode indices (n, m) given in Eq. (25), along with the polarization vectors p G n and p K n defined in Eq. (47), the coupling matrices acquire the simplified form:

The couplings U, V, W, and T are all real-valued and governed by specific Fourier components of Δϵ:

The C ₆ symmetry also simplifies the structure of the radiative Hamiltonian H ̂ rad . The diagonal elements are given by:

where γ ₀, γ ₁, and γ ₂ represent the radiative loss rates for the modes at Γ and K, respectively, under the assumption of a single radiation channel. These are governed by:

The off-diagonal elements are expressed as:

Finally, the full non-Hermitian effective Hamiltonians for guided resonances near the Γ and K points are constructed as:

and

3.4 Eigenmodes at Γ: emergence of symmetry-protected BICs

At the Γ point, H ̂ Γ ( k ) can be diagonalized analytically, yielding six eigenmodes Ω n = 1 , … , 6 ( Γ ) with eigenvalues:

From these expressions, we identify four BICs at the Γ point: Ω 1 ( Γ ) and Ω 4 ( Γ ) are non-degenerate, while Ω 2 ( Γ ) and Ω 3 ( Γ ) form a doubly degenerate pair. The remaining two modes, Ω 5 ( Γ ) and Ω 6 ( Γ ) , are leaky modes and also form a degenerate pair.

The corresponding eigenvectors, non-normalized, at k = 0 are:

Using Eqs. (22) and (23), the near-field distributions E Ω n ( Γ ) and H Ω n ( Γ ) can be computed analytically. Remarkably, these spatial field patterns are fully determined by the symmetry of the eigenvectors and are independent of the specific values of the coupling parameters U, V, and W. Figure 4 presents the calculated magnetic near-field profiles. Based on the spatial symmetry of these modes, we assign:

These modal patterns can also be classified according to the irreducible representations of the C ₆ point group symmetry [3], offering a clear group-theoretical interpretation of their polarization textures. We emphasize that the general near-field patterns predicted by our effective Hamiltonian model are in perfect agreement with those reported in the literature using full-wave finite-difference time-domain (FDTD) simulations [2], for the six photonic modes at the Γ point of a triangular lattice. This agreement not only validates the accuracy of our model but also underscores its ability to capture the essential physics of PhC slabs with high symmetry.

Figure 4:

General near-field patterns. Calculated magnetic near-field profiles H Ω n ( Γ ) for the eigenmodes at the Γ point. Arrows represent the electric field vectors E Ω n Γ .

Moreover, using Eq. (24), the radiation pattern of these six photonic modes can be computed. One may confirm that the farfield radiation of the monopolar mode corresponds to a polarization singularity of topological charge 1, while the topological charge of the hexapolar mode is −2, and the degenerated quadipolar modes are pinned at a polarization singularity of topological charge −2.

Interestingly, an accidental degeneracy between the quadripolar modes and the hexapolar mode occurs when the condition W = 2 U + V 3 is satisfied. Under this configuration, the resulting triple degeneracy at the Γ point consists of one quadratic band and a pair of Dirac cones. This accidental degeneracy, involving three BICs, is particularly relevant for applications such as zero-refractive-index metamaterials [40], [41] and scalable lasing [42], [43], both of which benefit from lossless Dirac cones at the Γ point.

3.5 Eigenmodes at K: symmetry-protected Dirac dispersion

In the vicinity of the K point, the non-Hermitian Hamiltonian H ̂ K ( k ) can be approximately diagonalized in the regime T ≫ γ ₁, γ ₂. This yields three eigenmodes Ω n = 1,2,3 ( K ) with eigenvalues:

The corresponding eigenvectors at k = 0 are:

The first two modes, Ω 1 ( K ) and Ω 2 ( K ) , are degenerate at k = 0 and split linearly with k, forming a Dirac cone centered at the K point. This Dirac cone is robust against variations in γ ₁, γ ₂, and T, as long as the condition T ≫ γ ₁, γ ₂ holds.

3.6 Effective theory versus numerical simulations near the Γ point

In this section, we focus on the eigenmodes in the vicinity of the Γ point, where various symmetry-protected and accidental BICs and EPs can emerge. The analysis of modes near the K point – characterized by distinct degeneracy lifting and topological transitions under C ₆ symmetry breaking – will be presented in the next section.

The simulated structures consist of air hole arrays in dielectric slabs with two lattice geometries: triangular and honeycomb. In both designs, the lattice constant is a = 400 nm, the air hole diameter is D = 0.35a, and the slab thickness is h = 100 nm. The refractive index of the slab is n = 2.0, and the structures are embedded in air.

3.6.1 Complex band structures and symmetry-protected BICs

We first show that the photonic band structures near the Γ point are accurately described by the effective non-Hermitian Hamiltonians introduced in Eq. (33). To validate this theory, we perform full-wave simulations using the finite element method (FEM) implemented in COMSOL Multiphysics. Floquet boundary conditions are applied in the in-plane directions, while perfectly matched layers (PMLs) along z model radiation into the far field. Complex eigenfrequencies Ω _n (k) are computed for a dense sampling of k-points near the Γ point, and are fitted using the analytical eigenvalues of the effective Hamiltonians. The corresponding fitting parameters are listed in Appendix G.

Figure 5a and d presents the real parts of the eigenfrequencies for the triangular and honeycomb lattices, respectively. The band structures show excellent agreement between numerical simulations and the analytical calculations, both for the real and imaginary parts of the eigenfrequencies, as evidenced by the quality factors plotted in Figure 5b and e. In particular, the expected scaling laws for BICs are recovered: Q ∝ 1/k ^2q , where q is the topological charge. The hexapolar mode with q = −2 follows Q ∝ 1/k ⁴, while the monopolar (q = +1) and each of the two quadrupolar modes (q = −1) exhibit Q ∝ 1/k ², fully consistent with theoretical calculations.

Figure 5:

Eigenmodes near Γ for triangular and honeycomb lattice design. a, d) Real part of the photonic band energies as a function of the in-plane wavevector. The zoomed-in insets highlight the crossing along ΓK(blue box) and anticrossing ΓM (green box) between the third and fourth bands of the triangular lattice. b, e) Quality factors of the photonic bands. Green and blue dashed lines indicate reference curves proportional to 1/k ⁴ and 1/k ², respectively. Red scatters represent numerically simulated photonic bands, while black lines show their corresponding analytical fitting using the effective theory. c, f) Far-field polarization textures (i.e., the orientation of radiated polarization) of the photonic bands hosting monopolar, quadrupolar, and hexapolar modes.

3.6.3 Emergence of chiral exceptional points

Beyond BICs, the effective theory also successfully predicts the emergence of EPs in the PhC slabs. As pointed out in Ref. [44], EPs are expected to appear near crossings of bands with opposite symmetry. To identify possible EPs, we examine the band structures along high-symmetry directions. In the triangular lattice, the third and fourth bands cross along the Γ → K direction (Figure 5a), suggesting the emergence of EPs in their vicinity. Due to the C ₆ symmetry, there are six equivalent crossing points. Without loss of generality, we focus on k _c = (0, k _c ). To probe the EPs, we map the amplitude and argument of the complex gap between the third and fourth bands around k _c . As shown in Figure 6a and b, two EPs are clearly identified by the vanishing of the gap amplitude and the presence of a singularity in the phase. The winding number of each EP, w = ± 1 2 , is computed from the gap argument as w = 1 2 π ∮ C d k ⋅ ∇ k arg ω 4 ( k ) − ω 3 ( k ) .

Figure 6:

Chiral EPs in triangular lattice design. a, b) Amplitude (a) and argument (b) of the complex gap between the third and fourth band in the vicinity of the crossing point of Figure 5a. Left panels are results obtained from numerical simulations, while right panels are from the analytical models. c) Ellipticity of the far-field polarization of the third and fourth bands in the vicinity of the crossing point k _c = (0, k _c ) of Figure 5a.

The phase map of the gap also reveals a bulk Fermi arc (BFA) connecting the two EPs [15], [44], characterized by a π jump in the argument (Figure 6b), marking a degeneracy of the real parts of the eigenfrequencies. The agreement between the effective theory and numerical simulations for both the gap amplitude and phase confirms the robustness of our model in capturing non-Hermitian degeneracies and their topological features.

Finally, we compute the ellipticity of the far-field polarization (see Appendix E for details) of the third and fourth bands near k _c . As shown in Figure 6c, the two EPs exhibit opposite handedness in their polarization textures, confirming their chiral nature. To the best of our knowledge, this is the first demonstration of chiral EPs in a triangular lattice without explicit symmetry breaking.

3.7 Hexagonal lattices with broken C ₆ symmetry

The C ₆ symmetry is broken either by using elliptical holes instead of circular ones in a triangular lattice, or by using two circular holes of different sizes in a honeycomb lattice (see Figure 7). In general, breaking the C ₆ symmetry lifts the degeneracy of the quadrupolar and dipolar modes at the Γ point, as well as the Dirac point degeneracy at the K point. However, depending on the specific geometry of the symmetry breaking, the form of the effective Hamiltonian will differ. In this section, we focus specifically on the band structure in the vicinity of the K point, where the lifting of Dirac degeneracies gives rise to rich topological and non-Hermitian phenomena [45], [46], [47].

Figure 7:

Hexagonal lattices with broken C ₆ symmetry. a) Triangular lattice with elliptical holes. b) Honeycomb lattice with two circular holes of different sizes.

3.7.2 Band structure: effective theory versus numerical simulations near the K point

To investigate the role of C ₃ symmetry, we designed a structure consisting of a triangular lattice of elliptical air holes with lattice constant a = 440 nm, slab thickness h = 180 nm, and refractive index n = 2.02, placed on a glass substrate with n = 1.46. The elliptical holes are defined by their semi-axes p and q, allowing controlled breaking of higher-order rotational symmetries. The results of three representative cases are presented in Figure 8, showing both the band structure and the associated quality factors.

Figure 8:

Triangular lattice with elliptical holes. Red scatters represent numerically simulated photonic bands, while black lines show their corresponding analytical fitting using the effective theory. The upper panel shows the band structures near the K point, while the lower panel depicts the corresponding quality factor for each band. (a) p = q = 60 nm. (b) p = 50 nm and q = 70 nm. (c) p = 40 nm and q = 80 nm.

In Figure 8a, we consider the high-symmetry case where p = q = 60 nm, which preserves C ₆ symmetry. In this configuration, a Dirac point is formed at the K point by the crossing of two upper bands, while the lowest band remains isolated. Interestingly, this lowest band exhibits a pronounced quasi-BIC character: its quality factor reaches a sharp peak (exceeding 10⁵) precisely at the momentum corresponding to the Dirac point. We then break the C ₃ symmetry while preserving inversion symmetry (C ₂) by elongating one semi-axis and reducing the other. Specifically, in Figure 8b, we increase p by 10 nm and decrease q by 10 nm. The Dirac point is no longer pinned to the high-symmetry K point but shifts along the K → Γ direction, appearing at 3(k _y − k _K )a/4π = −0.0075. Notably, the quasi-BIC peak in the lowest band follows this shift, indicating that the momentum-space location of maximal radiation suppression remains locked to the displaced Dirac crossing. This trend becomes more pronounced as the symmetry breaking increases. In Figure 8c, the semi-axes differ by 40 nm, and the Dirac point moves further to 3k _y a/4π = −0.015. Owing to the preserved C ₂ symmetry, this displacement is symmetric: if the major axis were instead aligned along the y-direction, the shift would occur in the opposite direction. These observations confirm that while the Dirac degeneracy persists due to inversion symmetry, its location in momentum space is no longer protected by C ₃ symmetry and becomes tunable through geometry.

Crucially, the momentum-dependent complex eigenfrequencies obtained from full-wave simulations are in excellent quantitative agreement with the calculations of the effective non-Hermitian Hamiltonian. This confirms that our analytical model faithfully captures both the band dispersion and the quasi-BIC behavior induced by symmetry breaking.

To investigate the role of inversion symmetry (C ₂) in honeycomb lattices, we consider a slab similar to the previous cases, but with a reduced lattice constant of a = 400 nm and air holes of different radii r ₁ and r ₂. When inversion symmetry is preserved (i.e., r ₁ = r ₂ = 50 nm), the structure exhibits a Dirac point at the K point, as shown in Figure 9a, consistent with the symmetry-protected degeneracy of the honeycomb lattice. However, when inversion symmetry is broken by introducing a small size asymmetry between the two sublattices (e.g., r ₁ = 50 nm and r ₂ = 55 nm), the Dirac point degeneracy is lifted, and a bandgap opens at the K point, as seen in Figure 9b. This gap becomes significantly larger with stronger symmetry breaking. In Figure 9c, a larger contrast between r ₁ and r ₂ results in a pronounced gap, demonstrating how geometric perturbations directly control the topological features of the band structure. Once again, the calculations of the effective non-Hermitian Hamiltonian show excellent quantitative agreement with the full-wave numerical simulations, confirming the accuracy and robustness of the theoretical model.

Figure 9:

Honeycomb lattice bands with different hole sizes. Red scatters represent numerically simulated photonic bands, while black lines show their corresponding analytical fitting using the effective theory. (a) r ₁ = r ₂ = 50 nm. (b) r ₁ = 50 nm and r ₂ = 55 nm. (c) r ₂ = 40 nm and r ₂ = 60 nm.

Compared to the case of triangular lattices with elliptical holes, where breaking C ₃ symmetry (while preserving inversion symmetry C ₂) causes the Dirac point to shift in momentum space without lifting the degeneracy, the honeycomb lattice exhibits a qualitatively different response: breaking inversion symmetry directly opens a bandgap at the K point. This contrast underscores the distinct roles of C ₂ and C ₃ symmetries in protecting Dirac points. Moreover, the energetic ordering and radiative properties of the bands also differ significantly between the two cases. In the triangular lattice, the singly degenerate band with quadratic dispersion lies below the Dirac point and exhibits a quasi-BIC character, with strongly suppressed radiation losses at the Dirac momentum. In contrast, for the honeycomb lattice, this quadratic band lies above the Dirac point, and all three bands near the K point exhibit significant radiation losses. No quasi-BIC behavior is observed in this case, reflecting the absence of symmetry protection and destructive interference mechanisms that suppress radiation. These differences further highlight how lattice geometry and symmetry breaking govern both the topological and radiative characteristics of the photonic band structure.

3.8 Parameter retrieval and computational efficiency

To determine the parameters of the analytical Hamiltonian, full-wave simulations are required only at a single high-symmetry point (e.g., Γ or K). From the complex eigenfrequencies at this point, we extract all model coefficients – including the coupling parameters U, V, W, the diagonal frequencies ω _Γ and ω _K , and the radiative loss rate γ ₀ (numerical values of parameters used in the results are reported in Appendix G. Once these coefficients are known, the effective Hamiltonian reproduces the full complex band structure in the vicinity of the high-symmetry point, including both the radiative linewidths and the far-field polarization textures, through instantaneous matrix diagonalization.

In our examples, subtle non-Hermitian features such as off-Γ bound states in the continuum and chiral exceptional points were first revealed by the analytical Hamiltonian. Only after their approximate locations were identified did we refine our full-wave simulations – using significantly increased mesh density and finer k-space sampling – to confirm these features numerically. This highlights the predictive power of the analytical model and its ability to guide full-wave solvers toward the relevant regions of parameter space.

For a representative hexagonal-lattice structure, a full-wave FEM sweep ( ∼ 50,000 mesh elements) required approximately 5 h on a standard desktop computer (AMD Ryzen 7 processor, 3.3 GHz; RAM 16 GB) to evaluate 1,000-points and 7 frequency samples. In contrast, once the Hamiltonian parameters were extracted, the analytical model generated the corresponding results within a fraction of a second on the same hardware. This demonstrates the substantial computational advantage of the proposed framework, particularly for broad parameter scans or for exploring high-Q resonances.