Second harmonic generation under doubly resonant lattice plasmon excitation

Second harmonic generation is enhanced at the surface lattice resonance in plasmonic nanoparticle arrays. We carried out a parametric investigation on two-dimensional lattices composed of gold nanobars where the centrosymmetry is broken at oblique incidence. We study the influence of the periodicity, the incidence angle and the direction of the linear input polarization on the second harmonic generation. Excitation of the surface lattice resonance either at the fundamental or second harmonic wavelength, achieved by varying the incidence angle, enhance the conversion efficiency. As a special case, we demonstrate that both the wavelengths can be simultaneously in resonance for a specific period of the lattice. In this double resonant case, maximum second harmonic power is achieved.


I. INTRODUCTION
Thanks to the fast development in manufacturing of subwavelength structures in a controlled and reliable way, metamaterials allow to manipulate the light-matter interaction beyond the possibilities of conventinal materials. Plasmonic metasurfaces [1] are bidimensional structures consisting of periodically arranged metallic nanoparticles. Plasmons are the collective oscillation of conduction band electrons in metals, and they have been largely studied due to intriguing properties, such as dramatic subwavelength localization and extreme electromagnetic field enhancement [2]. The optical properties depend on the electromagnetic mode supported by the elemental nanoparticles (localized surface plasmon resonance, LSPR) and on the long-range coupling determined by the periodic arrangement (surface lattice resonance, SLR) [3]. SLRs can emerge when the scattered field of a single nanoparticle is in phase with the emitted fields from the surrounding nanoparticles, exciting a collective plasmonic wave propagating along the lattice. Accordingly, the conditions for SLR are closely related to the presence of a diffraction mode propagating along the lattice, also known as Rayleigh anomaly (RA) [4,5]. The spectral position of the SLR depends on the nanoparticle periodicity and the surrounding refractive index, but can be tuned with the angle of incidence. At the SLR the oscillation damping is partially compensated [3,6], thus tuning the interplay between LSPR and SLR enables * These authors equally contributed to this work.
to design ultra narrow resonances [7][8][9]. The associated strong field enhancement indeed boosts nonlinear effects. In 1985 Couatz et al. showed that, as the incidence angle is changed, SHG undergoes strong enhancement in metallic mono-dimensional gratings when the impinging beam is resonant with the nonlocal surface plasmon [10], what is today known as SLR. However, in the same structure the SHG is also enhanced when the SLR condition is fullfield at the second harmonic frequency [11]. Similarly, for two-dimensional (2D) nanoparticle arrays the SHG is enhanced for the SLR at the fundamental frequency [12,13] and at the second harmonic [14]. The primary role of collective response yields surprising and counter-intuitive phenomena, such as larger SHG for less dense metallic gratings [15,16] or the strong dependence on the structure of the unit cell [17]. Even beyond the linear regime, plasmonic metasurfaces provide a versatile way to shape light [18]. The conversion process can be further improved by utilizing multiple-resonant nanoparticles, where the LSPR enhancement occurs simultaneously at the pump and signal wavelengths [19][20][21][22][23][24][25]. Recently, the relevance of double-resonance in the presence of SLRs along orthogonal directions has been stressed out theoretically [26] and was experimentally demonstrated [27]. In this paper, we perform a parametric investigation to study the influence of the SLR on the SHG, exploring various double resonant conditions, where the SLR of the fundamental and second harmonic are simultaneously excited. Therefore, multiple 2D rectangular arrays composed of centrosymmetric gold nanobars on fused silica are fabricated, whereby just the lattice period in one direction is varied while the LSPR value is kept constant. We measured the SHG from these metasurfaces at different angles of incidence and linear polarization directions. The required non-centrosymmetry for SHG is provided by the left-right symmetry breaking associated with a tilted incident wavefront. Relatively high average output powers are achieved, thus allowing the experimental investigation of the nonlinear effects using general purpose and cheap equipment such as CMOS cameras and com- pact spectrometers. Finally, we model our experimental results using the nonlinear inverse scattering method [28], which indeed connects the SHG to the the linear response of the structure at the two involved wavelengths [29].

