Multipolar origin of bound states in the continuum

TL;DR

This paper introduces a multipolar framework to understand and design bound states in the continuum within dielectric metasurfaces, enabling high-quality resonances for advanced flat-optics applications.

Contribution

It presents a novel multipole decomposition approach to explain the formation of bound states in the continuum in dielectric metasurfaces, considering symmetry conditions and spatial dispersion.

Findings

Bound states originate from specific multipolar modes.

Symmetry conditions determine the existence of bound states.

Metasurfaces can support bound states with wavevectors forming a line in reciprocal space.

Abstract

Metasurfaces based on resonant subwavelength photonic structures enable novel ways of wavefront control and light focusing, underpinning a new generation of flat-optics devices. Recently emerged all-dielectric metasurfaces exhibit high-quality resonances underpinned by the physics of bound states in the continuum that drives many interesting concepts in photonics. Here we suggest a novel approach to explain the physics of bound photonic states embedded into the radiation continuum. We study dielectric metasurfaces composed of planar periodic arrays of Mie-resonant nanoparticles ("meta-atoms") which support both symmetry protected and accidental bound states in the continuum and employ the multipole decomposition approach to reveal the physical mechanism of the formation of such nonradiating states in terms of multipolar modes generated by isolated meta-atoms. Based on the symmetry of the vector spherical harmonics, we identify the conditions for the existence of bound states in the continuum originating from the symmetries of both the lattice and the unit cell. Using this formalism we predict that metasurfaces with strongly suppressed spatial dispersion can support the bound states in the continuum with the wavevectors forming a line in the reciprocal space. Our results provide a new way of designing high-quality resonant photonic systems based on the physics of bound states in the continuum.

Figures

Equations

$${\mathbf{E}}^{in}(\mathbf{r})=\sum _{n=1}^{\infty }\sum _{{}_{p_i,p_r}^{m=0}}^n{E}_0[D_{p_i{p}_{r}mn}{\mathbf{W}}_{p_i{p}_{r}mn}^{(1)}(k_2,{\mathbf{r}}^{\prime })e^{i({\mathbf{k}}_b\cdot {\mathbf{r}}_j)}],$$

$$\mathbf{E}(\mathbf{r})=E_0\sum _{\mathbf{K},p_i,p_r,m,n}{\tilde{D}}_{p_i{p}_{r}mn}\frac{V_b{i}^{-n}}{2\pi k_1}\times$$

$$\times \underset{-\infty }{\overset{\infty }{\iint }}d{k}_{\parallel }\delta ({\mathbf{k}}_b-{\mathbf{k}}_{\parallel }-\mathbf{K})\frac{e^{i\mathbf{kr}}}{k_z}\left[{\mathbf{Y}}_{p_i{p}_{r}mn}\left(\frac{\mathbf{k}}{|\mathbf{k}|}\right)\right].$$

$$\mathbf{E}(\mathbf{r})=\frac{E_0{V}_b}{2\pi k_1{k}_{1z}}{e}^{i{\mathbf{k}}_1\mathbf{r}}\sum _{{}_{m,n}^{p_i,p_r,}}{i}^{-n}{\tilde{D}}_{p_i{p}_{r}mn}\left[{\mathbf{Y}}_{p_i{p}_{r}mn}\left(\frac{{\mathbf{k}}_1}{k_1}\right)\right],$$

$${\mathbf{M}}_{{}_1^{-1}mn}=\nabla \times (\mathbf{r}{\psi }_{{}_{-1}^{1}mn}),$$

$${\mathbf{N}}_{{}_{-1}^{1}mn}=\frac{\nabla \times {\mathbf{M}}_{{}_1^{-1}mn}}{k},$$

$${\psi }_{1mn}=\cos m\phi P_n^m(\cos \theta )z_n(kr),$$

$${\psi }_{-1mn}=\sin m\phi P_n^m(\cos \theta )z_n(kr)$$

$$\displaystyle{{\mathbf{\hat{G}}_{0}}}({\mathbf{r},\mathbf{r}^{\prime}})=\frac{ik_{1}}{4\pi}\sum_{n=1}^{\infty}\sum_{p_{r}=\pm 1}\sum_{m=0}^{n}(2-\delta_{0})\frac{2n+1}{n(n+1)}\frac{(n-m)!}{(n+m)!}\cdot\ \ \ \ \ \ \ \ \ \ \ \$$

$$\displaystyle\cdot\Bigl{(}[\mathbf{M}_{p_{r}mn}^{(3)}(k_{1},\mathbf{r})\otimes{\mathbf{M}}^{(1)}_{p_{r}mn}(k_{1},\mathbf{r}^{\prime})]+[{\mathbf{N}}_{p_{r}mn}^{(3)}(k_{1},\mathbf{r})\otimes{\mathbf{N}}^{(1)}_{p_{r}mn}(k_{1},\mathbf{r}^{\prime})]\Bigr{)},\ \ \ \ r>r^{\prime}$$

$$\mathbf{X}_{-p_{r}mn}\left(\frac{\mathbf{k}}{k}\right)={\nabla\times\left(\mathbf{k}Y_{p_{r}mn}\left(\frac{\mathbf{k}}{k}\right)\right)}\$$

$$\mathbf{Z}_{p_{r}mn}\left(\frac{\mathbf{k}}{k}\right)={i\frac{\mathbf{k}}{k}\times\mathbf{X}_{-p_{r}mn}\left(\frac{\mathbf{k}}{k}\right)}\$$

$$Y_{1mn}=\cos m\phi P_n^m(\cos \theta )$$

$$Y_{-1mn}=\sin m\phi P_n^m(\cos \theta )$$

$$\mathbf{E}(\mathbf{r})=\frac{k_1^2}{4\pi }\int d^3{\mathbf{r}}^{\prime \prime }\Delta \epsilon ({\mathbf{r}}^{\prime \prime }){\hat{\mathbf{G}}}_0(\mathbf{r},{\mathbf{r}}^{\prime \prime }){\mathbf{E}}^{in}({\mathbf{r}}^{\prime \prime })=$$

