Time-resolving the ultrafast H2 roaming chemistry and H3+ formation using extreme-ultraviolet pulses

The time scales and formation mechanisms of tri-hydrogen cation products in organic molecule ionization processes are poorly understood, despite their cardinal role in the chemistry of the interstellar medium and in other chemical systems. Using an ultrafast extreme-ultraviolet pump and time-resolved near-IR probe, combined with high-level ab initio molecular dynamics calculations, here we report unambiguously that H3+ formation in double-ionization of methanol occurs on a sub 100 fs time scale, settling previous conflicting findings of strong-field Coulomb explosion experiments. Our combined experimental–computational studies suggest that ultrafast competition, between proton-transfer and long-range electron-transfer processes, determines whether the roaming neutral H2 dynamics on the dication result in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm{H}}_3^ +$$\end{document}H3+ or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm{H}}_2^ +$$\end{document}H2+ fragments respectively.

T rihydrogen (H þ 3 ) is one of the most abundant ions in the universe with an important role in the formation of complex molecules in the interstellar medium (ISM) 1 . The formation of H 3 + and particularly its destruction in the ISM attracts a vast experimental and theoretical effort [1][2][3][4][5][6][7][8][9][10] . While H 3 + is typically formed by H 2 + + H 2 collisions 11 , it is a ubiquitous product of organic-molecule ionization processes such as electron impact 12 , fast ion bombardment 13 , multi-photon strong-field laser ionization 14,15 , as well as single extreme-ultraviolet (EUV) photon ionization 16 . However, the mechanisms for H þ 3 formation, as well as of the related H þ 2 product, are still poorly understood. Early fragment imaging experiments using strong-field laser ionization by Yamanouchi and co-workers, provided first evidence as to the time scale of the H þ 3 formation in ionization of methanol 14,15,17 . In contrast to products such as H þ and H þ 2 , in which the ejection directions exhibit a typical strong correlation with the strong-field laser polarization, H þ 3 ejection from ionization of methanol was observed to be isotropic, suggesting rotational depolarization due to a "much greater" than 1.4 ps lifetime of the parent methanol dication 17 . A particularly longtime scale is in accord with the significant structural rearrangement required for H 3 + formation from the parent methanol molecule. Early theoretical studies proposed an intricate mechanism involving formation of a neutral H 2 that roams on the doubly-ionized ground state potential surface of the organic molecule, culminating in proton abstraction and formation of H þ 3 18,19 The analogous roaming H-atom dynamics are known to exhibit prolonged roaming lifetimes ranging from hundreds of femtoseconds to thousands of picoseconds 20,21 .
Interpretation of multi-photon strong-field laser ionization experiments poses intrinsic challenges due to co-existence of direct and indirect ionization mechanisms. This intrinsic ambiguity has previously been the source of surprising and controversial results-for example concerning the double-proton transfer in a DNA base-pair model, which assignment as a sequential process was challenged due to the co-existing direct and indirect strong-field ionization mechanisms 22,23 . Furthermore, in the case of methanol, Itakura et al. demonstrated the non-linear sensitivity of the dication dissociation spectra and branching ratios to the exact strong-field laser parameters 24,25 . Nevertheless, in a recent pump-probe study implementing consecutive strong-field (~10 14 W/cm 2 ) laser pulses, Ekanayake et al. reported H þ 3 product signal appearance on an ultrafast (~100 fs) time scale 15 . However, in another recent pump-probe study, also implementing strong-field lasers with similar peak intensities albeit with shorter pulse lengths, Yamanouchi and co-workers were able to observe a~38 fs beating in the time resolved H þ 3 product signal, which time scale they assigned not to its formation dynamics on the final dication but rather to the vibrations of a singly-ionized intermediate 14 . Evidently, strong-field doubleionization can proceed via different competing direct as well as indirect mechanisms, involving intermediate state dynamics that may mask the time evolution of the making and breaking of chemical bonds 22,23 .
In addition to these experimental challenges, strong-field laser ionization is also very difficult to follow using computational simulations. Most theoretical modeling efforts do not address the strong fields directly. Recent ab initio molecular dynamics (AIMD) simulations on the methanol dication ground state showed~4% probability for the dominant H þ 3 formation channel 26 . However, those ground-state dynamics simulations could not account for the experimentally measured branching ratios or predict many of the experimentally observed products such as the C-O bond cleavage or proton migration to form H 2 O þ15, 26 . As we discuss below, these appear only when excited dicationic states and inclusion of non-adiabatic dynamics is allowed.
The purpose of this paper is to develop a combined experimental and theoretical approach that time-resolves the roaming H 2 chemistry responsible for the H 3 + formation. Taking advantage of high-order harmonic generation (HHG) of ultrafast EUV pulses 16,27,28 , we implement a single-photon double-ionization pump, followed by a time-delayed near-IR (nIR) probe. In this way the new approach not only bypasses the underlying uncertainty in strong-field laser experiments, but is also suitable for analysis using ab initio theoretical modeling 16 . We use XMS-CASPT2/(8e,8o)/ aug-cc-pVDZ/potential surfaces and nonadiabatic coupling terms calculated by the BAGEL code 29 , interfaced with a modified version Newton-X (v1.4.0) 30,31 . In our previous work we demonstrated that this level of theory, simulated for the singlet states manifold, allows for a detailed and multifaceted agreement with the Coulomb explosion (CE) product branching ratios, channel-specific kinetic energy release (KER) and three-body momentum correlation spectra, measured in time-independent experiments 28 . Here, the theoretical simulations have paved the way for the experimental timing of the roaming H 2 dynamics that forms H 3 + on the low lying states of the dication by excitation to higher lying states, which quench H 3 + formation and enhance three-body fragmentation. The combined pump-probe measurements and the nonadiabatic AIMD trajectory simulations allows us to unambiguously conclude that H þ 3 formation in double-ionization of methanol occurs on an ultrafast sub 100 fs time scale, in direct competition with the H 2 + formation channel, settling previous conflicting experimental findings 15,17 .

