BIC in waveguide arrays within a symmetry classification scheme

A counterintuitive concept of bound states in the continuum (BICs) was initially proposed by von Neumann and Wigner neumann as a mathematical construct represented by special states which lie in the continuum but remain localized without any radiation as a consequence of confining potential. By inspecting various scenarios of light scattering and their relation to location of the poles and zeros of the scattering matrix, BICs were identified as the limiting case of general Hermitian scatterer, when the pole and zero coalesce on the real axis with mutual coherent destruction monticone . BICs have been found in a wide range of different material platforms and have been experimentally observed in electromagnetic, acoustic and water waves Hsu -ysk .

In this Letter we employ a theoretical model PK_pra to investigate the coupled system consisting of an infinite array of single-mode optical waveguides, with two additional waveguides above and below the array (Fig. 1). Our configuration can be viewed as an extension of a photonic implementation of the Fano-Anderson model Mahan with a single additional waveguide ysk2 . This structure can be regarded as generalization of a famous finite system in which an optical symmetry-protected BIC was confirmed experimentally plotnik . Surprisingly, its comprehensive theoretical understanding has not yet been presented. We note that a similar structure consisting of a discrete linear chain, which in addition to Fano defects contains one $\delta$-like defect which gives rise to the Fano-Feshbach resonance, was considered in nguyen .

To explore physics and to provide a theoretical analysis of the formation of the symmetry-protected BIC, we determine analytically and calculate numerically spectrum of eigenmodes as well as scattering properties of the waveguide array. We establish the conditions, given by the relation between the self-coupling coefficients associated with side-coupled waveguides $\delta_{\pm 1}$ and the change of the propagation constant $\varepsilon$, that define regimes in which the governing coupled dynamic equations support the different eigenmode spectra. The eigenvalues $\varepsilon$ reveal new and interesting features, such as the level repulsion and exceptional points (EPs), that have been found by tracing their position within the complex plane with increasing difference between $\delta_{1}$ and $\delta_{-1}$. All the modes supported by the structure manifest themselves in the scattering spectra as the resonances that can be identified by means of the generalized Weierstrass factorization approach bonod1 ; bonod2 . We found that the resonance associated with quasi-BIC state possesses a Fano line-shape arising from the interference between quasi-BIC and another leaky mode. The transmittance near the resonance can be rewritten into the form of the Fano formula, where the shape parameter $f$ can be expressed in terms of the poles associated with the interacting modes. The knowledge of a set of discrete singular points associated with the scattering matrix, provides a widely applicable framework to interpret and to engineer the resonant properties of the nanostructures bonod1 ; bonod2 .

I. Model. To describe the coupled structure (Fig. 1) we use the standard approach based on the coupled mode theory (CMT) hardy ; shi . The total field is represented with slowly varying modal amplitudes $\psi_{m}(z)$ ($m\in Z$) and $\varphi_{\mp 1}(z)$ on the $m-$th waveguide in the array and on the upper/lower additional waveguide, respectively. The evolution of the amplitudes is described by a set of coupled equations

	$$\displaystyle i\psi_{m}^{\prime}=C\left(\psi_{m-1}+\psi_{m+1}\right),\;m\neq 0$$		(1)
	$$\displaystyle i\psi_{0}^{\prime}=\delta_{0}\psi_{0}+C\left(\psi_{-1}+\psi_{1}\right)+C_{1}\varphi_{-1}+C_{2}\varphi_{1}\;$$		(2)
	$$\displaystyle i\varphi_{-1}^{\prime}=\delta_{-1}\varphi_{-1}+C_{1}\psi_{0}$$		(3)
	$$\displaystyle i\varphi_{1}^{\prime}=\delta_{1}\varphi_{1}+C_{2}\psi_{0},$$		(4)