II. SAMPLE FABRICATION AND LINEAR CHARACTERIZATION
We study several samples, where gold nanobars are arranged into 2D lattices on a 1 mm thick fused silica substrate. The metasurface template is shown in Fig. 1a. In all the fabricated samples, the shape of the bars working as nanoantennas is preserved. The elemental antenna is 400 nm long (w x ), 300 nm wide (w y ), and 50 nm high. The period in the y-direction P y is fixed to 500 nm. We vary solely the period P x along the x-direction, the latter spanning the interval from 520 to 1200 nm with steps of 20 nm. The structured area of each lattice is 3x3 mm 2 . The wide area nanopatterning is accomplished using electron-beam lithography and standard metal liftoff technique (see Appendix VI A). The scanning electron microscope (SEM) image in Fig. 1b shows that our gold bars possess a large degree of uniformity, at the same time featuring very smooth edges on the nanometric scale.
To characterize the linear response of the fabricated metasurfaces, a supercontinuum laser is used as broadband light source in the setup (see Supplement 1, Fig. S1), featuring a beneficial larger spatial coherence [9,30]. The light is weakly focused on the metasurface to a beam diameter of approx. 300 µm with an achromatic 400 mm focal length lens. In all the experiments presented here the sample is rotated around the y-axis at the metasurface, changing the optical angle of incidence θ AOI on the structure. Thus, we are exploiting the angular dependence of the SLR only along the xdirection. Light propagates towards positive z-direction, and impinges the metasurface first before entering into the glass substrate.
The numerically computed LSPR for one isolated nanobar is at λ =1276 nm for light polarized along the x-direction (TM or p-polarization) and along the ydirection (TE or s-polarization) at λ =978 nm (see Supplement 1, Fig. S2). The corresponding bandwidth is very broad, approx. 700 nm. The experimental linear transmission spectrum as a function of the angle of incidence reveals the location of the plasmon resonances [6], exemplary shown for one lattice in Fig. 2a. The broad transmission minimum at about λ =1150 nm (for TE at λ =900 nm, not shown here) emerges from the hybridization between LSPR and the P y lattice mode [31]. Its position does not change with P x (see Supplement 1, Fig. S3). The narrow SLRs are visible as transmission anomalies [32,33], which closely follow the wavelength and the angle-depending RA conditions, represented by the superimposed white lines. The incidence angle to achieve the RA condition either for a reflection grating (light still propagates in air) or transmission grating (light propagates in glass substrate) is given by [6]: sin θ glass = ± n glass n air − mλ n air P x , where P x is the period along x, m the diffraction order, λ the vacuum wavelength, and n air and n glass are the wavelength-dependent refractive indices of the air and fused silica medium, respectively. At the RA, the diffraction order m is diffracted by ±90°and travels along the metasurface (see pictorial sketch in Fig. 2b). The asymmetric index environment for the metasurface introduces two families of RA curves, thus introducing an additional degree of freedom to achieve the double resonance. However, this comes at the price of a reduction of the SLR Q-factor [34]. Simulated spectra stemming from RCWA (Rigorous Coupled Wave Analysis) match very well with the experimental results (see Supplement 1, Fig. S4). The good agreement shows that small technological imperfections (surface roughness, particle shape deformation, edge and corner rounding, etc.) have a negligible effect on the long distance interaction [35] and that enough nanoparticles couple to each other to be approximated with an infinite array [36].
Beyond the RA condition stemming from the mo- mentum conservation, nonlinear effects are strongly dependent on the optical field distribution. We computed the optical field at the metasurface by using a commercial FDTD tool (see Appendix VI B). At normal incidence θ AOI = 0°, TM and TE polarized light excite an electric dipole in the nanobar along the x and the y-direction, respectively. A rotation of the metasurface around the y-axis leads to an asymmetric field distribution along x (Fig. 3). In the TE case, the dipole changes into an electric quadrupole distribution. Note that the field enhancement is stronger for TE than for TM. We confirmed the presence of a nonlocal plasmon oscillation with the field extending well beyond the unit cell [31] by measuring the power transported along x-direction on the metasurface plane depending on θ AOI . Peaks and discontinuities in the power are observed whenever a SLR is excited. We now focus on the two wavelengths corresponding to the fundamental (or pump, λ F F =1032 nm) and its second harmonic (λ SH =516 nm) in the frequency conversion experiment. Figure 4a-d show the simulated linear transmission at λ F F (top row) and at λ SH (second row) versus θ AOI and P x , for TM (left column) and TE (right column) input polarizations. Dips, peaks and jumps are mostly visible in proximity of the RA conditions, confirming that they originate from the SLR. For TE polarized light, the SLR in air is prevented by the reflection of the substrate [34]. Due to the smaller fill factor, the absolute transmission increases for longer periods. The average transmission at λ F F is lower than at λ SH due to the shorter spectral distance from the LSPR.

