Cizhe Fang, Qiyu Yang, Qingchen Yuan, Xuetao Gan*, Jianlin Zhao,Yao Shao, Yan Liu*, Genquan Han and Yue Hao

Keywords: all-dielectric metasurface; bound states in the continuum; optical nonlinearity; topological configuration

Introduction

In the last decades, research into optical metasurfaces1,2has been intensified because of their multifunctionality,which can be employed for many applications in photonics3-11, chemistry and biosensing12,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 reflection14-19. However, due to the Ohmic loss, the conventional plasmonic metasurfaces20,21suffer from high losses and heat dissipation, which hinders their application in various functional nanoscale devices22-25, especially lowloss meta-devices. To circumvent this constraint, alldielectric metasurfaces are of paramount interest in recent years due to their high efficiency25-27. For high-index dielectric material, individual nanoparticles can possess electric and magnetic Mie-type resonances25,28,29,leading to the electromagnetic field confinement and multiple interferences inside the particles30. However,nanoparticles made of conventional optical materials(e.g. silicon (Si)) can only support a relatively low quality factor (Qfactor)31. This broad spectral response limits their applicability to the devices requiring sharp spectral features such as lasing and sensing32,33. A recently emerged concept of BIC provides a new solution to overcome this problem31-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 mechanics34,35, and then BICs are observed in acoustics,optics, microwaves, and so on36-38. A true BIC is a mathematical abstraction with infiniteQfactor 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 highQfactor39, also known as a quasi-BIC mode. The recent work also demonstrates that there exists a direct link between this leaky mode and Fano resonances39. So far, tremendous work has been explored to distinguish different mechanisms of BICs36-41. Among them, symmetry-protected BICs show aQfactor dependence on the asymmetry parameter42, which provides a potential strategy to design all-dielectric metasurfaces supporting high-Qresonance.

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 theQfactor, which are further confirmed by the experiment. Moreover, it is proved that the high-Qresonance 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 processes43,44, but also for lasing and sensing applications.

