Photon acceleration and tunable broadband harmonics generation in nonlinear time-dependent metasurfaces

Time-dependent nonlinear media, such as rapidly generated plasmas produced via laser ionization of gases, can increase the energy of individual laser photons and generate tunable high-order harmonic pulses. This phenomenon, known as photon acceleration, has traditionally required extreme-intensity laser pulses and macroscopic propagation lengths. Here, we report on a novel nonlinear material—an ultrathin semiconductor metasurface—that exhibits efficient photon acceleration at low intensities. We observe a signature nonlinear manifestation of photon acceleration: third-harmonic generation of near-infrared photons with tunable frequencies reaching up to ≈3.1ω. A simple time-dependent coupled-mode theory, found to be in good agreement with experimental results, is utilized to predict a new path towards nonlinear radiation sources that combine resonant upconversion with broadband operation.

B efore the demonstration of the first laser by Theodore Maiman, light propagation was widely considered to be a linear process, with the photons not expected to interact with each other. This simple understanding of light−matter interactions was overturned in the early 1960s in second harmonic generation experiments by Franken et al. 1 . From this demonstration of the merger between two photons into a photon with doubled energy, the nonlinear optics was born. Subsequent realizations of the third 2 and higher-order 3 harmonics enabled efficient light sources 4 , high-resolution microscopy 5,6 , and produced some of the most sensitive optical characterization techniques [7][8][9] .
However, fundamental effects limit the efficiency and spectral range of the canonical nonlinear processes. Mainly, the very nature of the standard n-photon processes (where n ≥ 2 is an integer) dictates that a narrow-band laser pulse centered at the frequency ω L with the spectral width Δω L cannot produce upconverted photons with the frequencies Ω n outside of the jΩ n À nω L j≲ ffiffiffi n p Δω L spectral interval. While using highfinesse optical cavities or other resonant structures can enhance the efficiency of harmonics generation, they even further limit the spectral range of the nonlinearly generated photons. Thus, a fundamental challenge is to find a nonlinear optical process that enables efficient nonlinear frequency conversion without sacrificing the spectral bandwidth. In this article, we propose and experimentally realize one such process: upconversion of mid-infrared (MIR) light undergoing rapid blue-shifting-also known as photon acceleration 10 -in a resonant nonlinear metasurface.
The concept of photon acceleration (PA) was originally introduced in gaseous plasmas 11,12 as a process of frequency conversion that occurs when electromagnetic waves propagate in a medium with a time-dependent refractive index 13 . A reduction of the refractive index via free carrier (FC) generation results in a measurable blue-shifting regardless of whether the FCs were produced by the radiation itself 11,14 or by an auxiliary electromagnetic pulse 15 , as well as in the broadening of the spectrum 16 , which was demonstrated for harmonics generation as well 17 . PA in a solid (e.g., semiconductor) medium can be achieved at much lower laser intensities than in a gas because of the ease of FC generation 18,19 , and can be further enhanced in high-quality factor (high-Q) optical cavities. For example, loading photons into a ring microcavity 20 or a photonic crystal cavity 21 and subsequently generating FCs by an external pump while the photons are still in the cavity resulted in continuous near-infrared wavelength shifts of up to several nanometers 20,22 . However, no nonlinear manifestations were observed, primarily because of high sensitivity of high-Q resonators to high-intensity nearinfrared light. Therefore, new photon-accelerating platforms based on free-space light coupling are needed.
Recently, a new paradigm of regularly nanostructured surfaces-metasurfaces 23-25 -has been established for ultrathin nonlinear and active materials 26,27 . While metasurfaces share with optical cavities the attractive properties of high spectral selectivity and strong field concentration, their important feature is the strong coupling to free-space beams. A variety of metasurface designs have been implemented for applications as diverse as wave front manipulation (in both linear 28 and nonlinear 29 regimes), rapid amplitude and phase modulation [30][31][32][33] , as well as efficient harmonics generation [34][35][36] and all-optical modulation [37][38][39] . Of particular interest are semiconductor-based metasurfaces that utilize strong, geometry-dependent Mie-type localized modes 40 with high-Q-factors 41,42 . They have already shown record-breaking nonlinear optical performance on the nanoscale [43][44][45][46] , making them an attractive platform for observing PA and tunable and broadband optical harmonics.
Here, we design and experimentally realize a photon accelerating semiconductor infrared metasurface (PASIM) that undergoes rapid refractive index changes due to highly nonlinear photoinduced generation of free carriers (FCs) in silicon by a MIR pulse. The main idea of photon acceleration in a PASIM is given in Fig. 1a. Briefly, MIR photons interact with, and get trapped by, the metasurface. As FCs are generated by four-photon absorption (4PA), the resonant frequency of the metasurface blue-shifts, and the frequency of the trapped photons follows. Accelerated MIR photons then upconvert via the standardχ ð3Þ nonlinear process, resulting in the observed blue-shifting of the THG. The PASIM is designed to have a high-Q resonance at λ R ≡ 2πc/ω R = 3.62 μm that enables efficient four-photon FC generation at modest pulse intensities. This enables us to clearly observe the effect of PA on harmonics generation, with the peaks of the upconverted radiation appearing at frequencies of up to Ω n ≈ 3.1ω L , where ω L is the central frequency of the MIR pulses chosen at ω L ≈ ω R to maximize the effect. Moreover, we detect anomalous levels of nonlinearly generated signal with frequencies of up to ≈3.4ω L in the wings of the THG spectra, corresponding to the spectral density enhancement of ≈10 8 over the projected signal from an unstructured film. Such enhancement can only be explained by the frequency boost of the spectrally dense population of photons with initial frequencies around ω L , and their subsequent nonlinear conversion. An intuitive coupled-mode theory (CMT) model with time-dependent eigenfrequency ω R (t) and damping factor γ R (t) ≡ ω R (t)/2Q R (t) (where Q R (t) is a time-dependent quality factor of the metasurface) accurately captures most features of the experimental data and provides further insights into PA efficiency improvements, thus paving the way to future applications utilizing nonperturbative nonlinear nanophotonics.

