Femtosecond-laser induced dynamics of CO on Ru(0001): Deep insights from a hot-electron friction model including surface motion

A Langevin model accounting for all six molecular degrees of freedom is applied to femtosecond-laser induced, hot-electron driven dynamics of Ru(0001)(2 × 2):CO. In our molecular dynamics with electronic friction approach, a recently developed potential energy surface based on gradient-corrected density functional theory accounting for van der Waals interactions is adopted. Electronic friction due to the coupling of molecular degrees of freedom to electron-hole pairs in the metal are included via a local density friction approximation, and surface phonons by a generalized Langevin oscillator model. The action of ultrashort laser pulses enters through a substrate-mediated, hot-electron mechanism via a time-dependent electronic temperature (derived from a two-temperature model), causing random forces acting on the molecule. The model is applied to laser induced lateral diffusion of CO on the surface, “hot adsorbate” formation, and laser induced desorption. Reaction probabilities are strongly enhanced compared to purely thermal processes, both for diffusion and desorption. Reaction yields depend in a characteristic (nonlinear) fashion on the applied laser fluence, as well as branching ratios for various reaction channels. Computed two-pulse correlation traces for desorption and other indicators suggest that aside from electron-hole pairs, phonons play a non-negligible role for laser induced dynamics in this system, acting on a surprisingly short time scale. Our simulations on precomputed potentials allow for good statistics and the treatment of long-time dynamics (300 ps), giving insight into this system which hitherto has not been reached. We find generally good agreement with experimental data where available and make predictions in addition. A recently proposed laser induced population of physisorbed precursor states could not be observed with the present low-coverage model.


I. INTRODUCTION
Femtosecond-laser (FL) induced reactions of molecules at metal surfaces are interesting for applications but also for fundamental reasons [1,2]. Applications range from energy research over photochemistry to photocatalysis. An example for a fundamental aspect of interest is the demonstration that FL pulses act in a much more subtle way than simply heating the substrate. This leads often to larger reaction cross sections, new reactions and reaction pathways, large isotope effects, and product energy distributions different from those found in thermal reactions [1][2][3][4][5].
To gain insight into the molecular details of the FL induced photoreactions at surfaces, theoretical modeling is indispensable. A fully quantum dynamical approach which also accounts for hot-electron mediated excitation and electronic damping would of course be most valuable. However, this goal has been achieved, for metal surfaces at least, only for reduced-dimensional models using open-system density matrix theory or stochastic wave-packet methods [15]. An alternative approach makes use of classical trajectories and either nonadiabatic surface hopping [3,16], or ground-state Langevin dynamics in which electronic friction and hotelectron excitation are included via a damping term and stochastic forces, respectively. The (quasi)classical approach allows for the treatment of multidimensional dynamics. In its most simple variant, the electronic damping can be derived from the substrate electron density in which the adsorbate atoms are embedded [the so-called local density friction approximation (LDFA) [17]], and stochastic forces arising from a time-dependent electronic temperature which is determined from the two-temperature model (2TM) [18]. This methodology has been successfully used elsewhere to describe FL induced desorption of diatomic molecules from metal surfaces [19][20][21]. In Ref. [21], for example, all six molecular degrees of freedom of H 2 and D 2 desorbing from Ru(0001) were treated in this way, providing a detailed microscopic picture of the hot-electron mediated dynamics.
The same methodology, called molecular dynamics with electronic friction (MDEF) for short, will be applied here for FL induced dynamics of CO on Ru(0001). Previous work based on ab initio molecular dynamics (AIMD) has suggested that for this system, hot substrate (phonon) motion cannot be neglected to explain, e.g., details of FL induced desorption (and oxidation) [13]. Therefore, we extend our previously adopted MDEF methodology to also account for surface motion, which will be achieved by using the generalized Langevin oscillator (GLO) model [22][23][24]. The combined method, called MDEF-GLO in what follows, includes electronic friction, hot-electron driven processes, surface motion, and, importantly, all (six) molecular degrees of freedom of CO relative to the surface as dynamical variables. To account for the latter, we make use of a recently developed ground-state potential energy surface (PES) which is based on gradient-corrected density functional theory. Note that a similar multidimensional, DFT-based MDEF-GLO approach has recently been used to model FL induced desorption and dissociation of O 2 at a Ag(110) surface [25]. Using a precomputed PES and friction coefficients allows us to follow many trajectories over long times, achieving this way much better statistics and time-resolved insight than by, e.g., costly (and fixed-temperature) AIMD calculations as done previously [13]. Armed with these methods, the following questions shall be addressed in this paper: (i) What kind of reactions/molecular motions will be induced by femtosecond-laser stimulation of CO@Ru(0001)?
(iii) How do the individual processes behave with respect to the variation of pulse characteristics, notably the laser fluence?
(iv) Are our findings consistent with previous experiments and their current interpretation? Or are new findings/interpretations suggested by our numerical simulations?
To answer these questions, in the next Sec. II we first of all describe methods, models, and the potential function in some detail. Results will be presented in Sec. III. A final Sec. IV summarizes and concludes this work.

A. Laser-driven molecular dynamics with electronic friction: MDEF
The Langevin equation provides a useful frame for classical dynamics corrected by frictional and fluctuation forces [19]. Specifically, here we consider the laser-driven dynamics of a CO molecule on Ru(0001), solving a Langevin equation of the following form: m k d 2 r k dt 2 = −∇ k V (r 1 ,r 2 ) − η el,k (r k ) dr k dt + R el,k (t).
Here, k = 1,2 for the two moving atoms, C (k = 1) and O (k = 2), with r 1 = r C , r 2 = r O being atomic coordinates of the moving atoms, and m 1 = m C and m 2 = m O are carbon and oxygen atom masses, respectively.
The first term on the right-hand side of Eq. (1) is a conservative force acting on atom k and arising from the potential V (r 1 ,r 2 ). The latter is based on a six-dimensional (6D) parametrization of a PES as obtained from periodic density functional theory within the generalized gradient approximation (GGA), including van der Waals corrections (see below). Further, the second term on the right-hand side accounts for the effect of electronic excitations via a friction force −η el,k (r k ) dr k /dt. Here, η el,k (r k ) is the friction coefficient for atom k, calculated within the local density friction approximation (LDFA) [17]. In this model, η el,k is obtained in terms of the scattering of electrons by an atom inside a free-electron gas (FEG). It is given by (l + 1) sin 2 [δ l (ε f ,r s ) − δ l+1 (ε f ,r s )]. (2) In this equation, the Wigner-Seitz radius r s = ( 3 4πρ ) 1/3 is a measure of the electron density ρ of the FEG. δ l is the scattering phase shift at the Fermi level corresponding to the potential induced by the static atomic impurity in the FEG, which is calculated within DFT [26][27][28]. Within the LDFA at each point of the trajectory ρ is chosen as the electronic density of the bare surface at the position of atom k. Here, it is calculated on a grid from GGA-derived electron densities (RPBE exchange-correlation functional [29]) of the bare Ru(0001) surface, the latter modeled by a three-layer slab. For the MDEF simulations, the friction coefficients are fitted to the following analytic expression: The numerical parameters for C (k = C) and O (k = O) are given in Appendix A. In Eq. (3), we consider electron densities ρ only above 4 × 10 −4 a −3 0 (for C) and 7 × 10 −4 a −3 0 (for O), respectively. At very low densities (large r s at large Z values), we assume the friction coefficients to be zero.
Finally, at finite (electron) temperatures a random force arises in Eq. (1), R el,k (t), which is calculated as Gaussian white noise with the properties R el,k = 0 and Equation (4) results from the second fluctuation-dissipation theorem [30], which relates frictional and random forces, in this case through the electron temperature T el (t). k B is Boltzmann's constant. Details of the numerical realization of Eq. (4) are given in Ref. [21]. The time-dependent electronic temperature T el (t) was calculated with the two-temperature model [18] (2TM), where one solves two coupled equations for electron and phonon temperatures T el (t) and T ph (t): C el and C ph are electron and phonon heat capacities, κ the electron thermal conductivity, z the perpendicular position relative to the surface, g the electron-phonon coupling constant, and S(z,t) the absorbed laser power per unit volume. S(z,t) depends on shape, center wavelength (λ), and intensity (or fluence) of the applied pulse (see, e.g., Refs. [6,21]). The 2TM uses a number of material constants, and we employ those of Ref. [21]. An exception is the optical penetration depth ξ of ruthenium, for which we considered in addition to 800 nm (as in Ref. [6]), also the case with λ = 400 nm, which was also used in Ref. [7]. In this case, the penetration depth as one of the material parameters is ξ = 6.9 nm, rather than ξ = 16.2 nm for λ = 800 nm [6,7]. For completeness, we list all parameters entering the 2TM in Appendix A.
The method just outlined is the molecular dynamics with electronic friction (MDEF) model as mentioned above.

