Monte Carlo wave-packet approach to trace nuclear dynamics in molecular excited states by XUV-pump-IR-probe spectroscopy

Recent research interests have been raised in uncovering and controlling ultrafast dynamics in excited neutral molecules. In this work we generalize the Monte Carlo wave packet (MCWP) approach to XUV-pump–IR-probe schemes to simulate the process of dissociative double ionization of H 2 where singly excited states in H 2 are involved. The XUV pulse is chosen to resonantly excite the initial ground state of H 2 to the lowest excited electronic state of 1 (cid:2) u + symmetry in H 2 within the Franck-Condon region. The delayed intense IR pulse couples the excited states of 1 (cid:2) u + symmetry with the nearby excited states of 1 (cid:2) g + symmetry. It also induces the ﬁrst ionization from H 2 to H 2 + and the second ionization from H 2 + to H + + H + . To reduce the computational costs in the MCWP approach, a sampling method is proposed to determine in time the dominant ionization events from H 2 to H 2 + . By conducting a trajectory analysis, which is a unique possibility within the MCWP approach, the origins of the characteristic features in the nuclear kinetic energy release spectra are identiﬁed for delays ranging from 0 to 140 fs and the nuclear dynamics in the singly excited states in H 2 is mapped out.


I. INTRODUCTION
For molecular systems, due to interactions between nuclei and electrons, the electronic dynamics directly induced by the applied laser pulses is usually accompanied by nuclear dynamics. For example, dissociation of molecules can be induced by electronic excitation or ionization [1,2]. Many pioneering works have exploited ultrashort laser pulses to reveal and control the breakup dynamics including ionization and dissociation dynamics in molecules, where the applied laser pulses can be either single laser pulses or two laser pulses in a pump-probe configuration. For instance, 40-fs laser pulses with wavelengths tunable from 800 to 1850 nm were employed to probe the nuclear dynamics in H 2 + by recollision excitation after tunneling ionization from H 2 [3,4]. Under the same recollision excitation mechanism, carrier-envelope stabilized fewcycle near-IR laser pulses were exploited to achieve electron localization, which is a crucial precondition in realizing control of chemical reaction dynamics, in molecular dissociation of H 2 + [5]. More efficient control over electron localization in the same system was later realized in several dissociation channels involving recollision excitation and bond-softening mechanisms by using few-cycle [6] or 50-fs two-color [7] mid-IR laser pulses. In addition, by using IR laser pulses with varying wavelengths, intensities, or pulse durations [8][9][10], nuclear dynamics in H 2 + were resolved by studying the nuclear kinetic energy release (KER) spectra following double ionization, i.e., full breakup of H 2 after removal of two electrons.
Pump-probe spectroscopy is powerful in uncovering and controlling ultrafast dynamics taking place in molecules, especially since the emergence of attosecond laser sources. In general, the pump pulse is used to initiate molecular dynamics and the probe pulse is used to take snapshots of the dynamics in real time as a function of the controlled delay. Recently, there has been a considerable amount of work devoted to revealing and controlling ultrafast dynamics in molecules by using XUV-pump-IR-probe spectroscopy. For instance, control over electron localization in dissociation of H 2 + was achieved by changing the delay between an attosecond pulse train (APT) and a many-cycle IR laser pulse [11,12]. Instead of applying an APT, exploiting a single isolated attosecond pulse as the pump pulse was also able to realize control over electron localization in the same H 2 system [13]. Apart from controlling electron localization, this XUV-IR scenario is especially useful in monitoring the nuclear dynamics in molecular excited states. For example, this scenario was applied to characterize the singly excited molecular dynamics by using the ionization signal [14][15][16].
For most of the works in the XUV-IR scenario (e.g., Refs. [14][15][16]) only the single-ionization channel is open, which limits the information about the ultrafast dynamics obtained by the measurements. More information can be obtained if double-or multiple-ionization channels are open. For example, former studies of the double-ionization channel of H 2 due to interaction with intense IR laser pulses have illustrated that the purely repulsive character of the Coulomb potential can provide an unequivocal identification of the different dissociative ionization pathways [17]. Up to now, there have been few theoretical and experimental works applying the XUV-pump-IR-probe spectroscopy to study the dynamics involved during the double-ionization processes in molecules. For the theoretical part, the difficulty of treating molecular double ionization lies in solving the time-dependent Schrödinger equation (TDSE) to describe both nuclear and electronic dynamics in the presence of two fields and nonperturbative dynamics. Even for the simplest molecules, this task is beyond the capabilities of conventional computers due to the involvement of a huge number of continuum states in the process. As a result, it is very attractive to develop a theory that simplifies the computations while still being accurate. Moreover, ideally, such a theory should facilitate an identification of the main physical mechanisms. To this end, one can resort to the Monte Carlo wave packet (MCWP) approach, which, in the context of strong-field physics, was first presented by Leth et al. [10] for the description of nuclear motion in the presence of dissipative electron dynamics. Since then, this methodology has been validated in a number of cases. For example, a detailed account of the theoretical foundation of the MCWP approach and its limitations including its relation to the master equation of the Lindblad form for open quantum systems was presented in Ref. [18] (see also the original works connected to quantum optics problems [19][20][21]). Later, a very detailed comparison with experimental nuclear KER spectra following interaction with near-infrared pulses was carried out for different wavelengths including pump-probe schemes [22]. The agreement between experiment and theory was excellent. Apart from reproducing experimental results for near infrared wavelengths, this approach was also capable of performing an accurate description of the dissociative multiple-ionization process in the XUV regime [23].
In this work we resolve nuclear dynamics in the singly excited states in H 2 by generalizing the MCWP approach to simulate double ionization of H 2 within XUV-pump-IR-probe schemes. The process of double ionization considered in this work is as follows (Fig. 1). Starting from the ground state of H 2 , i.e., from the X 1 g + (v = 0) state, the nuclear wave packets in the excited electronic states of 1 u + symmetry are first efficiently launched via a one-photon absorption process caused by the interaction with an ultrashort attosecond XUV pulse. Then a delayed intense femtosecond IR laser pulse strongly couples these excited states of 1 u + symmetry to the nearby excited states of 1 g + symmetry. Thus, there are nuclear wave packets evolving along the potential energy curves for the excited states of both 1 g + and 1 u + symmetries (Fig. 1). At the same time, substantial ionization from H 2 to H 2 + , induced by the intense IR laser pulse, occurs from these excited states due to their small ionization potentials. After the first ionization, the remaining part of the IR laser pulse couples the electronic states in H 2 + , e.g, the 1sσ g and 2pσ u states, and simultaneously the second ionization from H 2 + to H + + H + occurs. All these dynamics, including the ionization processes, are modeled by the present approach. Our theory treats the double ionization in a consistent theoretical framework for this attosecond pumpprobe KER spectroscopy. By studying the nuclear KER spectra ( 1 u + ) symmetry in H 2 below the 2 g + (1sσ g ) curve in H 2 + . The two lowest potential energy curves in H 2 + , i.e., the 1sσ g and 2pσ u curves, and the 1/R Coulombic curve of H + + H + are also shown. The thin (purple) vertical arrow indicates the interaction between the XUV pump pulse and the ground-state molecule and the thick (red) arrows represent ionization from H 2 + to H + + H + induced by the IR probe pulse. There is ionization from H 2 to H 2 + induced by the IR field in spite of no arrow representing this process in the figure. of the repulsive protons in the double-ionization channel for different pump-probe delays, the nuclear dynamics in the excited states of H 2 can be mapped out. For example, characteristic features in the nuclear KER distribution as a function of the pump-probe delay are analyzed by conducting a trajectory analysis, where the trajectories that contribute the most to the observed features within our MCWP sampling are identified.
This work presents the consideration of utilizing the nuclear KER spectra following double ionization of the prototype H 2 to recover the nuclear dynamics in the singly excited states by using XUV-pump-IR-probe spectroscopy. We hope our findings can stimulate the experimental studies in this scenario. This paper is organized as follows. In Sec. II we describe the application of the MCWP approach and its validation within XUV-pump-IR-probe schemes. In Sec. III the nuclear KER distributions following double ionization of H 2 are obtained for three different pump-pulse durations. Through analyzing the characteristic features in the nuclear KER distributions, the nuclear dynamics in the excited states of H 2 is retrieved. Section IV summarizes. A detailed description of the sampling method is reported in Appendix A and a discussion of nuclear wave packets in H 2 + from B3 1 u + is given in Appendix B. Atomic units are used unless otherwise stated. 043426-2

