Coherent anti-Stokes Raman Spectroscopy of single and multi-layer graphene

We report stimulated Raman spectroscopy of the G phonon in both single and multi-layer graphene, through Coherent anti-Stokes Raman Scattering (CARS). The signal generated by the third order nonlinearity is dominated by a vibrationally non-resonant background (NVRB), which obscures the Raman lineshape. We demonstrate that the vibrationally resonant CARS peak can be measured by reducing the temporal overlap of the laser excitation pulses, suppressing the NVRB. We model the observed spectra, taking into account the electronically resonant nature of both CARS and NVRB. We show that CARS can be used for graphene imaging with vibrational sensitivity.

S ingle-layer graphene (SLG) has a high nonlinear third-order susceptibility: |χ (3) |~10 −10 e.s.u. for harmonic generation 1 and |χ (3) |~10 −7 e.s.u. for frequency mixing 2 , where one electrostatic unit of charge (1 e.s.u.), in standard units (SI) is 3 : χ (3) (SI)/χ (3) (e.s.u.) = 4π/(3 × 10 4 ) 2 . This is up to seven orders of magnitude greater than those of dielectric materials such as silica (χ (3) = 1.4 × 10 −14 e.s.u. 4 ). This property is due to optical resonance with interband electronic transitions 5 and has led to the observation of gate-tunable third-harmonic generation 1,6 and nonlinear four-wave mixing 2,7,8 (FWM, i.e., the third-order processes whereby an electromagnetic field is emitted by the nonlinear polarization induced by three field-matter interactions 3 ). FWM can be exploited for graphene imaging, with an image contrast of up to seven orders of magnitude 2 higher than that of optical reflection microscopy 9 . However, FWM-based imaging reported to date in graphene 2 lacks chemical selectivity and does not provide the same wealth of information brought about by the vibrational sensitivity of Raman spectroscopy 10,11 .
Coherent anti-Stokes Raman scattering (CARS) [12][13][14][15] is a FWM process that exploits the nonlinear interaction of two laser beams, the pump field E P at frequency ω P and the Stokes field E S at frequency ω S < ω P , to access the vibrational properties of a material. As shown in Fig. 1a, when the energy difference between the two photons matches a phonon energy (ℏω P − ℏω S = ℏω v ), the interaction of the laser pulses and the sample results in the generation of vibrational coherences 4 , at variance with impulsive anti-Stokes spontaneous Raman (IARS) which generates vibrational population [16][17][18][19] . While spontaneous Raman (SR) scattering is an incoherent signal 20 , since the phases of the electromagnetic fields emitted by individual scatterers are uncorrelated 20 , in CARS, atomic vibrations are coherently stimulated, i.e., atoms oscillate with the same phase 4 , potentially leading to a signal enhancement of several orders of magnitude depending on incident power and scatterer density 21,22 .
The same combination of optical fields used for CARS can generate another FWM signal, a nonvibrationally resonant background (NVRB) 2 , Fig. 1b. In both processes, the optical response consists of a field emitted at the anti-Stokes frequency 4 ω as = 2ω P − ω S . However, the interference of the two effects usually generates an additional contribution which is dispersive with respect to the emitted optical frequency, i.e., shaped as the first derivative of a peaked function (resembling the real part of the refractive index around a resonance), which introduces an asymmetric distortion of the Raman peak profile in the region 23 ω as = ω P + ω v .
In the biological field 21,24 , a wealth of studies has demonstrated the potential of CARS for fast imaging 21,22,25 , with pixel acquisition times as low as 24~0 .16 μs, thus allowing for video-rate microscopy 24 . By contrast, there are only a few reports to date of CARS imaging of micro-structured materials (such as polyethylene blend 26 , multicomponent polymers 27 , and cholesterol micro-crystals 28 ) and nanostructured ones (patterned gold surfaces 29 , single wall nanotubes 18,30 , and highly oriented pyrolytic graphite 31 ). Such studies, performed with pixel acquisition times down to 32~2 μs, have shown the ability of CARS to identify chemical heterogeneities on submicrometer scales and characterize single particles that are part of a larger domain, enabling, e.g., to visualize microscopic domains (polystyrene, polymethyl methacrylate, and polyethylene terephthalate) in the case of the above mentioned polymer mixtures 33 , or to provide detailed maps of microcrystal orientation in organic matrices (e.g., cholesterols in atherosclerotic plaques 28 ).
In graphene, despite the large 1,2 χ (3) , no CARS peak profiles, equivalent to those measured in SR, have been observed to date, to the best of our knowledge. We previously reported SR with single-color pulsed excitation 34 , using the same picosecond lasers usually adopted for CARS 24 . However, in order to measure CARS, a combination of pulses with different colors must be adopted 35 . By scanning the pulse frequency detuning in a twocolor experiment, a dip has been observed 36 in the third-order nonlinear spectral response of SLG at the G-phonon frequency. This was interpreted as an anomalous antiresonance and phenomenologically described in terms of a Fano lineshape 36 .
Here, we use two 1 ps pulses (see inset of Fig. 2) to explore FWM in SLG and few-layer graphene (FLG). We experimentally demonstrate and theoretically describe how the inter-pulse delay, ΔT (Fig. 1c, d) can be used to modify the relative weight of CARS and NVRB components that simultaneously contribute to the FWM, thus recovering the G-band Raman peak profile. We show that the dip in the nonlinear optical response around the vibrational resonance is due to the interplay of CARS and NVRB under electronically resonant conditions, which allows vibrational imaging with signal levels as large as those of the third-order nonlinear response.

