Introduction

Nonlinear plasmonics of metallic nanoantennas is highly promising for creating subwavelength coherent light1 and electron2 sources, and for ultrafast all-optical control3. Importantly, the geometry of such nanoantennas can be flexibly designed, allowing to tune their plasmonic modes to be resonant with the fundamental or higher harmonic frequencies of the driving laser4 or even to be doubly resonant5,6. The nonlinear efficiency of pure plasmonic nanoantennas, however, is limited by their intrinsically weak optical nonlinearities and the weak penetration of light into the metal. To further enhance nonlinear emission, hybrid metal–semiconductor nanostructures have been fabricated by combining the strong field enhancements of plasmons and the large nonlinear susceptibility of semiconductors7,8,9,10,11,12,13,14. The enhancement of nonlinear emission often results from the coupling of fundamental plasmonic fields to higher-energy, electronic7 or excitonic9,14,15 transitions of quantum emitters. Although many studies have been performed on various combinations of plasmonic nanoantennas and semiconductors in the frequency domain, considerably less is known about their ultrafast optical dynamics. Accessing such dynamics is not only important to understand the microscopic origin of the coupling and the mechanism of nonlinear signal enhancement, but also necessary for full coherent control of coupled nanosystems3. Such femtosecond dynamics have been probed on pure plasmonic nanostructures and have proven crucial in understanding plasmonic coherence properties16,17,18,19. Furthermore, most studies on the hybrid nanostructures focus on enhancing second- or third-harmonic generation that doubles or triples energetically degenerate photons, respectively. Other nonlinear wave-mixing channels, such as sum-frequency (SF) generation, that provide more possibilities and tunabilities to mix different colors of fundamental photons, have been rarely explored.

Here, by using few-cycle driving pulses, we record interferometric time-resolved electron and light emission from bare Au and hybrid Au/ZnO nanosponges. We show that the coupling between plasmonic hot spot fields and excitonic quantum emitters, via SF quantum channels, boosts nonlinear excitonic emission from ZnO. The results are valuable for improving and controlling optical nonlinearity via nonlinear plasmon–quantum emitter coupling20.

Results

Nonlinear light emission from localized plasmonic hot spots

For our studies of nonlinear plasmon–exciton coupling, we take gold nanosponges—porous half-spherical nanoparticles with diameters of several hundreds of nanometers—as a promising plasmonic nanoantenna21. Such nanosponges are perforated with 10–20 nm sized nanopores throughout the entire particle, forming a randomly disordered ligament network (see scanning electron microscopy (SEM) image in Fig. 1a, more SEM images can be found in Supplementary Fig. 1). Their unique structure allows us not only to efficiently couple far-field light to the dipolar mode of the particle, but also to funnel the electromagnetic energy into a series of localized hot spots via multiple coherent scattering of surface plasmons in the disordered structure22,23. Dark-field white light scattering spectra show a broad resonance superimposed with some random spectral modulations, reflecting the excitation of multiple hot spot modes superimposed to the dipole mode of single nanosponges (Supplementary Fig. 3). Since the scattering spectra are largely inhomogeneously broadened, they do not give direct insight into the plasmonic properties of single hot spots. Light scattering from individual hot spots has recently been studied using scanning near-field optical spectroscopy, which directly shows that these hot spots are highly localized on a 10-nm scale23. The narrow spectral line widths suggest long dephasing times T2 of ~20 fs of hot spots23. Their high quality factors (up to 40) and large field enhancements make them appealing for plasmon–emitter coupling. Particularly, different types of quantum emitters can be readily infiltrated into the nanopores and thus be placed in the hot spots, without the need for lithography methods. Furthermore, the broad spectral variation of the hot spot resonances offers an opportunity to probe the effect of detuning on the coupling between plasmons and emitters.

Fig. 1: Time-resolved nonlinear light emission from long-lived localized plasmonic hot spots in a gold nanosponge.
figure 1

a Scheme of the interferometric frequency-resolved autocorrelation (IFRAC) experiment. A broadband, phase-locked 8-fs pulse pair is used to excite and time resolve the plasmonic fields of a single gold nanosponge—a porous half-spherical nanoparticle with a randomly disordered ligament network (see scanning electron microscopy image, scale bar: 100 nm). The back-scattered nonlinear signal is detected in a spectrometer as a function of inter-pulse delay τ. b IFRAC trace of a single gold nanosponge. c Zeroth-order (DC) band of the Fourier transform of the IFRAC trace along τ. The clear stretching along the detection axis ωd indicates an inhomogeneous broadening of the plasmonic hot spots of the nanosponge. Right inset: fitted homogeneous dephasing time of hot spots of 16 fs as a function of detection frequency. d Simulated DC band by taking the incoherent sum of the second-harmonic (SH) emission from five hot spots, see details in Supplementary Fig. 13. e Top: excitation-polarization-resolved SH spectra at a fixed τ 75 fs; bottom: representative SH spectra taken at θ1 = 175° and θ2 = 320°, showing broadband, multi-peaked spectral structures. f Polar plot of the integrated SH intensity. A multipole feature is characteristic for SH emission from multiple randomly orientated hot spots.

To time resolve the surface plasmon fields of the nanosponge, we performed interferometric frequency-resolved autocorrelation (IFRAC) measurements (Fig. 1a)24,25,26. The sample is excited by a pair of broadband, phase-locked 8-fs pulses, time-delayed by τ. The pulses are centered at ~880 nm and have a nearly transform-limited time structure EL(t) (Supplementary Fig. 4). The resulting local plasmon field EP(t) generates a second-order nonlinear field ENL(t, τ)  (EP(t) + EP(t − τ))2. The signal at detection angular frequency ωd is resolved in a spectrometer and recorded as a function of τ, giving an IFRAC trace, \(I_{{\mathrm{IF}}}\left( {\omega _{\mathrm{d}},\tau } \right) = \left| {\smallint E_{{\mathrm{NL}}}\left( {t,\tau } \right){\mathrm{e}}^{{\mathrm{i}}\omega _{{\mathrm{d}}t}}\mathrm{dt}} \right|^2\). Figure 1b depicts such a trace from a single gold nanosponge, showing broadband nonlinear emission centered at ωd ~ 4.4 fs−1 (~430 nm) that is modulated periodically along the delay axis τ. We find a strong and long-lived (>25 fs) coherent second-harmonic (SH) signal, superimposed on a weak incoherent two-photon photoluminescence background (Supplementary Figs. 5 and 7)25.

