Modeling of second-order nonlinear metasurfaces

We present a frequency-domain modeling technique for second-order nonlinear metasurfaces. The technique is derived from the generalized sheet transition conditions (GSTCs), which have been so far mostly used for modeling linear metasurfaces. In this work, we extend the GSTCs to include effective nonlinear polarizations. This allows retrieving the effective nonlinear susceptibilities of a given metasurface and predict its nonlinear scattering responses under arbitrary illumination conditions. We apply this modeling technique to the case of metasurfaces made of a periodic arrangement of T-shaped gold nanoparticles. For verification, several metasurfaces are fabricated and a fair agreement is found when comparing simulated data and experimental results. The proposed model may thus serve as a design platform to implement complex nonlinear metasurface based applications.


Introduction
The invention of the laser in 1960 has revolutionized physics and opened up the avenue for completely new possibilities and applications in optics [1,2]. Among many was the ability to concentrate coherent light into optical beams with extremely high energy densities, which enabled the first experimental demonstration of second-harmonic generation (SHG) from crystalline quartz [3]. This swiftly transformed nonlinear optics from a theoretical field of study into a booming applied research topic [4][5][6]. For several decades, nonlinear optical effects have been obtained using special types of crystals exhibiting particular nonlinear responses, e.g. the non-centrosymmetric crystals of barium borate (BBO) [7,8].
With the technical advances and increased availability of nanofabrication processes of the last 30 years, new methods to induce nonlinear scattering from the surfaces of nanoparticles have been explored. The first and, probably, the most extensive studies of SHG were performed with plasmonic spheres [9][10][11][12][13][14][15]. These developments established the fundamentals of the theory of SHG from plasmonic nanoparticles. Further, the study of non-spherical geometries has revealed new possibilities to effectively route and enhance the second-harmonic light radiation [16].
Interestingly, an assembly of nanoparticles can bring about even more subtle control over linear and nonlinear scattering [17,18]. Specifically, for artificial 3D structures made of a periodic or pseudo-periodic arrangement, a strong magnetic response can be obtained, opening up the possibility to achieve negative refractive index within a certain wavelength range [19]. This ultimately led to the emergence of nonlinear metamaterials [16,[20][21][22][23]. However, large metamaterials suffer from important losses and strong temporal dispersion, due to their resonant nature. Moreover, they are also difficult to fabricate due to their intricate 3D structure. This explains the strong interests in metasurfaces-the flat counterparts of metamaterials-that are easier to fabricate and much less lossy [24][25][26][27][28][29][30][31][32][33]. With their ability to strongly confine light at the nanoscale, metasurfaces provide enhanced nonlinear responses several orders of magnitude stronger than those achievable with conventional nonlinear crystals without stringent phase-matching conditions. It is therefore not surprising that metasurfaces have been considered a paradigm shift in nonlinear optics [34][35][36][37][38][39][40]. This is supported by the countless number of nonlinear metasurface based applications that have been developed in recent years, such as optical switching [41], the realization of multiple nonlinear holograms [42,43] and the generation of magnetic nonlinear effects [44][45][46].
An essential aspect of metasurface design lies in the ability to properly model such complex electromagnetic structures in order to fully exploit the large number of degrees of freedom that they provide. In the linear regime, one of the common approaches to model a metasurface is the description of its electromagnetic responses using specialized transition conditions relating the metasurface effective parameters to the fields that interacts with it [47][48][49][50][51][52][53]. The strength of such an approach is that it provides a direct connection between the desired optical transformation that the metasurface must achieve and its effective material parameters, which may themselves be connected to the shape of scattering particles. It becomes therefore much easier to find the optimal metasurface design for applications such as polarization conversion [54], angular scattering control [55] and perfect refraction or reflection [56,57].
The purpose of this work is to show how this modeling technique, which has been so far mostly used in the linear regime, may be extended to the case of nonlinear metasurfaces. Specifically, we restrict our attention to the case of second-order nonlinear metasurfaces and on the generation of second-harmonic light. The proposed modeling approach is based on a frequency-domain description of the nonlinear interactions taking place on a metasurface and bears some similarities with the effective nonlinear susceptibility approach discussed in [58,59]. Assuming an undepleted pump regime, this model provides a direct relationship between the metasurface nonlinear effective susceptibilities, the linear pump fields that excite them and the second-harmonic fields produced. This work is based on our previously published theoretical developments in [60,61], which complement and extend the discussions already available in [62][63][64][65][66][67].
This paper is organized as follows. Section 2 describes the fundamental concepts of the proposed modeling approach. Section 3 provides an example of the application of this model by showing how the effective nonlinear susceptibilities of a metasurface made of T-shaped gold scattering particles may be computed. Then, section 4 provides a comparison between theoretical prediction and experimental measurements. Finally, section 5 concludes the discussion.

General modeling concepts
A metasurface-a planar artificial structure made of a subwavelength array of electrically small scattering particles-may generally be modeled as a homogeneous sheet with effective material parameters [47,48,[48][49][50]53]. In the context of a nonlinear optical system where new frequencies may be produced, such a model is only valid if the period of the metasurface array (P) is smaller than half the smallest wavelength of interest. It follows that a nonlinear metasurface designed for second-harmonic generation (SHG) should satisfy the condition P < λ 0 /(2 r ), where λ 0 is the pump wavelength in vacuum and r is the relative permittivity of the medium on either side of the metasurface [53]. With that condition satisfied 1 , an excitation laser at frequency ω can produce very directive second-harmonic emission, as illustrated in figure 1(a). This process can be accurately modeled by replacing the metasurface by a fictitious sheet of homogeneous linear and second-order nonlinear polarization densities, as depicted in figure 1(b). The validity of the model ensures that, for identical illumination conditions, the linear and nonlinear fields scattered by the metasurface are equivalent in figures 1(a) and (b).
A metasurface model typically consists in relating the metasurface effective material parameters to the fields that interact with it. For a metasurface made of deeply subwavelength elements, its macroscopic electromagnetic responses may be satisfactorily modeled by only considering electric and magnetic dipoles expressed in terms of effective susceptibility tensors [53]. In this case, the generalized sheet transition conditions (GSTCs) may be used to model a metasurface [47,48,53]. For a metasurface lying in the xy-plane at z = 0, the time-domain GSTCs read where ΔE and ΔH are field differences defined as ΔE = E + − E − and ΔH = H + − H − with E ± and H ± being the electric and magnetic fields just below/above the metasurface (at z = 0 ± ), and P and M are the metasurface electric and magnetic polarization densities. While the GSTCs in (1) and (2) may generally be used in the time domain, it is typically more convenient to express them in the frequency domain [61]. Indeed, this simplifies the expression of the polarization densities P and M, which should otherwise be expressed as time-domain convolutions between fields and susceptibilities [5]. For a second-order nonlinear metasurface, the model may be further simplified by assuming that: (1) the metasurface only produces fields at frequencies ω and 2ω, and (2) the pump remains undepleted by the nonlinear process. These two assumptions are supported by the fact that higher-order harmonics generated by such a system are typically negligible and that the efficiency of the SHG process in a metasurface is very low compared to the linear scattering [5,67]. These simplifications allow one to greatly reduce the complexity of the modeling technique since the linear and nonlinear regimes may be split and treated individually. The case pertaining to the linear regime is already thoroughly discussed in [53], we shall therefore next concentrate our attention on the modeling approach for the second-order nonlinear regime.
For second-harmonic waves at frequency 2ω, the frequency-domain counterpart of the GSTCs in (1) and (2) is given by 2ẑ where P 2ω and M 2ω are the metasurface polarization densities evaluated at 2ω and are expressed as In these expressions, the terms P 2ω lin and M 2ω lin describe the linear interactions of waves at frequency 2ω with the metasurface. On the other hand, the terms P 2ω nl and M 2ω nl correspond to nonlinear sources that produce second-harmonic waves when excited by a pump illumination at frequency ω. In the context of second-harmonic generation, the terms P 2ω lin and M 2ω lin are defined as where χ ee , χ mm , χ em and χ me are the metasurface effective electric, magnetic and magneto-electric susceptibility tensors, respectively, and E 2ω av = [E + (2ω) + E − (2ω)]/2 and H 2ω av = [H + (2ω) + H − (2ω)]/2 are the averages of the second-harmonic electric and magnetic fields at the metasurface. In a similar fashion, the terms P 2ω nl and M 2ω nl are given by In (7) to (10), we have considered the most general case of second-order nonlinear dipolar metasurface meaning that all possible types of field interactions are taken into account. Hence the reason for the presence of so many susceptibility tensors. Note that the tensor χ eem in (9) is in fact a combination of both χ eem and χ eme since it is practically impossible to differentiate an interaction due to E ω av H ω av from one due to H ω av E ω av . For the same reason, the tensor χ mem (10) is a combination of χ mem and χ mme . In summary, the general working principle of this model is as follows: (1) a pump incident field at frequency ω impinges on the metasurface whose linear scattering response may be modeled using the ω-counterparts of (3) and (4) along with (7) and (8) following the procedure described in [53]. (2) The high intensity of the pump fields triggers the metasurface nonlinear susceptibilities, which results in the generation of the second-order nonlinear sources in (9) and (10). (3) These nonlinear sources produce second-harmonic waves, whose linear interactions with the metasurface may be modeled using (3) and (4) along with (7) and (8) in a similar way that the metasurface linear scattering is modeled in step (1).

Nonlinear susceptibility retrieval method
We shall now detail how the nonlinear metasurface model presented in section 2 can be applied to compute the effective linear and nonlinear susceptibilities of a given metasurface. To illustrate this susceptibility retrieval method, we consider the case of a metasurface made of an array of T-shaped gold particles [61,68], as illustrated in figure 1(a). To understand the electromagnetic response of such particles, we investigate their electromagnetic behavior following a mode analysis approach similar to that used in [59,69,70].
In the context of nonlinear optics, a T-shaped structure has the advantage of providing an almost independent control of the absorption of the pump incident power and the generation of the second-harmonic light [61,68]. Indeed, considering the metasurface unit cell in figure 2(a), if the pump field is y-polarized, its interactions with the metasurface will be mostly controlled by the length L y whose corresponding metallic rod acts as a receiving dipole antenna. Due to the mirror symmetry of the T-shape structure along the xz-plane and its lack of symmetry along the yz-plane, the second-harmonic light is mostly polarized along the x-axis. Therefore, the length L x may be used to tune the generation of second-harmonic light as the corresponding metallic rod acts like a source dipole antenna. The T-shaped structure can therefore be seen as a double-antenna system that can be tuned to simultaneously maximize the absorption of the pump field and maximize the generation of second-harmonic light [71,72].
This behavior is demonstrated in figure 2, which plots the linear and nonlinear fields scattered by a uniform metasurface and the corresponding real values of the charge distributions on one of its T-shaped particle. The charge distributions at the linear frequency ω plotted in figure 2(b) clearly show that the y-polarized pump field mostly interacts with the y-oriented dipolar part of the T-shaped structure, especially at wavelengths near 900 nm. On the other hand, the second-harmonic light is indeed x-polarized as verified by the nonlinear charge distribution shown in figure 2(c).
We now describe the method to obtain the effective nonlinear susceptibilities of the metasurface. The general procedure consists in solving the system of equations, formed by combining together equations (3) to (10), for the nonlinear susceptibilities in (9) and (10). The other parameters in this system of equations are the linear susceptibilities in (7) and (8), which may be obtained following the technique described in [53], and the linear and nonlinear electric and magnetic fields, which can be computed by full-wave simulations. The main difficulty is that the GSTCs in (3) and (4) contain only 4 equations, whereas the 6 nonlinear third-rank susceptibility tensors in (9) and (10) contain a total of 162 unknown components making the given system of equations largely underdetermined. It may nevertheless be solved by considering several independent illumination conditions that, when combined together, lead to a fully determined system [61,73]. However, depending on the applications, computing all nonlinear susceptibility components may not necessarily be needed and often only a few of these components turn out to be really relevant.
For instance, the nonlinear metasurface discussed above is only illuminated with a y-polarized pump field and produces only x-polarized second-harmonic light. Moreover, all waves interacting with this metasurface are plane waves that propagate normally with respect to the metasurface plane. This implies that only a small subset of susceptibilities really plays a role in the metasurface scattering process. In our case, the linear fields E  (4) is reduced from 4 to 2. Therefore, we now have a system of 2 equations in 6 unknowns that may be solved with the 3 illumination conditions whose electric fields are defined by 3 The 6 nonlinear susceptibilities mentioned above are now obtained by substituting in (3) to (10) the simulated linear and nonlinear fields computed by sequentially illuminating the metasurface with the fields in (11) to (13). For this purpose, we used a homemade surface integral equation method that was specifically implemented for both linear and nonlinear periodic simulations [74][75][76]. The resulting normalized susceptibilities are plotted in figure 3(a). For such a simple structure, we see that the metasurface response is dominated by χ xyy eee with a minor contribution from χ yyx mem , while the 4 remaining susceptibilities are not plotted as their contribution is negligible. This is easily understandable as the T-shaped structure essentially behaves as a nonlinear x-oriented electric dipole excited by a linear y-oriented electric dipole, as expected from the results shown in figure 2. The presence of χ yyx mem is most likely a consequence of retardation effects due to the thickness of the particle that results in the generation of a nonlinear magnetic response.
As a comparison, we also consider the case of a metasurface made of double T-shaped structures with inverted T's, as shown in the inset of figure 3(b). The susceptibility retrieval method for this metasurface is identical to the one used for the structure in figure 3(a) due to the similarities in terms of shape and orientations of the particles. However, the retrieved susceptibilities are this time quite different than those of the single T structure. Indeed, the metasurface response is now dominated by the susceptibility χ yyy mee indicating that the double T-shaped structure mostly behaves as a nonlinear y-oriented magnetic dipole  resulting from a circulating nonlinear current distribution between the top and bottom T's excited by the two linear y-oriented dipoles retarded in phase with respect to each other. The retrieved susceptibilities presented in figure 3 clearly show that the proposed modeling technique yields results that are consistent with the expected electromagnetic responses of the metasurface. Further more in-depth demonstrations of the validity of this modeling technique are available in [53,61].
The results plotted in figure 3 demonstrate how electric and magnetic nonlinear responses may be achieved and controlled. This is of particular importance since the local tuning of the electric and magnetic responses of a metasurface is one of the main mechanisms to control the phase of the nonlinear scattered light, thus allowing applications such as nonlinear refraction, collimation and hologram generation. Moreover, a balance between the nonlinear electric and magnetic responses can be used for ultra directive nonlinear light generation by using the Kerker effect [61,77].

Experimental results
In order to verify the validity of the nonlinear GSTCs method, we have fabricated several different metasurfaces and compared their measured nonlinear scattering to the full-wave simulations following the procedure described above.
Each metasurface was fabricated by first using an evaporator to deposit a uniform 40 nm layer of gold and a 20 nm sacrificial layer of chromium on top of a float glass wafer. Then, 170 nm of HSQ was spin coated which, after Ebeam exposition, was developed in TMAH. An ion beam etcher was then used to etch through the chromium and gold layers. Finally, wet etching was used to get rid of the sacrificial layer of chromium, which removed the remaining HSQ. A total of 12 different metasurfaces were fabricated using this process and their corresponding dimensions are provided in table 1. To show the quality of the fabrication process, figure 4 presents small areas corresponding to 4 examples of the fabricated metasurfaces.
To measure these metasurfaces, we used a y-polarized high-power pulsed laser beam with a central wavelength of λ 0 = 800 nm. The laser beam was focalized on the samples with a high-NA oil immersion objective lens. For each metasurface, the nonlinear x-polarized backward propagating waves were measured with a spectrometer. These measurements were repeated 10 times to investigate their variance. The results are plotted in figure 5, where the boxplots correspond to the measured samples and the black circles to the corresponding normalized full-wave simulated data.
The comparison presented in figure 5 shows that the overall trend describing the change in the normalized nonlinear scattering intensity between the different samples obtained via full-wave simulations  is in a fair agreement with the measured data. The main reason that the match between theory and experiment is not better may be attributed for the most part to geometrical variations between the fabricated samples and the simulated structures, which strongly affect the metasurface response [79], especially the nonlinear one.

Conclusion
We have presented an overview of a frequency-domain modeling technique for second-harmonic nonlinear metasurfaces. The model is designed around the assumptions that only linear and second-harmonic fields are scattered by the metasurface and that the pump may be assumed undepleted by the nonlinear interactions. The validity of the model was demonstrated using full-wave nonlinear simulations and the ability to practically implement such nonlinear responses was investigated with a series of experimental measurements. The agreement that was found between simulation and measurement data further confirms the capability to finely control second-harmonic generation with metasurfaces and the proposed modeling technique may serve as an excellent theoretical platform to help achieve advanced design goals.