Generation of even and odd high harmonics in resonant metasurfaces using single and multiple ultra-intense laser pulses

High harmonic generation (HHG) opens a window on the fundamental science of strong-field light-mater interaction and serves as a key building block for attosecond optics and metrology. Resonantly enhanced HHG from hot spots in nanostructures is an attractive route to overcoming the well-known limitations of gases and bulk solids. Here, we demonstrate a nanoscale platform for highly efficient HHG driven by intense mid-infrared laser pulses: an ultra-thin resonant gallium phosphide (GaP) metasurface. The wide bandgap and the lack of inversion symmetry of the GaP crystal enable the generation of even and odd harmonics covering a wide range of photon energies between 1.3 and 3 eV with minimal reabsorption. The resonantly enhanced conversion efficiency facilitates single-shot measurements that avoid material damage and pave the way to study the controllable transition between perturbative and non-perturbative regimes of light-matter interactions at the nanoscale.

factor collective modes, which were demonstrated in Si metasurfaces 7 , or plasmonic field enhancement 5 .However, several challenges to achieving highly-efficient HHG in the strongfield regime assisted by spectrally-selective metasurfaces may be identified.First, narrow-and moderate-bandgap semiconductors, with bandgap energies Δ  that are not much larger than the laser photon energy ℏ, are damaged at moderate laser fluences due to multi-photon absorption followed by rapid free-carrier generation 15,16 .Moreover, the overabundance of free carriers can drastically reduce the quality () factor of a resonant metasurface 17 , thereby defeating its key purpose: the creation of resonantly-driven optical hot spots.Second, harmonics absorption by opaque materials reduces the HHG-emitting volume and dramatically decreases the HHG efficiency 18 .Finally, only a subset of harmonics (odd) can be produced by centrosymmetric materials.Currently, non-centrosymmetric materials enabling even-order harmonics 9,19 have not been utilized for nanostructure-based HHG: to date, high ( ≥ 4) harmonics have only been reported from nanostructures biased by an external dc field 8 or 2D semiconductors 20 .
Therefore, it is desirable to develop a photonic platform and an optical system providing both the access to non-perturbative physics (defined by a strong perturbation by a laser pulse of the electron/hole motion in their respective conduction/valence bands), as well as the ability to use HHG as a probe of the microscopic processes inside a crystal [20][21][22][23] .Such combination of a photonic platform and optical system must meet the following conditions: (a) the electronic bandgap of the constitutive material should be sufficiently large, so that multiple harmonic orders can be utilized; (b) the optical system should enable single-shot measurements that do not suffer from the inherent limitations of multi-pulse averaging, such as long-term damage [24][25][26] and measurement biases (e.g., produced by a single high-intensity outlier in a train of laser pulses); and (c) the photonic structure should enable the production of nanoscale regions of a strongly-driven material phase embedded inside a weakly perturbed phase, thus opening the possibility of studying nonlocal effects in condensed matter phase without confounding laser damage.
The transition to nonlinear carrier motion occurs when the momentum gained from the laser electric field over a single period exceeds the size of the Brillouin zone of a solid material.This condition is expressed as  ≡   /2 > 1, 14,21,27 where   = /ℏ is the Bloch oscillation frequency 28,29,30 ,  a crystalline period,  is the laser frequency, ℏ is the reduced Planck's constant, and  is the hot spot optical field.Concurrently, the injection of free carriers (FCs) into the conduction zone also takes place.The latter is governed by the dimensionless Keldysh parameter 31  =  √  * Δ g /, where  * is the effective electron mass.Approximately equal to the ratio of the carrier injection time to laser period, the Keldysh parameter characterizes electron tunneling across the bandgap.Therefore, highly-efficient non-perturbative (saturated) HHG requires that ,  −1 > 1.
The key challenge addressed by our work is finding the appropriate photonic platofrm for entering this new regime without producing large numbers of FCs that can blue-shift 32 or dampen 17 the metasurface resonance.As illustrated by Figure S1 (see Supplementary Information Section 1 for the calculation of strong-field-induced FC generation), our choices of the metasurface material and laser wavelength  = 2/ are strongly constrained if we are to access the non-perturbative regime of HHG in nanostructures.
Large electronic band gap (Δ g (dir) = 2.78 eV and Δ g (indir) = 2.24 eV ≫ ℏ) of GaP drastically reduces multi-photon absorption of MIR light (see Supplementary Information Section 1) and prevents visible HHG reabsorption up to the  = 7 harmonic frequency   ≡ .Finally, the non-centrosymmetric zincblende crystal structure of GaP enables generation of even-order harmonics from the bulk 9,23 .
This selection of the laser wavelength and the underlying metasurface material enabled us to produce record-breaking unsaturated conversion efficiencies into high harmonics even in the perturbative regime of moderate laser intensity  max (MP) ≈ 80 GW/cm 2 in the multi-pulse (MP) illumination regime.By employing single-pulse (SP) measurements, we avoid laser-induced damage and reach the non-perturbative regime of HHG for incident laser intensities as high as  max (SP) ≈ 480 GW/cm 2 .We observe a resonance-dependent saturation of the HHG at high estimated values of normalized Bloch oscillation frequencies ( ≈ 2), opening exciting new opportunities for non-perturbative light-matter interactions at the nanoscale.
The metasurfaces for enhanced HHG (Fig. 1a) were fabricated from 400 nm thick GaP films using thin film bonding, electron beam lithography and reactive ion etching (see Supplementary Information Section 2 for details).A scanning electron image of a typical metasurface sample is shown in Fig. 1b.The metasurfaces consist of densely packed domino-shaped dielectric resonant antennas (DRAs) supporting externally excited resonant electric dipole (ED) electromagnetic modes at the nominal resonant wavelength  res (0) = 3.95 µm.These modes were experimentally identified for several metasurfaces with varying dimensions (and, correspondingly, varying resonant wavelengths  res =  res (0) +  res ) using Fourier-transform infrared (FTIR) collimated beam spectroscopy 38 .At resonancesmanifested as the transmission dips in the experimental (Fig. 1d) and numerical (Fig. 1e) spectra due to the excitation of the ED modes of the DRAsmetasurfaces funnel the MIR radiation into the "hot spots" (see Fig. 1c    imaging or with a visible spectrometer; see Supplementary Information Sections 3-5 for details.
A typical HHG spectrum, with the luminescence background subtracted, is shown in Fig. 2b.
Even-and odd-order harmonics are observed in the near-infrared and visible ranges: from ℏ 4 ≈ 1.2 eV to ℏ 9 ≈ 3.0 eV (where   = 2/ is the N'th harmonic frequency).No detectible harmonic signal was observed from either unstructured GaP film of the same thickness, or the SiO2/Al2O3 substrate.The power of the 7 th harmonic (H7) emitted from the sample was calibrated using an external laser source of the known power and a similar wavelength; see Supplementary Information Section 7 for the calibration procedure details.The absolute conversion efficiency reaches a value of  7 ∼ 2 • 10 −9 for H7 at  = 80 GW/cm 2 , i.e. two orders of magnitude larger than the previous demonstration in a metasurface 7 and more than one order of magnitude larger than that in an epsilon-near-zero material 10 .Crucially, even-order (H4 and H6) harmonics were detected alongside the odd-order harmonics (H5, H7 and H9) because of the non-centrosymmetric (zinblende) crystal structure of GaP.Note that H8 was not detected in our experiment because of the combination of the indirect transitions at ℏ 8 = 2.28 eV (making GaP partially opaque at H8) and the inherently lower conversion efficiency of the even-order harmonics.The relatively low efficiency of the even harmonics can be attributed to unfavorable orientation of the GaP principal crystalline axes inside the DRAs; it can be improved by about two orders of magnitude by a judicious choice of the crystal axis orientation (see Supplementary Information Section 8 and Supplementary Fig. S6).
To validate the importance of the dipole-active metasurfaces resonances, we have investigated the dependence of the H7 conversion efficiency on the polarization of the MIR pulse.The non-resonant pump polarization along the short side of the metasurface DRAs results in the efficiency reduction by two orders of magnitude compared with the resonant one as shown in Fig. 2c.This implies that optical field enhancement inside the hot spot produced by the resonant laser polarization aligned with the dipole moment of the ED mode is essential for the high efficiency of HHG observed in our experiments.
Next, we have analyzed the polarization states of the odd-and even-order harmonics.
Specific examples for H7 and H6 harmonics are plotted in Fig. 2d for the (1, 0) diffraction order, as measured by BFP imaging.We observe that the odd harmonics (green dots) are co-polarized with the MIR pump (dashed lines).In contrast, the even harmonics (orange circles) are found to be elliptically polarized owing to the highly asymmetric structures of the even-order nonlinear susceptibility tensors  … () , 39 where the  th -order nonlinear polarization density of the medium is given by   () =  … ()   …   (see Supplementary Information Section 8 for details).For odd values of , the diagonal matrix elements of  … () dominate, and the th harmonic polarization is collinear with that of the MIR pump.In contrast, for even , the elements of the  … (𝑁) tensor are predominantly off-diagonal, thereby enabling polarization changes of the even-order harmonics.
To investigate whether the HHG in the multi-pulse (moderate peak power) regime obeys the perturbative scaling laws, we have plotted in Fig. 2e the dependences of the harmonic intensity  (𝑁) on the MIR intensity  MIR .The unsaturated dependences  () ~  MIR  are plotted as the guides for the eye.In striking difference with the previous findings of HHG in nanostructures 7,10,40 , the response of the GaP metasurface does not exhibit any appreciable saturation.We conclude that the perturbative regime of harmonics generation persists up to the maximum pump intensity ( MIR ≈  max (MP) = 80 GW/cm 2 ) used in these experiments, which is equivalent to the hot spot intensity  hs ≈ 0.7 TW/cm 2 inside the metasurface.This is in agreement with our estimates of  < 1 for this range of intensities (see Table S1 of the SOM).
Because metasurfaces subjected to multi-pulse trains were visibly damaged for incident laser intensities of order  MIR ≈ 200 GW/cm 2 , the only non-destructive pathway to accessing the nonperturbative regime of laser-matter interaction is to resort to single-pulse (SP) experiments.
Moreover, unlike MP averaging that may not provide the full picture of nonlinear processes, the SP exposures yield accurate relationships between the pulse energy, HHG signal, and the excitation site within the sample while avoiding the accumulation of MP damage [24][25][26] .In order to access the high-intensity regime (0.2 -0.6 TW/cm 2 ), we replaced the focusing optics and synchronized the elements of the setup.As schematically shown in Fig. 3a, the OPA triggers a mechanical shutter, which directs a single laser pulse to the sample and into the pick-off detection arm.The sample resides on a three-dimensional translation stage and is monitored by a visible-light imaging system (not shown).Each area of the sample is exposed to a single laser pulse by moving it out of the laser path by 50 μm after each shot.For each shot, the trigger starts the fast camera acquisition that records BFP images of the HHG pattern; a typical single-shot BFP image is shown in Fig. 3b. the farthest from (blue circles) to the closest to (purple triangles) the driver wavelength.Solid lines: best fits to the power law  (5) =   .Deviation from the expected  (5) ~  5 indicates the saturation of nonlinear response.Inset: power exponent  vs resonance wavelength λ res .The mask damage threshold and the metasurface damage threshold are shown for the most resonant metasurface  res =  res (0) = .
As an example, the zeroth diffraction order is plotted as a function of the field intensity in Fig. 3c for five different metasurfaces, from the one with the smallest detuning between the pump and the resonance (purple triangles: λ res = λ res (0) ) to the largest detuning (blue circles).The solid lines show the best fits to the power law  (5) ~   , where the exponent  is expected to be equal to 5 for the perturbative H5 process.However, in contrast with moderate-intensity data in Fig. 2e, the drastic reduction of the H5 exponent ( < 5) signifies the onset of the nonperturbative regime.The inset of Fig. 3c shows ( res ) as a function of the detuning between the incident pulse and the resonant wavelength  res .The exponent  varies monotonically between  = 2.8 for the least resonant metasurfaces to  = 0.9 for the most resonant metasurface.
The scanning electron micrographs (SEMs) of the degraded metasurfaces reveal two types of damage caused by the single pulses: the mask damage for  MIR >  max (HSQ) ≈ 280 GW/cm 2 (detachments of the HSQ cap from the GaP resonators) and the structure damage for  MIR >  max (GaP) ≈ 480 GW/cm 2 (removal of the GaP resonators from the substrate).Surprisingly, even though well-defined damage thresholds are identified by observing the metasurface degradation, no abrupt changes in HHG are experimentally observable at those threshold intensities  max (HSQ) and  max (GaP) (see Fig. 3c).The lack of any abrupt changes in the HHG dependences is attributed to the finite size of the beam: the HHG output is maintained at the beam's periphery even when the centrally positioned portion of the sample is damaged by a laser pulse.The estimated conversion efficiency of H5 in the single-pulse regime at  = 200 GW/cm 2 for sample #5 (resonant case) is  5 = 1.4 ⋅ 10 −6 , which is almost two orders of magnitude larger than that in the multi-pulse case.A comparison between various solid-state HHG sources, provided in Supplementary Table S2, shows that the GaP metasurface provides the largest specific (per unit length) conversion efficiency among all the materials provided.
One of the primary mechanisms contributing to the HHG in the non-perturbative regime is the generation of the nonlinear currents by the Bloch oscillations 29 of the FCs.The local (hot spot) field strength that does not destroy the most resonant GaP metasurface (corresponding to  MIR ≈  max (GaP) ) can be estimated to be  max (hs) ≈ 0.24 /Å (assuming a factor × 10 intensity enhancement at the hot spot), bringing the value of the Bloch frequency up to  B ≈ 2 • 10 15 s -1 .
The corresponding ratio of the Bloch frequency to the driving MIR laser frequency is  =   /2 ≈ 2.1, thus suggesting a transition to a non-perturbative response of the underlying GaP crystal (see Table S1).The anisotropic response of the electron subsystem suggests the importance of crystal lattice orientation, whereby one can tailor the contributions from different harmonics by engineering the crystal axes with respect to the nanostructure.These effects comprise an intriguing topic for future studies.
In conclusion, we have demonstrated efficient visible high harmonics generation using midinfrared resonances in ultra-thin gallium phosphide metasurfaces.Our approach provides recordhigh conversion efficiency at the nanoscale, enabled by the combination of strong hot spot enhancement of the optical field, high resilience of the underlying material to strong fields, and the low level of HHG reabsorption.Single-pulse illumination format enabled us to utilize much higher laser intensities than in the multiple-pulse format, thereby accessing the non-perturbative regime of HHG without confounding structural damage.The robustness of the metasurface to laser damage under ultra-intense illumination opens new routes to accessing strong-field regimes with tailored light fields and enables non-perturbative light-matter interactions on a chip.

Experimental parameter and ionization rates estimates
Following the estimates of the field intensity and the temporal local field enhancement factor  ≈ 3 (defined as  =  ext , where  is the peak value of the field within the GaP resonator,  ext is the external field), we can provide estimates of the peak fields within the structure, as well as the general metrics of the light-matter interactions in the metasurface.In Supplementary Table S1, three values of field intensity are given: the maximum multi-shot intensity used in Fig. 2e, minimum single-shot intensity used in Fig. 3   vac (TW/cm 2 )  (TW/cm 2 )  (V/Å)  B (10 To estimate the ionization rates in a material with the ionization potential Δ  by an ac field at frequency  and strength , we apply Perelomov-Popov-Terent'ev theory S1 which describes both the multiphoton and tunneling ionization processes.In this model, time-dependent Schrodinger equation is analytically solved in the quasi-static, single active electron approximation.We estimate the ionization rate using the following solution: where where  = 5.44 Å is the crystal lattice constant.The regime with  > 1, whereby Bloch oscillations can meaningfully contribute to the process of HHG, is denoted by the shaded pink area.We have also estimated the ionization probability that causes more than  crit = 10 19 cm −3 electron-hole pairs at  crit =  crit  3 /8 ≈ 10 −4 , where 8 stands for the number of atoms in the unit cell, as marked by the dashed line in Supplementary Fig. S1.The plasma as dense as  crit can cause significant free-carrier-induced absorption and deteriorate both the resonance of the metasurface and its HHG output.The circle indicates the estimated experimental conditions ( ≈ 2, ℏ/Δ  ≈ 0.12), which signify an important physical regime, whereby an ultra-intense laser pulse can lead to non-perturbative response without actively damaging the material or causing sub-optimal nonlinear response hindered by the generated free carriers.The family of dashed curves plotted for a narrower band gap of 1 eV show a considerably higher ionization, which may prevent efficient HHG in narrow-gap semiconductors.It is important to note that, even though the PPT model is not as accurate as the Keldysh model, S2 because it does not take into account any details of the electronic band structure, we have verified that the two model are in approximate agreement for the case of GaP, as the Keldysh model provides us with similar ionization rates for the parameters we used in our experiment.
Supplementary Fig. S1 | Single-pulse ionization probability estimates.A system with ionization potential   = 2.7  is excited by  = 100 fs pulse with a photon energy ℏ and peak electric field  for various values of  = 4, 2, 1, 0.5, 0.25 (blue, orange, red, green and purple solid curves, respectively).
The circle denotes the estimated experimental conditions, the shaded pink area denotes the region of strong-field ( > 1) excitation and the gray area denotes the dense photon-induced e-h plasma regime ( > 10 19  −3 ).The dashed lines show the same family of curves for   = 1.0 .

Sample fabrication
Crystalline GaP layer (∼400 nm) is first grown on a GaAs substrate with an AlGaInP buffer layer by metal-organic chemical vapor deposition (MOCVD).Then this structure is directly bonded to a sapphire substrate (150 μm) after depositing ∼2 μm SiO2 layers on top of both surfaces.The AlGaInP/GaAs substrate is then removed by wet etching.The fabrication of the GaP nanostructures starts with a standard wafer cleaning procedure (using acetone, isopropyl alcohol and deionized water in that sequence under sonication), followed by O2 and hexamethyl disilizane (HMDS) priming in order to increase the adhesion between GaP and subsequent spincoated electron-beam lithography (EBL) resist of hydrogen silsesquioxane (HSQ).After spincoating of HSQ layer with a thickness of ∼200 nm, EBL and development in 25% tetra-methyl ammonium hydroxide (TMAH) were carried out to define the patterns in the HSQ resist.Finally, inductively-coupled plasma reactive ion etching (ICP-RIE) with N2 and Cl2 gases was used to transfer the HSQ patterns to the GaP layer and generate the GaP nanostructures; see Fig. 1b for a scanning electron micrographs (SEM) of the sample's fragment.The orientation of the GaP crystal lattice with respect to the metasurface is visualized by orienting the [001] direction perpendicular to the plane of metasurface, and then tilting the normal to the metasurface' plane by 15° toward the [111] direction of the GaP crystal lattice.

High harmonic measurement setup
In Supplementary Fig. S2, a detailed schematic of the optical setup used for high harmonic generation is shown.The Extreme Mid-IR (EMIR) optical parametric amplifier (OPA) is a homebuilt KNbO3/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 was projected onto the sensor of a thermoelectrically cooled back-illuminated CCD camera (Princeton Instruments PIXIS 1024 BUV).To spectrally filter the individual optical harmonics for back focal plane imaging, a set of long-, short-and band-pass filters (Thorlabs FEL and FESH series, FGB37) was used.Supplementary Fig. S3 shows the back focal plane images of harmonics H4, H5, H6, H7, H9 after spectral filtration, exposing the diffraction patterns, as well as the incoherent luminescence that fills the whole aperture of the objective for some wavelength ranges.Exposure times were 3 s, 100 ms, 10 s, 3 s, 10 s, respectively.
For HHG spectra acquisition, the back focal plane of the objective lens was projected onto the entrance slit of a monochromator (Chromex 250SM scanning monochromator) coupled to the same CCD camera.In a typical raw spectroscopic image, see Supplementary Fig. S4, the left image shows the camera output in the spectral range capturing H6 and H7, where both can be discerned on top of the luminescence background.In order to subtract the incoherent background, for each wavelength, we fit the y-section of the image to a Gaussian near the zeroth order diffraction (middle 30 pixels).For Fig. 2 the camera after the monochromator (left).On the raw image, H7 and H6 are visible, both their (0; 0) and (0; ±1) orders, as well as the broadband luminescence spectrum.

Single-shot HHG and damage threshold measurements
To reach the non-perturbative intensities with our experiments, we added a functionality for the setup to be able to irradiate the sample with individual laser pulses the emitted harmonics generation at substantially higher intensities than those used in multi-pulse measurements.The focusing lens was changed to one with  = 50 mm, and the beam size was measured to be an ellipse with axes Δ FWHM = 53 μm and Δ FWHM = 43 μm by putting a 2D micro-bolometer array sensor (DataRay WinCamD-IR-BB, pixel pitch 17 μm) in the focal point of an attenuated beam.For the single-pulse acquisition, the repetition rate of the laser was lowered to 10 Hz, and the measurements were done in the back focal plane setting with triggered exposure.The software-controlled trigger from the laser was sent to the mechanical shutter (1/30 s opening time) and an oscilloscope that received the signal from an amplified PbS photodiode that detected the energy of a pick-off pulse.The diode was calibrated using a pyroelectric power meter (Gentek-EO QE-B), averaging over 5000 pulses for each power setting in the range from metasurface that has been exposed by single pulses in the energy ranges from 0.5 μJ to 3.8 μJ.
We have divided the outcomes of single-pulse irradiation into three scenarios: state '0' with no apparent damage done to the sample, state '1' with the HSQ mask getting detached (as supported by the SEM images) and state '2' with the GaP resonators partially removed from the substrate.
Supplementary Fig S5 shows the dependence of the outcome on the measured single-pulse energy.The importance of the single-pulse measurements is pinpointed by the fact that under multiple pulse irradiation, the sample gets severely damaged even at the periphery of the beam, where the intensity is very low: note the large crater at the bottom of Supplementary Fig S5, where the shutter was accidentally opened for several seconds, allowing about 30-50 pulses through at a moderate average pulse energy of about 2 μJ.

Conversion efficiency estimates
A source of cw radiation at   (1) = 532 nm, with a measured power of  cw (1) = 80 W right after the focusing lens, was attenuated by a stack of neutral density filters with the measured transmittance of 0.11 (OD1 filter), 0.01 (OD2 filter) and 7 ⋅ 10 −5 (OD4 filter), with the combined attenuation of  (1) = 7.7 ⋅ 10 −8 , yielding the overall power at the sample site of  cw (1)  (1) = 6.2 pW.This beam then passed through the rest of the setup and was detected by the camera sensor.Under the exposure time of  cw (1) = 1 s, the camera yielded  cw (1) = 1.4 ⋅ 10 6 counts at the zeroth diffraction order.Therefore, assuming the linearity of the signal over the exposure time, the sensitivity of the detection system can be estimated at  (1) =  cw (1)  cw (1)  (1)  cw

Tensor analysis of harmonics' properties
Owing to the stark difference in the nonlinear susceptibility tensor structures for the even-and odd-order harmonics, their emission efficiencies and polarization states can be strongly dependent on driver polarization and crystal orientation.Figure 2d of the main text shows the intensity of H7 and H6 as a function of analyzer angle for one of the diffracted orders, with the analyzer placed after the collection objective.While H7 is polarized along the pump radiation (90°/270° direction), H6 is elliptical, with the main semi-axis directed at around 45° with respect to the pump polarization.Also, Fig. 2b shows disproportionality between the even and odd harmonics' intensities.The difference in the polarization response and conversion efficiencies between the even-and odd-order harmonics can be qualitatively explained in terms of the nonlinear tensor symmetries.For the odd-order processes: .
However, due to the zincblende crystal structure of GaP, in the bulk for the even-order processes, = 0;  1 =  2 = ⋯ =  2 , meaning the main contributions to even-order harmonics will come from multiple off-diagonal components, opening potential to various polarization states of the output harmonics.We have analyzed the symmetry-enabled components of zincblende crystal nonlinear susceptibility tensors up to the 6 th order; see Supplementary File 1.In the configuration where the input field is in the form of  =  0

√2
(1, 1, 0) in the frame of the crystal structure, the nonlinear polarization for odd-order processes is in the form of  odd =  odd,0
In contrast, for even-order processes, the polarization is in the form of  even =  even,0 (0, 0, 1), for a numerical simulation).The metasurface was nominally designed to provide moderate | enhancement of the MIR radiation  =  res (0) .The most efficient excitation of an ED mode occurs when its spectral bandwidth matches that of the MIR pump shown in Fig. 1d in gray.

Figure 1 .
Figure 1.GaP metasurfaces for strong-field light-matter interactions in the mid-infrared | a. Illustration of the HHG process: resonant GaP metasurfaces show efficient even and odd high harmonic generation (up to order H9) due to the wide direct electronic bandgap, high refractive index and non-centrosymmetric lattice.The indirect band gap is not shown.b.Fabricated GaP metasurfaces: SEM images.c.Calculated local field map of the metasurface mode excited by a MIR pulse with  =  res (0) ; peak local field enhancement: | loc / ext | 2 ≈ 10 at resonant wavelengths.d.Collimated (normal incidence) FTIR transmission spectra of three samples with varying DRA sizes: largest (upper curve) to the smallest (lower curve) size.The second and third data sets are offset for clarity by +0.4 and +0.8, respectively.e. COMSOL simulations of d.The second and third data sets are offset for clarity by +0.5 and +1.0, respectively.Red stars indicate the estimated wavelengths of the maximum local field enhancement.

Figure 2a shows a
Figure2ashows a simplified sketch of the experimental setup for the detection and

Figure 2 .
Figure 2. High harmonic generation in the perturbative multi-pulse (MP) regime | a. Simplified schematic of the HHG detection setup, with the detection arm represented by either a spectrometer or a back focal plane (BFP) imager.b.MP-HHG spectra of the resonant sample at  MIR = 80 GW/cm 2 .The N = 8 harmonic is not observed due to the onset of indirect interband transitions in GaP.The arrows indicate the predicted HHG wavelengths.c.Polarization dependence of H5 shows two orders of magnitude contrast between the resonant (horizontal) and nonresonant (vertical) MIR polarizations with  MIR = 100 GW/cm 2 .d. Linear polarization of the odd-order (H7: green dots) and elliptic polarization of the even-order (H6: orange circles) harmonics.Dashed lines: MIR laser pulse polarization ( MIR = 80 GW/cm 2 ).e. Solid lines: HHG intensity as a function of the pump intensity for the N = 4 (dark-red), N = 5 (red), N = 6 (orange), N = 7 (green), and N = 9 (blue) orders.Dashed lines: corresponding guideto-the-eye power laws,  () ~  MIR  .

Figure 3 .
Figure 3. Single-pulse (SP) fifth harmonic generation reveals the non-perturbative regime and high damage thresholds of resonant metasurfaces | a.A setup for SP-HHG back focal plane (BFP) imaging.Single pulses pass through a mechanical shutter, split into the main beam (sample irradiation) and the pick-off beam (individual pulse power calibration).The diffracted harmonics are detected in the BFP configuration by triggered camera exposure.b.A typical BFP image of the H5 from the resonant sample at non-destructive intensities.c.Zeroth diffraction order intensity of the H5 as a function of MIR pump intensity for five different metasurfaces with resonances at  res , from of the main text, and the intensity corresponding to the middle point between the damage thresholds of the mask and the resonators.The columns show the following calculated quantities: vacuum intensity  vac , local intensity  =  2  vac , electric field strength  = √[TW/cm 2 ] 0.137 V/Å; the Bloch oscillation frequency  B = ||/ℏ, where  = 1.6 × 10 −19 C,  = 5.44 Å is the lattice constant of GaP, and ℏ = 1.05 × 10 −34 J ⋅ s is the reduced Planck's constant;  parameter,  =  B /2, where  = 4.77 ⋅ 10 14 s −1 is the pump frequency;  parameter,  = / crit , where the critical field  crit = Δ  / and Δ  = 2.78 eV is the Γ-point gap; the Keldysh parameter  =  √  * Δ g /, where  * = 0.09 is the Γ-valley effective electron mass,  = 9.1 ⋅ 10 −31 kg.
Fig. S1 as a function of the excitation field frequency for various values of  = /2ℏ, (b) of the main text, the amplitude of the Gaussian is plotted to separate the HHG signal from the luminescence background as a function of the wavelength.Supplementary Fig. S4 | HHG spectroscopy schematic.The back focal plane image of the sample was projected onto the entrance slit of the monochromator (right) resulting in the two-dimensional image at