Full-dimensional treatment of short-time vibronic dynamics in molecular high-harmonics generation process in methane

We present derivation and implementation of the Multi-Configurational Strong-Field Approximation with Gaussian nuclear Wave Packets (MC-SFA-GWP) -- a version of the molecular strong-field approximation which treats all electronic and nuclear degrees of freedom, including their correlations, quantum-mechanically. The technique allows, for the first time, realistic simulation of high-harmonic emission in polyatomic molecules without invoking reduced-dimensionality models for the nuclear motion or the electronic structure. We use MC-SFA-GWP to model isotope effects in high-harmonics generation (HHG) spectroscopy of methane. The HHG emission in this molecule transiently involves strongly vibronically-coupled $^2F_2$ electronic state of the $\rm CH_4^+$ cation. We show that the isotopic HHG ratio in methane contains signatures of: a) field-free vibronic dynamics at the conical intersection (CI); b) resonant features in the recombination cross-sections; c) laser-driven bound-state dynamics; as well as d) the well-known short-time Gaussian decay of the emission. We assign the intrinsic vibronic feature (a) to a relatively long-lived ($\ge4$ fs) vibronic wave packet of the singly-excited $\nu_4$ ($t_2$) and $\nu_2$ ($e$) vibrational modes, strongly coupled to the components of the $^2F_2$ electronic state. We demonstrate that these physical effects differ in their dependence on the wavelength, intensity, and duration of the driving pulse, allowing them to be disentangled. We thus show that HHG spectroscopy provides a versatile tool for exploring both conical intersections and resonant features in photorecombination matrix elements in the regime not easily accessible with other techniques.