Device design and theoretical analysis

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 inxyplanes with periodic boundary conditions while perfectly matched layers are used along thezaxis45. The plane waves withx-polarization are incident from thez-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 bys,l,t,p,λ, andθ,respectively. In all simulations,pandtare fixed at 0.6 μm and 0.23 μm, respectively. Figure 1(a) describes theQfactor of the mode for differentθandl. One can see that the mode has an infiniteQfactor (more than 108) at the normal incidence. With the increase inθ, theQfactor becomes finite and drops rapidly. By contrast, this mode is insensitive to geometrical change (i.e. the side length,l) and preserves high-Qmode (i.e. ～108) with the change inl, indicating that it is a symmetry-protected BIC mode32,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 distinctz-directed magnetic dipole (MD)features in one unit cell, manifesting themselves as circular currents46. Simultaneously, the intensity of the magnetic field decreases rapidly once away from the block since the intensity of the magnetic field atz=0.345 μm nearly decreases by an order of magnitude compared with that atz=0 μm. These modes do not radiate in the vertical direction along thez-axis because the normal radiative decay is compensated by driving terms arising from the local field at the position of the resonator47,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 theQfactor, the different defects are introduced in the blocks (see the bottom inset of Fig.1(c). The size of the defect is only controlled byasincebis fixed at 0.1 μm. As shown in Fig. 1(c), with the increase in the defect size, theQfactor decreases sharply.When the defect size continues to increase (i.e.a＞0.2 μm), theQfactor begins to increase. TheQfactor does not monotonically decrease with the increase in the size of the defect, which is not consistent with the previous report42. We ascribe this phenomenon to the existence of another BIC mode ata=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, theQfactor is independent ofland 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-Qresonances also originate from the symmetry-protected BIC modes. Whenθcontinues to increase, theQfactor 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 largeQfactor. It is worth noting that compared to the symmetry-protected BIC,the resonance-trapped BIC is relatively less sensitive toθand thus remains a highQfactor at aroundθ=19°32,40. It occurs in the metasurface due to the destructive interference38,40and it is sensitive to geometric dimensions (i.e.l). To prove it, theQfactor is modeled as a function oflatθ=19° (see the top inset in Fig. 1(d)). Compared with the influence ofθ, theQfactor decreases more dramatically oncelleaves the optimal value. It is also different from the insensitivity of theQfactor tolin the symmetry-protected mode. The further study of this mode is beyond the scope of this work.

Fig. 1 | (a) Simulated Q factor and resonant wavelength as a function of θ (top) and l (bottom), respectively. Inset: the schematic view of the model (a square block with l=0.4 μm). (b) Simulated transmission spectrum as a function of the angle and wavelength. Inset: the magnetic field profiles under the normal incidence condition, where white arrows correspond to electric field vectors and the black wireframe outlines the single unit. (c) Modeled the Q factor as a function of a for a broken-symmetry square block (l=0.4 μm). Bottom inset: the schematic image of the concave Si block. The top two insets show the magnetic field distributions at around a=0 μm and a=0.4 μm, respectively. (d) Simulated Q factor and resonant wavelength as a function ofθ(top) and l(bottom) for the BIC mode at a=0.4 μm. The top inset: simulated Q factor as a function of l at θ=19°. The bottom left inset: the schematic view of the model (a rectangle block with l1=0.3 μm). The bottom right inset: the magnetic field profile atθ=0°. The scaling parameter is set as 1.2 in (a-d).

Owing to the fact that there exist two different BIC modes (i.e.a=0 μm anda=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 largerQfactor. However, the behavior of theQfactor 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, theQfactor 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 theQfactors and field profiles(see the magnetic field distributions shown in Figs. 1(b)and 1(d)). For the case ofa=0.4 μm, there are two additional BICs (indicated by white circles) on thekyaxis,which contributes to the variation in theQfactor distribution alongkyaxis. 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 resonances49. That indicates the reason for different decay rates of theQfactor in two structures. As a result, the corresponding leaky modes (i.e.β≠0) are affected and it finally leads to the different dependence ofQfactor onβ, as shown in Fig. 2(a). Moreover, according to ref.50, the generation,evolution, and annihilation of BICs can be governed by varying structural parameters, which provides theoretical support to adjust the relationship betweenQfactor 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 theQfactor 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.

Results and discussion

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 differents. 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. Thes-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 thez-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 theyzplane 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 theQfactor. Due to the process tolerance, only the proposed metasurfaces witharanging 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 ata=0.02 μm and 0.32 μm. For comparison,Qfactors for identical structures are also calculated and shown in Fig. 3(c). We observe a very good agreement between the experimental and calculated results. TheQfactor decreases with the increase inaat first and then increases along witha.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 measuredQfactor is modeled as a function ofβfor two definitions. It can be seen that theQfactor increases with a decrease inβ. More importantly, there indeed exist two different relationships (i.e. the slopes are-0.35 and -0.69) between theQfactor 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 changinga, both structures can be realized. In this case, for the sames, the structure holds a smallerAfor definition 2. It means that the introduced perturbation is more easily affected by the non-ideal factor in the preparation process.

Fig. 2 | (a) Dependence of the Q factor on asymmetry parameter β for both designs (log-log scale). Inset: two definitions of the asymmetry parameter. (b) Simulated Q in momentum space at a=0 μm (left) and a=0.4 μm (right). Inset: the designed patterns shown in the xy plane for both cases. (c) Schematic image of a period metasurface with concave nanopillars. Inset: device fabrication process. (d) SEM image of the fabricated metasurface.

Next, we aim to discuss the implementation of BICsupporting metasurfaces. As shown in Fig. 3(c), the measuredQfactor is several times smaller than that calculated one. Four major factors that lead to the low-Qresonance are listed as follows. First, as an experimental sample, surface cleanliness is inevitable. It can introduce an inherent loss caused by scattering51, which finally leads to the attenuation in theQfactor. 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, theQfactor 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 one51. To verify this, we calculated and plotted the evolution of theQfactor and resonant wavelength as a function of the thickness of the substratetsub. From Fig. 3(d), we can see that theQfactor rapidly drops with the existence of the substrate accompanied by the redshift in resonant wavelength. As the thickness continues to increase, theQfactor becomes relatively insensitive totsub. 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 largerQfactor 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 theQfactor. For bettering understanding, we simulated the responses of the different size arrays and their correspondingQfactors (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-Qmodes 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. largerQfactor). Above all, the low-Qmeasured 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 highQfactor.

Fig. 3 | (a) Measured transmission spectrum at different s (20×20 array). Inset: The magnetic field distribution in the yz plane for the metasurface with the substrate. (b) Polarization dependence measurements of the transmission (Texp) and reflection (Rexp) spectra at a wavelength of 1.55 μm and their corresponding calculations (i.e. Tsim and Rsim). The arrow indicates the direction of the x-axis. (c) Top: Modeled the calculated and measured Q factor as a function of a(s=1.2 and b=0.1 μm). Inset: SEM image of the single unit at a=0.02 μm and a=0.32 μm, respectively. Bottom: The evolution of Q factor at different α for two definitions (log-log scale). (d) Numerically simulated dependence of the resonant wavelength and Q factor on tsub(tsub=0 μm means that there is no substrate). The top inset shows the comparison of the measured Q factor for different dimensions and array sizes. Structure 1: l=0.4 μm, a=0.15 μm, b=0.1 μm, and s=1.15. Structure 2: l=0.5 μm, a=0.1 μm, b=0.1 μm, and s=1.02. The bottom inset shows the SEM images of the prepared structures before and after the process optimization.

Fig. 4 | Magnitude plots of the on-resonance magnetic field for the isolated resonator and 3×3, 5×5, 7×7, and 9×9 arrays. The magnetic field distributions are obtained at the vertical mid-plane of the blocks and the magnitude is indicated by the height. The structural parameters are l =0.4 μm, a=0.15 μm, b=0.1 μm, s=1.2. The table on the right shows the corresponding Q factors.

As we mentioned above, an optimized metasurface was fabricated to improve theQfactor. 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 formula52:

whereεrepresents the scale of reduced energy,qis the asymmetry factor,ω0is the resonant frequency, andγrepresents the damping rate. By fitting the data with Eq.(1), we can obtain the fitting parameters ofω0=1547.85 andγ=0.219. AQfactor as high as 3534 can be obtained,which can be further improved by introducing smaller defects. As one of the most promising applications53, 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-Qresonances and strong field confinement21,54,which are considered as the dominant factors in nonlinear optics instead of phase matching at the nanoscale55.On the one hand, for this nonlinear response, the intensities of the THG scale approximately with the factor(Q/V)356, whereVis the mode volume. An enhancement inQfactor 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.

Fig. 5 | (a) Measured transmission spectrum of the proposed metasurface (200×200 array) and its fitting of Fano function. Optimized structural parameters are p=0.6 μm, l=0.5 μm, a=0.1 μm, b=0.1 μm, s=1.05. (b) Comparison of THG in the proposed structure (red line) and a flat Si film with the same thickness (black line). The intensity of the THG in the flat Si film is amplified 100 times. (c) The electric field profiles at the resonance (left) and the off-resonance (right), respectively. (d) The corresponding polarization dependence measurements (normalized to the peak value of the resonance).

Conclusions

In summary, we have successfully designed and fabricated Si metasurfaces with highQresonances, which have elements of broken-symmetry blocks. The MD modes of the blocks contribute to these resonances. Both simulation and experiment results reveal that theQfactor 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 highQresonance, 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-Qdynamic resonances, and they constitute a significant step towards the development of high-Qresonant photonic applications such as lasers, sensors, and filters.