It has previously been shown24 that the time structure of the field that is driving the nonlinear emission can be retrieved from a Fourier transform of the IFRAC signal along τ, giving a new delay frequency ωτ. This Fourier transform shows five distinct bands, a zeroth-order (DC) band around ωτ = 0, two fundamental (FM) bands around ωτ = ±ωd/2 and two second-order bands around ωτ = ±ωd (Supplementary Fig. 8)24,25. Both DC and FM bands contain the desired information about the time profile of EP(t) (ref. 24). Figure 1c shows the DC band of the IFRAC trace in Fig. 1b. It appears to be clearly asymmetric and stretched along the ωd axis. This is surprising since the SH emission from a single mode would result in a symmetric, Lorentzian line shape of the DC band with the same, homogeneously broadened linewidth of 4/T2 along both ωd and ωτ axes (Supplementary Figs. 8 and 9). Even in the case of an inhomogeneous ensemble of coherently interfering nonlinear emitters, the DC band should also be symmetric, with a width given by the inhomogeneous linewidth (Supplementary Fig. 10)27. This is in contrast to the distinct asymmetry of the DC band seen in Fig. 1c. To explain this observation, we consider that the SH polarization from plasmonic hot spots in disordered metallic nanostructures is usually generated from an Angstrom thin sheet of nonlinear dipole emitters, orientated perpendicularly to the complex-shaped and highly curved metal surface28,29. In such a case, one expects that SH emissions from neighboring hot spots do not interfere and that the total SH intensity is given as the incoherent sum of the SH intensities of each hot spot, \(I_{{\mathrm{IF}}}^\prime \left( {\omega _{\mathrm{d}},\tau } \right) = \mathop {\sum }\nolimits_j I_{{\mathrm{IF}},j}\left( {\omega _{\mathrm{d}},\tau } \right)\), with j being the hot spot index29,30.

When considering this in the modeling of the IFRAC traces, we can reasonably well account for the asymmetric DC band (Fig. 1d), assuming an inhomogeneously broadened distribution of several hot spot modes (nonlinear emitters) with comparatively long T2 times of ~16 fs (see simulated IFRAC and more details in Supplementary Fig. 13). The T2 time is obtained by fitting the DC band signal at a given ωd to a homogeneously broadened line shape model (Fig. 1c, details in Supplementary Fig. 11). The deduced T2 times are close to those deduced using near-field linear scattering spectroscopy of individual hot spots23.

To confirm this observation, we recorded SH spectra at a fixed τ ~ 75 fs for different polarization angles of the linearly polarized incident light (Fig. 1e, completed data in Supplementary Fig. 5). The SH emission shows a rich dependence on polarization angle, with very broad and multi-peaked spectra. The polar plot of the SH intensity is indicative of a multipolar emission character (Fig. 1f), as expected for the emission from multiple randomly orientated localized plasmonic hot spots at the surface of a single nanosponge. The dipolar mode of the nanosponges can couple efficiently to the far-field light because of their submicron size that matches to the incident wavelength. The dipolar mode then couples to and excites localized hot spot modes. The linear and nonlinear fields from the hot spots can then couple back to the dipolar mode and radiate into the far field. Since the hot spot fields are much stronger than that of the dipole mode, the nonlinear emission is mainly governed by the hot spots.

Photoelectron emission from individual plasmonic hot spots

The IFRAC measurements suggest that the nonlinear emission stems from multiple localized hot spots on the surface of the nanosponges, but are unable to resolve emission from a single hot spot. For this, we employed interferometric time-resolved photoemission electron microscopy (tr-PEEM)31,32,33,34,35,36. Similar to IFRAC, the sample is excited by a pair of few-cycle pulses, centered at ~800 nm. Photoemitted electrons are detected, providing sub-100 nm spatial resolution to isolate individual hot spots. Figure 2a shows a SEM image of the gold nanosponge used in this experiment. The tr-PEEM images measured at different τ are depicted in Fig. 2b, from which three individual hot spots are well resolved. Supplementary Videos 1 and 2 show τ -dependent photoelectron emission from the nanosponge in Fig. 2a and other nanosponges. We note that the PEEM image at τ = 0 does not allow us to clearly isolate the individual hot spots, because of the overlapping emission from the delocalized dipole mode. Electron emission from this short-lived dipole mode can only be observed for short time delays, whereas for long time delays only the localized hot spot modes survive. Interferometric autocorrelation (IAC) traces of the number of detected electrons recorded at each of those hot spots are shown in Fig. 2c, displaying long-lived coherent oscillations, persisting well beyond the time resolution of the experiment of 5–6 fs (ref. 32). The pronounced fringe contrast is consistent with a high order (n = 2.5−3) nonlinear photoemission process. The distinctly different oscillation periods are a clear signature of photoemission from single hot spots with different resonance frequencies33. The IAC traces are reasonably well modeled by assuming nonlinear photoemission from a homogenously broadened hot spot emitter with a T2 time of ~13 fs, superimposed by the weak and short-lived field associated with the excitation of the delocalized dipole mode of the entire nanosponge22. The fitted T2 times of hot spots are close to those deduced from IFRAC. The tr-PEEM results thus directly confirm that each hot spot acts as an independent emitter that can be selectively excited by tuning the inter-pulse delay.

Fig. 2: Interferometric time-resolved photoelectron emission from individual hot spots of a single gold nanosponge.
figure 2

a Scanning electron microscopy image of a single gold nanosponge, scale bar: 300 nm. b Photoemission electron microscopy images of the nanosponge with its contour indicated by dashed ellipses, recorded at different time delay τ of the few-cycle pulse pair. Electron emission from three well-defined localized plasmonic hot spots is observed, which can be selectively excited by tuning τ. c Experimental (blue circles) and simulated (pink lines) interferometric autocorrelation traces of the electron emission intensity at the three hot spots. Each localized hot spot is modeled as a Lorentzian oscillator that is excited by the laser pulses and a nth-order (n = 2.5 − 3) nonlinear photoemission process is assumed. For each hot spot, the local field is given as the superposition between the hot spot field and a short-lived (dephasing time T2,dp = 3 fs) dipole mode field, delocalized over the nanosponge22. The fitting gives the following parameters: resonance frequencies \(\omega _{{\mathrm{P}},j}\) = 2.02, 2.48, and 2.37 fs−1, dephasing times T2,j = 17, 10, and 13 fs, order of nonlinearity nj = 3, 2.5, and 2.7, for hot spots 1–3, respectively.

Enhanced excitonic emission from Au/ZnO hybrid nanosponges

