Monitoring long-term evolution of molecular vibrational wave packet using high-order harmonic generation

We present a method for probing ultrafast nuclear dynamics in molecules using transient enhancement of high-order harmonic generation (HHG). The method exploits strong dependence of the overall harmonic yield in a fixed spectral range on the nuclear separation in a molecule, which is shown to take place for both aligned and randomly oriented molecules. Our numerical simulations show that this method is capable of monitoring long-term evolution of the nuclear vibrational wave packets in molecules, even in light-weight ones, with very high time resolution. By the example of D2+ vibrational wave packet launched via tunnelling ionization of D2 and probed by a time-delayed 8 fs laser pulse with λ = 800 nm and the peak intensity of 1014 W cm-2, we show that the time-delay dependence of the high harmonic signal exhibits pronounced features, which are due to the collapses and revivals of the nuclear wave packet. The time-frequency analysis of the pump–probe signal reveals structures with a periodicity down to 6 fs, which correspond to the fractional revivals of orders 1/5 and 1/10. With the laser parameters used, the deuteron motion during the probe pulse is shown to have almost negligible effect on the resulting signal.


Introduction
Recent progress in a methodology and technology for the time-resolved measurement of femtoand attosecond dynamics of electrons and nuclei in molecules has marked the emergence of a new branch of ultrashort pulse laser science, molecular ultrafast dynamic imaging (UFDI) [1]. Among others, high-order harmonic generation (HHG) proves to provide useful tools for molecular UFDI. The first theoretical study [2,3] and proof-of-principle experiments [4,5] have demonstrated the feasibility of extracting information on the molecular structure and its temporal changes from the harmonic spectra in molecular gases. As for the limitations of the above-mentioned methods, the probing attosecond dynamics by chirp encoded recollision (PACER) method [3,5] can provide reconstruction of the nuclear motion in a time window restricted to only a few femtoseconds [3]; the method relying on an analysis of the interference minima in the harmonic spectra from the molecules [2,4] is expected to be hardly applicable to light molecules, according to the semiclassical numerical simulations [6]. In this paper, we present an approach which can be used to monitor long-term evolution of the nuclear vibrational wave packets in molecules, even in light-weight ones. This approach exploits strong dependence of the overall harmonic yield in a fixed spectral range on the nuclear separation in a molecule.

Numerical simulation
We have carried out a three-dimensional (3D) numerical simulation of HHG in a D 2 molecular sample illuminated with two ultrashort laser pulses. The first one serves as the pump pulse preparing a coherent superposition of the D + 2 vibrational states via tunnelling ionization of D 2 , whereas the second time-delayed pulse is the driver for HHG, probing the nonlinear dipole response of the system at different times (see figure 1).

D + 2 nuclear wave-packet evolution
We consider a D 2 molecule initially in the ground state. As a result of fast ionization by a broadband pump pulse, the nuclear wavefunction is transferred from the neutral D 2 ground vibrational state to the lowest potential energy curve of the D + 2 molecular ion. The timedependent nuclear wavefunction can then be represented as an expansion over the set of D + (1) Due to the large bandwidth of the pump pulse, many vibrational states are excited whose energy eigenvalues can be expanded about the central value n 0 of the vibrational quantum number n as The nuclear wave-packet free evolution is illustrated in figure 2. In the first stage, the wave packet behaves classically, oscillating like a particle between the turning points in the potential well with the classical period T cl = 2πh/|E (n 0 )|. However, high anharmonicity of the D + 2 potential energy curve leads to a very rapid loss of spatial localization (collapse) of the wave packet due to the dephasing of its components. Unlike in the case of a wavepacket free motion, this scenario is not preserved for a long time. On the contrary, the wavepacket evolution in an anharmonic potential is well known to exhibit features called revivals [7]. Near the revival time T rev = 2πh/(|E (n 0 )|/2) (e.g. near t = 1180 fs for D + 2 , see figure 2) the different components of the wave packet are in phase again. As a result, the wave packet regains its original form and its evolution again becomes particle-like, resembling that at the initial stage. Moreover, as shown by Averbukh and Perelman [8], more subtle features called fractional revivals also take place (of which the most clearly seen in figure 2 are the half-revival near t = 590 fs and the quarter-revivals near t = 295 and 885 fs). The fractional revivals occur at times t = (m/n)T rev , where m and n are mutually prime integers. Near these times, the wave packet relocalizes in the form of q = (n/4)[3 − (−1) n ] 'clones' of the original wave packet separated by multiples of T cl q, each moving with the classical period (for the details, see reviews [9,10]). For example, near half a revival time the wave packet is reconstituted in its original form but shifted by half a classical period; near (1/4)T rev and (3/4)T rev (quarter-revivals) the wave packet consists of a pair of copies of the original packet antiphased with respect to each other. In fact, the reconstruction of the wave-packet form near the full and fractional revival times is only approximate, due to the contribution of the higher-order anharmonicities.