Results
Sample preparation and SR characterization. SLG is grown on a 35 μm copper (Cu) foil, following the process described in ref. 37 . The substrate is heated up to 1000°C and annealed in hydrogen atmosphere (20 sccm) for 30 min at~200 mTorr. Then, 5 sccm of methane (CH 4 ) are let into the chamber for the following 30 min to enable growth 37,38 . The sample is then cooled back to room temperature in vacuum (~1 mTorr) and unloaded from the chamber. SLG is subsequently transferred on a glass substrate by a wet method. Polymethyl methacrylate (PMMA) is spin coated on the SLG/Cu and floated on a solution of ammonium persulfate (APS) and deionized water. When Cu is etched 38,39 , the PMMA Interaction with pulses ω P , ω S , results in blue-shifted (a) CARS and (b) NVRB contributions at ω as = 2ω P − ω S . Since in CARS a vibrational coherence is stimulated by two consecutive interactions with the pump and Stokes fields, their frequency difference must correspond to a Raman active mode, ω P − ω S = ω v . c, d Constraints for the temporal sequence of the field-matter interactions (represented by circles at the top of the pulse envelopes), for CARS and NVRB. In NVRB, the three interactions generating χ (3) happen within the few fs electronic dephasing time 20 . In CARS, the third interaction can occur over the much longer vibrational dephasing time (a few ps) 20 , within the pump pulse (PP) temporal envelope membrane with attached SLG is then moved to a beaker with deionized water for cleaning APS residuals. The membrane is subsequently lifted with the target substrate. After drying, PMMA is removed in acetone leaving SLG on glass. SLG is characterized by SR after transfer using a Renishaw InVia spectrometer at 514 nm. The position of the G peak, Pos(G), is~1582 cm −1 , while FWHM(G)~14 cm −1 . The 2D to G peak area ratio, A(2D)/A(G), is~5.3, indicating a p-doping after transfer~250 meV 40,41 , which corresponds to a carrier concentration~5 × 10 12 cm −2 . FLG flakes are produced by micromechanical cleavage from bulk graphite 42 . The bulk crystal is exfoliated on Nitto Denko tape. The FLG G peak is~1580 cm −1 . The D peak is negligible. The 2D peak shape indicates this is Bernal-stacked FLG 10,11 .