We now make use of this selectivity to study the coupling of individual hot spots to nonlinear quantum emitters, by depositing a thin layer of ZnO onto the nanosponge surface37. The material is chosen because of its large second-order susceptibility and sharp exciton absorption resonance that is spectrally overlapping with the broad SH emission from the ensemble of plasmonic hot spots. SEM and transmission electron microscopy (TEM) images of such a hybrid Au/ZnO nanosponge (Fig. 3a, more in Supplementary Fig. 2) show that the ZnO material nicely fills the interior of the nanopores and covers the entire gold ligament surface with a 10-nm thick layer.

Fig. 3: Enhanced coherent excitonic emission from a single hybrid Au/ZnO nanosponge.
figure 3

a Scanning electron microscopy (left, scale bar: 100 nm) and cross-sectional transmission electron microscopy (right, scale bar: 20 nm) images of a hybrid Au/ZnO nanosponge. A 10-nm thick layer of ZnO has been deposited onto their outer surface and infiltrated into the nanopores. b IFRAC trace of a single hybrid Au/ZnO nanosponge. Two distinct emission bands can be observed: a broad nonlinear plasmonic emission at 4.3 fs−1 and a distinctly enhanced excitonic emission at ZnO exciton frequency of 4.8 fs−1. c Experimental interferometric autocorrelation (IAC) traces integrated over the range of the plasmonic (bottom) and excitonic (top) emission from the IFRAC trace in b. The beating pattern showing destructive (purple region) and constructive (yellow region) interference at the excitonic emission band results from SF generation from a single hot spot. d IAC traces simulated by a nonlinear plasmon–exciton coupling model, agreeing well with the experiment.

An IFRAC trace of a single hybrid nanosponge is shown in Fig. 3b. A broad emission is observed around ωd,P ~ 4.3 fs−1 (~440 nm), very similar to that seen for the bare gold nanosponge in Fig. 1b. We thus assign this band mostly to plasmonic SH emission. Additionally, a distinctly enhanced emission at the exciton frequency of ZnO, ωX ~ 4.8 fs−1 (~390 nm)25,37,38, is observed. This band therefore reflects excitonic emission from ZnO and is almost one order of magnitude stronger than the emission from the bare gold nanosponge at the same frequency (Supplementary Figs. 6 and 16). Nonlinear optical spectra taken from many individual ZnO-infiltrated nanosponges confirm this enhancement effect37. Spectrally integrated IAC traces \(\bar I_{{\mathrm{IF}}}\left( \tau \right) = \smallint I_{{\mathrm{IF}}}\left( {\omega _{\mathrm{d}},\tau } \right){\mathrm{d}}\omega _{\mathrm{d}}\) from the plasmonic band between 4.1 and 4.65 fs−1 show coherent oscillations persisting for >25 fs (Fig. 3c, bottom), indicating SH emission from long-lived plasmonic hot spots. Interestingly, the IAC trace from the excitonic peak (4.65–4.95 fs−1) displays a characteristic beating pattern, reflecting a coherent superposition of at least two modes (Fig. 3c, top). Notably, the IAC trace of the excitonic emission is slightly longer-lived (>30 fs) than the nonlinear plasmonic emission, reflecting the coupling to the sharp ZnO exciton resonance, as will be further discussed later. The DC band of the IFRAC trace has a similarly stretched shape as that of the bare gold nanosponge, again a feature of nonlinear emission from multiple hot spots (Supplementary Fig. 12).

Quantum pathways of nonlinear plasmon–exciton coupling

The FM band of the IFRAC trace in Fig. 3b is shown in Fig. 4a, plotting the nonlinear emission intensity as a function of detection (ωd) and excitation frequency (ωτ), which is obtained by Fourier transform along τ. This two-dimensional (2D) FM spectrum allows us to visualize the effect of a selective excitation of different hot spots (vertical arrows) on the nonlinear emission (horizontal arrows). In the plasmonic emission range, the signal is distributed along the diagonal line, \(\omega _{\mathrm{d}} = 2\omega _{\uptau}\), with a slope of 2, indicating a SH process. Here, the selective excitation of hot spots with resonant frequencies ωP,j (j is the hot spot index) leads to SH emission at 2ωP,j.

Fig. 4: Two-dimensional spectra reveal quantum pathways of nonlinear plasmon-exciton coupling.
figure 4

a Fourier transformed fundamental band of the IFRAC trace in Fig. 3b, correlating fundamental excitation (vertical arrows) to nonlinear emission signals (horizontal arrows). The signal is symmetric with respect to the diagonal line with a slope of 2, indicating a second-order nonlinear process. The split signal at the exciton frequency shows that the excitonic emission is SF generated (\(\omega _{\mathrm{X}} = \omega _1 + \omega _2\)) with two fundamental modes at \(\omega _1\) 2.3 fs−1 and \(\omega _2\) 2.5 fs−1. The nonlinear plasmonic emission signal at the diagonal line corresponds to second-harmonic generation (\(\omega _{\mathrm{d}} = 2\omega _{{P},{j}}\)). b Simulated fundamental band using a nonlinear plasmon–exciton coupling model, reproducing well the features of the experiment. c Schematic illustration of nonlinear plasmon–exciton coupling. Few-cycle pulses excite several plasmonic hot spot modes |1P> with distinct resonance frequencies at the surface of the hybrid nanosponge (bottom inset: transmission electron microscopy image, scale bar: 20 nm), generating a nonlinear plasmonic polarization \(P_{{\mathrm{Au}}}^{(2)}\). In second-order coupling (Ω2), \(P_{{\mathrm{Au}}}^{(2)}\) are coherently Rayleigh scattered by the ZnO inclusion, giving a linear polarization \(P_{{\mathrm{ZnO}}}^{(1)}\). In first-order coupling (Ω1), the locally enhanced plasmonic fields drive the second-order nonlinear polarization \(P_{{\mathrm{ZnO}}}^{(2)}\) of the ZnO exciton |1X>. In all quantum pathways, the nonlinear fields can be generated via either second-harmonic or SF channels (marked by gray arrows with their lengths representing different energies). d Simulated electric fields of the three emission pathways for a coupled hot spot–exciton system from a single plasmonic hot spot at \(\omega _{\mathrm{P}}\) = 2.5 fs−1. After scaling, \(P_{{\mathrm{ZnO}}}^{(1)}\) is almost identical and in-phase with \(P_{{\mathrm{Au}}}^{(2)}\). In contrast, \(P_{{\mathrm{ZnO}}}^{(2)}\) is phase-shifted by about −π/2 and prolonged in time with respect to \(P_{{\mathrm{Au}}}^{(2)}\).