$$=\frac{k_1^2}{4\pi }({\epsilon }_2-{\epsilon }_1)\sum _j\underset{V}{\int }{d}^3{\mathbf{r}}^{\prime }{\hat{\mathbf{G}}}_0(\mathbf{r},{\mathbf{r}}_j+{\mathbf{r}}^{\prime }){\mathbf{E}}^{in}({\mathbf{r}}_j+{\mathbf{r}}^{\prime })$$

$$\mathbf{E}(\mathbf{r})=E_0\sum _j\frac{i{k}_1^3({\epsilon }_2-{\epsilon }_1)}{(4\pi {)}^2}\sum _{p_i,p_r,n,m}\sum _{p_i^{\prime },p_r^{\prime },n^{\prime },m^{\prime }}(2-{\delta }_0)\frac{2n+1}{n(n+1)}\frac{(n-m)!}{(n+m)!}{D}_{p_i^{\prime }{p}_r^{\prime }{n}^{\prime }{m}^{\prime }}{\mathbf{W}}_{p_i{p}_{r}mn}^{(3)}(k_1,\mathbf{r}-{\mathbf{r}}_j)\cdot$$

$$\cdot e^{i({\mathbf{k}}_b\cdot {\mathbf{r}}_j)}\underset{V}{\int }{d}^3{\mathbf{r}}^{\prime }[{\mathbf{W}}_{p_i{p}_{r}mn}^{(1)}(k_1,{\mathbf{r}}^{\prime })\cdot {\mathbf{W}}_{p_i^{\prime }{p}_r^{\prime }{m}^{\prime }{n}^{\prime }}^{(1)}(k_2,{\mathbf{r}}^{\prime })]$$

$$\mathbf{E}(\mathbf{r})=E_0\sum _{{}_j^{p_i,p_r,m,n}}{\tilde{D}}_{p_i{p}_{r}mn}{\mathbf{W}}_{p_i{p}_{r}mn}^{(3)}(k_1,\mathbf{r}-{\mathbf{r}}_j)\cdot e^{i({\mathbf{k}}_b\cdot {\mathbf{r}}_j)}.$$

$${\mathbf{M}}_{p_{r}mn}^{(3)}(k,\mathbf{r})=\frac{i^{-n}}{2\pi k}\underset{-\infty }{\overset{\infty }{\iint }}d{k}_{||}\frac{e^{i(k_{x}x+k_{y}y\pm k_{z}z)}}{k_z}\left[{\mathbf{X}}_{p_{r}mn}\left(\frac{\mathbf{k}}{k}\right)\right],$$

$${\mathbf{N}}_{p_{r}mn}^{(3)}(k,\mathbf{r})=\frac{i^{-n}}{2\pi k}\underset{-\infty }{\overset{\infty }{\iint }}d{k}_{||}\frac{e^{i(k_{x}x+k_{y}y\pm k_{z}z)}}{k_z}\left[{\mathbf{Z}}_{p_{r}mn}\left(\frac{\mathbf{k}}{k}\right)\right].$$

$${\mathbf{W}}_{p_i{p}_{r}mn}^{(3)}(k,\mathbf{r})=\frac{i^{-n}}{2\pi k}\underset{-\infty }{\overset{\infty }{\iint }}d{k}_{||}\frac{e^{i(k_{x}x+k_{y}y\pm k_{z}z)}}{k_z}\left[{\mathbf{Y}}_{p_i{p}_{r}mn}\left(\frac{\mathbf{k}}{k}\right)\right],$$

$$\mathbf{E}(\mathbf{r})=E_0\sum _{j,p_i,p_r,n,m}{\tilde{D}}_{p_i{p}_{r}mn}\cdot e^{i({\mathbf{k}}_b\cdot {\mathbf{r}}_j)}\cdot \frac{i^{-n}}{2\pi k_1}\underset{-\infty }{\overset{\infty }{\iint }}d{k}_{x}d{k}_y\frac{e^{i(\mathbf{k}(\mathbf{r}-{\mathbf{r}}_j))}}{k_z}\left[{\mathbf{Y}}_{p_i{p}_{r}mn}\left(\frac{\mathbf{k}}{k}\right)\right].$$

$$\sum _j{e}^{i\mathbf{k}{\mathbf{x}}_j}=V_b\sum _{\mathbf{K}}\delta (\mathbf{k}-\mathbf{K})$$

$${\tilde{D}}_{p_i{p}_{r}mn}=\frac{i{k}_1^3}{(4\pi {)}^2}(2-{\delta }_0)\frac{2n+1}{n(n+1)}\frac{(n-m)!}{(n+m)!}{D}_{p_i{p}_{r}mn}\sum _{p_i^{\prime },p_r^{\prime },n^{\prime },m^{\prime }}\underset{V}{\int }{d}^3{\mathbf{r}}^{\prime }({\epsilon }_2-{\epsilon }_1)[{\mathbf{W}}_{p_i{p}_{r}mn}^{(1)}(k_1,{\mathbf{r}}^{\prime })\cdot {\mathbf{W}}_{p_i^{\prime }{p}_r^{\prime }{m}^{\prime }{n}^{\prime }}^{(1)}(k_2,{\mathbf{r}}^{\prime })]$$

$${\tilde{D}}_{p_i{p}_{r}mn}^M=i{D}_{p_i{p}_{r}mn}^M{a}^2{\epsilon }_1[k_1^2{j}_{n-1}(k_{1}a)j_n(k_{2}a)-k_1{k}_2{j}_{n-1}(k_{2}a)j_n(k_{1}a)]$$

$${\tilde{D}}_{p_i{p}_{r}mn}^N=i{D}_{p_i{p}_{r}mn}^N{a}^2{\epsilon }_1(\frac{n+1}{2n+1}[k_1^2{j}_{n-2}(k_{1}a)j_{n-1}(k_{2}a)-k_1{k}_2{j}_{n-2}(k_{2}a)j_{n-1}(k_{1}a)]+$$

$$+\frac{n}{2n+1}[k_1^2{j}_n(k_{1}a)j_{n+1}(k_{2}a)-k_1{k}_2{j}_n(k_{2}a)j_{n+1}(k_{1}a)])$$

Full text

