Quasi-classical model of non-destructive wavepacket manipulation by intense ultrashort nonresonant laser pulses

A quasi-classical model (QCM) of nuclear wavepacket generation, modification and imaging by three intense ultrafast near-infrared laser pulses has been developed. Intensities in excess of 1013 W cm-2 are studied, the laser radiation is non-resonant and pulse durations are in the few-cycle regime, hence significantly removed from the conditions typical of coherent control and femtochemistry. The 1sσ ground state of the D2 precursor is projected onto the available electronic states in D2+ (1sσg ground and 2pσu dissociative) and D++D+ (Coulomb explosion) by tunnel ionization by an ultrashort ‘pump’ pulse, and relative populations are found numerically. A generalized non-adiabatic treatment allows the dependence of the initial vibrational population distribution on laser intensity to be calculated. The wavepacket is approximated as a classical ensemble of particles moving on the 1sσg potential energy surface (PES), and hence follow trajectories of different amplitudes and frequencies depending on the initial vibrational state. The ‘control’ pulse introduces a time-dependent polarization of the molecular orbital, causing the PES to be modified according to the dynamic Stark effect and the transition dipole. The trajectories adjust in amplitude, frequency and phase-offset as work is done on or by the resulting force; comparing the perturbed and unperturbed trajectories allows the final vibrational state populations and phases to be determined. The action of the ‘probe’ pulse is represented by a discrete internuclear boundary, such that elements of the ensemble at a larger internuclear separation are assumed to be photodissociated. The vibrational populations predicted by the QCM are compared to recent quantum simulations (Niederhausen and Thumm 2008 Phys. Rev. A 77 013404), and a remarkable agreement has been found. The applicability of this model to femtosecond and attosecond time-scale experiments is discussed and the relation to established femtochemistry and coherent control techniques are explored.


