Subwavelength nonlinear phase control and anomalous phase matching in plasmonic metasurfaces

Metasurfaces, and in particular those containing plasmonic-based metallic elements, constitute an attractive set of materials with a potential for replacing standard bulky optical elements. In recent years, increasing attention has been focused on their nonlinear optical properties, particularly in the context of second and third harmonic generation and beam steering by phase gratings. Here, we harness the full phase control enabled by subwavelength plasmonic elements to demonstrate a unique metasurface phase matching that is required for efficient nonlinear processes. We discuss the difference between scattering by a grating and by subwavelength phase-gradient elements. We show that for such interfaces an anomalous phase-matching condition prevails, which is the nonlinear analogue of the generalized Snell's law. The subwavelength phase control of optical nonlinearities paves the way for the design of ultrathin, flat nonlinear optical elements. We demonstrate nonlinear metasurface lenses, which act both as generators and as manipulators of the frequency-converted signal.

M etamaterials are a class of artificial materials whose optical properties can be tailored to exhibit phenomena not commonly found in nature, such as negative 1,2 or anomalous refraction 3,4 . These unique optical properties are frequently engineered by single-or multi-layered nanometric objects, often metallic, fabricated on the surface of 'classical' standard materials. Metasurfaces constitute a particularly interesting and attractive subset of such materials leading to the possibility of designing and creating, by means of modern nanolithographic fabrication techniques, flat and ultrathin optical elements [4][5][6][7][8] . Nano-plasmonic-based metallic elements are the commonly utilized building blocks, and their (linear) optical properties are quite well understood 6,7 . The significance of phase changes across a metasurface has been recognized early on 9 and the laws of refraction across such surfaces have been recently reformulated by Yu et al. 4 in terms of a generalized Snell's law. For metasurfaces with linear phase gradient [5][6][7] , it was shown that the transverse phase gradient dF/dx must be accounted for as light crosses the metasurface (ref. 4, equation (2)). Based on these principles, the same group 10 demonstrated achromatic lenses by proper design of the phase elements in metasurfaces.
Metal-based metasurfaces can enhance optical nonlinearities at plasmonic resonances by orders of magnitude 11 . Negative refraction 12 and zero index materials 13 were demonstrated, enhancement in clusters and Fano resonances 14 were shown, and the potential for super resolution in plasmonically enhanced four-wave mixing (FWM) was also discussed 15 . As for nonlinear harmonic generation, second harmonic generation (SHG) was studied [16][17][18] and recently Segal et al. 19 and Li et al. 20 discussed the SHG signal generated on metasurfaces and demonstrated beam bending and focusing of such light. Other works, Lee et al. 21 and Wolf et al. 22 , reported the use of metasurfaces to gain control over a SHG signal generated within the nonlinear substrate on which the metasurface is located. The next, third order, nonlinearity was investigated even less, although, in principle, it exists for any structure and for any surface symmetry. As an example, FWM was demonstrated for cavities perforated on gold films 23 . However, to the best of our knowledge, the fundamental issue of phase matching across metasurfaces has not yet been thoroughly addressed. This is, in part, due to the somewhat limited phase and amplitude control over the nonlinearities of individual plasmonic element, which in many cases drives researchers to resort to periodic structures (gratings) that impose specific angles of diffraction 24 .
Based on the linear results, one may anticipate that transverse phase gradients at the interface provide an additional momentum that must be included in any nonlinear phase-matching scheme. Here, we demonstrate full phase control over nonlinear optical interactions in plasmonic metasurfaces. This control is achieved by introducing a spatially varying phase response of a metallic metasurface consisting of subwavelength nonlinear nanoantennas designed specifically for the frequency of the nonlinear signal. For such metasurfaces, we derive a new, anomalous nonlinear phasematching condition that differs from the conventional phasematching schemes in nonlinear optics. The complete phase control over the nonlinear emission enables us to design flat nonlinear optical components, such as nonlinearly blazed elements and ultrathin frequency-converting lenses with tight focusing.