nanophotonics, multipolar decomposition, metasurface, photonic crystals, bound states in the continuum

I Introduction

The quest for compact photonic systems with high quality factor ($Q$ factor) modes led to the rapid development of optical bound states in the continuum (BICs). BICs are non-radiating states, characterized by the resonant frequencies embedded to the continuum spectrum of radiating modes of the surrounding space Hsu et al. (2016); Koshelev et al. (2018). The BICs first appeared as a mathematical curiosity in quantum mechanics Neumann and Wigner (1929). The discovery of BICs in optics immediately attracted broad attention (see, e.g., Refs. Paddon and Young (2000); Bulgakov and Sadreev (2008); Marinica et al. (2008)) due to high potential in applications in communications Dreisow et al. (2009); Gentry and Popović (2014), lasing Ha et al. (2018); Kodigala et al. (2017); Penzo et al. (2017); Rybin and Kivshar (2017), filtering Foley and Phillips (2015), and sensing Foley et al. (2014); Romano et al. (2018a, b). Recent achievements in the field of BIC are discussed in Refs. Carletti et al. (2018); Gao et al. (2016); Bulgakov and Sadreev (2014); Azzam et al. (2018); Jin et al. (2018); Chen et al. (2019); Abujetas et al. (2019); Kartashov et al. (2018); Romano et al. (2018c); Bulgakov and Sadreev (2017).

Decoupling of the resonant mode from the radiative spectrum, which is the basic idea behind the BIC, can be interpreted in several equivalent ways. Within the coupled-mode theory, it corresponds to nullifying the coupling coefficient between the resonant mode and all radiation channels of the surrounding space Wei Hsu et al. (2013). Alternatively, the appearance of BICs is explained as vanishing of Fourier coefficients corresponding to open diffraction channels due to the symmetry of the photonic structure. At the particular high-symmetry points of the reciprocal space, for example, at the $\Gamma$ point, the continuous spectrum is divided into the modes of different symmetry classes with respect to the reflectional and rotational symmetry of the photonic system. The bound states of one symmetry class can be found embedded in the continuum of another symmetry class, and their coupling is forbidden as long as the symmetry is preserved. Such kind of BIC is called symmetry-protected, and they also allow interpretation in terms of topological charges defined by the number of times the polarization vectors winds around the BICs presented as vortex centers in the polarization field Zhen et al. (2014). In contrast to the symmetry-protected BIC, the so-called accidental BICs Hsu et al. (2013); Zhen et al. (2014); Hsu et al. (2016); Yu et al. (2018); Ni et al. (2016) can be observed out of the $\Gamma$-point due to an accidental nulling of the Fourier (coupling) coefficients via fine tuning of parameters of the photonic system. Such a mechanism is also known as Friedrich-Wintgen scenario Friedrich and Wintgen (1985).

Despite the number of existing approaches to understand the nature of BICs, there is still a room for further development of the theory. During the previous few years the electromagnetic multipole theory Jackson (1999) has been extensively developed as a natural tool of nanophotonics dealing with the lowest (fundamental) resonances of the system. The main advantage of the multipole decomposition method (MDM) is that it provides a representation of an arbitrary field distribution as a superposition of the fields created by a set of multipoles Grahn et al. (2012); Kruk et al. (2017). Namely, the multipole expansion has been widely used to determine the polarization and directivity patterns of the scattered field of single particles and their clusters both plasmonic and dielectric Evlyukhin et al. (2010); Grahn et al. (2012) for a variety of applications such as polarization control device Kruk et al. (2017), dielectric nanoantenna Kuznetsov et al. (2016), light demultiplexing Panmai et al. (2018), and others. A number of novel optical phenomena have been explained within the MDM such as anapole effect Baryshnikova et al. (2018); Wu et al. (2018), optomechanical phenomena Hsu (2015); Poshakinskiy and Poddubny (2019), and Kerker effect Shamkhi et al. (2018); Liu and Kivshar (2018).

In this work, we extend the MDM approach for explaining both symmetry protected and accidental BICs. We provide a theory of BICs origin in terms of MDM for a general case of any periodic structure and develop an analytical method, which determines the contribution of the vector spherical harmonics to the far field (Section II). Working in the vector spherical harmonics (VSH) basis, we take the advantage of the internal symmetry and provide the group-symmetry approach to identifying the BIC formation in terms of the unit cell and lattice symmetries (Section III). We implement field multipole expansion of the eigenmodes of a periodic two-dimensional (2D) photonic structure supporting BIC. We illustrate the developed technique by considering a 2D square array of spheres and extending it to the case of a photonic crystal slab with a 2D array of cylindrical holes (Section IV). The developed approach can be easily extended even further to periodic structures with other types of the unit cell and other lattice symmetries. The proposed method both provides a deeper understanding of the photonic BIC physics and gives a tool for an effective designing of high-$Q$ resonant photonic systems.

II Multipolar approach

In this section we consider the modes of a two-dimensional periodic array of dielectric nanoparticles with arbitrary shape (see Fig. 1), and obtain an expression connecting the multipolar content of the field inside and outside the nanoparticles. We denote the VSHs as $\mathbf{M}_{p_{r}mn}(k,\mathbf{r})$ (magnetic) and $\mathbf{N}_{p_{r}mn}(k,\mathbf{r})$ (electric) relying on the definition presented in Ref. Bohren and Huffman, 1983. We introduce an additional notation $\mathbf{W}_{p_{i}p_{r}mn}$ for the both types of VSHs, where inversion parity index $p_{i}=(-1)^{n+1}$ for $\mathbf{M}_{p_{r}mn}(k,\mathbf{r})$, and $p_{i}=(-1)^{n}$ for $\mathbf{N}_{p_{r}mn}(k,\mathbf{r})$. Index $n$ is the multipole order, $m$ varies from 0 to $n$, and $p_{r}=1$ if $\mathbf{W}_{p_{i}p_{r}mn}$ is even under reflection from $y=0$ plane ($\varphi\rightarrow-\varphi$ in the spherical system), and $p_{r}=-1$ if it is odd Frizyuk (2018).