Introduction
The first femtochemistry experiments of Zewail and co-workers temporally resolved the fragmentation of iodocyanide [1] into its atomic constituents. Soon after, the disintegration of tetrafluordiiodethane (C 2 I 2 F 4 ) into tetrafluorethylene (C 2 F 4 ) and two iodine atoms was studied, allowing the asynchronous breaking of two molecular bonds to be observed for the first time [2]. These pioneering studies relied on the generation of ultrashort laser pulses, typically produced by chirped pulse amplification in titanium: sapphire [3] and the subsequent stable and repeatable creation of two pulses, a 'pump' and a 'probe'. The most significant parameter of these pulses was that they were ultrashort (typically hundreds of femtoseconds, 1 fs = 10 −15 s or even shorter, and thus comparable to vibrational periods in most molecules), of a wavelength corresponding to a sufficient photon flux and energy to allow the coupling of molecular states (Ti:S operates in the visible IR, hence photon energies of a few eV are common allowing excitation or ionization by the resonant absorption of few photons; the probability of electronic transition is high for typical target number densities and photon fluxes) and stable shot-to-shot over hours, hence statistically significant measurements were possible, even for low probability transitions. 3 The highly complementary field of coherent or quantum control goes beyond observing wavepacket motion, looking to actively control chemical reactions by influencing the evolution of a nuclear wavepacket [4]- [12]. As with any quantum wavepacket, a nuclear wavepacket is a coherent superposition of states, which can exhibit a time-dependent localization (or wavepacket revival) of amplitude in one or more degrees of freedom, here the nuclear co-ordinate. Driving a quantum system from an initial state towards a desired final state is most readily performed by suitably shaped light fields, whereby the wavelength, flux or intensity, polarization and relative phase can all influence how a molecule proceeds from the reactant, via transition state(s) to the products (see [13] and references therein). As highlighted in a recent review by Wollenhaupt et al [6], the ability to dynamically guide or direct quantum state(s) in molecules has applications from bio-and solid-state physics via atomic and molecular physics and photochemistry, to quantum computing or cryptography.
To quantify the observations of femtochemistry and coherent control, novel theoretical treatments were developed to accurately quantify the evolution of nuclear wavepackets created and controlled by the light-matter interaction. At low intensities (generally I <10 10 W cm −2 ), the evolution of a nuclear wavepacket can be described by the perturbation theory and the population of the initial state is assumed to be constant. The Brumer-Shapiro scheme [14] facilitates wavepacket control by changing the photon energy as the system evolves, thus interfering different states; the Tannor-Kosloff-Rice scheme [15] employs a phase change induced by time-delayed quantum paths. In the intermediate intensity regime (which we define as 10 10 W cm −2 < I <10 12 W cm −2 ), the population of the initial and final states must be considered; the intensity regime is also dependent on the coupling strength between the electronic states. A comprehensive review of wavepacket dephasing and revival, quantum and quasi-classical computational treatments is given in [16]. Bichromatic and pump-dump control, many-photon interference processes, control over collisional processes and chaotic systems, the influence of decoherence and decay are covered in Shapiro and Brumer [5]. The analytical solution of the Schrödinger equation is often impossible, and hence numerical integrations or alternative representations of wavepackets are required: quasi-classical trajectories are often employed, which may be intuitively understood. The recent review of Bonacic and Mitric [17] unifies the interdisciplinary aspects of quantum chemistry, molecular dynamics, wavepacket propagation and optical control as applied to atomic clusters by discussing a number of illustrative polyatomic systems.
Rather than relying on changing resonance conditions to cause wavepacket interference, a number of intermediate intensity schemes have been proposed, whereby the intensity of the applied field broadens or shifts the energy levels, modifying the wavepacket. The adiabatic passage family of control schemes rely on the generation of light-induced potential energy surfaces through the application of an external field with at least two time-varying parameters. Examples are chirp + amplitude (rapid adiabatic passage, RAP), amplitude + delay (stimulated Raman adiabatic passage, STIRAP) and additional field amplitude + STIRAP pulses (Starkshifted-chirped rapid adiabatic passage, SCRAP). Yatsenko et al [18] have discussed the topology of all these techniques.
In the strong-field regime, which we define as I >10 13 W cm −2 , multiphoton or tunnel ionization become highly probable processes in atoms, and photodissociation and Coulomb explosion become dominant in molecules. Studying these processes has led to a significant advancement in understanding the dynamic structure of small molecules in intense laser fields, as recently reviewed by Posthumus [19]. While the interference of quantum pathways occurs under such conditions, it is often more useful and intuitive to think in terms of laser-induced 4 dipole forces. The electric field induced by the laser field is significantly strong that atomic and molecular orbitals are polarized and impulsive Raman processes are possible [20]. The propagation of liberated electrons in the laser field should also be considered, particularly in the case of high harmonic generation (see [21] and references therein). Under strong-field conditions it might seem unlikely that the delicate finesse required for quantum control would be possible; however, by applying the control pulse for a time significantly shorter than the fastest vibrational or rotational motions, the electronic structure of the molecule can be modified. This alters the environment in which the nuclei move, changing the motion of the nuclear wavepacket. An additional advantage is that ionization and fragmentation become less likely with shorter control pulses.
With the advent of few-cycle near infrared (NIR) laser pulses, a number of experimental groups have carried out studies to characterize vibrational wavepackets in hydrogenic molecular ions [22]- [28], as recently reviewed by the authors in [29]. Following illumination by an ultrafast strong-field pump, tunnel ionization of neutral H 2 or D 2 molecules leads to the formation of molecular ions supporting vibrational wavepackets. Such targets are attractive as they are theoretically tractable, and hence the Schrödinger equation can be solved within the Born-Oppenheimer and dipole approximations. Furthermore, the excited electronic states are isolated from the ground state of the D + 2 ion, and hence electronic excitation is negligible. As is often true of vibrational wavepackets in diatomic molecules, the anharmonicity of the potential energy surface causes the wavepacket to dephase, dissipating the vibrational state components over a range of internuclear separations within a few classical periods (≈ 16 fs for H + 2 , ≈ 24 fs for D + 2 ). Theoretical predictions were recently confirmed experimentally, as reviewed by Calvert et al [29].
The duration of the NIR pulses employed to generate and image vibrational wavepackets in hydrogenic molecular ions have also been demonstrated to produce rotational wavepackets in the neutral precursor [30,31]. It might be expected therefore that a vibrational wavepacket in D + 2 would actually be in a coherent superposition of rotational and vibrational states (i.e. a rovibrational wavepacket); however, the signature frequency components are not present in recent observations by a number of experimental groups [24,28,30]. As a result, treating the vibrational wavepacket one-dimensionally is valid, as the final orientation of the molecular axis can be selected experimentally. This is possible by measuring the molecular orientation using cold-target ion recoil ion momentum spectroscopy (COLTRIMS) [32], whereby the momentum vector of the fragmenting ions lies along the molecular axis. Alternatively, as demonstrated by Bryan et al [30], traditional ion time-of-flight mass spectroscopy can be modified by limiting the angular acceptance to a few degrees thus only allowing the detection of molecules with a known orientation.
Just as femtochemistry led to quantum or coherent control, the observation of nuclear dynamics in strong-field few-cycle pulses also opens possibilities for analogous manipulation of the quantum state of molecules. The present work discusses this regime which is on a much faster time-scale, at intensities in excess of an order of magnitude stronger and without the requirement of resonant transitions. Niikura et al have employed NIR laser pulses to perform experimental investigations of the laser-induced dipole forces to influence vibrational motion [33]- [35] in H + 2 and D + 2 , focusing on measuring the ion fragmentation energy as the pump-control delay is varied. Thumm and co-worker have reported an accurate quantum mechanical treatment of a vibrational wavepacket [36] and the influence of an ultrashort control pulse [37,38]. The authors in collaboration with Murphy and McCann have used a similar 5 numerical treatment to investigate the heating and cooling of vibrational populations [39]- [41], recently leading to the prediction of the 'quantum chessboard' effect [42]. As will be discussed later in detail, the underlying process facilitating wavepacket control can either be considered as an impulsive Raman process or a dynamic Stark shifting of the potential energy surface induced by polarizing the molecular orbital. Such processes have recently been employed by Stolow and co-workers [43] to control the outcome of the laser-induced photodissociation of IBr. A control pulse is applied to modify the non-adiabatic evolution of a nuclear wavepacket as it traverses the resulting avoided crossing in the potential energy surface. As in the present work, the electronic states are isolated, and hence the process is reversible. An interesting application of wavepacket manipulation is employing the internal states of a molecule for quantum information processing [44]- [46].
Another topic attracting significant theoretical and experimental interest is the intentional localization of an electron to a chosen nucleus within a molecule [47]. This is directly related to the present work: while electron motion occurs in tens or hundreds of attoseconds, by applying a few-cycle carrier-envelope phase (CEP) stable control pulse, the electric field induced between the nuclei is sufficient to transfer electron density by overcoming the internal Coulomb barrier. Electron density is forced to oscillate at the frequency of the control pulse, and if the pulse is sufficiently short (sub-10 fs), the probability of the electron being localized to a chosen nucleus can be in excess of 80% [48]. By varying the CEP offset by π, the election is localized from one nuclei to the other. If we now consider the motion of the vibrational wavepacket, it is clear that the electronic and nuclear wavepackets will be coupled; if the control pulse is applied around the revival of the nuclear wavepacket, the fast oscillation of the electron density between the nuclei is modulated by the slow motion of the nuclear wavepacket. This effect is most dominant at the vibrational revival, as discussed by the authors in [48].
Considering electron localization naturally leads us to consider electronic processes on attosecond time-scales: UV-XUV pulses allow the generation of vibrational wavepackets in excited electronic states or molecular ions [49]. However, the intrinsic ultra-broadband nature of such radiation will make studying many-electron molecules troublesome, as all electronic states within the laser bandwidth can be populated. We address this point as we conclude, as nuclear dynamics may prove influential when identifying populated electronic states. We also refer back to femtochemistry and coherent control, as there is much to be learnt for applying trajectory approximations to many-electron and polyatomic molecules.