where the prime stands for the derivative according to the $z$, $C$ is the coupling coefficient in the array, $C_{1,2}$ are the coupling coefficients between the upper/lower additional waveguide and $0$-th waveguide in the array, and $\delta_{0}$ and $\delta_{\mp 1}$ is self-coupling coefficients in $0$-th waveguide in the array and in the upper/lower additional waveguide, respectively. The self-coupling coefficients may include the effect of nearest-neighbours due to the geometry as well as the possible external perturbations of $0$-th waveguide and the additional waveguides. We assume that all the parameters which enter Eqs. (1)-(4) are real.

The unperturbed array supports propagation of Bloch waves $\psi_{m}=A\exp\left(-ik_{x}am-i\varepsilon z\right)$, where $k_{x}$ is the Bloch wavenumber, $a$ is the period of the lattice and the “energy” (change of the propagation constant) $\varepsilon$ fulfils the well-known dispersion relation

$$\varepsilon=2C\cos\left(k_{x}a\right).$$

(5)

II. Scattering problem. The additional waveguides cause scattering of the Bloch waves. To evaluate the effect, we follow the standard approach, see, e.g., Weimann ; ysk2 , and represent the field in the left part of the periodic array in terms of incident and reflected waves

$$\psi_{m}=e^{-ik_{x}am-i\varepsilon z}+re^{ik_{x}am-i\varepsilon z},\;m<0$$

(6)

in the right part as the transmitted wave

$$\psi_{m}=te^{-ik_{x}am-i\varepsilon z},\;m>0.$$

(7)

and the fields in the additional waveguides as

$$\varphi_{\pm 1}=B_{\pm 1}e^{-i\varepsilon z}$$

(8)

After substituting Eqs. (6)-(8) into Eqs. (1)-(4) and few straightforward steps one finds that Eqs. (6), (7) are valid also for $m=0$ and obtains the analytical expressions for the amplitude transmittance $t$ and reflectance $r$ as

$$t=r+1=\frac{2C\sin(k_{x}a)}{2C\sin(k_{x}a)+i\mu}$$

(9)

where $\mu$ reads

$$\mu=\delta_{0}+\frac{C_{1}^{2}}{\varepsilon-\delta_{-1}}+\frac{C_{2}^{2}}{\varepsilon-\delta_{1}}.$$

(10)

Figure 2 demonstrates typical behaviour of the intensity transmittance $T=|t|^{2}$ evaluated for various structural parameters. Figs. 2 (a), (b) show results for symmetric structures where the self-coupling coefficients in both additional waveguides are the same $\delta_{1}=\delta_{-1}$. The most dominant feature is vanishing $T$ at single point $\varepsilon=\delta_{\pm 1}$ (for fixed $k_{x}$, $0<k_{x}a<\pi$) due to diverging parameter $|\mu|\rightarrow\infty$ in Eq. (9), see also Eq. (10). This resonance corresponds to existence of the symmetry protected BIC and, for varying $k_{x}$, is indeed located at positions given by the dispersion relation (5) [see dashed lines in Figs. 2 (a), (b)]. Note also, that the zero transmittance at $\varepsilon=\delta_{\pm 1}$ is double degenerate and identical with energies of localized states of the additional waveguides, i.e., when $\varepsilon$ is calculated from Eqs. (3) and (4) for $C_{1,2}=0$.

By breaking the vertical symmetry, the degeneracy is lifted as is illustrated in Figs. 2 (c), (d) and Fig. 3 for asymmetric structures defined as $\delta=-\delta_{-1}=\delta_{1}$, where the parameter $\delta$ describes an asymmetry strength. We observe two different zeros with energies equal to those of states of isolated additional waveguides, $\varepsilon=\pm\delta$. It will be shown further, that the resonances correspond to the excitation of two leaky modes of the composed structure, one of them being quasi-BIC, and the related spectral features may be explained by interference of the two leaky modes that leads to the Fano profile.