In homogeneous medium with permittivity $\varepsilon$, any solution of the Helmholtz equation with wavevector $k=\sqrt{\varepsilon}\frac{\omega}{c}$ can be expanded in terms of electric and magnetic spherical harmonics Stratton (2007). The field inside the medium of $j$-th cell of the array $\mathbf{E}^{in}$ can be written in terms of multipolar decomposition:

[TABLE]

where $k_{2}=\sqrt{\varepsilon_{2}}\frac{\omega}{c}$ is the wavevector in the material, the superscript $(1)$ stands for spherical Bessel functions in the radial part of the VSH, $\mathbf{k}_{b}$ is the Bloch vector, $\mathbf{r}^{\prime}=\mathbf{r}-\mathbf{r}_{j}$, and $\mathbf{r}_{j}$ is the position of a single sphere (see Fig. 1). $D_{p_{i}p_{r}mn}$ are the coefficients of the multipolar decomposition, which will be discussed in Section III. In similar manner, the field outside the array is expressed as follows (see Appendix B):

[TABLE]

Here $\mathbf{K}$ is the reciprocal lattice vector, $k_{z}\leavevmode\nobreak\ =\pm\sqrt{k_{1}^{2}-k_{\parallel}^{2}}$, $|\mathbf{k}|=k_{1}=\sqrt{\varepsilon_{1}}\frac{\omega}{c}$, $V_{b}$ is the volume of the first Brillouin zone. The spherical vector functions $\mathbf{Y}(\mathbf{n})$ depend on the spherical coordinates of a unity vector $\mathbf{n}$ and they are given in the Appendix A. Their symmetry coincides with the symmetry of $\mathbf{W}$-function with the identical indices. The relation between $D$ and $\tilde{D}$ is discussed in Appendix C, and for the case of an array of spherical particles, coefficients $\tilde{D}$ can be derived analytically.

The summation in Eq. (II) over the reciprocal vectors corresponding to open diffraction channels, i.e. the terms with real $k_{z}$, provides the contribution into the far field. For the frequencies below the diffraction limit, only the zero-order term with $\mathbf{K}=0$ gives non-zero contribution to the far field:

[TABLE]

where $\mathbf{k}_{1||}=\mathbf{k}_{b},$ and $k_{1z}=\pm\sqrt{k_{1}^{2}-k_{b}^{2}}$. According to Eq. (3), the contribution of the multipole with numbers $p_{i},p_{r},m,n$ into the far field in the direction defined by the wave vector $\mathbf{k}_{1}$ is proportional to the multipole expansion coefficient $\tilde{D}$ and the value of spherical vector function $\mathbf{Y}$ in the given direction. Equation (3) provides the correspondence between the radiation pattern of a single unit cell and the far-field properties of the whole infinite array allowing for interpreting the BIC in terms of MDM. In strong contrast to a single nanoparticle, where each multipole contributes to the far field, in case of a subdiffractive array there might be direction, where non of the multipole gives any contribution, or alternatively the non-zero contribution of different terms may eventually sum up to zero. The formulated alternative gives a sharp distinction between the symmetry protected and accidental BIC.

At-$\Gamma$-point BIC. The $\Gamma$-point BIC corresponds to the absence of the far-field radiation in the direction along the $z$-axis. Due to the structure of VSH, it appears that a number of multipoles do not radiate in the vertical direction along the $z$-axis. If the field inside a single unit cell consists only of such multipoles, there will be no total radiation in $z$-direction. This simple fact is illustrated in the upper panel in Fig. 1(b). Noticing that only $\mathbf{Y}_{p_{i}p_{r}1n}$ functions with m=1 are non-zero in parallel to the $z-$axis direction, we can conclude that at the $\Gamma$-point in the subdiffractive array all the modes which do not contain the harmonics with $m=1$ are symmetry-protected BICs. This fundamental conclusion lies in the basis of recent experimental demonstration Ha et al. (2018) of lasing with BIC in a 2D subdiffractive array of nanoparticles. The particular operational mode consisted of vertical dipoles oriented along the $z$-axis, thus, not contributing into the only open channel. There exist an approach Zhen et al. (2014); Hsu et al. (2016) that the eigenmodes at the $\Gamma$-point can radiate in the normal direction $z$ if their fields are odd under $C^{z}_{2}$ rotations, and do not have any other rotational symmetry of $C^{z}_{n}$-type. In terms of multipole moments, this follows from the fact that at the $\Gamma$-point any radiative mode should contain multipoles with $m=1$. On the other hand, in virtue of the symmetry, the even modes have zero radiation losses, i.e. infinite radiation quality factor, which are known as symmetry-protected BICs.

Off-$\Gamma$ BIC. Let us now turn to the case of accidental BIC description. In general case, coefficients $D$ are complex numbers. They define the amplitudes and phase delay between the multipoles. However, in accordance with Hsu et al. (2013), if the structure has time reversal and inversion symmetry, the eigenmodes must satisfy the condition $\mathbf{E}(\mathbf{r})=\mathbf{E}^{*}(\mathbf{-r})$. This fact imposes strict conditions on the multipoles’ phases, because some of them are even under inversion and some of them are odd, the coefficient $D_{-1p_{r}mn}$ before the odd ones must be imaginary. It follows that every term of the sum in Eq. (3) is purely real, and all multipoles are in-phase or anti-phase. All coefficients depend on k-vector and structure’s parameters and in a case of the off-$\Gamma$ BICs this sum turns to zero. In other words, for the particular $\mathbf{k}_{1}$ all vector harmonics add up to zero in the direction of $\mathbf{k}_{1}$, analogously to the anti-Kerker effect, because they are already in phase and only amplitudes are modulated while k-vector is changed [Fig. 1(b), lower panel]. One can intuitively understand accidental BIC (TM-polirized) via a toy dipole model composed of vertical electric and horizontal magnetic dipoles interfering destructively at some k-vector Doeleman et al. (2018).

The expansion coefficients $D_{p_{i}p_{r}mn}$ depend on shape of the nanoparticles, material parameters and symmetry of the lattice. Obviously, that the lattice symmetry should impose restrictions on $D_{p_{i}p_{r}mn}$ and some coefficients should vanish due to the symmetry. To explore these selection rules, we employ theoretical-group approach.

