Quantum-mechanical simulation and χ
(2) calculation

Starting from the general form of the dipole matrix formalism for the second-order nonlinear susceptibility32,38,46, the expression for the electron (e) and the heavy-hole (hh) susceptibility of the resonant tensor element can be expressed as

$$\begin{array}{l}{\chi }_{xzx,{\rm{e}}}^{(2)}({\omega }_{1}+{\omega }_{2},{\omega }_{1},{\omega }_{2})=\displaystyle\frac{{N}_{z}{e}^{3}{r}_{{\rm{e}},\mathrm{hh}}^{2}}{6{\epsilon }_{0}{\hslash }^{2}}\sum _{{k}_{\parallel }}\sum _{m,n}\sum _{l}\\\qquad\qquad\qquad\qquad\qquad\left(\displaystyle\frac{\langle {\psi }_{\mathrm{hh},m}| {\psi }_{{\rm{e}},n}\rangle \langle {\psi }_{{\rm{e}},n}| z| {\psi }_{{\rm{e}},l}\rangle \langle {\psi }_{{\rm{e}},l}| {\psi }_{\mathrm{hh},m}\rangle }{({\omega }_{\mathrm{hh},m}^{{\rm{e}},n}({k}_{\parallel })-{\omega }_{1}-{\omega }_{2}+i\Gamma )({\omega }_{\mathrm{hh},m}^{{\rm{e}},l}({k}_{\parallel })-{\omega }_{1}+i\Gamma )}\right)\end{array}$$

(2)

$$\begin{array}{l}{\chi }_{xzx,\mathrm{hh}}^{(2)}({\omega }_{1}+{\omega }_{2},{\omega }_{1},{\omega }_{2})=\displaystyle\frac{{N}_{z}{e}^{3}{r}_{{\rm{e}},\mathrm{hh}}^{2}}{6{\epsilon }_{0}{\hslash }^{2}}\sum _{{k}_{\parallel }}\sum _{m,n}\sum _{l}\\\qquad\qquad\qquad\qquad\qquad\left(-\displaystyle\frac{\langle {\psi }_{{\rm{e}},n}| {\psi }_{\mathrm{hh},m}\rangle \langle {\psi }_{\mathrm{hh},m}| z| {\psi }_{\mathrm{hh},l}\rangle \langle {\psi }_{\mathrm{hh},l}| {\psi }_{{\rm{e}},n}\rangle }{({\omega }_{\mathrm{hh},m}^{{\rm{e}},n}({k}_{\parallel })-{\omega }_{1}-{\omega }_{2}+i\Gamma )({\omega }_{\mathrm{hh},l}^{{\rm{e}},n}({k}_{\parallel })-{\omega }_{1}+i\Gamma )}\right)\end{array}$$

(3)

where [m, n, l] denote bound states in the conduction and heavy-hole bands, ψe,hh are the corresponding envelope wavefunctions, ω1,2 denote the input photon frequencies, \({\omega }_{\mathrm{hh};m,l}^{{\rm{e}};n,l}\) represent interband transition energies, Γ is the phenomenological broadening parameter (here set to 5 meV), k∣∣ denotes the in-plane momentum, Nz is the spin degeneracy, e is the elementary charge, ℏ is the reduced Planck constant and i is the imaginary unit. The quantity \({r}_{{\rm{e}},\mathrm{hh}}=\langle {u}_{{\rm{e}}}| r| {u}_{\mathrm{hh}}\rangle\) is the interband dipole matrix element where r is the position operator, and ue and uhh are the cell-periodic Bloch functions of the electron and heavy hole states. The electron and heavy-hole susceptibilities contain the intersubband matrix elements of the electrons and heavy holes, respectively.

The envelope functions of the coupled quantum wells were obtained from self-consistent Schrödinger–Poisson simulations performed with Nextnano47. The interband matrix element re,hh for GaAs was calculated using density functional theory in the Vienna Ab initio Simulation Package with the HSE06 hybrid functional. The summation over in-plane k states was evaluated by converting it to a two-dimensional integral over (kx, ky), truncated at one-tenth of the Brillouin zone, beyond which the contribution to χ(2) was found to be negligible34.

Material growth and substrate transfer

The III–V heterostructures were grown on semi-insulating GaAs(100) wafers using molecular beam epitaxy in a Varian Gen II system. The system was equipped with solid-source thermal effusion cells for Al and Ga and a solid-source valved cracker for As. The growth temperature was maintained at 600 °C, monitored by band-edge thermometry. AlGaAs was grown at 1.85 Å s−1 and GaAs was grown at 0.83 Å s−1 under a ×15 As overpressure. The heterostructure contains 16 periods, each comprising two AlGaAs barrier layers (shown in bold) and two GaAs quantum wells, following the layer sequence (in nanometres) of 18.2/7.1/1.8/2.9. Etch-stop layers were grown to enable subsequent flip-chip transfer to Al2O3 substrates. The substrate transfer was done by bonding the heterostructure to Al2O3 with 353ND, EPO-TEK epoxy resin, followed by mechanical lapping for etching of the GaAs substrate and wet etch removal of the etch-stop layers (Supplementary Fig. 6).

