Mixed quantum-classical approach to multiphoton dissociation of the hydrogen molecular ion

We present a mixed quantum-classical propagation method for the time-dependent dynamics of nuclear as well as electronic degrees of freedom for para-H2+ exposed to short intense laser pulses of 800 nm wavelength. Depending on the initial vibrational state, the angular distributions of photofragments show characteristic shapes in very good agreement with our full-dimensional nuclear dynamics quantum calculations. The results can be understood in terms of two-dimensional adiabatic Floquet surfaces, which depend on the internuclear separation and rotation angle, demonstrating that adiabatic light-dressed surfaces are also a useful concept for full-dimensional nuclear dynamics. Using kinetic energy release spectra, we are also able to extract the contributions of different photon channels in the quantum-classical case.

3 barrier suppression and partial trapping of the vibrational wave packet, respectively-can be understood with the two lowest Floquet states, as was shown in the case of 266 nm photons [13]. In the present case of 800 nm and including rotational dynamics, more light-dressed surfaces are involved and they are now intrinsically 2D, depending in addition on the alignment angle. Nevertheless, we will see that the angular distributions can still be qualitatively understood by analyzing the most important dressed surfaces. In [14], angular distributions of fragments for low-lying vibrational states were analyzed on the basis of kinetic energy spectra, rather than on wave packet propagation on 2D field-dressed potential surfaces, for a laser of 532 nm wavelength and 50 TW cm −2 maximal peak intensity, where the adiabatic Floquet surfaces exhibit a simpler angular dependence than that in our case. A 1D analysis of multi-photon processes based on a Floquet-like propagation method has been performed in [15] for a laser of 800 nm wavelength.
Acknowledging how useful the Floquet picture is for understanding the dynamics qualitatively, we present a novel scheme for extracting multi-photon branching ratios from our mixed quantum-classical NA-QMD approach.
This paper is organized as follows. In section 2, the wave packet evolution by using the two methods is described. In section 3, we present vibrationally resolved angular distributions for the dissociation of photofragments, contrasting the results of the nuclear full quantum with those of the mixed quantum-classical calculation. Section 4 provides an interpretation of the angular distributions, identifying the three most important Floquet surfaces corresponding to 0, 1 and 2 photon absorption. This interpretation can also be derived consistently from the quantumclassical calculation based on the extracted multi-photon branching ratios. The conclusions and an outlook are given in section 5, followed by relevant technical details in the appendices.

Theory
In the center-of-mass frame and for dipole approximation in length gauge, the motion of the electron and the two nuclei of H + 2 in an intense laser field (t) is described by the time-dependent Schrödinger equation (TDSE) (atomic units (au) will be used unless stated otherwise): whereT R = −∇ 2 R /(2µ) is the kinetic energy of the protons with reduced mass µ being half the proton mass andĤ is the electronic Hamiltonian. The electric field with amplitude (t) = 0 f (t) cos(ωt) (3) is assumed to be linearly polarized along the z-axis with photon frequency ω and pulse envelope f (t) = sin 2 πt T , where T = 2T FWHM is the full pulse length. The laser parameters entering (3) used in the remainder of this paper are λ = 800 nm, corresponding tohω ≈ 1.57 eV, T = 25 fs and 2 0 ≡ I = 5.7 × 10 −3 au, which reads I = 2 × 10 14 W cm −2 in familiar units for the intensity. As we will compare the results of our full quantum calculations for the dissociation of H + 2 on laser-coupled Born-Oppenheimer (BO) surfaces with the results of mixed quantumclassical calculations using the NA-QMD method [3,6], let us first briefly recall the theoretical formulation of both approaches.