A. The MCWP approach
The MCWP approach has been thoroughly described in Refs. [10,18,[22][23][24][25][26]. In this work we extend it to obtain the nuclear KER spectrum of the protons following double ionization of H 2 in the two-color XUV-pump-IR-probe scheme where singly excited states in H 2 are involved. When we use the MCWP approach to describe double ionization of H 2 , three systems are involved, i.e., the neutral H 2 , H 2 + , and H 2 2+ . For a specific system containing electrons, due to its interaction with the applied external field, ionization may occur which takes the specific system to the system with one electron less. Meanwhile, the electron removed from the system is now in the continuum states. The MCWP approach regards the integrated influence of the continuum states as the surroundings of a given system. Thus, the unidirectional loss of electrons from a given system to its surroundings can be treated by adding a non-Hermitian term to the system Hamiltonian. For a given system, the Hamiltonian is The equivalence of the constructed non-Hermitian Hamiltonian (1) to the Lindblad form master equation for an open system is shown in Ref. [18]. In Eq. (1), C m is a jump operator, which specifies the transition from the electronic state |φ el R,m in a given charge state to the electronic states |φ el R,n in the system with one electron less [10,18], i.e., where | R is the position eigenket of the nuclear coordinate and m ( R,t) is the ionization rate for the molecule when the electron (electrons) is (are) in the |φ el R,m state. The term n M n ( R,t)|φ el R,n in Eq. (2) represents the fielddressed ground state in the new charge state, which is a linear combination of the field-free electronic states |φ el R,n . For example, following a quantum jump from H 2 , the expression for M n ( R,t) in H 2 + can be found in Ref. [27]. The above-mentioned electronic states, which have a parametric dependence on the internuclear distance R, can be obtained by solving the time-independent Schrödinger equation of the electronic Hamiltonian H el at fixed R, i.e., H el |φ el R,m = E el m ( R)|φ el R,m . Since the influence of the continuum states is integrated in the non-Hermitian term in Eq. (1), the electronic states involved in the jump operator in Eq. (2) are bound states. The electronic Hamiltonian H el is related to the system Hamiltonian H s through H s = T N + H el + V L (t), where T N is the nuclear kinetic energy operator and V L (t) = − μ · F (t) is the light-molecule interaction operator. In this operator, μ is dipole operator of the electrons and F (t) = F XUV (t) + F IR (t) is the combination of the electric fields of the applied pump and probe pulses. In this work, both pulses are linearly polarized parallel to the molecular axis. For the XUV pulse, its vector potential is described by a sin 2 envelope, i.e., where T XUV is the XUV pulse duration and ω XUV is the central frequency of the XUV pulse. The instantaneous field strength of the XUV pulse is obtained from the relation F XUV (t) = −dA XUV (t)/dt, where A XUV0 ω XUV is its peak field strength.
For the time-delayed femtosecond IR laser pulse, the electric field profile is Gaussian, i.e., where F IR0 is the peak field strength, τ the time delay between the two pulses, ω IR the central frequency of the IR field, and T IR the FWHM duration of the temporal intensity profile of the IR laser pulse. In this work, both the first and second ionization events of H 2 are overwhelmingly induced by the IR laser pulse, i.e., the XUV pulse is relatively weak with its central photon energy well below the first ionization threshold of H 2 . For dissociative double ionization of H 2 , the three involved charge states correspond to H 2 , H 2 + , and H + + H + , respectively. For a specific charge state, the total wave function | (t) satisfies the TDSE i|˙ (t) = H (t)| (t) , with H (t) from Eq. (1) and | (t) expressed using the ansatz where χ m ( R,t) is the nuclear wave function corresponding to the electronic state |φ el R,m . Similar to the jump operator in Eq. (2), the total wave function in Eq. (5) is also a coherent superposition of bound electronic states. Substituting Eq. (5) into the TDSE and projecting on the electronic eigenstate φ el R,m |, we obtain a system of coupled equations for the nuclear wave functions where we have neglected nonadiabatic couplings. The last term in Eq. (6) describes the influence from the dipole couplings between the |φ el R,m state and the other states within a given system, e.g., in either H 2 or H 2 + . Equation (6) shows that the nuclear wave packet χ m ( R,t) evolves along its electronic potential energy surface E el m ( R) accompanied by coherent coupling with the other bound electronic states in the same charge state and ionization loss to the adjacent charge state.
We assume that the molecule is aligned along the linear polarization direction of F (t). This assumption reduces the three-dimensional coupled equations to one-dimensional equations and the equations of motion within a given charge state  can be written in matrix-vector form as · · · · · · · · · · · · · · · · · · − d 2 dR 2 + E el m (R) − i 2 m (R,t) · · · V mk (R,t) · · · · · · · · · · · · · · · · · · · · · V km (R,t) where each χ m (R,t) is expanded in a primitive basis for the R degree of freedom. To obtain the complete nuclear dynamics for double ionization of H 2 , we solve Eq. (7) in three charge states corresponding to H 2 , H 2 + , and H + + H + , respectively. It is the quantum jump that transfers nuclear motion in a given charge state to nuclear motion in the charge state with one electron less. Thus, in our model, we have the first jump from H 2 to H 2 + following the first ionization and the second jump from H 2 + to H + + H + following the second ionization. In Fig. 1 the relevant Born-Oppenheimer potential energy curves are plotted. The energies corresponding to the ground and excited states of the neutral molecule are taken from accurate correlated ab initio calculations [28,29]. The cation energies are obtained from exact simulations using a one-electron diatomic molecule code [30]. The figure shows the potential energy curves for the six lowest electronic states of 1 g + and 1 u + symmetries in H 2 , the potential energy curves for the two lowest electronic states in H 2 + , and the 1/R Coulombic curve of H + + H + . The dipole moment functions between electronic states entering the couplings V mk (R,t) have been obtained through a configuration-interaction calculation for H 2 by using a procedure similar to that described in Ref. [31] and taken from the literature for H 2 + [32]. To make Eq. (7) solvable, it is necessary to have an estimate of the state-resolved ionization rates m (R,t). In H 2 , the R-dependent ionization potentials I m (R) of the considered singly excited states are much smaller than the ionization potential of the ground state ( Fig. 1). Thus, in order to ensure ionization from the excited potential energy curves to work in the tunneling regime, the central wavelength of the IR laser pulse is chosen to be long enough to satisfy the corresponding ionization conditions, i.e., the IR photon energy ω IR should be smaller than the ionization potential for the highest-lying excited state in H 2 in Fig. 1 and the Keldysh parameter γ for this state should be smaller than 1 [33]. Apart from numerically solving the TDSE [34], methods such as the molecular Ammosov-Delone-Krainov theory [35,36] and the weak-field asymptotic theory [37][38][39][40][41] have been developed to obtain analytical formulas for the tunneling ionization rates for molecules. The tunneling ionization rates from the two lowest electronic states in H 2 + , i.e, the 1sσ g and 2pσ u states, were ab initio values taken from Ref. [42]. Because of the charge-resonance-enhanced ionization (CREI), for the peak probe intensity used in this work, the R-dependent ionization rates from the 2pσ u state are relatively large at internuclear distances around R = 7 and 11 a.u. Due to lack of more accurate published results, the tunneling ionization rates induced by the IR laser pulse from the excited states |φ el R,m in H 2 are approximated by the exponential factor in the tunneling We apply the split-operator method to the short-time propagator to solve Eq. (7) for H 2 and H 2 + as specified by the quantum jumps from H 2 to H 2 + and from H 2 + to H + + H + (see below). The practical application of the MCWP approach contains the following three steps. (i) In H 2 , starting from the ground state of H 2 (X 1 g + ,v = 0), we first evolve the initial nuclear wave packet along the potential energy curves in H 2 (we include 12 curves in this case) under the influence of the pump and probe pulses. The overall probability of staying in H 2 , i.e., P H 2 (t) = m dR|χ m (R,t)| 2 , decreases over time since the IR pulse induces ionization. The first jumps from H 2 to H 2 + are assumed to have a chance to occur at each time step and the ionization probability P 1 (t) within a small time interval t around t during the presence of the IR field is the drop in the probability P H 2 (t) at adjacent time steps, i.e., For a given first jump time t 1 , there are 12 individual pathways corresponding to jumps from the 12 involved potential energy curves to H 2 + and the relative probability from the potential energy curve for the |φ el Here g and u denote the 1sσ g and 2pσ u states in H 2 + , respectively, and M k (R) describes the field-dressed ground state in H 2 + and it is taken from Ref. [27]. Inserting this normalized initial nuclear wave packet into Eq. (7), we obtain the evolution of the nuclear wave packets along the 1sσ g and 2pσ u curves under the influence of the IR pulse. The ionization probability within a time interval t around t for the second jump from H 2 dR|χ k (R,t)| 2 is the overall probability of being in H 2 + at time t. For a given second jump time t 2 , the relative probabilities for the second jumps taking place from the two pathways, i.e., from the 1sσ g and 2pσ u curves to (iii) After the second ionization, since there is only the Coulombic curve in H + + H + , the term M m in Eq. (2) for the second jump operators is unity, thus the initial wave packet in H + + H + , corresponding to the case where the second jump takes place from the 1sσ g (2pσ u ) curve, is, before normal- ]. After removing two electrons from H 2 , the IR pulse has no influence on the evolution of the nuclear wave packets in H + + H + and the nuclear KER spectra are obtained by projecting the renormalized nuclear wave packets χ c (R,t 2 ) in H + + H + on the energy-normalized Coulomb waves K E (R), i.e., where m specifies that the first jump takes place from the |φ el R,m state in H 2 , k specifies that the second jump takes place from the 1sσ g or 2pσ u state in H 2 + , and it is implied that χ c (R,t 2 ) carries the jump history, i.e., the first jump time t 1 and the states |φ el R,m and |φ el R,k from which the first and second jumps take place. The total nuclear KER spectrum can be obtained through summing over the KER spectra describing jumps from different states in H 2 and H 2 + at different first and second jump times, i.e., In order to produce converged results, the time interval t should be small enough, i.e., t < π/(3 E) [43], where E can be approximated by the largest ionization potential in Fig. 1. Thus, there are hundreds of time steps when the simulation time is a few tens of femtoseconds. Accordingly, there are hundreds of first and second jumps in Eq. (9), thus without efficient sampling and simulation procedures, it would be very resource and time demanding to obtain the KER spectrum following double ionization of H 2 . In Appendix A we provide a detailed discussion and validation of our sampling method. For the parameters considered in Sec. III, the number of terms in Eq. (9) is reduced by a factor of 3.3 × 10 5 compared to the conventional sampling method, leading, of course, to a significant speedup of the calculations.