Electron microscopy characterization of MQW film

Using a focused ion beam, lamellae for STEM were prepared and subsequently fine-polished with Ar using a nano-mill. Using a double-aberration-corrected JEOL ARM-300CF microscope operating at 300 kV, STEM images were acquired. We performed the z-contrast HAADF-STEM imaging with a probe convergence semi-angle of 25.7 mrad and an inner collection angle of 53 rad.

EDS maps were taken with dual 100 mm2 Si drift detectors. The EDS maps were constructed by summing 100 drift-corrected scans, each with a dwell time and step size of 10 ms and 52 pm, respectively. Subsequently, the composition profiles of Ga and Al were fitted as piecewise linear functions in the growth direction to be used for the quantum-mechanical simulations of the bound-state wavefunctions and energy levels of the grown material.

Metasurface simulations

The electromagnetic simulations were carried out with rigorous coupled-wave analysis using GRCWA48. The structural parameters used in the simulations were a pillar height of 390 nm, a radius of 230 nm, an x periodicity of 891 nm and a y periodicity of 650 nm, with p-polarized incident light. The refractive index of the TiO2 was measured by ellipsometry of an amorphous TiO2 thin film deposited by atomic layer deposition under the same conditions as the pillars of the metasurface. The refractive index of the MQW stack was determined by ellipsometry of the equivalent GaAs/AlGaAs stack on a GaAs substrate. Measured refractive indices are in Supplementary Figs. 7 and 8. In the simulations, the MQW stack was treated as a single material layer with the measured refractive index. Equivalent simulations were completed for the pump and second-harmonic wavelengths, and the corresponding pairs were used to compute the modal overlap.

Metasurface fabrication

Undiluted positive-tone electron-beam lithography resist, ZEP 520A, was spin-coated at 3,800 rpm for 45 s. The sample was pre-baked on a hotplate at 90 °C for 3 min, followed by a hotplate at 180 °C for 3 min. A conductive polymer (Showka Denko ESPACER 300) was spun at 1,500 rpm for 45 s to avoid charging effects. The resist was patterned using electron-beam lithography (Elionix BODEN 150) with an acceleration voltage of 150 kV and a current of 1 nA. The sample was developed with o-xylene under gentle agitation. TiO2 was deposited by atomic layer deposition (Savannah, by Cambridge NanoTech), followed by reactive-ion etching (Oxford PlasmaPro 100 Cobra 300) of overgrown TiO2. Finally, the resist was removed with Remover PG. A schematic is in Supplementary Fig. 9.

Experimental characterization

Linear optical characterization of the unpatterned GaAs/AlGaAs heterostructure was done using a Cary 7000 Universal Measurement Spectrophotometer. Linear and nonlinear characterization of the metasurface device was done using a MenloSystems ELMO femtosecond erbium laser centred at 1,560 nm with an average output power of 330 mW, pulse duration of 70 fs and repetition rate of 100 MHz. The linear signal was collected with an Anritsu MA9710B optical spectrum analyser, and the nonlinear response was measured with an Andor SR-500i-B2-R spectrometer with an Andor Newton 971 EMCCD detector. The optical characterization set-up is shown in detail in Supplementary Fig. 2.

Measurement of \({{\boldsymbol{\chi }}}_{{\boldsymbol{x}}{\boldsymbol{z}}{\boldsymbol{x}}}^{({\bf{2}})}({\boldsymbol{\omega }})+{{\boldsymbol{\chi }}}_{{\bf{x}}{\bf{x}}{\bf{z}}}^{({\bf{2}})}({\boldsymbol{\omega }})\) by comparison with LiNbO3

To determine the magnitude of \({\chi }_{xzx}^{(2)}(\omega )+{\chi }_{xxz}^{(2)}(\omega )\) in the MQW, we performed relative measurements against a reference sample. A 0.60-μm-thick x-cut congruent thin film of LiNbO3 on sapphire was measured at normal incidence, and a 0.60-μm-thick MQW film on sapphire substrate was measured at 45° incidence. The samples were excited from the substrate side under equal pumping conditions, and the transmitted pump and generated second-harmonic spectra were recorded. Second-order nonlinear susceptibility \({\chi }_{zzz}^{(2)}(1,567\,{\mathrm{nm}})\) values measured close to a 1,567 nm pump wavelength for LiNbO3 are lacking; the measured value referred to most often in the literature is \({\chi }_{zzz}^{(2)}(1,064\,\mathrm{nm})\)= 54.4 pm V−1 (refs. 49,50,51).