A. Second harmonic generation
For the nonlinear investigations we used an ultrafast laser delivering 200 fs (FWHM) short pulses at a wavelength of 1032 nm with a repetition rate of 200 kHz and a variable average power. The linearly polarized beam is focused onto the structured area to a 250 µm beam diameter with a 400 mm focal length lens. The metasurface leads to a frequency conversion of the impinging light. We investigate the SHG light which propagates in the direction of the incidence light (0 th transmission order). We experimentally verified that SHG from the pure glass surface [37] is not detectable in the range of powers we employ (more information see Appendix VI C and Supplement 1). While the centrosymmetry of gold is broken at the nanoparticle surface, the SHG is still prohibited due to the spatial symmetry. However, oblique incidence angles break the spatial symmetry and allow SHG from the structure [38], see Fig. 3. The angular dependency of the SHG signal for different periods P x is shown in Fig. 4e-f. By a direct comparison with the linear response ( Fig. 4a-d), the strict connection between linear and nonlinear regime becomes immediately clear [14]. At the RAs, a fast change versus the incident angle is observed in the SHG signal, including even the 2 nd diffraction order. It is obvious that the SLRs have a significant impact on the frequency conversion. Pronounced SLR related anomalies in the linear transmission also cause higher SHG efficiency. The strong linear dependence on the input polarization is transferred to the SHG. In fact, the TE case is mostly sensitive to SLR in glass (dashed red line), whereas for TM inputs SLRs are excited on both surfaces. While lattice resonances of λ F F are more relevant in determining the SHG, surprisingly even the barely visible λ SH SLR in the linear regime has an influence on the frequency conversion strength for both input polarizations.
Especially interesting are the multiple intersections between the λ F F and λ SH SLRs, corresponding to the double resonance case. To see the behavior approaching the double resonance in more detail, we plotted in Fig. 5 the average SHG power versus the incidence angle for the period region of interest. For TM, the most prominent peak is at (P x =760 nm, θ AOI =20.5°) see Fig. 4e and Fig. 5a, corresponding to the double resonance. Close to this condition, the linear diffracted light in air travels for both wavelengths along the nanoparticle array, schematically shown in Fig. 2b. For slightly detuned periods, the two lattice resonances are clearly visible.The region from 520 nm ≤ P x ≤ 620 nm, θ AOI < 32° (Fig. 4e) shows a SHG enhancement with a peak at the SLR for λ F F in glass. For increasing periods the SHG peak undergoes a sudden increase reaching a local maximum in P x = 620 nm, θ AOI = 11.5° (Fig. 5b), correspond-ing to the double resonance with λ SH in air and λ F F in glass. The other double resonance of these two SLR at (P x = 1140 nm, θ AOI = 33.5°) leads also to a small signal peak (not shown). On the other hand, the double resonance with the λ SH SLR in glass as well as the simultaneous excitation of both λ F F SLRs does not show an enhancement. In the TE case for 520 nm < P x < 700 nm (Fig. 4f), the SHG rises for increasing periods, reaching a global maximum in P x = 620 nm, θ AOI = 11.5° (  Fig. 5c). The effect of the double SLR in glass manifests as a fast decay of the peak for longer periods. In the TE case there is another enhancement region P x ≥ 1080 nm (Fig. 4f), where two very close SLR peaks merge into a single peak for several periods after crossing with a third resonance (second order at λ F F in glass). Other double resonances show minor effects: near the fundamental SLR (P x = 860 nm, θ AOI = 13.5°), a relative SHG maximum is observed. The TE polarized light provides higher conversion efficiencies than TM polarized light, related to the stronger field enhancement. The highest conversion efficiency is achieved for TE polarized light at the double resonance. On the other side, for TM inputs the strongest signal is observed at large angles and short periods along the SLR in air ( Fig. 4e and 5b), where the emitted SHG is slightly lower than in the double resonant case. We ascribe the high efficiency to the larger E-field asymmetry and the longer coupling length between the incident beam and the metasurface. We observed a saturation at large incidence angles θ AOI >45°. At normal incidence, the signal is usually below the detectable limit, except when the SLR is close to normal incidence, see Fig. 4f at P x ≈ 700 nm. This is probably caused by a weak quadrupolar mode [13] or small fabrication imperfections [39] enhanced by the resonance. An analyzer at the output is used to determine the polarization of the SHG. Independently from the incident polarization, the SHG component along the x-direction (TM output) is always dominant. Indeed, the SHG polarization is determined by the spatial symmetry breaking, which occurs along the x-direction due to the tilted input wavefront [38], see Fig. 3.

