High-harmonic spectroscopy of low-energy electron-scattering dynamics in liquids

High-harmonic spectroscopy is an all-optical nonlinear technique with inherent attosecond temporal resolution. It has been applied to a variety of systems in the gas phase and solid state. Here we extend its use to liquid samples. By studying high-harmonic generation over a broad range of wavelengths and intensities, we show that the cut-off energy is independent of the wavelength beyond a threshold intensity and that it is a characteristic property of the studied liquid. We explain these observations with a semi-classical model based on electron trajectories that are limited by the electron scattering. This is further confirmed by measurements performed with elliptically polarized light and with ab-initio time-dependent density functional theory calculations. Our results propose high-harmonic spectroscopy as an all-optical approach for determining the effective mean free paths of slow electrons in liquids. This regime is extremely difficult to access with other methodologies, but is critical for understanding radiation damage to living tissues. Our work also indicates the possibility of resolving subfemtosecond electron dynamics in liquids offering an all-optical approach to attosecond spectroscopy of chemical processes in their native liquid environment.


I. INTRODUCTION
High-harmonic generation (HHG) is an extremely nonlinear process that occurs when a strong laser field interacts with gaseous, solid, or liquid targets.It results in an up-conversion of photon energies up to the tender X-ray regime [1], and has been well established as a highly versatile table-top source of attosecond pulses [2][3][4][5][6][7][8].HHG has also led to a new branch of attosecond spectroscopy that relies on extracting dynamical information directly from measured spectra, known as high-harmonic spectroscopy (HHS).This technique has already enabled imaging of molecular orbitals [9][10][11], the reconstruction of charge migration [12] and time-dependent chirality [13], as well as tunneling-ionization dynamics [14].While the majority of these applications are based on gas-phase HHG, solid-state HHG has recently attracted considerable attention because of its potentially higher efficiency and access to ultrafast dynamics and light-driven phase transitions in condensed matter [15][16][17][18].A prerequisite for an accurate interpretation of the underlying dynamics from the measured high-harmonic spectra lies in the * These five authors contributed equally formulation of a broadly applicable theoretical model.
In the gas phase, this understanding and modelling is often based on the semi-classical three-step model (TSM) [19], or its quantum-mechanical extension [20], which describes HHG as a set of electron trajectories initiated by a tunneling process.This approach is usually in good agreement with full ab-initio calculations, and allows interpretation of dynamical information in HHS [21,22].A hallmark of the model is that it correctly predicts the HHG cut-off and its quadratic dependence on, both, electric-field amplitude and wavelength, connected to the most energetic returning electron trajectory [23,24].For HHG in crystalline solids a similar trajectory picture, based on the material's band structure, can in principle be applied in momentum space (after applying the Bloch theorem) [16,[25][26][27][28][29][30][31].The predicted cut-off energy was shown to scale linearly with both field amplitude and wavelength [16,17,25,30,[32][33][34][35].However, there still remains some debate about the scaling based on the active HHG mechanisms in different solid systems [18,[36][37][38].Moreover, the direct comparison of crystalline and amorphous solids of the same composition (quartz vs. fused silica) has shown that the cut-off energy is much lower in the latter under the same driving fields [39].Whereas this observation remains to be fully explained, it points to the importance of long-range order in condensed-phase HHG, which was also explored for one-dimensional models [40,41].In contrast to gases and crystals, HHG in the liquid phase is far from being well understood.This is because it is challenging to: (i) experimentally measure HHG spectra from bulk liquids, (ii) numerically simulate strong-field physics in liquids, and (iii) there is still no intuitive model that describes non-perturbative lightdriven dynamics in liquids.Liquids, therefore, present a unique case where neither gas (single isolated particle approach) nor solid-state (Bloch theorem and periodic boundary conditions) approaches are strictly applicable.This gap in knowledge limits potential applications to ultrafast spectroscopy that are especially appealing in liquid targets.Only very recently, HHG in bulk liquids have been demonstrated beyond the visible domain utilizing the flat-jet approach [42][43][44].However, fundamental questions about the dominant HHG microscopic mechanisms, the cut-off scaling with wavelength or the macroscopic effects still remain unanswered.As the major biochemical processes take place in a liquid environment, detailed experimental results and the development of theoretical tools capable of describing the HHG process are crucial for understanding the electron dynamics in liquids.We note that our work addresses HHG at typical intensities of ∼ 10 13 W/cm 2 , which is the basis of HHS, in contrast to HHG in the coherent-wake-emission regime taking place at intensities beyond 10 17 W/cm 2 .The latter has been demonstrated on the surface of liquids [45].As a consequence of the broken inversion symmetry, both even and odd harmonics were observed in those experiments.Our present experiments, in contrast, probe the bulk of the liquid phase, such that no even harmonics are observed.
Here, we experimentally measure high-harmonic spectra from liquid water and alcohols over a broad range of laser wavelengths.We observe that the HHG cut-off energy (E c ), i.e. the energy marking the end of plateau region, as defined in the Lewenstein formalism [20], is wavelength independent in strong contrast with the semiclassical TSM for gases [24], as well as some models for solid-state HHG [16, 25-27, 29-31, 33].This implies that potentially new mechanisms are relevant in liquid HHG, and that the structural arrangement of the liquid (i.e. the lack of long-range order) might play a crucial role in the dynamics.We investigate this experimental result with a combination of newly developed ab initio techniques, and introduce a semi-classical model for HHG in liquids.Our proposed model takes electron scattering into account and successfully reproduces the observed wavelength-independence of E c .We identify a key parameter in HHG from the liquid phase -the effective mean-free path (λ MFP ) -which we extract from measurements using the extended semi-classical model.Our results shed light on fundamental strong-field-driven processes in liquids, that are distinct from what happens in either gas phase or solid-state environments, and form the basis of a first intuitive picture of HHG in liquids.

A. Effect of electron scattering on HHG
One noticeable difference between HHG in dilute gases and in condensed phases is the significance of electron scattering in the latter.We therefore start our analysis by formulating a semi-classical real-space trajectory picture similar to the TSM of gas-phase HHG, but include scattering from the beginning.Within this picture, harmonic photons are emitted as a result of electron trajectories that recombine with their parent ion.We assume that an electron may be photo-excited to the conduction band of the liquid at any time, t ion during the laser cycle.Following this, Newtonian equations of motion can be analytically solved to obtain the electron trajectory along the laser polarization axis (x(t), given in atomic units): x(t) = qE 0 mω 2 [cos(ωt) − cos(ωt ion ) + ω(t − t ion ) sin(ωt ion )], (1) where E 0 is the laser field's peak amplitude, q and m are the electron charge and mass, and ω is its angular frequency.Recombining trajectories arise by demanding that the electron returns to its initial position at time, t rec , i.e. x(t rec ) = 0.This generates a set of short and long electron trajectories that are initiated at t ion , recombine at t rec , and have an emission energy of Ω = I p + 0.5( dx/dt) 2 .The resulting cutoff in the absence of scattering is then obtained semi-analytically as Ω cutoff = I p + 3.17U p , where U p = E 2 0 /4ω 2 is the classical ponderomotive energy.In the following, we consider only the trajectories for which the electrons return to their parent molecule.We extend this model to include scattering processes of electrons with neighboring molecules: as a simple approximation, we assume that any trajectory that exceeds a characteristic excursion length (denoted as l max ) will likely scatter and therefore not contribute significantly to the HHG emission (see Fig. 1A).That is, we neglect the contribution of scattered electrons to HHG.The electron could also recombine to other centers (as reported in the case of solids [46]).However, semi-classical calculations show that this leads to a higher-energy cutoff and a different behavior under elliptically polarized light than measured experimentally (see Extended data Fig. 3), which allows us to discard this channel as the dominating one for our driving conditions.This effectively translates to the following constraint for linearly polarized light -only recombining trajectories that uphold |x(t ion < t < t rec )| < l max emit harmonics.The set of trajectories that fulfills these equations of motion are identical to those in the gas-phase case, except that trajectories extending beyond l max are not included.As a consequence, this introduces a mechanism that modifies E c and its scaling, and also naturally reduces E c compared to the gas phase by ∆E c .
Typical results for the model are presented in Fig. 1(C) where the cutoff energy follows the standard TSM pre-diction for short wavelengths (where both short and long trajectories do not surpass l max ), but rapidly saturates around 800 nm where the cutoff trajectories in the gas phase exceed a distance of few angstroms.We will show below that this simple model reproduces the main features of both measurements and ab-initio calculations.We note that the non-scaling of the cutoff with wavelength is reproduced by the semi-classical model for any choice of l max , which only changes the maximal E c (see Figure 1(B)).Moreover, this behavior does not depend on the laser intensity (see Supplementary Material (SM), Section S4.C and Fig. S6).We emphasize that this simple picture likely does not capture the full physics of strong-field light-matter interactions in the liquid state.Nonetheless, the fact that the very peculiar cutoff behavior (compared to other phases of matter) is reproduced, is encouraging.We also note that some of the approximations utilized here might not be accurate in the liquid phase (e.g. the strong field approximation (SFA) or neglecting multi-center recombination), but: (i) corrections accounting for these effects can conceptually be added, and (ii) the characteristic physical behaviour of the cut off is independent of these approximations, as demonstrated by the ab-initio results shown in Section C.

