Bound states in the continuum based on the total internal reflection of Bloch waves

ABSTRACT A photonic-crystal slab can support bound states in the continuum (BICs) that have infinite lifetimes but are embedded into the continuous spectrum of optical modes in free space. The formation of BICs requires a total internal reflection (TIR) condition at both interfaces between the slab and the free space. Here, we show that the TIR of Bloch waves can be directly obtained based on the generalized Fresnel equations proposed. If each of these Bloch waves picks up a phase with integer multiples of 2π for traveling a round trip, light can be perfectly guided in the slab, namely forming a BIC. A BIC solver with low computational complexity and fast convergence speed is developed, which can also work efficiently at high frequencies beyond the diffraction limit where multiple radiation channels exist. Two examples of multi-channel BICs are shown and their topological nature in momentum space is also revealed. Both can be attributed to the coincidence of the topological charges of far-field radiations from different radiation channels. The concept of the generalized TIR and the TIR-based BIC solver developed offer highly effective approaches for explorations of BICs that could have many potential applications in guided-wave optics and enhanced light–matter interactions.

As a platform for nanophotonics, photoniccrystal (PhC) slabs can guide light perfectly for optical modes below the light cone [33]. Above the light cone, guided modes become guided resonances since they are leaky [34,35]. BICs can ex-ist as isolated points on the bands of guide resonances [5][6][7][8][9][10][11][12][13]. From the far-field viewpoint, they can be interpreted as the vortex singularities of far-field polarizations with quantized topological charges [9]. These topological charges can be created, annihilated and merged in the Brillouin zone [12][13][14][15]. From the viewpoint of wave interference, some BICs in PhC slabs can be treated as the Friedrich-Wintgen type that originates from the destructive interference of two different guided resonances [36][37][38].
We proposed that the formation mechanism of BICs in a PhC slab can be further interpreted in terms of the interference of bulk Bloch states [10]. For a uniform dielectric slab, the formation of guided waves requires two conditions: a total internal reflection (TIR) at the interfaces between the slab and free space and that waves along the direction perpendicular to the slab are standing ones. The formations of BICs in a PhC slab must also satisfy these C The Author(s) 2022. Published by Oxford University Press on behalf of China Science Publishing & Media Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. two conditions. In a uniform dielectric slab, the condition of TIR is simply that the angle between the propagating direction and the slab-surface normal is greater than the critical angle. However, any optical mode supported in a PhC slab is the superposition of bulk Bloch waves of such an infinite PhC rather than a single plane wave. As a result, the TIR condition for a PhC slab becomes that of the TIR of constituent Bloch waves. If the total transmission for an optical mode consisting of multiple Bloch waves from the PhC slab side to the free-space side vanishes, TIR will occur. It is just the condition of the TIR of Bloch waves stemming from the multiple interference of the constituent Bloch waves. Therefore, the study of BICs can start from a basic problem: the diffraction of Bloch waves at a single interface.
The key point is that there may exist multiple reflected and refracted waves because of Bragg scattering [33,39]. At the interface between a uniform dielectric and free space, the wave vector component k parallel to the interface is a good quantum number due to the continuous translational symmetry at the interface. TIR can occur when |k | is larger than the free-space wave vector since the perpendicular component of the wave vector k ⊥ on the freespace side becomes imaginary. In PhCs, the continuous translational symmetry is broken. However, the discrete translational symmetry leads to the equivalence of k and k + nG, where n is an integer and G is a reciprocal lattice vector. This new degree of freedom renders the TIR of Bloch waves possible via a coherent way [10,21] to be discussed in detail later.
Here, the TIR of Bloch waves is fully investigated from the viewpoint of diffraction. The generalized Fresnel equations for Bloch waves are derived and formulas for the TIR of two Bloch waves with a very compact form are obtained analytically. For PhC slabs, the conventional conditions for the existence of waveguide modes can be directly generalized based on the TIR of Bloch waves and the solutions of the generalized conditions are exactly BICs. A BIC solver is therefore developed with low computational complexity and fast convergence speed, and can be used for the search and determination of BICs in a very large parameter space. Different from previous studies of BICs in PhC slabs, which are restricted to a single radiation channel, the generalized conditions can be also applied to the case of multiple radiation channels. Therefore, the BIC solver can find BICs for any number of radiation channels at any high frequency. Examples of BICs with two radiation channels are given and it is demonstrated that multi-channel BICs require the coincidence of the topological charges of far-filed radiations in all radiation channels.