B. SHG characterization
To verify that the origin of the detected signal is a second order nonlinear process, the emitted spectrum is analyzed via a compact back-thinned CCD spectrometer. In the bandpass filter limited range from 500 nm to 789 nm, only a narrow peak is observed around the second harmonic, as shown in Fig. 6a. Thus, two-photon luminescence is negligible in our case [40]. The center wavelength of the SHG (516.6 nm) is slightly red-shifted with respect to the half excitation wavelength, probably due to the asymmetry between forward and backward scattered light [41]. Furthermore, the central wavelength of the SHG does not change noticeably with the incidence angle or the lattice period P x (not shown here).  A typical angular dependence of the SHG spectrum is shown in Fig. 6b. Sharp discontinuities are clearly visible at the indicated RAs, their positions changing with both the angle θ AOI and the emitted wavelength λ SH . Noticeably, a spectrally-resolved measurement permits to reveal the sharpness of the transition. The sharp features get indeed dimmed when using a detector uncapable of distinguishing the different wavelengths.
The origin of the SHG signal is also verified by measuring the input-output power, as a double-log plot in Fig. 7a for the case of TM polarization and double resonance. The output power ramps up with the average input power, reaching a threshold around 180 mW (corresponding to a peak intensity of ≈15 GW cm −2 ) where the sample gets permanently damaged. The maximum average output power is 190 pW, corresponding to a maximum SHG conversion efficiency of 1.06 × 10 −9 . Taking the metasurface thickness into account, the maximum second order nonlinearity is around 1.0 pm V −1 (for more details, see Supplement 1), which is of the same order of magnitude with respect to similar plasmonic metasurfaces [11,27,42]. Polynomial best-fitting finds a linear slope in the double log scale of 2.35. The exact value of the slope depends on the incidence angle and on the period. Slightly steeper slopes than 2 have been already reported for nonlinear processes in metallic nanostruc-tures [27,43,44]. Such a deviation from a pure square dependence can originate from the spatial and temporal nonlocality associated with the excitation of electrons in an energetic band [45]. We finally investigated what happens when the sample is damaged. The recorded images in Fig. 7b reveal that above 180 mW average input power a black area is formed in the center of the SHG spot, where the nano-structure is permanently modified. Increasing the excitation power above this threshold leads to a wider modified area but also to a larger SHG spot. The output signal thus gradually decreases above the permanent damage threshold. Normalized SHG spectrum for Px=800 nm at different angles of incidence θAOI ; the superposed solid lines show the RA for the fundamental (red color) and the second harmonic (green color).

IV. INTERPRETATION OF THE RESULTS
The simplest way to describe SHG from metallic metasurfaces is the nonlinear scattering technique first described in Ref. [28]. Using the Lorentz reciprocity theorem, the electric field at the second harmonic E SH in the far field is given by the following surface integral extended to all the gold interfaces [29] is the nonlinear polarization at 2ω, supposed to be nonvanishing only at the six interfaces enclosing each nanobar. The factor γ is equal to cos θ AOI forl =x (TM output) and equal to unity forl =ŷ (TE output). Given we are interested in the transmitted SHG, the quantity E PW is the electric field generated on the sample by a plane wave which propagates from the observation point (coincident with the direction in the far field corresponding to the 0 th order) to the metasurface at 2ω with polarization parallel tol. The nonlinear polarization P NL depends on the linear coupling of the fundamental beam propagating in the z-direction and impinging on the sample, whereas E PW depends on the linear coupling at 2ω for waves propagating towards negative z. According to Equation (3), the overlap integral between these two quantities determines the amount of generated SH (for more information, see the Supplement 1). SHG is then maximized in correspondence to the linear resonances at ω and 2ω, in agreement with the coupled dipole approximation [14], with resonance at the fundamental intervening two times in the case of Type I SHG. Equation (3) permits to understand why in the experiments the maxima of SHG are observed always at some RA, but not all the RAs correspond to a peak of nonlinear emission (e.g. SLR in air for TE inputs) [19]. The SHG indeed depends also on the spatial overlap between the two fields in the linear regime, the overlap being moreover modulated by the nonlinear coefficient χ (2) of the nonlinear lattice. When comparing different lattices, maximum SHG is obtained when double resonance occurs, i.e., for periods where fundamental and second harmonic resonances are observed for the same incidence angle θ AOI .
A more physical interpretation can be given in correspondence to the lattice resonances. Owing to the momentum conservation, the pump beam efficiently excites a propagating delocalized surface plasmon at ω. The second-order nonlinearity of the metasurface then generates a plasmon at 2ω. The last step is the coupling of the plasmon at 2ω with the free space: this is maximized when a RA in the linear regime is observed at 2ω. We finally compare the predictions of the nonlinear inverse scattering model with the measured data. Specifically, in Fig. 8 the theoretical (red lines) and experimental (black lines) normalized SHG in correspondence to the calculated RA in air (Eq. 1) for |m| = 1 and λ F F = 1032 nm. For TM polarized pump (Fig. 8a), the experiment and simulation have a peak near the double resonant condition P x = 760 nm. For lower periods the SHG power decreases, and the experiment features a lower depression than predicted. The signal reaches another maximum at P x = 600 nm. For even shorter periods and higher incident angles the signal suddenly drops because the calculated RA deviates from the SHG peak. The experimental results in Fig. 5(b) show that the maximum signal stays constant for higher incidence angles. In the TE case plotted in panel Fig. 8(b), the experiment shows a broad peak due to the interaction between the double resonances at P x = 760 nm and the crossing with the SLR of the pump wavelength in glass at P x = 840 nm. Interestingly, the latter does not cause any visible effect in the TM case. The simulations just predict a narrower peak at P x = 840 nm, not the double resonance. Limitation of the simulation is the prediction of the SHG enhancement at the second harmonic SLR. Finally, we stress out that the numerical simulations also confirm that i) the main contribution of SHG is TM polarized, no matter what the polarization state of the pump is; ii) The SHG is stronger for TE than for TM polarization input.