Importantly, the zero positions observed in Figs. 2, 3 do not depend on the coupling coefficients $C_{1,2}$. On the other hand, $C_{1,2}$ have profound influence on the shape of the resonances: decreasing of $C_{1}/C$ and $C_{2}/C$ (weaker coupling) leads to narrower dip in the transmittance $T$, while asymmetry in $C_{1,2}$ strongly affects asymmetric behaviour of the two dips [compare Figs. 3 (a) and (b)].

III. Modes. The spectral features described above can be explained in terms of the modes of the composed structure. The modal field in the linear chain has the form

$$\psi_{m}=Ae^{-ik_{x}a|m|-i\varepsilon z},\\$$

(11)

while that in the additional waveguides is given by Eq. (8). By substituting Eqs. (11) and (8) into Eqs. (2)-(4) we obtain the following nonlinear eigenvalue problem

	$$\displaystyle\left[\delta_{0}-\varepsilon+2C\exp\left(-ik_{x}a\right)\right]A+C_{1}B_{-1}+C_{2}B_{1}=0,$$		(12)
	$$\displaystyle C_{1}A+\left(\delta_{-1}-\varepsilon\right)B_{-1}=0,$$		(13)
	$$\displaystyle C_{2}A+\left(\delta_{1}-\varepsilon\right)B_{1}=0.$$		(14)

The relation between the self-coupling coefficients $\delta_{\pm 1}$ and change of the propagation constant $\varepsilon$ plays the key role in a symmetry classification scheme, which defines the regimes with different families of the eigenmodes.

A. We first consider the case when $\varepsilon=\delta_{1}$ or $\varepsilon=\delta_{-1}$. According to Eqs. (12)-(14) a nontrivvial solution is not supported unless the structure possesses the vertical symmetry $\delta_{-1}=\delta_{1}$ when one obtains a doubly-degenerate bound state $\varepsilon=\delta_{\pm 1}$, which yields a vertically antisymmetric eigenfunction $C_{1}B_{-1}=-{C_{2}}B_{1}$, $A=0$ and represents the symmetry protected BIC provided $-2C<\delta_{\pm 1}<2C$.

B. Now we examine the case $\varepsilon\neq\delta_{1}$ and $\varepsilon\neq\delta_{-1}$. Eqs. (13)-(14) yield nontrivial solutions for the amplitudes associated with the attached waveguides in the form

$$B_{\mp 1}=\frac{C_{1,2}}{\varepsilon-\delta_{\mp 1}}A.$$

(15)

By substituting the both expressions into Eq. (12) one obtains the eigenvalue equation

$$2C\sin(k_{x}a)+i\mu=0.$$

(16)

It is important to keep in mind, that none of the solutions of Eq. (16) is a true BIC as we have $A\neq 0$, however, for the vertically asymmetric structures ($\delta_{-1}\neq\delta_{1}$), one of them may represent a quasi-BIC as will be shown below.

B.1 Symmetric states. The family of the symmetric states can be divided to several groups according to the vertical symmetry of the structure. B.1a We first consider the structure with the trivial vertical symmetry, $\delta_{-1}=\delta_{1}=\delta_{0}=0$, for which Eq. (16) provides 4 eigenvalues

$$\varepsilon=\pm\sqrt{2C^{2}\pm 2\sqrt{C^{4}+C_{\rm av}^{4}}},$$

(17)

where all combinations of the signs are allowed and $C_{\rm av}^{2}=(C_{1}^{2}+C_{2}^{2})/2$. The modes (in contrast to antisymmetric modes supported by the same vertically symmetric structures) exhibit symmetry

$$C_{2}B_{-1}={C_{1}}B_{1}.$$

(18)

