# Nonlinear metasurfaces governed by bound states in the continuum

###### Abstract

Nonlinear nanostructured surfaces provide a paradigm shift in nonlinear optics with new ways to control and manipulate frequency conversion processes at the nanoscale, also offering novel opportunities for applications in photonics, chemistry, material science, and biosensing. Here, we develop a general approach to employ sharp resonances in metasurfaces originated from the physics of bound states in the continuum for both engineering and enhancing the nonlinear response. We study experimentally the third-harmonic generation from metasurfaces composed of symmetry-broken silicon meta-atoms and reveal that the harmonic generation intensity depends critically on the asymmetry parameter. We employ the concept of the critical coupling of light to the metasurface resonances to uncover the effect of radiative and nonradiative losses on the nonlinear conversion efficiency.

The study of nonlinear metasurfaces has emerged recently as a new exciting direction of research, and it attracted growing interest in the photonics community due to its ability to provide a paradigm shift in nonlinear optics Li et al. (2017); Krasnok et al. (2018). Frequency conversion processes such as second- and third-harmonic generation are commonly realized in conventional nonlinear optics of macroscopic structures where high conversion efficiency is achieved by employing the concept of phase matching. Ultra-thin metasurfaces are capable to replace bulk elements yet providing reasonable conversion efficiencies Krasnok et al. (2018). Metasurfaces have the advantages over three-dimensional metamaterials in their lower optical losses and technical feasibility of fabrication, providing diverse functionalities such as efficient frequency conversion, fast optical switching, and modulation of light, and they can be employed for many applications in photonics, chemistry, material science, and biosensing.

The first study of nonlinear effects with metasurfaces revealed inherently small conversion efficiencies in the plasmonic metamaterials based on split-ring resonators. In such structures, nonlinear effects are driven by magnetic dipole resonances, but the first reported conversion efficiencies did not exceed , for the second-harmonic generation, and , for the third-harmonic generation Klein et al. (2006, 2007, 2008). The subsequent study of optically resonant dielectric metasurfaces driven by Mie resonances has improved substantially these numbers Yang et al. (2015); Tymchenko et al. (2015); Tong et al. (2016); Wang et al. (2018); Vabishchevich et al. (2018); Gao et al. (2018); Bar-David and Levy (2019) making nonlinear dielectric metasurfaces very attractive for practical applications.

All-dielectric resonant meta-optics has emerged recently as a novel research field driven by its exceptional applications for creating low-loss nanoscale metadevices Kruk and Kivshar (2017). The tight confinement of the local electromagnetic fields and multiple interferences available in resonant high-index dielectric nanostructures and metasurfaces can boost many optical effects and offer novel opportunities for the subwavelength control of light-matter interactions. In particular, recently emerged concept of bound states in the continuum (BICs) in nanophotonics enables a simple approach to achieve high-Q resonances (or quasi-BICs) for various platforms ranging from individual dielectric nanoparticles Rybin et al. (2017) to periodic arrangements of subwavelength resonators such as metasurfaces or chains of particles Hsu et al. (2016); Koshelev et al. (2018a); Koshelev et al. (2019). Moreover, very recently it was revealed Koshelev et al. (2018b) that metasurfaces created by seemingly different lattices of dielectric meta-atoms with broken in-plane inversion symmetry can support sharp high-Q resonances arising from a distortion of symmetry-protected BICs. For applications in nonlinear optics, with similarity to the recently studied isolated dielectric resonators Carletti et al. (2018), we wonder if such high-Q resonances can boost substantially nonlinear parametric processes in dielectric metasurfaces such as third-harmonic generation (THG), as shown schematically in Figs. 1(a, b). Indeed, we notice that some earlier studies suggest that the use of certain types of meta-atoms with broken symmetry, such as those shown in Fig. 1d Vabishchevich et al. (2018), Fig. 1d Gao et al. (2018) and Fig. 1e Tong et al. (2016) can lead to enhanced nonlinear effects in metasurfaces, often associated with Fano resonances. A rigorous theory of such asymmetric periodic structures developed in Ref. Koshelev et al. (2018b) revealed a close link between the Fano resonances and BICs, so we expect that the BIC-assisted resonances may provide a generic enhancement for nonlinear metasurface, including those discussed in Ref. Koshelev et al. (2018b).