V. CONCLUSION
We experimentally demonstrated the role of double surface lattice resonance in maximizing the SHG in plasmonic metasufaces. Our parametric study on centrosymmetric 2D nanobar arrays confirms [14] that local maxima of the SHG occur in correspondence to the SLR either at ω or 2ω, thus leading to a further enhancement of frequency conversion when both pump and emitted SLR get excited simultaneously. Such a result has been predicted based upon the discrete dipole approximations [26] and recently experimentally confirmed by tilting the sample along two orthogonal directions [27]. Here we created the double resonance by tilting the meta-surface with a specific period in just one direction. The close connection between the linear and the nonlinear response of the metasurface has been theoretically confirmed using the nonlinear inverse scattering approach, providing a generalization of the Miller's rule to the SLR case [29]. Finally, due to the high quality of our samples and the fact that we are exploiting the double resonance, we achieved a maximum SHG conversion efficiency of 1.06 × 10 −9 . Such a comparatively large efficiency, in combination with high average excitation power, allowed us the direct measurement of both the spatial profile and spectrum of the generated second harmonic.
Straightforward extensions of our work include nanoantennas of different shape, with specific attention to nano-gapped resonators. A deeper theoretical treatment will be pursued in order to fully understand the nonlinear coupling between the nonlocal surface plasmons propagating across the metasurface, for example concentrating on the role played by the multipolar response of the elemental antenna.

A. Fabrication
To fabricate 2D gold nanobar arrays, we have used a standard metal lift-off process. First, we took the clean fused silica substrate and pattern the nanostructures using electron-beam lithography on a bi-layer of positive tone resists (300 nm ARP617.06 and 100 nm ARP6200.04) using a 10 nm gold layer as a conductive surface on top. The sample was then developed in AR600-546 for 30 s to remove the ARP6200.04 resist layer. The remaing ARP617.06 resist was then developed in MIBK:IPA (1:1) for 80 s. A 50 nm Au layer was then deposited using an electron beam evaporator at 1 nm s −1 deposition rate at a pressure lower than 2 × 10 −5 mbar, using a 3 nm adhesion layer made of titanium. After the deposition process, a lift-off based on acetone and isopropanol is performed to realize the gold nanobar arrays on the fused silica substrate.

B. Numerical simulations
Linear optical transmission spectra of a 2D array of nanobars were simulated using a commercial Finite Difference Time Domain (FDTD) solver by Ansys Lumerical and Rigorous Coupled Wave Analysis (RCWA) [46]. The gold optical constants used for simulations were taken from Ref. [47]. Transmission computed from FDTD and RCWA results are in very good agreement with each other, regardless of the geometric parameters of the structure. Due to the larger computation time required by the FDTD simulations, RCWA is a better choice for investigating the parametric dependence of the optical response of the nanostructures. In RCWA-based simulations, 225 spatial harmonics were enough to obtain a good match with the experimental results for the 2D grating. After the simulation, the reflection from the backside of the substrate was taken into account in order to be comparable with the experiment. The map of the electric fields on the surface of the nanostructures is simulated using FDTD solver. FDTD solutions are also used in the computation of the overlap integral providing the SHG according to the nonlinear inverse scattering technique. In all FDTD simulations, Bloch boundary conditions were used in x and y-directions, whereas PML (Perfectly Matched Layer) boundary conditions were applied along z. The nanobar object was meshed 5 nm on the transverse plane xy-directions and 2 nm on the zdirection using mesh override region.

C. Nonlinear detection
The optical intensity at the sample plane is imaged onto a standard CMOS camera to investigate the emitted SHG power. Bandpass filters limit the detectable spectral range from 500 nm to 540 nm. The signal detected on the camera is calibrated versus the SHG of a nonlinear crystal, losses in the system were taken into account. The linear behavior of the camera in the range of interest is verified. The validity of the camera measurement at low signal is also cross-checked with a photo-multiplier tube. For comparison we verified that the input-power dependency of the SHG from a bulk nonlinear crystal (BBO) show a quadratic behavior (slope 2.04). We use a camera because several SHG output beams, parallel to the original 0 th transmission order, can arise at high incidence angles and large periods, due to multiple reflections inside the glass substrate of the sample. The camera allows an easy distinction of these few millime-ters separated replicas. In addition, the beam is laterally displaced when the whole substrate is rotated around the y-axis. Finally, for better visibility we only showed the response for positive angles, but we scanned from −70°t o 70°to ensure that the SHG response under the employed illumination conditions is equal for negative and positive angles within our experimental accuracy. Further details regarding the experiment can be found in the Supplement 1.

ACKNOWLEDGMENTS
The authors thank the Deutsche Forschungsgemeinschaft (DFG) for funding the project (project number 398816777) within the framework of the CRC 1375 NOA. We also acknowledge the valuable support of Werner Rockstroh, Natali Sergeev, Detlef Schelle, Holger Schmidt and Daniel Voigt in the fabrication of the nano-structured samples.
Supporting Information:

Second harmonic generation under doubly resonant lattice plasmon excitation
In this supplemental material we first provide more details on the linear response of the structure, including the setup and the comparison between experiments and numerical simulations. The second part contains more details about the setup we used for the measurement of the SHG, including how we calibrated our camera. In the last part we provide all the step in deriving the generalization of the nonlinear inverse scattering to the case of surface lattice resonances. Finally, the method to find the effective second-order nonlinear tensor is presented.

I. LINEAR TRANSMISSION SPECTRUM
A. Linear transmission spectrum setup Figure S1 shows the setup for the experimental linear transmission measurements. The beam from a supercontinuum laser (SuperK Extreme EXW-1 by NKT) is polarized, expanded to approx. 4 mm and slightly focused onto the sample with a 400 mm lens, similarly to the SHG characterization setup in Fig. S5. The transmitted 0 th order is coupled into a multimode fiber guiding it to a scanning spectrometer (AQ6370D by Yokogawa). Changing the incidence angle leads to a plan-parallel beam displacement by the 1 mm thick fused silica substrate, which also changes the coupling into the fiber. This is taken into account by the reference measurements were performed on the unstructured area of the sample.  Figure S2 shows the simulated extinction cross-section versus the wavelength at normal incidence of a single isolated nanobar. A peak is visible at the LSPR.  Figure S3 shows the simulated linear transmission spectra for the shortest and longest fabricated periods as computed from the RCWA. The broad transmission valley is independent from the periodicity, but changes with the polarization direction due to the nanobar shape. The absolute transmission decreases for smaller periods due to the larger fill-factor.  Figure S4 shows the experimental (left panel) and theoretical (right panel) transmission spectra versus the incidence angle for a period of P x = 800 nm. The simulated spectra stemming from Rigorous coupled-wave analysis (RCWA) match very well with the experimental results, regardless of P x .   Figure S5 shows the SHG characterization setup. An ultrafast laser (PHAROS-SP by Light Conversion) delivers 200 fs (FWHM) short pulses at a wavelength of 1032 nm with a repetition rate of 200 kHz and an average power of 150 mW. The input power and linear polarization direction are adjusted via a half wave plate and a polarizer. The beam is focused onto the structured area to a beam diameter of 250 µm using a plano-convex 400 mm focal length lens. Before the sample, parasitic optical radiation at shorter wavelengths is removed by a long pass filter (FELH0900 by Thorlabs). During the measurements the sample is automatically rotated along the y-direction at the focus position. We investigate the straight transmitted SHG light depending on the angle of incidence. Behind the sample, the illumination light is seperated from the weak SH light with a short pass filter (FESH0800 by Thorlabs), which is highly transmissive (>95 %) from 500 nm to 789 nm. In addition, a 520 nm band pass filter (FBH520-40 by Thorlabs) with a FWHM of 40 nm can be used to filter out ambient light. To image the sample with its emitting signal, a 100 mm bi-convex lens is placed at a distance of 370 mm behind the sample. An MgF 2 Rochon prism (RPM10 by Thorlabs) analyzer can be placed in front of the imaging lens to determine the polarization direction of the SHG. The sample is imaged onto a CMOS camera (UI-1490SE-M-GL by iDS). Once the background level is compensated for, the sum of the signal measured on each pixel overlapping with the observed SHG spot divided by the used exposure time is proportional to the average optical power. The calibration factor of the camera to determine the impinging SHG power was found in the following manner. For the power calibration we inserted a BBO crystal in the setup instead of a sample. The output power was then high enough to be measured with a calibrated photodiode (S121C by Thorlabs) at the position of the camera. In the following step several calibrated neutral filters were inserted to reduce the power to the range of the second harmonic signal emitted from the metasurface. Knowing the output power of the crystal and the optical density of the filters, it is then possible to connect the observed signal on the camera with the incident power at 2ω. The ratio between the calculated input power and this camera value is the calibration factor. From the captured images, the digital signal was integrated over the observed SHG spot, divided by the exposure time used and multiplied with the calibration value to get the optical power. In the procedure, the linear behaviour of the camera was ensured for different exposure times and digital signal levels. On the other side, to investigate the SHG spectrum, the output light is directly coupled into a compact spectrometer (Ocean HDX by Ocean Insight) with a back-thinned CCD. The image of spectrum shown in the main article are not subject to any processing, except for the subtraction of the constant background, the latter value being taken in region where no signal is measured. The error bars of Fig. 7(a) take into account the measured small laser fluctuations and their effect on the SHG, the measured background noise and the resolution of the measurement devices. We find that the main error source is the thermal noise from the pump power detector. The 3-10 % error of the SHG power is scarcely visible due to the logarithmic scale. On the other hand, the error bars in Fig. 8 correspond to the difference between the SHG signal between positive and negative angles of incidence.