B. Experimental results
The experimental setup is shown in Fig. 2(A).It consists of a laser system delivering ∼30-40 fs laser pulses with adjustable central wavelength (800-1800 nm) and a high-vacuum chamber containing the liquid flat-jet system and a flat-field imaging spectrometer.Further details are given in the Methods section and in the SM, Section S1.We have measured high-harmonic spectra of water (H 2 O) and ethanol (CH 3 CH 2 OH) from the liquid and gas phases of each species at different wavelengths.A typical background-corrected HHG spectrum of water is shown in Fig. 2(B).The liquid-and gas-phase spectra are recorded back-to-back to minimize drifts in the experimental parameters.Figure 2(B) presents HHG spectra for H 2 O with the top and bottom panels directly comparing the gas-and liquid-phase signals recorded with 800nm and 1500-nm drivers, respectively.The liquid-phase harmonics are roughly ten times brighter than the gasphase harmonics.Both spectra exhibit a distinct plateau, followed by a sharp cut-off region where the harmonic yield drops exponentially.Notably, the cut-off energy E c is around H9 in the liquid spectrum and H17 in the gas spectrum.To reduce uncertainty, we determine the cutoff energy following the formalism elaborated in section S3 of the supplement.In brief, as the harmonic yield in the cut-off region is expected to decay exponentially, the logarithmic value of the harmonic yield is fitted to a linear function of the harmonic energy.In contrast, the harmonic yield for the plateau harmonics remains constant as a function of energy.The intersection of these two lines (the linear fit of the log(yield) values as a func-tion of the harmonic energy in the cut-off region and the line indicating the average log(yield) value of the plateau harmonics) defines the cut-off energy or the end of the plateau region.We find that the liquid phase shows a much-reduced cut-off compared to the gas phase.For generation in H 2 O, the gas-to-liquid difference ∆E c is about 10 eV at 800 nm and about 26 eV at 1500 nm.This observation is a first hint at different dominant mechanisms in each phase of matter, which prevents the emission of higher energy photons from the liquid.
We next explore the wavelength scaling of E c .This basic property reveals information about the laser-driven electron dynamics in the liquid phase.Figure 2(C) and (D) show measured high-harmonic spectra from water and ethanol at two different wavelengths, respectively.All spectra display the characteristic envelope with a plateau and a sharp cut-off region.This allows us to define the cut-off energy E c as the intersection point of two lines that connect the plateau and cut-off intensities, respectively.The details of this procedure that is followed throughout this work are given in Section S3.For both liquids, all spectra share the same cut-off energy of the plateau, i.e.E c = 14.2 eV and E c = 11.4 eV in the case of water and ethanol, respectively.These results substantially differ from the gas-phase results, as well as the standard TSM, which both show that for the used laser intensities the cut-off should have extended by ∼25 eV between the 800 nm and 1500 nm drivers.
Notably, E c in ethanol is ∼3 eV smaller than in water.The difference in cut-off energies between these liquids is substantially larger than the difference in their band gaps (∼ 8 eV for H 2 O and ∼ 8.5 eV for ethanol [48,49]).This is a crucial point, since in the gas phase, and within the standard TSM, the cut-off should only vary by the difference of these value.The larger variation is indicative of the fact that the liquid structure, and more precisely, the electron dynamics in the liquid-phase, are playing an additional and yet to be specified role.
A critical aspect of HHG in liquids consists in ensuring that the measured signals originate from the bulk liquid phase.This requires explicitly excluding HHG emission from the evaporating gas phase, as well as HHG from the gas-liquid interface.A complete experimental separation of HHG emission from the gas and liquid phases has been achieved by using the wedge-like geometry in the upper part of the liquid jet, as shown in Extended Data Fig. 1.Additional experiments with the liquid jet placed at an angle of 45°with respect to the driver-beam propagation allow us to exclude HHG from the liquid-gas interface to contribute a measurable signal because of the absence of any even harmonics (see Extended Data Fig. 2).
We thus reach two main conclusions: (i) in the liquidgenerated high-harmonic spectra the position of the cutoff depends on the nature of the liquid sample, and (ii) the cut-off energy is wavelength independent, at least in water and ethanol.In what follows, we will show that these results are reproduced by ab-initio calculations.FIG. 3. TDDFT calculations in the supercell approach (A) A schematic representation of liquid water in supercell approach.There are 64 water molecules in a cubic cell at the experimental density of 1 gr/cm 2 ; periodic boundary conditions are implemented in 3D where each lattice vector is of length ∼ 12.43 Å. (B) HHG spectra calculated for liquid water using two different driving wavelengths.The peak intensity for all wavelengths is the same and it is equal to approximately 20 TW/cm 2 .These HHG spectra are averaged over 5000-5500 water molecules in the liquid phase (for details see the SI and Ref. [47]).TDDFT calculations in the cluster approach (C) Illustration of the computational approach with a cluster radius of ∼ 15.5 Å. (D) Wavelength scaling of high-harmonic spectra calculated for a constant peak intensity.

C. Numerical results: time-dependent density-functional theory
We now compare these experimental findings and the results of our simple model to two newly developed ab initio techniques for describing the strong light-matter response of liquids.Figure 3 presents simulated HHG spectra from liquid water that are based on a combination of well-established Car-Parrinello molecular dynamics (CPMD) [51] and time-dependent density-functional theory (TDDFT) [52] simulations in a periodic supercell including 64 water molecules at the experimental density of 1 g/cm 3 and temperature of 300 K (for details see Ref. [47] and the SM, Section S4.A).This is the first realistic, currently tractable description of HHG in liquids.TDDFT naturally includes mean-free path effects as it includes electron-electron and electron-ion scattering.Overall, very good agreement with the experimental results is observed, and most importantly, the cut-off energy and its wavelength independence are well reproduced in Figure 4.Moreover, a time-frequency analysis of the TDDFT results (see Extended Data Fig. 4) shows that only very short electron trajectories contribute to the HHG spectra, in agreement with our semi-classical model.Note that since the DFT-GGA underestimates the liquid-water band gap, the calculated HHG cutoff is about 1.5 eV lower than the experimental value.This numerical approach qualitatively reproduces the experimentally observed weak dependence of the cut-off on the laser intensity (see discussion in SM, Section S4.A).This further confirms that the above experimental findings are a signature of the microscopic mechanism in the liquid phase and not the result of macroscopic effects, which are absent in our theoretical modeling.This result is complemented by a second set of abinitio TDDFT calculations based on molecular clusters which employ some additional approximations (for details, see ref. [57]).The advantage of this approach is that it allows for faster calculations while still leading to accurate results; thus, it can be employed for a more detailed numerical study.Figure 3(B) shows simulation results of HHG in liquid water with the cluster approach for many wavelengths at a fixed laser intensity.Clearly, the same trend is observed and the cut-off is indepen-FIG.4. Scaling of the cut off with the density of the liquid (A) Measured high-harmonic spectra (1800 nm driver, 6 × 10 13 W/cm 2 ) as a function of vertical position on the liquid flat-jet, corresponding to the indicated change in density.The change in density has been determined from an absolute temperature measurement carried out by Raman thermometry [50], as shown in (B).(C) Spectra calculated with the ab-initio cluster approach (900 nm, 4 × 10 13 W/cm 2 ) for varying densities (∆ρ/ρ = -15%, -7%, 0%, 7%).The harmonics above the cut-off show a systematic decrease in yield with increasing density.(D) Maximal observed kinetic energy as a function of the inverse density for ethanol and the corresponding lmax obtained from the formula shown in the inset (derived in SM, Section S4.D).The vacuum electric field is corrected using the 1800-nm refractive index in ethanol, to obtain the peak electric field (E peak ) inside the liquid.
dent of the wavelength, at least in the range of 500 nm -1800 nm.In the SM, Section S4.B, we show that the cutoff with the cluster approach is similarly weakly dependent on the laser intensity (above a saturation intensity of 5 • 10 13 W/cm 2 ), and that the wavelength independence of the cut-off is maintained for other laser intensities, as well.With the cluster approach, we also performed calculations for two additional liquids (ammonia, NH 3 , which is polar, and methane, CH 4 , which is non-polar).In the SM, Section S4.B, we show that HHG calculations in liquid NH 3 and liquid CH 4 also predict the same wavelength independence of E c .These results, in combination with our measurements, lead to the conclusion that this characteristic non-scaling of E c is a fundamental general and unique property of the liquid-phase HHG, and applies both for polar and non-polar liquids.These accurate quantum-dynamical simulations reproduce and complement our experimental findings, which validates the broad applicability of our conclusions.