The full spectrum of the trivially symmetric structures, shown in Fig. 4(a),(b), consists of symmetry protected doubly-degenerate BIC within the continuum, two bound modes (B1,B2) with real energy, which bifurcate from the lower and upper edges of the continuum when the parameter $C_{\rm av}$ is increased, and two modes (L,L’) with pure imaginary energy. We note that Eq. (17) illustrates the general feature of Eq. (16): the energy eigenvalues are either real or occur in complex conjugate pairs. The states with $\mathrm{Im}(\varepsilon)<0$ (solid lines in Figs. 4 and 5) fulfil the conventional definition of leaky modes and satisfy outgoing boundary conditions $0<\mathrm{Re}(k_{x}a)\leq\pi$, states with $\mathrm{Im}(\varepsilon)>0$ (dashed lines in Figs. 4 and 5) correspond to ingoing boundary conditions $-\pi<\mathrm{Re}(k_{x}a)<0$ PK_pra .

B.1b For the structure belonging to the class characterized by the general vertical symmetry $\delta_{-1}=\delta_{1}\neq 0$, Eq. (16) yields 4 solutions, the modes are symmetric according to Eq. (18). The spectrum shown in Fig. 4(c),(d) also consists of the bound modes (B1, B2, BIC); the latter one enters the gap at $\delta_{\pm 1}/C=2$ when $\delta_{\pm 1}/C$ is increased. The L, L’ modes initially exhibit complex conjugate energies $\varepsilon$, the real component of their energy $\mathrm{Re}(\varepsilon)$ also enters the gap at $\delta_{\pm 1}/C=2.5$. The imaginary component of their energy $\mathrm{Im}(\varepsilon)$ remains nonzero in the range $2.5<\delta/C<3.1$ until the (L,L’) doublet reaches the EP beyond which the level repulsion effect occurs. Here the level repulsion, in contrast to the $\cal{PT}$-symmetric structures, arises from the coupling between two degenerate modes with positive energy bernier .

B.2 Structures with broken vertical symmetry (asymmetric states). Breaking the vertical symmetry, $\delta_{-1}\neq\delta_{1}$, yields coupling of the true BIC antisymmetric state to the continuum and the doubly degenerate mode is split into quasi-BIC (qBIC) and qBIC’ mode as shown in Fig. 5, which demonstrates a rich spectral behavior of the structure. Eq. (16) yields 6 solutions; instead of Eq. (18), the amplitudes $B_{\pm 1}$ satisfy the condition $C_{2}B_{-1}(\delta_{-1}-\varepsilon)=C_{1}B_{1}(\delta_{1}-\varepsilon)$.

The spectra for the structures with $C_{1}=C_{2}$, see Figs. 5(a),(b), reveal in comparison with those shown in Fig. 4 the EP which occurs when asymmetry strength $\delta/C=0.5$. The behavior of the EP differs from common one known in $\cal{PT}$-symmetric structures: the L, L’ and qBIC, qBIC’ doublets, initially (for $\delta/C<0.5$) with complex conjugate values of energy $\varepsilon$ and vanishing its real component, $\mathrm{Re}(\varepsilon)=0$, coalesce at the EP in L, qBIC and L’, qBIC’ pairs. In the range between $0.5<\delta/C<3.1$, the common value of $\mathrm{Re}(\varepsilon)$ for the L, L’ (or qBIC, qBIC’) doublet increases (decreases) and approaches the upper (lower) edge of the continuum, while $\mathrm{Im}(\varepsilon)$ for the L, qBIC and L’, qBIC’ pairs coalesce. Beyond another EP at $\delta/C=3.1$ the L, L’ and qBIC, qBIC’ doublets split into the singlets and transform into the bound states.

The structures with nonequal coupling coefficients $C_{1}\neq C_{2}$ exhibit, when increasing $\delta/C$ from zero thresholdless behavior, see Fig. 5(c),(d), i.e., the common value of $\mathrm{Re}(\varepsilon)$ for the L, L’ (or qBIC, qBIC’) doublets starts to increase (decrease) with $\delta/C$ and enters the gap above (below) the continuum. With further increasing $\delta/C$ the doublets L, L’ and qBIC, qBIC’ reach EPs at which they also transform into the bound states with $\mathrm{Im}(\varepsilon)=0$.