In this Letter, we suggest and develop a novel concept for enhancing and tailoring nonlinear response of dielectric metasurfaces with broken in-plane symmetry exploiting the powerful physics of bound states in the continuum. We apply our approach to silicon-based metasurfaces generating third-harmonic fields. We fabricate a set of metasurfaces varying the parameters of meta-atoms and analyse how the THG efficiency depends on the degree of asymmetry. We demonstrate that the interplay of radiative and nonradiative losses can control the intensity of the third-harmonic response. In particular, we reveal that tuning the metasurface parameters to the regime of the critical coupling, when the contributions of two loss mechanisms coincide, allows to achieve the maximum efficiency of the frequency conversion in metasurfaces.

We consider a meta-atom of the designed metasurface in the form of an asymmetric pair of bars, as illustrated in Fig. 2a; the rectangular bars have the widths and , respectively. The asymmetry of the unit cell is controlled by the difference in bar widths, and it is characterized by the asymmetry parameter . We perform numerical analysis of linear response and eigenmode spectra of the asymmetric metasurface using a commercially available software for full-wave electromagnetic simulations (see Methods below). We consider the metasurface as an infinite structure with perfect periodicity and all geometrical and material parameters taken from the experimental data below. A symmetric metasurface with supports a BIC at nm with infinite quality factor (the Q factor) protected by the in-plane symmetry of the unit cell Hsu et al. (2013). The BIC is not manifested in the transmission spectrum shown in Fig. 2d at due to the symmetry incompatibility with the modes of free space. For broken-symmetry metasurfaces, the transmission curves reveal a sharp resonance with a Fano lineshape associated with a quasi-BIC with high Q factor. Figure 2c demonstrates the near-field distribution for quasi-BIC resembling the interaction of two oppositely directed magnetic dipoles with slightly dissimilar amplitudes. The evolution of the quasi-BIC radiative Q factor on the asymmetry parameter of the meta-atom follows the inverse quadratic law for small Koshelev et al. (2018b), as shown in Fig. 2e,

(1) |

where is a constant determined by the metasurface design being independent on . For larger the decrease of is going faster because the deviation from the symmetric unit cell cannot be considered as a weak perturbation.

To demonstrate how the meta-atom asymmetry shapes the TH response, we fabricate a set of six silicon metasurfaces on a glass substrate with different asymmetry parameters using the electron beam lithography with the details provided in Methods. Figure 2b shows the scanning electron microscope (SEM) image of the fabricated pattern, and the inset shows high-resolution SEM image for a single meta-atom. The pattern size of each fabricated metasurface is m m, the period nm, height is nm, bar length is nm, larger bar width is nm, and the distance between bar centres is nm. The width of the smaller bar varies in the range , , , , , nm, which corresponds to the asymmetry parameter of , , , , , and , respectively. We pump the metasurfaces from air in the wavelength range from to nm with femtosecond pulses from a coherent optical parametric amplifier (see Methods). The generated TH signal is collected in the forward direction by a high numerical aperture objective.

Figure 3 illustrates our experimental results for the THG enhancement in the broken-symmetry metasurfaces driven by BICs. The measured transmission spectra demonstrate a quasi-BIC resonance with a Fano lineshape which resonant wavelength and linewidth strongly depend on . The mode positions in the experimental linear spectra agree well with the simulation (see Fig. 2d), though the Q factor and peak transmittance are reduced, which is expected due to non-radiative losses in the fabricated sample induced by surface roughness, structural disorder, and a finite size of the sample Sadrieva et al. (2017); Jin et al. (2018). We estimate the value of the nonradiative Q factor as by using the fitting procedure (see Methods). Figure 3b shows the evolution of the measured TH intensity with respect to the meta-atom asymmetry. The highest TH intensity is observed for the intermediate value of the asymmetry parameter , whereas for larger and smaller the nonlinear signal is weaker. The dependence of the output power on the pump power for the metasurface with is shown in Fig. 3c. For average input powers below mW, the TH power clearly follows the third-order power law, while for higher pump intensities the dependence shows a typical saturation behaviour.