Full nuclear quantum dynamics
To obtain the full quantum solution and to extract the contribution of individual photon channels, we expand (R, r, t) in diabatic field-dressed Born-Oppenheimer (FBO) states where φ ξ (r; R) are the bare parametrically R-dependent eigenstates of the electronic Hamiltonian and ξ,n (R, t) the nuclear wavefunction on BO surface V ξ (R) dressed by photon number n. Using expansion (4), the TDSE reads where D ξ η (R) = d 3 r φ * ξ (r; R)zφ η (r; R) are the dipole matrix elements and L is the angular momentum. Equation (5) solved numerically by expanding the nuclear wavefunction ξ,n (R, t) in a basis of plane waves and spherical harmonics Y m l using the WavePacket code [16]. At the center of interest in this paper is the angular distribution of the photofragments P n D (θ) corresponding to the n-photon channel, which is given by where θ is the angle between the laser polarization and the molecular axis and t f is a final time well after the laser pulse is over; R D denotes the lower boundary of the dissociated region. We stress that these quantities are probability densities. To extract the total angular dissociation probability density P D (θ ), the different channel contributions have to be summed over. Furthermore, for the total dissociation probability P D , the resulting quantity has to be integrated over the angle using an appropriate sin(θ) weighting.

Mixed quantum-classical dynamics
At the heart of this description is a classical treatment of the nuclear motion, while explicitly describing the electron dynamics quantum mechanically in the framework of NA-QMD [3]. This implies a basis expansion for the electronic wavefunction, which is guided by the classically determined location of the nuclei in time. Hence, it is convenient to keep the coordinates R 1,2 of the two nuclei separately, instead of working as usual in the center-of-mass frame with the relative coordinate R = R 1 − R 2 of the nuclei. In order to solve the electronic Schrödinger equation (r denotes the electronic coordinate) depending parametrically on the nuclear coordinates, we expand the electronic wavefunction (r, t; R 1 , R 2 ) into a large set of basis functions centered at the two nuclei and a hexagonal grid R A α ≡ R kl (for details, see appendix A). The grid functions are necessary in order to include the continuum of ionized electrons realistically. Using equation (8) the electronic Schrödinger equation readṡ with the overlap matrix S αβ = φ α |φ β , the Hamilton matrix H αβ = φ α |Ĥ |φ β and the nonadiabatic coupling B αβ = φ α | d dt φ β . In order to avoid non-physical reflections of probability, an absorber potential matrix V (abs) αβ is introduced [6]. For this purpose we introduce field-following adiabatic states χ (a) as eigenstates of the electronic HamiltonianĤ at time t, After the expansion of χ (a) in the set of basis functions used in (8) as where the coefficients f (E (a) ) determine the strength of the absorber at eigenenergy E (a) and are given explicitly in (A.3) of appendix A. Self-consistently with the electronic Schrödinger equation (9), we solve the Newton equations of motion for the nuclei The correction term D A αβ is given by with To calculate the angular distribution of photofragments within the mixed quantum-classical method, we take into account only those dissociated trajectories, denoted by the index j D , that end up in the corresponding angular interval Each of those trajectories has to be weighted, on the one hand, according to its final norm N j D (t f ), i.e. the norm of the electronic wavefunction of the respective trajectory j D , which may decrease due to the absorber. On the other hand, the initial angular weight due to initial uniform 6 angular sampling has to be taken into account by a factor sin(θ j D (t 0 )). The angular distribution of photofragments is then given by where denotes a normalization constant and the factor sin(θ k ) −1 is included in order to guarantee that the total dissociation probability is given by summation over the final angle θ k as P D =