The quasi-classical model (QCM) in one dimension
The QCM is broken into a series of steps. The creation of a coherent superposition of states is modelled by extending the non-adiabatic ionization model of Yudin and Ivanov [56] to a molecular system, allowing the vibrational state population to be evaluated. A non-interacting classical ensemble is then created, with weighted quantized vibrational levels to reflect the initial state of the molecule. The wavepacket evolution is then approximated by allowing the classical ensemble to propagate on the internuclear potential energy surface (PES). Now, applying a 'control' pulse to the propagating ensemble causes a time-dependent deformation of the PES [33]: the large gradient of the induced electric field causes a large AC Stark shift. Essentially, the interaction of the induced field and the electron orbital causes a force on the nuclei. The resulting deformation of the potential accelerates or decelerates components of the ensemble depending on their direction of motion, transferring energy into or out of the system.
Moving away from the classical interpretation, this transfer of energy is directly analogous to an impulsive Raman process, except that it is continuous (rather than discrete) and applies throughout the duration of the control pulse. The influence of this energy transfer is twofold, changing the relative phase of the ensemble components with respect to the unperturbed system, and transferring population between vibrational states. The transfer of population is considered to occur when the motion of an ensemble component is sufficiently perturbed by the control pulse that it takes on the characteristics of a higher-or lower-lying vibrational state (i.e. amplitude and frequency). Each of these processes will be discussed in detail.

