Giant enhancement of third-harmonic generation in graphene-metal heterostructures

Nonlinear nanophotonics leverages engineered nanostructures to funnel light into small volumes and intensify nonlinear optical processes with spectral and spatial control. Due to its intrinsically large and electrically tunable nonlinear optical response, graphene is an especially promising nanomaterial for nonlinear optoelectronic applications. Here we report on exceptionally strong optical nonlinearities in graphene-insulator-metal heterostructures, demonstrating an enhancement by three orders of magnitude in the third-harmonic signal compared to bare graphene. Furthermore, by increasing the graphene Fermi energy through an external gate voltage, we find that graphene plasmons mediate the optical nonlinearity and modify the third-harmonic signal. Our findings show that graphene-insulator-metal is a promising heterostructure for optically-controlled and electrically-tunable nano-optoelectronic components.

Nonlinear nanophotonics leverages engineered nanostructures to funnel light into small volumes and intensify nonlinear optical processes with spectral and spatial control. Due to its intrinsically large and electrically tunable nonlinear optical response, graphene is an especially promising nanomaterial for nonlinear optoelectronic applications. Here we report on exceptionally strong optical nonlinearities in graphene-insulator-metal heterostructures, demonstrating an enhancement by three orders of magnitude in the third-harmonic signal compared to bare graphene. Furthermore, by increasing the graphene Fermi energy through an external gate voltage, we find that graphene plasmons mediate the optical nonlinearity and modify the third-harmonic signal. Our findings show that graphene-insulator-metal is a promising heterostructure for optically-controlled and electricallytunable nano-optoelectronic components.
The strong light-matter coupling regime can be reached by concentrating light into nanometric volumes, opening a wide range of applications that extend from optical sensing 1 to quantum technologies. 2 To this end, metallic nanostructures that support plasmonic excitations -the coherent oscillations of conduction electronsare widely used to intensify electromagnetic near fields, and turn out to be particularly important for enhancing nonlinear optical processes on the nanoscale. 3 Although some aspects of plasmons can be modified by the geometry and optical properties of the host conductive media, 4 creating an actively tunable plasmonic platform remains an open challenge when relying on traditional plasmonic materials, such as noble metals. [5][6][7] In addition, plasmons in noble metals face intrinsic ohmic losses that limit both their lifetimes and their optical nonlinearities. 8 This has motivated the development of schemes to mitigate such losses by, for example, employing lattice resonances in engineered nonlinear metasurfaces. 9 Recently, a variety of nonlinear optical effects in graphene -including third-harmonic generation (THG), [10][11][12][13] four-wave mixing, 11,14 the optical Kerr effect, 15 and high-harmonic generation 16,17 -have been observed, demonstrating that graphene exhibits an intrinsically strong and actively tunable nonlinear optical response. Although, the optical nonlinearity of graphene is relatively efficient when normalised to the number of carbon atoms involved in the process, the atomically thin character of this material reduces the interaction volume. In this context, the use of photonic waveguides has been shown to substantially improve the nonlinear response of graphene. 18 Besides possessing an intrinsically large optical nonlinearity, graphene can sustain surface plasmon polaritons. These plasmons, which have exceptionally long lifetimes, are highly-confined and can be electrically tuned across a wide spectral range. [19][20][21][22] Graphene-insulatormetal heterostructures have been used to demonstrate strong optical field confinement, down to single-atom length scales, 23 as well as near-perfect absorption of impinging light beams. [24][25][26] It has also been argued that such systems can reach the strong-coupling quantum regime. 20,27,28 Nonetheless, in spite of the extensive theoretical work predicting coherent nonlinear plasmonic effects in graphene, 27,[29][30][31][32][33][34] very few experiments have experimentally investigated plasmonic excitations in the nonlinear response regime of this material. [35][36][37] Moreover, the existing observations have required nonlinear mixing, 35 THz radiation, 37 or directly patterned graphene 36 to excite plasmons.
Here, we demonstrate efficient THG in graphene assisted by metallic elements that enhance light coupling into plasmons of the monatomic carbon layer. We study heterostructures made up of nanometer-thick gold nanoribbons placed above graphene, from which they are separated by an insulating spacer layer. Metal nanoribbons play a double role to intensify the electric field of a far-field mid-IR incident light beam into graphene and launch graphene surface plasmons. 23 Our experiments demonstrate that gold nanoribbons serve as efficient nanoantennas, increasing the observed THG intensity by three orders of magnitude above that of bare graphene.
We confirm that THG originates in the graphene layer, which allows us to actively tune the enhanced nonlinear signal by controlling the Fermi energy of graphene using an externally applied voltage. This degree of tunability is simply not possible in conventional nobel metal plasmonics because the carrier density in those materials is too large to be substantially modified through gating. Additionally, our experiments reveal the role of graphene plasmons on the third-harmonic signal, emerging in particular witnessed when varying the carrier concentration. We further observe that when the incident photon wavelength is tuned, the observed plamonic feature appears at different carrier concentrations, as predicted by our simulations, in which acoustic graphene plasmons are excited by the nanoribbons and confined below the metal. Moreover, eliminating the plasmons from our model results in poor agreement with experimental data. These signatures of plasmon-enhanced and -suppressed thirdharmonic generation provide a new route toward the amplification and control of light at deep subwavelength lengthscales.