III. DISCUSSION AND CONCLUSIONS
We have so far demonstrated the wavelengthindependence of E c , both experimentally and theoretically, (Figs. 2 and 3) and we have shown that a scattering-limited trajectory model reproduces this behavior (Fig. 1).Importantly, a main conclusion arising from our results is that if E c is limited by the electron mean free path λ MFP in the liquid, then E c should scale with the density of the liquid.Figure 4 demonstrates that this is the case, both experimentally (A) and theoretically (C).
Very recently, some of us have reported the first measurements of the temperature of liquid flat-jets [50].The temperatures, measured by Raman thermometry under conditions identical to those of the present HHG experiments, range from ∼300 K at the top, to ∼255 K at the bottom of the first sheet, translating to a density variation of close to 5% (Fig. 4(B)) in the case of ethanol [56].Over this range of conditions, E c is found to decrease by ∼2 eV (panel (A)).The same trend is also observed in the calculations performed on liquid water (panel (B)).
Having measured E c over a range of densities, we can now verify how the maximal energy E k,max , gained by the electron from the driving laser field, scales with the density.Experimentally, we use E k,max = E c − E g , where E g is the band gap of the liquid.We find that E k,max scales linearly with the inverse density (blue symbols in Fig. 4(D)).This type of scaling precisely corresponds to the prediction of our simple trajectory-limited model (Fig. 1), because λ MFP =1/(nσ) ∝ 1/(ρσ).This conclusion is further supported by converting the measured E c to the corresponding maximal excursion length l max .As we show in the SM, Section S4.D, we find that E k,max = 3.73/4 * e * E * l max .A direct consequence of this relation is that it allows us to retrieve l max from the experimental spectra, provided that they were recorded under conditions where the wavelength-independence of the cut-off is observed, which is the case here (see Fig. 2D).The orange symbols in Fig. 4D show that l max also scales linearly with the inverse density.We therefore conclude that all experimental and theoretical evidence available at present agrees on the fact that E c is proportional to the maximal excursion length of the laser-driven electrons in the liquid phase.This suggests that it should be possible to accurately determine effective electron mean-free paths (MFPs) from liquid-phase HHG spectra.Electron MFPs play a very important role in describing electrondriven processes in the liquid phase, but they are notoriously difficult to measure and calculate, at least at low energies.The interest in developing new methods for accessing these quantities is therefore considerable and relevant for many physical processes.Here, we do not attempt to determine the MFPs with high precision because this would require a more sophisticated scattering model, including a large number of different scattering channels (see Ref. [58] and references therein).Instead, we aim at retrieving an effective MFP, which is best thought of as containing all types of scattering processes.Since the elastic scattering cross sections are by far dominant at the very low kinetic energies (∼10 eV) of interest here [53,[58][59][60], we compare our results to the elastic MFPs in Fig. 5.In this comparison, we use λ MFP = l max , taking into account that the electron travels up to the maximal excursion length before getting scattered.In the SM (Fig. S10), we show that this simple approximation is physically meaningful because replacing the "sharp" truncation of the trajectories (at the travel distance l max ) with an exponential distribution of path lengths (inherent to the definition of λ MFP ) leaves E c unchanged.Figure 5 compares the λ MFP values obtained from the HHG spectra (symbols) with the available literature values.In the case of liquid water, we are comparing to the most recent MFPs (blue dashed line), which were determined from a Monte-Carlo simulation of experimental liquid-microjet data using the most accurate ab-initio differential scattering cross sections available to date [53].In the case of the alcohols, liquid-phase MFPs have to our knowledge not been reported in the literature so far.We are therefore comparing our results to MFPs determined from the corresponding experimental gas-phase elastic scattering cross sections and the known number densities of the alcohols.The agreement is very good in all cases, confirming the possibility to retrieve electron MFPs from liquid-phase HHS.To summarize, we explored here the microscopic mechanisms responsible for liquid HHG with a combination of experimental and theoretical methods.Our measurements in water and ethanol show that contrary to crystals and gases, the cut-off energy in liquid HHG is mostly independent of the laser wavelength.Microscopic quantum mechanical calculations based on both supercells and clusters agree with this result, and show that it extends to other liquids and laser conditions.We proposed an extended semi-classical model for the electron dynamics in the liquid to explain this result.The model incorporates effects of ultrafast scattering of electrons off neighboring molecules, which are shown to reduce the HHG cut-off compared to the gas phase case through the spatial limitation of electron trajectories.The model reproduces well the wavelength independence of the HHG cut-off and highlights the importance of the electron mean-free path in liquids, indicating that this quantity is imprinted onto the high-harmonic spectra and can be retrieved.We also expect that our results would be highly relevant for HHG from amorphous solids [39,61].Our work paves the way to a deeper understanding of the strong-field dynamics in disordered condensed phases, and to resolving attosecond dynamics in liquids.

METHODS
The experimental set-up consists a 1 kHz Ti:Sapphire laser delivering ∼ 30 fs pulses at 800 nm.Driving wavelengths of 1500 nm and 1800 nm are obtained by optical parametric amplification (OPA) of the 800 nm pulses, respectively.The driver beams are focused with a spherical mirror on a sub-micron-thin liquid flat-jet target, further described in Refs.[42,43].The beam intensities are calculated from the harmonic cut-off energy of the gas-phase measurements, using the semi-classical TSM [19].The emerging high harmonics are analysed with a custombuilt XUV spectrometer consisting of an aberration-free flat-field grating (SHIMADZU) and a multi-channel plate (MCP) coupled with a phosphor screen.The phosphorscreen image is recorded by a CCD camera.Each spectrum is typically integrated for 10 -20 ms and measured 20 times.These spectra are averaged prior to all subsequent analysis.For the ellipticity dependent studies, an elliptically polarized 800 nm light wave generated using a combination of a rotating half-wave plate(HWP) and a fixed quarter-wave plate (QWP).The rotation of the HWP axis from 0 • to 22.5 • with respect to the QWP axis, changes the polarization of the input light from linear to circular keeping the axes of the polarization ellipse fixed.At a number of different ellipticities, ranging from = 0 for linear polarization to = 1, for circular polarization, the harmonic spectrum was measured for different liquids.Further details of the experimental and theoretical methods are described in the SM, Section S1 and the extended data Figure 1.
ture of surface-active salt solutions: photoelectron spectroscopy and molecular dynamics simulations of aque-ous tetrabutylammonium iodide, The Journal of Physical Chemistry B 108, 14558 (2004).

FIG. Extended data 1. (A)
Schematic of the IR beam optical path and the high-harmonic radiation in the top part of the liquid jet (see also Ref. [42]).The bottom-left inset of panel (A) shows the overview with regions labeled a-c.The highharmonic emission from the gas before the jet (region a) is absorbed in the first few layers of the liquid (∼10 nm [62]).Any XUV radiation generated in front of the liquid jet passes unrefracted through the liquid medium.This is because the index of refraction of liquid ethanol (as well as other alcohols and water) in the relevant XUV range (7-44 eV) is very close to 1 [62].The high-harmonic radiation generated in the bulk liquid (region b) is not refracted significantly at the exit liquid-gas interface.As a result, the XUV beam is refracted negligibly in comparison to the IR beam (the refractive index nIR is 1.35 at 1800 nm in ethanol [63]), as it exits the liquid phase into region c.In the schematic θi and θ i denote the angle of incidence of the IR beam at the entry and exit surfaces of the liquid jet.θr and θ r denote the angle of refraction of the IR beam at the entry and exit surfaces of the liquid jet and the lines T and N denote the tangent and normal to each surface.The strong refraction of the IR beam at the exit liquid-gas interface gives rise to the spatial separation of the liquid-phase harmonics (generated from the bulk liquid (region b)) and the gas-phase harmonics (generated from after the liquid jet (region c)) of ethanol at 1800 nm, on the micro-channel plate (MCP) detector.The green dashed box on top of the liquid-phase harmonics indicates the region of interest on the detector selected for extracting the gas-phase harmonic spectrum (C).Similarly, the purple dashed box indicates the region of interest selected for extracting the liquid-phase harmonic spectrum (D).The green dashed box below the liquid-phase harmonics indicates the region where the transmitted gas XUV (generated in front of the liquid jet (a)) is expected to be incident on the detector.The absence of these harmonics is explained by their expected absorption in the liquid jet.(C) We also studied the harmonic emission from 50-mmol Tetrabutylammonium iodide (TBAI) solution at 45°jet orientation.TBAI was chosen here because of its high surface affinity caused by the hydrophobic interactions of the TBA + cation and the large polarizability of the I − anion [64].This should further enhance the visibility of any even surface/interface harmonics arising from the breaking of the inversion symmetry of the p-polarizer driver by the liquid jet.We also compared this spectrum with that of pure water (D) for the same 45°jet orientation.However, even harmonics were not observed in any of these cases, indicating that the surface/interface contribution to the HHG response is negligible as compared to the bulk of the liquid.The corresponding HHG spectra are plotted for reference, along with the laser electric field (in solid white).The cutoff from the semi-classical model is indicated in dashed black.(B) Same as (A) but for an isolated molecule, where the response is orientation averaged to describe an un-oriented gas of water molecules.The cutoff in dashed black indicates the one obtained from SFA trajectories, and only SFA trajectories from the HOMO orbital are plotted (while some low yield higher energy emission still exists due to emission from the HOMO-1 and HOMO-2 orbitals).Both plots show the same kinetic energy scale and are calculated for 900-nm laser driving with an intensity of 5×10 13 W/cm 2 .The plot clearly indicates that in the liquid phase there is a dominant very-short trajectory contribution, whereas longer trajectories are strongly suppressed.The ab-initio cutoff for the liquid agrees very well with the prediction from the semi-classical model that incorporates the MFP as extracted from experiments.The time-frequency structure for the gas-phase similarly agrees with the expected structure, validating the approach.We note an additional emission contribution in the liquid phase at longer times that possibly corresponds to some other scattering process or a different emission channel that will be investigated in future work.The time-frequency analysis is obtained with a Gabor transform with a Gaussian window of T/3, and the isolated molecule calculations are obtained with the methodology described in ref. [57].
The experiments were performed with a 1 kHz regeneratively-amplified Ti:Sapphire laser (Coherent).From the laser output, 1.2 mJ beam was split and directed to a dedicated grating compressor delivering about 960 µJ, ∼30 fs pulses with a central wavelength of 800 nm.The linear and vertically (s) polarized driving beam was focused on a liquid flat-jet target with a pair of metallic mirrors (protected Ag) in a Z-shaped configuration, where the first mirror was flat and the second mirror was spherical with R = 800 mm.The whole set was mounted on a common base plate on top of a manual translational stage to allow fine tuning of the overlap between the laser focus and the flat jet.
For experiments with mid-IR wavelengths (1500 nm and 1800 nm), a commercial optical parametric amplifier (OPA, HE-TOPAS from Light Conversion) pumped with 6.5 mJ, 30 fs pulses at 800 nm was used.The signal or idler beam was separated from the residual wavelengths exiting the OPA using a pair of dielectric mirrors for the particular wavelength and focused with the same set of two metallic mirrors.

B. Flat Jet
The flat jet was created by two colliding cylindrical jets with about ∼54 µm inner diameter.The samples were pumped via an HPLC pump (JASCO) into the interaction region with a flow rate ranging from 3 ml/min to 5 ml/min.The established flat jet was about 1 µm thick [1].Due to evaporation, the liquid jet was surrounded by a region with high gas density, as compared to the backing pressure in the target chamber (∼5 • 10 −4 mbar), * These five authors contributed equally allowing a simultaneous probing of liquid and gas phases.All liquid samples have been used without further purification: pure Milli Q water with an electric resistivity of 18 MΩ and 99.8 % pure alcohols (ethanol, methanol, iso-propanol) from Sigma Aldrich.

S2. INTENSITY SCALING OF GAS-AND LIQUID-PHASE HHG SPECTRA
Figure S1 presents experimentally measured HHG spectra emitted from liquid water at 800 nm driving wavelength for different laser intensities in the range of 1.4 • 10 13 W/cm 2 to 5.8 • 10 13 W/cm 2 .At each intensity, we compare the gas-(Fig.S1A) and liquid-(Fig.S1B) phase spectra.For the gas phase there is a systematic linear increase of the cutoff with laser intensity (red dashed line), in accordance with the standard TSM.In contrast, the liquid HHG cut-off energy remains roughly constant (blue dashed line).This effectively demonstrates that in the given intensity range, the cut-off energy of the liquid spectra is independent (or weakly dependent) on the peak laser intensity.We have also verified that this is the case for liquid ethanol.
A. Scaling of the highest harmonic order vs.
scaling of the plateau's cutoff In relating the present results to the previous literature, it is important to distinguish the scaling of the highest emitted harmonic order (E max ) from the scaling of the plateau's cutoff (denoted E c , and discussed in the main text).Previous literature on condensed-phase HHG has mainly studied and discussed the scaling of E max .Most prominently, this is the case for Ghimire et al. for the solid state [2] and for Luu et al. for the liquid phase [3].In contrast, the overwhelming majority of the literature on gas-phase HHG has studied and discussed the scaling of E c , influenced by the definition provided by Lewenstein et al. [4].Here, we show that the scaling FIG.S1.For 800 nm driver, a comparisons between gas (A) and liquid (B) phase high-harmonic spectra of H 2 O is given for three intensities.As expected, the gas phase cut-off energy depends on intensity (red dotted line) but the liquid phase is intensity independent (blue dotted line).For 800 nm driver, a comparison between gas(A) and liquid (B) phase high-harmonic spectra of H2O is given for a range of intensities.
of these two quantities is actually different in the liquid phase and thereby show that the present conclusion regarding the scalings of E c are consistent, or a least in no contradiction, with previous literature, in particular with Luu et al. [3].
Figure S1C and D show the data in a form that highlights the scalings of both E max and E c .Panel C shows that, in the gas phase, both E max and E c scale linearly with the intensity.Panel D, in contrast, shows that the scalings of the two quantities are noticeably different.A non-linear least-squares fitting of E max to the functional form E max ∝ I y returns y = 0.53, consistent with Luu et al. [3].The independence of E c on I is visible, both in Fig. S1B and in Fig. S1D.This validates the conclusions reached in the main text and clarifies their relation with previous work.
We also performed ab-initio calculations to test the in-tensity dependence of the cutoff with both the supercell (Fig. S3) and the cluster (Fig. S5) approach for liquid HHG.Both types of calculations show that the cutoff is independent of the laser driving intensity, beyond a certain threshold.In each case the cutoff increases with the laser intensity until reaching a saturation point, where the increase stops.For the supercell approach, this saturation is found at ∼ 0.25•10 14 W/cm 2 , while for the cluster approach it is observed at ∼ 0.5 • 10 14 W/cm 2 .Variations between the two methods are a result of slightly different electronic structures (band gaps, in particular) obtained with the two approaches.Notably, the ranges of laser intensities for which the cutoff is intensityindependent correspond well with the experimental measurements.
Lastly, we also point out that the observed weak dependence of the HHG cutoff with respect to the laser intensity is also captured by our proposed extended semiclassical model.Figure S6 presents the calculated HHG cutoff vs. the peak laser intensity with the extended semi-classical model, assuming 800 nm driving, and intermolecular distances similar to those in liquid water.Initially, the HHG cutoff increases quadratically with the laser intensity just as in the standard TSM (because the trajectories are very short and do not extend beyond l max ).However, this dependence is reduced to a weak linear scaling in the rage of intensities > 0.5•10 14 W/cm 2 .In fact, over the intensity range of 0.5 • 10 14 W/cm 2 to 10 14 W/cm 2 , the HHG cutoff only increases by ∼ 2 eV.This result substantially differs from the gas phase, where the standard TSM predicts that the cutoff should increase by ∼ 9 eV in that region.Notably, a change of ∼ 2 eV in the cutoff energy at 800 nm driving would only move the HHG cutoff by approximately one harmonic order, which might be difficult to detect experimentally.Thus, we overall conclude that in our examined conditions, the liquid-HHG cutoff is weakly dependent on the laser driving intensity, an effect that is described remarkably well by our suggested semi-classical picture that includes scattering.We note that potential improvements to our extended semi-classical model (e.g.relaxing some of the utilized approximations) might also improve its correspondence to the measured and ab-initio calculated results.In order to accurately determine the cut-off from the experimental spectra we use the following method.The signal at each harmonic wavelength from the experimental spectra is acquired.It is known that the harmonic signal beyond the plateau region shows as exponential decay, which results in a linear decline in the logarithmic scale.Fig. S2 depicts this clearly, where the plateau region has a constant harmonic yield followed by a linear decay in the harmonic signal.The green line indicates the average signal of the plateau harmonics and the orange line indicates a linear fit to the intensities in the cut-off region.The intersection of the green and orange lines denotes the cut-off wavelength, which is converted to E c , the cut-off energy of the plateau.

S4. THEORETICAL DETAILS
A. Ab initio supercell calculations

Molecular dynamics simulations
In order to study theoretically the interaction of strong-field pulse with liquid water, the first step is simulating the structural and dynamical properties of liquid water.In this regard, we use Car-Parrinello molecular dynamics (CPMD) simulations in the canonical ensemble at the temperature of 300 K.The simulated system, as displayed in the main text, is a periodic cubic supercell, including 64 H 2 O molecules in the experimental liquid water density of 1 gr/cm 2 .Our computed radial distribution function between oxygen atoms (shown in Ref. [5]), which is averaged over a time length of approximately 20 ps, is in a good consistent with the experimental neutron scattering and x-ray diffraction measurements; it confirms the accuracy of our liquid water molecular dynamics simulation.Our CPMD simulations are performed using Quantum Espresso package [6].We use norm-conserving GGA-revPBE [7] pseudopotentials for oxygen and hydrogen atoms and Grimme-D2 dispersion correction to describe the non-local correlation effects.For more details please see Ref. [5].

Ultrafast dynamics simulations
In the next step, to study the interaction of ultrashort intense pulse with liquid target, we use real-time timedependent density functional theory implemented in Octopus package [8][9][10].Note that the ultrashort laser pulse stimulates the system in the time-scale of femtosecond while our CPMD simulation reproduce the macroscopic characteristics of liquid water in the time-scale of picosecond.So in order to obtain the converged response, we need to consider a sufficient number of different configurations of the liquid supercell system, denoted as 'subsystems'.The time interval between consecutive subsystems are in the order of picosecond, and the number of required subsystems for convergence depend on the target system as well as pulse characteristics (wavelength FIG.S3.High-harmonic spectra obtained from calculations with the supercell approach using a driving wavelength of 1500 nm. and intensity).In the supercell simulations, we consider the dynamics of nuclear based on Ehrenfest-TDDFT approach [11] during the interaction of ultrashort intense pulses with the liquid system.However, the effect of ionic current is negligible and the calculated HHG arises totally from electron dynamics.The exchange-correlation term in TDDFT calculations is described by the same GGA-revPBE pseudopotentials as those used in CPMD calculations.But, vdW interactions are not included, as they lead to no difference in the spectra.We consider a grid spacing of 0.3 bohr thorough TDDFT calculations, and a dense k-point grid of 5 × 5 × 5. Also, we use aetrs algorithm to approximate the evolution operator and the time step of 0.2 a.u.within our TDDFT calculations.

Pulse characteristics
In order to evaluate the impact of pulse wavelength and its intensity on HHG spectra, we consider the pulse characteristics corresponding to our experimental conditions.To see the impact of driving wavelength on HHG cutoff energy, we consider two different wavelengths of 800 and 1500 nm under the constant intensity of about 20 TW/cm 2 (shown in Fig. 3B of the main text).Additionally, to check the intensity scaling of HHG cutoff, the pulse intensity enhances from 25 to 100 TW/cm 2 , the pulse wavelength is constant and equal to 1500 nm (shown in Fig. S3).In all cases, the pulse duration at full width at half-maximum (FWHM) is equal to 18 fs with a sine-squared envelope shape for the vector potential.The numerical approach based on clusters that is utilized in the main text is thoroughly described in ref. [12].We provide here some additional technical details for the calculations presented in the main text.All calculations were performed using the real-space grid-based code, OC-TOPUS [8][9][10].Cluster geometries for liquid H 2 O were taken from ref. [13] using 43-molecule large clusters for the cut-off scaling calculations, and larger 54-molecule clusters for the calculations with varying liquid density that are more sensitive.Cluster geometries for NH 3 and CH 4 were taken from refs.[14] and [15], respectively, with 48-molecule and 40-molecule sizes, respectively.For each cluster, the ground state was obtained by DFT within the PBE [16] approximation for the XC functional with an added van-der-Waals correction term [17], with a realspace cartesian grid of spacing of 0.4 Bohr and spherical grid boundaries with a radius that extended 15 Bohr beyond the farthest atom from the cluster center.In the next step, the HHG response was calculated by propagating the KS states following the prescription in ref. [12] (by suppressing contributions from the surface states, employing orientation averaging, and utilizing the independent particle approximation (the XC potential is frozen to its ground-state form)).The laser pulse was taken to have a trapezoidal envelope in all calculations with a two-cycle long turn-on and turn-off, and a fourcycle-long flat top.The dipole response was filtered with a super-Gaussian prior to obtaining the HHG spectra.
Figure S4 show high-harmonic spectra obtained with this approach for NH 3 clusters (top) and CH 4 clusters (bottom).Both panels show a cut-off energy E c that is clearly independent of the wavelength.Ammonia was chosen as an example of a different polar molecule that also forms hydrogen bonds (like water and alochols, studied in the main text) and methane was chosen as an example of an apoloar molecule.The results of Fig. S4 thus show that the wavelength independence of E c is not specific to the liquids studied in the main text, but appears to be a general property of liquid-phase HHG.
FIG. S5.High-harmonic spectra obtained from calculations on water clusters using a driving wavelength of 900 nm.
To complement the data shown in the main text (Figs.3C, 3D, 4C), we show here in Fig. S5 the scaling of the high-harmonic spectra with the peak intensity of the driving field with a central wavelength of 900 nm.The figure nicely shows the saturation of the cut-off beyond a threshold intensity.

C. Extended semi-classical calculations
Based on our suggested extended semi-classical model (main text, Section II.A), we performed semi-analytical calculations for the HHG cutoff vs. the driving laser wavelength and intensity.We also performed calculations with varying l max parameters.In all cases, the Newton equations of motion for the electrons were solved analytically for a given time of ionization, and then represented numerically on a temporal grid to scan for recollisions.The solutions were numerically filtered for each re-colliding trajectory to test whether that particular trajectory extruded beyond l max or not, while the standard TSM was not filtered with this condition.In each case the cutoff was given by the most energetic recolliding trajectory with the added band gap of the liquid.The calculations in Fig. S6 confirm the linear scaling of the maximal kinetic energy with the peak electric-field amplitude of the driving pulse.The onset of saturation is observed at the highest field amplitudes for l max = 7 Å.

D. Semi-classical description of HHG and cutoff scaling including scattering
In this section, we aim at deriving an approximate fully analytical picture that describes the cutoff energy of HHG, as described by our extended semi-analytical semi-classical three-step model.We examine what is the analytical connection between the HHG cutoff and the maximum allowed excursion length of the electron trajectories in the extended model.Our starting point is the full semi-classical picture discussed in the main text.We consider the classical motion of an electron in an electric field E(t) = E cos(ωt), where ω is the frequency of the laser.The magnetic field, as well as the electric field due to the parent ion are neglected and the dipole approximation is employed.The electron is ionized at a time t i , and we consider the dynamics following this time.
where we imposed as initial condition that the trajectories start at x(t i ) = 0, with a vanishing velocity (v(t i ) = 0).The time at which the electron returns to its parent ion, referred here as the return time t r , is obtained by imposing that x(t r ) = 0.This gives a first relation, that defines trajectories returning to the parent ion: From the return time, we can obtain the kinetic energy for a trajectory starting at t i The solution of the Eq. ( 3) can be performed numerically, see Fig. S7, from which we obtain the well known gasphase result for the maximum kinetic energy of ≈ 3.17U p , corresponding to ωt i ≈ 0.305.

FIG. S7
. Evolution of the return kinetic energy versus the ionization phase, obtained from Eq. ( 4), solving numerically Eq. ( 3) for each value of ti.
Let us now consider the case that not all trajectories are allowed, as formulated in our extended model for the liquid.We first need to determine for a given trajectory, determined uniquely by t i and Eq. ( 3), what is the maximum excursion.The maximum excursion is reached when ∂ t x(t) = v(t) = 0.This gives us the condition defining the time of maximal excursion t * : which has a simple analytical solution ωt * = π − ωt i .We therefore have an exact expression for the excursion distance l exc for any given ionization time ) with e = |q| is the elementary charge.Solving this equation numerically for the most energetic trajectory (ωt i = 0.306), we find that without limiting the trajectories, the cutoff energy corresponds to an excursion of ≈ 1.145 eE mω 2 , which verifies the analytic analysis.
We next consider that the maximal excursion is limited to l max .We can distinguish two cases.As long as 1.145 eE mω 2 is smaller than l max , the kinetic energy gained is given by the usual formula ∆E K ≈ 3.17U p .Above this threshold, the value of the kinetic energy gain will be modified, as the most energetic trajectory is forbidden, and the cutoff will be reduced.Below, we assume that we are in this regime.The task is therefore to find the most energetic trajectory that returns to the parent ion, and whose excursion is smaller or equal to l max .This gives us the condition The function l exc (t i ) is a function that decreases monotonously from 2 to 0 over the interval of interest of the ionization time [0 : π/2].Therefore a value t * i exists such that the condition ( 7) is satisfied, for all t i > t * i (and smaller than π/2).Because we are in a regime where not all the trajectories are allowed, this time corresponds to a "short trajectory", meaning that the ionization phase is larger than 0.306, which is the time corresponding to the most energetic trajectory without scattering.Said differently, we only select the short trajectories that have an excursion smaller than l max .Now, we also know that the return kinetic energy decreases monotonously for the short trajectories (see Fig. S7).Therefore, the critical time t * i corresponds to the highest kinetic energy in our scattering model.This ionization time that corresponds to the most energetic trajectory with an excursion smaller or equal to l max is given by the relation where we define the dimensionless quantity L red = mω 2 lmax eE .We therefore have two defining equations, Eq. ( 8) and Eq. ( 3), where the first one determines the ionization time as a function of the reduced maximum excursion L red , and the second one that determines the corresponding return time.
We have solved numerically these two equations: for each value of L red , we can find the numerical solution of Eq. ( 8) and then use this value to numerically solve Eq. ( 3).The result of the kinetic energy (in units of U p ) as a function of L red is shown in Fig. S8.For not too large values of L red , this is well approximated by a linear function.In this regime, which is the one relevant for our experiment, we find that From this result, we clearly see the role of introducing a mean-free path in our model.The cutoff becomes independent of the wavelength and is linear with the electric field.The comparison between the full numerical result and the linear fit is shown for different laser parameters in Fig. S9, showing that for our range of parameters, the approximate linear scaling law matches very well the exact numerical result of our model.