III Symmetry approach

The group-theoretic approach is a powerful method which is widely used for analyzing the properties of periodic photonic system Hergert and Dane (2003); Ohtaka and Tanabe (1996a, b); Sakoda (2001); Yu et al. (2018). In this section, we apply this method to reveal which of the coefficients $D_{p_{i}p_{r}mn}$ are non-zero in accordance with the mode symmetry imposed by the symmetry of the periodic structure Hayami et al. (2018); Gelessus et al. (1995) and to explain the formation of BIC. To make further analysis more illustrative we will provide it by the example of a square periodic array of dielectric spheres shown in the inset in Fig. 2(a). Nevertheless, it is necessary to highlight that further considerations remain true for any kind of structures with square lattice with the point group symmetry $D_{4h}$.

Figure 2(a) shows the dispersion of eigenmodes along the $\overline{\Gamma X}$ direction for a square periodic array of dielectric spheres with permittivity ${\varepsilon}_{2}=12$ embedded in the air with permittivity ${\varepsilon}_{1}=1$. The dispersion is calculated numerically using COMSOL Multiphysics package. The radius of the spheres is $a=100$ nm and the period of the array is $d=600$ nm. By the analogy with a dielectric slab waveguide, the eigenmodes of the 2D array are split into transverse electric (TE) and transverse magnetic modes (TM), which have mainly the $z$-component of magnetic and electric field, correspondingly. Figure 2(b) shows the dependence of the Q-factor on the Bloch wavenumber $k_{x}$ for the modes, which are BICs at the $\Gamma$-point. One can see that additionally to at-$\Gamma$-BIC, TM3 mode turns into off-BIC in the middle of the Brillouin band (accidental BIC).

The eigenmodes of periodic structures have certain symmetry in accordance with the fact that every mode is transformed by an irreducible representation of the structure’s symmetry group Agranovich and Ginzburg (2013); Ivchenko and Pikus (1995). While Bloch functions are already the basis functions of the translation group irreducible representation labeled by $\mathbf{k}_{b}$, under point group operations the basis functions with different $\mathbf{k}_{b}$ transform through each other. They carry high-dimensional representations of the space group. We are interested in the multipolar content of the mode with particular $\mathbf{k}_{b}$, so we look at the point group of $\mathbf{k}_{b}$, i.e. the subgroup of the whole point group, which transforms the $\mathbf{k}_{b}$ into an equivalent. Note that symmetry of the unit cell should be also taken into account since it alters the point group of the full structure.

At the $\Gamma$-point, the group of $\mathbf{k}_{b}$ is the whole group $D_{4h}$, $\mathbf{k}_{b}=0$ and it is not transformed. At $\overline{\Gamma X}$-valleys the point group of $\mathbf{k}_{b}$ is $C_{2v}$, which consists of $\pi$-rotation around the $x$- or $y$- axis and two plane reflections at $z=0$ and at $y=0$ or $x=0$. These operations keep vector $\mathbf{k}_{b}$ invariant. Analogously, in the $\overline{\Gamma M}$-valley the group is also $C_{2v}$. Solutions with particular $\mathbf{k}_{b}$ are transformed by one of the $\mathbf{k}_{b}$-group irreducible representations. Thus, since the solution is transformed as a basis function of some particular representation, the multipolar content is strictly limited. Namely, all the multipoles with non-zero contribution must be transformed by the similar irreducible representation of the $\mathbf{k}_{b}$-group. We will use common notations for irreducible representations, listed, for example, in Ref. Goss, .

$\Gamma$-point. For example, we consider the TE1 mode of the square array, which is transformed by $A_{2g}$ at the $\Gamma$-point. Under the transformations of $D_{4h}$ group the only low-order multipoles transformed by $A_{2g}$ are magnetic dipole $\mathbf{M}_{-101}$, magnetic octupole $\mathbf{M}_{-103}$, and electric hexadecapole $\mathbf{N}_{-144}$. All of them are invariant under $C_{4}$ rotations, even under inversion and $z=0$-plane reflection and odd under other $D_{4h}$ transformations. Higher-order multipoles which behave in the same way are also presented in the multipolar content of this mode.

Analogously, we classify all possible multipoles at the $\Gamma$-point in accordance with their symmetry and provide the tables with multipolar content of the modes (Tables 1 and  2). Note, that TE2 and TE3 modes degenerate, and they transform through each other as two basis functions of the representation $E_{u}$.

$\overline{\Gamma X}$* valley*. After the symmetry reduction, when we step out the $\Gamma$-point into $\overline{\Gamma X}$ valley, some symmetry operations remain, e.g. mirror reflections in $z=0$ and $y=0$ planes and rotation by $\pi$ around the $x$-axis. Eigenmodes must behave in the same way under these symmetry operations, as at the $\Gamma$-point.

As an example, we consider TE1 mode, which transforms by $A_{2g}$ at the $\Gamma$-point, and TM3 mode ($A_{1u}$) at the $\overline{\Gamma X}$ valley. Using the compatibility relations Sakoda (2001), we obtain that the mode which transforms by $A_{2g}$ at the $\Gamma$-point is transformed under $B_{2}$ representation of the group $C_{2v}$ in the $\overline{\Gamma X}$ valley. The TM3-mode, which is $A_{1u}$ at the $\Gamma$, is transformed under $B_{1}$ in $\overline{\Gamma X}$. For the $A_{2g}$ mode at the $\Gamma$-point, which is odd under reflection in the $y=0$ plane and $\pi$-rotation around $x$ and even under reflection in $z=0$ plane, the only possible multipoles in $\overline{\Gamma X}$ valley should have the same symmetry properties. For the $A_{1u}$ mode the possible multipoles in the $\overline{\Gamma X}$ valley must be odd under reflection in the $z=0$ plane and $\pi$-rotation and even under reflection in the $x=0$ plane. Low-order possible multipoles in the $\overline{\Gamma X}$ valley are listed at the right column of Table 1 for TE modes, and for TM modes are easily derived by replacing $A_{2}\leftrightarrow A_{1},B_{2}\leftrightarrow B_{1}$, $u\leftrightarrow g,$  $p_{r}\leftrightarrow-p_{r},\ \mathbf{M}\leftrightarrow\mathbf{N}$. Analogously, for the $\overline{\Gamma M}$ valley we have specific symmetry in the $x=y$ or $x=-y$ direction and possible multipoles are the same as those which are transformed by the same representation in $\overline{\Gamma X}$ valley, but rotated by $\pi/4$ with help of Wigner D-matrixes Huayong and Yiping (2008); Aubert (2013).