To simulate the nonlinear response of the fabricated metasurfaces, we introduce artificial nonradiative losses characterized by to the numerical model (see Methods). The simulated TH spectra shown in Fig. 3d demonstrate good agreement with the experimental data. Considering the simulated directivity pattern of the TH wave, we conclude that only of the total TH intensity is emitted to the zero-diffraction order in the forward direction, while the rest of signal is generated in the backward direction and scattered in high-order diffraction channels. At the same time, the effective area of illumination is about which gives an estimate for the amount of the pump power coupled to the metasurface. Using these assumptions, we calculate the intrinsic experimental efficiency of THG for the sample with demonstrating the best performance. The estimated overall conversion efficiency is for the average pump power of mW, which is similar to the previous reports on TH efficiency in Si nanostructures in the vicinity of dipole modes, composite resonances or anapoles Yang et al. (2015); Tong et al. (2016); Bar-David and Levy (2019); Carletti et al. (2018); Chen et al. (2018); Shcherbakov et al. (2014); Grinblat et al. (2016). The maximum value of TH signal is affected dramatically by the non-radiative losses due to structural imperfections, disorder and finite size of the fabricated metasurface, which can be significantly decreased by enlarging the sample footprint Yang et al. (2014) and improving the fabrication quality.

The question about the interplay between radiative and nonradiative losses attracts a wide attention in various areas of optics, especially, for describing hybrid resonator-coupler systems, where the condition of equality of internal resonator losses and coupling losses, called critical coupling, is a fundamental feature Gorodetsky and Ilchenko (1999); Xu et al. (2000); Yariv (2002). The fulfillment of this condition is favorable since it allows for complete exchange of energy between a propagating mode in a coupler device and the given resonator mode. The critical coupling principle is studied in application to plasmonics, metamaterials, and other areas Bliokh et al. (2008), and here we employ it for describing nonlinear response of subwavelength metasurfaces governed by geometric and material resonances. We develop a model based on the temporal coupled mode theory (TCMT) Fan et al. (2003) and use it for describing the intensity of the third-order nonlinear signal.

The fabricated broken-symmetry metasurfaces support a quasi-BIC mode lying in the wavelength range from to nm (see Fig. 3a) characterized by the radiative and nonradiative Q factors, where depends on the asymmetry of the meta-atom according to Eq. (1). The total Q factor is also determined by the asymmetry

(2) |

Here, is the critical value of the asymmetry parameter

(3) |

The TCMT predicts that the amplitude of the quasi-BIC mode at the resonant wavelength is determined by the Q factor and the pump intensity Maier (2006)

(4) |

The maximum amplitude is achieved exactly at the critical coupling condition , when . The intensity of nonlinear TH signal is proportional to which at the pump resonance gives

(5) |

Equation (5) shows that the highest THG conversion efficiency is achieved in the critical coupling regime, and the maximum value of efficiency scales with respect to the magnitude of nonradiative losses as .

A comparison of experimental, numerical, and analytical results for the dependence of the Q factor on the asymmetry parameter is illustrated in Fig. 4a. The total Q factor is extracted from the experimental data by the fitting procedure, with an accuracy shown with the error bars. The value of nonradiative Q factor is extrapolated as . The numerical dependence is calculated by (see Methods). Using the eigenmode analysis for small and Eq. (1), we calculate the coefficient , which gives the analytical estimate . The analytical dependence based on Eq. (2) with given parameters and is shown with a black solid line. The green area illustrates the critical coupling condition corresponding to the range of asymmetry parameters , which is blurred due to an error in calculating . The critical coupling regime is achieved for the fabricated metasurface with which agrees well with the analytical estimate for .