HHG is a highly non-linear process [9,39], which often makes the interpretation of experimental results complex * Serguei.Patchkovskii@mbi-berlin.de and controversial. Similar to conventional spectroscopies, measurement of the isotope effects and quantum-path interferences in HHG has emerged as one of the most helpful techniques in the detailed analysis of high-harmonics spectra [15,[40][41][42][43][44][45][46][47][48][49]. In the form of the PACER (Probing Attosecond Dynamics by Chirp-Encoded Recollision [50]) experiment, isotope-dependent HHS promises direct access to the electronic and nuclear dynamics on a fewand sub-femtosecond time scale, including measurement of the phases of the nuclear wave packets [44,45].
For the electronically simple cases, where the Born-Oppenheimer (BO) variable separation applies, the theory underlying PACER experiments is well understood [40,51]. The relevant theory object, the short-time nuclear autocorrelation function [40], is readily available for most small molecules [52]. At the same time, the BO assumptions fail already for one of the first molecules examined with PACER spectroscopy. In methane, CH 4 , the key dynamics occur in the vicinity of a symmetry-required triple electronic surface intersection in a transiently-prepared CH + 4 cation. In larger, more electronically complex systems, which are beginning to be explored with HHS [53][54][55][56][57][58][59][60] (See Ref. [61] for a review), conical intersections and non-separable vibronic dynamics are expected to become ubiquitous [62,63].
The goal of this work is to extend the strong-field approximation (SFA) [39,98,99], a non-perturbative theory underlying much of attosecond science [2,100], to account for the short-time vibronic dynamics. The rest of this work is structured as follows. Section II derives the working expressions for the Multi-Configurational Strong-Field Approximation with Gaussian nuclear Wave Packets (MC-SFA-GWP). Section III gives technical details for the MC-SFA-GWP calculations of HHG emission in methane and deuterated methane, including the details of the electronic structure calculations, the diabatization procedure, and integration of the MC-SFA-GWP equations. Section IV presents the results of the numerical calculations. Specifically, section IV A discusses vibronic nuclear autocorrelation function in methane and analyses its features in terms of the wave packet composition and dynamics. Section IV B presents calculated HHG spectra for three wavelengths (800, 1600, and 2400 nm) and two intensities (300 and 1000 TW cm −2 ) of the driving field, and dissects the physical origin of the observed spectral features. Section IV C analyses calculated PACER signals, compares them to the available experimental data, and makes predictions for the previously unexplored range of the experimental parameters. Section V summarizes the results and offers an outlook for future developments. Finally, some of the technical aspects of the MC-SFA-GWP derivation and implementation are relegated to the appendices A-C.

II. THEORY
Our aim is to model high-harmonics emission due to the short-time, coupled electronic and nuclear dynamics of a molecular system under the influence of an intense, long-wavelength (near-to mid-IR) laser field. It is convenient to write the total HamiltonianĤ T as a sum of three terms, all of which may be time-dependent [101][102][103]: H T (t) =P 0Ĥ0 (t)P 0 +P 1Ĥ1 (t)P 1 +P 1VL (t)P 0 , (1) whereP 0 andP 1 are respectively projectors onto the neutral and singly-ionized spaces. The contributionĤ 0 (t) acts within the subspace of the neutral states. The termĤ 1 (t) acts within the subspace of singly-ionized states, while the interaction HamiltonianV L describes laser-induced transitions between the two spaces. In the dipole approximation and length gauge: where r i are coordinates of the i-th electron and F (t) is the laser electric field, which is is assumed to be slowlyvarying. N is the total number of electrons in the system. In Eq. 1, we have already neglected the possibility of multiple ionization, as well as the stimulated recombination process.
A formally exact solution of the time-dependent Schrödinger equation with Hamiltonian 1 is given by [99,102,103]: The propagatorsÛ 0 ,Û 1 , andÛ T are: where X = 0, 1, T. In Eq. 3, r and q are respectively electronic and nuclear coordinates. We assume that the initial wavefunction at time t 0 (Ψ ( r, q, t 0 )) corresponds to the neutral species.
Classically, the high-harmonics intensity is determined by the power spectrum of the time-dependent dipole d (t). Neglecting contributions arising entirely within the neutral and cation manifolds, the radiating dipole is given by: where the recombination dipole operator is: with sum running over coordinates of all electrons. Equation 5 can be brought into a computationally tractable form by introducing an identity operator1, given by a sum of identity operators within the neutral, singly-ionized, etc. spaces: The identity operator within the neutral space (1 0 ) is defined in terms of the neutral electronic states |Φ a (which depend on r and parametrically on q) and harmonic nuclear vibrational states |n (which depend only on q). Quantity n is a vector of non-negative integers, with each element n k defining the excitation level of the k-th normal mode of a reference potential v ref ( q), which is chosen for computational convenience [52,90].
The singly-ionized space identity operator (1 1 ) contains bound states of the (N −1)-electron ion |X b and the corresponding one-electron scattering states | k b with the asymptotic momentum k b . Functions | k b depend on the coordinate of the N -th electron r N . Formally, sum over b in Eq. 9 does not include autoionizing states |X b , which belong to the doubly-ionized continuum. Such states are unlikely to be populated by strong-field ionization, but may become important due to electron-correlation effects in the recombination step [30,104,105]. If the lifetime of such autoionizing state is long compared to the laser cycle, its effects can still be accounted for, by including the state in Eq. 9.
Both |X b and | k b depend parametrically on q. Op-eratorÂ 1 antisymmetrizes the product wavefunction |X b | k b with respect to the coordinates of the continuum electron [106,107]: whereÎ is the N -electron identity operator, and opera-torP iN permutes coordinates of the i-th and N -th electrons. Similar to the neutral-space identity1 0 , the nuclear wavefunction in1 1 is expanded in terms of the harmonic vibrational states |m of the reference potential v ref ( q). Because we are not interested in processes involving multiple ionization, the sum in eq. 7 is truncated after1 1 . Operator1, as well as the individual1 i operators are idempotent. Identity operators within different subspaces are assumed to commute: This condition is equivalent to the strong orthogonality assumption for all | k b : where |φ D ba is the Dyson orbital for the ionization of the neutral state |Φ a , forming cation |X b .
We now insert operator1 of the Eq. 7 to the left and to the right of each propagator and interaction operator in Eq. 5. Following the usual SFA assumptions [99,102,103], we also replaceÛ T (t, t 1 ) by a product of the (N − 1)-electron ion propagatorÛ I (t, t 1 ) and the Volkov propagatorÛ V (t, t 1 ): where A (t) is the vector-potential of the laser field, and Volkov states | k are taken in the length gauge. The explicit coordinate representation of | k normalized to δ k − k is given by Eq. 16 Neglecting contributions due to correlation-driven inelastic electron scattering in the continuum [106,107], we obtain: C an (t) = a n C ana n (t, t 0 ) C a n (t 0 ) , C a n (t 0 ) = n | Φ a |Ψ ( r, q, t 0 ) , In eq. 18, R bman k is the dipole matrix element for a bound to continuum vibronic transition. The first term in eq. 18 ( k| r|φ D ba ) arises due to the direct interaction between laser field and the active electron. The second contribution ( k| φ C ba ) is due to exchange-mediated interaction of the laser field with the inactive core electrons [108][109][110]. It is described by the 3-component "cradle" vector orbital of eq. 19 (after Newton's cradle, where a force acting on one ball in a multi-ball pendulum causes a different ball to swing [106]). Coefficients C an (t) describe the vibronic wave packet on the neutral surface at time t. Coefficients C ana n (t, t 0 ) and D b m bm (t, t 1 ) are respectively the matrix elements of the neutral and ionic propagatorsÛ 0 (t, t 0 ) andÛ I (t, t 1 ), where n is the energy of harmonic vibrational state |n . The explicit form of the arbitrary-order Taylor expansion of these matrix elements was given in Ref. [90] and need not be repeated here. Equivalently, these matrix elements can also be obtained with MCTDH time propagation [91][92][93][94][95]97]. The phase factors containing n in eqs. 22-23 arise to compensate for the presence of the vibrational phase factor in the definition of the wavefunction Ansatz in ref. [90] (cf. eq. 2 of Ref. [90]). In eqs. 22 and 23, we have extracted the rapidly-oscillating overall phase of the matrix elements, using E N and E I as the characteristic energy of the neutral and cationic manifolds, respectively. Finally, φ d (Eq. 24) is the rapidly-varying part of the HHG dipole phase accumulated from t 1 to t. The first term is due to the bound-state dynamics in the cation, with I p being the characteristic ionization potential; the second contribution is the Volkov phase of the continuum electron.
Applying the usual stationary-phase approximation [102,103] for the k integral in eq. 17, we obtain: where k s is the stationary electron momentum at the time of ionization t 1 .
Direct application of the stationary-phase approximation to the dt 1 integral in eq. 26 leads to a generally complex stationary ionization time t s , causing difficul-ties when the matrix elements of eq. 18 are only known numerically, on the the real axis. An alternative, real treatment is based on an observation that for all but the smallest ∆t = t − t 1 values, the phase φ d (Eq. 24) is dominated by the Volkov phase [98,102,103]. If the term containing the I p can be neglected and the field is linearly polarized, the stationary ionization time t s is then determined by the condition k s (t, t s ) = 0. Small, non-zero I p values can then be accommodated using a Taylor expansion around k s = 0 point [98]. For elliptically polarized fields, k s = 0 can not satisfy eq. 27, as a small transverse initial momentum is required to bring the stationary trajectory back to the origin. The general condition determining the stationary ionization time then becomes: Eq. 28 together with eq. 27 define the real ionization time t s .
An additional complication arises when the vibronic wavefunction is even with respect to the origin and the laser field is linearly polarized. Then, the matrix element R bma n k s in eq. 26 vanishes at k s = 0, and the next order in the power series determines the overall integral over dt 1 . Expanding the ionization dipole through the first order around k s and assuming that the remaining matrix elements in Eq. 26 vary slowly with k s , we obtain (also see Appendices A and B): where E a n and E bm are respectively energies of the neutral and cationic vibronic basis functions. In Eq. 30 we neglected the explicit time dependence of F , which is in the same order as other terms omitted in deriving eq. 30 (See Appendix A). We have also neglected the k s dependence of the recombination matrix elements R * b m a n . If necessary (for example, in the vicinity of a Cooper minimum), this dependence can be reintroduced in a similar manner.
Eq. 29 is our final working equation. A closely-related special-case result was previously obtained by Chirilǎ and Lein [65], and by Madsen et al. [86], which however did not consider the possibility of state crossings and non-adiabatic vibronic coupling. The special-case treatment of SF 6 by Walters et al [83] is similar in spirit, and in principle includes all effects considered here. Another closely-related result by Faisal [84] includes rotational motion (neglected here), but restricts vibrational dynamics to a single surface. In contrast, the treatment of Ref. [88], while superficially similar, completely neglects all subcycle effects arising due to the motion on the intermediate cationic surfaces. These effects are also completely neglected by Ref. [87], which however treats the vibronic dynamics in the neutral manifold. Finally, treatments of Patchkovskii and Schuurman [90], Varandas et al [91][92][93][94][95], and Arnold et al [97] fully treat cationic vibronic dynamics, but neglect neutral dynamics and their coupling to the HHG process.

A. Electronic structure calculations
All electronic structure calculations used correlationconsistent, valence triple-ζ basis set, augmented with diffuse functions ("aug-cc-pVTZ" [111,112]). Calculations of the neutral CH 4 geometry and vibrational Hessian used second-order Møller-Plesset (MP2) correlation treatment, with the 1s electrons on the carbon atom left uncorrelated. Optimized C-H bond length (r 0 = 1.110Å) and unscaled harmonic vibrational frequencies (ν 1 ...ν 4 = 3069, 1589, 3204, 1356 cm −1 ) are in a good agreement with the experimental equilibrium geometry and fundamental frequencies (r e = 1.094Å, ν 1 ...ν 4 = 2916, 1534, 3019, 1306 cm −1 [113]). The equilibrium neutral geometry corresponds to a symmetry-required triple degeneracy point of the 2 F 2 electronic ground state of the methane cation [114]. The electronic structure in the vicinity of this degeneracy point was explored by displacing all pairs of Cartesian coordinates by ±0.02Å, for the total of 450 unique distorted structures. The reference geometry was taken in the standard setting of the T d group, with the C 2 axes along the main Cartesian directions (See Fig. 1). Symmetry was used neither to reduce the number of the displaced geometries, nor in the electronic structure calculations at these geometries.
Electronic wavefunctions of the three lowest cation states correlating to the 2 F 2 state were calculated using multi-reference configuration interaction with single excitations treatment (MR-CIS). All configurations within the minimal complete active space with 7 electrons in 4 frontier orbitals [CAS (7,4)] were included in the reference set. Single-particle orbitals were optimized in a CAS (7,4) self-consistent field (CASSCF) calculation, with the energies of the three cation states correlating to 2 F 2 weighted equally. The neutral wavefunctions were determined within the CIS space, using the lowest closed-shell determinant constructed from the cation-optimized orbitals as the reference. The MR-CIS vertical ionization potential (MR-CIS: 13.14 eV) slightly underestimates the experimental value (Expt: 13.60 eV [115,116]). Including double excitations in the CI expansion (MR-CISD; Langhoff-Davidson +Q correction [117] not included) leads to an overestimation of the vertical IP by a similar amount (MR-CISD: 13.97 eV). Given the large number of the distorted geometries, further increase in the CI size is not practical. It is in principle possible to adjust the calculated vertical IP to match the experimental value. However, the experimental vertical ionization potential is somewhat uncertain, with a very very broad, complex peak in the 13.2-14.0 eV range [116,118]. Because all qualitative features of the energy surfaces are already adequately described at the MR-CIS level, all subsequent calculations use MR-CIS wavefunctions and energies. Due to the underestimation of the vertical IP, we expect ionization rates (Eq. 30) to be somewhat overestimated.

B. Diabatization procedure
Evaluation of matrix elements entering Eq. 29 requires that both the potential energy surfaces and and the wavefunctions are smooth and continuous in the vicinity of the expansion point. It is therefore necessary to diabatize the MR-CIS states. Construction of quasi-diabatic states is a long-standing, and still active area of research (See Ref. [63] for an overview), with many competing prescriptions available in the literature. An appealing, generallyapplicable approach for constructing local quasi-diabatic surfaces involves fitting adiabatic energy gradients and derivative couplings at selected points in the nuclear coordinate space [119][120][121]. Unfortunately, this technique does not offer an easy access to the transformed quasidiabatic wavefunctions, which are needed for evaluating the somewhat non-standard matrix elements Υ and R (Eqs. 18,30). Here, we retain the quadratic vibronic Hamiltonian form of Ref. [119], but revert to the older maximum-overlap approach [122][123][124] to directly determine the diabatic wavefunctions: where Ψ dia ja is the j-th quasi-diabatic state at geometry a, Ψ ad ia is the i-th adiabatic state at the same geometry. The unitary transformation U jia connects the two wavefunction spaces. At each displaced geometry (a), the optimal matrix U a is determined [122][123][124][125][126][127] by the overlap matrix S a between the adiabatic wavefunctions at the displaced geometry and reference wavefunctions Ψ ref k (see Fig. 1): For methane, reference wavefunctions Ψ ref are taken at the high-symmetry T d geometry, treated within the D 2 sub-group. Neutral manifold uses the 1 A 1 ground-state wavefunction, while cationic manifold is referenced to the lowest 2 B 1 , 2 B 2 , and 2 B 3 wavefunctions. Once the diabatization transformation U a is determined, the quadratic vibronic Hamiltonian is obtained with a least-squares fit to the diabatic Hamiltonian matrices at the displaced geometries. Electronic matrix elements are calculated directly from the transformed diabatic wavefunctions (eq. 34) as described elsewhere [110,128] (also see Appendix C). For the methane molecule in the vicinity of the equilibrium neutral geometry, the resulting vibronic Hamiltonian coincides with the result of [119], while the residual derivative couplings are found to be numerically negligible. Direct non-linear minimization of the quadratic vibronic Hamiltonian fit with respect to the diabatization parameters U a leaves the result unchanged, confirming that the transformation of Eq. 36 represents a local optimum.
It should be noted that the quadratic vibronic Hamiltonian of the 2 F 2 manifold of CH + 4 determined presently focuses on the small region of configuration space in the vicinity of the neutral equilibrium geometry. It is unbound from below, and corresponds to a dissociative state. As a result, it is not capable of describing longtime dynamics of the cation, and should not be compared globally to the more conventional, spectroscopic surfaces, such as used in Refs. [91,92].

C. Dyson orbitals and dipole matrix elements
Dyson orbitals (Eq. 13) corresponding to the groundstate ( 2 F 2 ) ionization channel of the CH 4 molecule are shown in Fig. 1. All three components are strongly 1electron allowed, with the Dyson orbital norm of 0.916. Visually, these Dyson orbitals are indistinguishable from the triply-degenerate HOMO of the neutral species. The dominant "cradle" terms (Eq. 19; not shown) are derived from the (s-like) 2a 1 molecular orbital of the methane molecule.
It is instructive to examine the electronic part of the photoionization matrix elements: (see Eq. 18) a little closer. Due to the high symmetry of the CH 4 molecule, it is sufficient to constrain the photoelectron observation direction to one of the Cartesian axes (See Figs. 1 and 2). The dominant first term in the matrix element is identical to the first Born approximation in photoionization, known to be qualitatively inaccurate close to the ionization threshold. Beyond 30 eV photon energy, however, the calculated cross-sections are in a (fortuitously) good agreement with experimental data [116,[129][130][131][132][133] and accurate calculations [134][135][136]. For example, Eq. 37 yields photoionization crosssection of ≈ 17 Mbarn at 30 eV photon energy (expt: ≈ 14Mbarn [136]), decreasing to ≈ 5 Mbarn at 60 eV The "cradle" terms become important [108] beyond ≈ 150 eV photon energy (See Fig. 2), leading to an ≈ 0.4 kbarn minimum in the calculated cross-sections at ≈ 196 eV. It is followed by a maximum (≈ 7 kbarn) at 295 eV, (again fortuitously) close to where the crosssection is expected to increase due to the K-edge intensity borrowing [104]. Unfortunately, we are not aware of reliable final-state resolved calculations or measurements in methane for photon energies above 100 eV, where ionization from the inner 2a 1 and 1a 1 shells dominates the overall cross-section [116]. It is therefore unclear whether the photoionization matrix elements in Fig. 2 provide a reasonable description of methane photoionization beyond 90 eV. Nonetheless, the minimum and the associated π phase jump at 196 eV serve to illustrate an important point (see below), which remains qualitatively valid even it occurs at a different energy in the actual CH 4 molecule than in our crude calculations here.

D. Nuclear autocorrelation and MC-SFA-GWP HHG calculations
Nuclear autocorrelation functions are a special case of the vibronic matrix elements D b m bm (t, 0) of Eq. 23, with b = b and m = m. Both D b m bm and C a n matrix elements are propagated in time numerically, using the 4-th order Runge-Kutta integrator with a uniform time step [137]. Calculations of the autocorrelation function used time step of 0.02 au[t] (≈ 0.48 as). Shorttime autocorrelation functions calculated presently numerically coincide with the analytical results [90], for all times where the power-series expansion in Ref. [90] converges.
Calculations of the time-dependent MC-SFA-GWP dipole d (t) (Eq. 29) used a time step corresponding to 1/3 of the Nyquist limit for the cut-off harmonics. The resulting time step ranged from 0.386 au[t] (≈ 9.3 as) for the 800 nm, 300 TW cm −2 driving field, to 0.037 au[t] (≈ 0.89 as) for the 1.6µm, 1 PW cm −2 driver. We did verify that all calculations are converged with respect to the time step. The high-harmonics spectrum was calculated from the Fourier transform of the time-dependent dipole [102,103]. No window function was applied.
The vibronic wavefunction expansion was performed in the basis of multi-dimensional harmonic oscillator functions of the neutral species. Rotational and translational modes were excluded. Calculations of the autocorrelation functions allowed up to 12 quanta in any of the remaining 9 normal modes. This choice is fully converged with respect to the basis set up to at least 4.5 fs evolution time, but is too computationally expensive for MC-SFA-GWP calculations, where the vibronic matrix elements need to be recalculated for each trajectory. Instead, MC-SFA-GWP calculations limited the vibronic basis to at most 4 quanta in any of the vibrational modes. This choice changes the autocorrelation function by less than 0.2% at 1.7 fs delay (800-nm cutoff trajectory), increasing to 3% at 2.6 fs (1200-nm cutoff) and 25% at 3.4 fs (1600-nm cut-off). Our MC-SFA-GWP results are therefore sufficiently converged at 800 and 1200 nm, but are only qualitative for the high-energy part of the spectrum at 1800 nm (beyond 160/500 eV at 300/1000 TW cm −2 ).
All MC-SFA-GWP calculations were performed for laser electric field polarized along the molecular X axis. Due to the high symmetry of the CH 4 molecule, we expect orientational averaging to play only a minor role for this system. Only the short trajectories were considered.

A. Free-cation dynamics and autocorrelation function
The simplest summary of the subcycle nuclear dynamics is provided by the nuclear autocorrelation function, shown in Fig. 3 for both CH 4 and CD 4 cations, using neutral vibrational ground state as the initial wavefunction at time zero. For comparison, we show autocorrelation functions for the same species and initial conditions from Ref. [52]. At very short times (up to ≈ 1.2 fs), all four autocorrelation functions are Gaussian in time (please note the linear-log scale of the Fig. 3a), which is a generic feature of short-time autocorrelation functions [52,138]. The agreement between the data sets is remarkably good, con- sidering the differences in the theoretical treatment. In Ref. [52], the autocorrelation functions are calculated on a single component of the 2 F 2 electronic surface, with the degeneracy lifted by applying an intense (≈ 2.6 VÅ −1 ) static electric field. Because the components of the 2 F 2 state are coupled by the electric field (MR-CIS transition dipole of ≈ 0.7 Debye), the surface shape is also modified, allowing the single-surface dynamics to mimic the more complex vibronic dynamics at early times.
At the same time, the single-surface treatment does not allow for the "orbiting" motion of the wave packet around the conical intersection point, precluding the possibility of revivals on a few-femtosecond time scale. Beyond ≈ 1.2 fs, a single, field-diabatized surface no longer adequately represents vibronic dynamics in this system. For CH 4 (CD 4 ), the population of the initial surface reaches a minimum of 52% (51%) at 1.5 fs (1.7 fs), then begins to increase again as a fraction of the wave packet completes a half-revolution around the CI. Correspondingly, the autocorrelation function changes sign (See Fig. 3b) and undergoes a half-revival at 2.15 fs (2.50 fs). For both CH 4 and CD 4 , the nuclear autocorrelation factor in HHG reaches 1.5% at the half-revival -more than an order of magnitude above the expected damping factor due to the Gaussian decay [52]. A full wave packet revival appears as a shoulder at the 0.5% level, at ≈ 3.6 fs (≈ 4.2 fs). Similar, non-Gaussian intermediate-time decay and unexpected revivals has been previously predicted in other non-adiabatic systems [89,90,97].
The short-time autocorrelation functions in Fig. 3 are qualitatively consistent with the previously reported calculations by Mondal and Varandas [91,92]. These authors find the first half-revival at ≈ 2.4 fs for CH + 4 (≈ 2.8 fs for CD + 4 ). The full revivals are also found at times comparable to our results. At the same time, the maximum of the CD 4 /CH 4 autocorrelation functions ratio calculated in Ref. [92] is much higher (≈ 6, compared to ≈ 2.7 here. It is also reached at a later time (≈ 1.85 fs, compared to ≈ 1.5 fs here). Both discrepancies are well within the range of available experimental resolution, and would be interesting to explore.
Furthermore, in contrast to Ref. [91,92], we cannot attribute the revivals to the dynamics reaching the Jahn-Teller minima on the cationic energy surface. Calculated expectation values of normal coordinates on each electronic surface evolve monotonically at least until 3.5 fs. The Cartesian displacements in the centre of mass coordinate system remain small (less than 0.02Å), with the structure remaining at the D 2d symmetry, indicating that the minimum is not reached at the times relevant for the PACER experiments at mid-IR wavelengths. It is therefore instructive to examine the origin of the oscillations in the short-time autocorrelation function in a little more detail. The square modulus of the autocorrelation function measures the population of the initiallypopulated vibronic basis function (D 1 |gs in Fig. 4). The strongest vibrational coupling of |gs on the same quasidiabatic electronic surface is to a vibrational basis function with a single quantum in the ν 2 vibrational mode (|ν 2 : e symmetry, 1589 cm −1 , H − C − H bending motion). The strongest coupling to the other two quasidiabatic surfaces (D 2 and D 3 ) involve single excitation of the ν 4 (|ν 4 : t 2 , 1356 cm −1 , H − C − H bend) and ν 3 (|ν 3 : t 2 , 3204 cm −1 , C − H stretch) normal modes. All of these couplings are of a similar magnitude. The number of significant couplings increases rapidly with the vibrational excitation level; once two or more vibrational quanta have been excited, the probability of the wave packet returning to the initial state D 1 |gs becomes negligible.
The short-time evolution of the population of the initially-populated and singly-excited vibronic basis functions in CH + 4 is illustrated in Fig. 4. One can immediately see that the populations of all five dominant singlyexcited functions (D 1 |ν 2 , D 2 /D 3 |ν 4 , and D 2 /D 3 |ν 3 ) oscillate in-phase at least up to 4.5 fs. The oscillations are out-of-phase with the autocorrelation magnitude (D 1 |gs population), indicating the back-and-forth population transfer. This behavior is easy to rationalize by considering a reduced-dimensionality model, consisting of just the six essential states listed above. The corresponding model Hamiltonian is: where all elements not in the first row or column, or on the main diagonal, are zero. As long as all d i ≈ d, this Hamiltonian is equivalent to a two-level system: where the transformed basis is: The φ 2 basis function corresponds to a vibronic wave packet component "orbiting" the conical intersection. In the nuclear coordinate space, |φ 2 forms a prolate spheroidal shell around the conical intersection. The oscillatory behavior in the autocorrelation function (ie the population of the |φ 1 basis function) is due to the interference between the two eigenstates of Eq. 39, while the overall decay of the signal in this model arises from the imaginary part of d ( (d) < 0). To the order O(d 2 ), the difference between the eigenvalues of H mod,2 is: Thus, the oscillation frequency of the short-time autocorrelation function is determined by the strength of the coupling between the initially-prepared component of the wave packet and its decaying part "orbiting" the conical intersection.
As can be clearly seen from Eq. 40, non-adiabatic vibronic coupling between the electronic surfaces plays a key role in the initial wave packet decay. If vibronic coupling is neglected, the initial decay is slowed down dramatically (data not shown), as was seen before for the benzene cation [90]. The resulting lower isotope effects are no longer compatible with experimental data (See Section IV C below). A similar observation was made in Ref. [86], which neglected the vibronic coupling.  Coordinate dependence of the tunneling-ionization matrix elements (eq. 30) can lead to substantial reshaping of the neutral vibrational wave packet upon ionization [139]. In methane, ionization primarily populates ν 4 (asymmetric bend) and ν 3 (asymmetric stretch) modes. Nuclear autocorrelation functions for single-quantum excitation of these two modes in CH 4 and CD 4 are shown in Fig. 5. These autocorrelations are qualitatively similar to the case of the vibrational ground state (Fig. 3). However, the first half-revival occurs at an earlier time, and has a substantially higher magnitude (CH 4 : 2.9/2.5% at 1.93/1.98 fs for ν 4 /ν 3 ; CD 4 : 2.8/2.7% at 2.23/2.25 fs for ν 4 /ν 3 respectively). The faster short-time dynamics of vibrationally excited states has been noted before [52,65,67,85,94]. This property will be important for understanding the HHG spectra at longer wavelengths (Section IV B below).
Overall, non-adiabatic nuclear autocorrelation factors suggest that HHG emission in methane should persist at much longer times than expected previously. Even at 4.3 fs, close to the short-trajectory cut-off for a 2µm IR driver, the autocorrelation factor remains above 0.2% (0.4% for CD 4 ) -within the possible measurement range. Furthermore, the non-adiabatic dynamics around the CI can be expected to lead to a rich PACER signal, including inverse isotope effects between 1.8 and 2.3 fs. A more detailed investigation of the HHG signal therefore appears justified, and is attended to in Sections IV B and IV C below.

B. High-harmonics generation
For calculations of high-harmonics spectra, we consider three driver wavelengths (800 nm, 1.2 µm, and 1.6 µm) and two peak intensities (300 TW cm −2 and 1 PW cm −2 ). For each combination of the wavelength, intensity, and isotopic species, we perform four simulations, namely: "F": In this simulation, the nuclei are "frozen", and are not allowed to move between ionization and recollision. To this end, the vibronic matrix elements of eqs. 22 and 23 are replaced by: All other matrix elements and the initial nuclear wave packet, represented by the C a n (t 0 ) coefficients in Eq. 20, remain unchanged. As the result, this simulation can still exhibit isotope effects due to the differences in the initial ground-state wave packet widths combined with coordinate dependence of the matrix elements.
"AC": In this "autocorrelation" simulation, the electronic part of the dipole is calculated exactly as in the frozen-nuclei simulation "F". For each trajectory, the dipole is multiplied by the magnitude of the nuclear autocorrelation function (section IV A). Thus, electronic continuum and vibronic dynamics are assumed to be factorized and mutually independent. This is the the approach typically taken for the analysis of PACER spectra.
"ND": In this "no-dipole" simulation, we use the full Eq. 29. However, the effects of the laser field are neglected when evaluating vibronic matrix elements of Eqs. 22 and 23. Thus, vibronic dynamics in the neutral species and the cation are taken to be field-free, and depend only on the initial composition of the wave packet and the elapsed time.
The laser field interacts with the continuum part of the wavefunction at all times, but with the rest of the molecule only at the moment of ionization, through the Υ bma n matrix element of Eq. 30. Electronic and nuclear dynamics are otherwise fully correlated.
"D": Finally, the "dipole" simulation uses Eq. 29, with no further simplifications to any of the matrix elements. The laser field can cause (subcycle) vibrational excitation in the parent neutral species, modify the nuclear wave packet through the ionization matrix elements Υ bma n , and cause vibronic transitions on the cationic surface. "D" is our besteffort simulation.
Calculated normal and deuterated methane HHG spectra for the 800-nm driving field are shown in Fig. 6. At the moderate, 300 TW cm −2 intensity of the driving field (panel a), harmonic cut-off is found at ≈ 74 eV photon energy (harmonic 47). For these field parameters, the frozen-nuclei spectra of the two isotopic species are indistinguishable. Close to the cut-off, nuclear motion leads to HHG suppression by a factor of ≈ 75× (CD 4 ) or 100× (CH 4 ). All three approaches we consider for treating the nuclear motion yield essentially the same results, indicating that the vibronic and recollision dynamics are largely uncorrelated, and the sub-cycle vibronic dynamics is not affected by the laser field.
The situation changes when the intensity is increased to 1 PW cm −2 (Fig. 6b). Now, the HHG cut-off extends to ≈ 206 eV (H133), past the minimum in the recombination cross-section (See Fig. 2). In the vicinity of the minimum, the finite width of the initial nuclear wave packet now leads to isotope dependence even in the absence of sub-cycle nuclear motion: HHG emission from CD 4 ("CD 4 F", solid black line) in the vicinity of the minimum is ≈ 55× stronger than for CH 4 ("CH 4 F", dotted black line). At a first glance, this result appears counter-intuitive: the broader ground-state vibrational wave packet in CH 4 is expected to sample a wider range of nuclear configuration, and thus better "fill in" the 196 eV structural minimum in the recombination matrix element. However, in the vicinity of the CH 4 harmonic minimum (≈ 194 eV), the phase of the recombination matrix element varies substantially over the characteristic extent of the ground-state nuclear wavefunction. The destructive interference then leads to a much stronger suppression in "CH 4 F" harmonic emission than might have been expected from the matrix elements at the equilibrium geometry alone. The narrower distribution of the nuclear positions in the heavier CD 4 reduces the extent of the destructive interference, leading to increased HHG intensity near the minimum.
The qualitative features of the frozen-nuclei highharmonics spectra remain essentially unchanged at longer wavelengths: the two isotopic species are indistinguishable away from the structural minimum in the recombination matrix elements. In the vicinity of the structural minimum, CD 4 shows much smaller suppression compared to the lighter isotopologue. As the result, we will neither show nor discuss the frozen-nuclei HHG spectra for longer wavelengths.
Returning to the 800-nm, 1 PW cm −2 case (Fig. 6b), the results obtained using different approaches for the treatment of sub-cycle nuclear motion remain similar within the harmonic plateau (below ≈ 150 eV), but start to differ closer to the structural minimum. As expected, the factorized "AC" approach faithfully reproduces the shape of the frozen-nuclei spectra (CH 4 : solid red line; CD 4 : solid blue line, Fig. 6b), and yields a pronounced minimum around 194 eV. However, the minimum is completely filled-in when using a more elaborate treatment ("ND" or "D"), reflecting the spreading of the nuclear wave packet between ionization and recombination. The magnitude of the isotope effects is also much reduced, and is possibly inverted (see Section IV C below). Furthermore, laser-field modification of vibronic dynamics may begin to play a role. Calculated HHG spectra for single-sided recollision at 1.2 µm are collected in Fig. 7. At 300 TW cm −2 , the cutoff is found at ≈ 145 eV (H140), well below the structural minimum. On the low-energy side of the plateau (below ≈ 80 eV), the three approaches to the treatment of nuclear motion ("AC", "ND", and "D") yield very similar results. This is not unexpected: this low-energy part of the spectrum corresponds to recollision time delays below 1.7 fs, which were explored by the 800-nm results above. At longer time delays, correlated vibronic calcu-lations (both "ND" and "D") now predict significantly higher HHG intensity than the direct-product "AC" approximation (but still two orders of magnitude lower than the frozen-nuclei results). The enhancement is due to the coordinate dependence of the strong-field ionization (Eq. 30). Close to the peak of the laser electric field, it leads to a substantial population of the the nuclear basis functions excited along the ν 3 (asymmetric stretch) and ν 4 (asymmetric bend) normal modes. For CH 4 and linear polarization along X (Fig. 1), the initial amplitudes of the singly-excited ν 4 (ν 3 ) on the dominant D 3 electronic surface reach 14% (33%) of the vertical-ionization amplitude at 300 TW cm −2 , increasing to 14% (36%) at 1 PW cm −2 . Ionization-induced nuclear wave packet reshaping is smaller, but still substantial for CD 4 : ν 4 (ν 3 ) relative amplitudes of 9% (19%) at 300 TW cm −2 , increasing to 10% (21%) at 1 PW cm −2 . The initiallyexcited component of the vibronic wave packet reaches the half-revival faster (after ≈ 2 fs delay), while the higher magnitude of the revival (Fig. 5) compensates for the reduced population relative to the vertical ionization. In the heavier CD 4 , the ν 4 /ν 3 half-revivals occur later (after ≈ 2.2 fs), so that the "ND" and "AC" spectra remain similar until much closer to the cut-off. Finally, the "ND" (no laser coupling in the vibronic dynamics) and "D" (full laser coupling) spectra begin to diverge close to the cut-off, indicating that field-induced bound-state dynamics becomes important beyond 2 fs delays.
At the higher 1 PW cm −2 intensity (Fig. 7b), the harmonics cut-off now extends to 443 eV (H429). The structural minimum in the recombination matrix elements is now well within the plateau region, inducing a false cutoff near 200 eV. The direct product form ("AC") remains a reasonable approximation early within the plateau (up to ≈ 1.5 fs; 160 eV emission energy). At higher emission energies, both the correlations between the vibronic and continuum dynamics and the field-induced vibronic transitions become important. The position of the structural minimum is shifted to higher photon energies (225-250 eV, depending on the species and the details of the treatment). A similar shift of an HHG feature was experimentally observed [43] in PACER experiments on H 2 , where however coordinate dependence of the recombination matrix elements moves the minimum to lower energies.
Vibronic dynamics "fills in" the structural minimum ("CD 4 ND": purple line with triangles; "CH 4 ND": orange line with inverted triangles). At the same time, field-induced vibronic transitions sharpen the minimum again ("CD 4 ND": teal line with crosses; "CH 4 ND": pink line with plus signs). The isotope effects in the high-energy part of the plateau thus reflect a number of dynamics (see IV C below).
Our final example uses a 1.6-µm driver (Fig. 8). Already at 300 TW cm −2 , (panel a) the cut-off extends beyond the structural minimum, to 244 eV (H315). The qualitative features of the calculated HHG spectra are similar to the 1.2-µm, 1 PW cm −2 case above: the sim- plified product representation ("AC") is adequate at low energies (short times); correlated vibronic dynamics and laser-induced vibronic coupling are increasingly important at higher energies, and especially in the vicinity of the structural minimum. At the 1 PW cm −2 intensity, the still higher harmonic cut-off (775 eV, H1000), leads to a particularly clean demonstration of a transition from largely pure, field-free vibronic dynamics to fullycorrelated dynamics under strong-field control (Fig. 8b).
On the left side of the plateau (before 170 eV/1.6 fs), the spectra are grouped by the isotopic species, with all three treatments ("AC", "ND", and "D") yielding very similar results for either species.
On the right-hand side of the plateau (beyond 450 eV/2.3 fs), the HHG spectra are instead grouped by the treatment of bound-continuum correlations and laser field effects in vibronic dynamics. For example, at 600 eV, the product form ("CD 4 AC", blue line with squares) underestimates the HHG intensity by a factor of ≈ 3× compared the field-free correlated treatment of the vibronic dynamics ("CD 4 ND", purple line with triangles). At the same time, the full treatment including field-induced vibronic transitions ("CD 4 D", teal line with crosses) predicts HHG yield at this energy ≈ 2.3× below the "AC" treatment, and ≈ 8× below the field-free dynamics prediction. For nuclear dynamics simulations neglecting field-induced vibronic transitions, isotopic effects in this region do not exceed 1.5×. However, the full treatment ("CH 4 D", purple line with plus signs) predicts a very strong suppression of harmonic emission beyond 500 eV, which is likely to appear as a false cut-off in experiment. Clearly, neglecting field-induced bound-state transitions in the cation is not a viable treatment in the high-intensity and multi-cycle limit [65,71].
Finally, in the transition region around the structural minimum, both the position and the depth of the structural minimum are highly sensitive to the isotopic species, treatment of non-adiabatic dynamics, and laserinduced bound-state dynamics. Again similar to the 1.2-µm, 1 PW cm −2 case, correlated vibronic dynamics shifts the apparent position of the structural minimum to higher energies and reduces its contrast, while inclusion of the laser-induced vibronic dynamics partially reverses the shift. In CH 4 , the minimum is found at 198 eV ("AC"), 236 eV ("ND"), and 219 eV ("D"). In CD 4 , the minima are shifted slightly to higher photon energies: 200 eV ("AC"), 239 eV ("ND"), and 220 eV ("D").

C. Isotope effects
Because PACER experiments are routinely interpreted in terms of ionization-recollision time delays, it is instructive to examine the isotopic ratios of the highharmonics spectra in Figs. 6-8 in a similar way. In the absence of continuum resonances, there exists a one-to-one mapping between the short-trajectory HHG spectrum and the ionization-recollision time delay [9]. This relationship breaks down close to sharp features in photorecombination matrix elements (See Section IV B and Refs. [102,103,105,140]). Although the emission time in these cases can still be recovered through the time-frequency analysis [141], there appears to be no unambiguous way of reconstructing the ionizationrecollision time delay close to resonances. We therefore choose the classical "simple man's" mapping [9] (SMM) for the ionization-recombination time delay, with the twin caveats: a) this mapping is known to be inaccurate for short trajectories, especially for the 800-nm driver [7,141]; and b) the mapping should be treated as undefined close to structural minima in the harmonic spectrum. The advantage of the SMM is that the range of possible time delays depends only on the wavelength of the driving laser, so that different intensities can be compared directly. The underlying assumption of the PACER method is that the isotope effects represented in this form are intensity-and wavelength-independent.
Calculated isotope effects for the 800-nm driving field are shown in Fig. 9. The available experimental data [41] (black error bars) was obtained at the estimated intensity of 200 TW cm −2 , and can be most directly compared to  (Fig. 3). Dash-dotted and dotted grey lines in panel a) represent the ratio of autocorrelation factors for ν3 and ν4 vibrationally-excited initial wave packets (Fig. 5), respectively. Panels a) and b) correspond to laser intensity of 300 and 1000 TW cm −2 , respectively. Experimental points in panel a) are from Ref. [41], measured at the estimated intensity of 200 TW cm −2 .
the numerical results of panel (a), calculated at a somewhat higher intensity of 300 TW cm −2 . As expected from the similarity of the calculated spectra for the three approximations we consider here (Fig. 6a), the calculated PACER ratios are nearly identical. The simple ratio of the autocorrelation functions (Fig. 3) yields nearly identical results, except very close to the harmonic cut-off (t > 1.5 fs). The agreement with the experiment is satisfactory, although the experimental isotope ratios are consistently slightly higher than the calculated values. A possible reason for the discrepancy is the vibrational excitation of the neutral molecule by the raising edge of the 8 fs pulse used in experiment [41], which is not included in the present single-cycle simulation. As can be seen from Fig. 5, population of the IR-active ν 3 and ν 4 vibrational modes in the initial wave packet is expected to increase isotope effects for the delays below 1.3 fs (dash-dotted and dotted lines in Fig. 9a). Both the intrinsic correlations and laser-induced vibronic dynamics become important for the isotope effects at the higher, 1 PW cm −2 intensity (Fig. 9b). Compared to the uncoupled approximation ("AC", blue line with diamonds), correlations between the continuum and vibronic dynamics in the cation ("ND", green line with triangles) substantially reduce the calculated isotope effects. Laser-induced vibronic dynamics ("D", red line with squares) partially counteracts this effect. Our best-effort calculation ("D") suggests that the observed isotope effects at this intensity should start decreasing beyond 1.2 fs delays, with inverse isotope effects predicted beyond 1.45 fs. We emphasize that the inverse isotope-effect in this case is due to the structural minimum at 196 eV in our photoionization matrix elements (Fig. 2) and the associated breakdown of the time-frequency mapping. An inverse isotope effect of a similar origin was previously predicted for the D 2 /H 2 pair [65,67]. A change in the position of the photoionization resonance will also change the apparent time delay where the inverse isotope effect is predicted.  Fig. 7. Also see Fig. 9 caption.
The simulated PACER results for the longer-wavelength, 1.2µm driver and moderate 300 TW cm −2 intensity are shown in Fig. 10a. Superficially, this PACER trace is remarkably similar to the "D" (correlated and laser-coupled) result at 800 nm and 1 PW cm −2 (Fig. 9b): the isotope effect first increases, then changes sense at longer delay times. However, the physics behind the trace is entirely different. At 1.2µm, the signal reflects the field-free vibronic dynamics in the transient cation, with neither vibronic-continuum correlations nor laser-driven vibronic dynamics qualitatively affecting the result. Thus, neglecting the laser coupling in the cation ("ND") and neglecting both the laser coupling and vibronic-continuum correlations ("AC") yield results very similar to the full simulation. All three simulated PACER traces are nearly on-top of the simple ratio of the autocorrelation functions (dashed line). At the 1.2µm wavelength, the laser cycle is long enough to allow vibronic wave packet in the cation to reach the halfrevival at the conical-intersection point (See Fig. 3 and Section IV A). The PACER trace containing a reversal of the isotope effect thus represents a true signature of the CI dynamics. The situation becomes more complex at the higher 1 PW cm −2 intensity (Fig. 10b). Now, the harmonic spectrum extends far enough to access the structural minimum in the recombination matrix elements. As the result, the apparent isotope effects at time delays corresponding to the structural minimum ("AC": ≈ 1.63 fs, "ND": ≈ 1.78 fs, "D": ≈ 1.73 fs) become very large, and sensitive to the details of the treatment ("AC": 14×; "ND": 3.9×; "D": 175×). At time delays unaffected by the structural minimum (t < 1.55 fs and t > 1.85 fs), the calculated PACER signal is qualitatively similar to the 300 TW cm −2 results. As was already seen above for the 800 nm case, the influence of continuum-vibronic correlations and laser-driven vibronic dynamics increases at the higher intensity.
Finally, simulated PACER traces at 1.6µm wavelength are shown in Fig. 11. Already at the moderate 300 TW cm −2 intensity, harmonic cutoff reaches beyond the structural minimum. The corresponding PACER trace (panel a) is representative of the intrinsic vibronic dynamics for the delays below ≈ 2.5 fs. At longer times, the structural-minimum feature dominates. Further increasing the intensity (panel b) leads to a very complex spectrum, which loses nearly all information on the intrinsic vibronic dynamics beyond ≈ 1.7 fs. If the intrinsic, field-free dynamics of the cation is of interest, extending the wavelength to 1.6µm therefore appears to offer little advantage.