Results
Metasurface design and fabrication. The metasurface was engineered to enhance the local fields, which is crucial for the efficiency of the nonlinear photon conversion and multiphoton absorption. In our design of the metasurface, we make use of high-Q collective resonances common in regular arrays of semiconductor particles 41,42,47 . The specific implementation of the PASIM comprised of nearly touching rectangular Si nanoantennas is shown in Fig. 1b. The quality factor of the resonance-and, hence, the local field enhancement-is controlled by the gap g between the rectangles (see Supplementary Note 1 for ultrahigh Qs of up to ≈10 4 ). For photon acceleration of broadband femtosecond laser pulses, the Q-factor of the PASIM was designed to be moderate: Q ≈ 72, as experimentally confirmed by Fouriertransform infrared (FTIR) spectroscopy measurements using collimated beams 41 , as shown in Fig. 1c. Using full-wave simulations, we determine the local MIR field enhancement to be up to |E loc /E 0 | 2 = 350, where E 0 is the amplitude of the incident electric field polarized along the short dimension of the rectangle. The sharpness of the resonance is attributed to the very small net polarization current integrated over one unit cell 48 as shown in the inset of Fig. 1c. The samples were fabricated according to a standard procedure described in Methods and Supplementary Note 2.
Upconversion spectroscopy. Intense ultrashort MIR laser pulses with a variable nondestructive fluence in the 1 < F < 6 mJ cm −2 range, with central wavelength λ L = 3.62 μm and time duration τ L ≈ 200 ± 30 fs (spectral FWHM Δλ L ≈ 150 ± 50 nm) were focused onto the PASIM to 5 < I < 30 GW cm −2 peak intensity (see Supplementary Note 3 for the details of the optical setup). The advantage of operating the PASIM in the MIR regime is that the refractive index change scales as Δn(t) ∝ − N FC (t)λ 2 , which will be crucial for PA; here, N FC (t) is the time-dependent FC density. Note that we neglect Kerr and thermal additions to the refractive index; corresponding discussion and estimates are provided in Supplementary Note 7. Kerr effect may play an important role in frequency conversion processes and may present an interesting avenue for further research in metasurfaces 49,50 .
The third-order nonlinear polarization, which is responsible for the third harmonic generation (THG), can be expressed in the perturbative regime as follows: whereχ ð3Þ ð3ω ¼ ω þ ω þ ωÞ is the third-order nonlinear susceptibility tensor of silicon, and E(ω) is the local electric field strength at the pump frequency ω. Although the resonant local field is confined to a very small volume of each nanoantenna, the cubic dependence of P (3) on the local field enables considerable THG boost in the metasurface as compared to the unpatterned silicon film. This is indeed experimentally observed in the transmitted THG spectrum plotted in Fig. 1d for the metasurface (blue curve) versus the unpatterned film of the same thickness (black curve) cases: an order of magnitude spectrally integrated THG enhancement is provided by the metasurface, and a maximum conversion efficiency estimated at around 10 −9 . While similar magnitudes of THG enhancement and conversion efficiencies have been observed in the past 51 , much more revealing is the spectrum of the THG that has never been collected from rapidly changing metasurfaces, and which reveals several unambiguous signatures of the PA process. We observe three features of the THG spectra that have not been previously experimentally observed or theoretically c Transmittance spectra of the metasurface measured using Fourier-transform infrared spectroscopy. Dashed line: central wavelength λ L of the MIR fs pulses. Inset: simulated electric field E (arrows) and intensity |E| 2 (color) inside a metasurface unit cell for the incident light at λ = λ L . The red region shows the Gaussian fit to the spectrum of the pump pulses. d Experimentally measured third harmonic generation (THG) in the metasurface (blue curve) and in an unstructured Si film (black curve) of the same thickness h. Both spectra are measured at the same MIR fluence of F = 2.3 mJ cm −2 (peak intensity I = 11 GW cm −2 ). A notable blue-shift of the upconversion peak from the metasurface with respect to the expected THG is a result of the plasma-induced acceleration of MIR photons. Inset: the same data plotted on the logarithmic scale. Gray region: the noise floor. The THG from the metasurface extends to λ = 1.06 μm, with intensities eight orders of magnitude higher than projected from the unstructured film (Gaussian fit given in dashed line): a clear indication of the emergence of the photons that are not present in the incident spectrum. e Spectra of MIR pulses transmitted through the metasurface for different input fluences revealing significant pulse self-modulation NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-09313-8 ARTICLE predicted for a solid-state medium. All three indicate that the THG takes place under nonperturbative conditions, when the optical properties of the metasurface are strongly modified while the laser pulse is interacting with it. First, the spectral peak is strongly blue-shifted with respect to its unperturbed 3ω L (corresponding to λ ð0Þ TH % 1:207 μm) spectral position measured with an unpatterned film, and with respect to the tripled frequency of the PASIM resonance at 3ω R . Second, because the spectral shift of the peak to ≈3.05ω L is larger than the full widths at half maxima (FWHMs) of the unperturbed THG spectrum and of the metasurface resonance, there exists a highfrequency spectral region of the THG (λ TH < 1.17 μm, on the blue side of the spectral peak) where the THG signal from the metasurface is more than two orders of magnitude higher than that from the Si film. This is a signature of the photons accelerated from the spectral peak of the incident pulse outside of its spectral width. In the inset of Fig. 1d, we show that signal with the wavelengths as short as λ TH~1 .06 μm (i.e. more than 100 nm on the blue side of the highest energy signal detected from the unstructured film) can be detected. For example, for two specific wavelengths, λ TH = 1.12 μm and λ TH = 1.06 μm, the measured THG intensities from the metasurface are, respectively, four and eight orders of magnitude stronger than the corresponding projections calculated from the THG spectrum that was collected from the unpatterned film (dashed line). Previously, FC-induced blue-shifting of SHG was observed in metal-based metasurface 52,53 , but PA was weak due to the low Q-factor of the plasmonic resonance.
Finally, the THG spectrum reveals another counter-intuitive property of a PASIM: it is possible to resonantly enhance a nonlinear process (THG) without sacrificing the spectral bandwidth. All three features are related to the emergence of the new, higher energy photons, and can be attributed to photon acceleration due to the dynamic multiphoton FC generation. The transmitted MIR spectra shown in Fig. 1e reveal strong power dependence, thus providing further confirmation of the nonperturbative nature of the PA. No such dependence on the MIR fluence F was observed in bulk silicon (Supplementary Note 4), thus validating the key role of the resonantly excited hot spots that enable FC generation through 4PA.
To quantify the combined process that manifests as blueshifted harmonics generation, we measured the NIR spectra as a function of the incident MIR fluence F, and compared them with the corresponding spectra generated in the unpatterned Si film (Fig. 2a). We observe that the spectral peak and width of the THG light generated in the PASIM can be controlled by incident fluence. The central THG wavelength can be blueshifted by more than 30 nm, as shown in Fig. 2c, enabling the THG with center frequencies of up to ≈3.1ω L . In contrast with the common belief that the resonant enhancement of the THG must be accompanied by spectral narrowing, our results plotted in Fig. 2d indicate the opposite. The upconverted signal has a spectrum that is up to 50% broader than that from the unstructured film, which is a clear fingerprint of PA. We therefore conclude that the perturbative approach expressed by Eq. (1) fails, suggesting the need for a more accurate model of harmonic upconversion.
Theoretical model. The observed blue-shifting, broadening and saturation of the THG spectra can be explained by a simple CMT involving a single cavity mode with a mode amplitude p(t), characterized by its time-dependent resonant frequency ω R ðtÞ ω 0 R þ Δω R NðtÞ=N max and damping factor γ R ðtÞ γ 0 R þ Δγ R NðtÞ=N max , and coupled to the incident optical fieldẼðtÞ according to refs. 54,55 : where κ is the coupling constant. Here, we assume a Gaussian incident laser pulse withẼðtÞ ¼ , where τ L = 105 fs and ω L ¼ ω 0 R , andĨ is the intensity of the pump. The model does not aim to reproduce the peak intensity of the resulting harmonics, which is affected by the absolute normalization of the MIR pulse intensity. The resonant frequency/   In Eq. (2), the unperturbed ω 0 R and γ 0 R are obtained by fitting the transmission spectrum obtained with FTIR (see Fig. 1c), the coupling constant κ ¼ ffiffiffiffiffiffiffi 2γ 0 R p 54 is calculated by neglecting nonradiative losses in the absence of FCs, and N max is set to be the maximum FC density achieved in the experiments. The carrier-induced shifts are assumed to be proportional to carrier density N(t) < N max produced via the 4PA: whereĨ max is the maximum intensity. We further assume that Δω R = 2×10 13 s −1 and Δω R = 2.5Δγ R ; these quantities are close to those obtained by pump-probe measurements (see Supplementary Note 5). The resulting dynamic frequency sweep profiles are shown in Fig. 3a for three characteristic intensities ofĨ ¼ 0:6Ĩ c ,Ĩ c and 1:3Ĩ c , corresponding to three different regimes of PASIM operation; the meaning of the characteristic critical intensityĨ c will be clarified below. After numerically solving Eq. (2) to obtain the enhanced nearfields inside the metasurface ∝p(t), we calculate their spectra |p(λ)| 2 (see Fig. 3b) and observe their significant blue-shifting and spectral broadening with the increasing incident intensityĨ. In the time domain, the latter feature translates into shorter and more intense bursts of the electric field in Si, thus giving rise to more intense nonlinear THG spectra I TH (λ TH ) given by where λ TH = 2πc/ω TH . The resulting upconversion spectra are displayed in Fig. 3c, d for several values of the input light intensityĨ. In Fig. 3c, we illustrate the three regimes of PASIM operation by comparing the spectra with those obtained for a set of fixed ω R . At low intensities, e.g.,Ĩ ¼ 0:6Ĩ c (blue curve), the shape of the upconverted spectrum does not differ from that produced by a time-independent metasurface with ω R ¼ ω 0 R . As the intensity increases, the upconverted light progressively blue-shifts, and for the critical intensityĨ ¼Ĩ c (green curve) we start observing PAinduced photons that cannot be produced by any static metasurface with time-independent resonances. The transition to the PA regime is best illustrated for the intensityĨ ¼ 1:3Ĩ c (red curve in Fig. 3c). It is instructive to attempt reproducing the nonlinear spectrum by using a set of stationary resonant frequencies/linewidths corresponding to different values of N j /N max , where 1 ≤ j ≤ 6. We observe that, without invoking PA, it is impossible to reproduce the red curve in Fig. 3c by any one of such static metasurfaces, or even by using a weighted average of their corresponding THG spectra. The reason for that is the   By comparing the theoretical spectra shown in Fig. 3d with their experimental counterparts shown in Fig. 2a, we find that the simple model of PA described by Eqs. (1-4) semi-quantitatively captures the key spectral features of the PASIM-based THG: simultaneous increase in intensity, blue-shift, and spectral width with increasing laser intensity. The model rather accurately captures the saturation of the THG as shown in Fig. 3e, the value of the blue-shift (~25 nm), and the spectral broadening from the initial linewidth of the resonance (~22 nm) up to~40 nm, as shown in Fig. 3f. Even though the dynamic resonance sweep is a very complex process, especially because of the highly localized nature of FC generation, we have demonstrated that our model can explain the main features of the experiment, and can be potentially used for optimizing future metasurface designs.