HHG in D +
2 molecular ion In our numerical study of the molecular dipole response we first considered that the nuclear positions do not change during the laser probe pulse. Using the split-operator fast Fourier transform technique [11] we calculated the time-dependent 3D wavefunction of the electron for the set of fixed values of the internuclear distance D and the molecular orientation angle θ with respect to the laser electric field with steps of D = 0.1 au and θ = 5 • . The laser electric field was considered linearly polarized along the z-axis. We numerically solved the time-dependent Schrödinger equation on a 3D Cartesian lattice of size 128 × 128 × 256 points with a uniform grid spacing x = y = z = 0.4 au. The time step was t = T /1440, where T is the laser period. We used the full two-centre molecular potential with Coulomb singularities. To treat the singularities, we took the grid points along the y-axis at half integer times y and put the molecular axis in the y = 0 plane. The nonlinear radiative response of aligned molecules was calculated via the time-dependent dipole acceleration expectation value [12] (atomic units are used throughout the paper). We have also carried out the calculations for the randomly oriented molecular ensemble according to From these quantities, the harmonic spectra were obtained by Fourier transformation into the frequency domain.
The probe pulse with λ = 800 nm was considered to have Gaussian shape with the peak intensity of 10 14 W cm −2 and full width at half maximum (FWHM) of 8 fs, which are the parameters accessible at present.  in each case the spectra calculated for the parallel (θ = 0 • ) and perpendicular (θ = 90 • ) orientations and randomly oriented molecular ensemble are plotted. As seen in figure 3, both the details of the harmonic emission spectra and the overall harmonic yield are very sensitive to the molecular orientation and size. The latter is explicitly demonstrated by figure 4, which plots the bond-length dependence of an integrated harmonic signal in the H19-H49 spectral range for various molecular alignment distributions.
The HHG process is well described in the frame of the three-step semiclassical model [13]. When an atom or molecule interacts with a strong linearly polarized laser field, the electron is optical field-ionized, accelerated by the ac electric field, and can recombine with the core to generate a high-energy XUV photon. Each step of this process contributes to a greater or lesser extent to the observed modulation of the harmonic signal with changing internuclear distance.
As figure 4 shows, the D-dependence of the harmonic yield is strongest for molecules aligned parallel to the laser field. Particularly, the contrast ratio between the strongest signal (at D > 3.5) and the weakest one (at D ≈ 2) is as high as about five orders of magnitude. This is mostly due to enhanced ionization, an effect which is essentially due to the fact that, in a certain range of the internuclear separation, the tunnelling of the electron through the internal barrier of the field-modified double-well molecular potential directly to the continuum dominates over the inter-well tunnelling [14]. The predominance of this effect in the range of D > 3 for θ = 0 • is confirmed by strong correlation between the bond-length dependence of the harmonic signal (figure 4, blue squares) and ionization probability by the end of laser pulse (figure 5). Since ionization provides the electrons necessary for high harmonic emission, strong variation of the ionization probability results in a corresponding change of the total harmonic yield (cf also figures 3(c) and (d)). Note that for the laser parameters used in our study, ionization induced by the probe pulse is well below the saturation level (even in the enhanced ionization case the ionization probability is about 0.1, see figure 5) and the bound-state depletion is insufficient to suppress harmonic generation. We also note that the saturation of the HHG signal for D > 3.5 is due to quite slow D-dependence of the ionization rate at 3.5 < D < 7 for the laser intensity I = 10 14 W cm −2 [15], rather than due to complete ionization of the system in this region.
The recombination step of the HHG process also influences the harmonic signal in our case through the two-centre interference effect [2]. As seen in figures 3(a) and (b) for the θ = 0 • case, the harmonics beyond the cutoff position (N > 35) and those around the 21th one are suppressed by an order of magnitude for D = 1.4 and 2, respectively. This suppression of the harmonic emission around some particular harmonic frequencies results from the destructive interference of the contributions to the emission from the two centres in a diatomic molecule. In the plane-wave approximation for the continuum electron, interference minima are expected for D cos θ = (2n + 1)λ/2 n = 0, 1, . . . , where λ = 2π/k is the de Broglie wavelength of the electron returning to the core. The frequency of the suppressed harmonic can be estimated as where the ionization potential I p is subtracted from the value determined by the simple-man model [13], = k 2 /2 + I p , in order to approximately take into account the influence of the attractive potential of the core. The positions of the interference minima observed in the 3D  numerical calculations agree well with (5) and (6) (see also figure 6, where the intensities of harmonics ranging from the 19th order to the cutoff in the spectrum are plotted).
The suppression of the plateau harmonics due to the two-centre destructive interference thus occurs in the range of internuclear distances between D ≈ 1.5 and 2.5, thereby diminishing the integrated harmonic yield for this interval of D, as observed in figure 4.
Importantly, for the randomly oriented molecular ensemble the D-dependence of the harmonic signal turns out to also be very strong (with a contrast of more than three orders of magnitude). This is because both enhanced ionization and two-centre interference effects, in fact, take place for all molecular orientations, except for the perpendicular one. For unaligned molecules the width of the minimum in the bond-length dependence of the signal is significantly larger (and, therefore, the region in which the harmonic yield grows dramatically with D is narrower), compared to the θ = 0 • case, since the interval of internuclear distances responsible for the destructive interference in the registration frequency window shifts towards larger D with increasing orientation angle.
As for the perpendicular orientation, the overall harmonic yield in this case, despite the absence of the aforementioned effects, also depends on the internuclear distance, although more smoothly. The efficiency of the generation of the plateau harmonics from the molecules aligned perpendicular to the electric laser field grows with increasing D in the range of D > 2.5. This behaviour can be explained by a reduced velocity of spreading of the electron wave packet during propagation in the continuum due to a larger delocalization of the initial wavefunction in a stretched molecule [16].