V. CONCLUSIONS AND OUTLOOK
In this work, we introduced MC-SFA-GWP -a version of the molecular strong-field approximation which treats all electronic and nuclear degrees of freedom, including their correlations, quantum-mechanically. The new technique allows, for the first time, realistic simulations of the nuclear motion effects on high-harmonic emission in polyatomic molecules without invoking reduceddimensionality models for the nuclear motion or the electronic structure.
We use the new technique to model isotope effects in methane. The intermediate 2 F 2 electronic state of the CH + 4 cation transiently accessed by the HHG process possesses a symmetry-required triple-state conical intersection at the Franck-Condon point. Simulations of the field-free vibronic dynamics of the 2 F 2 state in CH + 4 indicates that a fraction of the initially-prepared wave packet undergoes a half-revival, accompanied by a sign change of the vibronic wavefunction, within 2.2 fs. This time scale is well within the laser-cycle duration of near-IR light, and is accessible with HHG spectroscopy. The revival is associated with population of an intermediate vibronic wave packet, composed of the singly-excited ν 4 (t 2 ) and ν 2 (e) vibrational modes coupled to the degenerate components of the 2 F 2 electronic surface. In the space of nuclear coordinates, this wave packet forms a prolate spheroidal shell around the conical intersection. The revival times are determined predominantly by the strength of the vibronic coupling near the CI, thus permitting its direct experimental measurement.
A well-understood difficulty in PACER and other HHG spectroscopies of molecules is the strong (Gaussian in time) suppression of the HHG emission due to the loss of nuclear wavefunction overlap. We demonstrate that beyond 1.5 fs, nuclear motion in methane no longer causes a Gaussian suppression of the HHG signal, with the expected nuclear factor persisting at ≈ 1% level up to 3.5 fs -high enough to allow experimental detection. This unexpected persistence of the fraction of the HHG signal appears to be universal, and has been predicted for several other cations [89,97].
We analyze and identify a number of physical mechanisms which contribute to the isotopic PACER signal in methane. At intensities below 300 TW cm −2 and wavelengths below 1.2µm, the autocorrelation contribution is dominant. This contribution derives from the intrinsic vibronic dynamics around the conical intersection. It manifests as an inverse isotope effect at ≈ 2.2 fs ionization-recollision delay, and is wavelengthand intensity-independent. This effect has been predicted previously [91,92]. An experimental observation was claimed very recently [49].
The autocorrelation term is modified by a number of wavelength-and/or intensity-dependent contributions, including: coordinate dependence of strong-field ionization amplitudes [64]; laser-driven vibronic dynamics in the neural and the transient cation [71]; the coordinate dependence of the photorecombination amplitudes [43,65,67]. For wavelengths beyond ≈ 1.6µm or intensities above 300 TW cm −2 , the latter contributions may dominate the PACER spectrum (see also Ref. [65,67]), even for a single-cycle driving pulse.
Furthermore, we demonstrate that the PACER concept breaks down for harmonic emission close to resonances (constructive or destructive) in the recombination matrix elements. At these energies, no simple relationship exists between the harmonic photon energy and the emission time. The isotopic ratio is then no longer representative of the intrinsic vibronic dynamics in the cation. For any given combination of the wavelength and inten-sity of the driving field, the PACER spectrogram resulting from a resonance may be indistinguishable from the signal due to the intrinsic vibronic dynamics. It is therefore essential to perform PACER experiments at a number of wavelengths and/or intensities. Features which appear at a fixed, or nearly fixed emission energy are likely to originate from a resonance in the recombination matrix elements.
Although the present analysis focuses on one of the simplest polyatomic molecules, the methane, both conical intersections and resonant features in photorecombination matrix elements are ubiquitous in polyatomic molecules. High-harmonics spectroscopy thus provides a tool for exploring both, in the regime not easily accessible with other techniques.
The effects discussed presently are of the sub-cycle nature, and are already present for a single-cycle driving pulse. For longer pulses, additional physical mechanisms will come into play [64,68,71,77,79]. A careful examination of these multi-cycle effects and their interaction with the intrinsic sub-cycle dynamics is one of the possible future applications of the MC-SFA-GWP approach.
Appendix A: Evaluation of the dt1 integral in Eq. 26 We need to evaluate an integral in the form: dt 1 f k s (t 1 ) e −iφ d ( ks(t1),t,t1) , where k s is in turn a function of t 1 (Eq. 27) and φ d is defined by Eq. 24. We assume that f k s depends at most linearly on k s . Expanding φ d through the third [98] order in t 1 around t s , changing integration variable to τ = t 1 − t s , and extending the integration limits to infinity, we obtain: In Eq. A2, f , φ d , their derivatives, and F are evaluated at k s and t s pairs solving Eq. 28. For linearly-polarized driving field, ∂ 2 φ d ∂t 2 s (Eq. A6) and all but the last term in ∂ 3 φ d ∂t 3 s (Eq. A3) vanish. For low-frequency fields, these terms remain negligible for moderate non-zero ellipticities as well (high harmonic signal vanishes for large ellipticities [98]), so that: We now note that (Eq. 10.4.32 of Ref. [142]): where Ai is an Airy function. (Eq. A8 is obtained by differentiating Eq. A7 with respect to x.) Therefore, Eq. A2 becomes: I ≈ e −iφ d 2π f 2m where again all quantities are evaluated at pairs k s , t s satisfying Eq. 28. If multiple roots are present, summation over all roots is implied.
Appendix B: Derivation of Eq. 30 The presence of the vibronic state energy expectations E bm and E a n , rather than the manifold-average I p , in Eq. 31 is intuitively appealing, but requires some additional justification. Application of Eq. A9 to the integral of Eq. 26 leads to: Υ bma n = F (t s ) · R bma n k s 2π 2m e 2 2 F 2 (t s ) 1/3
If the HamiltoniansĤ I andĤ 0 are diagonally-dominant in the basis of corresponding vibronic product states, Eqs. B2-B3 reduce to: where E bm and E an are given by Eqs. 32-33. Substituting Eqs. B4-B5 into Eq. B1, we note that the first three terms are in fact the lowest-order contributions to the Taylor expansion of the first term in Eq. 30 for δE = (E bm − E a n ) − I p . Similar expansion for the second term corresponds to a term quadratic in τ , neglected in deriving Eq. A9. To within the accuracy expected from Eq. 30, we can therefore replace I p by (E bm − E a n ) in the last term of Eq. B1 as well, giving Eq. 31. The exponential part of Eq. 31 coincides with the result of the weak-field asymptotic tunneling theory (WFAT) [143].