Theory for the TIR of Bloch waves
The TIR of Bloch waves can be interpreted from the perspective of diffraction. We start from a simple planar grating shown in Fig. 1a. If there is only 0th diffraction order, the direct transmission is not zero in general. Therefore the simplest non-trivial case is that there are two propagating diffraction orders with wave vectors k and k − G. It is the discrete translational symmetry of the grating that leads to the Bragg scattering between k and k + nG, offering a higher degree of freedom to control the incident waves that is not restricted to a single plane wave with a fixed k. To be specific, if a zero total transmission can occur by introducing two incident plane waves with k and k − G, TIR is thus realized in a very unusual way even though the transmission for each of the incident wave is not zero. In Fig. 1a, diffraction of incident plane waves with wave vector k (purple arrows) and k − G (red arrows) is shown. For both cases, there are 0th and −1st (or 1st) diffraction orders. By definition, the 0th order remains the same wave vector as the incident one, so the order of the principal maximum and the secondary maximum is in fact switched for these two incident waves. If the two plane waves k and k − G are incident onto the grating simultaneously, the elimination of transmission shown in the right figure of Fig. 1a requires that the intensities of the diffracted waves satisfy the relation: where I k (m) represents the intensity of the mthorder diffraction for the incident wave vector k. However, according to the Fraunhofer diffraction from a diffraction grating [40], the principal maximum of m = 0 is usually the dominant maximum and Equation (1) cannot be satisfied generally. If a high-index (n h ) material is adopted as a substrate shown in Fig. 1b, we can possibly make the propagating diffraction order with wave vector k b − G evanescent. Here, k b is the wave vector on the transmission side with refractive index n b and k b, = k . The condition is that k satisfies n b k 0 < |k − G| < n h k 0 , where k 0 is the free-space wave vector. Under this condition, only one propagating diffraction order survives for both the incidence of k and k − G, and Equation (1) is relaxed greatly and reduced to I k (0) = I k−G (1). The destructive interference of the transmitted waves can be readily realized just by appropriately choosing the relative phase and amplitude of the two incident waves. Note that the essential point is that we have a sufficient degree of freedom for the incident waves to Natl Sci Rev, 2023, Vol. 10, nwac043 cancel out the transmission. Similar destructive interference was considered to achieve some unique phenomena such as complete reflections [41][42][43] and perfect anti-reflections [44]. In fact, the combination of a grating and a highindex material can be replaced by a PhC [33], as shown in Fig. 1c. We focus on a 1D semi-infinite PhC with a period of a in the x direction and uniform in the y direction. The alternating dielectric layers in the PhC have relative permittivity ε 1 and ε 2 , and thicknesses a-d and d, respectively. The background is chosen to be air with ε b = 1. Different from that in the Fraunhofer diffraction of gratings, we choose the incident waves as Bloch waves rather than plane waves since Bloch waves are eigenstates of the periodic structure and any optical modes supported are a superposition of these Bloch waves. A Bloch wave with wave vector k + nG is equivalent to that with k. Obviously, all of the arguments above for the existence of TIR can be applied to Bloch waves here.
The existence of multiple Bloch waves can be clearly seen in isofrequency contours. The dispersion relation, which relates the frequency ω, the normal wave vector k z and the Bloch wave vector k x , is given in the Supplementary Information. Figure 2 shows three examples with different frequencies ω. The isofrequency contours in air are shown by black lines, whereas those for the PhC are indicated by red lines. The red lines will be folded back when they go beyond the first Brillouin zone (see the dashed red lines) due to the periodicity in the x direction.
In Fig. 2a, at a low frequency, the isofrequency contours for the PhC and air in the first Brillouin zone are simply two circles without band folding, offering conventional refraction and transmission. The coefficient for the incident wave is denoted by a 1 , whereas that for the reflected and transmitted waves are denoted by r 1 and t 0 , respectively. The subscript in t 0 stands for 0th-order diffraction in air. In this case, the number of propagating Bloch waves in the PhC (N p ) and that of propagating diffraction orders in air (N r ) are equal to 1. The propagating diffraction orders can also be viewed as radiation channels. As frequency increases, band folding takes place and band gaps appear at the edges of the first Brillouin zone, as shown in Fig. 2b, and N p of the propagating Bloch waves is increased to 2, whereas N r of the radiation channels is still 1. Thus, for a single Bloch wave incident with a coefficient a 1 , in addition to the reflection r 1 (the same Bloch wave), an additional Bloch wave with a coefficient r 2 will also be excited. When frequency is further increased, more Bloch waves will be present as additional reflected waves, such as r 3 shown in Fig. 2c. Moreover, when frequency goes beyond the diffraction limit, the −1storder diffracted wave in air will change from evanescent to propagating, so that for the case in Fig. 2c, we have N p = 3 and N r = 2.
The formalism for the diffraction of Bloch waves incident from a PhC to air are outlined as follows.
Here, we only consider transverse electric (TE) Bloch waves (E = E yŷ , H = H xx + H zẑ ). Transverse magnetic (TM) Bloch waves (H = H yŷ , E = E xx + E zẑ ) are discussed in the Supplementary Information. Suppose that a series of TE Bloch waves with a fixed frequency ω and Bloch wave vector k x impinge on the PhC/air interface at z = 0. The electric field inside the PhC can be written as follows: where a n and r n are respectively the complex coefficients of the incident and reflected Bloch waves, k (n) z is the normal wave vector of the nth Bloch wave and u (n) (x) is the periodic-in-cell part of the nth Bloch wave function. Suppose we have N p incident propagating Bloch waves 1 ≤ n ≤ N p . Bloch waves with the order n > N p are evanescent waves with k z being purely imaginary. Physically, incident evanescent waves that increase away from the interface should be excluded for this semi-infinite PhC, namely a n>N p = 0. The transmitted wave in air can be expressed as follows: Natl Sci Rev, 2023, Vol. 10, nwac043 where k x,m = k x + mG and k z,m = k 2 0 − k 2 x,m . Here, t m is the complex transmission coefficient for the mth diffraction order. At the interface, a Fourier transform of the boundary conditions (the continuity of tangential E and H fields) gives: and where (T) m = t m , (A) n = a n , (R) n = r n , Note that the first N r elements of T correspond to the radiation channels in air. To solve Equations (4) and (5), the number of Fourier components (indexed by m) should be chosen to be the same as the number of Bloch waves (indexed by n). By eliminating T, the relation between the reflection and incidence reads: (6) Then, the transmission can be expressed as follows: where ↔ I is the identity matrix. Since we are only interested in the transmission of the radiation channels (denoted by T r ), Equation (7) can be reduced to: Equations (6) and (8) are in fact the generalized Fresnel equations for Bloch waves. Based on these two equations, the problem of the incidence of any number of Bloch waves can be solved. Obviously, the TIR condition of Bloch waves is also a direct consequence, given by: When the number of propagating Bloch waves is equal to that of radiation channels, namely N p = N r , a non-trivial solution of this condition requires that det( ↔ Mr) = 0, which is difficult to realize for a PhC. However, if N p > N r , a non-trivial solution of incidence A p always exists for Bloch waves.
Based on the TIR of Bloch waves, light can be further guided in a PhC slab with a finite thickness h. Distinct from the semi-infinite PhC, all evanescent Bloch waves are allowed in a PhC slab, with either positive or negative attenuation in the z direction. The origin of the z-axis is now set at the center of the PhC slab for convenience. Equations (6) and (7) are also slightly modified via replacing a n with a n e −i k (n) z h . Supposing that the TIR condition T r = 0 is satisfied at the upper interface for some properly initiated incidence a n , the reflected waves will then become the incident waves at the lower interface. In the case that the ratio r n /a n remains a constant for any arbitrary n, the TIR condition can be maintained at the lower interface.
However, the TIR condition is not the only condition for forming a waveguide mode. The phase accumulated after a round trip should be integer multiples of 2π , also called the guidance condition [45]. This guidance condition can be directly generalized just by counting the accumulated phase for each Bloch wave. At the interface of the PhC slab, a phase shift ϕ (n) r = arg(r n /a n e −i k (n) z h ) takes place for the nth Bloch wave. Note that the additional term e −i k (n) z h in the phase shift comes from the shift of the origin of z compared with the above semi-infinite PhC case. Similar to that for conventional waveguide modes, the total phase change for a round trip should be integer multiples of 2π for the nth Bloch wave, which can be simply expressed as follows: where m (n) is an integer. Equations (9) and (10) can be viewed as the generalized conditions for waveguide modes in a PhC slab, as summarized in Table 1.
Waveguide modes that satisfy the generalized conditions inside the light cone are precisely BICs. Combining Equations (9) and (10), we also obtain that r n /a n = ±1 for all Bloch waves, where the positive and negative signs correspond to even and odd symmetries in the z direction, respectively. The generalized conditions for waveguide modes can be used to efficiently determine BICs in the k x -ω space. In addition to propagating Bloch waves, evanescent waves with purely imaginary k z can exist near the interface of the PhC slab and should also be taken into account. Based on the generalized conditions for waveguide modes, a BIC solver has been designed [46] with the advantage of very low computational complexity and fast convergence speed. Since the Bloch waves we adopt form an appropriate basis set inside the PhC, the positions of BICs in the k xω space converge very quickly if only a few evanescent waves are considered, in addition to the propagating Bloch waves (see Supplementary Fig. 1). To be specific, we first ensure that the TIR condition is satisfied at one of the interfaces for every (k x , ω) point. The TIR condition requires that the number of propagating Bloch waves is larger than that of radiation channels (N p > N r ). The corresponding phase shift, ϕ (n) r , at this interface can be obtained by solving Equations (6) and (9) (see Supplementary Information for details). Second, we build a database of ϕ (n) r for a PhC in the whole k x -ω space. Finally, for any thickness h, the total phase of a round trip for the nth Bloch wave inside the PhC slab is simply k r . What the solver should do is to determine whether this phase is integer multiples of π . Therefore, the computational time is mainly spent on the construction of a reflection-phase database. With this database, the time to search BICs for different h values is negligible. We show an example of searching BICs in a range of k x -ω space with N p = 2 and N r = 1 in Supplementary Fig. 1, in full agree-ment with the results simulated by the finite element method.
It is worth mentioning that the algorithm based on the generalized conditions for waveguide modes can be applied to not only the k x -axis but the whole Brillouin zone. As a result, this BIC solver can work in the whole k -ω space, where k = (k x , k y ). However, BICs usually exist on high-symmetry lines. In Supplementary Fig. 2, we give another example of searching BICs in the k y -axis in the parameter space.