Results and discussion
Methanol dication dynamics. Figure 1 shows the potentials of the low-and the high-lying energy states of the methanol dication as a function of the CO stretch. The dynamics in the latter case typically results in rapid dissociation, accompanied by cleavage of the C-O bond as well as a three-body breakup 16 . In contrast, the dication formed in the low-lying states initially exhibits a potential charge imbalance that hinders direct CE due to the~3 eV potential barrier indicated in Fig. 1. As a consequence dynamics on the low-lying states are prolonged and highly complex, involving the emergence of a roaming neutral H 2 . Figure 2 compares the simulated branching ratios for the H þ 3 product of roaming H 2 , with the combined branching ratio for the three-body breakup channels, calculated for each initial Fig. 1 The adiabatic potential curves for the C-O bond breaking channel of CH 3 OH 2+ using MS-CASPT2/(12e,10o)/aug-ccpVTZ 16 . As the C-O strech coordinate is extended away from the Frank-Condon geometry of the neutral methanol ground state, shown by the black curve, the dication potentials exhibit a~3 eV barrier that prevents the C-O bond breaking on the low lying states. The potentials are calculated while keeping all other coordinates at their neutral methanol ground state value.
excitation of the methanol dication. Over 1/3 of the ground-state trajectories produce H þ 3 , this compared with~4% obtained in previous ground-state simulations that did not include the second order perturbation theory corrections 26 . The high H þ 3 formation probability on the ground state, drops for the higher lying states and is completely quenched once the CO bond cleavage becomes possible above the third excited state. In contrast, the three-body breakup exhibits an opposite trend, which increases for higherlying excited states.
These theoretical predictions provide a handle for the experimental time-resolved probing of the dynamics using a delayed nIR pulse, following excitation with the ultrafast EUV pulse. Where the time delayed probe will excite the transient dication to higher-lying states, consequently quenching H þ 3 formation and enhancing three-body breakup. However, once the excess internal energy is released in a successful CE, the product branching is expected to be less affected by the probe pulse.
In designing the probe pulse, we ensure that its peak intensity is kept well below the threshold for strong-field CE, such that at long negative time delays the branching ratios are identical to the ratios measured with the EUV pulse alone. In particular, the H þ 3 formation branching ratio is 6% and the three:two body ratio is 3 to 1. The effect of the nIR probe pulse delay (with respect to the EUV pulse) on the relative enhancement of the three:two body ratio is shown in Fig. 3a. The three:two body ratio increases by up to~8% at positive time delays, as the nIR probe arrives shortly after the dication formation. This effect decays as the probe pulse arrives and positive time delays longer than~70 fs. For comparison, Fig. 3b shows the enhancement of doubly-ionized Ne 2þ yield as a function of the nIR probe delay that reflects the instrumental response time. The full line in Fig. 3b represents a fit of the Ne 2þ yield, assuming photoionization of high lying Ne þ* cations by the time delayed nIR pulse, which rise time reflects the cross-correlation of the EUV and nIR pulses. The dashed red line represents the corresponding Gaussian cross-correlation function, in agreement with the <35 fs FWHM of our laser pulses. Figure 3c shows the time correlated relative change in the H þ 3 þ COH þ branching ratio, which exhibits upto~12% suppression. The full lines in Figs. 3b and 3c show a model trace including an exponential~70 ± 25 fs lifetime, convoluted with the instrumental time response directly determined based on the Ne 2þ data shown in Fig. 3b. While the three:two body ratio appears to return to its unperturbed value, the asymptotic H þ 3 formation remains suppressed by 2.5% even at long time delays. It should be mentioned that while the energy needed to dissociate the H þ 3 ground state is~4.5 eV 32 , nIR photodissociation of the highly vibrationally excited H þ 3 can still be expected 33,34 . We therefore assign the residual H 3 + depletion at long times to photodissociation of the vibrationally hot H þ 3 cations after the CE is completed. To provide additional insight into the nIR probe mechanism we compare, in Fig. 4, the shapes of normalized KER distributions, collected at different pump-probe delays. The red line shows the KER measured for negative times, where the nIR probe arrives over 150 fs before the EUV pulse. Similar to the branching ratios, the KER spectrum at negative times is identical to the one measured with the EUV pulse only. The blue line shows the KER at positive delays longer than 150 fs, while the open circles show the KER during the transient suppression times. Thus, while the branching ratios reveal clear time dependence, the KER spectra are not significantly affected by the nIR pulse. We therefore conclude that the field of the nIR probe acts as a switch between H 3 + formation and three-body breakup but does not significantly change transient dynamics leading to the specific channel, as suggested for other experiments implementing EUV pump and a  strong-field nIR probe 27,35 . Photo-excitation of the transient timeevolving dication before the CE is expected to promote also competing C-O bond breaking channels and enhances fragmentation as predicted by the AIMD simulations. Interestingly, the C-O bond breaking branching ratio does not exhibit a clear time evolution within the experimental error bars (not shown), possibly due to its competing enhancement on the higher lying states and suppression by three-body fragmentation 28 . For the less abundant CE channels of H + + CH 3 O + , H 2 + + CH 2 O + and H 2 O + + CH 2 + , statistical errors limit the determination of their individual time evolving branching ratios.
Further dynamical insight can be obtained from analysis of the AIMD trajectories culminating in the formation of H þ 3 . Figure 5a shows the time evolving inter-fragment velocity, corresponding to the time derivative of the distance between the H þ 3 and COH þ products. Before dissociation occurs, the inter-fragment velocity exhibits oscillations between positive and negative values, reflecting the roaming H 2 motion away from and towards the HCOH 2+ , prior to the formation of H þ 3 as can be seen in the typical AIMD movies provided (see Supplementary Movie 1). Interestingly, Palaudoux et al proposed a possible concerted mechanism for H þ 3 ejection from a CH 3 Cl 2þ dication 8 . However, none of the trajectories simulated here could be attributed to a concerted ejection of H þ 3 . The arrow in Fig. 5a indicates the dissociation time assigned to the highlighted black trajectory, where for each trajectory the dissociation time is defined as the last point of attraction between two dissociating fragments, after which their relative velocity is monotonically increasing until reaching the asymptotic KER. The bars in Fig. 3d show the histogram of the total of 66 trajectories resulting in H þ 3 þ COH þ dissociation, peaking at~100 fs. While explicit theoretical modeling of the nIR probe is beyond the scope of this paper, the experimental time resolved branching ratios can be compared with the simulated suppression time window that tentatively extends from the formation of each simulated dication until its dissociation time. The full line in Fig. 3d shows the average simulated suppression, constructed from the 66 trajectories that form H þ 3 and convoluted with the experimental instrumental response, in good agreement with the transient branching ratio measurements and an ultrafast sub 100 fs lifetime.
Previous simulations also showed similar roaming H 2 dynamics within the dication and ultrafast H þ 3 formation 18,26 . However, to understand the ultrafast lifetime of the roaming H 2 it is important to consider also its other decay channels. Earlier AIMD simulations using CISD and CASSCF electronic potentials reported unbalanced charge dissociation of the methanol dication ground state to form H 2 þ CHOH 2þ , with over 11% and 18% branching ratio respectively 18,26 . In contrast, the non-adiabatic AIMD simulations using CASPT2 potentials indicate that the neutral H 2 cannot escape, it is polarized and bound by the CHOH 2þ dication. Like the charge-transfer barrier preventing C-O bond cleavage with unbalanced charge, clearly visible in Fig. 1, the neutral H 2 cannot dissociate before charge is transferred and CE can proceed. In the H þ 3 formation mechanism the abstraction of a third proton from either methyl site or hydroxyl site allows dissociation by the transfer of a proton. This mechanism is in direct competition with the transfer of an electron that results in a H þ 2 þ CHOH þ breakup. Figure 5b represents typical AIMD trajectories evolving towards H þ 2 formation, showing the interfragment velocity between the H þ 2 and CHOH þ products. The highlighted trajectory exhibits the "inverse harpooning" mechanism 28 , observed in all the trajectories resulting in the H þ 2 þ CHOH þ breakup. The star labeled arrow indicates the time at which a neutral H 2 molecule begins to separate from the CHOH 2þ dication on the highlighted black trajectory. Although the relative velocity continues to be positive until the asymptotic dissociation limit is reached, the neutral H 2 is still bound by the CHOH 2þ dication. This is evident from the deceleration of the relative velocity. Although reaching long inter-fragment distances, as high as 9 angstroms, no neutral H 2 escape. Eventually, at a time indicated by the second arrow, a long-range adiabatic electron transfer from the neutral H 2 to the CHOH 2þ ignites a CE, producing H þ 2 . The charge transfer is evident both from the computed Mulliken charges as well as from the sudden transition to a long acceleration to a high asymptotic KER, typical of the long-range Coulomb repulsion. The open and full bars in Fig. 3e show histograms of the times assigned to the neutral H 2 separation and the inverse harpooning times respectively. Interestingly, the competing proton and electron transfer mechanisms that facilitate the release of the two molecular hydrogen ions proceed on comparable ultrafast time scales of~100 fs, in agreement with the measured H þ 3 suppression time window on the transient dication. This competition can be directly visualized in a selected trajectory Fig. 4 Measured H3 + + COH + KER spectra. Comparing spectral shapes for three t nIR windows: t nIR < −150fs, −60 fs < t nIR < 60 fs and 150 < t nIR , shown by the solid red line, black circles and blue solid lines respectively. (see Supplementary Movie 1), in which the roaming neutral H 2 significantly separates from the dication, typical of "inverse harpooning". Nevertheless, it is pulled back to abstract a proton from the hydroxyl and form H þ 3 . By using a low-field ultrafast EUV pulse to doubly ionize methanol and a time delayed nIR probe pulse it is possible to remove the intrinsic uncertainties of strong-field laser experiments 14,15,26 . A~70 fs lifetime is observed for the transient suppression of H þ 3 þ COH þ branching ratio, accompanied with a correlated enhancement of the three:two body fragmentation ratio. Furthermore, non-adiabatic AIMD simulations on the ground and excited CASPT2 potentials provide a detailed picture of the ultrafast roaming H 2 dynamics, culminating in the trihydrogen product on an ultrafast~100 fs time scale, in agreement with the experimental result. In contrast to earlier simulations that did not include second-order perturbation theory corrections and predict significant ejection of neutral H 2 26 , the simulated roaming neutral H 2 was found to exhibit ultrafast competition between proton abstraction resulting in H þ 3 and the "inverse harpooning", a longrange electron transfer that results in the H þ 2 þ CHOH þ CE channel. Further experimental and theoretical work on deuterated methanol as well as other organic systems will allow to provide a more detailed understanding of the different pathways for H þ 3 formation, explored so far only by strong-field laser experiments. Such ultrafast roaming H 2 chemistry, accompanied by competing proton and electron transfer dynamics described here are expected to occur also in other ionized systems produced by ionizing radiation damage, e.g., by cosmic radiation in planetary and interstellar environments or manmade light sources for single molecule crystallography experiments 19,[36][37][38] .

