Enhanced third-harmonic generation induced by nonlinear field resonances in plasmonic-graphene metasurfaces

Nonlinear metasurfaces offer new paradigm for boosting optical effect beyond limitations of conventional materials. In this work, we present an alternative way to produce pronounced third-harmonic generation (THG) based on nonlinear field resonances rather than linear field enhancement, which is a typical strategy for achieving strong nonlinear response. By designing and studying a nonlinear plasmonic-graphene metasurface at terahertz regime with hybrid guided modes and bound states in the continuum modes, it is found that a THG with a narrow bandwidth can be observed, thanks to the strong resonance between generated weak THG field and these modes. Such strong nonlinear field resonance greatly enhances the photon-photon interactions, thus resulting in a large effective nonlinear coefficient of the whole system. This finding provides new opportunity for studying nonlinear optical metasurfaces.

Third-order harmonic generation (THG), as a typical example of frequency conversion, is one of the most fundamental manifestations of nonlinear optical response, and offers opportunities for many applications in photonics, material science, all-optical signal processors and biosensing. A variety of nonlinear metasurfaces [25][26][27][28][29] have been demonstrated theoretically and experimentally to enhance the THG, including plasmonic metasurfaces, alldielectric metasurfaces and dielectric-metal hybrid structures. In an optical nonlinear system, including metasurfaces, the linear and nonlinear process happen at the same time, and they are in no particular order. The linear process, such as reflection, refraction or scattering, is the result of photon-electron interactions, while the nonlinear process is due to the photon-photon interactions that are intrinsically weak and can only be realized at high intensity of light.  , the interaction of such produced THG and the specific resonances and the associated feedback between linear and nonlinear processes lead to enhanced effective nonlinear coefficient.
So far, most of nonlinear metasurfaces for boosted THG response rely on the enhancement of linear local electric field (i.e., photon-electron interactions), and its typical physical process is schematically shown in Fig. 1 ff  . In this work, we design and study a plasmonic-graphene nonlinear metasurface, which can exhibit strong THG response not in this typical way. We contribute such strong THG response to nonlinear coefficient amplification originating from nonlinear field resonance. Figure 1(b) schematically illustrates its concept. Likewise, in order to obtain a strong interaction between light and matter, a linear metasurface system with a specific resonant frequency R f is designed. Distinguished from the above-mentioned typical way, the frequency of incident FF light is /3 ff  , far from the specific resonance. Due to the added nonlinear medium, the incident FF light with /3 FF R ff  will produce a weak THG with 3 TH FF R f f f  , which implies that the third harmonic frequency exactly coincides with the resonance frequency of the system. As a result, the interaction of such produced THG and the specific resonances leads to boosted THG. This seemingly simple operation is actually highly nontrivial because the THG field resonance greatly enhances the photon-photon interactions and effectively amplify the nonlinear coefficient of the entire system, (i.e., (3) (3) e   ). The two ways for strong THG in Fig. 1 seems to be the result of a simple interchange of the order of linear and nonlinear processes, while they differ fundamentally in physics. We will show that the THG based on our discussed way in some cases may be more efficient than that based on the linear local field enhancement corresponding to the typical way.