Molecular ionization in a few-cycle pulse
As with recent experiments [22]- [28], with this prototypical molecular ion the vibrational wavepacket is created by the multiphoton or tunnel ionization of the precursor neutral molecule by a few-cycle NIR (800 nm) pump pulse. There is no delineation between the ionization mechanisms; however, when the Keldysh parameter γ > 1 is generally accepted to be in the multiphoton regime and when γ < 1 is in the tunnel regime, where γ 2 = I p /2U p , where I p is the ionization potential of the electron under consideration and the pondermotive potential where ε 2 is the laser intensity and ω L the laser frequency; atomic units are assumed here. The MO-ADK (molecular orbital Ammosov-Delone-Krainov) formulation [51] has been relatively successful in predicting ionization rates; however, any theory that makes use of the original ADK treatment [52] is not ideally suited for a few-cycle pulse, being a quasistatic approximation. Instead, here we modify the recent work of Yudin and Ivanov [56] in which an analytical expression for ionization rate is found for arbitrary values of the Keldysh parameter. The flexibility of these expressions allows us to treat ionization between electronic states of the molecule efficiently.
The distribution of vibrational states depends on the probability of ionization and the overlap between the ground state wavefunction in the neutral and the available vibrational states in the ion. The recent experimental findings of Urbain et al [57] and theoretical discussion by Kjeldsen and Madsen [58] identified a significant deviation from the expected Franck-Condonlike distribution of vibrational states being populated by strong-field ionization in an ultrafast NIR laser pulse. This deviation was found to be dependent on peak laser intensity: on the leading edge of the pump pulse, the D 2 → D + 2 (1sσ g ) ionization rate, g D 2 (R) increases highly nonlinearly with laser intensity. The distribution of vibrational states is a direct result of the dependence of the ionization potential on internuclear separation, such that the maximum ionization rate (hence maximum population) does not necessarily occur at the optimum overlap of the ground state neutral wavepacket and the states in the molecular ion. Furthermore, once the intensity is sufficiently high to generate the molecular ion, there is also a significant excitation rate to the D + 2 (2pσ u ) dissociative state, u g (R). This is a loss channel and depopulates the number of D + 2 ions supporting vibrational wavepackets. An additional loss channel exists in D + + D + Coulomb explosion, CE g (R), still requiring a higher intensity. Figure 1(a) shows the result of numerically evaluating the coupled differential ionization rates and integrating over the pump pulse as a function of time. This result is comparable to that of Thompson et al [53] and indicates that a relatively narrow intensity range generates stable D + 2 ions. The distribution of vibrational population over this intensity range is found by calculating the rates between the neutral precursor and the 1sσ g PES as a function of R and scaling by the ground state nuclear wavefunction in the neutral shown in figures 1(b)-(g). In all cases, a 7 fs, 800 nm pump is simulated. The populations are found to vary subtly, most visibly for the , photodissociation (0,1) and Coulomb explosion (1,1) probabilities, an indication of the total number of ions generated at a particular intensity. The probability of surviving the pump pulse unaffected is also indicated (D 2 ). (b)-(g) Relative vibrational state populations as the intensity of the pump pulse is increased from (b) 5 × 10 13 W cm −2 to (g) 1 × 10 14 W cm −2 . The Franck-Condon distribution of states is indicated by the thin line. The largest discrepancy between Franck-Condon and the predicted distributions occurs at low pump intensity.
lowest-lying states where the potential energy gap varies most rapidly with internuclear separation, hence changing the rates g D 2 (R), u g (R) and CE g (R), in agreement with [50,58]. Clearly, it is vital that the distribution of laser intensity throughout the pump laser focal volume be known (as demonstrated in [54,55]), as this modifies the probability distribution in figure 1(a). A marked difference is found between the predicted distribution of vibrational states and the Franck-Condon distribution of states expected from charged-particle impact or long laser pulse ionization, also in agreement with [57] and [58]. 8 While this may not be an exact description of molecular ionization, it is a useful treatment as it is a reasonable approximation to the underlying physics, is mathematically straightforward, is applicable to both the multiphoton and tunnel ionization regimes and holds for ultrashort laser pulses. This is also the first use of Yudin and Ivanov's treatment of ionization for a molecular system. Furthermore, it should be pointed out that this model does not take into account the variation with ionization rate with angle with respect to the molecular axis; however, it hopefully conveys an insight into the underlying mechanism. As discussed later, the QCM can be extended to complex multi-electron systems on the condition that the PES of the neutral and ionic states are known along the coordinate frames of interest, such that the ionization potential can be calculated at all points. Furthermore, ionization by an attosecond UV/XUV pulse can also be included, as long as the active electronic transitions and associated cross-sections are identified across the full bandwidth of the pump pulse. The photon flux could then be converted directly to electronic state populations in the manner employed in synchrotron studies.

Propagation of ensemble components on the potential energy surface (PES)
In a manner directly comparable to the trajectory methods pioneered in femtochemistry and coherent control as highlighted in the introduction, we treat the evolution of the vibrational wavepacket in the molecular ion as an ensemble of classical particles moving on the bound PES. The initial internuclear positions of the elements of the ensemble are defined by the intersection of the ground electronic state of the molecular ion and the energy corresponding to the mid-point between successive vibrational states v = 0, 1, 2, . . .. For example, the (eigen)energy of v = 3 is −0.0782 au (where zero energy is defined by dissociation asymptote) and the (eigen)energy of v = 4 is −0.0719 au. The initial position of the element of the ensemble corresponding to v = 3 is R = 1.439 au, found by taking the mid-point of these two energies and finding the corresponding point on the PES.
The ionization event that launches the vibrational wavepacket enforces a coherence to the resulting motion, as the ionization rate is maximum for all vibrational states in the ion at the inner turning point. It is assumed that the unit mass of the classical particle is initially stationary following projection onto the ionic PES, which is reasonable considering that the tunnel time is much shorter than the characteristic time-scale for wavepacket propagation. We then treat the motion of each element of the ensemble as a Newtonian particle accelerated according to the differential of the PES. The presented discussion is confined to one dimension to be applicable to the hydrogenic molecular ion; however, as discussed in detail in Bonacic and Mitric [17], this holds for any number of degrees of freedom, and hence has significant promise for quantifying wavepacket motion in complex polyatomics exposed to strong-field pump pulses causing excitation or ionization.
It is assumed that the PES is smooth and well behaved, such that the numerical differential of a relatively coarse sample of the PES is an accurate representation of the continuous differential. In the one-dimensional case, a five-point stencil is employed; the key requirement is that the discrete numeric differential can be reliably and efficiently interpolated. In this work, a natural Spline function is employed; however, it is suggested that the PES be sampled to a sufficient spatial resolution so that the selection of the interpolation method is essentially arbitrary on the condition that a smooth continuous output results. Now, by discretely calculating the equations of motion of the unit mass as influenced by the differential of the PES, the ensemble propagation is simulated as a function of time.