THG IN GRAPHENE HETEROSTRUCTURES
Our samples are van der Waals heterostructures, consisting of a graphene sheet with a metallic nanoribbon array placed a few nanometers above it and separated by either an insulating Al 2 O 3 (3 − 20 nm) or a monolayer h-BN spacer, as depicted in Fig. 1a. All of the gold nanoribbons are 8 nm thick, with an additional 2 nm titanium adhesion layer between the spacer and the gold nanoribbons. In order to isolate the nonlinear signal from the heterostructure, we use a modified z-scan setup with a tight depth of focus (Fig. 1d). In our configuration, the sample is moved through the focal point of a fspulsed mid-IR incident light beam (with a wavelength of 5.5 µm nm or 3.9 µm), and a third-harmonic signal (at 1.833 µm or 1.3 µm, respectively) is measured in transmission (as detailed in the Methods section). All the measurements in this work are performed under ambient pressure and temperature conditions. A set of representative z-scan measurements is presented in Fig. 2a, showing that we only observe signals from bare graphene, and gold+graphene heterostructures (see Figure S2 of the Appendix for more detail). The spectrum of the nonlinear signal (with 3.9 µm incident light) is presented in Fig. 2c, showing a clear peak at the third-harmonic wavelength. The wavelength of the THG signal with 5.5 µm incident light is confirmed in Fig. S1 in the Appendix.
Our data clearly show that the THG signals from the heterostructures are greatly enhanced compared to bare graphene, and additional control experiments (see Figure  S2 of the Appendix) show that the metal structures alone do not produce a measurable THG signal. Moreover, the THG signal is maximised when the polarisation is perpendicular to the direction of the nanoribbons (red squares in Fig. 2d,e). Additionally, as shown in Fig.   2d, the THG signal of bare graphene is co-polarised with the incoming light (red triangles), and the THG signal of the heterostructures is perpendicular to the nanoribbons (red squares). In both cases, the strongly polarised signal indicates a coherent nonlinear process.
We quantify the enhancement and verify the thirdorder nature of our signal by measuring the power dependence of bare graphene (red triangles) and the gold+graphene heterostructures (squares), shown in Fig.  2b. The slope of the linear fits on a log-log scale is fixed to 3 and the y-intercepts are free parameters, which we use to calculate the enhancement of the heterostructures over bare graphene. Although at higher powers a small saturation effect can be observed in the gold+graphene data, a clear third-order power scaling is supported by the data. To explain the saturation, we model the effect of the increasing incident light power increasing the electronic temperature (see Methods section Electron temperature). The result for bare graphene, plotted as the dashed curve, fits our data well. We experimentally find the maximum enhancement for a 3.9 µm incident wavelength with a monolayer h-BN spacer and a ribbon width of W = 200 nm. Under these conditions, the heterostructures produce a THG signal that is 1600±800 times larger than that of bare graphene, which corresponds to a maximum THG conversion efficiency of 2 · 10 −7 % at 1.9 mW of incident power.
To understand the enhancement mechanism and the role of the metal, we perform rigorous coupled-wave analysis (RCWA) simulations that are presented in Fig. 3a, the details are presented in the Methods section. Notice that due to the form of the nonlinear field integral (Eq. 7), the analogous expression for higher-harmonics predicts an even larger enhancement. The simulations show a strong concentration of the electric field in the gap between the nanoribbons only when the polarisation of the incident field is perpendicular to the ribbons. In contrast, the bare graphene signal is independent of the incident polarisation (red triangles Fig. 2c). Note that the slight asymmetry is caused by a polarisation-dependent detection efficiency of our superconducting detector. Therefore, we conclude that the enhancement is mediated by the gold nanoribbons, which amplify the electric nearfield in the graphene layer. The simulations in Fig. 3a also show that the field strength in the gap depends on the width of the nanoribbons. To verify this experimentally, we performed a series of THG measurements for different nanoribbon widths, with a spacer thickness of s = 5 nm and an incident light wavelength of 5.5 µm (results for other spacers are shown in Fig. S4). From the THG signals we estimate an effective third-order susceptibility χ (3) as described in the Methods section. The result is shown in Fig 3a. The experimental data agree well with our RCWA simulations, which assume an uncertainty of ±20 % on the gap size (nominally set to 50 nm) caused by experimental imperfections. Figure 1. Gate-tunable graphene heterostructures. a. Graphene is encapsulated by a few-nm-thick Al2O3 or a monolayer h-BN film, setting the space s between the graphene and the gold nanoribbons. The gold nanoribbon arrays are characterised by the ribbon width W and the inter-ribbon gap g. Normally-incident light of frequency ω undergoes third-harmonic generation, which is collected in transmission. A gate voltage (V) tuned from −150 V to +150 V sets the graphene Fermi energy EF. b. The conical electron dispersion relation of graphene can be tuned in resonance with one, two, or three incident photons. c. Scanning electron microscopy image of one of our high-quality gold nanoribbon arrays. d. Sketch of our experimental setup. Difference-frequency generation (DFG) between signal and idler beams of an optical parametric oscillator (not shown) provides mid-IR ≈ 260 fs pulses. A half-wave plate (HWP), together with a polariser (pol), selectively rotate the linear polarisation of the incident light, which is then focused onto the sample. A second lens collimates the incident and outgoing third harmonic light. A band-pass filter (BP) isolates the third-harmonic signal, which is coupled into a multimode fiber and sent to a superconducting-nanowire single-photon detector (SNSPD). The sample is moved in the z direction (z-scan), in and out of the focal point of the beam.