IV Multipolar composition of the eigenmodes in periodic structures

IV.1 Multipole analysis of metasurfaces and photonic crystal slabs

Before applying the developed approach to certain structures, we would like to emphasize that the power of the group-theoretic approach and multipole decomposition method is that these methods applied equally for metasurfaces (arrays of meta-atoms) and photonics crystal slabs (slabs with holes). For example, a square array of dielectric spheres and dielectric slab with circular holes arranged in a square lattice are identical from the point of view of the group theory. Therefore, these structures have the same classification of eigenmodes and, that more important, the same multipole content (Mie-resonances). The difference will be only in amplitudes of the multipole coefficients.

At first glance it may seem unphysical to characterize the unit cell with a hole by Mie resonances. This contradiction can be eliminated by means of the proper choice of the unit cell. Indeed, Fig. 3 shows the modes of the photonic-crystal slab and metasurface corresponding to the same irreducible representation $B_{\text{1u}}$ (see Table 2). One can see that the mode structures are completely identical. However it is more physical to assign Mie resonances not the unit cell with a hole but to the cross-like shape unit cell (see the low panel in Fig. 3). In particular, it was shown in Yang et al. (2019) that the optical response of both photonic-crystal slab and metasurface is dominated by the electric and magnetic Mie-type dipole resonances, and both type of the structures provides $2\pi$ phase control of light.

IV.2 Periodic arrays of dielectric spheres

In order to understand how the interference of the multipole moments forms the BICs, we applied a method of multipole decomposition by expanding the eigenmodes fields inside the nanoparticles in terms of spherical harmonics Jackson (1999); Grahn et al. (2012). Figure 4 shows the numerical results for the multipole decomposition for TE1, TE4, and TM4 modes at the $\Gamma$-point and in a point of the $\Gamma$X-valley. At the high-symmetry $\Gamma$-point only a small fraction of all possible multipoles is presented, while after lowering the symmetry extra multipoles appear out of the $\Gamma$-point. Figure 5 shows the numerical results for the multipole decomposition for TM3 modes at the $\Gamma$-point, in a point of the $\Gamma$X-valley, and at the off-$\Gamma$ BIC.

At-$\Gamma$-BIC. One can see, that for at-$\Gamma$ TE${}_{1}(A_{2g}$) BIC-mode magnetic dipole along the $z$-axis ($\mathbf{M}_{-101}$) is mainly contributed [Fig. 4(a)]. Its directivity pattern restricts radiation along the $z$-axis, and the radiation patterns of all remaining multipoles at the $\Gamma$-point are also zero along the $z$-axis. Other directions of radiation are forbidden due to the subdiffraction regime. Similarly, at-$\Gamma$ TM3($A_{1u}$) BIC is formed by electric dipole along the $z$-axis ($\mathbf{N}_{101}$) mostly, prohibiting the radiation in vertical direction itself [Fig. 5(a)]. The TE4 ($B_{2g}$) and TM4 ($B_{1u}$) modes’ lowest multipoles are the electric quadrupole $\mathbf{N}_{-122}$ and the magnetic quadrupole $\mathbf{M}_{122}$ respectively. They have $m=2$, which also prohibits the radiation [Figs. 4(c) and 4(e)]. In contrast to BICs, the radiative modes ($E_{g},E_{u}$) degenerate at the $\Gamma$ point since they are transformed by the two-dimensional representations. From the symmetry-group approach, we know that TE2,3 and TM1,2 modes contain electric $\mathbf{N}_{\pm 111}$ and magnetic $\mathbf{M}_{\pm 111}$ spherical harmonics. The numerical multipole expansion shows that degenerated modes contain in-plane electric or magnetic dipole moments as the main contribution. These numerical results validate the symmetry-group approach for the system with square lattice, confirming that any symmetry-protected BIC is characterized by multipole moments with $m\neq 1$.

Off-$\Gamma$ BIC. Away from the $\Gamma$-point other multipoles appear in decomposition [Figs. 4(b), 4(d), 4(f) and  5(b)] and BIC is destroyed turning into a resonance state. As mentioned earlier in the Section II, the accidental off-$\Gamma$ BICs is formed due to either in- or anti-phase contributions of different multipoles. We obtain the same result for TM3 mode with off-$\Gamma$ BIC in the $\Gamma$X-valley [Fig. 5(c)]. This mode is transformed by $B_{1}$ representation and consist of the multipoles, which are odd under reflection in the $z=0$ plane and $\pi$-rotation around $x$-axis and even under $y=0$ reflection, eg. $i\mathbf{N}_{101},\mathbf{N}_{112},i\mathbf{N}_{103},i\mathbf{N}_{123},\mathbf{M}_{111},i\mathbf{M}_{122},\mathbf{M}_{113},\mathbf{M}_{133}$. They sum up into $TM$-polarized wave in the direction given by vector $\mathbf{k}_{1}$ and cancel each other in the off-$\Gamma$ point forming accidental BIC, shown in Fig. 2(b). In addition to all, it is well known, that the $z=0$ plane reflection symmetry of the structure is required to obtain the off-$\Gamma$ BIC Hsu et al. (2013). Indeed, in lack of such symmetry, each mode would contain both odd and even multipoles under reflection in the $z=0$ plane. To restrict the radiation both in the upper- and lower half-spaces, odd and even multipoles should be summed up into zero independently, while for the symmetric structure only one type of multipoles is presented for each mode, which makes it possible to achieve the BIC by tuning the structure parameters.

IV.3 Photonic crystal slab