IV. Factorization approach. The scattering properties of an arbitrary structure can be characterized by the scattering matrix, which for our structure has two eigenvalues PK_pra $S_{o}=r-t=-1$ and

$$S_{e}=r+t=2t-1=\frac{2C\sin(k_{x}a)-i\mu}{2C\sin(k_{x}a)+i\mu}.$$

(19)

The dependence $S_{e}(q)$, $q\equiv k_{x}a$, can be expressed through the Weierstrass factorization theorem bonod1 ; bonod2 which, for our structure, takes the simple form PK_pra

$$S_{e}(q)=\prod\limits_{n}{\frac{q-p_{n}^{\ast}}{q-p_{n}}},$$

(20)

where $p_{n}$ are poles of $S_{e}(q)$. It follows from Eq. (19) that the poles $p_{n}$ are solutions of the eigenvalue problem (16) in terms of $k_{x}a$. Therefore, the knowledge of eigenmodes is sufficient to restore $S_{e}(q)$ through Eq. (20) and to evaluate the scattering spectra as $r=(S_{e}-1)/2$ and $t=(S_{e}+1)/2$, where, in the reflectance $R=|r|^{2}$, each pole creates a peak of Lorentzian shape bonod1 . This way, we can interpret observed spectral features in terms of the eigenmodes. The procedure is illustrated in Fig. 6 where, in addition to the reflectance (solid line), we also display contributions of the individual modes (shaded areas) to the total spectrum. To demonstrate the effect of poles corresponding with all 6 previously described eigenmodes we extended the range of $k_{x}$ (B2’ refers to the pole B2 shifted by $-2\pi$). Clearly, B1 and B2 are responsible for peaks at band edges while L and qBIC for the peaks in the continuum. Even L’ and qBIC’ (or other poles such as B2’) may affect the spectrum in the physical range $0\leq k_{x}a\leq\pi$ provided that the overlap with their tails is sufficiently large.

The destructive interference mainly between L and qBIC modes leads to the asymmetric resonance associated with qBIC mode and minima $R=0$. This feature corresponds to the asymmetric dip and maxima $T=1$ in the transmittance, such as in Fig. 3 (b). In fact, it can be shown PK_pra that superposition of L and qBIC modes, playing the roles of the continuum and a narrow discrete state with poles $p_{c}$ and $p_{s}$, respectively, leads to the Fano formula

$$T(\Omega)=\left(\frac{1}{1+f^{2}}\right)\frac{\left(\Omega+f\right)^{2}}{\Omega^{2}+1},$$

(21)

where

$$\Omega\equiv\frac{q-\mathrm{Re}(p_{s})}{\mathrm{Im}(p_{s})},\;\;f\approx\frac{\mathrm{Im}(p_{c})}{\mathrm{Re}(p_{c})-\mathrm{Re}(p_{s})}.$$

(22)

In conclusion, we studied systematically spectral and scattering properties of a photonic analog of an extended Fano-Anderson tight-binding model. We classified the modes according to the relation between the self-coupling coefficients and the change of the propagation constant. The resonant features in the scattering spectra of the structures with broken symmetry can be interpreted through a generalized Weierstrass factorization. In particular, the Fano resonance associated with quasi-BIC arises from the interference between this state and another leaky mode. The transmittance near the quasi-BIC resonance can be rewritten into the form of the Fano formula, where the shape parameter $f$ can be expressed in terms of the poles of the interacting modes. Our work provides the theoretical framework which describes transformation of the symmetry protected BIC into a leaky mode and allows to interpret the resonant properties of the more complex systems. In addition through the engineering of zeros of the transmission matrix, it offers a possibility to investigate their nontrivial topological properties genack and may prove to be useful in various fields of optics as well as in cold matter and quantum confined systems.

We acknowledge financial support by the Czech Science Foundation (CSF) through the project 1900062S.