ELECTRICAL TUNABILITY
While using gold nanribbons of the appropriate width can greatly enhance the nonlinear response of the system, this width cannot be actively changed. In contrast, the optical nonlinearity in graphene can be electrically tuned, thus providing a practical route to active control. The optical response in graphene depends on the strength of the intraband and interband transitions, 11,13,34,38 which in turn depends on the ratio of the impinging light energy ω 0 to the graphene Fermi energy E F . The latter can be tuned in-situ by applying an external voltage to the graphene layer relative to the silicon substrate.
Conceptually, we can understand the effect of E Ftuning on THG as illustrated in Fig. 1b. By sweeping the gate voltage, one can match the Fermi energy to be resonant to an interband transition for one, two or three incident photons. Each transition results in a different resonance in the third-order nonlinear susceptibility given by, and H is the Heaviside step function. 39 However, the resonances enter the expression with different signs. Thus, for small Fermi energies, 2E F < ω 0 , although the three transitions can occur, the total nonlinear susceptibility nearly cancels out. For large Fermi energies, 2E F > 3 ω 0 , all three of these transitions are Pauli-blocked and there is only a non-resonant small intraband contribution. 11,13 At intermediate Fermi energies, ω 0 < 2E F < 3 ω 0 , it is possible to increase THG by, for example, Pauli blocking the one-photon and two-photon transitions, so that only Figure 2. Characterization of the third-harmonic signal. a. Z-scan measurements on four different regions of the sample (symbols). Gaussian fits to the data (curves) provide visual guides. b. Power scaling of the measured third-harmonic signal (symbols) in bare graphene and in gold-nanoribbon-graphene heterostructures for different nanoribbon widths with a monolayer h-BN spacer between the gold and graphene. The incident wavelength is 3.9 µm and the Fermi energy is EF ≈ 150 meV for this measurement. Linear fits to the data (solid lines) are used to determine the enhancement of the heterostructure relative to bare graphene. The dashed curve, which models the increase of the electron temperature with the incident light power, explains the observed saturation. c. Spectrum of third-harmonic generation from bare graphene, measured with 3.9 µm incident light. d. Third-harmonic signal of bare graphene and heterostructures measured as the input polarisation is rotated. Bare graphene (red) is isotropic with respect to the incident polarisation, while the gold nanoribbons (blue) result in a cos 6 θ dependence on the polarisation angle θ relative to the direction perpendicular to the ribbons (i.e., as expected from a third-order power scaling). e. Third-harmonic emission measured for fixed input polarisation when a polariser placed after the sample is rotated. For bare graphene (red), the third-harmonic signal is co-polarised with the input light, while for the heterostructures (blue) the third-harmonic polarisation is always orthogonal to the nanoribbons. In both cases, the polarisation of the third-harmonic signal is coherent relative to the incident light, which is indicated by the grey cos 2 θ line. For the data presented in panels d and e, the incident light wavelength is 5.5 µm and the Fermi energy is EF ≈ 150 meV. The studied heterostructure has an Al2O3 spacer of s = 5 nm, and ribbons with W = 200 nm, and g = 50 nm.
the three-photon transition is allowed and the other two transitions no-longer cancel it out. For low electron temperatures the gate response is expected to result in several sharp features as the system is tuned in and out of resonance. However, thermal broadening turns these features into broad shoulders. 11,13 The observed THG signal in bare graphene is plotted in Fig. 3b as a function of E F for four different incident wavelengths λ 0 = [5.5, 5.0, 4.5, 4.1] µm ([0.225, 0.248, 0.276, 0.302] eV). A prominent peak in the THG is observed that emerges at larger E F for shorter wavelengths and thus larger ω 0 . This is more clearly shown in Fig. S5d of the Appendix, where we plot the Fermi energy at which the THG is maximised versus the incident photon energy. As we show in the Appendix, the exact location of the maximum Fermi energy is affected by the electron temperature, which depends on the incident fluence. Importantly, these gating data show that we can actively modulate the THG signal in bare graphene. In particular, over the four data sets for bare graphene presented in Fig. 3b, we achieve an average intensity modulation of the THG by a factor ≈ 10 ± 3. As illustrated in Fig. S5 of the Appendix, the high electron temperature (which we estimate to be ≈ 1100 K in the Methods Section) is the main limitation on this modulation.
Similar gating measurements on the heterostructure for the geometry in which we obtain a maximum field x(ω) E x(3ω) /E0, relevant for third-harmonic generation. b. Gate-dependence of the third-harmonic signal in bare graphene for various incident light wavelengths. c. Gate-dependence of the third-harmonic signal in heterostructures with W = 200 nm, g = 50 nm and s = 5 nm, and incident light wavelengths of 3.9 µm and 5.5 µm (green and red symbols, respectively). The gray curve is the measured resistance of the the sample, the peak of which indicates the charge neutrality point. enhancement (W = 200 nm, g = 50 nm) are shown in Fig. 3c. Once again, we observe the expected shift as a function of the energy of the incident photons. We can modulate the THG by a factor of 7.4 ± 0.2 with a 3.9 µm incident light wavelength, and by a factor of 4.7 ± 0.2 with 5.5 µm. We stress that this active tunability comes from the unique combination of atomic thickness and linear electronic dispersion in graphene, which cannot be achieved in standard noble metal plasmonics. 3