II. SECOND HARMONIC GENERATION EXPERIMENT
SHG from the pure fused silica surface is detectable at high angles of incidence for average powers above the damage threshold of the metasurface. The maximum conversion efficiency is 3 order of magnitude lower than from the metasurface, therefore the susceptibility of the glass surface is neglected.

Nonlinear inverse scattering technique
Let us consider two electrical currents J 1 and J 2 , sources respectively of two electric fields E 1 and E 2 . When both dielectric and magnetic permittivity tensors are symmetric, the Lorentz reciprocity theorem in the harmonic regime (fields proportional to e iω t ) states (S1) We want to apply Equation S1 to predict the SH emission from a periodic metallic metasurface. Following [29], we consider J 1 to be originated from the nonlinear polarization at ω = 2ω. We thus set J 1 = 2iωP NL . The other source J 2 is instead a localized dipole J 2 (r ) = J 0 δ(r − r)ĵ, used to sample in each spatial point r the j-th component of the field E 1 emitted from J 1 . Equation S1 then provides [28] E emission (r, 2ω) ·ĵ = 2iω where we set E 1 = E emission and E 2 = E dipole . The field E emission is the second harmonic emitted from the metasurface, whereas E dipole is the field on the metasurface induced by the elemental dipole J 2 placed in r. As well known, the presence of nonlinearity breaks the reciprocity of the system [48], seemingly invalidating Equations S1 and S2. The point is that the derivation presented above assumes a fixed nonlinear polarization, thus treating it as a known forcing term. This approach holds valid in the non-depleted regime, where the number of photons in the fundamental beam is not significantly reduced by the SHG. Nonetheless, the nonlinear inverse scattering technique has been also generalized to more general conditions, including losses and saturation [49]. The last ingredient is then finding the correct form to model the nonlinear polarization. In our case the main contribution comes from the surface nonlinearity at the gold interfaces. Accounting for the isotropic response of gold, the associated elements of the nonlinear tensors are χ  [50], where the subscripts t and n refer to the components of the electric field tangential and normal to the interface, respectively. Equation S2 is already apt to describe the second harmonic emitted from our metasurface. But given the periodicity of our structure, it is more convenient from a numerical point of view to compute the electric field inside the integral by considering a plane wave excitation [29]. The electromagnetic field radiated from a dipole tends to a plane wave when the observation point is far enough from the source. Thus, we suppose the point r to be far enough from the metasurface, corresponding to compute E emission in the far field. The Green function G in free space for the magnetic potential vector A is the spherical wave G(r − r ) = µ 4π|r−r | e −ik|r−r | , where k is the wavevector in air. In the Lorenz gauge the electric field can be derived from the magnetic potential vector A using E = −2iωA − i/ (2ωµ ) ∇∇ · A, in turn providing in the far field where r is the distance from the source and θ is the angle between the dipole axis and the emission direction. Next, we need to find the relationship between the amplitudes of the dipole J 0 and the impinging plane wave a 0 . Such connection depends on the relative orientation between the dipole and the interface, i.e., the emitted polarizationĵ. For TE outputs (ĵ parallel toŷ) the dipole stays parallel to the interface regardless of the incidence angle θ AOI . Thus in Equation S3 it is always θ = π/2, in turn providing For TM outputs it isĵ =x = cos θ AOIt + sin θ AOIn , wheret andn are the unit vectors parallel and normal to the interface. The component of the dipole alongt behaves like the TE component and follows Equation S4 as well. The normal component of the dipole is equivalent to θ = π/2, thus radiative coupling with the metasurface is forbidden.
Remembering that θ = π/2 − θ AOI and that the effective dipole amplitude is J 0 cos θ AOI we find Equation S5 can also be found more directly from Equation S3 by considering the metasurface being rotated around the constant linear dipole of amplitude J 0 and accounting for the relation between θ and θ AOI . All the above polarization-dependent results can be gathered together defining the diagonal matrix in the laboratory framework xyz Substituting back to Equation S2 we find The new field E planewave is the field on the surface of the metasurface when it is illuminated with a plane wave of amplitude a 0 , frequency 2ω and wavevector k = cos θ AOIẑ + sin θ AOIx . With respect to the original formula derived in [29], we have an additional anisotropic factor modelled by the matrix Υ accounting for the skewed incidence angle on the sample.