Methods
Experimental. The single-photon CE imaging setup has been described earlier 16,27,28 . Briefly, in our experiments, a~7 mJ, <35 fs, 803 nm near IR (nIR) pulses from a Solstice 39 laser are split into a 2.1 mJ pump and 4.9 mJ probe pulses. The pump pulse is focused in a semi-infinite neon gas cell, for HHG of ultrafast EUV pulses at a 1 kHz repetition rate 40,41 . Where the EUV is spatially filtered from the higher divergence nIR fundamental 42 . The remaining~4.9 mJ nIR pulse is time delayed to probe the EUV initiated dynamics and is merged with the ultrafast EUV pulse at a small~1 degree angle at the center of a home built 3D coincidence imaging spectrometer 16,43 , where both beams cross a skimmed effusive beam of CH 3 OH. The nIR beam is mildly focused behind the spectrometer with a 610 mm lens,~203 mm behind the spectrometer, such that at the nIR beam is inside the EUV, as validated by using the spectrometer to image the parent ion birth positions. The cationic products are accelerated from the interaction volume towards a time and position sensitive detector, allowing 3D coincidence imaging of the ion recoil velocities [43][44][45] . Low count rate and center of mass momentum conservation are used to suppress random coincidence cation signal due to dissociative ionization of two different parent molecules. Any residual contributions from random coincidence are estimated and subtracted based on the measured single cation event probabilities. The time resolved Ne 2+ signal shown in Fig. 3b was used to characterize the temporal overlap of the EUV and nIR pulses, while singly ionized methanol signal was used to monitor the temporal overlap during the long acquisition times. Both 3-body enhancement and H þ 3 depletion were fitted together using the Ae À t τ 1 þ erf t ffiffi 2 p σ À σ ffiffi 2 p τ functional form 46 . Where the crosscorrelation width σ was fixed according to the Ne 2+ rise time analysis. The free parameters were the respective enhancement\ depletion amplitudes and a common lifetime τ for both processes. In addition to the transient depletion of H þ 3 , a second term with an infinitely long lifetime was added to describe the residual H þ 3 depletion at long time delays.
Theoretical. The initial configurations (geometries and velocities) of ground state methanol were sampled from a 300 K AIMD simulation of neutral methanol, calculated at the CASSCF level using the MOLCAS package 47 at the (14e,10o)/aug-ccpVTZ active space/basis-set level. Each of the 100 configurations was used for initiating non-adiabatic molecular dynamics calculations on each one of the seven lowest-lying electronic excited states of the dication. This generated 700 trajectories altogether. The ensuing non-adiabatic dynamics were approximated using surfacehopping molecular dynamics trajectories 48 generated at the XMSCASPT2/(8e,8o)/ aug-cc-pVDZ/density-fitting level using the BAGEL electronic structure package 29 within the so-called "SS-SR" contraction scheme 49 used for internally contracted basis functions in CASPT2, where a vertical shift was set to 0:2E h . We checked the basis set convergence by comparing aug-cc-pVDZ and aug-cc-pVTZ results for three single-point calculations: (a) double ionization energy at the Franck Condon point; (b) the excitation energy between the dication ground state and its excited states. (c) the barrier height on the dication ground-state. We found that basis set choice led to a very small relative difference of 1% in all cases. For example, the barrier height on the dication state is nearly 3 eV and the difference between the values given by the two basis sets is 0.02 eV. The BAGEL code was interfaced with a modified version Newton-X (v1.4.0) program 31 for carrying out the surface hopping dynamics 30 . To facilitate the trajectory calculations, the system in adiabatic state n is allowed to hop only to the state m nearest in energy above or below it (i.e. we neglect the nonadiabatic coupling terms τ nm unless n À m j j¼ 1). We modified the way Newton-X interfaces with the BAGEL code to enable this approximation. The time step for the NA-AIMD trajectories is 0.3 fs. The ab initio dynamics are typically followed until 300 fs or until the inter-fragment velocities are observed to reach an asymptotic monotonic behavior. At this stage, the effect of the residual long-range Coulomb repulsion on the final velocities is taken into account using the classical equations of motion applied to the center of masses of the cationic fragments. Once the asymptotic fragment identity is determined, the inter-fragment velocities are calculated as the time derivative of the distance between the fragment centers of mass.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Code availability
Simulations were performed using open source programs. Modified code files are available from the corresponding authors upon reasonable request.