Results for angle-dependent photodissociation
We first present the central numerical results on the total dissociation probability density P D (θ) = n P n D (θ ) with P n D (θ ) from (6) as a function of the rotational angle. Computational details of the two different approaches that we have used are given in appendix A. For the laser parameters used (see above), our NA-QMD results reveal that the ionization probabilities of each vibrational level do not exceed 20%, i.e. the photofragmentation dynamics is dominated by dissociation. Hence, our full quantum results, allowing only the description of dissociation, can serve to test the accuracy of the NA-QMD approach.
The total angular distributions of fragments according to the dissociation channel, for the initial vibrational levels from ν = 1 → 14 under investigation, are shown in figure 1, as well as the corresponding Franck-Condon average. As can be seen clearly, the mixed quantum-classical and the full quantum results show good qualitative agreement. For ν = 0 the NA-QMD results (not shown) display negligible dissociation, a fact that has already been observed in the laseraligned case [7]. Whereas for ν 4 the probability densities decrease as a function of angle, for intermediate ν a double humped distribution can be observed, with a tendency for almost uniform distributions, most clearly visible at the highest vibrational excitation ν = 14. In the FC case, low ν are dominating the average, leading to a forward peaked distribution. We will explain the reasons for the nontrivial dependence on the rotational angle of the probability densities for the specific ν values in the next section.
Before doing so, to compare with already published results in the aligned case we show the total dissociation probabilities for the initial vibrational levels, comparing full quantum results with the NA-QMD ones in figure 2. Again, the almost perfect agreement between the two results up to ν = 10 is striking. The deviations above ν = 10 are probably related to the Ehrenfest (mean-field) treatment, which can lead to an overestimation of the dissociation probability. Furthermore, large probabilities of about 0.8 are observed around ν = 10, which is unprecedented in the aligned case (θ = 0, see dotted lines in figure 2) [7], where maximal dissociation probabilities were about 60%.

Interpretation of the results
In the following, we aim to extract the multi-photon effects in fragment angular distributions from mixed quantum-classical calculations including nuclear rotation. To this end, we first   angular-resolved dissociation probability from the most important few-photon channels in quantum and NA-QMD description. How we extract the different photon channel contributions from the NA-QMD results is explained in appendix B.
For low initial vibrational levels, the angular distributions are peaked in the forward direction. For these initial states (ν = 1-5), the photodissociation occurs mainly via the 2ωdissociation channel. As can be seen in figure 3(c), in this case the relevant Floquet surface shows a ridge (indicated by the solid line) between θ = 30 • and θ = 60 • for R 4.0 au. This leads to an alignment behavior of the dissociating wavepacket on the 2ω-Floquet surface, exemplarily shown in figure 4 for initial ν = 4. At the beginning of the laser pulse and also for higher angles during the laser pulse, the dynamics is mainly governed by the diabatic pathways of the Floquet surfaces for low ν, which results from the fact that the avoided crossings in these cases are very small. Only at low angles is the Floquet barrier of the 2ω-surface lowered enough at high intensities, allowing for a dissociation on this surface (figures 4(c) and (d)), leading to a forward peaked angular distribution. In figure 4(b), the results of the mixed quantum-classical calculations reproduce these findings well.
For higher vibrational levels, 1-photon dissociation plays a dominant role, maximally present for initial ν = 9 (see figures 3(b) and 5). In these cases, the initial wave packet already    Starting with ν = 9, 0-photon dissociation represents an additional influence on the angular distribution of photofragments. As can be seen in figures 1 and 6, a broadening of the angular distribution due to the 0-photon contribution occurs. This results from the fact that a relevant part of the wavepacket is located at internuclear distances, which are far away from the avoided crossing that forms between the 0ω-and 1ω-surface (see figures 3(a) and (b)). For high ν the initial wavepacket extends beyond the beginning of the ridge on the 0ω-surface and part of the wavepacket can slide to higher angles, leading to a slight anti-alignment behavior (finite dissociation contribution at θ ≈ 90 • ), as shown in figures 6(a)-(d). Again, the mixed quantumclassical calculations reproduce these alignment features.
Finally, we stress that our 2D Floquet interpretation is based not only on kinetic energy spectra, but also on propagation on about 20 FBO potential surfaces, allowing a time-resolved analysis of dissociation dynamics, shown in the density plots in figures 4(c) and (d), 5(c) and (d) and 6(c)-(h), which is impossible on the basis of kinetic energy spectra [14].