9
The motion of the ensemble (and hence the approximation to the wavepacket motion) is then predicted by allowing the unit masses to propagate as a function of time. The result of propagating the QCM for a typical number of vibrational states (generally from v = 0 to v ≈ 24) is presented in figure 2, and a visual comparison made to the solution of the time-dependent Schrödinger equation (TDSE) for D + 2 within the Born-Oppenheimer and dipole approximations, as published in [26]- [28], identical to that presented in comparable work [22]- [25]. Clearly, as the QCM generates a series of trajectories and the quantum treatment results in a continuous wavefunction, a direct numerical comparison is difficult. However, from figure 2 it is apparent that the QCM captures the essence of the quantum wavepacket result. This is not surprising, as both motions are the direct consequence of the PES shape as a function of R; however, this comparison is intended to indicate the convergence of the classical trajectory calculation despite only requiring 7000 simulation steps to model the wavepacket motion over 700 fs. A typical split operator propagation of the TDSE (where a solution exists) requires of the order 30 000 steps, i.e. a time-step of 1 au (24.2 as). The classical trajectory model can therefore be executed on a desktop PC in seconds, and hence is highly attractive for immediate comparison with experimental results in the laboratory. Furthermore, it is straightforward to numerically integrate the wavepacket motion over the distribution of laser intensities in the laser focal volume within a minute. Figure 2 presents two subsets of the now familiar wavepacket dephasing and revival. Comparing the QCM trajectories and the predicted TDSE-derived wavepacket amplitude at 590-600 fs, the intense features in the false colour representation of the wavepacket are apparent in the overlapping or closely spaced QCM trajectories. The slightly imperfect rephasing of the wavepacket predicted by the TDSE solution is replicated in the QCM trajectories.

Quantifying deformation of the PES
Weak and intermediate field modifications of wavepacket motion rely on coupling vibrational states through other PESs, possibly with relatively low intensity (hence 1 eV) Stark dressing of the states. In the strong-field regime, the optically induced electric field is sufficient to not only distort the PES via the Stark shift, but also to polarize the molecular orbital(s), which the nuclei experiences as a time-dependent force, therefore changing the potential energy. The orbitals respond to the electric field (i.e. on an attosecond, or at least few femtosecond time-scale), but the nuclei are only driven by the cycle-averaged envelope; hence, if a control pulse is applied to a molecule supporting a vibrational wavepacket, the nuclear motion adjusts according to the external influence. Consequentially, our control pulse will modify the effective PES by the dipole force generated by polarizing the molecular orbital and subsequential Stark shifting.
Niikura, Villeneuve and Corkum were the first to explore this process in D + 2 [33]- [35], for reasons of accessibility and electronic simplicity as identified in this and other recent work. It is assumed that the structural dynamics of hydrogenic molecular ions can be accurately predicted by wavepacket motion on, and coupling between, the two lowest-lying potential energy surfaces (bound state 1sσ g and dissociative state 2pσ u ). This approach is not unique to the QCM, but is accepted because the 1sσ g and 2pσ u states are energetically isolated from the electronically excited states of the ion. Clearly, this will be a contentious issue when dealing with manyelectron molecules, as a large number of overlapping bound and dissociative states exist; this point is addressed later. Following Niikura and co-workers, we quantify the distortion of the 1sσ g state [59] as where = R|E(t)|/2, and we define E(t) as the time-dependent electric field envelope induced by the control pulse, and hence is a function of I control (t). The V (R) term allows the transition dipole to be estimated, and gives the dynamic Stark shift. Clearly, the influence of the transition dipole is maximum at small R and the Stark effect is dominant at large R. Intuitively, this can be understood by considering the polarization of the molecular orbital: the largest distortion of the orbital occurs well away from the nuclei (hence causing ensemble trajectory modification and loss through dissociation); however, the nuclear motion is influenced by comparatively small distortions all the way into in the vicinity of the PES minimum. As the modification of the PES is numerically propagated, it is trivial to include an arbitrarily complex control pulse. The only condition is that the temporal step size is sufficiently small that the envelope varies smoothly.
In figure 3(a), the field-free PES and the maximally laser-distorted PES are shown; such a dramatic deformation of the PES is only present for a fraction of the vibrational period. The distortion of the PES by the control pulse causes a time-dependent variation of dV (R)/dR, and hence the ensemble components experience an additional acceleration, the direction and magnitude of which depends on the direction the component is moving and its location on the PES. This causes a net increase or decrease of energy (analogous to an impulsive Raman process) depending on whether the control pulse dipole force acts to enhance or retard the ensemble elements, interpreted as a change of vibrational state in a quantum representation. This is illustrated in figure 3(b), whereby a 7 fs FWHM Gaussian control pulse of peak intensity I control = 3.5 × 10 13 W cm −2 is applied 24 fs after the wavepacket is created by the pump pulse (i.e. t = 0). The trajectories around R = 5 au at 24 fs illustrate this point well: both are accelerated to larger R under the influence of the control.
If the potential is too dramatically distorted, elements of the ensemble are no longer bound by a concave PES, and dissociation occurs (as seen in the case of the five highest vibrational states in figure 3(b)). However, the dynamic nature of the distortion should be considered: if the control pulse is shorter than approximately a quarter of the period of a particular vibrational beat frequency, that component of the ensemble is resilient to dissociation if it is in the vicinity of the inner turning point. Furthermore, the likelihood of dissociation also depends on the position of the components in R and motion along the potential 'trough' caused by the control pulse.