PLASMON MEDIATION
Interestingly, the THG signal for a larger range of E F and for structures with smaller nanoribbon widths reveals an intriguing gate response. We plot the THG gate response of a heterostructure with W = 55 nm and g = 45 nm for different incident wavelengths in Fig. 4a-c. When the incident wavelength is 5.5 um, two peaks and a dip are clearly present in the data (Fig. 4a). As expected, when the wavelength (energy) is decreased (increased) these features shift to higher Fermi energies. The dip is still visible with an incident wavelength of 5.0 um at 2E F / ω 0 ≈ 4 (Fig. 4b). However, it shifts beyond our accessible gating range at a wavelength of 4.5 um (Fig. 4c).
None of these features are evident in the bare graphene data presented in Fig. 3b.
To provide an interpretation, we note that the graphene third-order nonlinearity associated with THG is determined by the interplay between the third-order nonlinear graphene conductivity σ (3) 3ω and the nonlinear field integral L E x (ω) 3 E x (3ω)dx (i.e., the cubic of the linear field is associated with the THG current amplitude, while the field at 3ω represents the emission enhancement produced by the heterostructure at the THG frequency). Tuning the Fermi level to higher energies can lead to the excitation of acoustic graphene plasmons under the metal 23 that also affects the nonlinear response. The third-order conductivity σ (3) 3ω depends on the excitation frequency ω, Fermi energy E F , the electron relaxation time τ e , and electron temperature T e associated with incident light absorption; the latter quantity can reach high values relative to the ambient room temperature that cause anomalous behaviour, such as shifts in the features of the optical conductivity as a function of doping, compared to constant T e (see S4 of the Appendix). Similarly, these parameters affect the nonlinear field integral through the linear graphene conductivity σ (1) , although in a different manner, thus leading to different field dependencies on E F in the gap and below the metal  It is thus important to stress that, as we show in the Fig. 4a-c, the third-order conductivity alone, despite its rich dependence on these parameters, cannot fully explain the observed THG signal. This is more evident in Figure S5 of the Appendix, wherein we plot the thirdorder conductivity at various electron temperatures, observing no features qualitatively similar to the dips in our experimental data. However, the peaks and dips appearing in the response at specific values of 2E F /hν correlate well with the integral of the electric fields, stemming in part from plasmonic interferences, and driving the nonlinear response. This integral, shown as a dashed line in Fig. 4a-c, is found to exhibit the same trend as the data.
From our simulations, we conclude that the nonlinear response in our experiment is mainly driven by the field in the gap region where it takes comparatively larger values (see Fig 4e). In particular, with an incident wavelength of λ = 5.5 µm, W = 55 nm, g = 45 nm, and s = 5 nm, we find a dip between the 3 ω 0 and 4 ω 0 transitions, and a peak just below the 4 ω 0 transition, at 3.45 × 2E F / ω 0 and ∼ 3.9 × 2E F / ω, respectively. The dip in the data can be explained by the partial cancellation of positive and negative complex components of the incident field in the gap by the field under the gold nanoribbons, which produces a reduction of the observed χ (3) .
In a more intuitive picture, at the dip we excite acoustic graphene plasmons under the metal ribbons. 23,24 As shown in Fig. 4d, for short lifetimes τ e < 100 fs, exciting an acoustic plasmon decreases the net integrated nonlinear field (dashed maroon curves in panels a-c, labeled 'integral'). This manifests as a dip in the nonlinear signal. In Fig. 4a-c, the model with τ e = 25 fs fits best to our data, confirming that the dips at 1.8 and 3.45 × 2E F / ω in Fig. 4a correspond to the excitation of acoustic plasmons. These plasmonic features do not appear in the W = 200 nm heterostructure data presented in Fig. 3c because the width is too large, and exciting acoustic plasmons in these structures would require more doping or a longer incident wavelength. As an additional experimental confirmation, in panel b and c we repeat the same gate-dependent THG measurements for higher incident photon energies 0.25 eV (5.0 µm) and 0.28 eV (4.5 µm). This shifts the acoustic plasmon resonance to higher gate voltages, causing the dips to move to higher Fermi energies in our simulations. This shift is a key signature of graphene plasmons, which we confirm experimentally, observing that the dip disappears for high incident phototon energies.
Finally, when graphene plasmons are explicitly suppressed in our simulations (by decreasing the plasmon lifetime τ e ), we are unable to reproduce the dips and peaks in our data (see the purple curve in Figure 4a), further confirming the excitation of plasmons in our graphene heterostructures.