CARS measurements setup.
For CARS experiments, we use a two-modules Toptica FemtoFiber Pro source, with two erbiumdoped fiber amplifiers at~1550 nm generating 90 fs pulses at 40 MHz, seeded by a common mode-locked Er:fiber oscillator 43 , Fig. 2. In the first branch (FemtoFiber pro NIR), 1 ps pulses at 784 nm (pump pulse, PP) are produced by second-harmonic spectral compression 44 in a 1 cm periodically poled lithium niobate (PPLN) crystal. In the second branch (FemtoFiber pro TNIR), the amplified laser passes through a nonlinear fiber, wherein a supercontinuum (SC) output is generated. The SC spectral intensity can be tuned with a motorized Si-prism-pair compressor. A PPLN crystal with a fan-out grating (i.e., a poling period changing along the transverse direction) is exploited to produce broadly tunable (from 840 to 1100 nm) narrowband 1 ps Stokes pulses (SP), with a power <10 mW. A dichroic mirror is used to combine the two beams, whose relative temporal delay is tuned with an optical delay line. A long-working distance 20× objective (O, numerical aperture NA = 0.4) focuses the pulses onto the sample (Sa). The Stokes power is less than 2.45 mW during the scan, while the pump power is 2 mW for SLG and 6 mW for FLG. We note that no light-induced damage of SLG occurs up to 13.5 mW as verified by SR under the same repetition rate and similar pulse durations 34 . Thus, we rule out any structural degradation by the <6 mW pulses used here. The generated and transmitted light is collected by a condenser (C) and the PP and SP are filtered out by a short-pass filter (F). The total FWM signal is collected with an optical multichannel analyzer (OMA, Photon Control SPM-002-E). A dichroic mirror reflects the SP in order to measure its intensity (I s ) with a powermeter (P). FWM spectra are obtained by scanning ω S around ω P − ω v (at fixed ω P ) to probe the G-band phonon frequency, ω v = ω G . To assess data reproducibility we repeated the CARS measurements (scanning ω S ) finding no appreciable changes.
CARS by time-delayed FWM. Figure 3 displays the FWM intensity, normalized to I s , for different ΔT. In both SLG and FLG measurements, for ΔT shorter than the vibrational dephasing time τ~1 ps 45 , i.e., the characteristic time of coherence loss 20 , a Lorentzian dip at ω as = ω P + ω G appears on top of a background 36 . For ΔT > 2ps, while the total FWM signal decreases by nearly two orders of magnitude, the dip observed for FLG at ΔT 0 ps evolves into a Raman peak shape at the G-phonon energy.
No dispersive features are seen at any ΔT, unlike what normally expected for the interference between NVRB and CARS 23 . Here, we use pulses with duration δt~1 ps, since this allows us to scan the inter-pulse delay across the vibrational dephasing time τ to suppress the NVRB cross section more than the vibrational contribution, while minimizing the spectral broadening due to the finite pulse duration 1=δt $ 15 cm À120 .
Since both CARS and NVRB depend quadratically on the number of scatterers 33 , the SLG intensity is significantly reduced with respect to FLG (Fig. 3), with a lower signal-to-noise ratio, hampering the detection of peak-shaped vibrational resonances expected for ΔT > 1.4 ps.
The data in Fig. 3 can be qualitatively understood as follows. The anti-Stokes signal, I(ω as ), generated by CARS and NVRB, is proportional to the square modulus of the electric field emitted by the third-order polarization 4 , P (3) , as: CARS and NVRB are simultaneously generated by two lightmatter interactions with the PP and a single interaction with the SP, with different time ordering, Fig. 1a, b. Consequently, P ð3Þ / E 2 P E Ã S , where * indicates the complex conjugate. Therefore, the FWM signals are quadratic with respect to the pump power and linear with respect to the Stokes power. However, the temporal constraints for such interactions are significantly different for the two cases. As shown in Fig. 1c, d, in the case of NVRB the three interactions must take place within the dephasing time of the involved electronic excitation, which in SLG is~10 fs 46,47 , i.e., much shorter than the pulses duration (δt~1 ps). Hence, P ð3Þ NVRB (ω as ) is only generated during the temporal overlap between the two pulses P ð3Þ NVRB / E 2 P ðt À ΔTÞE Ã S ðtÞ (the three field interactions, in a representative NVRB event, are indicated by three nearly coincident dots in Fig. 1d). In CARS, the electronic dephasing time only constrains the lag between the first two interactions that generate the vibrational coherence (the two stimulating-field interactions are represented by the two nearly coincident dots in Fig. 1c). This can be read out by the third field interaction within the phonon dephasing time, τ~1 ps 45 (indicated, for a representative CARS event, by the third dot in Fig. 1d). Thus, P 48 . Therefore, ΔT can be used to control the relative weights of P Electronically resonant and non-resonant FWM. The system response can be evaluated through a density-matrix description 20 of P ð3Þ ðω; ΔTÞ. In the presence of a temporal delay between PP and SP, their electric fields can be written as 3 : E P ðt; ΔTÞ ¼ A P ðt; ΔTÞe Àiω P t and E S ðt; 0Þ ¼ A S ðt; 0Þe Àiω S t , where A P/S (t, ΔT) indicates the PP/SP temporal envelope. By Fourier transform, the fields can be expressed in the frequency domain asÊ P ðω; ΔTÞ ¼ R þ1 À1 E P ðt; ΔTÞe iωt dt and E S ðω; 0Þ ¼ R þ1 À1 E S ðt; 0Þe iωt dt, which can be used to calculate P ð3Þ CARS ðω; ΔTÞ as 20,53 where η CARS = n CARS μ ba μ cb μ cd μ ad , μ ij is the transition dipole moment between the i and j states, n CARS is the number of scatterers involved in the CARS process, ω ij ¼ ω ij À iγ ij ¼ ω i À ω j À iγ ij is the energy difference between the levels i and j, and γ ij ¼ τ À1 ij is the dephasing rate of the |i〉 〈j| coherence 20 ; a and c denote the vibrational ground state |g, 0〉, and the first vibrational excited level, |g, 1〉, with respect to the electronic ground state |g〉 (π band). In our experiments, c corresponds to the G phonon at q~0, b and d indicate the vibrational ground state, |e, 0〉, and the first vibrational excited level, |e, 1〉, with respect to the excited electronic state |e〉 (π * band).
Using the conservation of energy represented by the δdistribution in Eq. (2) and integrating over ω 2 , we get: Definingν ¼ ω=ð2πcÞ, the third-order nonlinear polarization can be expressed as a function of the Raman shiftν Àν P ð Þ as P ð3Þ ðω; ΔTÞ ¼ P ð3Þ ð2πcν; ΔTÞ.
The ω ca level in the denominator of Eq. (3) is the frequency of the Raman mode coherently stimulated in CARS, while ω ba and ω da are frequency differences between the electronic levels. In the case of real levels, resonance enhancement occurs 20 . In view of the optical nature of the involved phonons (q~0), and due to momentum conservation, only one electronic level must be included in the calculation and, consequently, the nonlinear response can be derived for ω ba = ω dc = ω P . In a similar manner, P ð3Þ NVRB can be expressed as 20 where η NVRB = n NVRB |μ ea | 4 , n NVRB is the number of scatterers involved in the NVRB process, and ω ea is the energy of the electronic excited level involved in the NVRB process. Since the cross section of third-order nonlinear processes in graphene is enhanced by increasing the photon energy 54,55 , we consider only the dominant case, i.e.,ν ea ¼ 2ν P . We describe the spectral FWM response assuming monochromatic fields with no inter-pulse delay:Ê P ðωÞ ¼ E P Á δðω À ω P Þ, E S ðωÞ ¼ E S Á δðω À ω S Þ. From Eqs. (3) and (4), the CARS and NVRB nonlinear polarizations can be expressed as 4,20 which can be used to calculate the total FWM spectrum according to Eq. (1). Figure 4a plots the electronically nonresonant case. The CARS polarization, defined by Eq. (5), is a complex quantity: the real part has a dispersive lineshape, while the imaginary part peaks at ω ca . The NVRB polarization, defined by Eq. (6), is a positive real quantity. Accordingly, the FWM spectrum in the electronically non-resonant condition, I(ω as ) NR , can be written as 20 Iðω as Þ NR $ jP and it can be significantly distorted by the third term in Eq. (7) depending on the relative weight of the two corresponding susceptibilities. The maximum of the signal, when the dispersive contribution is dominant, can be frequency shifted from the phonon frequency. This is the most common scenario, with the dispersive lineshapes hampering direct access to the vibrational characterization of the sample in terms of phonon frequencies  Þat different ΔT between the beams at tunable ω S and fixed ω P . In a, colored dashed lines are fits to the data using Eq. (1) and the nonlinear polarization obtained from Eqs. (19) and (20). Vertical black dashed lines indicate three energies ðν 1;2;3 Àν P ¼ 1545; 1576; 1607cm À1 Þ, taken as reference for the FLG CARS images in Fig. 5 and lifetimes. Such limitation is particularly severe when χ ð3Þ NRVB is comparable to χ ð3Þ CARS and the NVRB and CARS contributions have the same order of magnitude. This condition is common in the case of a weak vibrational resonant contribution μ ba μ cb μ cd μ ad jμ ea j 4