In the excitonic emission range, we observe a comparatively large intensity, about four times stronger than that in the 2D-FM spectrum of the bare gold nanosponge (Supplementary Figs. 7 and 16). Surprisingly, the signal splits into two side peaks, symmetrically positioned around the diagonal. Analytical analysis suggests that this split is the characteristic of SF generation, ωd = ω1 + ω2, i.e., the nonlinear mixing of different linear fields at ω1 and ω2 (see a completed analytical formalization in Supplementary Note 7). The horizontal cross-sectional spectrum at ωX ~ 4.8 fs−1 shows that the excitonic emission is mainly SF generated from fundamental modes at ω1 ~ 2.3 fs−1 and ω2 ~ 2.5 fs−1, whereas SH generation from 2.4 fs−1 mode is weaker. This splitting also explains the beating of the IAC trace at ωX in Fig. 3c, reflecting the coherent superposition of those two distinct fundamental modes. The 2D-FM spectrum thus nicely correlates fundamental excitation to second-order emission frequencies (see more detailed discussion and prove in Supplementary Note 7). The resonances and their line shapes in this map provide an opportunity to reveal the quantum pathways of the enhanced excitonic emission, which can either be excited through direct light–exciton coupling or driven by plasmonic near fields via linear or nonlinear plasmon–exciton coupling. We shall examine each of the three quantum pathways.

First, we assume that ZnO excitons are excited only by far-field light. In case of a narrowband laser excitation, efficient excitonic emission can only be achieved via SH generation when tuning the laser to half of the exciton frequency. In the 2D-FM spectrum, this pathway will show up as a peak at the diagonal, centered at (ωX/2, ωX), in stark contrast to the experiment. For a broadband, few-cycle laser excitation, excitonic emission can also be triggered by many SF channels, mixing a pair of colors ω1 and ω2 from the laser spectrum such that ω1 + ω2 = ωX. Hence, for each pair, we will see two split peaks at (ω1, ωX) and (ω2, ωX) in the 2D-FM spectrum. A simulation of such a 2D spectrum using a classical nonlinear oscillator model and the time profile of our laser is shown in Supplementary Fig. 17. We find that the cross-sectional spectrum I(ωτ) at ωd = ωX extends over the entire bandwidth of our laser and has a spectral width that is much wider than the experimental data. Far-field laser excitation therefore cannot account for the experiments, in agreement with observing no detectable signal from a pure, 10-nm thick ZnO layer under the same experimental conditions.

Secondly, we consider plasmon–exciton coupling. Here, the far-field laser first couples to the delocalized, dipolar plasmon mode of the nanosponge, which resonantly excites different plasmonic hot spots at the nanosponge surface22,23. This strong hot spot field, \(E_{\mathrm{P}}^{(1)}\) generates a second-order nonlinear polarization at the gold surface, \(P_{{\mathrm{Au}}}^{(2)}\), accounting for the nonlinear emission in the bare gold nanosponge (Fig. 4c, left). In the hybrid nanosponge, this nonlinear plasmonic field can be coherently scattered into the far field by the ZnO inclusion, inducing Rayleigh scattering, i.e., a linear dipole polarization \(P_{{\mathrm{ZnO}}}^{(1)}\) (via second-order coupling, Ω2, magenta arrows in Fig. 4c). This enhances the amount of light that is scattered into the far field. An estimate based on a simplified coupled dipole model (see details in Methods section and Supplementary Note 8) suggests a maximum intensity enhancement factor of \(\left| {\frac{{3\varepsilon (\omega )}}{{\varepsilon \left( \omega \right) + 2}}} \right|^2\). The dielectric function of ZnO, ε(ω) = εB + εX(ω), is dominated by a strong, almost frequency-independent background εB ≈ 4 − 5 and a much weaker resonant exciton contribution with a magnitude \(\left| {\varepsilon _{\mathrm{X}}} \right| \approx 1\) (ref. 39). Due to this strong background, the on-resonant enhancement factor of the scattering intensity at ωX should be only ~10% stronger than the off-resonant enhancement factor at ωd,P (Supplementary Note 8). That is, the enhancement is almost frequency independent. This contrasts with the experimentally observed ratio of on- and off-resonant enhancement factor of >350% (Supplementary Fig. 16). Therefore, the scattering of the nonlinear plasmonic field by a ZnO inclusion alone cannot explain the large resonant enhancement at ωX. Also, a change of the linear optical properties of the hot spot modes, specifically their linewidth and frequency11,13, due to the ZnO inclusion, does not account for this resonant enhancement since it would affect the on- and off-resonant nonlinear signal in exactly the same way13.

We are led to the third pathway: the linear plasmonic field instead couples to ZnO excitons via first-order coupling (Ω1, green arrows in Fig. 4c) by driving the exciton nonlinearly. Here, the part of \(E_{\mathrm{P}}^{(1)}\) that is spatially overlapping with the ZnO inclusion can off-resonantly drive a second-order nonlinear polarization from ZnO, \(P_{{\mathrm{ZnO}}}^{(2)}\). Earlier measurements have shown that the SH emission from ZnO thin films is strongly enhanced near the ZnO bandgap with a second-order susceptibility that is increased by more than a factor of five when approaching the exciton resonance40. Therefore, the nonlinear signal \(P_{{\mathrm{ZnO}}}^{(2)}\) will mostly be centered around ωX (refs. 25,38), showing excitonic enhancement that is consistent with the experiment. We therefore consider this first-order coupling as the main mechanism for the enhanced excitonic emission. The time structure of \(P_{{\mathrm{ZnO}}}^{(2)}\) will be given by a convolution of \(E_{\mathrm{P}}^{(1)}\) with the resonant excitonic response of ZnO. In comparison to \(P_{{\mathrm{Au}}}^{(2)}\), \(P_{{\mathrm{ZnO}}}^{(2)}\) thus displays a characteristic phase shift by ~π/2 arising from the resonant response of the exciton to the plasmonic driving field (see detailed phase analysis in Supplementary Note 9), and is considerably prolonged in time due to the persistent coherence of the excitonic transition. This is seen in the simulated electric fields for the excitation of a single hot spot at ωP = 2.5 fs−1 in Fig. 4d (blue and green curves). In contrast, the linearly scattered field \(P_{{\mathrm{ZnO}}}^{(1)}\) (pink curve in Fig. 4d) is in-phase with and shows almost the same dynamics as \(P_{{\mathrm{Au}}}^{(2)}\), resulting from the dominant dielectric background contribution as discussed above. The optical dynamics of the plasmon-driven nonlinear excitonic emission therefore depends on both the optical dynamics of the fundamental driving hot spots and the finite dephasing time of the exciton transition, as reflected by the elongated IAC trace of the excitonic emission compared to that of the plasmonic emission in Fig. 3c.

