Extracting sub-cycle electronic and nuclear dynamics from high harmonic spectra

We present a new methodology for measuring few-femtosecond electronic and nuclear dynamics in both atoms and polyatomic molecules using multidimensional high harmonic generation (HHG) spectroscopy measurements, in which the spectra are recorded as a function of the laser intensity to form a two-dimensional data set. The method is applied to xenon atoms and to benzene molecules, the latter exhibiting significant fast nuclear dynamics following ionization. We uncover the signature of the sub-cycle evolution of the returning electron flux in strong-field ionized xenon atoms, implicit in the strong field approximation but not previously observed directly. We furthermore extract the nuclear autocorrelation function in strong field ionized benzene cations, which is determined to have a decay of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _0 = 4 \pm 1$$\end{document}τ0=4±1 fs, in good agreement with the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \tau _0 = 3.5$$\end{document}τ0=3.5 fs obtained from direct dynamics variational multi-configuration Gaussian calculations. Our method requires minimal assumptions about the system, and is applicable even to un-aligned polyatomic molecules.


S1 Comparison of harmonic spectrum from normal and deuterated benzene
We observed no significant differences between normalized harmonic spectra from benzene and deuterated benzene. Figure S1 compares typical spectra from the two molecules with laser intensity 0.28 × 10 14 W/cm 2 . Any differences were at the 10% level or at the high photon energy range of the spectrum where small laser intensity fluctuations translate into large signal changes. 20 40 30 12 14 16 Figure S1. Measured harmonic spectrum from benzene (blue) and deuterated benzene (orange) with a laser intensity of 0.28 × 10 14 W/cm 2 and wavelength 1.8 µm.

S2 Fits to harmonic spectrum laser intensity scan
The measured laser-intensity-dependent harmonic spectra are fit to eq. (1) in the main text using using an iterative procedure. Each iteration consists of minimizing the error with respect to one of the three factors while holding the other two constant. The algorithm cycles through the three factors until convergence; this typically takes less than ten minimizations with respect to each factor. Numerically, the static and intensity factors, A(ω) and B(I L ) are sampled at the experimental values of ω and I L respectively. Sinceω takes a unique value for each (ω, I L ), C(ω) must be described using a basis to control the number of degrees of freedom. We use simple linear interpolation: C(ω) is stored on a uniformly spaced vector ofω, and linear interpolation used to evaluate it at the required values ofω. The number of samples N C is between 20-30. This procedure is similar to finding the best rank-1 approximation to the rank-3 tensor S(ω, I L ,ω) 1, 2 . We use stability of the results with respect to small changes in N C as a measure of the quality and robustness of the fit.
The fit is of high quality for all presented intensity scans. Figure S2(a) shows a representative measured dataset and Fig. S2(b) shows its corresponding fit for comparison. All significant features are captured by the fit. For clarity, the data are scaled by the product of the static (A(ω)) and intensity-dependent (B(I L )) factors, which means the experimental data at high photon energy but low laser intensity appears very noisy in Fig. S2(a). For a more quantitative comparison, Fig. S2(c) compares the experimental and fitted spectra (blue and red) at 18 TW/cm 2 , a mid-range intensity. The two curves are nearly indistinguishable -significant differences emerge only below 10 TW/cm 2 due to low signal. Figure S2 Although the focus of this article is on dynamics, it is useful to compare the ionization factor B(I L ) obtained from experimental and simulated data. This is shown in Fig. S3.