Figure 4b demonstrates the evolution of the peak values of the THG conversion efficiency for different asymmetries based on the experimental data, numerical calculations, and our analytical model. The experimental and numerical curves illustrate the values from Figs. 2(b, d), respectively, the analytical dependence is based on Eq. (5), and all data are normalized independently. The sample with demonstrates the best performance, as is predicted by the critical coupling model. The obtained results show high sensitivity of the nonlinear signal on the meta-atom asymmetry: nm change of results in a decrease of the THG efficiency by at least one order of magnitude.

In summary, we have developed a new approach for engineering the nonlinear response of dielectric metasurfaces composed of meta-atoms with broken in-plane symmetry closely associated with the physics of bound states in the continuum. We have employed our approach for silicon-based metasurfaces to tailor the THG efficiency depending on the degree of the unit cell asymmetry. We have demonstrated that performance of nonlinear metasurfaces can be improved by exploiting the interplay of radiative and nonradiative losses near the critical coupling condition. The similar approach can be applied to the case of the second-harmonic generation from broken-symmetry III-V semiconductor metasurfaces Vabishchevich et al. (2018) or nonlinear metasurfaces composed of arrays of zinc oxide nanoresonators designed for the nonlinear optical generation of VUV light Semmlinger et al. (2018). Thus, we believe that our general approach paves the way to smart engineering of sharp resonances in metasurfaces with inherent absorption losses and fabrication imperfections for nonlinear meta-optics and nanophotonics, towards efficient frequency conversion, fast optical switching, and modulation of light.

###### Acknowledgements.

The authors thank Sergey Kruk and Andrey Bogdanov for useful discussions and suggestions. G.L. was supported by the Guangdong Provincial Innovation and Entrepreneurship Project (2017ZT07C071), Applied Science and Technology Project of Guangdong Science and Technology Department (2017B090918001), and National Natural Science Foundation of China (11774145). Y.K. acknowledges a financial support from the Australian Research Council and the Strategic Fund of the Australian National University. K.K. was supported by the Russian Science Foundation (17-12-01581) and the Foundation for the Advancement of Theoretical Physics and Mathematics BASIS.## Methods

### Sample fabrication.

For fabrication of metasurfaces composed of silicon bar pairs on the glass substrate we first deposit thin-films of hydrogenated amorphous silicon (a-Si:H) with a thickness of nm, using plasma-enhanced chemical vapor deposition (PECVD) at a temperature of C on standard microscope coverslips. Next, the substrates are spin-coated with the positive-tone electron-beam resist (ZEP520A from Zeon Chemicals). The pattern of bar pairs is then defined by electron-beam lithography (EBL). A -nm thick Al layer was subsequently deposited by e-beam evaporation (Temescal BJD-2000), accompanied by a lift-off process in which the sample is soaked in a resist remover (ZDMAC from ZEON Co.). An array of remaining rectangular Al patterns was used as the etch mask to transfer the designed pattern into the a-Si:H film through inductively coupled plasma-reactive ion etching (Plasmalab System 100, Oxford). As etch gases, we used SF6 ( sccm) and CHF3 ( sccm). Etching was performed at C with mTorr at W induction power and W bias power. Finally, residual Al was removed using Al wet etchant (H3PO4:HNO3:CH3COOH:H2O).

### Optical experiments.

