Formation mechanism of guided resonances and bound states in the continuum in photonic crystal slabs

We develop a formalism, based on the mode expansion method, to describe the guided resonances and bound states in the continuum (BICs) in photonic crystal slabs with one-dimensional periodicity. This approach provides analytic insights to the formation mechanisms of these states: the guided resonances arise from the transverse Fabry–Pérot condition, and the divergence of the resonance lifetimes at the BICs is explained by a destructive interference of radiation from different propagating components inside the slab. We show BICs at the center and on the edge of the Brillouin zone protected by symmetry, BICs at generic wave vectors not protected by symmetry, and the annihilation of BICs at low-symmetry wave vectors.

Scientific RepoRts | 6:31908 | DOI: 10.1038/srep31908 in a periodic inhomogeneous backgrounds, we find symmetry-protected BICs both at the center and on the edge of the Brillouin zone, in addition to those inside the Brillouin zone where they are not protected by symmetry. Although we only show examples for systems with one dimensional periodicity, our mode expansion method can be extended to systems that are periodic in two dimensions.

Methods
For simplicity, we consider TM modes (H x , H y , E z ) in a PhC slab that is periodic in y and uniform in z (Fig. 1a); the generalization to TE modes (E x , E y , H z ) and to PhC slabs with two-dimensional periodicity is straightforward. We consider structures that are mirror-symmetric in the normal direction x (this symmetry is necessary for reducing the number of radiation channels 6 ) and where the slab permittivity ε 1 (y) is uniform in x (common in most fabricated structures of PhC slabs). Here we consider two cases: the permittivity of the surrounding medium ε 2 (y) is a constant (Fig. 1b) or is periodic with the same period a as the slab (Fig. 1c). Inside (|x| < 0.5h) and outside (|x| > 0.5h) the slab, the structure is uniform in x, so the fields can be expanded in the eigenmodes of ε 1 (y) and ε 2 (y) with a sinusoidal dependence along x; the expansion coefficients can be determined by continuity at the slab surface |x| = 0.5h. Specifically, with an outgoing boundary condition in x, an even-in-x TM mode with wave vector k y can be written as 35 where h is the slab thickness, the cosine inside the slab guarantees the even-in-x symmetry, and the complex exponential outside the slab guarantees the outgoing boundary condition. TM modes satisfy the wave equation where k 0 = ω/c, ω is the frequency and c is the vacuum speed of light; inserting Eq. (1) into this wave equation, we find that the eigenfunctions and propagation constants u m (y), ϑ m (y), β m and γ m inside and outside the slab satisfy for l = 1, 2 are the Hermitian operators governing the wave equation for the layers inside and outside the slab, subject to periodic boundary condition u m (y + a) = u m (y), ϑ m (y + a) = ϑ m (y). For a given frequency and k y , there will be a finite number of eigenmodes with β > 0 ) that are evanescent in x. Similar expansion methods were used previously for water waves 21 and for quantum waveguides 36 . Odd-in-x TM modes can be written similarly by replacing the cosine in Eq. (1) with sine, and TE modes can be treated by replacing Ĥ l with that for H z . C m and T m are coefficients of the eigenmode expansion, and they can be determined via the continuity of E z and ∂ E z /∂ x at x = 0.5h, which requires The standing waves inside the slab [the cos (β m x) in Eq. (1)] are superposition of waves propagating in + x and in -x directions, so one can interpret the in-slab fields as waves circulating within the slab with reflection and transmission at the two slab surfaces. The transverse phase shift for every propagating component is an integer multiple of 2π after a round trip with two reflections, which is the same resonance condition as the Fabry-Pérot resonances in uniform dielectric slabs; the difference is that here multiple propagating components are coupled due to the periodicity in y.
While the two sets of eigenmodes {ϑ m } and {u m } each form an infinite-dimensional basis, the high-order ones correspond to fast oscillating fields that are negligible at low frequencies. Therefore, in our calculations, we truncate down to M terms by expanding the eigenmodes in an M-dimensional basis (details below). In the truncated basis, the transformation between the two bases {ϑ m } and {u m } is given by to an . In this way, the continuity condition Eqs (4-5) can be written in matrix form as where Note that the matrix B is purely real since β m 2 is real for all m. Eq. (4) is a linear equation group for vectors T and C. Substituting Eq. (6) into Eq. (7) yields (iγP + PB)C = 0. Non-trivial solutions exist when y where ||*|| denotes the determinant of the matrix. Therefore, solving for Eq. (8) for a given k y yields the dispersion relation ω (k y ) for arbitrary resonances and BICs, as well as regular bound states. The vector C corresponding to the zero determinant and the associated vector T from Eq. (6) yield the field profile as given in Eq. (1). At frequencies in the continuum spectrum of the surrounding medium (ω ε > k c/ y B for a homogeneous medium), some of the γ m 's are real, and f (k y , ω ) is generally complex-valued; in such region finding the zeroes of f (k y , ω ) requires searching for solutions on the lower half of the complex-frequency plane with ω r = Re(ω) and ω i = Im(ω) being the parameters. The imaginary part of the frequency is the decay rate, and the quality factor of the resonance is The transformation matrix P is deduced from Eqs (2-3). Inspired by ref. 24, we expand both sides of Eqs (2-3) in Fourier series and truncate to Fourier orders from -N to N (with a total of M = 2N + 1 terms). Then, Eqs (2-3) can be written as matrix equations The m-th column of Φ (or Θ) contains the -N-to-N Fourier coefficients of u m (or ϑ m ). Φ and Θ are transform matrices connecting {u m } and {ϑ m } to the same basis with plane-wave elements, hence P = Θ −1 Φ. Note that when the permittivity is real and mirror symmetric in y, , the matrices H l are real symmetric, so β 2 and γ 2 are real; moreover, Φ, Θ, and P can be chosen to be purely real.
BICs arise when the decay rate γ = − 2ω i of a resonance becomes zero, or equivalently when the amplitudes of the radiating waves vanish: . From Eq. (6), T m′ is given by where ′ P m m , is the (m′ ,m)'s element of matrix P, which is the m′ -th component of the in-slab eigenmode u m in the basis of {ϑ m′ }. Equation (11) reveals that the radiating wave of port m′ comes from interference of all the contributions from {u m } to the radiating mode ϑ m′ . Therefore, T m′ = 0 is a result of destructive interference from the in-slab eigenmodes to the radiation modes.
Numerically it can be ambiguous to determine whether ω i is very small or identically zero. Therefore, we use a slightly modified scheme to look for exact BICs. Define f ′ as : the propagation constants of the radiating waves γ m′ are artificially set to zero; this function f ′ is purely real for a lossless dielectric structure that is symmetric in y (where H l is real symmetric). A BIC not only satisfies f = 0; it also satisfies f ′ = 0 since T m′ = 0 for a BIC, and from Eq. (5) it can be seen that setting γ m′ to zero does not change the solution. Therefore, to search for BICs, we first solve the real-valued equation f ′ (k y , ω ) = 0 at each k y for a real-valued frequency ω ; the solution also provides a mode profile given by T and C. However, such a mode profile will only satisfy the continuity condition, Eqs (4 and 5), if T m′ = 0. In this work, we study the frequency range where there is only one leaky channel (only one m′ with γ ′ > 0 m 2 ), and we perform a root finding with k y being the free parameter to look for solutions of f ′ = 0 where the amplitude of this radiation channel vanishes, T m′ = 0. Once found, such a solution will be a true bound state at a purely real frequency and with no radiation.