Discussion
We have estimated the energy portion of the upconverted radiation that is unreachable by a passive frequency sweep, η PA ≈ 22%, by calculating the shaded area in Fig. 3c relative to the area covered by the gray curves. Comparing to previous demonstrations of PA in microcavities 20,21 , which yielded up to η PA ≈ 50%, we note that the value of η PA demonstrated here is achieved at much larger relative wavelength shifts (about 2.7% of the central wavelength versus <0.5% shown before), and without any external pump. Moreover, the ultrathin nature of the PASIM potentially enables tunable high-harmonics generation despite their finite absorption in Si. While comparable PA-associated blue-shifts of the THG have been observed 17 in ionizable gases, the required laser intensities were in the PW cm −2 range, i.e., almost five orders of magnitude higher than those used in our experiments.
As an application example, we will discuss how the PA mechanism enables a new approach to improving photons' capture and nonlinear conversion: utilizing optical pulse shaping/ chirping 56,57 to engineer the right-and left-hand sides of Eq. (2). We demonstrate that significant enhancement of the bandwidth and intensity of the THG can be achieved by matching the instantaneous frequency of a chirped laser pulse to the timedependent resonant frequency of the metasurface. Figure 4 shows the calculated THG spectra produced by the interaction of a time-dependent PASIM with ω R (t) = ω R (0)(1 + αt) (where α = 0.05 ps −1 ) with an incident chirped laser pulse whose electric field is given by where the last term describes a linear frequency chirp with the normalized rate of change δ ch . By choosing the FWHM of the pulses' intensity spectrum Δω L , and the instantaneous width of the metasurfaces' resonance Δω R in such a way that Δω L ) Δω R , we demonstrate that the benefits of the high-Q metasurface and a broadband incident laser pulse can be combined to achieve broadband THG with high conversion efficiency. Four cases are considered in Fig. 4a: cases 1-2   evolving resonance frequency ω R (t) = ω R (0)(1 + αt) excited by either a chirped or a transform-limited (compressed) MIR pulse, and cases 3-4 of a static metasurface with ω R (t) = ω R (0) excited by the same pulses. The corresponding parameters of the pulses and metasurface are listed in the caption of Fig. 4. The chirped and compressed pulses are chosen to have identical spectral intensities. The predictions of the CMT calculations for these four cases are plotted in Fig. 4d. The PASIM excited by both the chirped pulse (case 1) and the compressed pulse (case 2) produces broadband THG, with former being much more efficient than the latter. The nonlinear response of the PASIM strongly depends on δ ch , including its sign. Specifically, the enhancement by the PASIM is the highest if the chirp parameter δ ch matches the evolution rate α of the PASIM's resonance. This can be clearly observed in Fig. 4e, where the THG spectra from the PASIM are plotted for different values of δ ch . To our knowledge, the PASIM concept is the first proposal for a chirp-sensitive nonlinear metasurface. In contrast, the passive metasurface (cases 3 and 4) produces narrow-band THG with low efficiency because they do not utilize the portion of the optical bandwidth that is outside of the metasurface bandwidth.
The importance of utilizing the entire bandwidth of the laser pulse can be most easily appreciated from the following observation. We predict that a compressed THG signal from a PASIM pumped by an optimally chirped laser pulse can be as short as τ TH ≈ 80 fs despite the long lifetime of the mode τ R ≡ 1/γ R ≈ 1 ps. This is accomplished by using an additional passive dispersive element (e.g., a pair of compressing gratings) that turns a long chirped THG pulse emerging from the PASIM into a transformlimited short pulse while preserving its bandwidth. The comparison between the four cases discussed above is shown in Fig. 4f. It is apparent that only in case 1 an intense short THG pulse is produced. These results establish that a high-Q photonaccelerating metasurface can exhibit simultaneously efficient and broadband response, and thus provides a path toward surpassing the time-bandwidth limit found in systems with resonances 58 .
As an outlook, here we speculate on the possibilities of scaling our approach to other wavelength ranges, as well as outline several potential applications of PASIMs. Many applications benefit from efficient frequency conversion in the near-infrared and the visible. In order to scale the design down by a factor of 6, so as to push the PASIM operation down the visible, the smallest feature (the gap size) should be on the order of about 35 nm at the same height-to-gap aspect ratio of 2:1. Current e-beam resist technology can produce sub-20-nm resist feature sizes 59 that, in appropriate etching conditions, will result in the desired pattern. As far as the specific semiconductors, GaP is arguably the best candidate for applications of PA in the visible, for it has a large refractive index (n = 3.2-4) and band gap (2.25 eV) required for resonant operation using collective modes used in this work. The main obstacle to large frequency shifts at shorter wavelengths is the λ 2 -scaling of the Drude term. However, this is not the only known term that affects the refractive index of semiconductors. In our previous work 39 , GaAs was used near the band gap, where the refractive index was affected by band filling and band shrinkage effects 60 . Therefore, by judicious choice of materials and band gap engineering in ternary semiconductors, PA may be extended to near-IR and visible frequencies.
Below we briefly discuss a potential application of PASIMs: filling the spectral gap between high optical harmonics, with the objective of generating satellite-free isolated attosecond pulsessomething that is currently accomplished using bulky optical components 61 . We propose to utilize the process of highharmonic generation (HHG) by an MIR pulse which nonlinearly interacts and gets spectrally broadened by a time-dependent PASIM. The role of the PASIM is to sufficiently broaden the spectrum of the pulse so that adjacent (N′th and N + 1′th, where N ) 1) harmonics spectrally overlap. A similar PASIM based on a high-mobility semiconductor (GaAs) with Q 10 3 would be utilized, and its resonance wavelength swept by Δλ R /λ R~0 .1, where λ L ≈ λ R = 3.6 μm and τ comp L $ 170 fs are the wavelength and FWHM duration of the transform-limited intense MIR pulse used for metasurface-enhanced HHG. Such sweep of the metasurface resonance frequency requires FC concentrations on the order of N FC~3 ×10 18 cm −3 due to the low effective mass of electrons in GaAs 60 . In Supplementary Fig. 8, we show that the spectrum of the N = 15th harmonic generated from a PASIM is sufficiently broad to fill the spectral gap between the Nth and N + 1st harmonics. Because the PASIM-induced spectral broadening of the Nth harmonic is proportional to N, and high harmonics up to N = 32 have already been produced in solids 62 , we anticipate that PASIMs could play a pivotal role in obtaining resonant solid-state HHG continua for the generation of attosecond pulses in the extreme UV.
To summarize, we have provided the first demonstration of photon acceleration in ultrathin semiconductor metasurfaces by observing a blue-shifted THG with central frequencies of up to 3.1ω L . Relative wavelength shifts as high as 2.7% have been observed under moderate laser intensities owing to excitation of collective high-Q metasurface resonances. Using a CMT with time-dependent mode parameters, we have validated our experimental findings and estimated the overall photon acceleration efficiency at around η PA ≈ 22%. In the measured spectra, new frequencies of up to ≈3.4ω L have emerged, with the spectral intensity of up to eight orders of magnitude higher than the projected intensity from an unstructured silicon film. These findings indicate that photon-accelerating nanostructures represent a novel time-dependent nonlinear photonic platform that can find various applications in novel pulsed light sources.