With this, we can now successfully understand the experimental results. To describe the nonlinear emission from the entire hybrid nanosponge, we assume several well-separated hot spots at its surface with different resonance frequencies (ωP,j = 2.0 − 2.6 fs−1) and dephasing times of ~16 fs. Direct evidence for the existence of these hot spot resonances, their long dephasing times and small spatial localization lengths of ~10 nm is given by scanning near-field optical spectroscopy23 and PEEM (Fig. 2) measurements. Each of these hot spots is driven by the laser field and coupled to neighboring ZnO excitons. We consider all three contributions to the nonlinear emission depicted in Fig. 4c (details of the model in Methods section). The simulated IFRAC trace (Supplementary Fig. 19), 2D-FM spectrum (Fig. 4b), and IAC traces (Fig. 3d) can reasonably well account for the experiments, as long as we assume that plasmon-driven nonlinear excitonic emission, \(P_{{\mathrm{ZnO}}}^{(2)}\) dominates the emission around ωX. Indeed, as in the experiment, the model predicts a splitting of the 2D-FM spectrum near ωX into two excitation peaks around ω1 = 2.3 fs−1 and ω2 = 2.5 fs−1, suggesting that SF generation dominates the nonlinear response in this frequency range.

Discussion

In principle, for a single plasmonic eigenmode, such a splitting is unlikely to occur for a smooth laser spectrum, which would favor SH generation from hot spot fields resonant at ωX/2. This SH channel, however, is not efficiently excited by our laser, which shows a slight dip in the spectrum at ωX/2 ~ 2.4 fs−1 (Supplementary Figs. 4 and 20), thus suppressing the coupling of this hot spot field to the ZnO exciton. The spectral structure of our laser therefore results in nonlinear excitonic emission that is off-resonantly driven by two different frequency components ω1 and ω2 of a hot spot field such that ω1 + ω2 = ωX. As seen in the 2D-FM spectra, this new, SF channel is efficient if the hot spot resonance is reasonably close to ωX/2. This SF channel is particularly helpful to enhance the overall yield of nonlinear emission in realistic samples with a large number of hybrid plasmon–exciton antennas, where disorder-induced fluctuations vary the resonance conditions from one antenna to another. In such systems, efficient nonlinear emission can be achieved via SF generation by excitation with ultrashort broadband laser pulses, even if the plasmon resonance is detuned from ωX/2, for which SH generation with narrowband excitation is inefficient.

We note that for more detuned hot spots the coupling to the ZnO excitons is reduced, as seen in Fig. 4a. This arises because the two fields for SF mixing are most likely created by the same hot spot eigenmode, and thus only one of the two fields, ω1 ≈ ωP, will be resonantly enhanced (Supplementary Fig. 20). To maximize the efficiency of SF generation, one could use multifrequency antennas with two or more energetically detuned linear plasmonic eigenmodes that are spatially partially overlapping with the same hot spot. Here, more frequency components for SF mixing can be resonantly excited and can coherently interact with the same nonlinear emitter placed in the hot spot. Typical sample geometries of multifrequency plasmonic resonators are two-arm antennas with different arm lengths41, cross antennas with two pairs of nanorods perpendicular to each other11,42, and other tailored sample geometries6,43. These structures support spatially overlapping and thus coherently interacting plasmonic modes that can lead to efficient nonlinear SF coupling to excitons. SF generation thus allows for highly tunable nonlinear plasmon–exciton coupling in a variety of hybrid nanostructures. Our result suggests new SF quantum channels for nonlinear plasmon–exciton coupling. Compared to commonly adopted SH quantum channels, SF generation offers much wider spectral tuneability for the coherent manipulation of the linear and nonlinear fields, opening up new routes for optical control in such a nonlinearly coupled nanosystem. The results also opens up the possibility to explore difference frequency generation \(\left| {\omega _1 - \omega _2} \right|\) of two plasmonic fields for coupling to low-frequency infrared or terahertz vibrations, for example, for plasmon-enhanced stimulated Raman scattering44.

In summary, we have spatially and temporally resolved individual plasmonic hot spot modes at the surface of a gold nanosponge, and studied their nonlinear coupling to ZnO excitons by probing the ultrafast optical dynamics in a hybrid Au/ZnO nanosponge. We demonstrate how fundamental plasmonic fields couple to higher-energy excitons via SF mixing, enhancing coherent excitonic emission. The much prolonged excitonic fields may provide a valuable tool for coherent control via nonlinear plasmon–emitter coupling20. Our results not only allow us to unravel the different microscopic contributions to the nonlinear signal enhancement, but may also prove helpful for the design of a new class of hybrid nanostructures with coupled, phased arrays of plasmonic antennas15, and quantum emitters with strong optical nonlinearities.

Methods

Preparation of nanosponges

Gold nanosponges were synthesized in two steps21,45. First, gold–silver alloy nanoparticles with (sub)micron diameters were fabricated by solid-state dewetting of a bimetallic Au/Ag (8 nm/20 nm thick) thin film on a SiO2/Si substrate. Subsequently, the nanoparticles were dealloyed in which silver was chemically removed in a solution of HNO3 (65 wt%), leading to nanoporous structure with 10-nm sized channels perforating the entire particle22. To prepare the hybrid nanosponge, a ZnO layer of 10-nm thickness was deposited into the nanopores of bare nanosponges using plasma-enhanced atomic layer deposition at 150 °C. The samples were annealed at 500 °C in Ar for 1 h to improve the crystallinity of ZnO. The nanosponges were transferred to indium tin oxide substrate for measurements.

IFRAC microscopy

We used a home-built nonlinear microscopy setup38 combined with broadband, 8-fs pulses from a 80-MHz Ti:Sapphire oscillator (Venteon: ultra) centered at ~880 nm. A phase-locked pulse pair is prepared in a dispersion-balanced Mach–Zehnder interferometer with inter-pulse time delay τ controlled by a piezo scanner (Physik Instrumente P-621.1CD PIHera). The laser is focused onto the sample to a diffraction-limited spot size of 1.0 µm through an all-reflective Cassegrain objective with a numerical aperture of 0.5, preserving the time structure of the pulses46. The laser power at the focus is ~0.45 mW at τ =0. The pulse in the focus has been characterized by IFRAC measurements on a 15-µm thick beta barium borate crystal and a retrieval algorithm was used to retrieve the electric field of the pulse (Supplementary Note 2), showing a nearly transform-limited time structure. A half-wave plate is used to tune the polarization of the incident pulse. The back-scattered light is spectrally filtered by a dichroic mirror and then short-pass filtered to remove the linear scattering signal. The nonlinear signal is spectrally dispersed in a monochromator (Princeton Instruments Acton SpectraPro-2500i) and detected with a liquid nitrogen cooled CCD camera (Princeton Instruments Spec-10).