Probing the nuclear wave-packet evolution
To calculate the transient harmonic emission signal W(t) from a sample of D + 2 ions driven by the probe pulse delayed by time t relative to the pump pulse, we average the harmonic yield W(D) (see figure 4) over the nuclear probability distribution whose time evolution is shown in figure 2. The results for various orientations of the molecular axis are shown in figure 7.
As can be seen, for example, for the case of molecules aligned parallel to the electric field (figure 7(a)), the time-delay dependence of the harmonic signal plotted for the range of delays extending up to 1.3 ps exhibits intervals of quasi-periodic behaviour. The mean positions of these intervals and the corresponding instantaneous carrier frequencies agree well with those observed for the mean D value (figure 8). Similar time structures are observed for randomly oriented molecules (figure 7(c)). The above-mentioned features can be interpreted as arising from the revivals of the nuclear wave packet. Indeed, since near the fractional revival the wave packet consists of q equidistant mini-packets, the probability of the wave-function localization in the region of enhanced ionization oscillates with a period equal to T cl /q. The same should be  true for the harmonic signal. This conclusion can be verified by the Fourier and time-frequency analysis of the time-delay dependence of the HHG yield. Figure 9 shows the Fourier spectrum of the pump-probe signal where W(t) is the integrated harmonic signal in the H19-H49 spectral range sampled in the interval of delays 0 t T (T is the sampling time which we chose to be 1.3 ps). The spectrum was calculated for the unaligned sample of D + 2 ions. Several groups of lines localized approximately near the multiples of the classical vibrational frequency (in our case, ν cl ≈ 33 THz) are clearly seen in figure 9. The lth group arises from beats between vibrational states |v and |v − l (see, e.g., the analysis in [17] where the first experimental observation of the fractional revivals in a molecular system has been reported). Numerical ratios in figure 9 provide the identification of various groups of lines with certain fractional revivals near which these beats occur. The information on when the particular frequency arises can be extracted from the time-frequency analysis of the signal, as demonstrated for the case of fractional revivals in [17,18].  Figure 10 presents the result of the time-frequency analysis of the pump-probe signal W(t) for unaligned D + 2 . This analysis was made by means of the wavelet transform. We used the analysing function T t 0 ,ω (t) = √ ω (ω (t − t 0 )) with the Morlet wavelet (η) = 1/ √ τ exp (iη) exp −η 2 /2τ 2 (see, e.g., [19]) chosen as the mother wavelet; t 0 and ω are the variables used respectively to translate and scale the wavelet window. We chose the parameter τ equal to 5π/ √ 2 ln 2 ≈ 13.34, so that the analysing function FWHM comprises 5 oscillations. Figure 10 shows the contour plot of the wavelet energy density function (scalogram) W (t)T * t 0 ,ω (t) dt 2 , which in our case gives the contribution to the harmonic signal at the specific beat frequency ω and delay time t 0 . Vertical lines above the plot show the expected locations of various fractional revivals, as follows from an analysis in the spirit of the theory developed in [8]. One can easily see a series of patches centred about the frequencies relevant to various fractional revivals. The locations of these patches in time are in good agreement with the predictions of the theory by Averbukh and Perelman [8]. For example, the intense maxima centred about the frequency ω ≈ 33 THz corresponding to the full revival (near t = 1180 fs) and half-revival (near t = 590 fs) are clearly seen. Higher-order fractional revivals (see the patches around the frequencies 66, 99, 132 and 165 THz) are also easily identified in the scalogram. These observations confirm the interpretation of the features seen in the time-delay dependence of the HHG signal we mentioned above. We note that every patch in figure 10 is, in fact, split in time into the two ones, in agreement with the fact that the signal in figure 7 is suppressed at times half a classical period after the full and fractional revivals. This property, which can be used to identify accurately the revivals from the inspection of the signal, can be explained by the less synchronized harmonic emission near these times because of the higher delocalization of the outer tail of the nuclear probability distribution (the part of the wave packet in the shaded area of figure 11), as compared with that at the adjacent times of the wave packet return to the vicinity of the outer turning point (cf figures 11(a) and (b)).
The maximum central frequency of the patches observed on the time-frequency plane in figure 10 agrees with the estimated maximum number of resolvable fractions of the wave packet, q < L/ x, where L is the size of the classical orbit and x is the initial width of the packet [9].  Figure 11. Snapshots of the D + 2 nuclear wave packet at t = T rev + T class /2 (a) and t = T rev + 3T class /2 (b). Note that the observation of the fractional revivals of orders 1/5 and 1/10 for the D + 2 system, which has a classical period of approximately 30 fs, means that we resolve the time structures with the period T q about 6 fs (see the inset in figure 7(c)), i.e. shorter than the probe pulse duration. Given that at least four points are necessary in order to resolve one oscillation with the period T q , we can conclude that, in fact, 1-2 fs resolution of the method is demonstrated. This is because, due to the high nonlinearity of the HHG process, the time window to which the HHG is confined is very short compared with the duration of a laser pulse used as a probe.
Finally, we studied the impact of the nuclear motion during the probe pulse on the accuracy of the method. Since the pulse peak intensity is rather weak and the pulse duration is very short, we neglected the laser-induced changes of the positions of the nuclei, i.e. only the free evolution of the packet during the probe pulse was taken into account. The results are shown in figure 12.
It follows from the comparison of figure 12 with figure 7 that with the laser parameters used, the nuclear motion does not have any serious effect on the resulting signal. This can be explained by the fact that, because of a rapid increase of the efficiency of harmonic generation with the nuclear separation, the main contribution to the harmonic signal comes from the wave packet outer turning point, i.e. from the region where nuclei are at rest. Furthermore, the HHG process in our case lasts for a very short time as compared with the wave-packet classical oscillation period.