Results and discussions
Let us start with the considered nonlinear metasurface, as shown in Fig. 2(a). The designed metasurface is a periodic array of gold (Au) cylinder with a radius R, laid on a graphene/dielectric/metal substrate along the y direction. As the designed metasurface is 2D structure, the cylinders are infinite along x direction. Figure 2(b) shows the details of a unit cell with a period p. The optical property of Au at terahertz regime is described by 22 / , where 0  and (3)  are the linear term and the third-order nonlinear term, respectively [31]. For the linear term 0  , when /2 F E   , the interband transition of electrons in graphene is forbidden and only the intraband transition is left to dominate its optical response with incident light. As a result,   is the angular frequency, e is the electron charge, is the reduced Planck constant, and  is the relaxation time which is estimated as 10 −13 s [32]. For the nonlinear term (3)  , its exact expression is taken from Ref. [33]. A transverse-magnetic (TM) polarized FF light (i.e., the magnetic field only along the x direction) is normally incident on this nonlinear metasurface from air. COMSOL Multiphysics is used to study the nonlinear effect of the proposed metasurface, and the nonlinear graphene is modeled by a nonlinear surface current with where FF E and TH E are the electric fields of the FF and the TH, respectively. Based on the simulation scheme used in Ref. [28], we can obtain THG radiation in COMSOL Multiphysics.
Intuitively, due to localized surface plasmons (LSPs) [35], the electromagnetic (EM) energy of incident light can be guided into the bottom of each Au cylinder, leading to a strong local field near it. Rm   are numerically found to design the metasurface with specific resonances in this region. Clearly, owing to the LSP resonances (see the field pattern in the right panels), a resonance peak appears in each absorption spectrum (see the left panel). According to the typical way, this strong local field will significantly enhance the interaction between light and nearby nonlinear graphene, producing a strong THG. However, the understanding of THG enhancement effect found in our designed nonlinear metasurface is not so straightforward. Figures 2(f) . The corresponding electric field patterns of THG are plotted in the right panels. Most energy of the THG is bounded at the surface of the graphene/PMMA, featured with a standing wave with nine nodes. Note that although the FF of THG resonance is very close to the LSP resonance frequency and the difference is very small, the THG peak in the nonlinear spectrum is not caused by the LSP resonance. We will show later that such frequency difference is ubiquitous, and can be altered further by changing the geometric parameters, such as the thickness of dielectric layer and the radius of the metal cylinders (discussed later in Fig. 5 and Fig. 6). To uncover the mechanism of seemly unusual THG, the photonic band diagram of the considered structure was numerically calculated, as shown in Fig. 3(a). Here we take the case of 1 dm   as an example for illustrations. In calculations, we only considered the optical properties of the purely linear system, thus removing the nonlinear term of graphene, as the nonlinear effect is weaker compared to the linear one. Owing to the interaction between the metal cylinder array and the dielectrics, a plasmonic-photonic band is plotted in Fig. 3(a), and we call it hybrid guided mode for convenience because these hybrid modes can propagate freely in the deep subwavelength dielectric layer along the y direction. The eigenfrequency at  point ( 0 y k  ) is about 18.0 THz and the corresponding eigenmode profile is shown in Fig. 3(b), which is exactly consistent with the field distribution of THG at 6.0 FF f THz  in Fig. 2(g). In addition, we also calculated the Q factor of this band, as shown in Fig. 3(c). The Q factor at ky = 0 is the largest, approximately Q = 200, which is almost the same as the Q value calculated from THG spectrum in Fig. 2(g). Therefore, these results of photonic band diagram clearly reveal that the physical mechanism of the significant THG in the designed nonlinear metasurface is attributed to hybrid guided mode resonance of THG, rather than the enhancement of the linear local field. In physics, the THG effect depends on the term ( 11 3.7 10   ), respectively. In fact, the linear field enhancement in both cases is comparable. But in the case of 6.0 FF f  THz, the nonlinear field resonance stemming from hybrid guided mode resonance of THG greatly enhances the photon-photon interactions, thus resulting in an effective third order susceptibility greater than the intrinsic one of graphene, i.e., (3) ( . Therefore, by designing mode resonance of third harmonic frequency in a linear system, the boosted THG could be obtained as well. Alternatively, as our designed metasurface has mirror symmetry with respect to z axis, the complete decoupling between incident wave and the eigenmode of odd symmetry can happen. As a result, the quasi-BIC mode or BIC mode can exist in our metasurface system, which can greatly boost THG due to their extremely strong resonances [29]. To find it, we numerically analyzed the photonic band structure of the considered structure. As shown in Fig. 4(a), it is a band of odd mode (see Fig. 4(b)), where the mode at 0 y k  indicated by a red circle is quasi-BIC mode with eignefrequency at about 4.5 THz. This point is confirmed by calculating the quality factor Q of the whole band structure. As shown in the Fig. 4(c), the Q factor reaches an extremely large value at 0 y k  . Note that if the Ohmic's loss in metal is removed, i.e., 0   , the obtained Q value goes to infinity as 0 y k  . Now let us explore what happens to THG caused by BIC modes in the current case. Based on the two different ways for THG shown in Fig. 1, two resonant  f  THz is caused by the generated THG field resonance with above mentioned LSPs. Specifically, the peak around 1.6 THz, which is very close to 1.5 TT f  THz (i.e., 4.5/3 THz), is caused by the resonance between THG field and BIC mode. While the peak around 4.8 THz, which is slightly off the eigenfrequency 4.5 THz of BIC, is caused by the local field enhancement resulting from the resonance between the incident light and BIC. Owing to the nature of BIC itself, a small deviation is seen in both cases, which is consistent with the results observed in [29], that is, the strongest nonlinear effect does not occur exactly at BIC point in the spectrum, but at a position near it. As we know, the Q factor of BIC is infinite, which means that the incident EM energy cannot effectively couple to the optical system through BIC resonance. Strictly speaking, the energy coupled to the system is zero, which results in extremely weak THG. Therefore, in this sense, a true BIC is not an ideal choice for enhancing the interaction between light and matter. By contrast, the quasi-BICs with a moderate Q factor are more favorable to enhance the interaction between light and matter. Moreover, one can see that in the current case, the THG output power at 1.6 TT f  THz (the corresponding conversion efficiency is about 12

10  
) is larger than that at 4.8 TT f  THz. It implies that the THG based on our proposed way may be more efficient than that based on the typical way. The THG resonance frequency can be adjusted by changing the thickness of the dielectric layer or the radius of the metal cylinders, because the hybrid guided modes and BIC modes depend on the geometry of the considered structure. Figure 5(a) numerically show the eigenfrequency (the red circles) of hybrid guided mode and the resonance frequency of THG (the black five stars) for the thickness d ranging from 0.6 to 2.0 m  . It is clearly seen that as d increases, both frequencies decrease monotonically, and both results exactly agree with each other, which further prove the validity of our proposed way for THG. Figure 5(b) shows their corresponding Q factors, which undergo an increased and decreased process as d increases. To further confirm that the enhanced THG does not stem from LSP resonances, we calculated the resonance frequency of LSP-based absorption for different thickness, which is displayed by the blue circles in Fig. 5(c). For comparison, the fundamental frequency (the red circles) corresponding to THG resonances in Fig. 5(a) Fig. 5(d) shows the relationship between the THG output power and the thickness of the dielectric layer (the blue curve), as well as the ratio (the red curve) of this output power to THG output power of nonlinear metasurface without Au cylinder array. The maximum value of ratio is about 3 1.8 10  , occurring at 1.1 dm   , which is mainly attributed to the quality of the hybrid guided modes. We also carried out a study on the influence of the radius of the metal cylinders. The calculated results are shown in Fig. 6, which are similar with those in Fig. 5. For instance, as R changes, the eigenfrequency (the red circles) of hybrid guided mode and the resonance frequency (the five stars) of THG exactly agree with each other (see Fig. 6(a)). The FF of LSP-based absorption (the blue circles) and the FF corresponding to THG resonances (the red circles) also follow different variation trends, with the intersection point at 17.6 Rm   . Likewise, for BIC modes, the calculated results (not shown here) demonstrate that the THG resonance frequency shifts with the change of geometric size.

Conclusion
In conclusion, we have examined and discussed a new way to obtain enhanced nonlinear response by designing and studying a plasmonic-graphene nonlinear metasurface at terahertz regime. We have shown that due to the existence of hybrid guided modes and quasi-BIC modes with resonance frequency R f in the designed structure, a significant THG can be observed at /3 FF R ff  , and the maximum conversion efficiency can be enhanced by three orders. Different from the typical way of local field enhancement, the proposed mechanism for enhanced THG relies on strong resonance of the generated THG field with these specific modes, which greatly magnify the effective nonlinear coefficient of the whole system via photonphoton interactions. With respect to the conversion efficiency based on our proposed way, we think it should be discussed case by case, as the conversion efficiency is highly depending on the features of resonance modes and can be improved by carrying out the optimized studies. In addition, owing to the lower incident frequency, our proposed way makes the nonlinear elements further miniaturized, showing some advantages in the integrated application of photonic devices. Similar strategy can be extended to study other nonlinear effect, such as second-order harmonic generation.