CONCLUSION
Recent studies have reported a wide range in the estimate of χ (3) using third-harmonic generation. [10][11][12][40][41][42] By calculating an effective nonlinearity for the W = 200 nm heterostructure with a monolayer h-BN spacer, E F = 0.45 eV, irradiated with 5.5 µm light, and an electron temperature of approximately 1100 K, we find χ (3) exp ≈ 3.4 × 10 −6 esu, which is an order of magnitude larger than that of bare graphene, for which we obtain a maximum value of χ (3) exp ≈ 3.9 × 10 −7 esu (with E F = 0.39 eV and irradiated at 5.5 µm). Moreover, our all experimental measurements agree well with simulations based on the third-order nonlinear conductivity taken from References . 34,43 Unlike in metal plasmonics, we can actively modulate the nonlinearity of our graphene heterostructures by controlling E F with an external gate voltage. Graphenebased linear optical devices have already been shown to operate at GHz frequencies 44 and hence, our system provides a new route toward high-speed nonlinear optoelectronic switches and frequency converters. Additionally, our measurements reveal intriguing plasmonic effects supported by simulations in which graphene surface plasmon polaritons appear to directly modify the nonlinear optical response of our structures. These plasmonic excitations potentially provide a novel approach to the manipulation and amplification of light at subwavelength scales. In the present work, acoustic graphene plasmons seem to modulate the nonlinear response of graphene, while our simulations suggest that improving the plasmon lifetime by a factor of five would already increase the nonlinear response by one order of magnitude. We also observed that the maximum field enhancement (and hence the maximum nonlinearity) was obtained in a device that did not support acoustic graphene plasmons. However, one could in principle design different metal nanostructure geometries to simultaneously enhance the field and launch graphene plasmons. Engineering smaller nanostructures also has the potential to excite plasmons at shorter wavelengths. Our findings suggest that graphene plasmonic devices could provide unprecedentedly strong nonlinearities, potentially resulting in nonlinear optical effects at the single-photon level. . It reflects only the author's view, the EU Agency is not responsible for any use that may be made of the information it contains.

