In the last decades, research into optical metasurfaces^{1, 2} has been intensified because of their multifunctionality, which can be employed for many applications in photonics^3-11, chemistry and biosensing^{12, 13}. Metasurfaces are ultrathin planar arrangements consisting of subwavelength spaced scatterers with different shapes and sizes. They can modify the amplitude, phase, and polarization of propagating optical waves upon transmission or reflection^14-19. However, due to the Ohmic loss, the conventional plasmonic metasurfaces^{20, 21} suffer from high losses and heat dissipation, which hinders their application in various functional nanoscale devices^22-25, especially low-loss meta-devices. To circumvent this constraint, all-dielectric metasurfaces are of paramount interest in recent years due to their high efficiency^25-27. For high-index dielectric material, individual nanoparticles can possess electric and magnetic Mie-type resonances^{25, 28, 29}, leading to the electromagnetic field confinement and multiple interferences inside the particles³⁰. However, nanoparticles made of conventional optical materials (e.g. silicon (Si)) can only support a relatively low quality factor (Q factor)³¹. This broad spectral response limits their applicability to the devices requiring sharp spectral features such as lasing and sensing^{32, 33}. A recently emerged concept of BIC provides a new solution to overcome this problem^31-34.

BICs represent a wave phenomenon of modes, which have the energy lying in the delocalized states inside the continuum. This concept firstly appears in the quantum mechanics^{34, 35}, and then BICs are observed in acoustics, optics, microwaves, and so on^36-38. A true BIC is a mathematical abstraction with infinite Q factor and the access to radiative channels is completely forbidden. In fact, for a practical device (i.e. finite-size structure), all possible BICs turn into the leaky mode with a high Q factor³⁹, also known as a quasi-BIC mode. The recent work also demonstrates that there exists a direct link between this leaky mode and Fano resonances³⁹. So far, tremendous work has been explored to distinguish different mechanisms of BICs^36-41. Among them, symmetry-protected BICs show a Q factor dependence on the asymmetry parameter⁴², which provides a potential strategy to design all-dielectric metasurfaces supporting high-Q resonance.

In this paper, we introduce the defects in the unit cell of the periodic metasurfaces based on Si. The quasi-BIC modes are emerged, which manifest themselves as Fano resonances. By employing the finite element method (FEM), the physical mechanism is studied. The theoretical results also show the relationship between the asymmetry parameter and the Q factor, which are further confirmed by the experiment. Moreover, it is proved that the high-Q resonance can dramatically enhance the intensity of the third harmonic generation (THG) signal. We believe that these results can be not only useful for improving the nonlinear conversion processes^{43, 44}, but also for lasing and sensing applications.