S5. INTENSITY DEPENDENCE OF LIQUID-HHG CUTOFF
In the main text we concluded that the liquid-HHG cutoff was independent of the driving wavelength.However, since the HHG measurements for different wavelengths slightly vary in the driving peak laser intensity (because the threshold intensity of plasma generation depends on the laser wavelength), this conclusion is only valid if the liquid-HHG cutoff also weakly depends on the driver peak intensity in these conditions.We here show that this is indeed the case with both measurements and calculations, which validate that the HHG spectra from liquids is independent of the driving wavelength.

S6. COMPARATIVE STUDY OF CUT-OFF ENERGY WITH DIFFERENT TRAJECTORY-TRUNCATION METHODS
To study the dependence of the harmonic energy on the implementation of the electron mean free paths, we compare the cut-off energies obtained from different truncation methods.Fig. S10 shows the variation of harmonic yield (i.e. the number of returning trajectories) as a function of the kinetic energy.The analysis is run for different models to find the harmonic yield dependence on the harmonic energy due to the energy dependent λ MFP .The simulations are performed with a 1500 nm light of 1.6V/ Å electric field amplitude.The purple solid line denotes the gas-phase classical SFA model where all trajectories upto the cut-off return with equal probability.The red circles indicate the trajectory model where returning trajectories with excursion lengths greater than 3.1 Å are rejected.The maximum kinetic energy is therefore reduced to ∼2 eV as opposed to 25 eV (SFA).However, in reality trajectories with trajectory lengths extending beyond the λ MFP are expected to return but with exponentially decaying probability.Indeed, the probability that a particle is scattered between the travelled distance l and l + dl is given by where λ MFP is the mean free path of the particle.From this definition, we can compute the probably that a particle is getting scattered after a travelled distance L, given by the travelled distance at the returned time.This is given by As L grows, the probability converges to one, as one would naturally expect.In this treatment, we have assumed that the mean free path is independent of the energy of the particle, and therefore the scattering probably only depends on the travelled distance.
It is however possible to go beyond this and to take into account the whole variation of the mean free path along the path of the electron, as we know analytically the trajectories.For this, we convert the integral over the travelled distance into an integral over the time, using the fact that the travelled distance at a point in time t is given by l(t) = t ti |v(t )|dt .
(12) The yellow line indicates further sophistication of the trajectory model as given by equation (11).This model considers exponential damping by a fixed λ M F P value of 3.1 Å. Finally the blue solid line indicates the returning probability for trajectories calculated considering energy-dependent MFP given by equation (12).We observe that even with the statistical Monte-Carlo-like energy-dependent MFP model (solid blue line), the harmonic yield (returning probability) varies slowly up to ∼ 2 eV falling to only 0.5 the value with respect to the hard-cut off and the SFA model.Beyond ∼ 2 eV the harmonic yield falls sharply indicating the onset of the cut-off region.This is in considerable agreement with the simplest hard-cut-off model.As a result, the harmonic yield below cut-off does not change considerably to extract the energy-dependent λ M F P (E) with reasonable accuracy taking into account the effect of experimental parameters.Beyond the cut-off energy there is a significant dependence of the harmonic yield on the harmonic energy resulting from λ M F P (E).However, the above cut-off harmonics also have contributions from back scattering and nearest-neighbor recombination as discussed previously.Thus isolating the effect of the energy-dependent λ M F P only for above-cut-off harmonics is beyond the scope of this manuscript.