B. Substrate phonons: MDEF-GLO model
To account for surface phonons in an approximate way, the GLO model [22][23][24] is adopted as in our previous work [25,31]. Accordingly, surface motion is described in terms of a three-dimensional (3D) harmonic oscillator of mass m s with coordinates r s and associated diagonal 3 × 3 frequency matrix . Dissipation and thermal fluctuations are modeled with the help of a ghost 3D oscillator with coordinates r g , which is subject to friction and random fluctuation forces that model the energy exchange between the surface atoms and the bulk thermal bath. The mass and the associated frequency matrix for the ghost oscillator are the same as for the surface oscillator, and both are coupled by a matrix gs . The equations of motion for surface and ghost oscillators are then m s d 2 r s dt 2 = −∇ s V (r 1 − r s ,r 2 − r s ) − m s 2 r s + m s gs r g and m s d 2 r g respectively. Here, the friction force −η ph dr g dt models energy dissipation from the interacting surface atoms to the bulk thermal bath. The random force R ph models the heating of the surface atoms due to the thermal motion of the bulk atoms. The random force R ph is Gaussian white noise calculated similarly as above for the electronic case. In particular, Eq. (4) is used, however, with η el and T el (t) replaced by η ph and T ph (t), respectively. T ph is the surface (phonon) temperature which is calculated in the 2TM model. System-specific parametrizations of the friction coefficient for phonons η ph , frequency matrix , coupling matrix gs , and some further details are given in Appendix A. When combined with the electronic friction model as described above, we call the GLO model the MDEF-GLO model in what follows.

C. A six-dimensional potential energy surface for CO adsorbed on Ru(0001)
The interaction of CO with a rigid Ru(0001) surface is described by the ground-state PES V (r 1 ,r 2 ) = V (r C ,r O ) of above. A global, six-dimensional ground-state PES has been developed in Ref. [32], which we adopt here. Rather than using Cartesian coordinates of individual C and O atoms, x k , y k , and z k , we occasionally express the PES as a function of the three center-of-mass (COM) coordinates of the CO molecule relative to the surface, X,Y,Z, as well as the interatomic distance r, the polar angle θ ∈ [0,π ] and the azimuthal angle φ ∈ [0,2π ]. The relations between curvilinear and Cartesian coordinates are For the construction of the PES in Ref. [32], a 2 × 2 unit cell with a single CO molecule was chosen corresponding to Ru(0001)(2 × 2):CO and a coverage of 0.25. Using a three-layer slab, more than 90.000 points were calculated with periodic density functional theory, using the RPBE functional [29] and D2 dispersion corrections [33] (RPBE-D2). Based on these points, a symmetry-adapted corrugation reducing procedure [34] (CRP) was used to interpolate the PES. In the current paper, a slightly modified CRP parametrization compared to Ref. [32] was adopted, which accounts for the local C 3v symmetry and allows to distinguish small energetic differences for adsorption in fcc and hcp sites. Namely, using the same computational setup and grid representation of the PES as described in Ref. [32], we calculated additional DFT points, at high-symmetry sites top, fcc, bridge, and t2f [see Fig. 1(a)]. The C 3v irreducible wedge used in the symmetry-adapted Fourier interpolation is represented at the six adsorption sites shown in Fig. 1(a). We further extend the 165447-3 range of validity of the PES to lower Z values, that is, down to Z = 0.9Å.
The PES predicts that CO adsorbs in upright position atop a Ru atom, C down, with an adsorption energy E ads of 1.69 eV. The surface-molecule distance Z eq is 2.59Å and r eq = 1.17Å. The latter value is slightly longer than that of gas phase CO (1.15Å). Energetic and geometric data obtained are in good agreement with experiment, giving Z eq = 2.63 ± 0.05Å, r eq = 1.17 ± 0.05Å, and (for coverages comparable to here) E ads = 1.65 eV [35]. The PES predicts also that dissociation of the molecule is facilitated on the surface compared to the gas phase, however, dissociation barriers remain high, in the order of 4-5 eV (depending on orientation of the molecule). Diffusion, on the other hand, of the entire molecule is possible at much lower energies, with minimum barriers for top-to-top diffusion of around 240 meV. In the experimental reference [36], activation energies between 480 and 270 meV were reported for diffusion, depending on coverage. We shall return to this issue below. We also note that bridge, hcp, and fcc sites are local maxima within our potential with energies of about 260, 270, and 280 meV, respectively, higher than top. More details on the PES and further comparison to experimental data can be found in Ref. [32]. (There, however, an error was made in Fig. 5 where bridge and hcp energies compared to top were obtained from an unconverged steepest-descent calculation.) In Fig. 1, we show a contour plot of the CRP-RPBE-D2 PES in the (X, Y ) space, with other (four) coordinates optimized.

III. RESULTS
In the following MDEF(-GLO) simulations, often a large number of trajectories were propagated over long time periods. Also, a large range of laser parameters (wavelength, fluence, pulse length) or other conditions (initial temperature) were considered, to make contact to different experiments as reported in the literature. The computational settings/parameters will be described in the respective subsections, and are summarized in Appendix A for the reader's orientation.