Green function approach
The nonlinear inverse scattering technique described in the previous paragraph can be interpreted directly in terms of antenna theory. The crucial point is that we are assuming the nonlinear polarization P NL to be fixed and determined by the the interaction of the fundamental beam with the metasurface in the linear regime. For a distribution J of currents, the emitted electric field is the superposition integral between J and the dyadic Green function Γ = −iω I + k −2 ∇∇ G int [51], where G int is the Green function accounting for the presence of the interface glass-air and the metallic nanostructures. For nonlinear thin materials we find E emission (r; 2ω) = −2iµω Γ(r , r; 2ω) · P NL (r ; 2ω)dV . (S8) The equivalence between Equation S2 and Equation S8 can be established by accounting for the symmetry properties of the tensorial Green function stemming from the reciprocity theorem. As a matter of fact, reciprocity imposes on the dyadic Green function the constraint Γ mn (r , r) = Γ nm (r, r ): exchanging the position of the field probe and of the emitting dipole does not induce variations in the electromagnetic response, if the respective polarization directions are simultaneously interchanged. When the j-th component of E emission is calculated, the elements of the dyadic product k Γ jk P NL k are selected. Accordingly, in Equation S2 the emitting current is taken parallel to the direction j, thus confirming the swapping of the indices in the Green function. From a physical point of view, this alternative interpretation based upon the Green function stresses out that Equation S7 accounts for the exact Green function G int of the structure, including both the reflection from the glass-air interface and the interaction with the metallic nano-structures. Direct usage of Equation S8 requires to find numerically the far field expression for Γ in each point r on the elemental unit. In the nonlinear inverse scattering approach discussed in the previous section, the numerical effort instead consists in the solution of the associated linear problem for only the direction of interest. Thus, in our case the nonlinear inverse scattering technique minimizes the demand of computational time.

Evaluation of the effective metasurface nonlinearity
The unidirectional second harmonic generation in the Type I case is macroscopically described by [52] where d eff is the effective second order nonlinearity for the given input and output polarizations, * stands for complex conjugation, and ∆k is the phase mismatch. In the case of a metasurface placed in z = 0, we can assume a delta-Dirac functional form for the nonlinearity seen on a macroscopic scale, d eff = d surf δ(z). Accounting for the metasurface thickness h = 50 nm, we can set d surf = hd meta where d meta is defined to account for the overall average longitudinal response of the metasurface. Under the no-depletion hypothesis, integration of Equation S9 along z yields for z > 0 E(2ω) = − ik 0 h n(2ω) d meta E 2 (z = 0, ω). (S11) By rewriting the electric fields in term of the respective intensity I = 0.5cn 0 |E| 2 |d meta | = κ(ω) h I(2ω) I 2 surf (ω) where we defined κ(ω) = c 0n(2ω)n 2 (ω) 2k 2 0 , I(2ω) is the second harmonic intensity emitted in the forward direction and I surf (ω) is the intensity effectively coupled to the metasurface. Given that the refractive index of a metasurface is not well defined, for simplicity we use a refractive index of n(ω) = n(2ω)=1 for both frequencies.
From an experimental point of view, the next step is to connect the intensity I surf to the impinging intensity I 0 . The connection is not straightforward, as shown in the seminal paper by Chen concerning the surface-enhancement of SHG [53]. Given we are impinging from the air side, in a first approximation the effective field on the metallic nano-structure is E inc [1 + R(θ AOI ; j)], where R(θ AOI ; j) is the Fresnel coefficient for plane waves polarized alongĵ at the interface air-glass. In this approximation we find I surf = |1 + R(θ AOI ; j)| 2 I inc . Equation S12 turns into |d meta | = κ(ω) h|1 + R(θ AOI ; j)| 2 I(2ω) I 2 inc (ω) The previous approach is valid in the limit of plane waves and monochromatic waves. For a proper comparison with our experiments, we need to account for the pulse duration, beam width and repetition rate. The instantaneous intensity can be factorized out as i inc (t) = i peak u(t)f 2 (x, y). We then set f (x, y) = e −(x 2 +y 2 )/w 2 inc , u(t) = e −2t 2 /τ 2 and i peak = 2 √ 2E/ π 3/2 τ w 2 inc , where w inc is the beam width at the metasurface, τ is the (1/e 2 ) pulse duration, and E is the pulse energy. Strictly following the concept of root mean square, the equivalent monochromatic plane wave intensity E eq (ω) is found by a spatio-temporal average across the electromagnetic pulse P eq d eff = E 2 eq (ω) ≡ where erf(x) = (2/ √ π) x 0 exp −t 2 dt is the standard error function, and erf(x) ≈ (2/ √ π) x for small enough x. For a train of pulses, the average power P is related to the pulse energy E via the repetition rate f rep , P = Ef rep . From the experimental side, the known quantities are the emitted average second harmonic power P measured (2ω) and the measured beam width w measured , the latter being measured with the camera at large enough generation. Substituting Equation S15 into Equation S13, we finally find the effective nonlinear coefficient of the metasurface |d meta | = κ(ω) h|1 + R(θ AOI ; j)| 2 π 2 3/2 w 2 inc τ f rep w measured P measured (2ω) P 2 inc (ω) Using a pump pulse with a wavelength of 1032 nm, τ = 200 fs, f rep = 200 kHz, P inc = 180 mW, ω inc = 125 µm and R(20.5°, T M ) = 0.166 leads to a SHG signal with ω measured = 40 µm and P measured = 190 pW. Therefore, the h = 50 nm thick metasurface has a second order nonlinearity d meta of 1.0 pm V −1 .