PEEM measurements and fitting

Broadband 6-fs laser pulses centered around 800 nm are amplified in an optical parametric chirped pulse amplification laser system at a repetition rate of 200 kHz. They are split into a pulse pair in a dispersion-balanced Mach–Zehnder interferometer and focused onto the gold nanosponge sample in an ultra-high vacuum chamber, using a lens with a focal length of 20 cm under an angle of 65° to the surface normal (p-polarized)32. The laser pulse energy is attenuated such that only single electrons are emitted from the sample per laser shot. The emitted electrons are collected with the extractor lens of a PEEM (Focus GmbH) and imaged with a spatial resolution of ~50 nm. The time profile of the laser pulses EL(t) is characterized by the dispersion scan and frequency-resolved optical gating techniques close to the vacuum chamber with a replica of the chamber window and the focusing lens in the beam path. The extracted spectral phase from both techniques is very similar and confirms a nearly transform-limited pulse duration of 6.0 fs (full-width at half-maximum of the intensity profile).

To model the IAC traces of the electron emission shown in Fig. 2c, the local electric field at the position of each hot spot ES,j(t) is given as the convolution of the laser pulse and the local response function, which is taken as the sum of a long-lived hot spot mode (subscript j being the hot spot index) and a short-lived dipolar mode (subscript dp)22,

$${\it{E}}_{{\mathrm{S}},{\it{j}}}\left( {\it{t}} \right) = {\it{E}}_{\mathrm{L}}({\it{t}})\otimes \left\{\Theta {({\it{t}})({\it{a}}_{\it{j}}\, \mathrm{e}^{-{\mathrm{i}}{{\omega }}_{{\mathrm{P}},j}t} \mathrm{e}^{-\frac{{{t}}}{{{{T}}_{2,{{j}}}}}}\, +\, {{a}}_{{\mathrm{dp}}}\, \mathrm{e}^{ -{\mathrm{i}}{{\omega }}_{{\mathrm{dp}}}\, {{t}}} \mathrm{e}^{ - \frac{{{t}}}{{{{T}}_{2,{\mathrm{dp}}}}}})} \right\}$$
(1)

where ω and T2 are the central frequency and dephasing time, respectively, of the plasmonic modes (hot spots or dipole mode). Θ(t) is the Heaviside function. The amplitudes of the hot spot mode and dipole mode are represented as aj and adp, respectively. The IAC trace of each hot spot is then modeled as

$${{C}}_{{j}}\left( \tau \right) = \smallint \left| {\left( {{\it{E}}_{{\mathrm{S}},j}\left( {\it{t}} \right) + {\it{E}}_{{\mathrm{S}},j}\left( {t - \tau } \right)} \right)^{{\it{n}}_{\it{j}}}} \right|^2{\it{\mathrm{d}t,}}$$
(2)

where nj is the order of nonlinearity for photoemission.

In the experiment, we can clearly observe a change of the central wavelength when varying the time delay. At early time delays, both delocalized dipole mode and localized hot spot modes contribute, whereas at later time delays only the hot spot modes survive. Since the PEEM signal depends in a highly nonlinear manner on the local electric field amplitude (IPEEME6), both fields must be of very similar amplitude (adp ≈ aj) to make both contributions significant to the signals at early time delays. For fitting the experimental IAC trace, the amplitude of the dipole mode adp is kept the same for each simulation of ES,j of the three hot spots in Fig. 2. We find that the amplitudes aj of the hot spots 1–3 are 1.1, 0.84, and 0.83 as large as adp, respectively.

The fitting gives the following parameters: ωP,j = 2.02, 2.48, 2.37 fs−1, T2,j = 17, 10, 13 fs, and nj = 3, 2.5, 2.7 for hot spots 1–3, respectively. For the dipole mode, ωdp = 2.35 fs−1 and T2,dp = 3 fs. Supplementary Figure 21 shows the response functions of the dipole mode and the three hot spot modes, and their convoluted fields with laser in time and frequency domain. We note that the observed order of nonlinearly is slightly lower than estimated from a photon picture needed for the multiphoton-emission of electrons from gold. This may be explained by (i) the applied DC bias voltage bending the vacuum level and lowering the work function47, (ii) emission from local surface (adsorbate) states closer to vacuum level, and (iii) emission from a hot, nonequilibrium electron system with an effective electron temperature well above room temperature47.

Model of nonlinear plasmon–exciton coupling

We use a classical nonlinear oscillator model to simulate the second-order nonlinear emission from the sample by modeling plasmonic hot spots and exciton resonances as Lorentzian oscillators48. We consider three distinct quantum pathways that contribute to the nonlinear emission from hybrid Au/ZnO nanosponges.

The first pathway is the nonlinear plasmonic emission from each of the plasmonic hot spots, which can be excited by an external laser field EL(t),

$${\ddot {x}} + 2{\it{\upgamma }}_{\mathrm{P}} {{\dot{x}}} + {\it{\omega }}_{\mathrm{P}}^2{{x}} + {\it{a}}_{\mathrm{P}}{\it{x}}^2 = - eE_{\mathrm{L}}({\it{t}})/{\it{m,}}$$
(3)

where the displacement of the induced collective charge oscillations is denoted as x, and −e and m are electron charge and mass, respectively. The resonance frequency of the plasmonic hot spot is ωP and its damping rate γP = 1/T2,P, with T2,P being the dephasing time of the hot spot. The strength of the second-order nonlinearity of the hot spot is given as aP which relates to the second-order nonlinear susceptibility \(\chi_{\mathrm{P}}^{(2)} \propto a_{\mathrm{P}}\). We solve Eq. (3) in two steps using perturbation theory48,

$$\ddot x^{(1)} + 2{\it{\upgamma }}_{\mathrm{P}}{{\dot{x}}}^{(1)} + {\it{\omega }}_{\mathrm{P}}^2{\it{x}}^{(1)} = - eE_{\mathrm{L}}({\it{t}})/{\it{m,}}$$
(4)
$$\ddot x^{(2)} + 2{\it{\upgamma }}_{\mathrm{P}}{{\dot{x}}}^{(2)} + {\it{\omega }}_{\mathrm{P}}^2{\it{x}}^{(2)} + {\it{a}}_{\mathrm{P}}({\it{x}}^{(1)})^2 = 0.$$
(5)