Results
Photonic crystal slab in a homogeneous background. In this work we study two systems. The first one is a layered slab in a homogeneous medium (Fig. 1b) and Θ being an identity matrix. In Fig. 2a, the region with one leaky channel (one real γ m′ ) is shaded in yellow. In the slab, the Fourier coefficients of the permittivity is ξ 1 (n) = (d/a)(ε A − ε B ) sin c (nd/a) + ε B δ n . For the basis truncation, we take M = 21 Fourier terms and eigenmodes (N = 10), which is enough for the results to converge within the frequency range we consider. As an example, we take ε A = 4.9 (this is a reasonably small dielectric constant and is close to many common optical materials such as silicon nitride, zinc oxide, gallium nitride, indium tin oxide, and diamond), ε B = 1, d = 0.5a, and h = 1.4a. By solving Eq. (8), we obtain the band structure Re(ω) and the quality factor  Fig. 2a,b. As a validation, we also perform finite-difference time-domain (FDTD) simulations for the same structure. The FDTD simulations are carried out at a spatial resolution of 32 2 points per area of a 2 (high enough for the calculated Q to converge) and take significantly longer than our method; the results, shown as circles in Fig. 2a,b, are in perfect agreement with our method.
The quality factor diverges at the BICs. Non-symmetry-protected BICs occur at k y = 00.3156 (2π/a) for even-in-x modes (blue curve) and at k y = 0.1640 (2π/a) for odd-in-x modes (green curve). Symmetry protected BICs can be found at the Γ point (k y = 0), where radiation vanishes because E z is odd in y for the resonance but even in y for the radiating wave. Such symmetry-incompatibility mechanism also holds at the edge of the Brillouin zone (k y = π/a), but at the zone edge of this system modes are either regular bound states (below the yellow shaded area) or states with multiple leaky channels (above the yellow shaded area) for which BICs are harder to come by. We will show zone-edge BICs in the second system. Fig. 2c shows the field profiles of the four BICs.
In Fig. 2d, we plot the coefficients C and T of the even-in-x non-symmetry-protected BIC. Inside the slab (upper panel), the amplitudes C are dominated by two propagating modes (shown in red). Outside the slab (lower panel), there is only one propagating mode, and its amplitude vanishes at the BIC. The amplitude of the propagating mode outside the slab is the transmission from the in-slab modes, as shown in Eq. (11). Therefore, the disappearance of radiation arises from destructive interference of the transmission from the in-slab modes, which, as shown in the upper panel, primarily consist of two propagating modes.
Photonic crystal slab in a periodically-modulated background. The second system we consider is a PhC slab surrounded by a periodically-modulated background (Fig. 1c). The out-of-slab region has the same period as the slab but with a different filling fraction. We consider ε A = 4.9, ε B = 1, d 1 = 0.5a, d 2 = 0.2a, h = 1.4a. Figure 3a shows the band structure obtained from Eq. (8), with yellow shading over the region with one leaky channel (γ > 0 1 2 , γ < 0 2 2 ); the corresponding quality factor is shown in Fig. 3b. Results from FDTD simulations are shown as circles in Fig. 3a,b; again the simulation results quantitatively agree with our method.
In this structure, symmetry-protected BICs can be found both at the Γ point (k y = 0) and at the zone edge (k y = π/a) inside the yellow-shaded region. The zone-edge BICs are possible because the periodic modulation in the background breaks the degeneracy of the propagating waves on the zone edge and opens up a finite yellow-shaded region where only one leaky channel is present. On the zone edge and within this region, the leaky wave is even under mirror flip around y = 0, but the BICs are odd (as can be seen from the mode profile in Fig. 3c) so they decouple from radiation. Meanwhile, non-symmetry-protected BICs still exist, as marked by red dashed lines in Fig. 3b and with the mode profiles shown in Fig. 3c. Scientific RepoRts | 6:31908 | DOI: 10.1038/srep31908 In this system, we can observe an interesting phenomenon that two BICs annihilate at low-symmetry k points. On the odd-in-x band, as we vary the filling fraction in the background from d 2 /a = 0.23 to d 2 /a = 0.25, we observe that another non-symmetry-protected BIC emerges from the Γ point, moves along the k y axis, and then annihilates with the other non-symmetry-protected BIC near k y = 0.14 (2π/a). The annihilation removes both BICs and leaves behind a finite peak in Q, as shown in Fig. 4a. The annihilation can be understood from the radiation coefficient T 0 , which we plot in Fig. 4b. Each zero-crossing of T 0 corresponds to a BIC. Since T 0 is expected to change continuously, two adjacent zero-crossings must have opposite slope and will cancel each other when they meet. This is an example of the topological charge of BICs 18 ; here the two neighboring BICs have opposite charges and can annihilate with each other. Note that near the annihilation of the two BICs, quantitative prediction of the quality factor using FDTD becomes exceedingly hard due to an increased sensitivity on structural variations (which requires an unusually high spatial resolution in FDTD); nonetheless, we can still calculate Q efficiently with our method.