Conclusions and outlook
By comparing with the 'full' quantum results on nuclear dynamics on FBO surfaces, we have shown that a mixed quantum-classical NA-QMD treatment allows for a good description of angle-resolved dissociation probability densities, as well as of total dissociation probabilities of H + 2 driven by a fs laser pulse, as a function of initial vibrational excitation. Furthermore, we were able to resolve the different multi-photon contributions to the probability densities, both fully quantum mechanically using 2D Floquet BO surfaces, and in the quantum-classical NA-QMD approach, using an assignment of classical trajectories together with kinetic energy release spectra. Details of the latter methodology are presented in appendix B. Using these approaches, we were able to show that for different initial vibrational excitations of the molecular ion, different Floquet channels play a dominant role. Whereas for low initial ν the 2ω surface is needed for explaining the dissociation dynamics, for intermediate ν the 1ω surface and for high ν the 0ω surface play a decisive role.
Including the rotational degree of freedom enhances the total dissociation probabilities with respect to the aligned case for ν 6 [7]. The large dissociation probabilities we observe for high initial excitation hints at the fact that ionization is suppressed by the rotational motion in these cases. We will elaborate on this topic in a separate publication.
Finally, we stress that the NA-QMD method is also applicable to polyatomic many-electron systems [5].

Appendix A. Computational details
In this appendix, we give computational details that are used for the dressed state dynamics and motivate our choice of the basis in the NA-QMD case.

A.1. Dressed state quantum dynamics
In order to solve the coupled channel TDSE (5), we use 7σ g and 10σ u dressed state BO surfaces, which was checked to be sufficient to obtain converged angular distributions of photofragments. For the radial grid we use 2000 equally spaced grid points ranging from 0.2 to 200.0 au, whereas angular momenta up to l = 59 are used. The time step is chosen as 0.01 fs ≈ 0.4 au. We take the initial state for the propagation as a product of the vibrational eigenstate ν and the spherical harmonic corresponding to the rotational ground state, starting from the σ g + 0ω surface: The onset of dissociation has been arbitrarily set to R D = 10 au, a value close to what has been used previously by other authors [17].

A.2. Mixed quantum-classical dynamics
For mixed quantum-classical NA-QMD calculations, we use a basis set consisting of 125 basis functions in (8), where 33 Gaussian basis functions (for the functional form, see (14) in [18]) were centered at each nucleus (parameters are collected in table A.1) and 59 s-type Gaussians with width σ = 5.156 au, which are centered at R A α ≡ R kl on a hexagonal grid, defined by where the distance is d = 3.869 au and e 1,2 = ± 1 2 e y + √ 3 2 e z . The absorber function f in (11) needed for the description of ionization is chosen as using a typical lifetime of τ min = 4.76 au and a reference energy of E ref = 0.3 au. Further details of the absorber can be found in [6]. In order to solve the full NA-QMD equations (9) and (12) for the H + 2 problem, including electronic and nuclear motion, the initial nuclear separations and vibrational momenta are sampled according to a microcanonical distribution attributed to the vibrational level ν. The angle θ is sampled uniformly in the interval [0, π 2 ]. The initial angular momentum is set to zero. The time step for the classical motion (which takes place on a larger time scale than the electronic motion) is chosen as 1 au. In order to split the angular distribution of photofragments into multi-photon channels within the mixed quantum-classical method, where the nuclear motion is governed by the Ehrenfest (mean-field) force equation (12), the occupation probabilities of the bare BO, i.e. σ g and σ u , states are calculated for each trajectory. If the final population of the bare BO state differs by at most 10% from the total final population N j D (t f ), then the trajectory indexed by j D is regarded as of σ g -, respectively σ u -type. An exemplary visualization of this procedure is given in figure B.1, showing three different trajectories in the upper panel. Mean field (MF) here corresponds to a trajectory that cannot be assigned to one specific surface, whereas this was possible for the trajectories ending on the σ g , respectively σ u surface. The assignment of the photon number is done using the nuclear kinetic energy spectrum (KER), which can be extracted from NA-QMD analogous to (16)