A. Vibrational lifetimes
The electronic friction model accounts in an approximate way for the coupling of vibrational modes to electronic degrees of the metal surface. By running Langevin trajectories for molecules which are initially excited along a (normal) mode i with vibrational energy ω i initially, this mode will be damped accordingly. The energy loss to the electronic bath is in good approximation exponential E i (t) = ω i e −t/τ i , from where a vibrational lifetime τ i of mode i can be estimated. For each degree of freedom (DOF) we start a Langevin trajectory from the minimum of the PES (perpendicular on top, see above), with all but one (the ith one) of the initial velocities being zero. The initial velocity for mode i is chosen such that the kinetic energy of the adsorbate equals one vibrational quantum of the respective mode, obtained from classical normal mode analysis in Ref. [32]. These modes can be roughly categorized as vibrational motions along X,Y , Z, r, θ , and φ or combinations thereof, called T modes (X,Y , doubly degenerate), R modes (θ,φ, doubly degenerate), S mode (Z), and CO stretch mode (r), respectively. For instance, for the molecule's center-of- mass motion along Z the harmonic frequency is 414 cm −1 , corresponding to 51.3 meV [32], and the CO stretch mode has 2097 cm −1 (260.0 meV). Corresponding experimental frequencies are 445 and 2010 cm −1 , respectively. The initial kinetic energy is then obtained via 1 2 m i v 2 i = ω i , where m i is the vibrational mass of mode i.
In Fig. 2, it is demonstrated that the CO stretch mode (r) and the CO surface mode (Z, the S mode) decay on time scales of 5.5 and 7.5 ps, respectively, within the MDEF model. Lifetimes of other modes R and T (not shown) are similar.
To the best of our knowledge, lifetimes for vibrations of CO on Ru(0001) are not reliably known experimentally, but lifetimes in the order of several ps for Z and r have been observed for similar CO/metal surface systems, both experimentally and theoretically [37]. Typically, however, there are larger differences between lifetimes for r and Z modes than here. Our lifetimes obtained for the internal (r) mode are probably overestimated in the independent-atom LDFA approach [38], while those of the Z mode are probably somewhat underestimated [37]. Frictional approaches beyond the original independent atom LDFA have been suggested [38][39][40][41]. Still, the LDFA accounts for multimode electronic damping in a reasonable and parameter-free way. This is in stark contrast to one-dimensional Arrhenius-Kramers-type rate models which are frequently used to model FL induced desorption [1,7,21,42]. These models also account for electronic friction through an effective, coordinate-independent damping parameterη el = 1/τ , however, in a rather empirical way. To fit FL-desorption data for CO/Ru(0001), for example, in Ref. [7] values of τ in the order of several 100 fs had been used (see below). Assuming Z being the dominant reaction coordinate for desorption, we have to compare this number to the τ Z lifetime, which is about 7.5 ps according to the present model. We will show later that the multimode dynamics employed here accounts also for observed desorption data, despite the large difference in Kramers model and LDFA η values.