Discussion
In principle, more BICs can lie above the yellow shaded area in Figs 2a and 3a, where there are two or more radiation ports. In this case, three or more independent equations have to be satisfied simultaneously: f ′ = 0, T 0 = 0, T 1 = 0 … It will require root-finding with more variables than k y and ω , which means the structure itself has to be fine-tuned to find BICs; such structure-sensitive BICs are less practical as experimental realization will be harder.
Our analysis of PhC slabs with 1D periodicity can be considered an extension of the previous topological vortex work 18 to the 1D parameter space. As shown in Fig. 4b, here BICs correspond to nodal points where the radiation amplitude crosses zero (instead of vortex centers), which are manifestations of topological charges in 1D.
The preceding examples concern structures where the dielectric is real and symmetric in y, for which the matrices H l are real symmetric and so the function f ′ is real-valued. However, BICs can exist in even more general systems. As long as ε ε − = ⁎ y y ( ) ( ) l l , H l is real (although not necessarily symmetric and not necessarily Hermitian). In such PT-symmetric systems, if the non-Hermiticity is below the PT-breaking threshold, the eigenvalues and eigenvectors of H l can still be real, and the function f ′ can still be real-valued. Such systems can also support BICs. However, if the introduction of gain leads to lasing, one will need to account for the nonlinearity resulting from gain saturation [37][38][39] , which is beyond the linear model considered in this work.

Conclusion
We have presented a mode expansion method that can efficiently and quantitatively describe guided resonances and BICs in PhC slabs, and the method also reveals their underlying formation mechanisms. We find symmetry-protected BICs at the Γ point and at the zone edge, as well as BICs not protected by symmetry. The formalism is easily extendable and applicable to a wide range of structures. This is an attractive approach for the study of guided resonances and BICs in periodic structures.