For linear the optical measurement, a halogen white light is focused onto to the silicon metasurface after passing through a Glan-Thompson linear polarizer. The transmitted light with both co-polarization and cross-polarization is then spectrally resolved by using NIR Andor spectrometer (KYMERA-328i). In the nonlinear optical measurements, we use a spectrally tunable fs laser source (repetition frequency: MHz, pulse duration: fs). The considered spectral range of the pump laser is from to nm. The linearly polarized fundamental wave with a spot size of m in diameter was normally incident onto the silicon metasurface after passing through an objective lens (NA = ). The THG signals in transmission direction are collected by an infinity-corrected objective lens (NA= ) onto an Andor spectrometer (SP500i) with a EMCCD detector. The origin of the TH signal is verified by the direct measurement of its spectrum and its power dependence which is in a cubic manner.

### Numerical calculations.

For numerical simulations of the linear and nonlinear response, we used the finite-element-method solver in COMSOL Multiphysics in the frequency domain. For simulations of the eigenmode spectra, we used the eigenmode solver in COMSOL Multiphysics. All calculations were realized for a metasurface on a semi-infinite substrate surrounded by a perfect matched layer mimicking an infinite region. The simulation area is the unit cell extended to an infinite metasurface by using the Bloch boundary conditions. All material properties including absorption losses are extracted from the ellipsometry data. The incident field is a plane wave in the normal excitation geometry polarized along the long side of the bars. For linear calculations the non-radiative losses were limited to the actual absorption losses. For nonlinear calculations, the additional non-radiative losses corresponding to were artificially introduced via correction to the extinction coefficient , where is the refractive index. The nonlinear simulations of TH response we employed within the undepleted pump approximation, using two steps to calculate the intensity of the radiated nonlinear signal. First, we simulated the linear scattering at the pump wavelength, and then obtained the nonlinear polarization induced inside the meta-atom. Then, we employed the polarization as a source for the electromagnetic simulation at the harmonic wavelength to obtain the generated TH field. The nonlinear susceptibility function was considered as a tensor corresponding to the cubic crystallographic point group with m/V. The extraction of Q factors from the experimental transmission spectra was based on the single-peak fitting to a Fano lineshape using the Levenberg-Marquardt algorithm.

## References