S7. ANALYTICAL TREATMENT OF THE PROPAGATION EFFECT
Our goal is to describe how the laser-induced current density lead light emission at a given frequency propagate in a medium, cross the surface of the material and reach the detector.We first calculate the E and B fields that result from given source terms.Our approach is quite similar to the approach of Ref. 18, that we generalize here to not only treat the polarization, but also the free charge current and densities, denoted respectively ρ(r, t) and J.In the context of high-harmonic generation in solids, it is mandatory to keep these terms, as the freecharge current (J) describes the intraband contribution while the polarization P, describes the interband contribution.
We begin with the macroscopic Maxwell equations for an homogeneous system where The total macroscopic polarization, P tot , and magnetization, M tot , have in frequency space the following expression where the macroscopic susceptibilities χ e and χ b describe the linear response of the medium, and are directly related to the macroscopic dielectric and the permeability tensors of the system Note that as we are interested here in systems which are homogeneous at equilibrium, these quantities do not depend on space.The remaining polarization P nl (r, ω) (magnetization M nl (r, ω)) term corresponds to the nonlinear contributions to the total polarization (magnetization).
Rewriting the Maxwell equations only in terms of E and B fields, we obtain where we defined the total macroscopic nonlinear current and charge densities as These are the source terms induced by the strong incoming laser, which in term induce the harmonic fields.There are many situations where the spatial variations of the current cannot be neglected.Due to the multiphoton absorption, the pump laser is depleted, which implies a spatial variation of the laser intensity, and hence of the generated current along the propagation direction.As the laser-matter interaction is highly nonlinear in the case of high-order harmonic generation, many effects have to be taken into account when describing the propagation of the pump laser in the matter, such as self-focusing effects [19], which is caused by the third-order non-linearity of the medium.Moreover, due to surface roughness, nano-structuration of the surface, or possible metalization of the first few atomic layers at the incident surface, the incoming laser can excite surface currents and surface-plasmon-polaritons when the field crosses the interface [20].It is clear that all these phenomena will induce a spatial variation of the pump laser will propagating in the matter, even at normal incidence, and will thus induce a spatial dependence of the generated nonlinear current.However, while our treatment is general and does no require any approximation of the spatial dependence of the current, we will in the following assume that we have a film which is thin compared to the Rayleigh length of the driving pulse, and neglect the spatial dependence and propagation effects of the incoming laser field.We now proceed, following the approach of Ref. 18 and we first give the solution of the homogeneous problem defined by Eqs. 16, without considering the presence of the surface.In the view of treating a surface, we construct plane-wave solutions of Eqs.16 which have a real momentum in the (xy) plane of the surface.[18] The z-axis is chosen to be perpendicular to the surface plane.
There is therefore two possible wave-vectors associated where κ is the in-plane unit vector of the form κ = cos(φ)x+sin(φ)ŷ, φ being the azimuthal angle, and q ⊥ (ω) is the out-of-plane part of the wave-vector.Here q(ω) 2 q ⊥ being chosen to have Im[q ⊥ ] ≥ 0 and Re[q ⊥ ] ≥ 0 if Im[q ⊥ ] = 0.This convention allows us to be sure that the wave associated with q + is upward propagating and that the one associated with q − is downwardpropagating.[18,21] We thus refer to the solutions of the homogeneous problem associated with q + (respec.q − ) as the upward-propagating (respec.downwardpropagating) solutions.The treatment of the problem represented by equation Eqs.16 if greatly simplified if we introduce the s and p polarizations, defined by the following unit vectors This is illustrated in Fig. S11(b).For conciseness, we have omitted the frequency dependence here.
One can easily check that we have the following relations ŝ × q± (ω) =p ± (ω) , q± (ω) × p± (ω) =ŝ , p± (ω) × ŝ =q ± (ω) , where q± (ω) = q ± (ω)/q(ω).Thus (ŝ; q+ ; p+ ) and (ŝ; q− ; p− ) are direct basis.Therefore, the wave associated with q+ (respectively q− ) only propagates along ŝ and p+ (respectively p− ).Therefore, the upward-propagating solutions of the homogeneous Maxwell equations is of the form The downward-propagating solution as the same form and is not reported here for conciseness.Knowing the solution of the homogeneous problem, we can now search for the solution of the full problem, described by Eqs.16.The search of the solution of Eqs.16 proceed in two steps.We first find the solution associated with source terms of the form [18] with r || = (x, y).This corresponds to a sheet of current and density located in the plane z = z 0 and having a in-plane wave-vector q || (ω).The physical solution of the Maxwell equations, for such a source term located at z = z 0 (which excludes exponentially diverging waves) still fulfill the homogeneous Maxwell equations for z = z 0 and therefore reads [18] E(r; ω) = E + (r; ω)Θ(z − z 0 )e −iq ⊥ (ω)z0 +E − (r; ω)Θ(z 0 − z)e iq ⊥ (ω)z0 and a similar expression for the magnetic field, with here Θ(z) being the Heaviside function.
The determination of the coefficients E s ± (ω) and E p ± (ω) of Eq. ( 18) proceed exactly as explained in details in Ref. 18, and is therefore not reported here.After some algebra, we obtain that These coefficients determine completely the electric and magnetic fields induced by the source terms of Eqs.(19).
The total electric field induced by one sheet of polarization located at z 0 , and evaluated at the position r, with a given in-plane wave-vector q || , denoted E q || ,z0 , is thus given by A similar expression holds for the B field.The first term of this expression corresponds to the upward propagating field radiating from the considered source term sheet.
The second term is the downward propagating field, and the last one is the local contribution that does not radiate [18].
Exploiting the linearity of Maxwell equations (Eqs.16), we can obtain the total electric field originating from any general source terms n tot (r, ω) and j tot (r, ω) by integrating over the corresponding in-plane wave-vector and z 0 , the latter ranging from z b to z t .The general fields E(r; ω) are obtained by summing the individual contributions from all the source terms sheets Similarly, one can also defined the total upward propagating electric field and the magnetic field.After some algebra and using Eq. ( 19), we obtain Similar expressions are obtained for the downward propagating electric field and the magnetic upward/downward propagating fields.In the following, we assume that we can neglect the in-plane position dependence of the electronic current.Thus A. Measured intensity for thick films From the previous results, we can now obtain the expression of the radiated harmonic field in the entire sys-tem.However, we are only interested in the transmitted harmonic field, which corresponds to the to the electric field in vacuum, measured at some value of z such that z > z t .As we are interested by a macroscopic treatment, the surfaces are modeled by abrupt surfaces with effective Fresnel coefficients.These effective coefficients can be always chosen such that they include the microscopic details of the surfaces.At each interface, a Maxwell saltus occur and only the in-plane component of the electric field's wave-vector is conserved.Therefore, choosing the polarization and wave-vector of the incoming field determines uniquely the fields in vacuum.We first neglect the multiple reflection, assuming a thick slab.The case of a thin-film with multiple reflections is discussed in the next section.Knowing the expression of the upward-propagating harmonic fields, we can now look at the expression of the transmitted harmonic light in vacuum.The harmonic electric field emitted by the crystal is considered to be weak enough to not induce new nonlinear effects, and to have its propagation well described by linear optics laws.At the interface with the vacuum, the field is therefore modified according to the Fresnel coefficients.The harmonic electric field that reaches the detector is thus where E + (z − t ; nω l ) is the upward propagating field computing just below the surface.
If we assume no pump depletion and neglect nonlinear effects in the pump propagation, we can drop position dependence in source term, the electronic current, and replace it by the power n of the phase of the incoming laser, which leads to where the last term is the phase-matching term and L = z t − z b .This implies here that we are in the perturbative regime, i.e., that the n-th harmonic scales as the n-th power of the electric field.The modification of the refractive index at the incoming frequency, due to the Kerr effect for instance, can still be taken into account in q ⊥ (ω l ).Note that experimental values points toward a dependence in the electric field of the n-th harmonic much lower that the n-th power of the electric field, and one might replace nq ⊥ (ω l ) by n eff q ⊥ (ω l ) where n eff is an effective number for all harmonic order.From this expression, the measured intensity of the harmonic field is given by the usual relation I(nω l ) = c 0 2 |E out (nω l )| 2 .One special case of interest is normal incidence.As we assumed a cubic, non-magnetic crystal, a s(p)-polarized driving field lead to a s(p)-polarized current, and thus I pp (nω l ) = I ss (nω l ) = I(nω l ), with This expression is the generalization of the typical phase-matching expression that one would obtain from perturbative non-linear optic to the n-th order harmonic generation from a slab of matter.It contains the propagation of the incoming light pulse from z b to z and then the subsequent propagation of the emitted light from z to z t , for all possible values of z.