TIR of two propagating Bloch waves
The simplest case for TIR of Bloch waves is that there are only two propagating Bloch waves in the PhC (N p = 2) and one radiation channel in air (N r = 1), as shown in Fig. 1c. In this case, an analytical solution can be obtained. We assume that only propagating Bloch waves are considered and other evanescent waves are neglected. According to the above analysis, we only need to achieve the TIR of two Bloch waves at a single interface and then adopt the generalized guidance condition to fix BICs. The generalized Fresnel equations for Bloch waves can be simplified considerably and a concise form for the relative coefficient of the incident waves can be directly obtained when TIR occurs at the interface. For TE waves, it can be expressed as (see Supplementary Information for details): where Z n = k z,−1 /k (n) z has a similar form of relative surface impedance [47]. Here, k z,−1 is the normal wave vector of −1st-order diffracted wave in air, which is purely imaginary. Moreover, the reflection coefficients at the interface are as follows: which takes a similar form of the reflection coefficient in the conventional Fresnel equations. Note that Equation (12) holds only when the TIR of two Bloch waves occurs. The TIR condition becomes slightly complicated for two TM Bloch waves since the electric field is a vector in nature for the TM case but is a scalar for the TE case [48]. An approximate form of the Fourier transform of ε(x) is used: ε −1 (x) ∼ κ 0 + κ 1 e i G x + κ −1 e −i G x . When the TIR of two TM Bloch waves occurs, similar forms of the relative incidence and reflection coefficients can be obtained as those in Equations (11) and (12). However, the definition of Z n should be modified and expressed as follows (see Natl Sci Rev, 2023, Vol. 10, nwac043 Supplementary Information for details): . The TIR condition of two Bloch waves can be realized via Equations (11) and (12) with appropriate definitions of Z n for TE and TM waves. The phase shift of TIR for the nth Bloch wave is as follows: Note that diffracted evanescent waves are neglected in the above TIR condition, so Equations (11)(12)(13)(14) work for the case when the index contrast is not too large, namely = |ε 2 − ε 1 |/ε 1 1. Strikingly, even when → 0, i.e. the index contrast is vanishingly small, Z n approaches a constant for any (k x , ω) point. Therefore, two important conclusions can be drawn. First, BICs obtained from Equations (10) and (14) approach a series of fixed points in the kω space [21]. Generally, the band of guided resonances can be regarded as the folded band of the waveguide modes in an effectively uniform waveguide. The existence of discrete BIC points in the limiting case manifests the non-trivial physical consequence that continuous translational symmetry is broken into a discrete translational symmetry even if → 0. Second, it is known that the introduction of a substrate can destroy BICs [49]. This is because the TIR conditions at the upper and lower interfaces are different; the combination of TIR at a single interface and guidance condition in Table 1 cannot restore the waveguide mode after a round trip. Or, in other words, the TIR conditions at the two interfaces contradict each other if there is a substrate.