Applying Miller’s rule52 to this value predicts \({\chi }_{zzz}^{(2)}(1,567\,\mathrm{nm})\)= 51.9 pm V−1. Careful measurements by Shoji et al. found \({\chi }_{zzz}^{(2)}(852\,\mathrm{nm})\) = 51.4 pm V−1, \({\chi }_{zzz}^{(2)}(1,064\,\mathrm{nm})\) = 50.4 pm V−1 and \({\chi }_{{zzz}}^{(2)}(1,313\,\mathrm{nm})\) = 39.0 pm V−1 (ref. 53); however, the latter value deviates considerably from the Miller’s rule expectation. The highest value, \({\chi }_{zzz}^{(2)}(1,064\,\mathrm{nm})\) = 83.4 pm V−1, was reported for stoichiometric LiNbO3, which contains a higher lithium fraction than congruent LiNbO3 (refs. 50,54). For consistency with the current body of literature, this work uses \({\chi }_{zzz}^{(2)}(1,567\,\mathrm{nm})\)= 51.9 pm V−1.

To account for reflections at the sample interfaces and thin-film interference, we model propagation into our experimental structures, including the respective incidence angle, using finite-difference time-domain simulations (Flexcompute Tidy3D; https://github.com/marcus-o/linbo_mqw_comparison). To model the experimental pump-pulse spectrum, we Fourier transform the incident laser pulses’ spectrum, inject the resulting electric field using a time-dependent source and record the polarization-resolved time-dependent electric field Ex(t) and Ez(t) at the centre position of our thin films.

The use of a time-dependent electric field E(t) inherently accounts for the pulsed nature of the excitation and avoids the need for a continuous-wave approximation. In particular, the time-dependent second-order nonlinear polarization P(2) is computed as P(2)(t) ∝ χ(2)E2(t), such that the temporal profile of the pulse, including its peak intensity, is explicitly captured in the simulation. As the extraction of χ(2) is performed through a relative calibration against a LiNbO3 reference measured under identical excitation conditions, the absolute amplitude of the electric field cancels. Consequently, the extracted χ(2) does not depend on whether the pump is expressed in terms of peak power or average power.

Whereas the frequency dependence of the nonlinear response of LiNbO3 is small in our laser’s wavelength range, we account for the heterostructure’s frequency dependence by modelling χ(2) in the time domain as the product of two independent and exponentially decaying (lifetime or dephasing time, 30 fs) oscillators. The coherence lengths for second-harmonic generation from 1,550 nm to 775 nm wavelength are 1.8 μm in the heterostructure and 9.6 μm in LiNbO3. As the samples are considerably thinner, phase matching and pump depletion are negligible, and we assume that sum-frequency radiation builds up coherently along the sample. We correct for the sum-frequency propagation and absorption (n780nm = 3.39 + 0.11i from ellipsometry measurements) outside of the samples using a finite-difference time-domain simulation. Under the assumption that the laser pulses’ temporal/spectral phases do not possess strong second-order (chirp) or higher-order phase components in the samples, the simulated and experimentally measured sum-frequency spectra match55 (Supplementary Fig. 7).

Finally, we scale the heterostructure’s \({\chi }_{xzx}^{(2)}(\omega )+{\chi }_{xxz}^{(2)}(\omega )\) so the simulated and experimentally measured sum-frequency fluxes match for LiNbO3 and the heterostructure (Supplementary Fig. 10). This relative approach removes the need for an absolute pump intensity calibration and enables a direct comparison between the MQW and LiNbO3 under consistent experimental and numerical conditions.

Measurement of resonant second-harmonic generation enhancement

To measure the resonant enhancement by the GMR, we compare the second-harmonic flux generated by the metasurface sample at a 0.3° incidence angle with the sum-frequency generated by the bare MQW film at 45° incidence angle using otherwise identical excitation conditions.

Experimentally, for the metasurface sample at 0.3° incidence angle, we observe no sum-frequency radiation generated by mixing pump photons from different branches of the GMR (which would be observable as additional peaks in the sum-frequency spectrum; Fig. 4a) and no sum-frequency radiation generated by mixing one GMR photon with a non-resonant photon (which would be observable as a broad sum-frequency background; Fig. 4a). Therefore, we fit the generated second-harmonic radiation using the sum of two squared Lorentzian profiles (Fig. 4c) to extract the resonance frequencies and widths of the GMR branches.

We account for the effects of the resonant metasurface in the simulation procedure detailed above by applying a Lorentzian filter to the electric field in the nonlinear layer before calculating sum-frequency generation. The filter uses the experimentally determined centre frequency and width. We then scale the Ex(t)Ez(t) product in the heterostructure until the simulated and experimentally measured sum-frequency spectra match those of the metasurface-covered and bare heterostructures, removing the need for an absolute pump intensity calibration.