B. Validation of the MCWP approach in the XUV-IR scenario
In this section we validate the MCWP approach in the XUV-IR scenario. A complete test of the performance of the MCWP approach in simulating the double-ionization process in the XUV-IR scenario is impossible due to two facts. First, the process considered is not readily accessible by other theoretical approaches. Second, there are no experiments available for this process. This method can, however, be partially validated by comparing with available experimental or ab initio results for the single-ionization channel in the XUV-IR scenario. For instance, we have run a simulation under the same pulse parameters as in Ref. [44]: The two pulses have sin 2 envelops, the central photon energy of the pump (probe) pulse is 10.885 (4.63) eV, the pump (probe) pulse has a duration of 2 fs (4 fs), and the peak intensity of the pump (probe) pulse is 10 9 W/cm 2 (10 12 W/cm 2 ). The ionization probability following single ionization of H 2 is presented in Fig. 2 as a function of delay between the pump and probe pulses together with the ab initio and model calculations from Ref. [44]. Good agreement between the results for the three cases is observed, thus validating the application of the MCWP approach to the XUV-IR scenario. In our MCWP calculation, the ionization rates from the excited states induced by the probe pulse are calculated using the formula developed for the photoelectric effect [45].

III. RESULTS AND DISCUSSION
In this section, the nuclear KER spectra after double ionization of H 2 in an XUV-pump-IR-probe setup are obtained and analyzed. An isolated ultrashort XUV pulse is first applied to pump the ground state of H 2 to the excited states u + symmetry is unavoidable due to the bandwidth of the XUV pulse. This is the reason why we include six singly excited states of 1 u + symmetry and six of 1 g + symmetry in H 2 in our present model. The central wavelength of the IR pulses is 2400 nm, i.e., the IR photon energy is ω IR = 0.019 a.u. This energy is smaller than the ionization potential from H 2 to H 2 + for the highest-lying state in H 2 , thus the IR field can induce tunneling ionization (γ < 1) from excited states in H 2 when its the peak field strength is F IR0 = 0.04 a.u. The peak field strength of the XUV pulse is A XUV0 ω XUV = 0.005 338 a.u. To explore the influence of the XUV pulse duration on the excited nuclear dynamics in H 2 , we consider three different cases where the XUV pulse durations T XUV are one, three, and five optical cycles, corresponding to T XUV = 334, 1002, and 1671 as, while the pulse duration of the IR laser pulses is fixed at three optical cycles (T IR = 24 fs). We have found that a time step of t = 1 a.u. and a spatial step of R = 0.02 a.u. ensure convergence of the results presented here (see Appendix A).

