Terahertz bound state in the continuum in dielectric membrane metasurfaces

Mie-resonant metasurfaces composed of subwavelength dielectric resonators enable an efficient route for electromagnetic wave manipulation. Among these manipulations, a localized mode with a high-quality factor coexisting with a continuous spectrum of radiating waves termed bound state in the continuum (BIC) can arouse many exotic applications in photonics. Here, we demonstrate the terahertz BIC in a dielectric membrane metasurface and analyze its resonant nature based on Mie-resonant multipoles and vector spherical harmonics. The intrinsic splitting of the resonances under oblique incidence is also explored, in which the conversion of multipole radiation patterns versus the oblique angle will drive the resonances from BIC to leaky modes or vice versa. Both Γ and off-Γ point BICs could be identified as the superposition cancellation of vector spherical harmonics for both p-wave and s-wave. Our research not only provides a novel perspective for exploring the essence of BIC metasurfaces in the terahertz regime, but also points new opportunities for achieving terahertz BIC metasurfaces with ultra-high quality factors.


Introduction
The confinement or localization of light is critical for the control of phase, amplitude, and polarization [1,2], which has promising applications in chemical/biological sensing [3,4] and lasing process [5]. The bound state in the continuum (BIC), a sort of non-radiating states locating in the continuum of a system, has been proved to be an efficient way to achieve light confinement. BIC is a general wave phenomenon and has been identified in electromagnetic waves [6,7]. The typical characteristic of BIC is its unparalleled light confinement capability and theoretically infinite quality factor, which attract many theoretical and experimental investigations in lasing, nonlinear generation and sensing [3,5,[8][9][10][11]. Several concepts and physical mechanisms have been demonstrated to understand the essence of the BIC. One is the symmetry-protected BIC, which comes from the symmetry prohibited coupling to the leaky modes [12][13][14][15]. Another is the parameter-tuned BIC, where destructive interference between two or more radiating modes cancel each other, and this mechanism is known as the Friedrich-Wintgen scenario [7,[16][17][18]. A true BIC is a mathematical object with an infinite Q factor and exists only in ideal lossless structures. In practice, quasi-BICs arising from metasurfaces with broken in-plane symmetry have been demonstrated, and they also support a sharp high-Q radiative resonance [19][20][21][22].
To realize BIC or quasi-BIC devices, the intrinsic material loss is a crucial topic and need to be suppressed as low as possible. Due to the low intrinsic materials loss, dielectric metasurfaces enable a vast variety of efficient flat functional devices for optics [23][24][25][26][27][28]. Especially for high-index dielectric metasurfaces, electric and magnetic Mie-resonances dominate the interference and scattering of the incident light. By tuning the Mie multipoles, dielectric metasurfaces can accomplish complex light shaping and manipulation including Huygens' surfaces and wave plates [29][30][31]. Most dielectric metasurfaces need a low-index host or substrate to sustain the structures. The index differences from the host will affect the resonance and overall efficiency of the resonators. It has been demonstrated that this limitation can be avoided in a high index dielectric hole array formed in a free-standing membrane metasurface [32,33], which resembles a photonic-crystal slab, except that the thickness is deep subwavelength and its functionalities come from Mie-resonant multipoles [34][35][36].
Both the symmetry-protected and parameter-tuned BIC can be supported in the dielectric slab that is governed by the Bragg scattering from a periodically varying refractive index [36]. Compared with the dielectric slab having multiple waveguide modes, the benefit offered by a free-standing membrane also lies in the fact that the membrane metasurface with a sub-wavelength thickness has fewer electromagnetic modes, and the scattering and resonance could be interpreted via the Mie theory [32]. Based on the multipoles and vector spherical harmonics analysis, the existence of BICs in dielectric Mie-resonant metasurfaces or isolated meta-atoms originating from the symmetries of both the lattice and the unit cell has been demonstrated [37][38][39][40]. Particularly, for the terahertz community, BICs show a massive potential for communication, nonlinear generation, filtering and sensing. However, analysis of the essence of BICs/quasi-BICs in terahertz metasurfaces from the perspective of the Mie-resonant multipoles is still blank. Compared to the visible light domain, the strategies for realizing terahertz BIC metasurfaces are few and poor. Here, we investigate the Γ point and off-Γ point BICs in a dielectric membrane metasurface and explore their resonant nature based on Mie-resonant multipoles and vector spherical harmonics analysis. The intrinsic nature of the resonances under the p-wave and s-wave is also analyzed. The constructive or destructive interference of the multipole radiation patterns determines the far field behavior of the resonances. With an increment of the oblique angle, the multipole radiation pattern controls the BICs for both the modes under the p-wave and s-wave, which drives the resonances from BICs to leaky modes or vice versa. Our research presented here provides a deeper understanding of BICs in the terahertz regime and this investigation promises a new route for designing high-quality terahertz BIC metasurfaces.  Figure 1(a) schematically illustrates the design of the membrane metasurface. Here, the silicon membrane (permittivity ε = 11.9) is suspended in air and perforated with circular holes, and the thickness of the membrane is chosen to be T = 75 μm. A plane wave with wave vector k in the x-z plane is incident on the membrane metasurface. First, we start with the p-wave and the oblique incident angle α ranges from 0 • to 10 • . To show the spectra of the proposed metasurface, we extract the transmission of the silicon membrane via the frequency-domain analyzer of CST microwave studio. In this work the dielectric hole array was built by high resistivity silicon (ρ > 10 4 Ω cm). Compared with other electromagnetic domains, high resistivity silicon has a more special significance for the terahertz band. Because of the negligible absorption and little dispersion at terahertz frequencies, high resistivity silicon is an excellent raw material for the design of terahertz metasurfaces. The high resistivity silicon hole array avoids the Ohmic loss to achieve a higher Q factor, making this terahertz BIC metasurface more valuable in the applications such as high sensitivity terahertz sensing and large data capacity terahertz wireless communications. The adoption of high resistivity silicon also enables us to use the loss free model in the simulation without producing experimentally observable differences in the terahertz band. Besides, the possible fabrication process could be seen from references [32,33]. Figure 1(b) displays the transmission responses for the p-wave versus the incident angle. For the normal incidence, there are two obvious resonances (1 and 2) locating around 0.7 and 0.9 THz, which mainly come from the electric and magnetic dipoles (MDs), respectively. With an increment of the oblique angle, two more individual resonances (resonances 3 and 4) will arise around resonance 2. It is worth noting that resonance 2 displays a redshift while a blueshift is observed for resonances 3 and 4 as the incident angle increases. However, resonance 1 remains the same for all oblique angles.

BIC under p-wave
Next, the eigenmodes calculation of the membrane metasurface was performed to further analyze these resonances. Figure 1(c) shows the Q factor profiles of resonances 1 to 4 along the Γ-X valley. At Γ point, the Q factors of resonances 3 and 4 tend to infinity, corresponding to the case of light being perfectly confined in the membrane metasurface. With the increment of incident angle, Q-factors decrease rapidly with the broken of the symmetry. Therefore, these two symmetry-protected BICs could be proven by the lifetime and near-field profile at the Γ point, as shown in figure 2. Here, the electric and magnetic field distributions of resonances 3 could not be characterized with a rigorous classification of modes, which corresponds to four pairs of MD oscillating out of phase. From figure 2(b), the phase of each pair of MD are opposite, so the collective distribution of these dipole pairs will cancel each other. As a result, zero leakage of electromagnetic waves introduces high Q factor of resonance 3. For resonance 4 shown in figures 2(c) and (d), the electric and magnetic fields at inner and outer rings also oscillate out of phase. At this time, the inner and outer rings contribute exactly the opposite to the far field, which also leads to the Q values of resonance 4 tending toward infinity at the Γ point. When the symmetry is broken, the equilibrium of the out-of-phase oscillation of these modes will be broken, the Q factors of these two resonances drop quickly, corresponding to resonances 3 and 4 being turned to leaky modes with a high Q factor. Especially, the Q factor of resonance 4 rises again with the increment of the oblique angle, and the corresponding off-Γ point BIC will arise. The Q factor of resonance 1 keeps still at all angles.

Multipole and far-field radiation pattern analysis
To demonstrate the resonant nature of the BIC, we first start by considering the Mie-resonant multipoles of the membrane metasurface via multipole decomposition with spherical coordinate representation. The total scattering cross-section of the silicon membrane is given as [41] where a E and a M are the coefficients of the electric and magnetic multipoles, respectively; (l, m) is the order of the electric and magnetic multipoles. Figure 3(a) shows the calculated scattering cross-section of the electric and magnetic multipoles under normal incidence. Electric dipole (ED) to magnetic octupole (MO) in the legend of figure 3(a) denote the ED, MD, electric quadrupole (EQ), magnetic quadrupole (MQ), electric octupole (EO) and MO, respectively. Here, MD and EQ (even component) dominate resonance 1, and ED (odd component) supplies a small section. Where the even component is much stronger than the odd one, the first Kerker's condition is not fulfilled here. In fact, both even and odd components radiate in the vertical direction along the z-axis, which contributes to the resonance of resonance 1. For resonances 3 and 4, EO and MQ supply the most scattering contributions. Due to the zero scattering of these two multipoles along the z-axis, there will be no total radiation in the z-direction, which generates these symmetry-protected BICs. However, for resonance 2, the strongest contribution comes from ED, whereas only a small proportion of scattering is attributed to MD, EQ and MQ, resulting in an observable z-direction scattering. Furthermore, the scattering cross-section, and multipolar content under oblique incidence were also calculated. Figures 3(b) and (c) display the scattering cross-sections of the multipoles at 5 • and 10 • incidence, respectively. Because the modes under the p-wave will not affect the modes under the s-wave, it can be seen that the scattering of MD at 0.7 THz remains the same, which makes it constant at different oblique angles. In the case of resonances 2, 3 and 4, the contributions of MQ, EQ and MD all increase around 0.9 THz, respectively. With an increment of the oblique angle, the scattering of these multipoles becomes much stronger than ED and EO that dominate the response at normal incidence (over more than one order of magnitude at 10 • incidence).
To explore more details of electromagnetic multipoles, the far-field distributions of multipoles can be also categorized into magnetic and electric groups and characterized by vector spherical harmonics [37,42]: where M ml (magnetic multipoles) and N ml (electric multipoles) are the vector spherical harmonics. Equation (2) provides the correspondence between the contribution of the Mie-resonant multipole and the far-field radiation patterns. Therefore, the vector spherical harmonics analysis could be employed to demonstrate the performance of the BIC. It is worth noting that the magnitude of the scattering cross-section of the vector spherical harmonics determines the amplitude of the radiation patterns. We further expand the scattering fields at 5 • incidence in terms of spherical harmonics in figures 3(d)-(f), which are analyzed as an example of different incidence. Here, magnetic multipole M 02 plays the main role in resonance 2, electric multipoles N 22 and N 12 for resonance 3, and resonance 4 is magnetic multipole M 01 , respectively. M 02 plays a significant role in the scattering of resonance 2, whose zero radiation will cancel the scattering along the propagation direction gradually with an increment of the incident angle. Notably, for resonance 2, MQ/M 02 dominates in the scattering of the resonance at 10 • incidence because contributions of other multipoles could be neglected. The conversion of MQ will enhance the Q factor of resonance 2. For resonances 3 and 4 at oblique incidence, the scattering of N 22 and M 01 invalidates the zero radiation in the propagation direction, which destroys the formation conditions of the Γ point BIC. Then, the growing contributions of EQ and MD further turn the symmetry-protected BIC into leak modes with a high Q factor.

BIC under s-wave
So far, we have finished the multipole analysis of the BIC under the p-wave. In the following, a similar analysis of the s-wave is also performed. Figure 4(a) shows the schematic map of the s-wave incidence. In figure 4(b), it can be seen that resonance 6 will split into three individual resonances too. Resonance 7 transforms from a BIC to leaky modes with an increment of the oblique angle while resonance 6 has the opposite trend. Meanwhile, resonance 8 remains the same all the time. To understand the performance of s-wave incidence, the Q factors of the resonances were also calculated via the eigenmode analysis. At Γ point, similar to the resonances 3 and 4, the Q factors of resonances 7 tends to infinity. However, the Q factors of resonance 6 tends to infinity away from Γ point. It is obvious that the s-wave displays a similar tendency to the p-wave case except for the resonant frequencies. Moreover, the near field distribution of resonance 7 at the Γ point is calculated and displayed in figure 4(d). Both electric and magnetic modes of resonance 7 oscillate out of phase, which destroys the scattering of the resonance 7. Figures 4(e)-(g) illustrate the multipole analysis of s-wave at incident angles of 10 • . For the normal incidence, resonances 5 and 6 come from the superposition of ED, MD and EQ with MD playing a major role, and the BIC (resonance 7) arise from the destructive interference between MD and EQ with the same amplitude. From the figure 4(e), when the oblique angle increases, the enhanced contribution from MQ and ED will destroy the interference of the other modes. Figure 4(g) shows due to the broken symmetry, the enhanced scattering of electric multipoles N 11 and N 01 will destroy the balance in the destructive interference of the BIC (resonance 7), which induces the appearance of leaky modes. Meanwhile, figure 4(f) shows the rise of magnetic multipoles M 22 in resonance 6 will cancel out the scattering in the propagation direction gradually with an increment of the oblique angle, which drives the transfer of resonance 6 from a leaky mode to an off-Γ point BIC. The rest of the resonances under the s-wave share identical conclusions with the p-wave case and will not be displayed here.

Discussion and conclusion
Decoupling between the resonant modes and the continuum spectrum from the BIC can be interpreted in several ways. It could be explained as destructive interference when the coupling constants with radiating waves vanish via continuous tuning of the parameters or broken symmetry. Moreover, within the frame of the coupled-mode theory, it corresponds to nullifying the coupling coefficients between the resonant mode and all radiation channels of the surrounding space. As a basic tool for Mie resonance, multipole decomposition and vector spherical harmonics could analyze an arbitrary field distribution as a superposition of the multipoles field profiles. Similarly, the scattering of the BIC also can be interpreted as the interaction of a set of Mie-resonant multipoles.
In summary, we propose and illustrate the terahertz BIC in a membrane metasurface, and the resonant nature of the BIC is demonstrated via the multipoles and far-field radiation pattern analysis. Γ point BICs could be considered as the superposition cancellation of the multipole radiation patterns for both the p-wave and s-wave. With an increment of the oblique angle, the scattering variation of each multipoles components could change the total radiation pattern of the resonances, which drives the conversation from BICs to leaky modes. The multipole analysis of the BIC gives a new sight on the essence of BICs and could provide additional advantageous functionality for designing high-quality terahertz devices.