Variation of vibrational state phase and population
As touched upon at the beginning of this section, applying a control pulse causes a timedependent distortion of the PES resulting in a dipole force acting along the internuclear axis. As the wavepacket (here the ensemble of quasi-classical particles) propagates, this force transfers energy into or out of the system, analogous to an impulsive Raman process. Clearly, the discrete vibrational states have an influence when considering the quantum nature: a change of vibrational state only occurs when sufficient energy has been gained or lost to instigate a change in state.
Recent time-dependent solutions of the Schrödinger equation [36], [38]- [41] quantify this process by calculating the dipole coupling between the 1sσ g and 2pσ u PESs in the coupling operator. Here we use a similar approach by allowing the components of the ensemble to rapidly adjust to the modified PES. Clearly, the diabaticity of this process is important; however, we have already assumed the nuclear motion is driven by the control pulse envelope rather than the few-femtosecond oscillation of the associated electric field. Recent experimental investigations by the authors support this assumption [60], where we report that the presented model accurately reproduces the modification of the recorded photodissication yield. As is illustrated in figure 4, the time-varying distortion of the molecular PES by the control pulse causes the component trajectories to vary from the unperturbed motion. If additional energy is transported into a trajectory, the amplitude of the oscillation increases, along with the vibrational period, and hence the trajectory takes on the characteristics of a higher-lying state. We propose that this increase in energy is equivalent to population being transferred from the a lower to higher vibrational state. Similarly, the control pulse can simultaneously remove energy from the system, causing a decrease in vibrational state.
Intuitively, the transfer of population between vibrational states can be understood by considering whether the action of the control pulse is to impede or encourage the component. In figure 4, the control pulse is applied at such a time that the unperturbed trajectories of v = 18 and 7 are just before and after arriving at the outer turning point, and hence have velocities in opposite directions. The force resulting from the modification of the PES acts away from r = 0 in figure 4; hence, following the control pulse, the trajectory of the component initially in v = 10 has the characteristics of v = 7 (i.e. is vibrationally excited). The opposite is true for the component initially in v = 13 as compared to v = 18.
As the control pulse can up-and down-shift vibrational population, it is possible that multiple ensemble components can end up in the same vibrational state, and the phase and relative population of each state must be considered. Variations in the phase and population are quantified by making a numerical comparison between the perturbed (e.g. figure 3) and unperturbed (figure 2) trajectories. Importantly, this comparison must be made once the intensity of the control pulse has dropped to zero, such that the trajectories are only defined by the  static PES. In practice, it is important to ensure that the control pulse(s) does not suffer from significant pedestals on the femtosecond and picosecond time-scale. Therefore a third-order autocorrelation contrast measurement is required when comparing with experimental results, in case an underlying PES distortion exists. Using the unperturbed trajectories as a reference (figure 2), the final vibrational state and phase are found by a global form of least-squares fitting which is unbiased towards the initial state of the trajectory. The square of the difference between unperturbed and perturbed motions is calculated, defining a quality of fit parameter and a time offset is introduced to the unperturbed trajectories and varied over a range greater than the period of the highest vibrational state. The minimum values of the fitting parameter therefore returns the distribution of final vibrational states and time offset that best represent the ensemble after the action of the control pulse.
The shifting of population is not treated in a totally discrete manner; rather, if a final state is equally represented by two states, the initial population will be shared according to the ratio of the fitting parameters. To determine the final distribution of vibrational states, the populations in each state are summed; this is illustrated in the vibrational population matrix shown in figure 4(b). The redistributed initial population (caused by a pump intensity of 10 14 W cm −2 ) by the control pulse is apparent. The vibrational population matrix is a convenient visualization of the influence of the control pulse. With no control pulse present, the vibrational population matrix would be a diagonal line as v in = v out .
On transferring between states, we place no constraint on the relative phase between the unperturbed and perturbed motions of the trajectories. So, the resultant phase with respect  [38] are reproduced for comparison with the QCM, which reproduces the periodicity and relative shift of vibrational population as the delay between the pump (6 fs, 1 × 10 14 W cm −2 ) and the control is varied.
to the natural motion can be varied depending on when a vibrational state undergoes state transfer, hence dependent on the intensity and delay of the control pulse. If a number of ensemble components contribute, a weighted mean is calculated, depending on the relative initial population. There is a clear phase difference between the vibrational motions highlighted in figure 4(a) following the control pulse. Depending on the intensity and arrival time of the control, this phase difference can be positive or negative. The phase change caused by the control pulse is a natural result of the fitting technique where the time offset is directly related to the phase via the known periods of the unperturbed motions.