A. Kinetic energy release spectra
We first show the nuclear KER distribution after double ionization of H 2 when the pulse duration of the XUV pulse is 043426-5 Ratios between the signals on the energy decreasing branches, which are apparent for delays larger than 50 fs with KER values decreasing from around 2 to around 1 eV, and the total yield for the three XUV durations as a function of delay. one optical cycle, i.e., T XUV = 13.7 a.u. = 334 as, in Fig. 3(a). When the delay is increased from 0 to around 10 fs, there are several branches with KER values decreasing from around 10 to around 5 eV. This behavior results from the nuclear wavepacket motion along the two potential energy curves in H 2 + , i.e., the 1sσ g and 2pσ u curves, moving from smaller to larger internuclear distances. Such a manifestation is also consistent with the fact that the proton yields are much smaller for smaller delays due to smaller ionization rates out of the 1sσ g and 2pσ u curves at smaller internuclear distances. For this small delay range (τ < 10 fs), the XUV pulse is within the FWHM region of the IR pulse where the IR field is strong, thus there is nearly no time for the nuclear wave packets in H 2 to move since they are ionized to H 2 + very quickly. For delays larger than 10 fs, the nuclear wave packets in H 2 have time to move along the excited potential energy curves in H 2 before substantial ionization from H 2 to H 2 + takes place. For the (10-30)-fs delay range, the KER values of the protons decrease from around 5 to around 4 eV, due to the spread of the nuclear wave packets in H 2 . Further increasing the delay from 30 to around 40 fs, there are significant enhancements of the proton yields for KER values around 3 eV. These enhancements result from the fact that the nuclear wave packets in H 2 + after the first ionization have time to reach internuclear distances around the CREI [46] position of R ≈ 11 a.u. where ionization takes place with a large probability [42]. The enhancements of the proton yields for KER values increasing from 3 to 5 and then decreasing back to 3 eV within the (40-60)-fs delay range is related to the nuclear wave packets moving from larger internuclear distances to smaller internuclear distances and then back to larger internuclear distances. When the delay is larger than 60 fs, the enhancements of the proton yields at KER values around 3 eV, i.e., in the (70-80)-fs range and in the (90-110)-fs range, have a similar origin as the above-mentioned enhancements of the signal at KER values around 3 eV within the (30-40)-fs delay range. This periodic occurrence implies that the vibrational period of the contributing nuclear wave packets is around 30 fs. Due to the spread of the nuclear wave packets in H 2 , the structures of the enhancements of the signal at KER values around 5 eV for delays larger than 60 fs are, however, not easy to observe. The energy decreasing branch with KER values starting from around 2 to around 1 eV for delays larger than 50 fs comes from the nuclear wave packets reaching internuclear distances larger than R ≈ 6 a.u. in H 2 , which, as we will see, are originating from the dissociating part of the nuclear wave packets evolving along the B 1  Fig. 3(a). This similarity is due to the fact that the XUV pulse is extremely short for the three cases and a very short XUV pulse creates nuclear wave packets in the excited states that resemble the initial ground vibrational state in the Franck-Condon approximation. When the XUV pulse duration is increased, this approximation breaks down and the nuclear KER spectra dramatically differ from Fig. 3. For the three panels in Fig. 3, there are small differences appearing at small delays (τ < 30 fs). For example, the periodic enhancements and suppressions of the proton yield every half of an IR laser cycle for τ < 30 fs become much easier to recognize when increasing the duration of the XUV pulse. Similar to those in Ref. [47], the periodic enhancements and suppressions of the proton yields originate from the population transfer between the ground state and the dominant excited state, i.e., the B 1 u + state, by the XUV pulse and its dependence on the IR field strengths around the instants of the short XUV pulse. If the IR field is weak (strong) around the instant of time when the XUV pulse is on, the Stark shift of the excited state by the IR field is small (large), thus the population transfer is large (small) due to resonant (off-resonant) coupling between the two states. When the population transfer from the ground state to the excited states in H 2 is large, the proton yields following double ionization are enhanced; otherwise they are suppressed. The reason why the periodic enhancements and suppressions become clearer for larger XUV durations is that the smaller bandwidth of the longer XUV pulse makes the population transfer more difficult for the same Stark shift. The main difference, appearing at large delays in Figs. 3(a)-3(c), is that the low-energy decreasing branch becomes less apparent when increasing the duration of the XUV pulse, which is 043426-6 because the smaller bandwidth makes direct excitation from the ground state to the B 1 u + state (responsible for this branch) less likely. We plot in Fig. 3(d) the ratios between the signals in the energy decreasing branches and the total proton yields for the three XUV pulse durations, which clearly show that the energy decreasing branch is relatively less important for a larger XUV pulse duration.