B. Laser induced, hot-electron mediated surface diffusion
Next, various femtosecond-laser induced, hot-electron mediated dynamical processes of CO on Ru(0001) are considered. The first process we wish to study is single-pulse induced surface diffusion of CO on the surface. For this purpose, single laser pulses were adopted within the 2TM, with a pulse width (FWHM) of 50 fs (peaking at t = 0), wavelength λ = 400 nm, and fluences in the range between 100 and 200 J/m 2 . The initial temperature was 300 K (cf. Appendix A). The same parameters were chosen in the experimental reference [8] (for an absorbed fluence F = 140 J/m 2 ), where also 2TM modeling was employed in the context of one-dimensional Kramers rate expressions. The T el (t) and T ph (t) curves arising from the 2TM with two different fluences (140 and 200 J/m 2 ), over short-and long-time intervals, are shown in Fig. 3(a). There, the usual behavior is found, i.e., large peak electronic temperatures of several thousand K are observed several tens of fs after the pulse maximum, followed by quasiequilibration with phonons on the time scale of a few ps. The phonon temperature T ph reaches maximal values after a bit more than a ps (of about 1800 K for 140 J/m 2 ), decaying slowly afterwards, on a tens-of-ps time scale.
Adsorbate DOF are excited after a laser pulse has been applied, and the adsorbate is rapidly set in motion by fluctuating forces. While hot-electron excited adsorbate modes are usually not necessarily thermally equilibrated in particular at short times (see, e.g., Ref. [43], but also Ref. [44]), we can nevertheless define an adsorbate temperature from the molecular dynamics with electronic friction (MDEF) as Here, E kin is the total kinetic energy of the adsorbate atoms obtained from MDEF or MDEF-GLO, and applicability of the equipartitioning theorem for each of the six DOFs has been assumed. In Fig. 3(b), we show two different T ads (t) curves, both for the case of a laser pulse with 140 J/m 2 . One is for MDEF only (frozen surface approximation), the second one for the MDEF-GLO model in which the surface is allowed to move. For both curves, the negative times refer to an equilibration period of 10 ps to reach an initial target temperature of 300 K (which in practice is not quite achieved within that time). The pulse form (source term in the 2TM) is Gaussian and peaks at t = 0. Here and everywhere in the paper, we start the MDEF(-GLO) trajectories from the on-top site, which is the global minimum of the PES (see above) and we propagate for further 300 ps after the initial equilibration time. It is seen that the MDEF derived adsorbate temperature rises until about 13 ps, reaching a maximum of ∼1250 K. The two curves (MDEF and MDEF-GLO) are practically identical to each other at longer times (>20 ps), however, at shorter times the MDEF-GLO model predicts a larger peak temperature (up to ∼1500 K), which also occurs earlier, at around 5 ps. The higher The initial temperature is 300 K, the pulse is maximal at t = 0. The inset shows the short-time behavior. (b) Adsorbate temperature T ads over time after excitation with the laser pulse with F = 140 J/m 2 , computed from MDEF and MDEF-GLO kinetic energies and Eq. (9) (see text). The laser pulse in the 2TM has its maximum again at t = 0 fs. The equilibration phase of T ads (t) starts at −10 ps and lasts for 10 ps. (c) T el (t) (solid lines) and T ph (t) (dashed lines) for a Ru surface driven with two laser pulses of 800 nm, with a fluence of 125 J/m 2 and a pulse width of 130 fs each. The first pulse peaks at t = 0, the second one after a delay time of t = 20 ps. The base temperature was 200 K in this case. peak temperature is due to the fact that the surface phonons act now as an additional energy source for the adsorbate, at least 165447-5 at early times. The faster response time of the adsorbate is, on the other hand, a bit of a surprise given the generally accepted notion that phonon degrees couple more slowly than electrons to the adsorbate. In passing, we note that during the action of the pulse some molecules desorb (see below), and we have not counted those for temperature analysis.
The most ubiquitous dynamical process driven by the hot electrons is lateral motion of a CO molecule in the (X,Y ) plane. This can most easily be visualized with a single, "representative" trajectory as shown in Fig. 4, in this case within the MDEF model. (Note that we work with one molecule per 2 × 2 cell which formally corresponds to a coverage of 0.25. However, by using periodic boundary conditions CO molecules in different cells move simultaneously, which corresponds more closely to a lower-coverage situation.) In Fig. 4, the motion of the molecule along the Z coordinate (upper, narrow panels), in the surface (X,Y ) plane (middle right) or along Y (middle left) is shown, as well as the motion along the tilting angle θ (lower left). It is seen that after the thermal preequilibration phase in the time interval [−10,0] ps, the action of the laser pulse sets in, setting the molecule in motion. In the first few ps, however, the gained energy of the molecule is not sufficient to overcome diffusion barriers: the molecule merely oscillates in the adsorption well. Later, the motion along Y (and X) shows up as a terraced structure. Here, vibrations inside a potential well alternate with hops between different adsorption (on-top) wells. Naturally, most hops lead the molecule to a neighboring well, but also long-distance hops over distances of up to 1 nm can be found. Hops occur quasirandom, thus showing typical signatures of (thermal) adsorbate diffusion on the surface.
Before analyzing "surface diffusion" further, we note that also motion along Z is appreciable during and after the action of the laser pulse, and so does the motion along the tilting angle θ . That is, the molecule becomes "hot" also in these nonlateral modes.
Returning to surface diffusion, we compute the meansquared displacement in the surface plane as a quantitative measure for lateral motion, Here, is the starting point of all trajectories. This quantity can be related to a diffusion coefficient D(t), which is time dependent due to the time dependence of the electron temperature. The time-dependent electron temperature can be seen to "translate" via Eq. (9) into an adsorbate temperature T ads (t), which drives the surface diffusion. A diffusion coefficient for two-dimensional diffusion can be defined in the usual way from the slope of R 2 (t) via [This equation gives the usual Fig. 5, the D(t) curves obtained are shown for various laser fluences. To obtain these curves, 20 000 trajectories were run for each fluence, each 300 + 10 ps long (where 10 ps refer to the initial thermalization time to approach the base temperature of 300 K). This setting was used in the rest of this section. Both results for MDEF (a) and for the MDEF-GLO model (b) are shown in the figure. Apart from the noise, it can be seen that D(t) still resembles reasonably well the general shape of the adsorption temperatures T ads (t) of Fig. 3(b). It is further noted that in the first few ps when the adsorbate is still quite "cold," the D values are small while at around t = 10-20 ps, when the T ads (t) values are large, D has its largest values, too. For MDEF, the maximum diffusion coefficients are in the order of 3-9 ×10 −4 cm 2 /s, depending on fluence. For MDEF-GLO these values are further enhanced by about 30%. The FL-derived diffusion coefficients are up to five orders of magnitude larger than purely thermal diffusion coefficients for CO/Ru(0001), measured at 300 K [36]. For example, at a coverage of 0.27 (comparable to ours), D ∼ 3 × 10 −9 cm 2 /s was reported in Ref. [36].
Since at each time t when D(t) is evaluated can be associated with an adsorbate temperature value T ads (t) through curves of the type shown in Fig. 3(b), we can also compute D as a function of the adsorbate temperature. Assuming Arrhenius-type behavior, an Arrhenius plot lnD versus 1/T ads gives a prefactor D 0 , and an "activation" energy E d for diffusion. In  Table I for different fluences and show little scatter around mean values, as it ideally should be. We also note from the table that there is little difference between MDEF and MDEF-GLO. This does not mean, however, that FL induced diffusion is insensitive to surface motion as had been demonstrated in Fig. 5. The diffusion rates in the MDEF-GLO model are larger than for MDEF, due to the higher adsorbate temperatures which can be achieved when phonons participate.   (0001), for laser pulses of 50 fs FWHM, λ = 400 nm, with three different fluences, and a base temperature of 300 K. Averages over all trajectories are shown (20 000 when analyzing desorption, 20 000 minus desorbing trajectories when analyzing diffusion.) E d and D 0 are Arrhenius activation energies and prefactors for surface diffusion, respectively, P des are percentages of desorbed CO molecules. Results for the MDEF and MDEF-GLO models are shown. For F = 200 J/m 2 , we give in brackets also the desorption probability for another condition: λ = 800 nm, FWHM = 130 fs, and a base temperature of 200 K, which fits better to experiments presented in Ref. [6], while the previous conditions fit to Ref. [8]. (Note that in Ref. [ As mentioned above, the lowest RPBE-D2 derived diffusion barrier of CO on Ru(0001) is around 240 meV, in good agreement with the data obtained from the Arrhenius fit. From an Arrhenius analysis of diffusion, in Ref. [36] an activation energy of 270 meV has been given for a coverage of 0.56, and a prefactor of D 0 = 0.06 cm 2 /s. This activation energy fits nicely to our theoretical values. It must be said, however, that according to Ref. [36], E d is strongly coverage dependent, giving E d = 480 meV at a coverage of 0.27. According to Ref. [36], there is also a clear coverage dependence of the prefactor D 0 , changing from D 0 = 0.06 cm 2 /s at coverage 0.58 to D 0 = 0.38 cm 2 /s at coverage 0.27. Our D 0 values derived from FL-diffusion data are smaller, for reasons not known so far. It should be said that both DFT-GGA as well as experimental uncertainties [45] (in particular also for prefactors) can be responsible for observed deviations.
The picture of isotropic diffusion implicitly suggested by our analysis so far is somewhat oversimplified. This is demonstrated in Fig. 7, where, for the MDEF-GLO model, the distribution of (X,Y ) coordinates of all (20 000) trajectories at various times are given, for two fluences: 100 and 140 J/m 2 . There, we note that before the action of the laser pulse, at t = −5 ps (i.e., in the middle of the equilibration phase, frames in upper row), the CO molecules reside in their on-top positions. Shortly after the pulse (at t = 0.5 ps, frames in the second row) the CO molecules start to explore larger regions of the (X,Y ) space on the surface. At intermediate times (t = 5 and 7 ps, frames in the fourth and fifth rows of the figure), the distribution outside the top position becomes nonuniform, and CO molecules tend to reside apart from on-top sites, also in fcc and hcp sites. This preference, however, is only clearly seen for the higher fluence F = 140 J/m 2 , and less clearly for F = 100 J/m 2 . This implies that different laser fluences can lead to different lateral distributions of CO molecules on the Ru(0001) surface. As mentioned above, our PES has no minima at the fcc and hcp and therefore the preferred population of these sites may come as a surprise. We assign this behavior to a dynamical trapping effect, which is also supported by the observation that at long times (at t = 30 ps, frames in the lowest part of the figure) the local population maxima at fcc and hcp sites are diminished again, even for F = 140 J/m 2 .
In the context of Fig. 7, we note that in the experimental work [9], also two different outcomes of FL-pulse excitation (50 fs FWHM, 400 nm) of CO/Ru(0001) were observed for F = 100 and 140 J/m 2 . In that reference, for the low fluence the dominant creation of "hot" CO molecules (with average distance closer to the surface than initially) was postulated, while at the higher fluence the population of weakly bound precursor states farther away from the surface was assumed to take place. We will return to the question of physisorbed precursor states later. We mention here already, however, that different lateral displacements for different fluences in principle offer a possible alternative to the interpretation given in Ref. [9].

C. Femtosecond-laser induced formation of "hot" adsorbates
As already indicated in the single-trajectory case of Fig. 4 (for 140 J/m 2 ), FL stimulation of CO/Ru(0001) not only induces lateral motion on the surface and diffusion associated with it, but also motion along other modes such as Z and θ are excited. This is further demonstrated in Fig. 8, where the distribution of coordinates Z and θ as a function of time, together with the average value for all trajectories, are given for F = 200 J/m 2 . Results for both the MDEF-GLO and MDEF-only models are shown. From inspection of the Z coordinate, we note that the molecules move at early times (after laser excitation) on average slightly towards the surface because they tend to drift away from the top sites to spots where a shorter molecule-surface distance is energetically more favorable according to our RPBE-D2 PES [32]. Note, however, that the distribution is broad and some molecules desorb, particularly at early times. Inspection of the θ coordinate proves that θ > 0 is acquired in the course of FL induced dynamics, also with a rather broad distribution. Also, this observation is in line with findings of Ref. [32], reporting a tilt of up to 23 • away from perpendicular once the molecule is displaced from on top. This way, the strength of the C-Ru bond is maximized. Acquiring large θ values with up to θ ∼ 60 • at early times, however, is energetically costly, which may also explain why molecules desorb relatively early.
In summary, both the Z and θ coordinates indicate that the CO molecules are "hot" after FL excitation. We also see no large differences in this respect when comparing MDEF and MDEF-GLO models.

D. Femtosecond-laser desorption of CO from Ru(0001): Single pulses
As seen from Fig. 8, within the first picoseconds after laser pulse excitation, desorption is a possible reaction channel. This is further demonstrated in Table I, where also desorption probabilities are listed. Trajectories were counted as desorbed if Z exceeded 5.1Å at the end of the trajectory. Most desorption probabilities in the table are given for the three pulse fluences considered so far, with laser parameters λ = 400 nm, FWHM = 50 fs, and a base temperature of 300 K. For a selected fluence, F = 200 J/m 2 , also alternative conditions: λ = 800 nm, FWHM = 130 fs at a base temperature of 200 K, were considered (cf. Appendix A). These fit better to an older experimental reference [6], where the FL desorption of CO from Ru(0001) has been studied in some detail. In contrast, the previous set of parameters fits better to the newer experiments done in Ref. [8]. From Table I, we already note that (i) the desorption probability is clearly non-negligible (in the order of a few percent) for "typical" conditions matched by current FL experiments, and (ii) choosing different laser parameters affects the desorption probability. Specifically, there are moderately strong dependencies on laser wavelength and pulse width (with longer wavelengths and/or longer pulses at a given fluence leading to less desorption). The main dependence arises from a strong increase of desorption yield with laser fluence. This latter observation is a well-known effect for FL induced desorption [1], often empirically modeled as a power-law fit where n is a characteristic exponent. The latter can be related roughly, but not directly, to the number of excitationdeexcitation cycles in a DIMET (desorption induced by multiple electronic transitions) scenario, as theoretically demonstrated elsewhere [43]. As a third feature from Table I, we see that (iii) lattice motion (considered via the GLO model) favors desorption considerably. In fact, P des increases easily by an order of magnitude when phonons are accounted for, in this case from about 3.6% to about 25%. (iv) Finally, one finds from a plot of the averaged desorption rate dP des /dt versus time, that in the MDEF-GLO model the desorption rate is maximal at earlier times compared to MDEF, namely at around 15 ps for MDEF-GLO compared to slightly above 20 ps for MDEF (not shown). This may come as a surprise since phonon response is usually considered to be slower than pure electronic response, however, this picture is consistent with the T ads (t) curves presented in Fig. 3(b), and also with the two-pulse correlation traces to be presented in Sec. III E. Next, by using the "new" parameter set (λ = 800 nm, FWHM = 130 fs, and a base temperature of 200 K), we computed the desorption probability of CO from Ru(0001) by systematically increasing the laser fluence. In our simulations, nine different fluences in the range [175, 375] J/m 2 were chosen. For each fluence, 5000 trajectories were calculated, all being 300 + 10 ps long as before (including 10-ps equilibration phase). In Fig. 9, we plot desorption probabilities as a function of laser fluence, for the MDEF and MDEF-GLO models comparing to experimental data from Ref. [6], where a similar range of fluences had been considered. A . Results from MDEF-GLO and MDEF models are shown, and experimental data from Ref. [6]. The laser wavelength was 800 nm, FWHM =130 fs, and the base temperature was 200 K. The experimental data were represented such that the desorption yield for F = 305 J/m 2 corresponds to 0.2, the value stated in Ref. [6] for this particular fluence.
double-logarithmic representation is chosen to obtain from a linear fit the exponent n according to Eq. (13). From the figure, the following observations can be made.
(i) Theoretical desorption yields follow approximately, but not exactly, a power law of the form of Eq. (13), with exponents n ∼ 5.2 (MDEF-GLO) and n ∼ 7.0 (MDEF), respectively. In particular for MDEF-GLO the desorption yields increase slower than exponential at least at large fluences, where P des approaches 1 and saturates.
(ii) MDEF-GLO desorption yields are larger than MDEF yields, notably at smaller fluences where the MDEF-GLO values are typically an order of magnitude higher than MDEF. This indicates once more a clear participation of phonons, consistent with the findings in Table I. (iii) The comparison with experimental values from Ref. [6] is fair: experimental yields are somewhat lower than MDEF-GLO (the currently most accurate theoretical model), and an exponent n ∼ 4.8 is found experimentally [6], not far from our MDEF-GLO value of 5.2. Note that, in experiment, slightly different conditions had been used compared to here, e.g., a coverage of 0.68 [rather than the formal 0.25 or lower value (see above) as in our model]. Note also that the damage threshold of a Ru crystal is near 270 J/m 2 (for 800 nm, 130 fs pulses) [46].
In this context, we mention that similar experiments as those in Ref. [6] have recently been done in Ref. [7], by Gladh et al.
Like in the present work (cf . Table I), the shorter wavelength leads at a given constant fluence to larger yields according to Ref. [7]. This enhancement effect is due to the smaller penetration depth ξ of light with shorter wavelengths [7], which concentrates energy near the surface. Quantitatively, Ref. [7] suggests a factor of 3-4 yield enhancement. Our enhancement factor is only slightly larger: at 200 J/m 2 , it is about 5 according to Table I. However, the comparison is not entirely straightforward, as for the shorter wavelength we also employed a shorter laser pulse and higher initial temperature (cf. Table III in Appendix A), which both enhances desorption as well.
In further agreement with the present theory, Gladh et al. emphasize the role of both electronic and phononic contributions to desorption. As said earlier, however, their interpretation is based on (electronic) friction coefficients in the order of (100 fs) −1 . The latter should be disputably considered according to the analysis in Sec. III A.

E. Femtosecond-laser desorption of CO from Ru(0001):
Two-pulse correlation By using two pulses rather than one laser pulse, desorption yields can also be enhanced. In addition, two-pulse correlation (2PC) experiments in which a second pulse is delayed from a first one by a delay time t give insight into mechanistic details of desorption, e.g., electronic versus phonon mechanisms [1]. Dominantly electronic mechanisms are generally accepted to be characterized by short correlation times (∼1 to a few ps), while phonon participation results in correlation times ∼10 to several tens of ps. A "correlation time" can loosely be defined as the HWHM (half-width at half-maximum) of a yield versus delay time plot of the desorption probability, vide infra.
Here, we have calculated 2PC spectra for desorption, using the following settings (cf. Appendix A): two identical pulses, each with λ = 800 nm and with a fluence of 125 J/m 2 , pulse width 130 fs; the base temperature was 200 K, trajectories were 300 + 10 ps long (with 10-ps equilibration phase), and 5000 trajectories were run per pulse delay. Eleven pulse delays between 0 and 150 ps were considered. In Fig. 3(c) we show the time evolution of T el (t) and T ph (t) obtained from the 2TM, when the pulse delay was t = 20 ps. Due to the usage of two identical pulses, P des ( t) values in 2PC spectra are gerade functions of t, i.e., P des ( t) = P des (− t). Here, negative delays refer to the second pulse preceding the first.
Again, there are two sets of experiments we can compare with, one from Ref. [6] and one from Ref. [7]. We shall refer to data from Ref. [6], but very similar conclusions could be drawn from the data of Ref. [7]. In Ref. [6], experimental settings were as follows: two pulses with λ = 800 nm, total fluence 250 J/m 2 , however, with slightly different individual pulse fluences (ratio F 1 : F 2 = 52 : 48), pulse lengths FWHM =130 fs; base temperature 200 K. In that case, 2PC spectra are slightly unsymmetric.
From Fig. 10(a), where P des ( t) curves are shown for MDEF and MDEF-GLO models, we note first of all that desorption yields are indeed grossly enhanced when using two pulses rather than one. The maximum desorption yield occurs at around t = 3 ps [and not t = 0, for reasons which have been explained elsewhere (see, e.g., Ref. [In Ref. [6], no absolute experimental values are given for 2PC yields. We estimate them by deriving from the reported, experimental single-pulse yield for F = 305 J/m 2 of 0.2 (see above), the absolute single-pulse yield at F = 250 J/m 2 (0.071), and equaling this value to the yield for two strongly overlapping pulses ( t ∼ 0) having the same total fluence.] Each theory data point was obtained with 5000 trajectories. Only the positive delay times were calculated, data with negative t were simply obtained by mirroring the data for positive delays. Lines between data points are drawn to guide the eye. Dashed vertical lines in (a) indicate correlation times (of about 18 and 29 ps, respectively), after which desorption yields dropped to 1 2 of their maximum value (HWHM values).
for longest delay times. (A pure single-pulse scenario with this fluence gives a desorption probability of less than 0.02%.) A second feature emerging from the figure is the rapid decay of the 2PC trace with increasing t. The correlation time, estimated from the HWHM of the two curves (the black and red vertical dashed lines in the figure), is around 29 ps (MDEF) and 18 ps (MDEF-GLO), respectively. Our theoretical model gives larger correlation times than the experiment, although they are still on the same order of magnitude. This is demonstrated in Fig. 10(b), where the MDEF-GLO model is compared with experimental data from Ref. [6], plotted on a half-logarithmic scale. It is seen that the falloff of the experimental curve gives a HWHM of around 10 ps, roughly half the value from our theory. The same order of magnitude for the correlation time has also been found in the newer experiments of Ref. [7], both at 400 and 800 nm. Correlation times in the range of a few tens of ps indicate, according to current interpretation, both electronic and phononic contributions to desorption [7]. Note, however, that in our theory the correlation time for the MDEFonly model is even longer than that of the MDEF-GLO model according to Fig. 10(a). This suggests that the connection of correlation times in 2PC spectra with mechanistic (electron versus phonon) interpretation is not entirely straightforward.

F. Physisorbed precursor states?
In recent experiments it has been suggested that FL stimulation of CO/Ru(0001) could not only lead to desorption and (lateral) diffusion, but to another interesting process as well [8,9]: After applying a 140-J/m 2 , 400-nm, and 50-fs laser pulse, on a time scale of >2 ps, the population of a physisorbed molecular state farther away from (rather than closer to) the surface at Z ∼ 4Å was suggested to take place.
[At shorter times (up to ∼1 ps), the formation of a higher-metal coordinated adsorbate was postulated [9]]. This conclusion was based on (i) real-time observation by x-ray absorption (XAS) and emission (XES) spectroscopy on the one hand, and (ii) on computed potential of mean force (PMF) curves along the reaction path (here Z, the molecule-surface distance) on the other. Referring to the PMF, the latter is defined as [8] at different temperatures T . Here, V 0 (Z) is the minimal value of the multidimensional PES at a certain Z value and other coordinates optimized, and where A is an arbitrary factor (which cancels), and V (Z,q) = V (Z,q) − V 0 (Z) is a potential energy difference. Assuming for now a rigid surface, q are the other five degrees of freedom of CO relative to the substrate, and V (Z,q) the 6D PES of CO interacting with a Ru(0001) surface. In practice we choose coordinates u = X − Y/ √ 3, v = 2Y/ √ 3, r, θ , and φ, and the volume element is dq = sin(θ ) dX dY dr dθ dφ. In Ref. [8], one of the coordinates (r) was frozen, and several model assumptions were made, including separability of translational and rotational degrees of freedom. Using the so-called BEEF-vdW DFT functional [47], in Ref. [8] PMF curves were then found to exhibit the usual chemisorption well at short molecule surface distances (Z ∼ 2.5Å), and at higher temperatures (around 300 K and above) a physisorption well at around Z ∼ 4Å in addition. The latter was around 80 meV deep, and separated by a clear barrier (of slightly above 0.3 eV at T = 2000 K) from the chemisorption well. According to Ref. [8], about 30% of the CO molecules are trapped in the physisorption precursor state after 140 J/m 2 laser-pulse excitation.
In our present MDEF simulations, which account for nonthermal multidimensional dynamics, no such trapping of molecules in a physisorption well could be found. In none of the trajectories such an event occurred. This is also evident from Fig. 8 (upper panels), where no persisting population around Z ∼ 4Å is found. Rather, particles strongly excited along Z either desorb (going to Z > 5.1Å) or return quickly to the chemisorption well. Also, using our fully coupled, six-dimensional RPBE-D2 based CRP PES we recalculated PMF curves along Z. That is, all (other) five degrees of freedom were considered, and no separability of degrees of freedom was assumed. The five-dimensional integration in Eq. (15) was done numerically for various (50) Z points in a range [0.9,6.0]Å and for temperatures in the range of [0, 4000] K. Equation (15) was integrated numerically, using around 4 × 10 7 points of the V (Z,q) fit of our RPBE-D2 PES, 1 for each Z. For comparison, we also did the same calculation with the pure RPBE potential, i.e., without vdW corrections. The results are shown in Fig. 11(a) (for RPBE-D2) and 11(b) (for RPBE), respectively.
Considering the van der Waals corrected RPBE-D2 PMF curves [ Fig. 11(a)], we note only extremely shallow physisorption wells (at most ∼20 meV deep) and almost no barrier between physisorption and chemisorption wells. Therefore, no non-negligible barrier between physisorption and chemisorption exists according to our PMF, and thus no clear precursor state can be identified in which CO molecules could be trapped. When "switching off" the van der Waals correction [RPBE PMF curves in Fig. 11(b)], a barrier towards desorption emerges at higher temperatures, as in Ref. [8], however then, as expected, there is no physisorption well any more at all.
The discrepancy between results of Ref. [8] and ours can have different reasons. Note that both in Refs. [8,32], three-layer slabs of a Ru(0001)(2 × 2):CO (formal coverage 0.25, see above) were adopted, and in both cases no phonon contributions to the PMF have been considered. There are some differences between RPBE-D2 as used here and BEEF-vdW as adopted in Ref. [8], however. Namely, our PES gives 1 Grid points for PMF calculations at various Z values were chosen equidistantly for four out of five coordinates: N r = 50 for r ∈ [0.9,2] A, N u = N v = 25 for u,v ∈ [0,L], and N φ = 35 for φ ∈ [0,2π ]. L = 2.725Å is the lattice constant in a 1 × 1 cell. The polar angle θ ∈ [0,π ] was represented on N θ = 37 Gauss-Legendre points. The PMF calculations were done with the (C 6v ) potential fit of Ref. [32], not with the slightly improved (C 3v ) fit of this work; we expect only small differences. approximately the mentioned, correct experimental adsorption energy of about 1.65 eV [35] {the RPBE-D2 potential well [and W (Z) at T = 0 K] is 1.69 eV deep, the zero-point energy of the CO-surface bond is around 0.025 eV, see above}, while DFT calculations of Ref. [8] seem to underbind chemisorption somewhat: the PMF (T = 0) potential well is only about 1.4 eV deep according to Fig. 3 of Ref. [8]. Further, we cannot exclude that our RPBE-D2 may somewhat underbind the van der Waals interaction in the physisorption region. On the other hand, the PMF (T = 0) in the region relevant for the possible existence of the barrier of Ref. [8] is found to be very similar to our RPBE-D2 derived PMF (T = 0) curve. So, one may attribute differences to the entropic contributions, arising 165447-12 from different assumptions made for the PES. For the latter, we adopt here a fully coupled, anharmonic 6D PES, while in Ref. [8] translational and rotational degrees of freedom were assumed to be separable. In Appendix B, we shall give evidence that neglecting mode coupling can be problematic, and that artificially neglecting mode coupling leads to clear barriers also with our RPBE-D2 potential.
Beyond the PMF picture, the more striking argument against the occurrence of precursor states, for the low-coverage adsorbate model at least, is anyway the multidimensional dynamics. The latter is free of the effective one-dimensional, free-energy path approximation and also accounts for phonon effects through the GLO model. Having said this, we are aware of the fact that also our PES can only be approximate, and that the inclusion of phonons via the GLO is quite rudimentary. Further, real-world conditions arising from sample preparation (with defects, coverage variations, clustering of CO molecules) are not captured by the model, and might lead to the occurrence of precursor states which could be populated in experiment. Another possibility for the interpretation of the experimental data of Ref. [8] is that alternatives/additional contributions may exist which have not been considered so far. It remains to be seen if such alternative scenarios, e.g., the preferential population of different lateral areas on the surface (as a function of time and/or fluence) (cf. Fig. 7) could also explain the characteristic changes in XAS and XES spectra in experiment [8].

IV. SUMMARY AND CONCLUSIONS
In summary, a Langevin model accounting for all six molecular degrees of freedom of CO on a Ru(0001) surface (at formal coverage 0.25) was applied to femtosecond-laser (FL) induced dynamics of the adsorbate on this substrate. The model is based on first-principles (RPBE-D2) potentials, accounts for electronic friction and hot-electron mediated excitation within the local density friction approximation (LDFA), and via the generalized Langevin oscillator (GLO) model also for substrate motion. No adjustable parameters enter the model, only system-specific parameters for the GLO and two-temperature models (2TM) were adopted. Since the PES and electron densities were precomputed, a large number of trajectories could be run over long time periods and analyzed with small statistical errors. As such, our treatment allowed us to go beyond previous theoretical work of FL induced dynamics of CO/Ru, which was based either on empirical Arrhenius-Kramers-type kinetics, thermodynamic analysis [e.g., via potential of mean force (PMF) arguments], or ab initio molecular dynamics (AIMD) over short times and at constant temperature. The model allowed us to gain (within its limits, which clearly exist) detailed insight into FL induced dynamics of CO/Ru(0001). Returning to our original physical questions posed in the Introduction, the following observations were made: (i) FL stimulation of CO/Ru(0001) leads to various dynamical processes, such as (lateral) diffusion on the surface, hot adsorbate formation, and molecular desorption. Within the present model, no indication of the population of physisorbed precursor states was found. On the other hand, the lateral distribution of CO on the surface can lead to favorable patterns (e.g., population of fcc and hcp sites) as a function of time and laser fluence. Apart from this "nonisotropic" diffusion behavior, diffusion coefficients can be obtained from an Arrhenius-type analysis of the MDEF(-GLO) data which are grossly enhanced compared to those found in thermal diffusion. Also, desorption yields (and rates) obtained after single-or double-pulse stimulation are much larger compared to thermal processes.
(ii) As far as microscopic mechanisms are concerned, we note that diffusion is somewhat affected by phonons, in addition to hot-electron mediated dynamics. The effect is relatively moderate, however, leading to only slightly enhanced diffusion coefficients within the MDEF-GLO model, compared to MDEF. In contrast, phonons affect the hot-electron mediated desorption dynamics considerably, leading to desorption yields easily an order of magnitude larger in the MDEF-GLO model, compared to MDEF. Perhaps surprisingly, the phonons act on similar time scales as the electron-hole pairs, at least in the sense that they enhance the energy transfer to the molecule from the surface on a few-ps time scale. It must be said, however, that our Langevin treatments of electronic friction and phonon motion both do not account for memory effects and may therefore not fully capture time-scale separations of electrons and nuclei. Here, Langevin methods which take memory into account [48] and/or account more explicitly for substrate phonon motion [49] should be consulted in the future to make a more definitive statement.
(iii) For desorption, we find the usual nonlinear increase of the yield with laser fluence, ideally modeled as a power-law dependence. For diffusion, increasing the fluence increases the (diffusion) rates, too, presumably in a more complicated fashion. The nonlinear (power-law) increase of desorption yields is the basis for the distinct enhancement of desorption in 2PC scenarios, for short-time delays. Further, while we did not analyze this in great detail, there is a dependence of the dynamics on laser wavelength and pulse length. For instance, shorter wavelengths (400 versus 800 nm) lead to larger desorption yields due to the smaller penetration depth of the laser into the solid. It will be a valuable task for the future, to exploit the possibilities of incoherent laser control of hot-electron mediated reactions at metal surfaces, using MDEF(-GLO) models.
(iv) Our theoretical findings are generally in qualitative, partly in quantitative, agreement with experimental observations. In addition, our approach allows for detailed timeand space-resolved insight into the ongoing dynamics. One discrepancy to an existing experiment or, more precisely, its current theoretical interpretation, is the lack of FL induced trapping of CO molecules in physisorbed precursor states. Of course, this statement is subject to the possible limitations of the current model, which still exist. It will be interesting to see if alternative explanations of the mentioned XAS and XES results obtained for CO/Ru are possible. Also, a more precise modeling of experimental conditions (notably the coverage conditions used there) might alter the picture somewhat, also within the MDEF model.
Work along the mentioned lines is currently being performed in our laboratories. In particular, the modeling of XAS and XES spectra for this system, using a frictional trajectory-based approach, is underway.  Table II. Two-temperature model. Parameters and derived quantities entering the 2TM equations (5) and (6) are as follows [6,7]: The heat capacity of the phonons C ph is calculated from the Debye model of the solid (see equations in Ref. [6]); necessary parameters there are the atom density ρ of the solid and its Debye temperature T D , for which we took ρ = 7.4 × 10 28 atoms/m 3 and T D = 600 K, respectively. Further, the source term in Eq. (5) is calculated from where R is the substrate reflectivity, I (t) the laser intensity profile (assumed Gaussian), and ξ the penetration depth, which is ξ = 6.9 nm for a laser wavelength λ of 400 nm, and ξ = 16.2 nm for λ = 800 nm [6,7]. The GLO model for substrate motion. In the GLO model to treat substrate motion of Ru(0001) in an approximate way, the friction coefficient for phonons η ph is derived from the Debye frequency ω D of the solid as η ph = m s πω D /6, ω D being taken as ω D = 1.758 × 10 −3 E h . Oscillator frequencies entering [25] are chosen as ω x = ω y = 8.45 × 10 −4 E h and ω z = 7.35 × 10 −4 E h for the Ru(0001) surface [50]. Further, m s is the mass of a single Ru atom. Coupling constants in gs are obtained from the phonon frequencies as described in Ref. [23].
MDEF protocol and parameters for propagation and laser pulses. In the calculations of Secs. III B-III E, in all  Table III.

APPENDIX B: POTENTIAL OF MEAN FORCE (PMF) ANALYSIS
As noted in the main text, for the calculation of the potential of mean force (PMF) the authors of Ref. [8] use the so-called BEEF-vdW DFT functional, while our PES is based on the RPBE-D2 DFT functional. The use of two different functionals results in the following differences in PMF at temperature T = 0 K (compare Fig. 11 of the paper and Fig. 3 of Ref. [8]): (i) In our case, the adsorption energy (∼1.7 eV) is larger than the BEEF-vdW adsorption energy (∼1.4 eV) and thus closer to experiment. (ii) In the range Z = 3.0-3.5Å, both van der Waals corrected DFT functionals give quantitatively very similar energy values (which we checked by comparing PMF curves at 0 K of Ref. [8] to ours). This Z range is the one in which the transition state energy barrier at higher T was found in Ref. [8], while it is negligible in our calculations.
[The barrier from physisorbed to chemisorbed state is <5 meV according to our RPBE-D2 PMF (T = 2000 K), but more than 300 meV according to the BEEF-vdW PMF (T = 2000 K).] The existence of a non-negligible barrier is crucial in order to have a well-defined precursor state in which molecules could be trapped. (iii) In the range Z = 3.5-5.0Å, the BEEF-vdW functional gives smaller energies compared to the RPBE-D2 DFT functional. The difference is around 60 meV. (Our RPBE-D2 PMF gives a physisorption minimum ∼20 meV deep almost independent of temperature, while the BEEF-vdW functional suggests a physisorption minimum ∼80 meV deep, also almost independent of T .) Based on facts (ii) and (iii), one may speculate that the transition state energy barrier that exists in Ref. [8], but not in our calculations, is due to the different approaches in calculating the PMF and not due to the chosen DFT functional itself. The main difference between our approach and that of Ref. [8] is that in evaluating PMFs we use a fully coupled 6D PES, while in Ref. [8] it was assumed that the potential can be separated in translational and rotational contributions. The aim of this appendix is to show that this separability assumption is probably not well justified, and that this is the reason for the appearance of a non-negligible transition state energy barrier at higher T in Ref. [8].
For this purpose, we will show that such a barrier also appears with the RPBE-D2 DFT functional if one uses the procedure of Ref. [8], i.e., if one makes the separability approximation. The information about the PMF calculation of Ref. [8] can be found in their accompanying Supplemental Material and some more details are also available in Ref. [13]. Unfortunately, complete detailed information how the PMF was obtained is not presented, so in some cases we have to resort to some reasonable assumptions. Fortunately, these TABLE III. Parameters for propagation and laser pulses used in the various sections of this work. For T (0), λ, F , and t see text; in addition, N traj is the number of trajectories per fluence F or delay time t, respectively, and t max the maximum propagation (after the initial equilibration phase of 10 ps).
where V 0 is the minimum of the PES for a given Z and where V (Z,q) = V (Z,q) − V 0 (Z) and dq is the volume element of the configurational space of other degrees of freedom (q), different from Z. [The prefactor A in the corresponding equation of the main text has been left out here, as it cancels in Eq. (B1).] In Ref. [8], separability of V (Z,q) is assumed such that V trans is the translational potential energy, V rot1 is the cartwheel rotational potential energy, and V rot2 is the helicopter rotational potential energy. The latter is considered only in case Z > 4Å. The vibrational component is ignored. Inserting Eq. (B3) in (B2), one obtains Note that in Refs. [8,13], g(Z) has been defined as instead, which is incorrect in view of Eq. (B3). Nonetheless, using either of the equations leads to qualitatively similar PMF curves as will be shown below. The translational potential energy in Ref. [8] was assumed to be representable by the following function: where d is the distance between two top sites. Further, the cartwheel rotational potential energy was represented as The helicopter rotational potential energy was represented as where C is the difference between the highest and the lowest potential energy for that degree of freedom. As discussed in the original work Ref. [8], the helicopter mode is not a normal mode at Z < 4.0Å. Due to this at Z < 4.0Å both rotations are represented by the cartwheel mode V rot1 so that g(Z) = g 2 trans (Z)g 2 rot1 (Z).
Details on the energy points used for the fitting of coefficients a i , b i , and C were not given in Refs. [8,13]. However, evaluation of our PES is very fast and we have checked that using a large number of points or just a few does not change the results noticeably. Also, it is not known whether in Ref. [8] some of the degrees of freedom were relaxed when the molecule was displaced from the minimum position for given Z. Therefore, we have used the following procedure. First, to obtain V 0 (Z) for each Z the molecule was positioned on the top site of the surface with the C atom closer to the surface and the molecular axis perpendicular to the surface and then relaxed. Afterwards, the molecule is displaced from the V 0 position along the particular degree of freedom (X,θ,φ) while keeping all the rest of degrees of freedom fixed. This gives us the energy points to which we fit the coefficients a i (Z), b i (Z), and C(Z).
Originally, in Ref. [8] the coefficients a 0 , a 1 , and a 2 were fitted to a set of energies for two molecules in a 2 × 4 unit cell by laterally moving one of the molecules along (probably) the X direction. In our PES, we do not have a second molecule so we will set a 2 to a reasonable value by hand. We will show that this choice does not affect the conclusions either.
Once we obtain a i , b i , and C for several Z values, we can evaluate the integrals in Eq. (B4). In practice, we will calculate W such that g(Z = ∞) = g(Z = 4.75Å). In the original paper, the authors use g(Z = ∞) = g(Z = 6.75Å).
The obtained PMF is shown in Fig. 12(a). Clearly, this PMF looks similar to the one obtained in the original paper in the barrier region. In particular, at T = 2000 K, a barrier of 0.6 eV arises. Quantitatively, this barrier is almost 0.3 eV larger than the one from Ref. [8] (being slightly higher than 0.3 eV, see above). However, if we use Eq. (B5) for g(Z), in the same way as in Ref. [8], one obtains the results shown in Fig. 12(b). In this case, the barrier for T = 2000 K is around 0.3 eV, very close to the value obtained in Ref. [8].
In both cases, the physisorption well is not visible in our PMF. This is due to the fact that BEEF-vdW is more attractive in the region Z = 3.5-5.0Å, as discussed above. Additionally, we set V 0 (Z = 4.75Å) = 0 and thus we lose some of the attractive long-range van der Waals contribution to the energy (see Fig. 11 in the main text).
Let us now discuss the importance of the parameter a 2 which we have added "by hand" and which could in principle also be a reason why we disagree with Ref. [8] in the barrier issue. In Fig. 13(a), we show the PMF in case a 2 = 0, so that the interaction with a second molecule (particle in a box) is neglected. The barrier at T = 2000 K decreases by 0.2 eV, compared with Fig. 12(a). Therefore, in case the value we choose by hand is overestimated, one would obtain somewhat smaller energy values.
In Fig. 13(b), we show the PMF in case a 2 = 100 times the a 2 value of Fig. 12 so that the interaction with the second molecule is overestimated. The barrier only slightly increases. All in all, it can be stated that the differences induced by changing the a 2 coefficient are not large enough to invalidate our main conclusions.
In conclusion, employing the same (or a very similar) procedure as given in Ref. [8], but using potential values obtained with the RPBE-D2 functional, we also obtain a clear barrier of several tenths of an eV in the PMF for higher temperatures. Because such a barrier is not (or almost not) present when the full, coupled potential is used instead of a separable one, we can conclude that the barrier is an artifact of the separability approximation, at least when we accept the picture that the RPBE-D2 method accounts reasonably well for mode couplings.
In fact, assuming such a separable potential seems not to be a well-controlled approximation according to details of the coupled, 6D RPBE-D2 potential. For the latter, we observe that in case of CO/Ru(0001), lateral X,Y motion and rotation along the polar angle θ are strongly coupled as the molecular axis tends to point to the closest Ru atom in order to preserve the C-Ru bond. In Ref. [8], precisely these degrees of freedom were uncoupled in the final expression for the potential.
In summary, we have shown that a separation between physisorption and chemisorption well by a non-negligible barrier (several tenths of an eV) in PMF curves at higher temperatures is a consequence of neglecting mode coupling, which most probably does not hold in practice. As a consequence, no distinct precursor state exists in the coupled situation (and also stronger van der Waals forces will not change this picture), and trapping in such a precursor state appears to be unlikely: the molecule will either desorb or fall back into the chemisorption well, depending on its energy, in full agreement with our multidimensional dynamics.