Vibrational Excitation of H2 Scattering from Cu(111): Effects of Surface Temperature and of Allowing Energy Exchange with the Surface.

In scattering of H2 from Cu(111), vibrational excitation has so far defied an accurate theoretical description. To expose the causes of the large discrepancies with experiment, we investigate how the feature due to vibrational excitation (the "gain peak") in the simulated time-of-flight spectrum of (v = 1, j = 3) H2 scattering from Cu(111) depends on the surface temperature (Ts) and the possibility of energy exchange with surface phonons and electron-hole pairs (ehp's). Quasi-classical dynamics calculations are performed on the basis of accurate semiempirical density functionals for the interaction with H2 + Cu(111). The methods used include the quasi-classical trajectory method within the Born-Oppenheimer static surface model, the generalized Langevin oscillator (GLO) method incorporating energy transfer to surface phonons, the GLO + friction (GLO+F) method also incorporating energy exchange with ehp's, and ab initio molecular dynamics with electronic friction (AIMDEF). Of the quasi-classical methods tested, comparison with AIMDEF suggests that the GLO+F method is accurate enough to describe vibrational excitation as measured in the experiments. The GLO+F calculations also suggest that the promoting effect of raising Ts on the measured vibrational excitation is due to an electronically nonadiabatic mechanism. However, by itself, enabling energy exchange with the surface by modeling surface phonons and ehp's leads to reduced vibrational excitation, further decreasing the agreement with experiment. The simulated gain peak is quite sensitive to energy shifts in calculated vibrational excitation probabilities and to shifts in a specific experimental parameter (the chopper opening time). While the GLO+F calculations allow important qualitative conclusions, comparison to quantum dynamics results suggests that, with the quasi-classical way of describing nuclear motion and the present box quantization method for assigning the final vibrational state, the gain peak is not yet described with quantitative accuracy. Ways in which this problem might be resolved in the future are discussed.

Introduction

In scattering of a diatomic molecule from a metal surface, vibrational excitation may be intimately linked to the molecule’s dissociative chemisorption, as bond stretching is involved in both cases. Given the importance of elementary molecule–metal surface reactions to heterogeneous catalysis[1,2] and the observation that vibrationally inelastic scattering can probe the barrier region of reactive potential energy surfaces (PESs),[3−6] it is not surprising that vibrationally inelastic scattering of molecules from metal surfaces has become a subject of intense study. Systems that have been studied experimentally include H2 + Cu(111),[3,4,7] H2 + Cu(100),[8,9] NO + Au(111),[10,11] NO + Ag(111),[12] N2 + Pt(111),[13] HCl + Au(111),[14] and CO + Au(111).[15]

There is considerable evidence that vibrationally inelastic scattering of molecules other than H2 from metal surfaces is governed by an electronically nonadiabatic mechanism.[10−15] However, the H2 + Cu(111) system has often been viewed as a system in which vibrational excitation happens in a mostly adiabatic mechanism, in competition with dissociative chemisorption and resulting from the stretching of the molecule as it approaches its transition state.[3,4,16] Features have been identified in the PESs for H2 interacting with low-index Cu surfaces that are thought to promote vibrational excitation in an electronically adiabatic manner.[5,6,17] These considerations would seem to make vibrational excitation in H2 + Cu systems a phenomenon that should be straightforward to model and understand, but as will now be discussed, this is not true.

Vibrationally inelastic scattering of H2 from copper surfaces has been studied experimentally for H2 + Cu(111)[3,4,7] and H2 + Cu(100).[8,9] In a detailed study on H2 + Cu(111),[7] on which we will focus, a high-energy molecular beam was scattered from Cu(111) at an incidence angle that was slightly off-normal (θi = 15°), with the [12̅1] azimuth selected as the incidence plane. The amount of H2 molecules scattered to the (v = 1, j = 3) state (v is the quantum number for vibration and j the quantum number for rotation) was determined in a time-resolved manner, using resonance enhanced multiphoton ionization (REMPI) and time-of-flight (TOF) techniques. With the setup that was used, vibrational excitation from several (v = 0, j) states to (v = 1, j = 3) is evident from a peak occurring at short times in the time-of-flight spectrum (see Figure ). By reference to the TOF signal of a freely moving H2 beam, Rettner et al. called this peak the “gain peak”.[7] At longer times another broad peak was evident (Figure ), which reflects rotationally inelastic scattering within v = 1 as well as loss of (v = 1, j = 3) H2 due to dissociative chemisorption and vibrational de-excitation. This peak was therefore called the “loss peak”.[7] Using an equation similar to the one we will be using below in attempts to reproduce this experiment, Rettner et al. were able to extract quantitative information on vibrational excitation on the assumption that j should be conserved in the vibrational excitation process (i.e., vibrational excitation probabilities P(v = 0, j = 3 → v = 1, j = 3)).[7] Theoretical work later indicated that this assumption is not justified, because j is not conserved in vibrational excitation and because the total vibrational excitation probability to v = 1 depends rather strongly on the initial value of j in the initial (v = 0, j) state of H2.[18] Finally, experiments performed for surface temperatures (Ts) of 400 and 700 K suggested that raising Ts weakly promotes vibrational excitation[7] (see also Figure ), which has been attributed to a mechanism involving surface phonons.[16]

Figure 1

Time-of-flight spectrum of (v = 1, j = 3) H2 scattering from Cu(111) at a surface temperature (Ts) of 400 K (black dots) or 700 K (red dots) or in a freely traveling molecular beam (blue dots) at a position corresponding to the same flight length as traversed by H2 in the scattering experiment. In all cases the data points were connected by lines, which merely serve to guide the eye. The data were taken from refs (7) and (16).

Prior to 2009, attempts to describe vibrationally inelastic scattering of H2 from Cu(111) were severely hampered by the accuracy in the potential energy surfaces (PESs) used in the dynamics calculations. By this, we mean that discrepancies between calculations and experimental observations could always be blamed on inaccuracies in the PES used in the calculations. This changed to a large extent when in 2009 a chemically accurate PES became available for H2 + Cu(111), from a semiempirical implementation of density functional theory (DFT).[19] With this PES, sticking probabilities of H2 and D2 on Cu(111),[19] the influence of the initial vibrational and rotational states of H2[19] and D2[20] on reaction on Cu(111), and rotational excitation of H2 scattering from Cu(111)[19] could all be described with chemical accuracy. The PES also allowed an accurate description of the rotational quadrupole alignment parameter of D2 desorbing from Cu(111) in two rovibrational states.[21] The good performance of the PES for a variety of reactive scattering experiments (in addition to rotationally inelastic scattering) suggests the PES (or the accompanying density functional) should also be good for studying vibrationally inelastic scattering, which occurs in competition with reaction and in the same energy regime.[7] Nevertheless, quantum dynamics (QD) calculations carried out within the Born–Oppenheimer static surface approximation (BOSS model) and using this PES underestimated the gain peak in the TOF spectrum of Figure by about a factor of 3. By this, we mean that the calculated vibrational excitation probabilities had to be multiplied by a factor of 3 to reproduce the TOF spectrum.[16,22]

Further analysis suggested that this failure should be primarily due to the failure of the dynamical model (i.e., the BOSS model) rather than to the PES used.[16] Specifically, the analysis showed that if seemingly plausible assumptions were made about how vibrational excitation should change from the hypothetical Ts effectively used in the theory (0 K) to its experimental value (400 K), and about the size of energy loss to the surface, the discrepancy between theory and experiment could be reduced to a factor of 2.[16] It was suggested that the absence of phonons in the dynamical model should be primarily responsible for the discrepancy between theory and experiment.[16] Another suggestion was to examine this further with ab initio molecular dynamics (AIMD) calculations, if possible with the method extended in such a way that electron–hole pair (ehp) excitation could also be modeled with an electronic friction approach to examine its role. The expectation was formulated[16] that the quasi-classical treatment of nuclear motion should not represent a severe limitation, as the quasi-classical trajectory (QCT) method should already be reasonably accurate for describing vibrational excitation at the collision energy (80 kJ/mol)[23] at which the contribution of vibrational excitation to the gain peak in Figure peaks.

Meanwhile, it has become possible to perform AIMD calculations with an electronic friction description of ehp excitation (AIMD with electronic friction, or AIMDEF).[24−27] However, the method is still quite expensive, especially if the goal is to obtain scattering probabilities with high statistical accuracy for a large range of incidence energies and initial states, as required for the accurate simulation of the measured TOF spectrum[7] (the gain peak in Figure ). Here, we will use the AIMDEF method to benchmark a computationally much cheaper to use method incorporating the effects of phonons and electronic friction, i.e., a generalized Langevin oscillator method incorporating electronic friction (GLO+F).[28,29] We will show that, compared to AIMDEF, the GLO+F method accurately describes vibrational excitation of H2 up to and including the incidence energy (Ei) most relevant to the simulation of the experiment, while exhibiting still reasonable accuracy for higher Ei.

The goals of the present work are as follows: We will explore whether the QCT method used in GLO, GLO+F, and AIMDEF calculations is capable of yielding quantitatively accurate results for the simulated TOF spectrum exhibiting vibrationally inelastic scattering, through comparison with quantum simulations. Next, we will use the QCT and the GLO+F methods to explore whether previous speculation of how vibrational excitation probabilities should depend on the incidence angle[16] was correct. This is useful knowledge as the need for performing quantum calculations for off-normal incidence would make a QD approach much more computationally expensive. Next, GLO and GLO+F calculations are used to explore how introducing surface phonon motion and ehp excitation into the dynamical model affects calculated vibrational excitation probabilities and whether their inclusion improves the agreement between the simulated and experimental TOF spectra, as speculated earlier.[16] We will also investigate whether the promotion of vibrational excitation through increased surface temperature[7] is due to heating the surface phonons, as assumed earlier,[16] or to heating the metal electrons. The calculations will reveal that the measured TOF spectrum is highly sensitive to how quickly the vibrational excitation probabilities rise above their threshold. We will therefore also perform some simulations to establish how shifting quantum dynamically calculated vibrational excitation probabilities along the energy axis, and uncertainties in the origin of time in the experiments (i.e., the beam chopper opening time), affect the computed TOF spectrum and its comparison to experiment.

This paper is set up as follows: Section describes the methods used, their implementation, and numerical details. Specifically, section describes how to simulate the experimental TOF spectrum exhibiting vibrational excitation. Sections , 2.3, 2.4, and 2.5 describe the QCT, the GLO, the GLO+F, and the AIMDEF methods used. Section discusses how the molecule–surface interaction and, from these, the forces are obtained in the present work. Section discusses several details of the implementation of the methods, such as operational definitions of scattering probabilities, the generation of initial conditions, the running of trajectories, and numerical details. Section presents the results of benchmarking the GLO+F method against the AIMDEF method. Section presents QCT, GLO, and GLO+F calculations for normal incidence, also benchmarking the QCT method against QD. Section investigates the effect of the incidence angle. Like the preceding section, this section also discusses the effects of introducing phonons and electronically nonadiabatic effects in the dynamics calculations, and which of these effects promote vibrational excitation if the surface temperature is raised. Section presents the results on how shifting vibrational excitation probabilities along the energy axis, and the time-origin in the experiments, affect the simulated TOF spectrum. Section discusses how the theoretical description of vibrationally inelastic scattering of H2 from Cu(111) (or from metal surfaces in general) might be improved in the future. Conclusions are presented in section .

Method

Modeling the Experiment

In the experiments we simulate (see Figure 1 of ref (7)), a chopped molecular beam of H2 molecules travels toward a Cu(111) crystal and is scattered from it. The incident beam makes a polar angle θi = 15° with the surface normal, and the incidence is along the [12̅1] azimuth.[7] The amount of molecules scattered from the surface in the (v′ = 1, j′ =3) state is measured in a time-resolved manner with REMPI and can be described according to[7,16]In eq , t measures the time from the opening of the chopper and vi and vs are the velocities of the incident and scattered molecules, respectively, with both depending on the initial (v, j) state of H2 and t through energy conservation, as described in ref (7). However, in some cases we take into account that the scattered molecule incurs a loss of a fraction fl of the kinetic energy that would be available to it if it were to scatter from a static surface in an electronically adiabatic manner:Here, m is the mass of H2, Ei its incident translational energy, and E its internal rovibrational energy depending on v and j. Furthermore, v0 is the stream velocity of the molecular beam (4115 m/s), and α(vi) a width parameter that can take on one of two values depending on vi (1358 or 2379 m/s for vi < v0 or vi > v0, respectively).[7] Also, N is a normalization factor, and c defines an offset; xi and xs define the distance traveled by the incident and scattered molecules, and xt defines their sum (values are given in ref (7)). The Boltzmann population of the initial (v, j) state in the beam divided by the Boltzmann population of the initial (v = 1, j = 3) state in the beam is given by the weight factor w. For the nearly effusive beam used in the experiments, these populations may be calculated assuming that both the rotational and the vibrational temperatures were equal to the nozzle temperature Tn of 2000 K.[7]

The use of a nearly effusive beam with a high value of Tn (i.e., with a broad energy distribution) in combination with the chopper and the detecting laser allowed the experimentalists to obtain vibrational excitation probabilities for very high Ei, which are not accessible in ordinary supersonic molecular beam experiments on H2. Although, at the time of writing, the experiments were done almost 25 years ago, they are quite well documented and yet to be surpassed in accuracy and information content when it comes to experiments on vibrational excitation of H2 scattering from metal surfaces.

For the simulation of the TOF signal described by eq , the crucial inputs from dynamics calculations are the state-to-state probabilities P(v, j → v′ = 1, j′ = 3) and the fractional energy losses fl (eq ), which may be taken to be dependent on v and j. Upon scattering, fractional energy losses may occur that reflect energy loss to phonons (in GLO, GLO+F, and AIMDEF calculations), to ehp’s (in GLO+F and AIMDEF calculations), and to translational motion perpendicular to the scattering plane (in all calculations performed for off-normal incidence).

Quasi-Classical Trajectory (QCT) Method

In the calculations with the QCT method,[30] the momenta and positions of the H atoms labeled by i and j are evolved according to Hamilton’s equations of motionIn eq , V(r, r) is a six-dimensional potential energy surface describing the interaction of the two H atoms at positions r and r with one another and with the static Cu(111) surface, and in eq , the velocities of the atoms are represented by v. The QCT method differs from the ordinary classical trajectory method in that in the QCT method zero-point vibrational energy (possibly added to extra vibrational energy if the molecule is vibrationally excited initially) is always imparted to the molecule at the start of the trajectories.

Generalized Langevin Oscillator (GLO) Method

In the GLO method,[31−33]eq is rewritten asIn eq , rs is the position of the surface atom nearest to H2, which is taken to move as a three-dimensional harmonic oscillator with mass ms (here taken as the mass of a surface Cu atom). The momentum of the surface atom obeyswhere Ω2 is the diagonal 3 × 3 frequency matrix associated with the surface harmonic oscillator and Λgs is a diagonal 3 × 3 matrix that couples the motion of the surface to a ghost oscillator with position rg. In turn, the momentum of the ghost oscillator is given byThe third term on the right-hand side (rhs) of eq models energy dissipation from the surface to the bulk of copper through a friction coefficient ηph, which is computed fromwhere ωD is the Debye frequency of the solid.[32] The randomly fluctuating force Rph is modeled as Gaussian white noise with variance[34]In eq , Δt is the time step used in the integration of the equations of motion and kB is the Boltzmann constant. The linking of the randomly fluctuating force with the phonon friction coefficient ensures that the fluctuation–dissipation theorem is obeyed,[34] so that thermal equilibrium can be restored after the direct scattering event (i.e., in the present case of H2 + Cu(111)). The elements of the diagonal matrix Ω2 are equal to 2ω2, and the elements of the diagonal coupling matrix Λgs are equal to ω2, with ω being the surface phonon frequencies (i = x, y, z).

Generalized Langevin Oscillator with Electronic Friction (GLO+F) Method

In the GLO+F method,[28,34]eq is extended further toIn eq , the effect of energy transfer involving ehp’s is modeled with molecular dynamics with electronic friction (MDEF)[35] using the local density friction approximation (LDFA)[36] in the independent atom approximation (IAA).[36] The friction coefficients used in the LDFA have been successfully applied to calculate the stopping power of atoms and ions by metal solids and surfaces[37−40] and in the modeling of scattering of H atoms from Au surfaces.[41] When using the IAA to apply the method to molecules, the assumption is made that the electronic friction is independent of the electronic structure of the molecule, so that electronic friction forces can be specified through atomic friction coefficients ηel,, as done in eq . In the LDFA, the atomic electronic friction coefficients depend on the electronic density of the bare metal surface at the position of the atom relative to the surface.[36] The LDFA-IAA method has now been used to study the effect of electronic excitations on the dynamics of molecules scattering from metal surfaces in several applications,[28,29,36,42,43] including the H2 + Cu(111) system.[44] In eq , the randomly fluctuating force Rel represents the nonadiabatic scattering of thermal surface electrons from the molecule. To ultimately enable descriptions in which the molecule becomes equilibrated to the surface, the fluctuation–dissipation theorem is taken into account[45] by modeling this force as Gaussian white noise with variance[34]

Ab Initio Molecular Dynamics with Electronic Friction (AIMDEF) Method

In the calculations with the AIMDEF method,[24−26] essentially quasi-classical calculations are carried out for the nuclear dynamics, with the forces calculated on the fly from DFT. Of course, we also model ehp excitation, by adding electronic friction forces and a randomly fluctuating force (second and third terms on the rhs of eq ) to the adiabatic forces on the H atoms in the simulations. The motion of not just the impinging H2 molecule, but also the Cu atoms in the upper layers of the Cu(111) slab is simulated. The Cu(111) slab is thermalized prior to the scattering calculations at the experimental Ts. Since the scattering and reaction of H2 on Cu(111) occur in a direct manner (without the molecule performing several bounces on the surface), there is no need to thermalize the atoms in the layers in which they are allowed to move while the collision proceeds. Therefore, we simply have one layer of stationary Cu atoms at the bottom of the slab in the simulations we carry out, and the GLO formalism is not applied to the motion of these atoms.

Molecule–Surface Interaction

In the QCT, GLO, and GLO+F calculations, the original specific reaction parameter (SRP) PES for H2 + Cu(111) was used. The SRP functional used effectively to generate the PES[19,22] is a weighted average of the revised Perdew–Burke–Ernzerhof (RPBE) functional[46] (mixing coefficient 0.43) and the Perdew–Wang 1991 (PW91) functional[47] (mixing coefficient 0.57). Further details of its calculation can be found in refs (19) and (22), and elbow plots of two-dimensional cuts through the PES are shown in Figure S1 of ref (19). The SRP barrier geometry and height are provided in Table for two geometries, i.e., the minimum-barrier geometry in which H2 impacts on a bridge site, with the H atoms moving to hollow sites (bth), and a geometry in which H2 impacts on a top site, with the H atoms moving to bridge sites (ttb). The first geometry is most relevant for reaction at low Ei, while the second geometry is thought to be most important to vibrational excitation, as the features of the elbow cut are thought to be conducive to vibrationally inelastic scattering (reaction path with large curvature in front of an especially late barrier;[5,6,17] see also Figure S1E of ref (19)).

Table 1

SRP Minimum-Barrier Geometry and the SRP and SRP48 Values of the Molecule–Surface Interaction Energy E at These Geometriesa

geometry	db (bohr)	Zb (bohr)	E (eV), SRP	E (eV), SRP48
bridge-to-hollow	1.95	2.20	0.628	0.628
top-to-bridge	2.64	2.62	0.891	0.876

In all cases, H2 is parallel to the surface. db is the H–H distance and Zb the molecule–surface distance at the barrier.

For reasons related to the need for being able to work with a variable (i.e., Ts-dependent) lattice constant, which are discussed in detail in ref (21) and its accompanying Supporting Information, the SRP functional had to be reparametrized to enable AIMD (and, here, AIMDEF) calculations.[21] With the reparametrized (i.e., SRP48) functional (which is a weighted average of the RPBE functional[46] (mixing coefficient 0.48) and the PBE functional[48] (mixing coefficient 0.52)), the molecule–surface interaction at the SRP minimum-barrier bth geometry is reproduced to within better than 1 meV, while the molecule–surface interaction at the SRP ttb barrier geometry is reproduced to within 15 meV (see ref (21) and Table ). Therefore, the use of a density functional in the AIMDEF simulations somewhat different from the one implicit in the use of the SRP PES in the QCT, GLO, and GLO+F calculations should, taken by itself, only lead to small discrepancies between the calculated vibrational excitation probabilities.

Numerical Details and Implementation

Operational Definition of Reaction and of Rovibrationally Inelastic Scattering

In all calculations a similar operational definition is used for reaction and for rovibrationally inelastic scattering. Reaction is defined to occur once the H–H distance in a trajectory becomes larger than 1.6 Å in the AIMDEF calculations (2.2 Å in all other calculations). Scattering is defined to occur once the molecule–surface distance becomes larger than 6.1 Å, with the velocity of the molecule pointing away from the surface in the AIMDEF calculations (9 Å in all other calculations).

The assignment of the final rovibrational state is done as follows: The classical analogue of the rotational quantum number is computed aswhere jc is the classical rotational angular momentum. Next, the rotational quantum state j is assigned by binning jq to the nearest odd value of j, keeping in mind the conservation of parity in H2 and our interest in scattering to an odd j state within v = 1 (i.e., (v = 1, j = 3);[7] see also section ). The vibrational state v is then assigned by computing the classical rovibrational energy of the molecule and comparing it to the quantum mechanical vibrational energies within the j ladder and assigning v to describe the (v, j) state with the nearest rovibrational energy in that ladder. Probabilities (whether for reaction or rovibrationally inelastic scattering) are simply computed by dividing the number of trajectories resulting in the outcome of interest by the total number of trajectories. In the AIMDEF calculations, a total of 1100 trajectories were run for each Ei. Much larger numbers of trajectories were computed in the QCT, GLO, and GLO+F calculations, i.e., 200000 trajectories for each v, j state at each Ei.

Initial Conditions

H2 was initialized with is center of mass 6 Å away from the surface in the AIMDEF calculations (9 Å in all other calculations), with a velocity directed toward the surface according to the Ei simulated. Depending on the simulation, either the incidence direction was normal to the surface or the incidence angles were taken equal to the experimental values (see below). A Monte Carlo integration was performed over randomly selected impact points on the surface and initial orientation angles and directions of rotational velocities, in accordance with the initial value of j and the initial magnetic rotational quantum number m (the projection of j on the surface normal) as described in, for instance, ref (49). In the AIMDEF calculations, a uniform sampling of m was performed by running equal numbers of trajectories for −j ≤ m ≤ j. In all other calculations, instead a random sampling was performed of the orientation of the classical angular momentum for each j. Initial values of the H–H distance d and its conjugate momentum were selected by performing a uniform sampling in time of these values along a one-dimensional quasi-classical trajectory run for isolated H2 for the quantum mechanical energy computed for the relevant initial (v, j) state with the Fourier-grid Hamiltonian method.[50]

In the GLO and GLO+F calculations, the initial position of the surface atom, rs, and its conjugate momenta are sampled through a conventional Monte Carlo procedure in such a way that they correspond to the experimental surface temperature. In the AIMDEF simulations, the initial coordinates and velocities are sampled from pre-equilibrated four-layer Cu slabs, in which the atom positions and velocities in the upper three layers are representative of a Cu(111) surface at the experimental Ts. The procedure used has been described in ref (21). The AIMDEF calculations used a value of the surface lattice constant that corresponds to a bulk lattice constant of 3.698 Å, based on the calculated SRP48 value of the static lattice constant of bulk copper and the experimentally determined thermal expansion of copper at 400 K[51,52] as described in ref (21). Initial positions and velocities were sampled from 1000 different snapshots from 10 equilibrated surfaces (10000 snapshots in total). The 10 surfaces, from which coordinates and velocities were sampled, are characterized by an average surface temperature of 399.5 K. The distribution of surface temperatures characterizing the 10 surfaces exhibited a standard deviation of 60 K.

Trajectory Calculations

In all methods used, the equations of motion were solved with the Beeman algorithm[53] as implemented in refs (33) and (54). This has the advantage that the coordinates and velocities are available to high accuracy at the same points in time. This is relevant to, for instance, the accurate calculation of the electronic friction forces, which require the coordinates and velocities to be available to high accuracy at the same point in time (because friction forces are based on atomic velocities and positions, where the latter determine the friction coefficients).

The AIMDEF calculations were performed with a user-modified version of the ab initio total energy and molecular dynamics program VASP (version 5.4) developed at the Universität Wien.[55,56] In the DFT calculation of the forces, the ultrasoft pseudopotentials[56,57] were used in combination, with which the SRP48 functional was parametrized for H2 + Cu(111).[21] The AIMDEF calculations on H2 + Cu(111) used a time step of 0.25 fs. All other details of the AIMD and DFT calculations are the same as described in ref (21).

Other Numerical Details

In the GLO calculations, the surface phonon frequencies ω were taken equal to 14 meV for i = x, y, z, where 14 meV is equal to the surface Debye frequency of Cu(111),[58] as used before in refs (22) and (58). The same value of the (surface) Debye frequency was used in eq to calculate the phonon friction coefficient ηph (see section ).

In the GLO+F calculations, the electronic density of the bare Cu(111) surface, which is needed to compute the friction coefficients, was calculated from a three-dimensional cubic spline interpolation using a previously prepared spline fit. In the AIMDEF calculations, the electronic density was calculated instead from the self-consistent density of the entire H2 + Cu(111) system with displaced Cu atoms, with subtraction of the densities due to the H atoms using a Hirshfeld partioning scheme,[26,27,59] which was also used in ref (60). Electronic friction coefficients for H atoms were obtained in the usual way[37,38,61] by computing the phase shifts of Kohn–Sham orbitals at the Fermi momentum for a proton embedded in a free electron gas for different values of the electronic embedding densities. The friction coefficients were parametrized throughusing a fit expression used earlier in ref (25), and employing atomic units. eq accurately describes the friction coefficient for values of the free electron radius r ranging from 1 to 10 bohrs, which covers the range relevant to our calculations.

Results and Discussion

Comparison of the GLO+F and AIMDEF Methods

Vibrational excitation probabilities P(v = 0, j = 5 → v = 1, j = 3) computed with the AIMDEF method and the GLO+F method for H2 + Cu(111) and Ts = 400 K are compared in Figure and Table . As is the case for all other AIMDEF results shown in this paper, the results are obtained for the conditions relevant to the experiments (i.e., the quoted value of Ts and θi = 15° with incidence along the [12̅1] azimuth[7]), and j = 5 was selected because the initial (v = 0, j = 5) state makes the largest contribution (see Figure 3 of ref (16)) to the gain peak (see Figure ) in the TOF spectrum. The GLO+F results reproduce the AIMDEF results rather closely for the incidence energies Ei ≤ 0.83 eV most important to describing the gain peak. For flight times corresponding to Ei > 0.83 eV, the “blue” (high-energy, short-time) tail of the gain peak drops off rather quickly, due to the exponentially decreasing amount of molecules present in the incident beam with these high Ei values.[7] The somewhat diminished accuracy with which the GLO+F method describes vibrational excitation at these Ei values should therefore be less relevant, and we conclude that the GLO+F method may justifiably be used to explore the effect of a range of factors on the computed TOF spectrum. The diminished accuracy of the GLO+F method for higher E might be due to this method being less capable of describing the more elaborate surface deformation that could occur at such energies. AIMDEF is intrinsically capable of describing surface deformation involving more than one atom, whereas the GLO+F method is not.

Figure 2

P(v = 0, j = 5 → v = 1, j = 3) as computed with the QCT method, with the GLO+F method for 400 and 700 K, and with AIMDEF for 400 K. The error bars on the AIMDEF results denote 68% confidence intervals.

Table 2

AIMDEF and GLO+F Results for Off-Normal Incidence at Ei = 0.829 eV and Ts = 400 K

	AIMDEF	GLO+F
P(v = 0, j = 5 → v = 1, j = 3)	0.037 ± 0.006	0.0354
Et of (v′ = 1, j′ = 3) H2 (eV)	0.417 ± 0.019	0.429
Etp of (v′ = 1,j′ = 3) H2 (eV)	0.399 ± 0.019	0.407

The energy lost by scattered H2 is relevant to the interpretation of the experiments on vibrational excitation to (v = 1, j = 3),[7] because it determines the velocity with which H2 travels through the detection zone,[16] where H2 is laser-excited with REMPI. The final translational energy Et of (v′ = 1, j′ = 3) H2 obtained upon scattering of (v = 0, j = 5) H2 at Ei = 0.829 eV is shown in Table , comparing AIMDEF and GLO+F results for Ts = 400 K. As can be seen, the AIMDEF and GLO+F results for Et are in excellent agreement with one another for this Ei. The same is true for the final translational energy in the scattering plane, Etp (see also Table ). In Table , we compare not only the calculated values of Etp, but also the standard deviations associated with the distributions of these energies for two different values of Ei. The GLO+F method correctly predicts not only the average Etp (Tables and 3), but also the standard deviations of the distributions of these energies (Table ), suggesting that the GLO+F method is capable of correctly predicting the distributions of the final translational energies in scattering and therefore also the loss of energy to the surface phonons and to ehp’s.

Table 3

AIMDEF and GLO+F Results Compared for Off-Normal Incidence and Scattering from (v = 0, j = 5) to (v = 1, j = 3) at the Values of Ei Indicated and Ts = 400 Ka

	AIMDEF	GLO+F
Etp (eV) at Ei = 0.6 eV	0.33 (0.07)	0.33 (0.07)
Etp (eV) at Ei = 1.05 eV	0.48 (0.17)	0.49 (0.14)

The first number presented is the average Etp value, and the second number (in parentheses) is the standard deviation of the distribution of Etp (see the text for its definition).

Dissociative chemisorption probabilities computed with the AIMDEF and GLO+F methods for (v = 0, j = 5) H2 + Cu(111) and Ts = 400 K are compared in Figure . As can be seen, compared to AIMDEF, the computationally much less expensive GLO+F method yields quite accurate results for the reaction. Including the effects of phonon motion as well as ehp excitation leads to a much larger broadening of the H2 reaction probability curve (relative to the QCT static surface result) than obtained when QCT static surface results are compared to AIMD results for D2 + Cu(111),[20] where the effects of phonons should actually be larger for the heavier D2. This is important, as the latter broadening was observed to be much smaller[20] than measured experimentally[62−64] for D2 + Cu(111). The above results suggest that ehp excitation should be taken into account for a correct description of the broadening effect of increasing Ts on the reaction probability curve for H2/D2 + Cu(111) and that this effect can be obtained just as well with GLO+F as with AIMDEF.

Figure 3

Reaction probability as computed for off-normal incidence with the QCT method, the GLO+F method for 400 K, and AIMDEF for 400 K for (v = 0, j = 5) H2 + Cu(111). The error bars on the AIMDEF results denote 68% confidence intervals.

Finally, probabilities P(v = 0, j = 5 → v′, j′) are shown for several v′ and j′ states for Ei = 0.829 eV in Figure . As can be seen, the GLO+F results are in quite good agreement with the AIMDEF results for this Ei. Together, Figures and 4 help to understand why the GLO+F method is quite accurate for computing state-to-state probabilities for scattering to the (v = 1, j = 3) state: the GLO+F method is also quite capable of describing the competition of reaction and of scattering to other rovibrational states with scattering to (v = 1, j = 3).

Figure 4

P(v = 0, j = 5 → v′, j′) as a function of j′ for E = 0.829 eV as computed with the GLO+F and AIMDEF methods for 400 K and for v′ = 1 (upper panel) and v′ = 0 (lower panel). The error bars on the AIMDEF results denote 68% confidence intervals.

TOF Spectra Simulated with Scattering Calculations Performed for Normal Incidence

To investigate the reliability of the quasi-classical approximation for computing TOF spectra for comparison with the experiments on vibrational excitation,[7] the TOF spectrum computed with the QCT method is compared to that computed with QD in Figure . To account for the fact that the experiments were done for off-normal incidence, in the simulation of the TOF spectra, the assumption was made that the vibrational excitation probabilities only depend on the total translational energy, and not on the incidence angle (total energy scaling, or TES). This approximation allows the dynamics calculations to be performed for normal incidence only, which was favorable for the earlier QD, time-dependent wave packet (TDWP) calculations.[16] While the TDWP method allows results to be obtained for a range of incidence energies in just one wave packet propagation,[65] this only applies to a calculation with a fixed initial parallel momentum,[66] making the simulation of TOF spectra based on QD calculations performed for off-normal incidence expensive.

Figure 5

TOF spectra simulated from calculations performed for normal incidence with the TDWP method (QD), the QCT method, and the GLO+F method for Ts = 400 and 700 K, assuming total energy scaling (TES) of vibrational excitation. No energy loss to the surface was taken into account. The black dots denote the experimental TOF spectrum at Ts = 400 K.[7]

We first note that the QCT method quite well reproduces the long-time loss peak in the TOF spectrum calculated with QD. This suggests that the QCT method correctly accounts for the loss of intensity in the TOF spectrum of (v = 1, j = 3) H2 due to reaction and vibrational de-excitation from v = 1 to v = 0 and that it correctly describes rotationally inelastic scattering within v = 1. However, the QCT calculations do not reproduce the QD short-time gain peak, which reflects vibrational excitation. Figure shows that this is due to the QCT method overestimating P(v = 0, j → v = 1, j = 3) at low Ei for j = 1, 3, and 5. As a result, the height of the gain peak increases considerably at long flight times (low Ei). While this brings the peak height into better agreement with experiment, the peak position is shifted to too long flight times (to too low Ei).

Figure 6

P(v = 0, j → v = 1, j = 3) as computed for normal incidence with the TDWP method (QD), the QCT method, and the GLO+F method for Ts = 400 and 700 K and for j = 1 (lower panel), 3 (middle panel), and 5 (upper panel).

The level of agreement achieved between QCT and QD for the P(v = 0, j → v = 1, j = 3) is disappointing in view of the expectation voiced in ref (16) that the vibrational excitation probabilities should be reasonably well described with QCT for Ei ≥ 0.8 eV. This expectation was based on calculations performed on H2 + Cu(110) using a PES calculated with the PW91 functional.[23] Being calculated with the PW91 functional, this PES most likely[19] contains a too low barrier for dissociation. Already near Ei = 0.8 eV, the vibrational excitation probability P(v = 0, j = 0 → v = 1, j = 0) computed in ref (23) took on much higher values (close to 0.1) than observed here for vibrational excitation to (v = 1, j = 3) (see Figure ), making it easier to reproduce QD results with the QCT method (the QCT method generally performing better for larger probabilities). The rather poor performance of the QCT method for calculating probabilities P(v = 0, j → v = 1, j = 3) for low Ei for our system and PES represents a setback, as it impairs our capability of accurately simulating the “red” (i.e., low-Ei, long-time) side of the gain peak in the TOF spectrum. However, we may still hope that the effects of the incidence angle, and of allowing energy transfer to surface phonons and ehp’s, are reasonably well described with quasi-classical mechanics, taking the QCT result for normal incidence as the reference. In doing so, we keep in mind that the QCT calculations (and most likely also the GLO+F calculations) are likely to overestimate the gain peak at its low-energy side. In the future, it should be worthwhile to investigate whether better results can be obtained with more sophisticated binning methods for assigning the final vibrational state[67] than used here in the QCT calculations.

To investigate the role of Ts and of allowing energy exchange with the surface, the TOF spectra computed with the GLO+F method are compared with the spectrum computed with the QCT method in Figure for Ts = 400 and 700 K. It is gratifying to see that the GLO+F results qualitatively reproduce the experimental finding[7] that the height of the gain peak increases with Ts increasing from 400 to 700 K (see also Figure ). This finding can be explained by the P(v = 0, j → v = 1, j = 3) calculated with GLO+F being higher at low Ei for Ts = 700 K than for Ts = 400 K for j = 1, 3, and 5, especially for j = 1 and 3 (see Figure ).

However, it is disturbing to see that the gain peak obtained with the GLO+F method for Ts = 400 K is much lower than that computed with the QCT method, which uses the static surface approximation. In ref (16), the effect of Ts was estimated by assuming that the static surface results obtained with TDWP calculations should approximately equal the results for a 0 K surface in which the surface atoms are allowed to move. To estimate how results for a 400 K surface should differ from results obtained with the static surface approximation, the assumption was made that the Ts dependence of the vibrational excitation probabilities, as inferred from the differences in the observed peak height for 400 and 700 K, could be extrapolated from 400 to 0 K. More specifically, the observation of a 20% increase in the gain peak on going from Ts = 400 K to Ts = 700 K was interpreted as evidence that the gain peak should rise by 27% on going from the static surface (0 K) result to the 400 K result.[16] The results of Figure show the opposite trend: the peak height is decreased on going from the static surface QCT result to the 400 K GLO+F result. This result suggests a rather limited potential of surface motion to account for the observed discrepancies between the experimental and previous theoretical results for vibrational excitation of H2 scattering from Cu(111). A caveat here is that the result obtained of course depends on at which Ei the QCT and GLO+F vibrational excitation probabilities cross. We have assumed that this crossing point is reliably obtained by taking the QCT result as the standard to measure the GLO+F results against for taking into account the effects of Ts and of allowing surface motion. Since we know that the QCT method does not accurately reproduce the QD results at low Ei, these assumptions may not be entirely correct and require testing in future work.

Figure shows the effect of taking into account that the (v′ = 1, j′ = 3) H2 molecules lose energy to the surface phonons and through ehp excitation, which ensures that the molecules fly less quickly through the detecting laser beam, making detection more likely, in the GLO+F calculations. For Ts = 400 K and Ei = 0.829 eV, depending on j, energy losses were in the range of 18–26% of the available energy (Ei – Eth), where Eth is the threshold Ei for vibrational excitation to the (v = 1, j = 3) state of H2 in the BOSS model (see also eq ). This ensures that for the GLO+F calculations the gain peaks in Figure are somewhat higher than in Figure . Whereas the energy loss percentages calculated with GLO+F are somewhat smaller than assumed in a previous analysis (i.e., 30%),[16] the results of Figure show that allowing energy loss to “surface modes” such as phonons and ehp’s does lead to a modest increase in the height of the gain peak attributed to vibrational excitation. This partly accounts for the differences previously observed between QD and experimental results for vibrational excitation of H2 scattering from Cu(111) (see also below for the results for off-normal incidence).

Figure 7

TOF spectra simulated from calculations performed for normal incidence with the TDWP method (QD), the QCT method, and the GLO+F method for Ts = 400 and 700 K. Energy loss to the surface was taken into account in the GLO+F results. The black dots denote the experimental TOF spectrum at Ts = 400 K.[7]

TOF Spectra Simulated with Scattering Calculations Performed for Off-Normal Incidence

The experiments were performed for θi = 15° with incidence along the [12̅1] azimuth.[7] However, previous QD calculations simulated TOF spectra from vibrational excitation probabilities computed for normal incidence. This required assumptions about whether at off-normal incidence the vibrational excitation probabilities should scale with the normal or total incidence energy or whether an intermediate scaling would be obeyed.[16] To investigate whether the assumptions made were correct, here we also performed QCT calculations for off-normal incidence for the experimental conditions. Figure shows that the gain peak obtained in this way is much closer in height to the result obtained from QCT normal incidence results assuming normal energy scaling (NES) than that obtained assuming total energy scaling (TES).

Figure 8

TOF spectra calculated with the QCT method for off-normal incidence (15° off-normal, QCT OFF) and with the QCT method from normal incidence results assuming total energy scaling (QCT TES) or normal energy scaling (QCT NES). The black dots denote the experimental TOF spectrum at Ts = 400 K.[7]

Figure shows that at low Ei the P(v = 0, j → v = 1, j = 3) values computed for off-normal incidence from normal incidence results assuming TES are larger than the P(v = 0, j → v = 1, j = 3) values computed directly for off-normal incidence (OFF) for j = 1–5. This explains why the gain peak computed from normal incidence results assuming TES is too high compared to the QCT result obtained directly for the actual incidence conditions. A previous quantum dynamical result obtained directly for the experimental off-normal incidence condition suggested that, in the range of Ei of interest, the vibrational excitation probabilities should fall midway between the probabilities computed from normal incidence results assuming TES and NES (“intermediate scaling”).[16] However, the off-normal incidence result was obtained for only one value of Ei (0.83 eV) and for one value of j (j = 3). The present work suggests that the assumption about the scaling being intermediate between TES and NES does not necessarily hold and that at low Ei it might be better to assume NES than intermediate scaling as done earlier.[16] Relative to the earlier work, this should lead to increased discrepancies between theory and experiment for the gain peak; i.e., the computed gain peak, which was already too small, should be further reduced by a factor of 1.25.[16] The present work also suggests that it should be better to obtain QD results for off-normal incidence. This could be done by performing TDWP calculations for several values of the initial momentum parallel to the experimental [12̅1] incidence plane[7] and interpolating these results to get results for the range of Ei of interest to the experiments. With present day computational resources, this should now be feasible.

Figure 9

P(v = 0, j → v = 1, j = 3) for off-normal incidence as computed with the QCT method directly (QCT OFF) and from normal incidence QCT results assuming total energy scaling (QCT TES) or assuming normal energy scaling (QCT NES) for Ts = 400 K and for j = 1 (lower panel), 3 (middle panel), and 5 (upper panel).

The TOF spectrum obtained from off-normal incidence results computed with the GLO+F method for Ts = 400 K is compared to the QCT static surface result and the GLO+F spectrum computed for 700 K in Figure . In the computation of the gain peak in the TOF spectrum based on the GLO + F results, all energy losses out of the translational motion in the incident plane were taken into account, now also including loss to translational motion out of the scattering plane. Depending on the initial value of j, energy losses ranged from 21% to 36% (see Table ), in reasonable agreement with the value of 30% assumed in earlier work.[16] Energy losses to translation outside the scattering plane contribute 5–14% to these percentages. Note that we have verified that the fractional energy loss to phonons, ehp’s, and motion out of the scattering plane is, to a reasonable extent, independent of Ei.

Figure 10

TOF spectra calculated from results for scattering at off-normal incidence with the TDWP method (QD TES, from normal incidence results assuming TES), the QCT method (QCT OFF), and the GLO+F method for Ts = 400 and 700 K (400 K OAL and 700 K OAL). All energy losses (to the surface and to motion out of the scattering plane) were taken into account in the GLO+F results. The black dots denote the experimental TOF spectrum at Ts = 400 K.[7]

Table 4

Energy Losses for Off-Normal Incidence GLO+F Calculations for Ts = 400 and 700 K and Ei = 0.829 eVa

j	Ei – Eth	EQCT	EGLO+F, 400 K	energy loss (%)	EGLO+F,700 K	energy loss (%)
1	0.2447	0.402	0.351 (0.364)	21 (16)	0.371	13
3	0.3175	0.445	0.368 (0.395)	24 (16)	0.388	18
5	0.4456	0.518	0.407 (0.446)	25 (16)	0.421	22
7	0.6246	0.669	0.480 (0.545)	30 (20)	0.487	29
9	0.8461	0.875	0.588 (0.689)	34 (22)	0.595	33
11	1.0471	1.115	0.743 (0.857)	36 (25)	0.740	36

EQCT is the total final translational energy obtained with QCT calculations. EGLO+F, 400 K is the translational energy for motion in the scattering plane as obtained with GLO+F for Ts = 400 K, and EGLO+F,700 K is the same obtained for Ts = 700 K. The energy loss (%) is obtained by taking the difference of EGLO+F and EQCT and dividing by (Ei – Eth). All energies are in electronvolts. Values in parentheses are for GLO calculations for 400 K.

As also found for normal incidence (Figure ), the GLO+F gain peak for Ts = 400 K is lower than the QCT gain peak, even though energy losses are taken into account in the GLO+F calculations but not in the QCT calculation of the TOF peak. The reason is the same as for normal incidence: the GLO+F vibrational excitation probabilities are smaller than the QCT probabilities for the initial j values that are important to the calculation of the gain peak and the low Ei values relevant to the gain peak (see Figure ).

Figure 11

P(v = 0, j → v = 1, j = 3) for off-normal incidence as computed directly with the QCT method (QCT) and with the GLO+F method (left) or the GLO method (right) for Ts = 400 and 700 K and for j = 1 (lower panels), 3 (middle panels), and 5 (upper panels).

It is of interest to see whether modeling ehp excitation (using the GLO+F rather than the GLO method) leads to an increase in the gain peak in the TOF spectrum, which one would then normally attribute to increased vibrational excitation. As can be seen from Figure , this is not the case: for the value of Ts at which most experiments were done (400 K), modeling ehp excitation leads to a lower gain peak. The reason is that for the important initial j states (j ≤ 5) and the Ei values most relevant to the calculation of the gain peak (≤0.9 eV) the computed vibrational excitation probabilities decrease if friction is introduced (see Figure ). This effect is more important than the extra translational energy loss to friction, which leads to increased detection because the vibrationally excited H2 flies more slowly through the detection zone (see Table , noting that the translational energy loss is larger in GLO+F than in GLO for Ts = 400 K). Apparently, the effect that the molecule has already lost some energy to ehp excitation when it hits the surface leads to decreased vibrational excitation at the lower Ei, and the ensuing effect on the gain peak is larger than the increased detection following from the energy loss to ehp excitation.

Figure 12

TOF spectra calculated from the results for scattering at off-normal incidence with the TDWP method (QD TES, from normal incidence results assuming TES), the GLO+F method (400 K F), and the GLO method (400 K NF) for Ts = 400 K. Energy loss to the surface and to motion out of the scattering plane was taken into account in the GLO and GLO+F results. The black dots denote the experimental TOF spectrum at Ts = 400 K.[7]

Figure 13

P(v = 0, j → v = 1, j = 3) for off-normal incidence as computed directly with the GLO+F and GLO methods for Ts = 400 K and for j = 1 (lower panel), 3 (middle panel), and 5 (upper panel).

One might also ask what the effect on the gain peak (and on vibrational excitation) is of raising Ts and whether this depends on whether ehp excitation is modeled. The answer to this question is rather surprising. Even though the effect of ehp excitation is to decrease the gain peak and vibrational excitation at lower Ei for Ts kept fixed, raising Ts only leads to a considerable increase in the height of the gain peak and to increased vibrational excitation if ehp excitation is included! This can be seen by comparison of Figure (showing TOF spectra computed with the GLO method) with Figure (showing TOF spectra computed with the GLO+F method). Just heating the phonons from 400 to 700 K hardly raises the gain peak (Figure ), whereas heating the phonons and the electrons does (Figure ). This GLO+F effect of raising Ts does not come from more efficient detection due to greater loss of translational energy; in fact, the percentages of energy loss to the surface are smaller for 700 K (Table ). Instead, the effect arises from the P(v = 0, j → v = 1, j = 3) calculated with the GLO+F method being larger for 700 K than for 400 K for low j and the relevant low Ei, especially for j = 1 and 3 (see Figure ). In contrast, the P(v = 0, j → v = 1, j = 3) values calculated for low initial j with the GLO method are more or less the same for 400 and 700 K (see also Figure ). We attribute the effect of raising Ts on the vibrational excitation probabilities computed with GLO+F to the concomitant greater random force exerted by the electrons on the H nuclei in the region of configuration space where the electronic friction coefficients are large (see eq ). The rise in the gain peak seen in the GLO+F calculations (Figure ) is in qualitative agreement with experimental observations[7] (see Figure ), and our analysis suggests that the experimentally observed increase of the gain peak with Ts can be attributed to an electronically nonadiabatic mechanism. However, it would be good to check in future research whether a phononic contribution to the rise of the gain peak with Ts would result from AIMD calculations, which treat the phonons in a more sophisticated way. Such calculations would probably need to employ a much larger number of AIMD trajectories than used here, to make the statistical error bars small enough to enable the detection of potential, reasonably sized phononic contributions to the rise of the gain peak with Ts.

Figure 14

TOF spectra calculated from the results for scattering at off-normal incidence with the TDWP method (QD TES, from normal incidence results assuming TES) and with the GLO method for Ts = 400 and 700 K (400 K NF and 700 K NF). Energy loss to the surface and to motion out of the scattering plane was taken into account in the GLO results. The black dots denote the experimental TOF spectrum at Ts = 400 K.[7]

Effects of Uncertainties in Computed Vibrational Excitation Probabilities and Time Origin in Experiment on the TOF Spectrum

In section we saw that uncertainties in computed vibrational excitation probabilities (in this case, due to the use of the classical approximation in the QCT calculations) can lead to large changes in the computed TOF spectrum. For this reason, we decided to explore two factors that might affect the TOF spectrum.

The first effect we explored concerns the calculations. It is conceivable that, due to errors in the PES (for instance, in the reaction barriers), the vibrational excitation probabilities should be shifted to lower incident translational energies. To explore this, we took the original quantum dynamical vibrational excitation probabilities (see, e.g., Figure ) and shifted them to lower energies by 1 kcal/mol (∼43 meV). Figure shows how this alters the gain peak in the TOF spectrum relative to the original QD results, in both cases assuming that at off-normal incidence the vibrational excitation probabilities obey TES. The shifts along the energy axis lead to a substantial increase in the height of the gain peak (Figure ), which is due to the strong dependence of the TOF signal on the incident velocity distribution (see eq ). If additionally an energy loss of 30% is assumed (as done originally in ref (16) and justified to a large extent by our GLO+F results; see section and Table ), the peak is further increased (Figure ). If additionally the gain peak is multiplied by a factor of only 1.5, already quite good agreement is obtained with experiment for the height of the gain peak (see also Figure ). To be fair, it should be noted that our present calculations suggest that it should be better to assume NES than TES and that this should change[16] the multiplication factor required for good agreement of the peak height with experiment to a factor of 2.2.

Figure 15

TOF spectra calculated with the TDWP method on the basis of the original QD vibrational excitation probabilities (QD TES), on the basis of the QD vibrational excitation probabilities shifted along the incident energy axis by −1 kcal/mol (QD SHF), additionally assuming 30% energy loss (f(K) = 0.3), and additionally multiplying the shifted probabilities by a factor of 1.5 (fm = 1.5). The black dots denote the experimental TOF spectrum at Ts = 400 K.[7]

The gain peak obtained by shifting the QD vibrational excitation probabilities by −1 kcal/mol is not yet in the right position; i.e., it occurs at a too low energy (too long time of flight), as may be seen from Figure . However, the experiments also contain an uncertainty of about 1 μs, which is due to a lack of perfect timing in comparing reference and scattered TOF distributions.[7] To see the possible effect of such a shift, we additionally shifted the time origin by −1 μs. The effect of this additional shift is shown in Figure . As may be seen, the position of the peak is now also in the approximately correct place. In other words, shifting the vibrational excitation probabilities by −1 kcal/mol in incident translational energy and shifting the time origin by a value within the experimental error (i.e., by 1 μs) leads to a TOF spectrum with the gain peak occurring at almost the correct flight time (corresponding Ei). However, the peak still needs to be multiplied by a factor of 2.2 (1.5), assuming vibrational excitation at off-normal incidence obeys normal (total) energy scaling. Nevertheless, these TOF spectrum simulations show that the position of the gain peak in the TOF spectrum is quite sensitive to both the dependence of the vibrational excitation probabilities on incidence energy and the origin of time in the measurements. The comparison of theory to experiment would clearly benefit from both a higher time resolution in the experiments and, possibly, improved accuracy of the description of the dependence of the computed vibrational excitation probabilities on the incident translational energy.

Figure 16

TOF spectra based on the original QD vibrational excitation probabilities (QD TES), on QD vibrational excitation probabilities shifted along the incident energy axis by −1 kcal/mol (QD SHF), additionally applying a shift of the time origin by −1 μs (QD SHFT), and additionally assuming 30% energy loss and multiplying the shifted probabilities by a factor of 1.5 (fm = 1.5). The black dots denote the experimental TOF spectrum at Ts = 400 K.[7]

Discussion: Remaining Issues and Future Improvements

The comparison between GLO+F and AIMDEF results (Figure ) suggests that, for modeling the TOF spectrum exhibiting vibrational excitation in the gain peak, it should be enough to use the GLO+F method and model the surface through a single oscillating atom connected to a ghost atom. In turn, this suggests that it might be possible to model the effect of (allowing) surface motion on vibrationally inelastic scattering of H2 from Cu(111) with QD calculations in which the motion of only one surface atom connected to a bath of oscillators and ehp’s is considered. This can perhaps be done in the spirit of the GLO method, for instance using a density matrix formalism[68−70] or using stochastic wave function approaches.[71−74] Before such QD calculations are carried out, it might be useful to test whether it is necessary to retain all three degrees of freedom of the surface atom or whether it might be enough just to model its motion perpendicular to the surface. Finally, the GLO+F method clearly suffices for modeling the effect of surface phonons and ehp excitation on dissociative chemisorption (Figure ) and to calculate energy loss to the surface (Tables and 3).

The present comparison between QCT and TDWP results clearly shows that QCT with box quantization of the final vibrational state, as employed here, is not of sufficient accuracy to reliably compute the gain peak representing the effect of vibrational excitation in the TOF spectrum (Figure ). Future work should test whether better results can be obtained if more sophisticated methods, such as Gaussian binning techniques,[67] are used to assign final vibrational states to the outcome of quasi-classical trajectories. Alternatively, it might be possible to analyze the final vibrational state of the scattered molecule using a so-called adiabatic switching procedure.[75] This obviously represents an important issue, as the present work suggests that the vibrational excitation probabilities computed with methods using the quasi-classical approximation (including QCT, GLO+F, and AIMDEF) are much too high at the incidence energies relevant to the simulation of the gain peak (Figure ). To enable a quantitatively accurate simulation of the TOF spectrum, this problem should be solved, or one should revert to the use of QD and solve the problem of how to model the effect of energy transfer involving surface phonons and ehp’s with a quantum dynamical approach.

Our present QCT results for off-normal incidence strongly suggest that, within a QD approach, it should not be an accurate approximation to assume that vibrational excitation obeys NES or TES or a scaling midway between these two limits characterized by a single mixing coefficient for all times of flight (Figures and 9). Rather, the calculations suggest that, in a QD approach, the calculations should be performed for off-normal incidence. This will make the QD simulation of the TOF spectrum much more expensive, as TDWP calculations will have to be performed for several values of the initial translational momentum in the plane of incidence. However, with present-day computational resources, this should now be possible, at least within the BOSS model.

Our TOF simulations based on QD vibrational excitation probabilities also show that the simulated TOF spectrum is quite sensitive to the origin of time in the experiments (Figure ) and to the exact dependence of the vibrational excitation probabilities on Ei (Figures and 16). Concerning the latter, it should be possible to remove some remaining uncertainties regarding the PES by reparametrizing an SRP functional for H2 + Cu(111), using a correlation functional approximately describing the attractive van der Waals interaction with the surface,[76−79] as done successfully for CH4 + Ni(111).[80] There are indications that such a functional should exhibit a somewhat higher, and later, barrier to reaction.[81] Both the higher barrier and the somewhat increased velocity toward the surface resulting from the slight acceleration in the van der Waals well (with an experimental depth of about 30 meV[82]) could lead to increased vibrational excitation and to energy shifts of the vibrational excitation probabilities, on which the simulated TOF spectra show a strong dependence (Figure ). The use of such a reparametrized functional might also lead to an even more accurate description of the orientational dependence of reaction of D2 on Cu(111) than the one already obtained with the SRP48 functional.[21] Additionally, with such a functional, one could also test whether the experimentally measured effect of selective adsorption resonances on scattering of H2 from Cu(111) in the low-incidence-energy regime (<50 meV)[83] could be accurately described. Calculations using correlation functionals that describe the van der Waals interaction in at least an approximate way suggest that this should be possible,[84,85] and an empirical potential describing scattering in the van der Waals energy regime and reproducing the resonances is available from potential inversion.[83] Furthermore, an investigation[81] of H2 + Cu(111) that considered a few experiments on H2 + Cu(111) also addressed with the SRP and SRP48 functionals showed that these experiments were equally well described with the optimized Perdew–Burke–Ernzerhof van der Waals density functional (optPBE-vdW-DF).[86] This latter functional was not used in the present study as it has not yet been used for quantitative comparison with the same wide range of experiments as the SRP and SRP48 functionals.

Concerning new calculations, it might also be of interest to investigate how the use of tensorial friction coefficients[87,88] might alter the results compared to the present use of the LDFA-IAA method. Note, however, that tensorial frictions should be calculated in the exact quasi-static limit, which to our knowledge has not been accomplished yet.[89]

Finally, it would certainly be advantageous if additional experiments were to become available for comparison. For instance, it would already be helpful if, in experiments similar to the one we now use for comparison, the origin of time would be much better defined than the present value (of 1 μs),[7] as the TOF spectrum is quite sensitive to this value (Figure ). As already pointed out in ref (16) it would also be advantageous to know the nozzle temperature used in the original molecular beam experiment with high accuracy: the original work stated its value (2000 K), but not its uncertainty.[7] TOF spectra simulated on the basis of TDWP calculations suggested a strong dependence of the gain peak on the experimental nozzle temperature (see Figure 3 of ref (16)). It would also be nice to have TOF spectra available for a large range of surface temperatures. This might help to determine through dynamics calculations whether the observed dependence of vibrationally inelastic scattering on this parameter can indeed be explained with an electronically nonadiabatic mechanism, as our calculations now suggest. Finally, it would be even better to be able to compare directly to individually measured state-to-state vibrational excitation probabilities P(v = 0, j → v = 1, j′). Then an eventual agreement between theory and experiment could no longer be due to error cancellation between values obtained for different j values, as might have occurred in the present and earlier[16] work for the simulated TOF spectrum.

Conclusions and Outlook

Previous work has shown that, with the SRP-DFT functionals now available for H2 + Cu(111), the sticking of H2 and D2, the dependence of associative desorption on the final rovibrational (v, j) state, rotationally inelastic scattering, and even the orientational dependence of reaction can all be described quite accurately for this system. In contrast, vibrationally inelastic scattering of H2 from Cu(111) has so far defied an accurate theoretical description. Previous work on vibrational excitation in this system had raised several questions regarding its dependence on the incidence angle, on the surface temperature, and on allowing energy exchange with the surface and regarding the applicability of classical mechanics. Here we have tackled these questions using the QCT, GLO, GLO+F, and AIMDEF methods, employing the SRP functionals available for H2 + Cu(111) to eliminate uncertainties due to possible inaccuracies in the molecule–surface interaction as much as possible.

The results of the present work strongly suggest that, to model the feature in experimental TOF spectra due to vibrational excitation (the so-called gain peak; see Figure ) with reasonable accuracy, it should not be necessary to use AIMDEF to describe surface deformation on impact. Rather, it should suffice to treat energy exchange with and dissipation to the surface within a GLO+F formalism. The GLO+F calculations also accurately describe dissociative chemisorption and the competition with scattering to other rovibrational states than probed in the experiments we have used for comparison. Importantly to the simulation of the TOF spectra exhibiting the effects of vibrational excitation, compared with AIMDEF, the GLO+F calculations accurately describe energy loss to the surface, which is relevant to the detected TOF signal. The GLO+F results for energy loss (about 30%) are in good agreement with assumptions made earlier[16] about the size of the energy loss to the surface accompanying vibrational excitation of H2 in scattering from Cu(111).

The present comparison between QCT and QD results suggests that the QCT method accurately describes the so-called loss peak in the TOF spectra, which reflects reaction and vibrational de-excitation and rotationally inelastic scattering within the vibrationally excited state. Unfortunately, the feature in the TOF spectra most relevant to this work (the gain peak due to vibrational excitation) is not described with quantitative accuracy using quasi-classical mechanics, because the QCT method overestimates the vibrational excitation probabilities for the relevant initial rotational states and (low) Ei values. However, it is possible to derive a number of important qualitative conclusions from the GLO and GLO+F results by adopting the QCT results as a reference, keeping in mind that some of these conclusions may require further checks through QD calculations, as discussed above in several places.

The GLO+F calculations reproduce the experimental finding that raising Ts from 400 to 700 K promotes vibrational excitation. Surprisingly, the comparison with GLO results for these temperatures suggests that the effect is due to an electronically nonadiabatic mechanism, in which the randomly fluctuating forces due to the hotter electrons promote vibrational excitation at the higher Ts. Earlier work had assumed that the effect should be due to a mechanism involving surface phonons.[16] The fluctuating forces referred to are based on friction coefficients, suggesting that, through comparison with experiments on nonreactive or weakly reactive systems exhibiting strongly increasing vibrational excitation with increasing Ts,[13−15] calculations might be able to test different friction models aiming to describe the effects of ehp excitation.

Importantly, the present work shows that at moderate Ts (400 K) the effect of allowing energy transfer to the surface phonons (as evident from the GLO calculations) and to ehp’s (as evident from the comparison of GLO to GLO+F calculations) is to reduce vibrational excitation. Thus, on a 0 K surface, the mobility of the surface atoms and the possibility of ehp excitation ensure that the vibrational excitation is reduced relative to that found with the Born–Oppenheimer static surface (BOSS) model. This highlights the disagreement already found earlier between theoretical work and experiments on vibrational excitation, where the earlier work assumed that allowing energy transfer from the surface phonons should lead to increased vibrational excitation for the experimental 400 K surface. However, this work assumed that vibrationally inelastic scattering from a 0 K surface should essentially be equal to that occurring over a hypothetical static surface. The present work strongly suggests that this should not be the case.

The present work also suggests that, at off-normal incidence, vibrational excitation probabilities cannot be accurately computed from normal incidence calculations assuming total or normal energy scaling of these probabilities, nor can one assume a constant (i.e., independent of Ei) mixing of TES and NES to obtain results for off-normal incidence. At the Ei most relevant to the simulation of the gain peak, the scaling of the vibrational excitation probabilities is closest to NES, and this further highlights the quantitative disagreement found earlier between theoretical work and the TOF experiments on vibrational excitation.

Simulations performed on the basis of earlier QD calculations of vibrational excitation probabilities show that the gain peak is quite sensitive to a constant energy shift of calculated excitation probabilities and to the exact chopper opening time in the experiments. Applying shifts that are reasonable (by −1 kcal/mol for the vibrational excitation probabilities and −1 μs for the chopper opening time) leads to considerably improved agreement of the theory with the experiment.

Given the comparative ease with which energy transfer to the surface can be modeled with classical mechanics, we recommend that future research be directed at more accurate schemes for assigning final vibrational states in quasi-classical calculations.[67,75] It is also worthwhile to test QD schemes for incorporating energy transfer between the molecule and the surface phonons and ehp’s modeled as a bath.[68−74] Finally, additional experiments, in particular aimed at obtaining fully rotationally resolved state-to-state vibrational excitation probabilities, could be quite helpful for improving the theory. The present TOF experiments on vibrational excitation reflect scattering from several initial j states within v = 0 as well as energy losses to the surface, making the attribution of the causes of the disagreement between theory and experiment somewhat muddled.