METHODS
Experimental details-We carry out our THG measurements using a modified z-scan setup, where the TH signal is measured while the sample is moved along the z axis through the focus of the laser beam (see Fig. 1d). Our incident light beam consists of linearlypolarised pulses of ∼ 260 fs duration with a tunable carrier wavelength of 5.5 µm (0.225 eV) at a 80 MHz repetition rate, which we create by performing Difference-Frequency Generation (DFG) between the signal and the idler beams of an Optical Parametric Oscillator (OPO). The typical acquisition time for a gate-dependent THG measurement is around 30 min; over this time we observe < 2 % power fluctuations. We use a Half-Wave Plate (HWP) to tune the polarisation of the incoming beam to that set by the polariser (pol). By rotating the HWP and the polariser together we can rotate the incident light polarisation relative to the orientation of the nanoribbons. Then, a lens with a 5.26 mm focal length focuses the incident light down to a waist of ≈ 20 µm for 5.5 µm nm light and ≈ 13 µm for 3.9 µm light. When the sample is moved parallel to the incident light beam (along the z axis), the nonlinear emission occurs most efficiently where the fluence is maximum (at the focal point). Afterwards, a lens with a 11 mm focal length collimates the beam, which is then sent through a bandpass filter (BP) to separate the THG signal from the excitation beam. Finally, the signal is coupled into a multimode fiber and sent to a largearea Superconducting Nanowire Single-Photon (SNSP) detector, with a ≈ 20 % detection efficiency at the thirdharmonic wavelength, 5.5/3 µm = 1.833 µm.
We verify the wavelength of this signal by removing all spectral filters, keeping the sample in the focus, and using a NIREOS GEMINI interferometer to perform Fourier transform spectroscopy on the signal. Sample fabrication and electrical doping-Here we use a heterostructure consisting of wet-transferred chemical-vapor deposition (CVD) graphene stripes shaped by dry ion etching (RIE) with a UV resist positive mask. The standard CVD graphene used in these samples has typical mobility below 1000 m 2 /Vs and ≈ 25 fs plasmon lifetimes, which is much shorter than those found in exfoliated graphene, where it can reach up to ≈ 500 fs at room temperature. 45 An insulator, which can consist of either wet-transferred CVD h-BN or Al 2 O 3 deposited by atomic layer deposition (ALD), defines the space between the substrate and the nanometer-thick gold nanoribbons, which are deposited by thermal evaporation onto a positive polymethyl methacrylate (PMMA) mask shaped using e-beam lithography. Note that for adhesive purposes, there is a 2 nm-thick Ti layer below the gold nanoribbons. These nanoribbons concentrate the electric field of a far-field incident light into graphene. 23 We perform the gating measurements applying a backgate voltage with a sourcemeter that allows us to monitor the current leakage between the Si and graphene layers, while increasing the applied voltage. A 1 mV voltage between the source and drain allows us to measure the graphene resistance, which we then use to estimate the induced E F . In order to calculate the induced E F given an applied voltage V , we consider the SiO 2 and graphene layers to behave as a parallel-plate capacitor, with the SiO 2 dielectric with relative permittivity d ≈ 4 and thickness d d in between. The conversion is given by with n = C · ∆V /e, where C = 0 d /d d is the surface capacitance and ∆V = V − V CNP is the difference between the applied voltage and the voltage at which the charge neutrality points (CNP) is found.
Extracting the third-order susceptibility-Experimentally we extimate χ (3) exp starting with the expression of the input (i) and output (o) average power as a function of the field as, where we assume laser pulses with repetition rate f , duration τ , waist W on the sample and Gaussian profile, and n ωi,o = 2.4 is the refraction index, and 0 and c are the permittivity and speed of light in vacuum. 11 Additionally, we can write the THG process as a function of the input and output fields as follows and solve for χ exp : where ω i = ω o /3 and d gr = 0.33 nm is the effective thickness of graphene. Electron temperature-In order to estimate the electron temperature in our samples, we have performed gating measurements on bare graphene for different incident wavelengths (see Fig. 3b) and determined the Fermi energy at which the maximum third-order susceptibility is found. As shown in Fig. S5d, this Fermi energy at which the maximum χ (3) is found shifts linearly with the photon energy of the excitation light. We fit the electron temperature of the simulations to best match the experimental data and determine a value of 1100 K, which is in good agreement with the simulations shown in Fig. S7. Note that, as mentioned in Ref., 11,13,46,47 the electron temperature changes with both the excitation wavelength and the Fermi energy. However, due to the small Fermi energy shift of about 80 meV, this effect was neglected in our simulations. For the power dependent measurements, it is necessary to describe the electron temperature as a function of fluence. To a good approximation, we assume a chemical potential µ > k B T e , which allows us to describe the electron temperature as follows: 48 T 0 is the ambient temperature, F 0 is the energy per pulse of the incident light and γ is the amount of energy absorbed that leads to hot electrons, which we have considered to be 3.5 × 10 −3 , which is consistent with previously reported values. 48 Considering the different electron temperatures for different excitation powers as discussed in Ref., 46 we are able to explain the power dependence of the TH signal (see Fig. 2b).
Simulating the third-order susceptibility-We write the third-harmonic susceptibility as where α 3ω is third-order polarisability given by where σ 3ω is the analytical third-order conductivity in graphene derived by Mikhailov in Ref., 43 and L is the length of the simulated region used for integration along the direction perpendicular to the nanoribbon width. Here, η 3 ω (x) is a dimensionless quantity representing the enhancement in electric field amplitude acting on the graphene layer relative to the incident field amplitude, so when multiplied by E 3 (ω i ) in Equation (4) it yields the actual field amplitude as a function of position x; likewise, η 3ω (x) represents the factor by which the THG field reaching the detector is modified by the presence of the structure surrounding the graphene. Equation (7) represents the contribution of the THG current to the far field. In particular, the η 3 ω (x) factor times σ 3ω is the THG current, which we represent as a polarisation density at the emission frequency 3ω through the continuity equation. Because of reciprocity, η 3ω (x) is exactly given by the enhancement amplitude in the near field relative to free space when the structure is illuminated with 3ω radiation. More precisely, we calculate η 3ω (x) as the complex factor of enhancement relative to the incident field in the field acting on the graphene layer when it is illuminated by a 3ω plane wave impinging from the detector direction. The electric field enhancements η ω (x) and η 3ω (x) entering the above expression are obtained using a RCWA Matlab script 49 and adapted to include graphene as an interface material adopting the nonlocal 2D linear optical conductivity of graphene σ(Q, ω) 20 that depends on the chemical potential µ and electronic temperature T e . We note that these enhancement factors depend on both the geometry of the heterostructures and the T eand µ-dependent linear conductivity. It requires periodicity of the structure perpendicular to the layered dimension. Here we consider the influence of electronic heating by the incident light pulse in the optical response of the graphene-metal hybrid system. In particular, following the procedure in Ref., 50 we make use of the implicit relation between E F , T e , and µ obtained from conservation of doping charge, along with the graphene heat capacity where F is the energy of the excitation pulse absorbed (i.e., F = η a F 0 where η a is the absorbed fraction of pulse energy F 0 ). The absorbed energy into graphene was obtained from the linear RCWA simulations at an ambient temperature of 300 K and β = 2 π ∞ 0 dxx 2 e x+µ/kBTe + 1 −1 The linear conductivity can thus be determined directly and the temperature dependent nonlinear conductivity can be compute by using Maldague's identity: .       The simulations in Fig. 3a and 4 of the main text are achieved with the Matlab script from Ref.. 49,51 The geometry of the 2D simulation is given by the layer stack of the original chip, which is composed by a 525 µm-thick silicon chip coated with a 285 nm-thick SiO 2 layer above and below. The monolayer graphene placed on the front side of the chip is modelled by a surface nonlocal conductivity, given in Ref.. 20 The spacer layer built by h-BN or Al 2 O 3 with thickness s was placed on top of the graphene. Finally a 10 nm-thick gold nanoribbon with width W was placed on top of the spacer. In a good approximation, the electric permittivity for silicon was assumed to be 12 and wavelength independent. The wavelengths depending dielectric permittivity of SiO 2 , h-BN and gold are found in Ref., 52 Ref. 53 and Ref., 54 respectively. As mentioned in the Methods, the simulation script requires a periodic structure in the lateral dimension. Here the periodicity is given by the sum of the nanoribbon width and the gap in-between. The electric field of the excitation incident light (THG signal) in the graphene layer was calculated for a normal incident of a monochromatic plane wave polarized parallel to the x axis and angular frequency ω (3ω). From the electric fields E x (ω) and E x (3ω), |χ sim | was calculated as described in the Methods section. For accessing the electron temperature T e in the graphene layer, we calculate the amount of absorbed energy in the graphene. To do so, we run the simulations with the RCWA, considering an ambient temperature of 300 K (see Fig.  S7a-c. Solving the implicit Eq. 8 − 10 of the main text leads to the electron temperature used in Fig. 4. These temperatures were considered in the third-order nonlinear conductivity and are shown in Fig. S7d-f. During the gating measurements, when the voltage is applied between the graphene and silicon substrate, a leakage current through the SiO 2 substrate can occur. Because of that, the charge densities in the graphene layer may differ slightly from those determined via the capacitor model introduced in the Methods section. This difference depends on the measurement time and the leakage. In order to compensate for this effect, the gating voltage of the experimental data shown in Fig. 4a-c of the main text has been fitted accordingly.