We extend our approach beyond the 2D array of spheres and apply it to the photonic crystal slab of the same symmetry. We consider a dielectric slab with a square array of cylindrical holes studied in Hsu et al. (2013). Both at-$\Gamma$ (symmetry-protected) and off-$\Gamma$ (accidental) BICs appear in the lowest TM band referred to as TMh. This mode has the field profile of the same symmetry as the TM4 mode of the array of the spheres and is transformed by $B_{1u}$ representation.

The description of the far field defined by the Eq. (3) is still valid, but the coefficients $\tilde{D}_{p_{i}p_{r}mn}$ in general cannot be expressed analytically through $D_{p_{i}p_{r}mn}$ (see Appendix C). However, only the multipoles which are presented in the field inside the slab contribute into far field, $\tilde{D}_{p_{i}p_{r}mn}=0$ if $D_{p_{i}p_{r}mn}=0$. The multipolar content for the considered slab is the same as in the periodic array of spheres because the modes of the slab have the same symmetry as the modes of the array. The multipole decomposition of TMh mode [Fig. 6(a)] reveals that a magnetic quadrupole $\mathbf{M}_{122}$ and electric octupole $\mathbf{N}_{123}$ make the major contribution to the at-$\Gamma$ BIC, as well as to the TM4 mode of the array of the spheres. However, the $B_{1u}$ mode of the photonic crystal slab is the lowest-energy TM mode while for the array of spheres it has the highest energy among modes under diffraction limit. Due to the variational principle Joannopoulos et al. (2008), for the mode of such symmetry, the electric field is more concentrated inside the high-index material in case of photonic crystal slab, minimizing the energy of the mode. Although dispersion curves $\omega(k)$ of the modes TM4 and TMh behave completely differently, multipole decomposition proves a common origin of them. At the $\Gamma$-point, we obtain contributions of multipoles only with $m=2,6,10,$ etc. for both modes, and none of these multipoles contribute in the far-field. At the $\overline{\Gamma X}$-valley, the multipolar content of the TM4 and TMh is similar [Figs. 4(e), 4(f), 6(a), and 6(b)]. However, for the TMh it is possible to obtain an accidental off-$\Gamma$ BIC in the $\overline{\Gamma X}$-valley [Fig. 6(c)]. Away from the $\Gamma$-point, multipole contributions change smoothly keeping the multipolar content invariable, and at a particular wave vector $k_{x}$ multipoles interfere destructively forming the accidental BIC.

V Summary and outlook

Importantly, our approach based on the multipole decomposition analysis of individual meta-atoms not only explains clearly and in simple physical terms the origin of both symmetry protected and accidental bound states in the continuum, but it has also a prediction power and may be employed for both prediction and engineering different types of BICs. As an example, we consider a metasurface consisting of meta-atoms packed in a subwavelength 2D lattice, which are polarized purely as octupoles, for example, $\mathbf{M}_{-103}$ (see Fig. 7). Each octupole of this type has a nodal cone and, therefore, we can expect that in a periodic subwavelength array of such meta-atoms, BICs form a line in the reciprocal space. However, to observe this phenomenon, the effective polarizability of the unit cell accounting for the interaction between all meta-atoms should not depend on the Bloch wavenumber or have very weak dependence. In other words, the line of BIC could be observed in metasurfaces with suppressed spatial dispersion that is still a challenging problem. Interestingly, such a kind of BIC will be observed for the same directions independently on the lattice symmetry of the metasurface.

We can expand this discussion even further, by adding $\mathbf{M}_{-101}$ dipole mode into consideration. To obtain the Q-factor of the at-$\Gamma$ BIC, we should find the dependence of energy loss rate on $\mathbf{k}-$vector in Eq. (3), assuming stored energy is almost constant when $\theta_{k}\rightarrow 0$. The asymptotic behavior for $\theta_{k}\rightarrow 0$ of $\mathbf{Y}_{p_{i}p_{r}mn}\left(\frac{\mathbf{k}}{|\mathbf{k}|}\right)$ functions is defined by $\frac{dP_{n}^{m}(\cos\theta_{k})}{d\theta_{k}}$. So, $Y_{1-103}\left(\frac{\mathbf{k}}{|\mathbf{k}|}\right)\propto(-6k_{x}+17k_{x}^{3}/2)\mathbf{e}_{\varphi}$, and $Y_{1-101}\left(\frac{\mathbf{k}}{|\mathbf{k}|}\right)\propto(k_{x}-k_{x}^{3}/6)\mathbf{e}_{\varphi}$. We can manage their relative contribution, and at particular point when coefficient before dipole is six times larger, we obtain that linear terms cancel each other and field asymptotic is proportional to $k_{x}^{3}$, so the quality-factor growth is proportional to $k_{x}^{-6}$. Similar effect was observed in Jin et al. (2018) for the photonic crystal slab. However, considering realistic situations, we should take into account all possible multipolar contibutions, including terms with $m=1$, where the multipole contribution growth rate plays role but not the asymptotic behavior. Thus, the multipole origin of BIC makes a new query for metasurfaces with a suppressed spatial dispersion.

In summary, we have demonstrated that symmetry-protected bound states in the continuum in dielectric metasurfaces and photonic crystal slabs at the frequencies below the diffraction limit are associated with the multipole moments of the elementary meta-atoms which do not radiate in the transverse direction. For any type of metasurfaces, the symmetry-protected bound states in the continuum can be observed only if there exist no multipoles with the azimuthal index $m=1$ in the multipole decomposition. The symmetry approach allows to determine which multipolar content the lattice eigenmodes have, and it can be analogously applied to the structures with different symmetries, for example, hexagonal lattices or arrays of nanoparticles of an arbitrary shape and in-plane broken symmetry. Similarly, we have revealed that the accidental bound state in the continuum corresponding to an off-$\Gamma$ point in the reciprocal space is formed due to destructive interference of the multipole fields in the far-zone. We have provided the general tools for the analysis of bound states in the continuum based on the irreducible representation of the appropriate photonic band. We believe that our results will provide a new way for designing high-quality resonant photonic systems based on the physics of bound states in the continuum.

Acknowledgements.

The authors acknowledge useful discussions with I.D. Avdeev. Y.K. acknowledges a support from the Strategic Fund of the Australian National University. Z. S. acknowledges support from the Foundation for the Advancement of Theoretical Physics and Mathematics ”BASIS” (Russia). Z.S. and K.F. contributed equally to this work.