Here, superscripts (1) and (2) represent the linear and nonlinear electronic motion, respectively. Equation (4) gives the linear polarization \(P_{{\mathrm{Au}}}^{\left( 1 \right)}(t) \propto x^{\left( 1 \right)}(t)\) and Eq. (5) generates the second-order surface polarization \(P_{{\mathrm{Au}}}^{\left( 2 \right)}(t) \propto x^{\left( 2 \right)}(t)\). For gold plasmons, the first-order perturbation is resonant at the frequency of the fundamental plasmon mode (ωP ~ 2.3 fs−1) and the second-order perturbation is off-resonant. This is different from the infiltrated ZnO, for which the polarization is resonantly enhanced at frequencies around the ZnO bandgap (~5.0 fs−1). Therefore, fields emitted from ZnO in the frequency range of the fundamental plasmon mode are comparatively weak and can safely be neglected. Also, the back-coupling of the resonant emission around the ZnO bandgap to the plasmonic hot spots is neglected since it is largely detuned from the plasmon resonance.

The second pathway is the nonlinear excitonic emission from ZnO. ZnO is a direct bandgap semiconductor with a bandgap of ~3.3 eV and a quite complex optical response that has been studied in detail during the past decades. For mid-bandgap optical excitation, as in the present experiments, the response is usually dominated by excitonically enhanced SH emission showing a pronounced resonant enhancement near the bandgap38,40. Even though more complex responses have been observed for strong optical driving49, we assume that this second-order nonlinear excitonic emission from ZnO is the dominant contribution to the detected signal in our experiment. It can be excited by (i) the incident far-field laser light EL(t) or (ii) the local linear and nonlinear optical polarization of plasmonic near fields x = x(1) + x(2) of neighboring plasmonic hot spots, with coupling strengths Ω1 and Ω2, respectively, defining the relation between plasmonic polarization and the field it creates at the position of the ZnO inclusion. The equation of motion follows

$$\ddot y^{(1)} + 2\gamma _{\mathrm{X}}{{\dot{y}}}^{(1)} + \omega _{\mathrm{X}}^2{\it{y}}^{(1)} = - e\left( {E_{\mathrm{L}}\left( {{t}} \right) + \Omega_{1}{{x}}^{\left( 1 \right)} + \Omega _2{\it{x}}^{\left( 2 \right)}} \right)/{\it{m,}}$$
(6)
$$\ddot y^{(2)} + 2\gamma _{\mathrm{X}}\dot y^{(2)} + \omega _{\it{{\mathrm{X}}}}^2y^{(2)} + a_{\mathrm{X}}\left( {y^{(1)}} \right)^2 = 0,$$
(7)

where y(1) and y(2) represent linear and second-order nonlinear electronic motion of ZnO, and ωX and γX are resonance frequency, and damping rate of the ZnO exciton, respectively. In contrast to gold plasmons, the linear excitation of ZnO (Eq. (6)) induced by EL(t) and the linear plasmonic field \(\Omega _1x^{\left( 1 \right)}\) is off-resonant, but is resonantly enhanced near the bandgap by the nonlinear plasmonic field \(\Omega _2x^{\left( 2 \right)}\) via second-order coupling, leading to the Rayleigh scattering as will be discussed in detail in the third pathway later. Importantly, however, the nonlinear excitonic oscillation y(2) is mainly induced by the locally enhanced linear plasmonic field \(\Omega_1x^{\left( 1 \right)}\) via first-order coupling. As discussed in the main text, laser excitation EL(t) cannot explain the experiment and can be excluded (see simulation result in Supplementary Fig. 17), and the effect of \(\Omega_2x^{\left( 2 \right)}\) is completely off-resonant (note the square in Eq. (7)). The oscillation y(2) induces a nonlinear excitonic polarization \(P_{{\mathrm{ZnO}}}^{\left( 2 \right)}(t)\) which emits into the far field, interfering with \(P_{{\mathrm{Au}}}^{\left( 2 \right)}(t)\) from the same hot spot. As discussed above, the back-coupling of this field to the plasmonic field can be neglected.

The third pathway is the Rayleigh scattering of nonlinear plasmonic field by ZnO. Here, the electromagnetic fields from the ZnO inclusions can also be emitted by resonantly scattering the local nonlinear plasmonic fields \(\Omega_2 x^{(2)}\) at the position of the ZnO inclusion. Knowing that the plasmonic fields are emitted from localized hot spots with diameters 2rP of ~10 nm (ref. 23), we model each hot spot as a small dipole emitter with nonlinear dipole moment \(p_{{\mathrm{Au}}}^{(2)} = P_{{\mathrm{Au}}}^{(2)}A_{{\mathrm{Au}}}\), where AAu is the surface area of the hot spot. This dipole moment generates a near field at frequency ω of \({\mathbf{E}}\left( {{\mathbf{r}},\omega } \right) = \frac{{\omega ^2}}{{\varepsilon _0{\it{c}}^3}}\mathop {{\mathbf{G}}}\limits^ \leftrightarrow \left( {{\mathbf{r}},\omega } \right){\it{p}}_{{\mathrm{Au}}}^{\left( 2 \right)}(\omega )\), with \(\mathop {{\mathbf{G}}}\limits^ \leftrightarrow\) denoting the near-field term of the Greens’ propagator. We model the ZnO inclusion as a small, spherical dielectric nanoparticle with radius rX. In quasi-static approximation, the linear polarizability of the ZnO particle is then given as

$${\it{\upalpha }}({{\omega }}) = 4\pi {\it{\varepsilon }}_0{{r}}_{\mathrm{X}}^{3}\frac{{{{\varepsilon }}\left( {{\omega }} \right) - 1}}{{{{\varepsilon }}\left( {{\omega }} \right) + 2}}$$
(8)

Here, the frequency dependent dielectric function of ZnO, ε(ω) = εB + εX(ω), is dominated by a strong, almost frequency-independent, real-valued dielectric background εB ≈ 4 − 5, and a much weaker resonant exciton contribution,

$$\varepsilon _{\mathrm{X}}\left( \omega \right) = \left| {\varepsilon _{\mathrm{X}}} \right|\frac{{\gamma _{\mathbf{X}}}}{{\omega - \omega _{\mathrm{X}} + {\mathrm{i}}\gamma _{\mathrm{X}}}},$$
(9)

with a magnitude \(\left| {\varepsilon _{\mathrm{X}}} \right| \approx 1\) (ref. 39). Taking the free-space Greens’ propagator, the maximum local electric field at the center of the ZnO nanoparticle created by the nonlinear plasmonic dipole moment is