First, we consider a design of the periodic structure with in-plane symmetry to support BIC modes (see the inset of Fig. 1(a)). The mode properties are studied by utilizing the FEM. For simplicity, a single unit in free space is simulated in xy planes with periodic boundary conditions while perfectly matched layers are used along the z-axis⁴⁵. The plane waves with x-polarization are incident from the z-axis and the refractive index of silicon is set as 3.48. The scaling parameter, the side length of the square block, the thickness of the block, the period, incident wavelength and angle are denoted by s, l, t, p, λ, and θ, respectively. In all simulations, p and t are fixed at 0.6 μm and 0.23 μm, respectively. Figure 1(a) describes the Q factor of the mode for different θ and l. One can see that the mode has an infinite Q factor (more than 10⁸) at the normal incidence. With the increase in θ, the Q factor becomes finite and drops rapidly. By contrast, this mode is insensitive to geometrical change (i.e. the side length, l) and preserves high-Q mode (i.e. ~10⁸) with the change in l, indicating that it is a symmetry-protected BIC mode^{32, 40}. For better understanding, the transmission spectrum dependence on angle of incidence is plotted in Fig. 1(b). It can be seen that the resonance (manifested as dips in the transmission) narrows and then vanishes when the incident angle approaches 0°, accompanied by the blue shift in the resonant wavelength. It corresponds to the decline in radiation loss and the emergence of the symmetry-protected BIC. To gain a deeper insight into the physics of this decline, the magnetic field profiles at θ=0° are displayed in the inset of Fig. 1(b). The electric field presents distinct z-directed magnetic dipole (MD) features in one unit cell, manifesting themselves as circular currents⁴⁶. Simultaneously, the intensity of the magnetic field decreases rapidly once away from the block since the intensity of the magnetic field at z=0.345 μm nearly decreases by an order of magnitude compared with that at z=0 μm. These modes do not radiate in the vertical direction along the z-axis because the normal radiative decay is compensated by driving terms arising from the local field at the position of the resonator^{47, 48}. After introducing the defect, this compensation is broken and the ideal BIC transforms into the leaky mode with the sharp resonance. The radiation powers of the MD moments are the dominant factor for these resonances. As can be noticed from the top two insets of Fig. 1(c), a distinct MD can be observed in the block with broken in-plane symmetry. To interpret the influence of the detect size on the Q factor, the different defects are introduced in the blocks (see the bottom inset of Fig. 1(c). The size of the defect is only controlled by a since b is fixed at 0.1 μm. As shown in Fig. 1(c), with the increase in the defect size, the Q factor decreases sharply. When the defect size continues to increase (i.e. a>0.2 μm), theQ factor begins to increase. The Q factor does not monotonically decrease with the increase in the size of the defect, which is not consistent with the previous report⁴². We ascribe this phenomenon to the existence of another BIC mode at a=0.4 μm (see the bottom left inset in Fig. 1(d)) and eigenvalue simulations are carried out to verify this speculation. Figure 1(d) shows that the mode is also extremely sensitive to θ as its quality factor drops sharply for a small θ. However, the Q factor is independent of l and remains high. Meanwhile, the electromagnetic field is perfectly localized in the block (see the bottom right inset of Fig. 1(d)), implying that these high-Q resonances also originate from the symmetry-protected BIC modes. When θ continues to increase, the Q factor begins to increase along with θ. It is believed that this phenomenon results from the existence of a resonance-trapped mode at θ=19° (indicated by the red dotted circle), which also shows a large Q factor. It is worth noting that compared to the symmetry-protected BIC, the resonance-trapped BIC is relatively less sensitive to θ and thus remains a high Q factor at around θ=19°^{32, 40}. It occurs in the metasurface due to the destructive interference^{38, 40} and it is sensitive to geometric dimensions (i.e. l). To prove it, the Q factor is modeled as a function of l at θ=19° (see the top inset in Fig. 1(d)). Compared with the influence of θ, the Q factor decreases more dramatically once l leaves the optimal value. It is also different from the insensitivity of the Q factor to l in the symmetry-protected mode. The further study of this mode is beyond the scope of this work.

Owing to the fact that there exist two different BIC modes (i.e. a=0 μm and a=0.4 μm), we define the asymmetry parameter β for both cases as seen in the inset of Fig. 2(a). It can be seen that for small values of β, a smaller asymmetry parameter corresponds to a larger Q factor. However, the behavior of the Q factor on β is different, where the slopes of the fitting lines are -0.61 and -1.8 for definition 1 and definition 2, respectively. This phenomenon can be explained by the topological configuration of BICs. Taking into account that the leaky mode arises when the BIC mode couples to far-field radiation, we focus on the model supporting the ideal BIC mode (i.e. β=0). The models for definitions 1 and 2 are shown in the inset of Fig. 2(b), corresponding to the square block (see the inset of Fig. 1(a)) and rectangle block (see the bottom left inset of Fig. 1(d)), respectively. To verify this, the Q factor distributions around the Γ point are calculated and plotted for both scenarios. A comparison is shown in Fig. 2(b), where BICs are generated at the Γ point for both cases. It arises because coupling constants with all radiating waves vanish due to symmetry. At this point, the light becomes perfectly confined in the metasurfaces, as is evident from the Q factors and field profiles (see the magnetic field distributions shown in Figs. 1(b) and 1(d)). For the case of a=0.4 μm, there are two additional BICs (indicated by white circles) on the k_y axis, which contributes to the variation in the Q factor distribution along k_y axis. The BIC at off-Γ corresponds to the resonance-trapped mode in Fig. 1(d). Since the BICs behave like a topological defect, it causes the difference in the topological configuration between both structures. It is known that the topological configuration of BICs can control radiative losses of nearby resonances⁴⁹. That indicates the reason for different decay rates of the Q factor in two structures. As a result, the corresponding leaky modes (i.e. β≠0) are affected and it finally leads to the different dependence of Q factor on β, as shown in Fig. 2(a). Moreover, according to ref.⁵⁰, the generation, evolution, and annihilation of BICs can be governed by varying structural parameters, which provides theoretical support to adjust the relationship between Q factor and β. Above all, the modulation of the topological configuration of BICs is also achievable in the metasurface. Therefore, this mechanism provides a simple but important tool to engineer the topological configurations and thus the behavior of the Q factor on β, which also provides a platform for developing the realistic BIC devices upon metasurface.

To confirm such theoretical findings, metasurfaces with different dimensions using electron-beam lithography (EBL), followed by inductively coupled plasma (ICP) etching, were fabricated, as shown in Fig. 2(c). Figure 2(d) presents the scanning electron microscope (SEM) images of the fabricated device. It should be pointed out that the dimensions of the blocks are always smaller than their design values, resulting from the fabrication tolerance, and it leads to a blue shift in resonant wavelength.

In this section, we focus on the discussion of measured results. The linear optical characteristics of the fabricated samples were measured using a home-built setup, where a continuous-wave (CW) laser was focused onto the sample, and the transmitted or reflected light was monitored using a spectrometer. Figure 3(a) depicts the transmission spectra for different s. One can observe that transmission spectra show sharp resonances and the resonance positions cover the spectral range from 1.2 μm to 1.5 μm. The s-dependence on resonant wavelength is also valid beyond this spectrum range. It is noted the measured resonance appears as a peak in the transmission spectrum, different from the simulation results. We attributed this phenomenon to the result of the combined action of the Fano resonance and the substrate. It is known to us that the quasi-BIC manifests itself as the Fano resonance in the transmission. For the prepared metasurface, there exists a substrate. On the one hand, the substrate introduces the perturbation (i.e. asymmetry along the z-axis). At around the resonant wavelength, the mode leaks to the substrate (i.e. the incident direction), which increases the transmission. For better understanding, the magnetic field distribution in the yz plane is plotted in the inset of Fig. 3(a). It is clear that at the resonant wavelength, there is more leakage into the substrate. On the other hand, the substrate suppresses the transmitted light. As a result, the “dip” is suppressed and the Fano resonance exhibits as a “peak”. We also perform the polarization dependences of transmission and reflection spectra. As the polar plots are shown in Fig. 3(b), both transmission and reflection spectra show a distinct polarization dependence, consistent with the results of the simulation. Now, we turn toward the influence of the defect size on the Q factor. Due to the process tolerance, only the proposed metasurfaces with a ranging from 0.08 μm to 0.3 μm can meet the design requirements. The insets of Fig. 3(c) show that the shapes of the unit cells are greatly changed at a =0.02 μm and 0.32 μm. For comparison, Q factors for identical structures are also calculated and shown in Fig. 3(c). We observe a very good agreement between the experimental and calculated results. The Q factor decreases with the increase in a at first and then increases along with a. Notably, there is a large difference between the experimental and calculated results. It will be explained in the next paragraph. As shown in the bottom figure in Fig. 3(c), the measured Q factor is modeled as a function of β for two definitions. It can be seen that the Q factor increases with a decrease in β. More importantly, there indeed exist two different relationships (i.e. the slopes are –0.35 and –0.69) between the Q factor and asymmetry parameter β, which confirms the results of the eigenvalue analysis. It is noted that the slopes of the measured results are different from that of the calculated results, especially for definition 2. We attribute this difference to the non-ideal factor in the experiment. According to the definition, the structures remove and add the perturbation for definitions 1 and 2, respectively. For convenience, we designed the pattern by removing the perturbation in the middle of the structure (see the bottom inset of Fig. 1(c)). By changing a, both structures can be realized. In this case, for the same s, the structure holds a smaller A for definition 2. It means that the introduced perturbation is more easily affected by the non-ideal factor in the preparation process.

Next, we aim to discuss the implementation of BIC-supporting metasurfaces. As shown in Fig. 3(c), the measured Q factor is several times smaller than that calculated one. Four major factors that lead to the low-Q resonance are listed as follows. First, as an experimental sample, surface cleanliness is inevitable. It can introduce an inherent loss caused by scattering⁵¹, which finally leads to the attenuation in the Q factor. The second one is fabrication tolerance, which has a great impact on the shape and size of the unit cell. As shown in the bottom inset of Fig. 3(d), after the process optimization, the surface is cleaner and the edges of the blocks are smoother in the structure. Importantly, the Q factor can be improved from 250 to 413. Through process optimization, it is possible to minimize the influences of the surface cleanliness and process preparation. Another long-standing factor is the substrate. The existence of the substrate will generate a field penetrating, which turns a BIC mode into the resonant one⁵¹. To verify this, we calculated and plotted the evolution of the Q factor and resonant wavelength as a function of the thickness of the substrate t_sub. From Fig. 3(d), we can see that the Q factor rapidly drops with the existence of the substrate accompanied by the redshift in resonant wavelength. As the thickness continues to increase, the Q factor becomes relatively insensitive to t_sub. Besides, the array size also plays a crucial role in achieving a sharp resonance. As depicted in Fig. 3(d), the prepared metasurface with larger array sizes shows a larger Q factor for both structures. It can be explained that with the increase in the array size, a finite array of the unit cell becomes an infinite periodic structure, which corresponds to an ideal BIC with an infinite value of the Q factor. For bettering understanding, we simulated the responses of the different size arrays and their corresponding Q factors (see Fig. 4). The magnitude of the magnetic field increases along with the array size and the innermost blocks of the array show larger field enhancement than that at the edge of the array. As we mentioned above, these high-Q modes hardly radiate in the vertical direction because the normal radiative decay is compensated by driving terms arising from the local field at the position of the resonator. However, for a finite array, the unit cells at the edge can only be partially compensated, leading to extra leakage. With the increase in the array size, the proportion of the unit cell at the edge can be reduced, corresponding to a less leakage (i.e. larger Q factor). Above all, the low-Q measured results are caused by the existence of the substrate and the finite array size, which are inevitable in practical application at the nanoscale. Therefore, process optimization and reasonable design can be more effective solutions to achieve a high Q factor.

As we mentioned above, an optimized metasurface was fabricated to improve the Q factor. Due to the resolution of the spectrometer, a narrowband tunable continuous-wave (CW) laser was utilized, and the transmitted optical powers were monitored using a telecom-band photodiode. Figure 5(a) illustrates the transmission spectrum at a wavelength of around 1.55 μm. The transmission peak is fitting by a Fano formula⁵²:

$ {F(ε ) = \frac{{{{(q + ε )}^2} + {\gamma ^2}}}{{1 + {ε ^2}}}\;,} $ (1)

$ {ε = \left( {\omega - {\omega _0}} \right)/\gamma \;,} $ (2)

$ {Q = {\omega _0}/2\gamma\; ,} $ (3)

where ε represents the scale of reduced energy, q is the asymmetry factor, ω₀ is the resonant frequency, and γ represents the damping rate. By fitting the data with Eq. (1), we can obtain the fitting parameters of ω₀=1547.85 and γ=0.219. A Q factor as high as 3534 can be obtained, which can be further improved by introducing smaller defects. As one of the most promising applications⁵³, we consider nonlinear effects in the proposed metasurface. For the nonlinear measurements, the metasurface was pumped by a pico-second pulsed laser, and the transmitted THG signals were monitored using a spectrometer. As we can see from Fig. 5(b), the intensity of THG signals in the proposed structure is 368 times larger than that in the flat Si film. This enhancement can be ascribed to high-Q resonances and strong field confinement^{21, 54}, which are considered as the dominant factors in nonlinear optics instead of phase matching at the nanoscale⁵⁵. On the one hand, for this nonlinear response, the intensities of the THG scale approximately with the factor (Q/V)³⁵⁶, where V is the mode volume. An enhancement in Q factor at the fundamental frequency can be further amplified in the third-order nonlinear interaction. On the other hand, strong light confinement means that the photons are localized in a small space for a long time, which can enhance various light-matter interactions. Also, the electric field distributions at both the resonance and the off-resonance are plotted in Fig. 5(c). As we can see, there is a significant attenuation in the intensity of the electric field once away from the resonant wavelength, corresponding to the decrease in the intensity of the THG. Figure 5(d) compares the polarization dependences of both THG. The pronounced THG from the flat Si film (black line) has another polarization dependence (indicated by the black arrow). This polarization dependence results from the THG in the sapphire substrate. It is worth noting that this polarization dependence disappears in the proposed structure (red line). It is because Si blocks possess the strong electromagnetic field, resulting in the significant enhancement in the intensity of the THG. It provides further evidence of the improvement in the THG signal caused by the proposed structure. The nonlinear response can be further enhanced and tailored by varying the dimensions of the metasurfaces.

In summary, we have successfully designed and fabricated Si metasurfaces with high Q resonances, which have elements of broken-symmetry blocks. The MD modes of the blocks contribute to these resonances. Both simulation and experiment results reveal that the Q factor increases when the size of the defects decreases and this relationship can be altered by changing the dimensions of the structure. As an example of the applications of the high Q resonance, high-efficiency THG is implemented from the silicon metasurface, which presents an enhancement factor of 368. The results provided by this work pave a way to manipulate BICs and realize high-Q dynamic resonances, and they constitute a significant step towards the development of high-Q resonant photonic applications such as lasers, sensors, and filters.