B. Trajectory and wave-packet analysis
Now we concentrate on identifying the origin of the aboveoutlined characteristic structures in Fig. 3. We first conduct a trajectory analysis for a specific characteristic feature in Fig. 3(a) as an example to show the strength of the MCWP approach in obtaining physical insights, including the stateand time-resolved information, of a given feature in the nuclear KER spectrum. Conducting a trajectory analysis is a unique possibility in the MCWP approach and it can tell us the individual contributions from all the involved trajectories to the signal for a specific KER value at a specific delay in the nuclear KER distribution. Here a trajectory means the complete dynamics leading to the final nuclear wave packet along the Coulombic curve. A trajectory is specified by four parameters: the first and second jump times t 1 and t 2 and the states from which the first and second jumps take place. Thus, by performing the trajectory analysis, we can, for any given feature in the KER spectra, identify the propagation in the neutral H 2 , the first jump time from the H 2 to the H 2 + system, the propagation in H 2 + , and the final second jump from H 2 + to the Coulomb explosion channel. It is interesting to analyze the nuclear KER spectra for τ ≈ 71.5 fs in Fig. 3(a) since at least two important structures are apparent for this value of the time delay τ , e.g., the enhancements of the proton yields at KER values around 3 eV and low-energy decreasing branch starting from τ ≈ 50 fs. When we carry out an analysis of the contributions to the signal at a given KER value in a nuclear KER spectrum, it is possible to obtain the origin of the signals at other KER values at the same time. This is because the signals for all the KER values in the same nuclear spectrum are formed from the same collection of trajectories with each trajectory giving its individual contribution to a particular feature. Thus, we present in Figs. 4(b)-4(d) the results obtained after a trajectory analysis to the signal for a KER equal to 3 eV at τ = 71.5 fs in Fig. 3(a). Figure 4(b) shows the respective probabilities for the first jump taking place from the states in H 2 for the considered energy (3.0 eV). We can see that the first jumps from the HH 1 g + and B3 1 u + states make the largest contributions to the signal. In Figs. 4(c) and 4(d), the contributions from the trajectories for the first jump taking place from the HH 1 g + and B3 1 u + states and the second jumps from the 2pσ u state are presented as a function of t 1 and t 2 , respectively. Contributions from the trajectories for the second jumps taking place from the 1sσ g state are not presented since they give much smaller contributions because the ionization rates from the 1sσ g curve are about two orders of magnitude smaller than those from the 2pσ u curve.
For a fixed first jump time, a clear periodic dependence of the dominant trajectories on the second jump time t 2 is observed in Figs. 4(c) and 4(d). These dominant trajectories correspond to ionization events taking place at instants around   Fig. 3(a) are mainly from ionization events occurring at the CREI position of R ≈ 11 a.u. from the 2pσ u curve in H 2 + . These evolution times in H 2 + imply that the contributing nuclear wave packets move slowly along the 2pσ u curve. This can be verified by the relatively small kinetic energies for these two trajectories before the second ionization taking place, which is the difference between the final kinetic energy and the Coulomb repulsion energy at around R ≈ 11 a.u., i.e., 3  To confirm this speculation, we show in Figs. 5(a) and 5(b) the evolution of the nuclear wave packets along the 1sσ g and 2pσ u curves in H 2 + when the first jump takes place at t 1 = 2193 a.u. (53 fs) from the HH 1 g + state. As seen from Fig. 4(b), the B3 1 u + state is as important as the HH 1 g + state. The nuclear wave packets motion along the two H 2 + curves for the first jumps occurring from the B3 1 u + state at t 1 = 2193 a.u. can be obtained and analyzed in the same way as the wave packets from the HH 1 g + state. We have found that the main features of the evolution of the nuclear wave packets in H 2 + are very similar irrespective of whether the intermediate HH 1 g + and B3 1 u + states are used and therefore here we only show the results associated with the HH 1 g + state in Figs. 5(a) and 5(b) (see Appendix B for the results associated with the B3 1 u + state). Now let us return to Fig. 5(b), where we can clearly see that the nuclear wave packet moving along the 2pσ u curve, starting from around R ≈ 6.2 a.u., reaches internuclear distances around the CREI position of R ≈ 11 a.u. at around t 2 = 3039 a.u. (74 fs).
We can also see from Fig. 5(b) that the nuclear wave packet at t 2 = 3039 a.u. is distributed over a very broad range of internuclear distances, i.e., from R = 5 to 15 a.u. Thus, apart from the above-mentioned nuclear wave packet at around R = 11 a.u., the other part of the nuclear wave packet would result in different structures in the nuclear KER spectrum. This means that a detailed analysis of this trajectory is worthwhile. Because of the R-dependent ionization rates out of the 2pσ u curve in H 2 + , nuclear wave packets along the 2pσ u curve located at internuclear distances with larger ionization rates are more likely to be ionized to H 2 2+ . Thus, for a given trajectory, the initial nuclear wave packet in H 2 2+ is dependent not only on the nuclear wave packet in H 2 + at its second jump time t 2 but also on the instantaneous ionization rates from H 2 + to H 2 2+ at t 2 , as discussed in Sec. II. In Fig. 6 we show the initial nuclear wave packets in H 2 2+ at different second jump times for the second jumps occurring from the 2pσ u state in H 2 + , when the first jump occurs from the HH 1 g + state in H 2 at t 1 = 2193 a.u. We can see that for the most dominant trajectory discussed above, i.e., for the trajectory where the first jump takes place from the HH 1 g + state at t 1 = 2193 a.u. and the second jump from the 2pσ u state at t 2 = 3039 a.u., there are enhancements of the nuclear wave packet densities at around R = 5, 10, and 15 a.u. By applying the reflection approximation, i.e., KER ≈ 1 R , such enhancements would result in KER signals peaking at around 5.5, 3, and 2 eV for this considered delay, respectively, which are consistent with the features in the nuclear KER spectrum for the 71.5-fs delay in Fig. 3(a).
To obtain knowledge of the initial wave packets in H 2 + we plot the evolution of the nuclear wave packets in H 2 along the HH 1 g + state in Fig. 5(f). The two regions with large probabilities at around t ≈ 1500 a.u. (37 fs) and t ≈ 2000 a.u. (49 fs) are from large dipole couplings between this state and the adjacent B 1 u + state. The evolution of the nuclear wave packets in the B 1 u + state is plotted in Fig. 5(e). We can see that the population of the nuclear wave packets in the HH 1 g + state is out of phase with that in the B 1 u + state for these two regions, i.e., t ≈ 1500 and 2000 a.u., directly showing a time-dependent population transfer between these two states. That is, at certain instants of time when there is a maximum of population in the Applying the reflection principle, the nuclear wave packets at around R = 5.5 a.u. along the 2pσ u curve in H 2 + [ Fig. 5(b)] at around t 2 = 3000 a.u. (73 fs) result in the enhancements of the proton yields at KER values around 5 eV for the 71.5-fs delay in Fig. 3(a), i.e., 1 5.5 × 27.2 ≈ 5 eV. In fact, for delays larger than 30 fs in Fig. 3(a), it is also the nuclear wave packets at around R = 5.5 a.u. that result in the enhancements of the proton yields at KER values around 5 eV. Moreover, the structure for KER increasing from 3 eV to 5 eV and then decreasing back to 3 eV in the (40-60)-fs delay range in Fig. 3(a) Fig. 3. As a result, one can benefit greatly from conducting a trajectory analysis to a given feature in the nuclear KER spectrum in uncovering its origin and the origins of other features in the considered nuclear KER spectrum. By conducting an analysis similar to that for the other characteristic structures in the nuclear KER distribution in Fig. 3, a thorough understanding of the formation of the structures can be achieved.