- Li et al. (2017) G. Li, S. Zhang, and T. Zentgraf, Nature Reviews Materials 2, 17010 (2017).
- Krasnok et al. (2018) A. Krasnok, M. Tymchenko, and A. Alu, Materials Today 21, 8 (2018).
- Klein et al. (2006) M. W. Klein, C. Enkrich, M. Wegener, and S. Linden, Science 313, 502 (2006).
- Klein et al. (2007) M. W. Klein, M. Wegener, N. Feth, and S. Linden, Optics Express 15, 5238 (2007).
- Klein et al. (2008) M. W. Klein, M. Wegener, N. Feth, and S. Linden, Optics Express 16, 8055 (2008).
- Yang et al. (2015) Y. Yang, W. Wang, A. Boulesbaa, I. I. Kravchenko, D. P. Briggs, A. Puretzky, D. Geohegan, and J. Valentine, Nano letters 15, 7388 (2015).
- Tymchenko et al. (2015) M. Tymchenko, J. S. Gomez-Diaz, J. Lee, N. Nookala, M. A. Belkin, and A. Alù, Physical review letters 115, 207403 (2015).
- Tong et al. (2016) W. Tong, C. Gong, X. Liu, S. Yuan, Q. Huang, J. Xia, and Y. Wang, Optics express 24, 19661 (2016).
- Wang et al. (2018) L. Wang, S. Kruk, K. Koshelev, I. Kravchenko, B. Luther-Davies, and Y. Kivshar, Nano letters 18, 3978 (2018).
- Vabishchevich et al. (2018) P. P. Vabishchevich, S. Liu, M. B. Sinclair, G. A. Keeler, G. M. Peake, and I. Brener, ACS Photonics 5, 1685 (2018).
- Gao et al. (2018) Y. Gao, Y. Fan, Y. Wang, W. Yang, Q. Song, and S. Xiao, Nano letters 18, 8054 (2018).
- Bar-David and Levy (2019) J. Bar-David and U. Levy, Nano letters 19, 1044 (2019).
- Kruk and Kivshar (2017) S. Kruk and Y. Kivshar, ACS Photonics 4, 2638 (2017).
- Rybin et al. (2017) M. V. Rybin, K. L. Koshelev, Z. F. Sadrieva, K. B. Samusev, A. A. Bogdanov, M. F. Limonov, and Y. S. Kivshar, Physical review letters 119, 243901 (2017).
- Hsu et al. (2016) C. W. Hsu, B. Zhen, A. D. Stone, J. D. Joannopoulos, and M. Soljačić, Nature Reviews Materials 1, 16048 (2016).
- Koshelev et al. (2018a) K. Koshelev, A. Bogdanov, and Y. Kivshar, Science Bulletin (2018a).
- Koshelev et al. (2019) K. Koshelev, G. Favraud, A. Bogdanov, Y. Kivshar, and A. Fratalocchi, Nanophotonics (2019).
- Koshelev et al. (2018b) K. Koshelev, S. Lepeshov, M. Liu, A. Bogdanov, and Y. Kivshar, Physical review letters 121, 193903 (2018b).
- Carletti et al. (2018) L. Carletti, K. Koshelev, C. De Angelis, and Y. Kivshar, Physical review letters 121, 033903 (2018).
- Hsu et al. (2013) C. W. Hsu, B. Zhen, J. Lee, S.-L. Chua, S. G. Johnson, J. D. Joannopoulos, and M. Soljačić, Nature 499, 188 (2013).
- Sadrieva et al. (2017) Z. F. Sadrieva, I. S. Sinev, K. L. Koshelev, A. Samusev, I. V. Iorsh, O. Takayama, R. Malureanu, A. A. Bogdanov, and A. V. Lavrinenko, Acs Photonics 4, 723 (2017).
- Jin et al. (2018) J. Jin, X. Yin, L. Ni, M. Soljačić, B. Zhen, and C. Peng, arXiv preprint arXiv:1812.00892 (2018).
- Chen et al. (2018) S. Chen, M. Rahmani, K. F. Li, A. Miroshnichenko, T. Zentgraf, G. Li, D. Neshev, and S. Zhang, ACS Photonics 5, 1671 (2018).
- Shcherbakov et al. (2014) M. R. Shcherbakov, D. N. Neshev, B. Hopkins, A. S. Shorokhov, I. Staude, E. V. Melik-Gaykazyan, M. Decker, A. A. Ezhov, A. E. Miroshnichenko, I. Brener, et al., Nano letters 14, 6488 (2014).
- Grinblat et al. (2016) G. Grinblat, Y. Li, M. P. Nielsen, R. F. Oulton, and S. A. Maier, Nano letters 16, 4635 (2016).
- Yang et al. (2014) Y. Yang, I. I. Kravchenko, D. P. Briggs, and J. Valentine, Nature communications 5, 5753 (2014).
- Gorodetsky and Ilchenko (1999) M. L. Gorodetsky and V. S. Ilchenko, JOSA B 16, 147 (1999).
- Xu et al. (2000) Y. Xu, Y. Li, R. K. Lee, and A. Yariv, Physical Review E 62, 7389 (2000).
- Yariv (2002) A. Yariv, IEEE Photonics Technology Letters 14, 483 (2002).
- Bliokh et al. (2008) K. Y. Bliokh, Y. P. Bliokh, V. Freilikher, S. Savelâev, and F. Nori, Reviews of Modern Physics 80, 1201 (2008).
- Fan et al. (2003) S. Fan, W. Suh, and J. D. Joannopoulos, JOSA A 20, 569 (2003).
- Maier (2006) S. A. Maier, Optics Express 14, 1957 (2006).
- Semmlinger et al. (2018) M. Semmlinger, M. L. Tseng, J. Yang, M. Zhang, C. Zhang, W.-Y. Tsai, D. P. Tsai, P. Nordlander, and N. J. Halas, Nano letters 18, 5738 (2018).