Methods
Sample fabrication and characterization. Samples of silicon metasurfaces were fabricated at the Cornell Nanoscale Facility (CNF) from a silicon-on-insulator wafer (600 nm undoped Si device layer on top of a 460 μm top-grade sapphire from University Wafer) using the following recipe. The substrate was cleaned with acetone, isopropanol and O 2 plasma; PMMA 495 was spun to form a 400-nm-thick layer and baked for 15 min at 170°C; PMMA 950k was spun to form a 100-nmthick layer and baked for 15 min at 170°C; E-spacer conducting layer was spun at 6000 rpm; the pattern was exposed at 1000 μC cm −2 (JEOL 9500FS) and developed in MIBK:IPA 1:3 solution; a 60-nm-thick Cr mask was electron-beam-evaporated and lifted off in sonicated acetone for 1 min; the pattern was transferred to the silicon layer through HBr reactive ion etch (Oxford Cobra). Finally, Cr was removed with the commercially available Cr wet etchant.
Infrared spectroscopy. Bruker Vertex 80 FTIR spectrometer was upgraded to an external transmittance spectroscopy setup as described elsewhere 41 that enables collimated MIR spectroscopy, with the MIR beam focused to a spot size of about 300 μm in diameter using a pinhole imaging technique. The transmitted beam was sent to the detector and Fourier-analyzed by the spectrometer. Normalization was done using the signal from the clear sapphire area.
Nonlinear optical measurements. In Supplementary Fig. 3, a schematic of the optical setup used for nonlinear measurements is shown. The Extreme Mid-IR (EMIR) optical parametric amplifier (OPA) is a homebuilt KNbO 3 /KTA 3-crystal/ 3-pass OPA. EMIR is pumped by The Ohio State University's GRAY laser, a homebuilt 80-fs Ti:Sapphire chirped pulse amplification system with a central wavelength of 780 nm and 4 mJ per pulse. The repetition rate of EMIR can be varied nearly continuously between 1 and 500 Hz using an external Pockels-cellbased pulse picker. EMIR was used to generate 200-fs mid-IR pulses with up to 40 μJ per pulse. The output wavelength of EMIR can be varied continuously from λ = 2.7 to 4.5 μm. For the experiments, the MIR (idler) beam was fixed at λ = 3.62 μm. The 780 nm NIR, MIR, and λ = 1 μm (signal) output beams are separated spatially, with the 1-μm signal being dumped and 780 nm pump being retained for use in pump-probe experiments. The residual NIR and MIR beams are roughly collimated to a size of about 2.5 mm. The NIR pulses were found to have a pulse duration of 200 ± 15 fs. Output modes were characterized for several different wavelengths using a WinCamD-FIR2-16-HR 2 to 16 μm Beam Profiler System. Residual NIR pulse length was characterized using a BBO crystal-based near-IR autocorrelator. MIR pulse duration was measured using an AGS-crystal-based MIR autocorrelator for 3 and 3.6 μm.
MIR spectra were obtained using a homebuilt spectrometer based on a ThorLabs GR1325-30035 blazed ruled diffraction grating with a blaze wavelength of 3.5 μm as the dispersion element and the beam profiler sensor as the detector array. On the setup schematic, an inset demonstrates a typical image of the diffracted MIR beam. The spectrometer was calibrated with an A.P.E. Wavescan USB MIR spectrometer, which, due to a low sensitivity and operation speed, could not be used for the routine MIR spectroscopy.
For pump-probe and upconversion spectroscopy, the horizontally polarized MIR pulses first pass through a waveplate-polarizer assembly for precise energy control. The pulses travel through a variable delay line after which they are recombined with the NIR pulses via a dichroic mirror. The NIR pulses follow a separate but similar path. The collinear beams are focused using a CaF 2 f = 100 mm plano-convex lens. In the sample plane, the spot sizes were found to be 300 μm FWHM for MIR and 400 μm FWHM for NIR. Both spots fit within the 500 × 500 μm 2 structured area of the metasurfaces. The relative delay between MIR and NIR pulses was controlled dynamically using either the manual MIR delay line or the electronically controlled NIR delay line with sub-ps resolution. For selftuning of the resonance, the NIR beam is blocked with a beam block. MIR fluences were varied from 1 to 6 mJ cm −2 (see Supplementary Table 1) and NIR fluences were varied from <1 to 4 mJ cm −2 for the experiments. As a control method, an Si wafer was pumped in place of the sample. With NIR beam blocked, contamination of scattered light from the NIR and 1 μm signal was measured at the sample location. The NIR content was found to be 0.5 pJ per pulse and 1 μm signal was estimated to be of order 1 pJ per pulse. These pulse energies were determined to be insignificant to affect the sample during the experiment.
Upon transmission through the sample, the MIR beam and any upconversion signal were collected with a CaF 2 f = 50 mm bi-convex lens. Any residual NIR was filtered using an Si window. In one configuration, a blazed grating/MIR camera combination is used as a high-resolution MIR spectrometer. In another configuration, a commercial Ocean Optics NirQuest spectrometer (900-2500 nm) is used for detection of the upconverted radiation. THG signal was powercalibrated using the signal beam from the OPA at λ = 1.2 μm that had a known power, after being attenuated by a set of neutral density filters with a known (measured) transmittance at this wavelength. By dividing the mean power of the THG beam by the mean power pump beam, an estimate maximum conversion efficiency of 10 −9 was obtained.
Finite element simulations. We used COMSOL Multiphysics to model the response of the metasurfaces (no free carriers) by defining the computational domain as a slab with the dimensions of p x × p y × 3 μm, where p x = 2.1 μm and p y = 2 μm. Periodic boundary conditions were used for the domain boundaries parallel to the x−z and y−z planes, and perfectly matched layers were used for the domain boundaries parallel to the x−y. The dimensions of the metasurface were chosen to match those obtained from the SEM image. Wavelength-independent refractive indices of n Si = 3.45 and n Al 2 O 3 ¼ 1:7 were used for Si and sapphire, respectively.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.