$${{E}}\left( {{\omega }} \right) = \frac{1}{{2{\pi}\varepsilon _0}}\frac{{{{p}}_{{\mathrm{Au}}}^{\left( 2 \right)}\left( {{\omega }} \right)}}{{{{\bar{r}}}^3}} = \Omega _2{{x}}^{(2)},$$
(10)

with \(\bar r\) being the distance between the center of the nonlinear plasmonic dipole and that of the ZnO nanoparticle. This field is now driving the ZnO oscillator, creating a linear dipole moment

$${\it{p}}_{{\it{{\mathrm{ZnO}}}}}^{\left( 1 \right)}\left( {{\omega }} \right) = {\it{\upalpha }}\left( {{\omega }} \right){\it{E}}\left( {{\omega }} \right) = \frac{{2{{r}}_{\mathrm{X}}^3}}{{{{\bar{r}}}^3}}\frac{{{{\varepsilon }}\left( {{\omega }} \right) - 1}}{{{{\varepsilon }}\left( {{\omega }} \right) + 2}}{{p}}_{{\mathrm{Au}}}^{\left( 2 \right)}({{\omega }}).$$
(11)

The time profile of the scattered field \(p_{{\mathrm{ZnO}}}^{\left( 1 \right)}(t)\) is obtained by Fourier transform \(p_{{\mathrm{ZnO}}}^{\left( 1 \right)}\left( t \right) = \smallint p_{{\mathrm{ZnO}}}^{\left( 1 \right)}\left( \omega \right)e^{ - {\mathrm{i}}\omega t}\mathrm{d}\omega\). This contains two contributions, \(p_{{\mathrm{ZnO}}}^{\left( 1 \right)}\left( t \right) = p_{\mathrm{B}}^{\left( 1 \right)}\left( t \right) + p_{\mathrm{X}}^{\left( 1 \right)}\left( t \right)\), in which the scattered field \(p_{\mathrm{B}}^{\left( 1 \right)}(t)\) by \(\varepsilon _{\mathrm{B}}\) is essentially in-phase and has the same time dynamics as that of \(p_{{\mathrm{Au}}}^{\left( 2 \right)}(t)\) because εB is a real-valued constant. Instead, the field linearly scattered by exciton polarization \(p_{\mathrm{X}}^{\left( 1 \right)}\left( t \right) \propto y^{\left( 1 \right)}(t)\) is ~π/2 phase-shifted with respect to \(p_{{\mathrm{Au}}}^{\left( 2 \right)}(t)\) and with elongated time dynamics given by the linear convolution between \(p_{{\mathrm{Au}}}^{\left( 2 \right)}(t)\) and the response function of the exciton resonance. The total scattered field \(p_{{\mathrm{ZnO}}}^{\left( 1 \right)}(t)\) is, however, dominated by the large εB with much weaker contribution from εX(ω). This is seen in the simulated field shown as the pink curve in Fig. 4d, which has almost the same profile as \(p_{{\mathrm{Au}}}^{\left( 2 \right)}(t)\) (blue curve).

Equation (11) also allows us to estimate the enhancement factor of nonlinear emission near the ZnO bandgap resulting from Rayleigh scattering of \(p_{{\mathrm{Au}}}^{\left( 2 \right)}\) by the ZnO inclusion. Optimally, \(p_{{\mathrm{Au}}}^{\left( 2 \right)}\) and \(p_{{\mathrm{ZnO}}}^{\left( 1 \right)}\) align in parallel, resulting in a total dipole moment \(p_{{\mathrm{tot}}}\left( \omega \right) = p_{{\mathrm{ZnO}}}^{(1)}\left( \omega \right) + p_{{\mathrm{Au}}}^{(2)}\left( \omega \right)\) and an enhancement factor of the emission intensity \(F_{{\mathrm{enh}}}(\omega ) = \left| {p_{{\mathrm{tot}}}\left( \omega \right)} \right|^2/\left| {p_{{\mathrm{Au}}}^{(2)}\left( \omega \right)} \right|^2\). To maximize the enhancement, we choose \(\bar r = r_{\mathrm{X}}\), i.e., place \(p_{{\mathrm{Au}}}^{\left( 2 \right)}\) at the surface of the ZnO nanoparticle, and get

$$F_{{\mathrm{enh}}}({{\omega }}) \approx \left| {\frac{{3{{\varepsilon }}({{\omega }})}}{{{{\varepsilon }}\left( {{\omega }} \right) + 2}}} \right|^2.$$
(12)

The enhancement factor thus simply depends on the dielectric function of ZnO. Using the refractive index \(n(\omega )\) from ref. 39 and \(\varepsilon \left( \omega \right) = n(\omega )^2\), the on-resonant enhancement factor \(F_{{\mathrm{enh}}}^{{\mathrm{on}}}(\omega = \omega _{\mathrm{X}})\) at exciton frequency is estimated to be ~10% larger than the off-resonant enhancement factor \(F_{{\mathrm{enh}}}^{{\mathrm{off}}}(\omega = \omega _{{\mathrm{d}},{\mathrm{P}}})\) at the main plasmonic emission peak \(\omega _{{\mathrm{d}},{\mathrm{P}}}\) ~ 4.3 fs−1, see details in Supplementary Note 8.

In this phenomenological model, the total nonlinear emission from a single hot spot, \(I_{{\mathrm{NL}},j}(\omega )\), in the hybrid nanosponges thus stems from the above discussed three distinct pathways. These are added coherently to give \(E_{{\mathrm{NL}},j}\left( t \right) = a \cdot x^{\left( 2 \right)}\left( t \right) + b \cdot y^{\left( 2 \right)}\left( t \right) + c \cdot y^{\left( 1 \right)}\left( t \right)\) and \(I_{{\mathrm{NL}},j}\left( \omega \right) = \left| {\mathop {\smallint }\nolimits E_{{\mathrm{NL}},j}\left( t \right)\mathrm{e}^{{\mathrm{i}}\omega t}\mathrm{d}t} \right|^2\). The total emission from an individual hybrid nanosponge is then obtained by incoherently adding up the nonlinear emission from all hot spots, \(I_{{\mathrm{NL}}}^\prime \left( \omega \right) = \mathop {\sum }\nolimits_j I_{{\mathrm{NL}},j}\left( \omega \right)\), as discussed in the main text. The amplitude coefficients \(a\), \(b\), and \(c\) are taken as free scaling parameters to achieve acceptable agreement between simulations and experiments.