S3 High-harmonic generation simulation
We use a form of the SFA 3 extended to high-harmonic generation in molecules by summing coherently over multiple cation states [4][5][6][7][8] . For completeness, our model also includes laser-induced electronic dynamics in the cation between ionization and recombination 9 . While it is in principle possible to include a full description of the vibronic dynamics of the cation in an SFA calculation 10 , such calculations are at present prohibitively expensive.
In the following, boldface symbols represent spatial vectors, and a quantity with an overbar indicates it is a vector or matrix running over the cation states. For a particular molecular orientation, the high-frequency part of the dipole response at recombination time t r is The equation is an integral over ionization times t i and the integrand, a matrix product, can be read from right to left. The factor describing the angular dependence of the ionization rate is a column vector over cation states the trajectory duration τ = t r − t i , and E E E(t) and A A A(t) are the laser electric field and vector potential respectively. The τ here is the same as the Methods section of the main text, as the time window over which the cation evolves is given by the trajectory duration. The next term in the integrand of (1) is a matrix describing evolution of the cation states under the influence of the laser field. It is the solution of the Schrödinger equation where the diagonal matrixĪ p contains the ionization potentials I p , . . . on its main diagonal, ∆μ µ µ is a matrix of cation state dipole couplings (the dot product E E E(t) · ∆μ µ µ runs over the spatial coordinates) andĪ is the identity matrix. The recombination step is described in (1) byd d d r (v v v), a column vector of recombination dipoles in the same format asd d d i (v v v), and † denotes the conjugate transpose. Finally the scalar exponential contains the standard SFA action 2 2 dt (6) and the τ −3/2 factor describes spreading of the continuum electron wavepacket. The ionization dipole amplitude is the length-gauge plane-wave approximation where |ψ n is the Dyson orbital corresponding to cation state n. The projection operator in (7) arises from our use of the dressed form of the SFA, as described by eq. (33) of Smirnova et al. 11 . This form is necessary to ensure invariance with respect to translation of the molecule. For recombination we use the plane-wave approximation We have verified that substituting spectrally-filtered photoionization cross-sections does not significantly affect the retrieved parameters (data not shown). Our calculations use the five lowest lying cation states. These are well approximated by a single hole in the neutral (the Koopmans' picture 12,13 ), and therefore throughout this work we equivalently refer to them by their corresponding molecular orbital. For the purposes of HHG, we checked this equivalence directly by comparing the spectra produced using Dyson and natural orbitals. The Dyson orbitals were calculated within the static exchange (i.e. Hartree-Fock) model 14 . The calculation included a diffuse basis of continuum orbitals with angular momentum up to L = 6. The inclusion of the diffuse continuum functions has proved important for an accurate description of the exponential tail of the Dyson orbitals 15 . The natural orbitals (orbitals diagonalising the 1-electron density matrix) were calculated with the complete active space self-consistent field (CASSCF) electronic structure method state-averaging over the two equally-weighted lowest-energy cationic eigenstates. We used an active set consisting of the 6 π orbitals and 6-31G* basis set. We found insignificant differences between the harmonics calculated using the two sets of orbitals. The factor g(τ) is the nuclear autocorrelation function which here takes Gaussian form as described in the main text. Because the solution of equation (4) must be evaluated for every pair of ionization time t i and observation time t, an efficient numerical implementation is essential. Because laser-induced cation state coupling plays a small role here and the laser is linearly polarized, we approximate the matrix exponential solution to (4) using a split-operator formula 16 Here,Ī is the ionization potential operator corresponding to half the time between ionization and observation andê is a unit vector parallel to the laser electric field. The matrix of dipole components along the electric field directionê e e · ∆μ µ µ is independent of t and t i . Therefore it can be diagonalized once for each orientation of a given molecule and its exponential evaluated quickly using only scalar exponentials. A higher-order splitting would increase the accuracy but we found that equation (10) was sufficient for convergence. The diagonal elements of the cation dipole matrix ∆μ µ µ lead to dynamic Stark shifts. These elements are equal to the difference between the dipole moments of the cation states and the neutral [17][18][19] . The off-diagonal elements lead to laser-driven transitions between the cation states. In the Koopmans' approximation, the transition dipole ∆μ µ µ is given by the single-particle matrix element: ∆μ µ µ m,n = ψ m | x x x |ψ n .
The transition dipoles from (12) are all within 0.04 au of transition dipoles calculated from the fully-relaxed cationic wavefunctions.
We have tested our model for invariance with respect to translation of the molecule 11 ; to this end inclusion of dynamic Stark shifts is essential. We take the centroid of the benzene ring as the origin in all calculations.
Macroscopic trajectory selection is modelled by coherently averaging across the plane of the laser focus, which is taken as spatially and temporally Gaussian. For the temporal integration we use the adiabatic approximation 20 . The macroscopic harmonic amplitudes are taken as the amplitude of the emission on-axis (k ⊥ = 0) and at the centre of each harmonic (ω = qω 1 ), which strongly favours the short trajectories. In this approximation the spectra are independent of the laser waist size and pulse duration except for a global scale factor.