Appendix A Basic definitions

Vector spherical harmonics are defined as Bohren and Huffman (1983)

[TABLE]

where

[TABLE]

$z_{n}(kr)$ can be replaced by spherical bessel function of any kind, and $P_{n}^{m}(\cos\theta)$ are associated Legendre polynomials.

Vacuum dyadic Green’s function expansion in terms of vector spherical harmonics Tai (1994); Li et al. (1994):

[TABLE]

where superscripts $(1)$ and $(3)$ appear, when we replace $z_{n}(\rho)$ by spherical Bessel functions, and the spherical Hankel functions of the first kind, respectively, $\delta_{0}=1$ when $m=0$, and $\delta_{0}=0$ when $m\neq 0$.

Spherical vectors $\mathbf{Y}_{p_{i}p_{r}mn}$ denote two types of functions, $\mathbf{X}_{p_{r}mn}$ and $\mathbf{Z}_{p_{r}mn}$, defined as Alaee et al. (2019)

[TABLE]

[TABLE]

where

[TABLE]

$p_{i}$ is $(-1)^{n+1}$ for $\mathbf{X}$, and $p_{i}=(-1)^{n}$ for $\mathbf{Z}$. Note that the transformation behavior is similar for $\mathbf{W}$ and $\mathbf{Y}$, $\mathbf{X}$ and $\mathbf{M}$, $\mathbf{Z}$ and $\mathbf{N}$, $\psi$ and $Y$, respectively.

Appendix B Lattice sums of the spherical harmonics

We assume that multipolar content of the mode is already known, and coefficients in the formula (1) are given. With help of vacuum dyadic Green’s function $\hat{\mathbf{G}}_{0}$, we express the field outside the array

[TABLE]

where $k_{1}=\sqrt{\varepsilon_{1}}\frac{\omega}{c}$ is vacuum wavevector, $V$ is the single nanoparticle’s volume.

Green’s function can be also expressed in terms of vector spherical harmonics (see Appendix A). Using the property $\hat{\bf G}_{0}({\mathbf{r},\mathbf{r}_{j}+\mathbf{r}^{\prime},\omega})=\hat{\bf G}_{0}({\mathbf{r}-\mathbf{r}_{j},\mathbf{r}^{\prime},\omega})$, and substituting (1) into (B), we obtain

[TABLE]

here superscript $(3)$ stands for the outgoing spherical Hankel wave. Now we dwell on the case when the array is composed of spherical nanoparticles. Exploiting the VSH orthogonality properties, the integral over sphere can be taken analytically Stratton (2007). It is proportional to Kroneker delta $\delta_{p_{i}p^{\prime}_{i}}\delta_{p_{r}p^{\prime}_{r}}\delta_{mm^{\prime}}\delta_{nn^{\prime}}$ which removes one summation. Combining all coefficients including the intergral into $\tilde{D}$ we obtain that the field outside the array is expressed with the formula:

[TABLE]

Note, that coefficient $\tilde{D}_{p_{i}p_{r}mn}$ is non-zero only if the harmonic with such numbers is presented in the field expansion inside the sphere. If the fields created by each unit cell are already known, we can also start the considerations from this formula.

To obtain formula (II) we exploit the Weyl identity for VSH expansion through the plane waves in case when $z_{n}(kr)$ is replaced by spherical hankel functionsWittmann (1988); Stout (2012):

[TABLE]

[TABLE]

Sign before $k_{z}$ is depends on z- sign. Redefining the harmonics, we have

[TABLE]

where $k_{z}=\sqrt{k^{2}-k_{x}^{2}-k_{y}^{2}}$. Substituting this expansion into (15), we get

[TABLE]

This expansion helps us to apply the summation formula

[TABLE]

where $\mathbf{K}$ is reciprocal lattice vector and $V_{b}$ is the volume of the Brillouin zone ($\frac{(2\pi)^{2}}{d^{2}}$ for square lattice), and substituting (20) into (19), we obtain (II).

Appendix C Relation between the coefficients $D$ and $\tilde{D}$

The coefficients $D$ in (B) and $\tilde{D}$ (15) are connected by the formula

[TABLE]

This formula describes both array of nanoparticles, and any photonic crystal slab, but in case of the array of spheres, the integral can be easily taken analytically.

Remark Here we give the dyadic Green’s function ${\mathbf{\hat{G}_{0}}}(\mathbf{r},\mathbf{r}^{\prime})$ only for the case when $r>r^{\prime}$. In case of photonic crystal slab or non-spherical particles, we also need the part of Green’s function when $r^{\prime}<r$ Li et al. (1994) to obtain the near-field. This will alter the intermediate calculations, since we have to compute the lattice sum for VSHs with spherical bessel functions. Nevertheless, the answer will have the same form.

We apply the orthogonality properties of vector spherical harmonics Stratton (2007) [p.418], and consider integrals of magnetic and electric harmonics separately. Implementing the angular integration, we reduce the integral to the integral of r-dependent Bessel functions products, which also can be computed analytically. For magnetic harmonics we have

[TABLE]

where $a$ is nanoparticle radius, and for electric

[TABLE]

Note that this expression turns to zero at some frequencies, so we can have zero $\tilde{D}$ when $D$ is non-zero. This refers to the anapoles of the spherical nanoparticles. The frequency where the anapole appears is the same as for single isolated nanoparticle.

If we have other type of surface, for example, photonic crystal slab with holes, or array of the cylinders, the orthogonality property can’t be applied and the integral $\int\limits_{V}d^{3}\mathbf{r}^{\prime}\Delta\varepsilon[\mathbf{W}^{(1)}_{p_{i}p_{r}mn}(k_{1},\mathbf{r}^{\prime})\cdot\mathbf{W}^{(1)}_{p_{i}^{\prime}p_{r}^{\prime}m^{\prime}n^{\prime}}(k_{2},\mathbf{r}^{\prime})]$ will mix some harmonics. However, all the harmonics, which admix, are already presented in the expansion of the field inside the cell. This will just alter the coefficients before the outgoing multipoles, but not the multipolar content.