B. Measured intensity for thin films
In the case of a thin-film, it is not possible to neglect the multiple reflections of the harmonics inside the ma-terial.In this case, Eq. 26 has to be modified, to account for internal reflections.The resulting total outgoing field is composed of the transmitted upward propagating field, plus the multiple reflection of the upward and downward propagating waves.
where L = z t − z b is the thickness of the thin-film and r± = [ŝr s mv (nω l )ŝ + pm ± (nω l )r p mv (nω l )p m ∓ (nω l )].Following the same approach as in previous section, we now neglect the position dependence of the current and assume a perturbative dependence on the propagation phase of the incoming laser pulse.In the case of normal incidence, assuming again a cubic non-magnetic material, the s and p polarization give the same result and the intensity of the radiated field is then given by Equation 30 formula is the equivalent of Eq. 28 for the case of multiple reflections.We found numerically that in our experimental conditions, the multiple reflections do not play any important role.

C. Results for the propagation effects
We now employed the above model to simulate the effect of light propagation in water jet, and in particular to simulate the effect of re-absorption of the harmonic light by the jet itself.We considered a jet thickness of 1 mum.We assumed that the laser generate an harmonic spectrum as in the gas phase at each possible position along the slab, and we employed the experimental di-electric function of water to model light-propagation effects [22].Due to the reabsorption, we found that the total measured signal, that reaches the detector, builds up within the 200nm close to the exist surface.However, our results, shown in Fig. S12 shows that the propagation effects, including phase-matching effects and light re-absorption, cannot explain the change in the energy cutoff.We also found (not shown) that the multiple reflections of either the incoming light or of the emitted Effect of propagation for a 1µm thick slab of liquid water.The emission is assumed to be given by the gas phase for the same laser intensity as for the liquid phase measurement.The effect of harmonic light propagations is computed here for (a) 800 nm driving field, and (b) 1500 nm driving field.In both cases, the change in the energy cutoff cannot be explained by the laser reabsorption.harmonic light are not changing significantly the effect of the propagaition.