Multi-channel BICs
When frequency increases, more than one propagating diffraction order (i.e. radiation channel) in air appears, shown in Fig. 2c. The construction of BICs is more subtle since all radiation channels should be closed. Note that multi-channel BICs occurring at k x = 0 or π /a were discussed in Ref. [50]. However, since these multi-channel BICs appear at the high-symmetry points in the Brillouin zone, the corresponding radiation channels are not completely independent. Strikingly, the above generalized conditions for waveguide modes can be directly applied to the case with multiple radiation channels. The BIC solver we designed can thus work well to determine multi-channel BICs. Two different types of BICs with two radiation channels are taken as examples and shown in Fig. 3a and b, which consist respectively of three and four propagating Bloch modes. These two multi-channel BICs appear on the TE (1) 0 and TE (−2) 0 bands, as highlighted by red dots in Fig. 3. Here, TE (m) 0 represents the fundamental TE mode with m being the index of the band folding in the reduced-zone scheme. It is known that BICs interpreted by topological vortexes can exist robustly in the parameter space [9]. The robustness of BICs should be reexamined for multi-channel BICs since they only exist for some specific thicknessesfor example, h BIC = 1.948a and 1.968a in Fig. 3a and b, respectively. The Q factors of the guided resonance modes near the multi-channel BICs are plotted in Fig. 3c and d for different h values. It can be clearly seen that the Q factor diverges only when the thickness is equal to h BIC . This divergence behavior disappears as long as the thickness is slightly varied away from h BIC , which is distinct from robust BICs below the diffraction limit. The divergence rates are also plotted in Fig. 3e and f, which are Q ∼ 1/δk 2 x and Q ∼ 1/δh 2 (inverse square law), respectively.
It has been demonstrated that BICs below the diffraction limit are vortex centers of the polarization directions of far-field radiations [9], characterized by topological charges and robust then in parameter space. However, for multi-channel BICs, an increased number of radiation channels can make an essential difference and the topological nature is manifested in other ways. To reveal the topological nature, the far-field polarization states are investigated for each radiation channel (see Supplementary Information for details). The far-field polarization states displayed in Fig. 4c correspond to the multi-channel BIC shown in Fig. 3a. There are two radiation channels in air and three propagating Bloch waves in the PhC, as shown in Fig. 4a. The total Q factor, defined by Q = (1/Q 0 + 1/Q −1 ) −1 , takes into account the radiative losses from the 0thorder diffraction (Q 0 ) and −1st-order diffraction (Q −1 ). In the upper and lower panels of Fig. 4, Q 0 and Q −1 are plotted separately as purple and green lines, respectively, and the polarization states of the 0th-and −1st-order diffraction are also shown correspondingly. Since Q diverges at the multi-channel BIC for the thickness h = 1.948a, both Q 0 and Q −1 have to diverge simultaneously. First, this implies that there is one topological charge (marked by the black dot) in both two polarization maps as shown in Fig. 4c; second, the two topological charges coincide with each other in momentum space, giving rise to a multi-channel BIC without any leakage. Note that the topological charge is defined by v m = (1/2π ) L d k · ∇ k φ m (k ). Here, L is a closed loop in momentum space surrounding the It is worth emphasizing that these two topological charges come from the same eigenstate with fixed k and ω but belong to different radiation channels (i.e. the propagating diffraction orders with k and k −G). Therefore, they are independent and will not merge or annihilate each other in momentum space. This topological property is distinct from that of merging BIC [13,31,38], wherein the topological charges are linked to the same radiation channel.
When the thickness of the PhC slab is slightly varied away from h BIC , the multi-channel BIC no longer exists, as shown in Fig. 4b and d, and Q 0 is bounded, whereas Q −1 still diverges at a certain k x . The topological charge v 0 = +1 for the 0th-order diffraction splits into two half-integer charges of 1/2 with the total topological charge conserved and each being circularly polarized. Because of the y-mirror symmetry of the system, the two circularly polarized states are symmetric about the k x -axis and carry the same charge with different handedness (or chirality). The states with right-handed circular polarization (RCP) and left-handed circular polarization (LCP) are indicated by red and blue dots, respectively, in the upper panels of Fig. 4b and d. The splitting of an integer charge into two half-integer charges here comes only from the change of thickness and the symmetry of the system is perfectly maintained. Note that below the diffraction limit, the breaking of the C 2 symmetry is necessary in order to observe this kind of splitting [12]. This non-symmetry-breakinginduced splitting manifests the unusual topological nature for multi-channel BICs. For the −1storder diffraction, the topological charge persists and slightly moves along the k x -axis (see the lower panel in Fig. 4). The dotted arrows in Fig. 4 are a guide for the eyes and indicate the evolution of topological charges. The half-integer charge of RCP (red point) passes through the k x -axis from positive to negative k y , while the one of LCP (blue point) passes through the k x -axis from negative to positive k y . The two half-integer charges meet each other at the k xaxis. Multi-channel BICs lying in the k x -ω space with only three propagating Bloch modes can be understood as the coincidence of two integer charges in momentum space: one coming from the merging of two half-integer charges and the other being a stable integer charge moving on the k x -axis slowly.
Multi-channel BICs can even manifest a different topological nature if they lie in the region of the k x -ω space with different numbers of propagating Bloch modes. Another example, the multi-channel BIC marked in Fig. 3b, is demonstrated by showing the far-field polarization states of the 0th-and −1storder diffraction separately in Fig. 5c. Similarly, both Q 0 and Q −1 diverge at this BIC point and the two topological charges coincide with each other in momentum space so that leakage is eliminated for these two radiation channels. Note that the two topological vortexes defined in the two radiation channels can either exhibit the same amount of charge, as shown in Fig. 5c, or different amounts of charge, as shown in Fig. 4c. Furthermore, both integer charges in Fig. 5c will split into a pair of half-integer charges of 1/2 with opposite chirality when the thickness of the PhC slab is slightly varied from h BIC , as shown in Fig. 5b and d. This non-symmetry-breaking-induced splitting of an integer charge into two half-integer charges is a generic phenomenon since it happens in both radiation channels. In short, multichannel BICs in the region with N p = 4 and N r = 2 can also be interpreted as the coincident point of two integer charges in momentum space, both of which result from the merging of two half-integer charges.

CONCLUSION
In summary, we derive the generalized Fresnel equations for the Bloch waves at a PhC/air interface, from which the TIR condition of Bloch waves are obtained. For a PhC slab, by combining the TIR of Bloch waves and the guidance condition, the generalized conditions for waveguide modes are given, with solutions being precisely the BICs. Distinct from BICs below the diffraction limit, multi-channel BICs with frequencies beyond the diffraction limit are found which can only exist for some specific geometric parameters of the PhC slab. They possess a quite different topological nature stemming from the coincidence of two integer charges in the polarization maps of two different radiation channels. Integer topological charges can split into two half-integer charges even without breaking any symmetry, which is generic for multi-channel BICs. Our BIC solver with the generalized conditions for waveguiding in PhC slabs incorporated offers a powerful tool for readily finding BICs at any frequency in momentum space. The distinct topological nature revealed in multi-channel BICs from conventional ones may render new opportunities in designs and applications of BICs possible in nanophotonics and enhanced light-matter interactions as well.