Discussion and conclusions
Although the approximation adopted in the computation for the moving-nuclei case and the explanation of the results shown in figure 12 seem reasonable, this reasoning can be in principle checked with more rigorous computations, i.e. via solving full TDSE, which is nowadays possible [15,20]. The calculations of this type could, for instance, help to clear up whether the probe field distorts the molecular potential via a bond-softening mechanism, which is the case [21] when the molecular dynamics is characterized using the time-resolved Coulomb explosion imaging (CEI) [22,23]. A similar question arises concerning the dissociation induced by the probe pulse. These points will be discussed in detail elsewhere. Anyway, since the laser peak intensity required for HHG is much lower than that used for CEI, the probing method discussed in this paper is expected to be much less perturbing than that based on Coulomb explosion. We should also mention that, while CEI provides richer information on the molecular dynamics [23], the HHG-based method has the advantage that the registration technique is much simpler than the advanced multiparticle detection technologies used for CEI [24]. We also note that, to the best of our knowledge, no clear experimental evidence for high-order revivals of vibrational wave packets in hydrogenous molecules has yet been reported in the literature.
Of interest is the question how the signal due to the effects discussed in the paper can be separated from the contribution from the molecules remaining nonionized after the pump is turned off. The answer depends on the parameters of the pump pulse. In the case of not too intense a pump, the neutrals will remain in the ground vibrational state and, therefore, will contribute only to the time-independent background which can be subtracted from the total signal. For shorter and more intense pulses the wave packet in the electronic ground state of the neutral molecule can be created due to the internuclear-distance dependent depletion of the initial wavefunction, the effect named 'Lochfrass' [25]. The wave packet created due to this effect was shown theoretically and experimentally to consist mainly of the v = 0 and 1 vibrational states [25,26]. The resulting periodic motion of the wave packet would lead to the contribution to the harmonic signal oscillating at a frequency equal to the difference between the v = 0 and 1 states of the neutral. This contribution can be eliminated using spectral filtering of the total signal. On the other hand, the observation of this frequency component in the signal could provide evidence of 'Lochfrass' in the pump-probe experiment using HHG for probing. The time-frequency analysis is also useful for distinguishing between the contributions due to the 'Lochfrass' effect and the high-order revivals of the wave packet in the molecular ion, since in the former case the oscillations are constant in time, while for the latter case the signal at a particular frequency occurs near particular instants of time. For even stronger pump pulses the contribution from 'Lochfrass' should drop because the number of nonionized molecules decreases with increasing intensity.
In conclusion, strong bond-length dependence of harmonic yield in molecules can be exploited for monitoring long-term evolution of molecular vibrational wave packets by means of HHG. In the case of a parallel orientation of the molecule with respect to the laser electric field, the harmonic signal is most sensitive to the molecular configuration. However, we have shown that the method can work efficiently even without molecular alignment. For the D + 2 molecular ion probed with 8 fs, 10 14 W cm −2 pulses, nuclear motion during the probe pulse is expected to play a minor role. 13 Finally, we should point to the experimental work [27] where the sensitivity of the HHG signal to the internuclear distance was demonstrated by comparing HHG from bound and dissociated iodine. We also note that in a recent experiment [28] the modulation of the pumpprobe signal measured for a particular (39th) harmonic was observed in SF 6 due to the vibrations excited by Raman scattering. Similar observation of a periodic enhancement of HHG due to laser-induced vibrations in argon clusters was reported in [29].