S8. EXPERIMENTAL RESULTS FOR DENSITY AND PROPAGATION EFFECTS
To isolate the effects of propagation and density liquid generated high harmonic spectra we experimentally measured jet thickness dependence of the harmonic cutoff and harmonic yield.If the cut-off energy is dependent on the propagation length through the liquid, a change in the jet thickness would show a change in the cut-off energy.The thickness for a flat-jet varies in the vertical direction.However, we have seen for Ethanol liquid (Figure 4(A) of manuscript) that the temperature also changes along the vertical direction.To isolate the effect of only the thickness variation we specifically choose H 2 O and the vertical positions of(0.66-1.26mm from top of jet).The thickness has been measured using a few-cycle 800 nm XUV-XUV high harmonic interference and verified through white light interferometry [23].Figure S13(A) shows the fractional density change in percentage as a function of the jet position and the corresponding jet thickness.The density change has been calculated from the temperature variation as a function of the jet position, obtained from the calculations of [24] for H 2 O flat-jet formed using 40 um nozzle diameter(used for the thickness measurements here).[24].The corresponding water density has been interpolated from Table 17 of [25].It is observed that in this vertical scan range, while the thickness changes from 1000 nm to 500 nm the fractional change in density with respect to the 1000 nm data point is only about 0.04%, thus for H 2 O and the given jet positions we can clearly say, that whatever cut-off energy shift,in any, we observe should solely be because of the jet thickness variation.Figure S13 (B) shows the harmonic spectrum of H2O measured for different jet thicknesses using 800 nm driving wavelength.Each harmonic spectra has been normalized to the 5th harmonic (energy 7.76 eV) of 800 nm.The Figure S13 (B) clearly indicates that the harmonic spectra are almost identical, irrespective of the jet thickness.To observe this carefully we plot Figure S13 (C) that shows the ratio of each harmonic signal (Normalized yield summed over an energy bin of 0.5 eV centered around each harmonic peak).The top x axis of Figure S13 (C) indicates the fraction change in density in percentage with respect to the data point at 1000 nm thickness (vertical position 0.66 mm).It is observed that the harmonic ratio varies negligibly as a function of the jet thickness.This observation is in contrast to what we see for Ethanol: Figure 4(A) (Figure S14): where the fractional density is significantly large and we see a clear variation of the harmonic ratio as a function of the fractional density change.This further supports the microscopic origin of the cut-off energy shift, that is now clearly isolated to be a density effect as opposed to an absorption effect of the liquid medium.
We compare these results to the harmonic yields observed for different spectrums presented in Fig. 4(A) of the manuscript, where we have normalized each spectrum to the 11 th harmonic of 1800 nm. Figure S14 (A) plots the normalized harmonic yields as a function of harmonic energy for the four different spectrums of Fig. 4(A) where the first four energies correspond to the harmonics below cut-off (H11 H13 H15 H17) and the two highest energies(H19 and H21) correspond to harmonics above cut-off.The harmonic yield has been computed by a summation of the harmonic signal within an energy bin of 0.5 eV centered on each harmonic peak.The measured signal is then normalized to the signal of harmonic 11. Figure S14 (B) shows the harmonic ratios obtained for different density variations from the cluster calculations done for a laser wavelength of 900 nm.Comparing Figure S14 (B) with the Figure S14 (A) we see that the harmonics below the cut-off behave similar to what we expect from trajectory simulation, where reducing the fixed λ M F P values (equivalent to increasing density) systematically suppresses the harmonic yields, which essentially emphasizes the validity of the fixed MFP assumption.However for harmonics beyond cut-off (H19 and FIG.S13.Effect of jet thickness on the cut-off energy of liquid harmonics of H2O measured at 800 nm wavelength.(A) Shows the fractional density change in percentage as a function of the jet position and the corresponding jet thickness obtained from [23], [24] and [25].(B)The harmonic spectrum of H2O measured for different jet thicknesses using 800 nm driving wavelength.Each harmonic spectra has been normalized to the 5th harmonic (energy 7.76 eV) of 800 nm.(C) Ratio of each harmonic signal (Normalized yield summed over an energy bin of 0.5 eV centered around each harmonic peak).The top x axis of Figure S13 (C) indicates the fraction change in density in percentage with respect to the data point at 1000 nm thickness (vertical position 0.66 mm).H21) the trend seems to change increasing in yield for higher densities.This might involve contribution from the energy dependent MFP but harmonics beyond cutoff also have contributions from backscattering and nearest neighbor recombination as discussed previously.Thus isolating the contribution of just one effect is beyond the scope of this manuscript.