IV. CONCLUSION
We have extended the MCWP approach to simulate the fragmentation upon double ionization of H 2 in an XUV-pump-IR-probe setup where the ultrafast dynamics of singly excited states is involved. We have applied the MCWP approach to this XUV-IR scenario after providing a validation of this method in this scenario. This theoretical study accounts for the ionization induced by the intense IR pulse from the excited states in H 2 , reached from the ground state of H 2 by an XUV pulse, and again the ionization from H 2 + to H + + H + . We have obtained the nuclear KER distributions for three different pump-probe schemes using XUV pulses of durations of one, three, and five optical cycles, respectively, to excite H 2 for the pump-probe delays ranging from 0 to 140 fs. The nuclear dynamics in molecular excited states is reflected in the characteristic features in the nuclear KER distributions and, for instance, the periodic nuclear motion in the B 1 u + state is mapped out by trajectory and wave-packet analysis. Apart from that, the Stark effect can play a role in determining the efficiency of population transfer from the ground state to the excited states in H 2 for time delays τ < 30 fs. The present MCWP approach can be applied to XUV-pump-IR-probe attosecond spectroscopy of any molecule to resolve and control ultrafast dynamics in excited states through laser-induced double-ionization processes, provided the Born-Oppenheimer potential energy curves in different charge states, the dipole moment functions between different states, and the state-resolved ionization rates are available. In this appendix we give a detailed description and validation of the MCWP sampling method for XUV-pump-IRprobe spectroscopy. Within the MCWP approach, in order to reproduce experimental results, the time step t and the space step R should be chosen small enough. We have found that t = 1 a.u. and R = 0.02 a.u. give converged results for the processes studied in this work. Apart from that, in order to capture the features of the nuclear motion at large internuclear distances, such as the CREI [46] regions for the 2pσ u state, and also to reduce the reflection from the boundaries, the simulation box should be large enough, e.g., R box 40 a.u. For the approach in this work, the first jumps from H 2 to H 2 + take place at every time step and the second jumps take place from H 2 + to H + + H + at every time step after the first jumps. In addition, the first jumps have 12 pathways and the second jumps have two pathways, corresponding to the number of states included in H 2 and H 2 + , respectively. Thus, for an IR pulse with a duration of 24 fs (992 a.u.), there are at least 12 × 2 × 992 3 /3 terms in Eq. (9). The reduction of the huge number of terms in Eq. (9) to a computationally more feasible level is the task of our sampling method.