Comparison with established results
While figure 2 has demonstrated that the QCM can successfully reproduce the unperturbed motion of a vibrational wavepacket by propagating the ensemble, it could be argued that any treatment relying on the differentiation of the field-free PES will give a similar result, assuming sufficient numerical accuracy. It is therefore vital to quantify the ability of the QCM to accurately describe the action of the control pulse; to this end, we compare the output of the QCM to established theoretical results. In figure 5, systematic scans of the control delay and intensity with respect to the pump pulse are presented for direct comparison with the reproduced results of Niederhausen and Thumm, figure 5 [38]. In both cases, the control pulse is 6 fs in duration and the initial vibrational state distribution is calculated using the tunneling theory. Importantly, the initial distribution of vibrational states in [38] was calculated using the ADK treatment [52]. The TDSE was then solved within the Born-Oppenheimer approximation, and the control pulse causes Raman transitions between vibrational states. Furthermore, following the control pulse, the wavepacket is projected onto the Coulomb explosion PES. Despite these differences in method, a good agreement is found (especially considering how disparate the two numerical techniques are), indicating that the QCM accurately captures the modification of the vibrational population as revealed by a full quantum mechanical model. Approximating the evolution of a wavepacket under the influence of an ultrafast strong-field control pulse as a quasi-classical ensemble is therefore justified.

QCM prediction of vibrational population and phase
In figure 6, the vibrational population matrices, initial and final vibrational state populations and final phases are shown as a function of control pulse intensity. The initial vibrational distribution corresponds to a pump pulse intensity of 1 × 10 14 W cm −2 , and the control delay is fixed at 25 fs. At the lowest intensity of 1 × 10 13 W cm −2 , the near diagonal distribution of vibrational population indicates the initial and final states are similar, as the control pulse perturbation is minimal. Nonetheless, the phase of the ensemble is altered, inducing a step change in phase between v = 4 and 7, which shifts to lower v-state as the control pulse intensity is increased. At the highest control intensity of 5 × 10 13 W cm −2 , the matrix indicates a dramatic redistribution of population, skewing the maximum in the population from v = 1 to v = 3 and 4. As is apparent from the matrix, initial states with v > 12 are dissociated by the relatively intense control pulse, and the phase is maximally distorted at high v. Figure 7 shows how the vibrational population matrices, initial and final vibrational state populations and the final phases vary as the delay between the pump (1 × 10 14 W cm −2 ) and control (3 × 10 13 W cm −2 ) pulses is varied in 4 fs steps. Even with a fixed control intensity, the degree by which the final vibrational population and phase can be varied is dramatic. Changing the control arrival time from 24 to 28 fs then to 32 fs induces a significant shift, which we suggest should be experimentally observable.
The result of applying a control pulse can be resolved by imaging the wavepacket (i.e. by destroying the molecular bond). As with previous studies, the temporal evolution of the wavepacket is revealed by an intense probe pulse which initiates photodissociation or Coulomb explosion of elements of the wavepacket at large R. The overall photodissociation yield is therefore directly proportional to the (amplitude) 2 of the wavepacket. Furthermore, measuring the kinetic energy of fragmentation reveals the location of the wavepacket in R when the probe pulse arrives, allowing an additional comparison with theory. As discussed earlier, an intensity of greater than 2 × 10 14 W cm −2 would be optimal to access the majority of bound states. Varying the delay between the pump and control changes the distribution of states in the superposition. Varying the delay between the control and probe, then Fourier transforming the resulting time-dependent fragmentation yield should reveal the influence of the control pulse.
In the final stage of the QCM, predicting the fragmentation yield can be done by modelling the distortion of the PES by the probe pulse, then identifying the vibrational states that do not survive, i.e. are dissociated. If the probe intensity is high enough, this will be almost all states. However, this calculation is time consuming: a far more efficient method has been recently demonstrated [48]. The dissociation yield can be accurately estimated for a chosen pulse intensity and duration by identifying the elements of the ensemble that are above some critical Initial vibrational state, v in (a) 25fs 1x10 13 (b) 25fs 2x10 13 (c) 25fs 3x10 13 (d) 25fs 5x10 13 Figure 6. Population transfer matrices, vibrational populations (thick solid lines) and phase relative to the unperturbed motion as the intensity of the control pulse is varied. The control arrives 25 fs after the ensemble is generated, and is 7 fs in duration. The vibrational population generated by the 1 × 10 14 W cm −2 pump pulse is indicated by the thin line. On the phase plots, the unperturbed wavepacket has a relative phase of zero, the resulting phase is the population weighted mean and is plotted twice separated by 2π to aid visual interpretation.
internuclear separation during the probe pulse, referred to as the 'critical R cutoff' method. In the case of the hydrogenic molecular ions, an R cutoff of 2.9 au is employed, corresponding to a probe intensity of 10 14 W cm −2 . and phase relative to the unperturbed motion as the arrival time of the control pulse is varied. The control pulse has an intensity of 3 × 10 13 W cm −2 , and is 7 fs in duration with a centre wavelength of 800 nm. The vibrational population generated by the 1 × 10 14 W cm −2 pump pulse is indicated by the thin line. On the phase plots, the unperturbed wavepacket has a relative phase of zero, the resulting phase is the population weighted mean and is plotted twice separated by 2π to aid visual interpretation.