<<1
, as in the case of low concentrations of oscillators n CARS n NVRB

<<1
. Hence, this produces an intense NVRB signal and reduces the vibrational contrast, hindering the imaging of electronically nonresonant samples. This is the case for cells and tissues which need to be excited in the near infrared to avoid radiation damage 56 .
For SLG, the linear dispersion of the massless Dirac Fermions makes the response always electronically resonant. In the case of FLG, absorption has a complex dependence on wavelength, as well as on the number of layers and their relative orientation, exhibiting, e.g., a tunable band gap in twisted bilayer graphene 57 . This is also reflected in the resonant nature of SR 58,59 . However, approaching visible wavelengths, the absorption spectrum flattens above~0.8 eV and it is $ ð1 À πe 2 =2hÞ N for Bernal-stacked Nlayer graphene 60 . Our exfoliated FLG are Bernal stacked, as also confirmed by the measured 2D peak shape in SR 10,11 . Accordingly, at the typical CARS wavelengths used here (784 and 894 nm), SLG and Bernal FLG are electronically resonant, unlike the situation for most biological samples 56 . This results in an opposite sign in the CARS response, i.e., a spectral dip in =ðχ ð3Þ CARS Þ, related to two additional imaginary unit contributions in the denominator of Eq. (5), ðω P À ω ba Þ and ðω À ω da Þ, wherein the iγ ba , iγ da components dominate. Further, the −i contribution from ð2ω P À ω ea Þ results in a NVRB dominated by the imaginary part, as illustrated in Fig. 4d-f.
Thus, the third term in Eq. (7) must be replaced with the contribution from the interference of the spectral dip =ðχ , which depends only on the material under examination and not on the pulses used in the experiment. Such a qualitatively different interplay between NVRB and CARS, compared with the experimental lineshapes for ΔT = 0 in Fig. 3, unambiguously indicates the presence of electronic resonance in SLG and Bernal FLG. For a given material, the relative weight of the two FWM contributions can be modified by using pulsed excitation and tuning the temporal overlap between PP and SP fields 48 , i.e., changing ΔT. The experimentally observed evolution of the FWM signal in FLG as a function of PP-SP delay in Fig. 3 has a trend similar to that shown in Fig. 4e, f as function of η NRVB /η CARS , validating the resonance-dominated scenario.
A more quantitative picture can be derived from Eqs. (3) and (4), where the PP and SP temporal profiles are taken into account,  7) and (8). In c, f, selected spectra corresponding to three η NVRB /η CARS from the colormap are reported matching those retrieved from the experimentally measured autocorrelation (Fig. 2).
As model parameters for the FLG we use the experimental SR ν ca ¼ 1580 cm À1 , with fitted τ ca = 1.1 ± 0.1 ps 10,45 (corresponding to FWHM(G) = 10 cm −1 ) and τ da = τ ba = τ ea = 10 ± 2fs in agreement with the value measured for SLG 47 . The ratio between NVRB and CARS contributions η CARS η NVRB ¼ ð3:0 ± 0:7Þ 10 À5 is obtained by fitting the data in Fig. 3 with Eqs. (1), (19), and (20). The resulting spectra (colored dashed lines in Fig. 3), evaluated by tuning only the PP-SP delays, are in good agreement with the experimental data, with η CARS η NVRB as the only adjustable parameter. This ratio, combined with Eqs. (5) and (6), allows us to extract the ratio between the third-order nonlinear susceptibilities for CARS and NVRB, The dependence of our spectral profiles on ΔT indicates that the peculiar FWM lineshapes for SLG and FLG originate from the interference between two electronically resonant radiation-matter interactions (NVRB and CARS) rather than from a matter-only Hamiltonian coupling the electronic continuum and a discrete phonon state (implying a resonance between the corresponding energies), resulting in the Fano resonance 61 suggested in ref. 36 .
Coherent vibrational imaging. In the electronically nonresonant case, CARS provides access to the real part 23 of χ (3) . However, due to the dispersive nature of the χ (3) real part 23 , it distorts the phonon lineshapes 3 , unlike SR. In SLG and FLG the FWM signal arises from the imaginary (non dispersive) CARS susceptibility, and it is amplified by its NVRB (third term in Eq. (8)). Then, the signal can be used for vibrational imaging, unlike the nonresonant case, for which the real part of the CARS susceptibility is involved, generating spectral distortions. Thus, coherent vibrational imaging can be performed without suppressing the NVRB, achieving signal levels as large as those of FWM, while preserving the Raman information.
The vibrationally resonant contribution I can be isolated by subtracting from the I 2 FWM signal atν 2 Àν P $ν G , the NRVB obtained by linear interpolation of the spectral intensities measured at the two frequencies at the opposite sides of vibrational resonance where the indexes i = 1, 2, 3 refer to data atν 1 ¼ν P þ 1545 cm À1 , ν 2 ¼ν P þ 1576 cm À1 ,ν 3 ¼ν P þ 1607 cm À1 (i.e., withν 2 near to the G-phonon frequency and jν 1;3 Àν G j greater than two halfwidths at half maximum of the measured profiles, as shown in Fig. 3). This combination of electronically resonant NVRB and CARS nonlinear responses gives CARS images (i.e., retaining vibrational sensitivity) with signal intensities comparable to those of NVRB, for which sub-ms pixel dwell times have been demonstrated with the use of a point detector, e.g., a photomultiplier 2 . In our case, the images in Fig. 5 are obtained with a pixel dwell time~200 ms using a Si-charge-coupled device array, aiming at a complete spectral characterization, and scanning the sample at fixed ΔT with stepper-motor translation stages. Figure 5a-c displays nonlinear optical images measured at two different ω S , corresponding to vibrationally nonresonant and resonant conditions. Extracting for each image pixel I 1 (Fig. 5a), I 2 (Fig. 5c), and I 3 , required for Eq. (9), we obtain an image with suppression of the signal not generated by FLG, as in Fig. 5e.
To obtain a quantitative comparison of the different images, we plot the pixel intensity histogram in Fig. 5b, d, f. This gives a bimodal distribution: one peak corresponds to the most intense pixels, associated with FLG (I g ) and the other is related to the weaker substrate signal (I s ). The ability to discriminate sample from substrate can be quantified in terms of 1) I g compared to I s , evaluated as the difference I g − I s , and 2) the proximity of I s to I = 0 in the histogram (dashed black line in Fig. 5b, d, f). These two features can be quantified by the contrast in order to compare the images: C ¼ I g ÀI s I s , where I g and I s are the mean FLG and substrate intensities, corresponding to the local maxima in the histograms in Fig. 5b, d, f. In Fig. 5g, h, we plot the intensity profiles along two scanning paths, one inside (dashed) and the other adjacent to (full line) the FLG flake.
Comparing the three histograms (Fig. 5b, d, f), the vibrationally off resonant FWM image (NVRB only, Fig. 5b) has the highest I g . The visibility of the flakes is limited by the noise of the detector and by the substrate χ (3) . NVRB, lacking vibrational specificity, can also originate from the glass substrate outside the FLG flake I s ) 0 À Á , as indicated by the scanning profile in Fig. 5h (red line). This may become a critical limitation in those substrates with χ (3) much larger than Si (|χ (3) |~2.5 × 10 −10 e.s.u. 33 ), such as Au (χ (3) = 4 × 10 −9 e.s.u. 62,63 ). Similarly, the vibrationally resonant FWM, I 2 , originating from concurrent CARS and NVRB processes (Fig. 5d), has a I s ) 0, related to NRVB. The depth of the FWM dip ( Fig. 5f) is related to the CARS signal intensity, and its vibrational sensitivity brings a substantial contrast increase, as demonstrated by the close-to-zero average value of the (green) scanning profile in Fig. 5h.
In summary, by using an experimental time-delayed FWM scheme, CARS peaks equivalent to those of SR were obtained from graphene. By explaining the physical mechanism responsible for the FWM signal, we demonstrated that the spectral response can be described in terms of joint CARS and NVRB contributions concurring to the overall signal. Unlike nonresonant FWM, where dispersive lineshapes hamper vibrational imaging of biological systems, the resonant nature of FWM in graphene, which can be traced back to its peculiar electronic properties, mixes CARS and NVRB, resulting in Lorentzian profiles which are either peaks or dips depending on their relative strength. We also demonstrated that CARS can be used for vibrational imaging with contrast equivalent to spontaneous Raman microscopy and signal levels as large as those of the third-order nonlinear response.

Data availability
All relevant data and Matlab codes are available from the authors.