ACKNOWLEDGMENTS
In fact, the sampling method in our previous works involving only IR fields [10,18,22,[24][25][26] is very resource and time saving, since it assumes that the first jumps occur at the field 043426-9 extrema of the IR laser pulse and the second jumps take place at every tenth time step after the first jumps. In those works, laserinduced coupling between the ground and excited states in H 2 was neglected, since more than eight photons were required to resonantly excite the ground state to the lowest excited state in H 2 for the 800-nm wavelength. Such coupling became even more unlikely for mid-IR wavelengths. Thus it was acceptable to only include the ground-state potential energy curve, i.e., the X 1 g + curve, in H 2 in the calculations. The assumption for the second jumps works well since its sampling points are dense enough to capture the dynamics due to the relatively low velocities of the nuclear motion. The assumption for the first jumps is reasonable due to two facts. First, the nuclear wave packet along the X 1 g + potential energy curve in H 2 does not change significantly during the evolution. Such a stability of the nuclear wave packet results from the initial wave packet in H 2 being centered around the equilibrium internuclear distance, i.e., R ≈ 1.4 a.u., of the X 1 g + curve. Second, the very small ionization rates from the X 1 g + curve increase when increasing the field strength from zero to the peak field strength of the applied laser pulses. Thus, within half an optical cycle where the field strength increases from zero to an extremum and then decreases back again, the first jump occurring at the field extremum describes the events with the largest first ionization probability. Thus, under these assumptions, the number of terms in Eq. (9) is reduced to around 2 × 6 (first jump times) × 99 (second jump times) from 2 × 992 3 /3 for a 2400-nm IR laser pulse with a duration of 24 fs.
The above assumption for the first jumps, however, no longer works in the present XUV-pump-IR-probe case, because of the spread of the nuclear wave packets in the excited states of H 2 . Now the ionization probability P 1 (t) within a time interval t depends not only on the instantaneous field strength at time t but also on the spatial distribution of the nuclear wave packets in H 2 at t, due to the relatively rapid variations of the large ionization rates out of the excited potential energy curves in H 2 with respect to R. Thus, the first jump at a given field extremum does not necessarily have the largest probability among the jumps occurring within the same half of the optical cycle, as shown in Fig. 7. In Fig. 7(a), the population in H 2 and the ionization probability within a time interval t during the presence of the IR field, i.e., P H 2 (t) and P 1 (t), are plotted as a function of t for an illustrative case: The central frequency of the XUV (IR) pulse is ω XUV = 0.46 a.u. (ω IR = 0.019 a.u.), the duration of the XUV (IR) pulse is T XUV = 1 × (2π/ω XUV ) = 13.7 a.u. = 334 as [T IR = 3 × (2π/ω IR ) = 992 a.u. = 24 fs], the peak field strength of the XUV (IR) pulse is A XUV ω XUV = 0.005 338 a.u. (F IR0 = 0.04 a.u.), and the pump-probe delay is τ = 33.5 fs. For better observation of the dependence of P H 2 (t) and P 1 (t) on the instantaneous external field strength, the electric field of the IR pulse F IR (t) is plotted in Fig. 7(b). The t axis in Fig. 7(a) is confined from 0 to 1000 a.u., since the population in H 2 , P H 2 (t), starts to converge to 0.999 55 for t larger than 1000 a.u. in spite of the presence of the IR pulse. This means that nearly all the excited molecules are ionized by the IR field before t ≈ 1000 a.u. In our sampling method, we first determine the instants when the ionization probability P 1 (t) is locally maximized or minimized, which correspond to the instants for the local maxima and minima of the black line in Fig. 7(a). The first jumps are assumed to take place at the instants for the local maxima of the black line in Fig. 7(a). For a given first jump taking place at the instant for a given local maximum, the instants for the adjacent left and right local minima t 1L and t 1R are used to define its ionization probability under our sampling method. The new ionization probability P 1new (t 1 ) at a given first jump time t 1 now becomes the sum over the ionization probabilities within a time interval t at each time step between t 1L and t 1R , i.e., P 1new (t 1 ) = t=t 1R t=t 1L P 1 (t). Sampling of the second jump times is the same as in our previous works, i.e., assuming them to occur at every tenth time step. Instead of summing over the ionization probabilities within t for the ten time steps around t, i.e., t+10 t t P 2 (t; t 1 ), the new ionization probability for the second jump occurring at t is chosen the same as before, i.e., P 2new (t; t 1 ) = P 2 (t; t 1 ). This assumption is reasonable since the ratio between P 2 (t 2 ; t 1 ) and P 2 (t 2 + 10 t; t 1 ) is nearly the same as the ratio between t 2 +10 t t 2 P 2 (t 2 ; t 1 ) and t 2 +10 t t 2 P 2 (t 2 + 10 t; t 1 ), since the positions of the nuclear wave packets in H 2 + almost remain the same and the external field does not change significantly in ten time steps. Thus, compared with summing over the ionization probabilities within t for ten time steps as the ionization probabilities for the second jumps, our choice of the second ionization probabilities only influences the absolute values of the nuclear KER yields but keeps the shape of the nuclear KER spectra unchanged. Nuclear KER spectra obtained through applying the original sampling method (solid black curve) and the sampling method presented in this work (dashed red curve). The P E on the vertical axis labels the probability to produce a signal at a given energy. The laser parameters are identical to those for Fig. 7(a) except that a 0.5-fs pump-probe delay is used.
In order to validate our sampling method, we have conducted comparisons between results obtained from this and the original sampling method. As an example, the nuclear KER spectra obtained by these two sampling methods are presented in Fig. 8. The laser parameters are identical to those in Fig. 7(a) except that the pump-probe delay is τ = 0.5 fs. The excellent agreement between the two curves in Fig. 8 clearly shows that the sampling method presented herein works very well. Our method can reduce the number of terms in Eq. (9) from about 12 × 2 × 992 3 /3 to about 12 × 2 × 10 × 99. Thus, we apply the MCWP approach with this much more efficient sampling method to obtain the results presented in Sec. III. In this appendix we show in Fig. 9 the nuclear wave packets along the 1sσ g and 2pσ u states in H 2 + when the first jump occurs at t 1 = 2193 a.u. from the B3 1 u + state in H 2 . We can clearly see that the evolution of the nuclear wave packets for the first jump occurring from the B3 1 u + state (Fig. 9) is very similar to that for the first jump occurring from the HH 1 g + state [Figs. 5(a) and 5(b)]. This means that at a given second jump time t 2 , the trajectory whose first jump takes place at t 1 = 2193 a.u. from the B3 1 u + state would result in a nuclear KER spectrum similar to that associated with the trajectory whose first jump takes place at t 1 = 2193 a.u. from the HH 1 g + state.