Applicability of the QCM to complex molecules
As has been demonstrated, when interpreting coherent control and femtochemistry experimental results (see [4]- [12] and references therein), starting from a simple system and then considering multiple active electronic states and more than two nuclei leads to rather more computational complexity. Treating the motion of polyatomic molecules on a known PES is well understood, with the key consideration being the reduction to the most significant nuclear coordinates.
Assuming the vibrational wavepacket propagates on a single known PES, applying the QCM to a polyatomic molecule would rely on an accurate PES, which can be calculated with modern ab initio quantum chemistry software [46]. Naturally such a calculation will be multidimensional, and the efficiency of which should be judged against experimental results. The first test of the application of the QCM to a polyatomic molecule will be the observation and prediction of an unperturbed wavepacket. An important benchmarking of the QCM is also conceivable against time-dependent density-functional theory (TD-DFT) [61]- [63]. While the non-time-dependent formulation of this numerical technique has been successful in predicting the dynamics of large molecules, its applicability to strong-field ultrafast processes is still developing (see for example [64]); quasi-classical wavepacket dynamics and density functional theory could potentially be developed in parallel.
As has been apparent for a number of years, vibrational wavepacket propagation in a resonant external field often results in the coupling of a number of electronic states, making simulation more complex. Such systems require the consideration of conical intersections, which are the realm of few-cycle strong-field NIR and attosecond soft x-ray or XUV pulses [21]. The use of such pulses allows the temporal tracking of nuclear and electronic dynamics through such non-adiabatic transitions; however, this is at the expense of narrow-band selectivity. A coherent attosecond pulse can only be supported by an ultrabroadband photon field; as a result, a wide range of electronic and vibrational states will be populated following the pump pulse interaction.
We propose an interesting solution, or at least a viable approach: it is possible to exploit the varying coupling between bound and dissociative states under the influence of the control pulse to carry out exactly the state identification required by such investigations. It is highly unlikely that the energy difference between the available bound and dissociative electronic states would be identical over all internuclear separations. Therefore, by applying a control pulse of a known intensity, the different degree of PES distortion will vary the amount of upand down-shift of elements of the ensemble depending on the electronic state. By then Fourier transforming the probe delay-dependent fragmentation yield, the electronic states are potentially identifiable [25].
Clearly the optimum solution would be to balance selectivity of the photon energy with temporal resolution. Ultrashort (few to tens of fs) UV-XUV photon sources are currently being explored, either generated by monochromating the coherent radiation produced by highharmonic generation, for example the Astra Artemis Facility at the Rutherford Appleton Laboratory [65], or in free-electron lasers operating in the x-ray region, such as FLASH in Hamburg or LCLS in Stanford. It is under such conditions that new science will emerge and, along with more computationally demanding quantum mechanical methods, the QCM will allow the quantification of such new frontiers.

Conclusion
A QCM has been proposed that allows the quantification of wavepacket dynamics modified by an ultrashort strong-field non-resonant laser pulse. We have discussed how the vibrational phase and population are adjusted by the control pulse, and a comparison has been made with established theoretical predictions. Systematic predictions of wavepacket dynamics as a function of pump intensity and control delay and intensity have been presented. Such results will be of interest to groups attempting to experimentally detect the manipulation of a wavepacket. Finally, the application of the QCM to polyatomic systems is discussed, and the reader is referred to the wealth of compatible studies and techniques in coherent control. Attosecond studies in complex many-electron molecules is discussed, and a method for identifying the populated electronic states suggested, whereby the coupling of electronic states in a strongfield distorts the vibrational wavepackets generated, thus revealing the nuclear and electronic dynamics.