Results
Finite-differences time-domain calculations. A coherent wavemixing process obeys the phase-matching condition Dk ¼ 0, where Dk is the vectorial sum of the momenta of all photons, incoming and generated, participating in the nonlinear mixing. This condition determines the direction of the coherent emission.
In FWM, for collinear input beams in a homogeneous and non-dispersive medium, the phase-matched generated beam propagates in the same direction, and this is also the case for the quasi-phase-matching scheme 25,26 , where a nonlinear material is designed with periodic reversal of the sign of the nonlinearity in the propagation direction. For non-collinear input beams, a more elaborate phase-matching scheme is required, often based on birefringence or temperature tuning. For our metasurfaces, the nonlinear phase gradient imposed by the plasmonic antennas determines the phase-matching conditions.
Optical nanoantennas, like other driven oscillators, reradiate the incoming light at the same frequency but with a shifted phase. In this work, we use rectangular nanocavities in thin gold films as optical antennas, and to convey the new physical principles more clearly, the aspect ratio (AR) of the rectangles was maintained as our single tuning parameter.
In Fig. 1, we present nonlinear finite-differences time-domain (FDTD) calculations (details given in the Methods section) performed for individual cavities of varying ARs. These separate calculations are combined to depict the transmission and nonlinear interactions in phased arrays of such cavities. In Fig. 1a,b, we show the calculated linear transmittance spectrum (intensity and phase) for light polarized along the short axis of the rectangles for a set of rectangular nanocavities of different ARs within a free-standing 250-nm thick gold film. Two distinct cavity resonances are seen 27 , and the correlation between the intensity and the phase shift acquired by the transmitted wave is clear.
Consider now a FWM configuration where two transformlimited laser pulses (with o j , k j and E j r; The third-order polarization induced at position r is given by 28 : Where w (3) is the third-order susceptibility of the metal and the fields E i (r,o i ) are the position-dependent electric fields. An antenna much smaller than the wavelength can be approximated by a point dipole (eliminating the position dependence within the antenna) and this leads to effective fields 29 is the field enhancement (a real quantity) and F(o i ) is the phase response. Therefore, we can rewrite equation (1) as Where S (3) is the effective third-order nonlinear susceptibility. The nonlinear FWM signal carries the frequency response at the fundamental frequencies through the phase factor e i 2F o 1 ð ÞÀF o 2 ð Þ ð Þ , which changes sharply for excitation close to the nanoantenna resonance. The nonlinear phase response of the nanoantennas can be directly calculated using full-wave nonlinear finitedifferences time-domain (NL-FDTD) calculations.
In Fig. 1c,d, we show the generated field at o FWM for a set of nanorectangles with varying AR. In the NL-FDTD calculation, two 60-fs long transformed-limited co-propagating pulses, with centre frequencies o 1 ¼ 800 nm and o 2 ¼ 1,088 nm, respectively, are temporally and spatially overlapped at a single nanocavity. Figure 1d depicts the linear phase accumulated by individually propagating waves at 800, 1,088 and 633 nm, and the nonlinear phase accumulated by the generated FWM beam at 633 nm. The calculated F FWM does not fit the phase of the 633-nm wave. The calculated phase difference, 2F 1 À F 2 , provides a much better estimate for the phase of the generated FWM signal, but the fit is still not perfect.
For a better description of the generated phase two additional factors need be included. The first, and less critical one, is integration over the laser bandwidth, which for these ultrashort pulses is not negligible. The more important factor is the fact that  the FWM signal is generated gradually over the length of propagation through the metasurface, so that the phase accumulated by this wave as it is building up should also be included, and can be incorporated as an effective dielectric constant e eff that appears in the nonlinear wave equation As mentioned, our direct NL-FDTD calculation provides the phase F FWM , and therefore we will use results such as of Fig. 1d for the design of our metasurface.
Experimental arrangement and observations. The experimental arrangement is the standard, forward propagating FWM configuration described in an earlier publication 27 . Two 60-fsec beams, one from a Ti:Sapphire and one from an Optical Parametric Amplifier serve as our inputs-their input angles and time delays are individually controllable. In Fig. 2a, the two ultrashort pulses, with wavevectors k 1 and k 2 , respectively, are now spatially and temporally overlapped and focused on the phase-gradient metasurface to generate a FWM signal at o FWM ¼ 633 nm and k FWM . After the sample, the fundamental beams are filtered and the FWM signal is imaged on a CCD camera which records its k-space information. We measure the FWM from two different metasurfaces-each consisting of four rectangles 450 nm apart, in the uniform case all rectangles are with AR ¼ 1.9, and in the phase gradient case the AR covers the range of AR ¼ 1.1-2.9.
In both cases we measure the FWM output angle as a function of the input angle y 800 of the o 1 ¼ 800 nm beam for normal incidence of the o 2 beam. There is an B10°difference in the output angle from the two different metasurfaces, indicating a different phase-matching condition. This new phase-matching condition stems from the additional momentum provided by the metasurface, along the gradient direction: k new FWM ¼ k FWM þ Dk x , where k FWM ¼ 2k 1 À k 2 is the conventional phase-matching condition. The net momentum provided by the metasurface to the FWM signal, Dk x ¼ (dF FWM /dx)u x , is transferred by the metasurface to the beams participating in the nonlinear conversion process.
Thus, the new anomalous phase-matching condition for FWM assumes the form: Equation (4) combined with the NL-FDTD calculation for F FWM provides a framework for the design of phase-gradient nonlinear metasurfaces. In Fig. 2e, the orange curve represents a nonlinear fit to equation (4), from which we extract a value for the phase gradient provided by the metasurface over a unit cell, in this case With the proper choice of dF FWM /dx, one can control the beam steering of the FWM emission. In Fig. 3, we show a NL-FDTD calculation, using parameters similar to the experimental, of the angle dependence of the FWM signal for different phase-gradient metasurfaces. The line fits to equation (4) using the nonlinear phase gradient shown in Fig. 1 (blue lines) are also plotted and are in good agreement to the NL-FDTD calculations.
Beyond the phase gradient unit cell, if several (many) cells are arranged in a periodic manner, they form a blazed grating. The analysis of the anomalous phase matching from such a grating should include the general theory of diffraction 24 in the linear case, or the Raman-Nath diffraction 30 in the nonlinear one. These blazed gratings, unlike gratings resulting from alternating positive-and negative-phased antennas in uniform unit cells 19,31 , are not symmetric in terms of 'positive and 'negative' transverse directions, resulting in much lower diffraction efficiency of the negatively diffracted orders (m ¼ À 1, À 2y). In Fig. 4, the scattering from blazed gratings is depicted. The spots seen on the CCD images for collinear, normal excitation are explained by the different orders of diffraction of the blazed grating. The angle of diffraction of the different diffraction orders is determined by the grating period, and agrees with the Raman-Nath diffraction formula sin y m ¼ ml FWM /L, where L is the period of the grating. High-order diffraction modes are also seen, but with weaker intensity compared with the zeroth and first order, also in accordance with the Raman-Nath diffraction theory.
FWM metasurface lenses. To further illustrate the power and flexibility of these nonlinear phase-gradient metasurfaces we designed nonlinear metalenses that focus the wavelength of choice in a specific FWM configuration. The ultrathin lenses operate by imposing a radially dependent relative phase shift at Here f is the desired focal distance of the lens and l 0 is the free-space wavelength. Metalenses with different focal lengths are shown and discussed in Fig. 5. They are made of concentric rings of phased gratings consisting of rectangular nanoantennas, and the FWM signal is a tightly focused, nearly diffraction limited Gaussian spot. Note that the experimentally observed focal spots are not as tight as the calculated ones. This discrepancy may be attributed to less than perfect fabrication accuracy, and input beams which are not fully collimated. Both factors will give rise to different focal parameters.

Discussion
In the present work, we demonstrate full control over the nonlinear phase on the subwavelength scale in phase-gradient metasurfaces. As in the linear regime, where Snell's law had to be modified, the phase control over the nonlinear nanoantennas leads to a modified manifestation of the laws governing nonlinear phenomena, such as NL scattering, NL refraction and frequency conversion. Related phase control had been previously reported for the linear case and analysed in terms of the Berry phase 9 , and more recently for nonlinear harmonic generation 19,20   subwavelength phase gradients enables treatment of any polarizations and for any input frequency. Quite a few ultrathin optical components 6,7,34-36 have been proposed and used for beam steering applications. These are based on the abrupt phase changes experienced by light on propagation through a properly designed ultrathin layer of metamaterial. The beam bending (at a specific angle) results from scattering by phase gratings generated by the proper design of the nanoantennas. These nanoantennas are generally uniform across a unit cell, and their phase changes abruptly (typically by p) at a predesigned periodicity. We have carried this concept one step further, designed and fabricated blazed metasurfaces, where the scattering is from the phase gradient across the unit cell, and carefully analysed their nonlinear optical properties. We show that FWM from such metasurfaces reveals a new feature: the scattering from a phase gradient unit cell enables anomalous phase-matching condition for FWM from such metasurfaces. This phase-matched FWM is efficient; its phase-matched direction of propagation may be controlled by proper design of the phase gradients, it enables beam bending at any angle and thus FWM focusing, and it differs from the scattering by phase gratings, as discussed in ref. 24. In all these cases, the experimental results were compared with numerical solutions (NL-FDTD, Lumerical solutions package 37 ) of the wave equations with nonlinear terms added to them. The agreement is generally very good, and whenever relevant, comparison with analytical expressions is added.
We designed and fabricated ultrathin FWM lenses, and demonstrated tight focusing with focal lengths of several microns. These FWM lenses do not have any restrictions on the symmetry of the design which is characteristic to elements based on SHG, and can be integrated in light detectors based on frequency conversion to provide more sensitive detection.
For pedagogical reasons we have limited our design to rectangles and the tuning parameter to the AR, but more generally parameters such as area of the hole, or shapes presenting multiple resonances such as V-shaped antennas can and will be used. Interestingly, even though the linear rectangular nanoatennas do not cover the full range of 2p phase shift, it is possible to obtain a 2p phase shift in the nonlinear case due to multiples or combinations of linear phases imprinted on the fundamental beams.
As is well known, metal structures are lossy, especially in the visible range, and thus the search is going on for nonmetallic replacements. This endeavour, however, is hampered by the lack of the electromagnetic enhancement readily available for metals, and presently combined structures consisting of nonlinear materials with plasmonic metallic structures seem to offer a pathway towards higher efficiencies 21 .
In conclusion, we demonstrate full control over the nonlinear phase in phase-gradient metasurfaces. We show that in such metasurfaces a new phase-matching condition applies, which differs from the phase-matching schemes known in nonlinear optics, and have designed and built ultrathin elements such as blazed elements and lenses that are based on nonlinear wave mixing. The phase control of nonlinear nanoantennas will enable the design of ultrathin nonlinear metamaterials, which can generate and control the wavefront of light, and may have implications for the next generation of efficient devices for spectroscopic (CARS or SERS) measurements.

Methods
Sample fabrication. The samples were fabricated by focused ion beam (FIB) milling on a high quality free-standing gold film. The procedure for fabrication of free-standing gold films is described in details in ref. 27 Briefly, using an e-beam evaporator, we deposited a 10-nm thick Cr adhesion layer and a 250-nm thick gold layer on a polished silicon wafer. On the opposite side of the wafer, we had previously grown a circular Si 3 N 4 mask by plasma-enhanced chemical vapour deposition. The mask was chemically etched using KOH and the remaining freestanding metallic area was etched with HCl to remove the adhesion layer.
Linear FDTD simulations. The transmittance spectrum and the relative spectral phase response of the rectangular metallic nanocavities were calculated using the commercial software Lumerical FDTD solutions. The values of the dielectric constants were taken from the data table of Gold from Palik 38 .
Nonlinear FDTD simulations. The nonlinear phase was calculated using the nonlinear material implementation of Lumerical. The base material is Palik gold, which is assumed to have instantaneous (non-dispersive) third-order nonlinearity w (3) ¼ 10 À 18 m 2 V À 2 . As input light sources we used two temporally overlapped transform-limited plane wave sources centred at o 1 ¼ 800 nm and o 2 ¼ 1,088 nm, with pulse duration 60 fs, propagating parallel to the z direction. The polarization of both pulses is perpendicular to the long axis of the rectangles. The y component of the real and imaginary parts of the electric field (that is, the propagating waves) of the FWM signal is recorded on a y normal plane spanning the whole simulation area. The dimensions of the mesh were set to dx ¼ dy ¼ dz ¼ 5 nm and perfectly matched layers were added in all dimensions.
k-Space analysis. For the k-space analysis of the FWM signal from phase-gradient antennas, the same simulation parameters used in the phase response were kept, except now the o 1 ¼ 800 nm source propagates in the x-z plane with a variable incidence angle y 800 with respect to the normal. The exit fields on the opposite side of the metasurface are recorded in a z-normal plane and projected to the far field, where the angle y FWM of the FWM signal at o FWM ¼ 633 nm is calculated as a function of y 800 .
K-space measurements. In the FWM experiments, we used the set-up described in details in ref. 27 An optical parametric amplifier, pumped by a 1-kHz amplified Ti:Sapphire laser, was used as the light source for the o 2 ¼ 1,088 nm pulses, while the pulses of the Ti:Sapphire laser that pumped the optical parametric amplifier were used as the fundamental o 1 ¼ 800 nm beam. Both o 1 and o 2 pulses have the same pulse duration of 60 fs. The beams travel two distinct optical paths, where the intensity and polarization of each individual beam could be controlled by a set of half-wave plate and polarizer to avoid optical damage to the samples. Both beams are focused and overlapped in the sample, by an objective lens of numerical aperture (N.A.) ¼ 0.42 (Mitutoyo M Plan Apo 50X Infinity-Corrected). The incident angle y 800 of the o 1 beam can be varied by controlling the lateral displacement of the beam with computer-controlled translation stage supporting a beam splitter whose primary role was to merge the optical path of two beams. The o 2 beam at 1,088 nm was always kept normal to the surface. The temporal overlap between the two beams as the input angle was varied, was monitored and controlled by properly delaying the o 1 beam. The FWM signal centred at o FWM ¼ 633 nm is collected by an objective lens with N.A. ¼ 0.42 (Mitutoyo M Plan Apo SL50X Infinity-Corrected) and focused by spherical lens of f ¼ 75 mm onto an EMCCD camera (Andor iXon DV885) The fundamental beams are filtered by a pair of shortpass filters (Thorlabs FES0700 and FESH0750). In the angular measurements in Figs 3 and 4, the focusing objective with N.A. ¼ 0.42 was replaced by a lens of focal length f ¼ 50 mm to illuminate a larger area.
Nonlinear lenses. For the design of our lenses we kept the same parameters used in the calculations of the NL phase response. However, to decrease the simulation time, we used symmetric boundary condition in the x dimension and anti-symmetric in the y dimension. The dimensions of the fine mesh around the lenses were set to dx ¼ dy ¼ 10 nm and dz ¼ 5 nm for the f ¼ 5 and 10 mm lenses and dx ¼ dy ¼ 15 nm and dz ¼ 5 nm for the f ¼ 30 and 60 mm lenses.
Nonlinear lenses measurements. The experimental set-up for measurements of the nonlinear lenses is a modification of the k-space set-up described above. Both o 1 and o 2 beams are normally incident. The focusing objective is replaced by a spherical lens with f ¼ 100 mm. The FWM signal out of the lenses is collected by the imaging (NA ¼ 0.42) objective and is imaged directly onto the EMCCD in real (physical) space. The imaging objective is supported on a computer-controlled translation stage that can vary the focal plane of the objective and record three-dimensional tomographic images of the focal region.