S9. COMPARISON OF CUT-OFF ENERGY OBTAINED FROM THE SEMI-CLASSICAL SCATTERING MODEL FROM NEAREST NEIGHBOR INTERACTION AND ONSITE INTERACTION
From the perspective of a semi-classical trajectory based model, the energy cutoff obtained from the re- (A) Normalized harmonic yields as a function of harmonic energy for the four different spectrum of Fig. 4(A) where the first three energy points correspond to the harmonics below cut-off (H11 H13 H15 and H17).(B) The harmonic ratios obtained for different density variations from the cluster calculations done for a laser wavelength of 900 nm and intensity of 4× 10 13 W/cm 2 .combination to the neighboring site is significantly larger than for a wave packet recombining to its parent ion, because of a higher kinetic energy at the collision time.Thus, if recombination channels on neighboring molecules were playing a significant role the cutoff energy would have been almost double than that of the onsite recombination, which contradicts our observations.The figure below shows the comparison of the maximum kinetic energy as a function of laser wavelength (Figure S15 A) and as a function of the laser field (Figure S15 B) for the on-site recombination and nearest neighbor recombination as calculated from the semi-classical scattering model.In addition to the nearest neighbor interaction having higher cut-off than the onsite interaction, we also observe that for wavelengths below 900 nm the nearest neighbor interaction should ideally give a cut-off higher than the SFA (for the gas generated harmonics).But from the data presented in Figure S1 of the supplement, this is definitely not the case.Moreover we see that for the nearest neighbor interaction the maximum kinetic energy has a stronger dependence on the electric field than the onsite interaction, which is in contrast to our observations in Figure S1(D).where we imposed as initial condition that the trajectories start at x(t i ) = y(t i ) = 0, with a non-vanishing velocity along the minor-axis direction (y), that guaranties that the electron returns to its parent ion (y(t r ) = 0) at the return time.As in linearly-polarized case, the return time is given by the condition x(t r ) = 0.This leads to v 0y = qE .From these equation of motions, we considered only returning trajectories which travel a distance d = tr ti |v(t)|dt smaller than an effective mean-free-path distance, and computed the corresponding gain in kinetic energy at the return time.This directly allows us to compute the change in cutoff energy for our model with respect to the driver ellipticity, as shown in the main text Figure ED1.We found that the energy cutoff decreases with increasing ellipticity.We note that if we instead would start from a zero velocity, we also obtain a decrease of he energy cutoff with the increasing ellipticity for small ellipticities, and an increase for large ellipticities, because of the unphysical starting point.If we consider instead of the recombination to the parent ion the recollision with the nearest neighboring molecule, we obtain that the increasing ellipticity leads to an increase of the energy cutoff, which is the opposite trend as observed experimentally.

FIG. 1 .FIG. 2 .
FIG.1.Effect of electron scattering on HHG spectra in liquids (A) Schematic illustration of the extended trajectorybased semiclassical model.An electron (green) is ionized by the laser field, accelerated and either recombines directly with its parent ion (solid green arrows), or scatters off another molecule (dashed green arrows).(B) The returning electron trajectories from the standard TSM within an 800-nm driving laser field of 1.9 V/ Å.The dots on the electric field represent the ionization times of the electrons (the colour of dots corresponding to the colour of the respective trajectories).The grey dashed line denotes the limited excursion length (lmax) imposed by scattering.(C) Wavelength-scaling of Ec in the absence (standard TSM: black dashed line) and presence (lmax-limited: colored lines) of scattering for a laser intensity of 5×10 13 W/cm 2 .The observable manifestation of scattering is a decrease in the cut-off by ∆Ec.

FIG. 5 .
FIG. 5. Comparison of the determined electron meanfree paths with literature data.Mean-free paths for electron scattering in the liquid phase determined from the experimentally observed cut-off energies (Ec) as a function of laser wavelength.The dashed blue line indicates the elastic meanfree path for liquid water [53] at E k = 4.5 eV.The arrows indicate the value of the mean-free paths obtained from the integral elastic scattering cross sections of the corresponding alcohols in the gas phase [54, 55] at E k of ∼ 3-4.5 eV, using the number densities obtained from the densities of different liquids at 20 • C [56].
FIG. Extended data 2. (A)Schematic showing the flat-jet orientation with respect to the incident laser direction used for excluding the presence of surface/interface-generated harmonics.The p-polarized (perpendicular to the plane of the figure) 800-nm laser beam is made incident on the flat-jet.The flat-jet is rotated such that its target normal is oriented at an angle θ with respect to the incident laser direction.If harmonics are generated at the interface where the beam exits the flat-jet, the liquid medium would break the symmetry of the light field (and therefore the electron trajectories) for the negative and positive half cycle of the electric field for any θ > 0, which would lead to the generation of even harmonics.(B) The raw MCP image for θ = 0 o (top panel) and the comparison of the harmonic spectra (bottom panel) taken for θ = 0 o (blue line) and θ = 45 o (orange line), for 800 nm laser beam incident on an ethanol flat-jet.No even harmonics are detected for the tilted jet geometry, which suggests that the detected signals are dominated by the bulk liquid with negligible interface contributions.(C)We also studied the harmonic emission from 50-mmol Tetrabutylammonium iodide (TBAI) solution at 45°jet orientation.TBAI was chosen here because of its high surface affinity caused by the hydrophobic interactions of the TBA + cation and the large polarizability of the I − anion[64].This should further enhance the visibility of any even surface/interface harmonics arising from the breaking of the inversion symmetry of the p-polarizer driver by the liquid jet.We also compared this spectrum with that of pure water (D) for the same 45°jet orientation.However, even harmonics were not observed in any of these cases, indicating that the surface/interface contribution to the HHG response is negligible as compared to the bulk of the liquid.
FIG. Extended data 3. (A) Experimental setup for studying the ellipticity dependence of high-harmonic generation from liquids.The rotating half-wave plate and fixed quarter-wave plate geometry is used to keep the axes of the elliptically polarized driving field fixed as the light changes from linear to circular polarization.(B) Harmonic spectra obtained from liquid water at 800 nm for different driver ellipticities.The pink dashed line in each sub-figure represents the cut-off energy for linearly polarized light.The blue dashed line indicates the actual cut-off for each of these spectra.It is observed that as the light changes towards circular polarization ( =1), the cut-off shifts to lower energies.(C) Comparison of the experimentally determined change in cut-off and that calculated from the MFP-limited scattering model which considers only recombination with the parent molecule.Details of the theoretical scattering model are given in the supplementary section S10.(D) experimental data showing high-harmonic yields vs. ellipticity.(E) Data from cluster calculations showing the high-harmonic yields vs.ellipticity.These calculations were performed using a 900-nm driving field with a peak intensity of 4×10 13 W/cm 2 .
FIG. S2.Plot of harmonic signal as a function of harmonic wavelength.The cut-off is determined by the intersection of the green and orange lines.

B
FIG. S4.High-harmonic spectra obtained from calculations on ammonia (NH3) clusters (top) and methane (CH4) clusters (bottom) over a broad range of driving wavelengths with a peak intensity of 50 TW/cm 2 .
FIG. S6.Scaling of the maximal electron kinetic energy with the peak electric field of the driving pulse.
FIG. S8.Evolution of the return kinetic energy versus L red .The linear and cubic fits are also shown.
FIG. S9.Evolution of the return kinetic energy versus electric field (left panel) and wavelength(right panel).The parabola correspond to the result without scattering, whereas the solid lines correspond to the result including scattering, for values of lmax ∈ {2, 4, 6, 8, 10, 12, 14} Bohr.The dashed lines correspond to the same results, but obtained using the simpler formula ∆EK = 3.73UpL red .
FIG. S10.Shows a comparison of the returning probability of trajectories with different kinetic energies for different variations of the trajectory based model.The purple solid line denotes the gas-phase classical SFA model where all trajectories upto the cut-off return with equal probability.The red circles indicate the trajectory model where returning trajectories with excursion length greater than 3.1 Å are rejected.The maximum kinetic energy is therefore reduced to 2 eV as opposed to 25 eV(SFA).The yellow line indicates further sophistication of the trajectory model where each returning trajectory is modulated by a scattering amplitude of exp(-lexc/3.1 Å).This model considers exponential damping by a fixed λMF P value of 3.1 Å. Finally the blue solid line indicates the returning probability for trajectories calculated considering energy dependent MFP.In this model at each time step of each trajectory, based on their instantaneous kinetic energy, a scattering amplitude of exp(-lexc/λMF P (E)) is applied to the returning probability.
FIG. S11.(a)The experimental geometry considered in this work.The incoming laser field E laser is at normal incidence.(b) Coordinate systems used to describe propagation of harmonics.The unit vectors and the link between the different coordinate systems are given in the text.The sketch is presented here for an azimuthal angle φ = 0. Adapted from Ref. 21 FIG. S12.Effect of propagation for a 1µm thick slab of liquid water.The emission is assumed to be given by the gas phase for the same laser intensity as for the liquid phase measurement.The effect of harmonic light propagations is computed here for (a) 800 nm driving field, and (b) 1500 nm driving field.In both cases, the change in the energy cutoff cannot be explained by the laser reabsorption.

FIG
FIG. S14.(A) Normalized harmonic yields as a function of harmonic energy for the four different spectrum of Fig.4(A) where the first three energy points correspond to the harmonics below cut-off (H11 H13 H15 and H17).(B) The harmonic ratios obtained for different density variations from the cluster calculations done for a laser wavelength of 900 nm and intensity of 4× 10 13 W/cm 2 .
FIG. S15.Comparison of maximum kinetic energies for onsite recombination and nearest-neighbor recombination: (A)As a function of wavelengths, calculated at an intensity of 5×10 13 W/cm 2 .(B) As a function of driving field for 1000 nm driving wavelength.