Vibrational response and motion of carbon monoxide on Cu(100) driven by femtosecond laser pulses: Molecular dynamics with electronic friction

Carbon monoxide on copper surfaces continues to be a fascinating, rich microlab for many questions evolving in surface science. Recently, hot-electron mediated, femtosecond-laser pulse induced dynamics of CO molecules on Cu(100) were the focus of experiments [Inoue et al. , Phys. Rev. Lett. 117 186101 (2016)] and theory [Novko et al. , Phys. Rev. Lett. 122 , 016806 (2019)], unraveling details of the vibrational nonequilibrium dynamics on ultrashort (subpicoseconds) timescales. In the present work, full-dimensional time-resolved hot-electron driven dynamics are studied by molecular dynamics with electronic friction (MDEF). Dissipation is included by a friction term in a Langevin equation which describes the coupling of molecular degrees of freedom to electron-hole pairs in the copper surface, calculated from gradient-corrected density functional theory (DFT) via a local density friction approximation (LDFA). Relaxation due to surface phonons is included by a generalized Langevin oscillator model. The hot-electron induced excitation is described via a time-dependent electronic temperature, the latter derived from an improved two-temperature model. Our parameter-free simulations on a precomputed potential energy surface allow for excellent statistics, and the observed trends are conﬁrmed by on-the-ﬂy ab initio molecular dynamics with electronic friction (AIMDEF) calculations. By computing time-resolved frequency maps for selected molecular vibrations, instantaneous frequencies, probability distributions, and correlation functions, we gain microscopic insight into hot-electron driven dynamics and we can relate the time evolution of vibrational internal CO stretch-mode frequencies to measured data, notably an observed redshift. Quantitatively, the latter is found to be larger in MDEF than in experiment and possible reasons are discussed for this observation. In our model, in addition we observe the excitation and time evolution of large-amplitude low-frequency modes, lateral CO surface diffusion, and molecular desorption. Effects of surface atom motion and of the laser ﬂuence are also discussed.

, unraveling details of the vibrational nonequilibrium dynamics on ultrashort (subpicoseconds) timescales. In the present work, full-dimensional time-resolved hot-electron driven dynamics are studied by molecular dynamics with electronic friction (MDEF). Dissipation is included by a friction term in a Langevin equation which describes the coupling of molecular degrees of freedom to electron-hole pairs in the copper surface, calculated from gradient-corrected density functional theory (DFT) via a local density friction approximation (LDFA). Relaxation due to surface phonons is included by a generalized Langevin oscillator model. The hot-electron induced excitation is described via a time-dependent electronic temperature, the latter derived from an improved two-temperature model. Our parameter-free simulations on a precomputed potential energy surface allow for excellent statistics, and the observed trends are confirmed by on-the-fly ab initio molecular dynamics with electronic friction (AIMDEF) calculations. By computing timeresolved frequency maps for selected molecular vibrations, instantaneous frequencies, probability distributions, and correlation functions, we gain microscopic insight into hot-electron driven dynamics and we can relate the time evolution of vibrational internal CO stretch-mode frequencies to measured data, notably an observed redshift. Quantitatively, the latter is found to be larger in MDEF than in experiment and possible reasons are discussed for this observation. In our model, in addition we observe the excitation and time evolution of large-amplitude low-frequency modes, lateral CO surface diffusion, and molecular desorption. Effects of surface atom motion and of the laser fluence are also discussed. DOI: 10.1103/PhysRevB.100.245431

I. INTRODUCTION
Femtosecond-laser (FL) induced dynamics of molecules on metal surfaces are interesting for materials processing, time-resolved spectroscopy, surface photochemistry, and photocatalysis. More importantly, sometimes completely new phenomena emerge beyond what is usually found under "ordinary," i.e., thermal equilibrium conditions. For example, "hot electrons" (or electron-hole pairs) are created in the metal, which efficiently couple to adsorbate vibrations. This can lead, compared to thermal reactions, to enhanced reaction cross sections, new reaction pathways, nonthermal product energy distributions, and sometimes entirely different products [1][2][3][4][5]. Here in particular, CO-covered transition metal surfaces were often in the focus, because of their relevance * Corresponding author: peter.saalfrank@uni-potsdam.de in heterogeneous catalysis, and as microlabs for fundamental surface science.
By combining femtosecond laser pulse excitation with subsequent, time-delayed probe pulses detailed insight can be gained into the microscopic, nonequilibrium dynamics of adsorbates on ultrashort (subpicoseconds) timescales. For instance, the desorption dynamics of CO on Ru(0001) was recently studied in real time, by probing the FL-driven, "hot" adsorbates with x-ray pulses [6,7]. Theoretical models were developed in Ref. [6] and elsewhere [8] to interpret the experimental data.
Here we concentrate on another popular system of nonadiabatic surface dynamics, CO adsorbed on Cu(100). In particular, the main motivation for the present work are recent experiments and simulations done by Inoue et al. [9], who studied FL-induced induced dynamics of CO on Cu(100) at a coverage of about 1/2 (one molecule CO per √ 2 × √ 2 cell), using vibrational sum frequency (VSF) generation in included via damping terms and stochastic forces, respectively. The method is also known as the molecular dynamics with electronic friction (MDEF) approach, and has been applied previously for modeling FL-induced desorption [8,13,[24][25][26][27]. As elsewhere [8,26,27], here we consider all six (coupled) molecular degrees of freedom of the diatomic adsorbate. The CO-surface interaction is described by the SAP potential, and surface phonons are included in an approximate way, by the generalized Langevin oscillator (GLO) model [28][29][30]. We treat electronic friction/damping in a firstprinciples manner, adopting the so-called local density friction approximation (LDFA) [31,32]. Stochastic forces arise, as in previous work [8,9], from a time-dependent electronic temperature which is determined from the two-temperature model [12]. Being widely used, we also critically assess the validity of different variants of the TTM. The MDEF calculations are supported by a few ab initio molecular dynamics including electronic friction (AIMDEF) trajectories [33,34].
In the end, we wish to study which reactions and molecular motions (including vibrations) are induced by the femtosecond laser, and what underlying mechanisms are, e.g., phonon vs electron mediated. To address these issues, in Sec. II we describe models and methods adopted in this work, and in Sec. III which properties were analyzed in which way. Results will be presented in Sec. IV. A final Sec. V summarizes and concludes this work.

A. Laser-driven molecular dynamics with electronic friction (MDEF)
A molecular dynamics with electronic friction (MDEF) approach is used here for the laser-driven dynamics of a CO molecule on a Cu(100) surface. In all cases below, a coverage of 1/2 is considered as in the recent experiment [9], with one CO molecule per √ 2 × √ 2 cell, i.e., CO:Cu(100)-( √ 2 × √ 2), as shown in Fig. 1. Here a word of caution is in order. In the (AI)MD simulations with periodic boundary conditions, all CO molecules move in concert and not independent of each other. In a sense, therefore, the coverage 1/2 is somewhat formal: While correct for the interaction potential (see below), the dynamics itself is effectively for lower coverage, formally even → 0. Furthermore, due to some desorption, also the experimental coverage is probably not exactly 1/2 [9].
Specifically, we solve a Langevin equation with friction treated in the independent atom approximation (neglecting two-body scattering corrections [35]), Here k = 1, 2 for the two moving adsorbate atoms, C (k = 1) and O (k = 2), with r 1 = r C , r 2 = r O denoting atomic coordinates, and m 1 = m C and m 2 = m O being corresponding atom masses, respectively. The first term on the right-hand side (r.h.s.) 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 6D parametrization of a potential energy surface (PES) obtained from periodic DFT within the generalized gradient approximation including dispersive interactions (see below).
The second term on the r.h.s. causes vibrational damping via a friction force, where η el,k (r k ) is the friction coefficient for atom k, calculated with the local density friction approximation (LDFA) [8,31]. Details on the parametrization of η el,k are given in Appendix A.
The last term in Eq. (1) accounts for hot-electron mediated excitation, by random forces R el,k (t ), calculated as white noise from the friction coefficients and the time-dependent electronic temperature T el (t ) as explained in Appendix A and Ref. [8]. T el (t ) is calculated from the TTM [12] to be discussed shortly.

B. Two-temperature model and its modification
As in previous work [8], the creation of hot electrons in the metal after laser pulse excitation is modeled with the help of the TTM [12]. Accordingly, time-dependent electron temperatures T el (t ) and phonon temperatures T ph (t ) are calculated from two coupled thermal diffusion equations which read (when neglecting phonon heat conductivity) S(z, t ) is a laser source term, where z is a vertical position within the metal, relative to the surface. We assume a Gaussian temporal pulse shape, characterized by a full width at half maximum (FWHM), a laser wavelength λ, and a laser fluence F , as described in Appendix B. In agreement with experiments done in Ref. [9], here we choose λ = 400 nm, FWHM = 150 fs, and we consider two fluences, F = 50 and 170 J/m 2 , respectively. For a critical look on (further) approximations made within the TTM, which we test here, a few details should be mentioned. Notably, the C el,ph in Eqs. (2) and (3) are electron and phonon heat capacities, κ el is the electronic heat conductivity, and g is the electron-phonon-coupling constant. The choice of the latter is described in Appendix B. The phonon heat capacity C ph is calculated using the Debye model [36] for the Cu crystal, also described in Appendix B. The T el dependence of the electronic heat capacity C el is usually approximated in a free-electron, Sommerfeld model as linear, where γ el is the specific heat capacity, also given in Appendix B. As shown there, this linear approximation does not work particularly well for Cu, at the electronic temperatures to be considered here. This is why in the following, a more general, nonlinear model developed in Ref. [37] is used instead, which takes band structure effects into account and which is also described in the Appendix. Using this nonlinear C el (T el ) model, which we refer to as the "semilinear model", we obtain indeed, lower electronic temperatures from the TTM as demonstrated in Fig. 2(a) (blue curve).
In the context of the TTM, another free-electron-like, linear approximation pertains to the electronic heat conductivity κ el in Eq. (2), where κ RT is the heat conductivity at room temperature (300 K), with a numerical value given in Appendix B. Also this approximation becomes questionable at high electron temperatures, due to enhanced electron-electron scattering. To improve on that, Petrov et al. developed a nonlinear, firstprinciple version of Eq. (5), details of which are described in Appendix B. Using the nonlinear model the electronic temperatures increase now noticeably in comparison to the old, linear model for κ el , as one can see in Fig. 2 (a). As shown there, taking both nonlinearities into account, those in C el and those in κ el , the electronic temperature reaches higher peak values (close to 7000 K instead of about 6000 K in the example shown), and stays "hot" for a longer time, compared to the usually adopted, linear models [8,9]. This will affect the hot-electron driven Langevin dynamics. The improved, nonlinear TTM with electronic and phonon temperatures obtained for the two fluences considered here, as shown in Fig. 2(b), will be used throughout this work.

C. Generalized Langevin oscillator model
In previous work [8,26,27], the generalized Langevin oscillator (GLO) has been employed for FL-induced desorption dynamics to account for interactions with lattice vibrations (phonons). The GLO model is a very approximate model to include surface phonons. Accordingly, surface motion is described in terms of a three-dimensional (3D) harmonic oscillator. Dissipation and thermal fluctuations are modeled via a ghost 3D oscillator coupled to the surface oscillator, which feels friction and random fluctuation forces to model energy exchange between the surface atoms and the bulk thermal bath. Details of the model and GLO parameters for CO:Cu(100) are given in Appendix C. Below, to disentangle the effect of surface atom motion from hot-electron driven dynamics, MDEF calculations have been performed with and without the GLO, i.e., MDEF-GLO and MDEF, respectively.

D. The six-dimensional SAP potential
For V (r 1 , r 2 ) in Eq. (1), the SAP-PES of Ref. [11] for CO:Cu(100)-( √ 2 × √ 2) (coverage 1/2) was employed. This Here, and in the following, the time origin t = 0 is chosen to coincide with the maximum of the laser pulse. In the usually adopted "linear model" (black curve), we use linear approximations according to Eqs. (5) and (4) for both κ el and C el (T el ). In the "semilinear model" (blue curve), the temperature is calculated with κ el linearly approximated according to Eq. (5), but for C el (T el ) a nonlinear model is used [37] as described in Appendix B. Finally, in the "nonlinear model" (red curve) both C el (T el ) and also κ el are calculated by nonlinear models, the latter from Eq. (B4) with parameters from Ref. [38] as described in Appendix B. (b) T el and T ph curves, obtained with the nonlinear model, for the two fluences used in this work. The temperatures of the nonlinear model were used for all dynamics simulations below.
PES is an analytical fit to periodic DFT-GGA (Perdew-Wang) calculations which reproduces experimental vibrational frequencies and adsorption energies of CO:Cu(100)- . Basic features of the SAP-PES are summarized in Appendix D.
Here we only show in Fig. 3(a) the potential V (X, Y ), with all other coordinates (r, θ , φ, Z) being optimized, measured with respect to the minimum of the potential. The relation between atomic Cartesian coordinates x k , y k , and z k (k = C, O) and the coordinates X, Y, Z, r and angles θ ∈ [0, π] and φ ∈ [0, 2π ] defined in Fig. 1, are also summarized in Appendix D. In Fig. 3 (a) it is shown that adsorption on top of Cu is energetically most favorable, while bridge and hollow positions are transition states and local maxima, respectively. The adsorption energy in on-top position is 0.76 eV, with CO upright and C pointing down.

A. Equilibration and statistical sampling
In what follows, MDEF trajectories were run for laser pulses with λ = 400 nm, FWHM = 150 fs, and two laser fluences (50 and 170 J/m 2 ). Calculations were done with (MDEF-GLO) and without the generalized Langevin oscillator model (MDEF), respectively. For each case 20 000 trajectories were run, i.e., 80 000 trajectories in total. For test purposes, MD calculations were performed without electronic friction and electron temperatures artificially set to zero. An initial temperature of 100 K was chosen in agreement with experiment [9], and all trajectories started with CO in its most stable adsorption geometry, i.e., upright with C down on top of a Cu atom, cf. Fig. 3 (a). Each trajectory was first equilibrated to reach the starting temperature using the Langevin algorithm for 15 ps. After that time, a laser pulse was applied, modeled by the Gaussian intensity profile Eq. (B1) defined in Appendix B, with a FWHM = 150 fs, again in agreement with Ref. [9]. The pulse intensity was maximal at 300 fs after the end of the equilibration phase (the time origin used here). Postequilibration trajectories were run for up to 5.0 ps after the pulse maximum.

Probability distributions and coordinate expectation values
Trajectories were analyzed with respect to the ensemble average values of individual coordinates Q (t ). In particular, ensemble average values for Z and θ will be considered, along with respective 1D probability distributions P(Z; t ) obtained from histograms. We further will show lateral probability distributions P(X, Y ; t ) of CO molecules on the surface. For details of their calculation, see Appendix E. Additionally, we also show distribution functions P(Z, θ; t ) of CO molecules there.

Kinetic energies and adsorbate temperatures
We further determine mode-resolved, ensemble-averaged kinetic energies of the moving CO molecules, the latter given as E kin (Q) = 1 2 m Q ( dQ dt ) 2 , where m Q is the reduced mass of mode Q, and · · · denotes an ensemble average. While FL-induced processes are not necessarily expected to create equilibrium energy distributions within the molecule, it is still useful and in some cases quite appropriate [39] to define timedependent, mode-resolved adsorbate temperatures according to We will also use mode-averaged temperatures if several modes are combined in a set comprising n individual modes. One may check the validity of the concept of a temperature for adsorbate modes, by testing if the kinetic energy distributions P(E kin ) fit to a Maxwell-Boltzmann distribution.
For the cases below where we show adsorbate temperatures, this was indeed the case to a good approximation.

Desorption and diffusion
For many trajectories the motion of CO after FL-pulse excitation was rather "moderate" in the sense that no largeamplitude motions were enforced within the simulation time leading to vibrational frequency shifts only. In other cases, however, lateral diffusion on the surface or even desorption took place.
To quantify the latter, the desorption probability was defined as the fraction of trajectories which reached a Z value of 6 Å or larger. To characterize the former, as in previous work [8], lateral, time-dependent diffusion coefficients D(t ) were calculated from the mean squared displacements in the surface plane as described in Appendix F. For diffusion, only those trajectories were included in the analysis which did not lead to desorption.

C. Analyzing vibrational frequencies
Time-resolved vibrational frequencies of the C-O stretch vibration of adsorbed CO, after photostimulation, were the major target of Ref. [9]. To obtain those from our MDEF simulations, two different measures were used in this work.

Instantaneous normal mode analysis
First, following a procedure similar to Ref. [9], instantaneous normal mode (INM) analyses were performed based on the SAP PES, to calculate the vibrational frequency of the C-O stretch mode, along r. The INM analysis can be done in a static fashion, giving the C-O stretch frequency as a function of fixed lateral positions X and Y of CO, i.e., ω r (X, Y ), or it can be evaluated dynamically as a function of time and as an average over a swarm of trajectories ω r,av (t ). Details of how this was done are described in Appendix G. In Fig. 3(b) we show (static) vibrational frequencies ω r (X, Y ), expressed as wave numbers (cm −1 ) for the r mode obtained in this way. It is seen that when CO adsorbs on top of Cu (in upright position, C down), the frequency is 2060 cm −1 . This value is not far from the experimental value (2079 cm −1 ), and also close to the anharmonically calculated value when using the SAP potential of 2053 cm −1 [11] (see also Table II below). When CO moves away from the top site, the C-O vibration can become "harder" or "softer," respectively. It is the "hardest" if CO is around X = Y = d/4 (modulo sign), i.e., halfway between top and hollow positions, where ω r ∼ 2150 cm −1 . It is the "softest" at the hollow positions around X = Y = d/2 (modulo sign), where ω r ∼ 1850 cm −1 . These frequency shifts correlate with (slightly) shortened C-O bonds at around X = Y = d/4, and elongated C-O bonds around the hollow sites, respectively, as can be seen from Fig. 3(b) in Appendix D. They do not correlate strictly with the potential energy itself [Fig. 3(a)], however, which follows (with its circular shape) trends for the optimal molecule-surface distance Z and, to a lesser extent, the tilting angle θ , respectively [see Figs. 3(b) and 3(d)]. In fact, simple σ → donation/→ π * backdonation arguments which may imply a correlation between adsorption energy and C-O frequency, are only strictly applicable at the top site, and not where CO interacts with several Cu atoms and can adopt tilted conformations.

Time-resolved Fourier transformation of the velocity-velocity autocorrelation function
Another way to obtain information on the vibrations during dynamics is from the Fourier transform of the velocityvelocity autocorrelation function. The latter is calculated in this work, in normalized form, as from an average over all (nondesorbed) trajectories. Here is the 6D velocity vector of the CO molecule (obtained with the Cartesian atomic coordinates for C and O). Furthermore, t is a reference time, chosen from a time grid along the trajectories, and τ is the (ordinary) propagation time.
A time-resolved vibrational spectrum (vibrational density of states) is then obtained from a Fourier transform as Since the velocity vector contains all six coordinates of the CO molecule, Eq. (9) gives frequencies for all six modes. They are calculated in a time-resolved fashion, and anharmonically, however, in classical approximation.

A. Hot-electron mediated heating of the adsorbates
When applying laser pulses with fluences of 170 and 50 J/m 2 to our model of CO:Cu(100)-( √ 2 × √ 2) as described above, the metal electrons and phonons are heated according to the TTM, as demonstrated in Fig. 2(b). According to Eq. (A2), then, random forces set the CO molecule in motion.
In Appendix H, Fig. 12, we show for two different models, MDEF-GLO and MDEF, the time evolution of the total kinetic energy of CO, as a result of energy uptake from the hotelectron bath, for both fluences. The figures give energy-and time-resolved distributions, and the time-resolved averaged kinetic energy E kin (t ). From there we see that before the laser pulse, when T = 100 K, narrow kinetic energy distributions around 3k B × 100 K = 0.026 eV are found, with k B being Boltzmann's constant. Thereafter, the kinetic energy increases on a picosecond timescale before slowly decreasing again, with energetic distributions which first broaden and then narrow. This is a consequence of hot-electron mediated excitation and the subsequent vibrational relaxation due to electronic friction. The maximum average kinetic energy is around 0.25 eV for MDEF/F = 170 J/m 2 at around 2.5 ps [Fig. 12(b)], corresponding to an energy uptake of about 0.22 eV. For the lower fluence the corresponding maximum average kinetic energy is around 0.10 eV, i.e., the uptake only roughly 0.07 eV [Fig. 12(d)]. In the MDEF-GLO models [Figs. 12(a) and 12(c), respectively], the energy gain is smaller than for the MDEF models, because the surface phonons, described by the GLO, reabsorb energy from the hot-electron stimulated adsorbate.
In Fig. 4 we plot, for the same cases as in Fig. 12, how the energy uptake translates into "adsorbate temperatures" according to Eqs. (6) and (7). Considering first Fig. 4(b) (MDEF, F = 170 J/m 2 ), where only electronic excitation and decay channels were included, we see that all modes are heated quite substantially after laser impact, on average up to about 1150 K at around 3 ps after the laser pulse maximum. This temperature is still much lower than the corresponding maximum electron temperature [∼7000 K, see Fig. 2(b)], and occurs at later times. It is higher, however, than the corresponding phonon temperature, which reaches a value of below 500 K after ∼4 ps. Note that there are also some differences between various modes. In particular the C-O stretch (r) reaches a higher peak temperature than the other modes (∼1300 K). This reflects the mode-dependent coupling of the CO vibrations on Cu(100) to electron-hole pairs.
Comparing the MDEF and MDEF-GLO models at the larger fluence [Figs. 4(a) and 4(b), respectively], we note that the r mode is rather unaffected, while all other modes are cooled down quite a bit, notably the Z mode. This reflects the more efficient coupling to the surface phonons, as described by the GLO, which mainly acts as an energy sink for these latter modes. Quite similar behavior is found at the lower fluence [Figs. 4(c) and 4(d), respectively], however, in this case the acquired temperatures are substantially lower.

B. Hot-electron mediated adsorbate motion
By the laser pulses the CO molecules are set in motion. In Furthermore, in Fig. 11 of Appendix E, we show the corresponding two-dimensional probability distribution P(Z, θ ) for different times relative to the laser pulse for F = 170 J/m 2 within the MDEF-GLO model. From Fig. 11 we see that the P(θ, Z ) distribution broadens with time and shifts towards lower Z and larger θ values, respectively. This finding is in qualitative disagreement with the reduced-dimensional Langevin dynamics of Ref. [9], where only these two modes had been considered. According to Ref. [9], both Z and θ shift to larger values (and broaden) with time, because no lateral motion away from the top site (where Z becomes smaller) is possible in the 2D model. In Fig. 6 we plot the distribution of CO molecules in the (X, Y ) plane, P(X, Y ), at different times. It is seen that the mentioned shift in the P(Z, θ ) distribution (Fig. 11 of Appendix E), appears to be a consequence of the molecules moving away from the on-top position (X = Y = 0) in the course of time. In fact, at 0.5 ps after the laser pulse maximum, for example, the probability to find the CO molecules at the bridge positions (at X = ±d/2, Y = 0 or X = 0, Y = ±d/2), is higher than finding CO on top. By computing the probability to find CO in a circle of radius of 0.7 Å around the top as a function of time, one finds that molecules start moving out of the top region already after 0.1-0.2 ps, with the probability dropping to a value of about 1/3 of its original value after 0.5 ps (see Fig. 13 in Appendix I). This implies that on the short timescales found in the experiment in which a C-O redshift is observed [9], the lateral motion of the molecules is non-negligible. From Fig. 6 we see that also the area around hollow positions (at X = ±d/2, Y = ±d/2) starts to be populated in the course of time. We note further that at very long times after the pulse (e.g., 4.5 ps), in particular the hollow sites start to be depopulated again. We remark that large populations at around hollow and bridge sites are also a consequence of the slower motion of CO molecules there, due to a relatively high potential energy. As will be shown below, the lateral motion of CO molecules has implications for the expected, time-resolved C-O stretch frequencies.
Before considering vibrational frequency shifts during dynamics, we analyze lateral diffusion of CO in more detail. For this purpose we show in Fig. 7 diffusion coefficients as calculated from Eq. (F2), for both fluences, using the MDEF-GLO model in each case. As a result of time-dependent temperatures, the diffusion rates are also time dependent. The largest diffusion coefficients are found for F = 170 J/m 2 , at around 2.5 ps after the pulse, with D ∼ 6 × 10 −4 cm 2 /s. For F = 50 J/m 2 , the maximal D is around 1.5 × 10 −4 cm 2 /s. In both cases, the diffusion coefficients are much larger than typically found for CO diffusion on metals under thermal (room temperature) conditions: For CO/Ru(0001), for example, D ∼ 3 × 10 −9 cm 2 /s was reported in Ref. [40] (for a coverage of around 1/4).
In analogy to Ref. [8], for each time t when D(t ) is evaluated, we can also determine the adsorbate temperature T ads according to Eq. (7), using n = 2 and considering the Q 1 = X and Q 2 = Y modes. The (averaged, lateral) adsorbate temperatures are shown as green curves in Figs. 4(a) and 4(c), for the MDEF-GLO model corresponding to both fluences. Assuming a simple Arrhenius expression for diffusion coefficients, from Arrhenius plots lnD vs 1/T ads , we can then determine the prefactor D 0 and the activation energy E a for diffusion, respectively. The Arrhenius plots, together with linear regression curves, are shown in Fig. 7. The data are quite noisy, however, from the linear regression we still get D 0 and E a values which are listed in Table I. As seen from there, the Arrhenius fits lead to a spread in prefactors and activation energies, which is partially due to the oversimplifications made by the Arrhenius model, and partially due to noisy data. Still, we can extract an activation energy in the order of 60-90 meV. As described in Appendix D, the SAP-PES gives a classical barrier for the top→bridge→top diffusion of about 60 meV [11]. Here we find an activation energy around or somewhat above this value, due to the multidimensional dynamics. As noted in Appendix D, the experimental measured diffusion barriers for CO on Cu(100) are controversial, ranging between 30 meV according to Ref. [41] and about 130 meV according to Refs. [42,43] (for diffusion via bridge). Given that our MDEF(-GLO) calculations are purely classical, without zeropoint energies which lower the barrier (see Appendix D), the dynamically determined barrier is probably too large. That would support the lower barrier found in Ref. [41]. On the other hand, barriers obtained from GGA-DFT are typically too low and, therefore, we cannot exclude that barriers on the order of 100 meV are more realistic.
Besides laser-induced diffusion, also laser-induced desorption of CO is observed in our simulations, at least at the higher fluence. According to Table I, the MDEF-GLO model predicts a desorption probability for F = 170 J/m 2 of ∼0.2% within the simulation time, i.e., a low value. Within the MDEF (only) model, it increases to ∼7.4%. That is consistent with the different kinetic energy gain predicted by both models (cf. Fig. 12), the GLO tends to absorb energy which is then not available to overcome the binding energy of ∼0.76 eV. According to Table I at the lower fluence of F = 50 J/m 2 , desorption can be entirely neglected.

C. Vibrational dynamics: AIMDEF calculations
As mentioned above, recent pump-probe experiments [9] concentrated on the C-O vibrational frequency and its time evolution following the excitation with 400 nm pump pulses. To set the stage and consider more general aspects of (dissipative) vibrational dynamics, we first of all demonstrate the performance of AIMD calculations including electronic friction (the AIMDEF model mentioned above [33,34]). The approach is computationally demanding, which is why we present only a few "proof of principle" trajectories here.
As for MDEF, we consider a CO:Cu(100)-( √ 2 × √ 2) model with coverage one half (for a slab with six copper layers). Furthermore, the revPBE exchange-correlation functional [44] as in Ref. [20], a slightly modified friction model, and the nonlinear TTM for 170 J/m 2 as before were used. We ran AIMDEF trajectories with a modified version of the Vienna atomistic simulation package (VASP) [33,34,45,46], with either only the electron temperature included and the surface kept fixed, or with both temperatures T el and T ph included and (most) Cu atoms being allowed to move also. Further details are described in Appendix J.
Already with single AIMDEF trajectories basic features of the hot-electron (and hot-phonon) driven vibrational dynamics FIG. 6. Probability distribution functions P(X, Y ) at different times, averaged over 20 000 trajectories, for the MDEF-GLO model and a fluence of F = 170 J/m 2 . The contours enclose regions with probabilities to find CO with 25% probability (1st Q), 50% probability (2nd Q), and 75% probability (3rd Q). emerge. In Fig. 8(a) we show how, for the T el -only case, in response to the elevated electronic temperature, the CO molecule acquires kinetic energy (upper panel) and starts to move. This can be seen, for example, as an oscillation of the z coordinate of the C atom (second panel). Taking simply the z component of the velocity of C, v C z (t ), Fourier transforming it with sliding Hann functions (and squaring it, see Appendix J), a simplified version to the full autocorrelation function approach of Eqs. (8) and (9) can be realized, giving information on which vibrational frequencies are excited at which time. This information is shown in the lowest panels of Fig. 8.
For the T el -only case, the time-energy analysis of v C z in Fig. 8(a) exhibits characteristic frequencies (wave numbers) around 2000-2100 cm −1 and around 300 cm −1 , respectively, which we interpret as the C-O stretch and the CO-surface stretch (S) and/or frustrated rotational (R) vibrations. In fact, according to Table II, the experimental and theoretical fre-quencies, from the SAP potential according to Ref. [11], for these modes fall in the respective ranges. (We note that the revPBE frequencies are expected to be slightly different from the SAP frequencies; in fact, a normal-mode analysis with VASP of CO:Cu(100)-( √ 2 × √ 2)/six fixed Cu layers gives 2039 cm −1 for the C-O stretch, 321 cm −1 for S, 277 cm −1 for the R, and between 53 and 57 cm −1 for the T modes [18].) We further see that intensities and frequencies are time dependent, e.g., the C-O stretch mode gains intensity, with a maximum about 1200 fs after the pulse maximum, and then loses intensity again, albeit nonmonotonic in this case. From the decay a "lifetime" may be deduced which is between about 1 to 3 ps for several (five) trajectories studied, which is due to vibration-electron coupling and in reasonable agreement with experimental lifetimes [19]. We also observe dynamical frequency shifts for the CO mode from AIMDEF, which will further be analyzed in Sec. IV D. Also for the low-frequency portion of the spectrum, intensities and frequencies shift with time. Careful inspection reveals the occurrence of "overtones" at about twice the S/R frequencies (around 600 cm 1 ), and of very weak signals at 100 cm −1 , corresponding to the frustrated translational mode T.
If phonon modes and T ph are included, we obtain Fig. 8(b), for a representative trajectory. In this example, the kinetic energy (upper panel) gained is about twice that of T el -only case [ Fig. 8(a)], and the motion of the C atom is of larger amplitude and more irregular, due to coupling to moving Cu atoms (middle panel). That more rather than less energy is gained is not contrary to the MDEF+GLO of Fig. 12, since there only the CO molecule was analyzed. The time-frequency plot in Fig. 8(b) reveals excitation, decay, and reexcitation of the CO stretch, and a rich, complicated pattern in the lowfrequency region at around 500 cm −1 and below. This pattern is due to the low-frequency adsorbate modes, phonons, and coupled modes, respectively, which are not further analyzed here.

D. Vibrational frequency shifts from MDEF calculations
As mentioned above, a striking result of recent pumpprobe experiment [9], are rapid frequency shifts ω r of the internal CO stretch mode. Specifically, for the larger pump fluence F = 170 J/m 2 a decrease by 15 cm −1 was observed within less than 500 fs, followed by a phase in which the old frequency was recovered. In Ref. [9] this relaxation phase was interpreted as a renarrowing of the nonequilibrium distribution along θ and Z (in particular of the θ mode), which was supported by the accompanying two-mode MDEF calculations. For the F = 50 J/m 2 pulse, the redshift of the reversible CO frequency shift was smaller ( ω r ∼ 5 cm −1 ), occurring on a slightly longer timescale.
To unravel details of vibrational response of CO adsorbed on Cu(100), the MDEF-GLO trajectories were analyzed for both fluences, using the same 20 000 realizations and protocols of above. In Fig. 9 we plot the time-dependent frequencies as obtained with the time-dependent/windowed Fourier transformation of the full velocity-velocity autocorrelation functions according to Eq. (9). The time resolution is determined by the step size separating the reference times in Eq. (9), which is 5 fs in all cases shown. Figures 9(a) and 9(c) give the entire frequency range (from 0 to 3000 cm −1 )  Fig. 9, center values at t = −1 ps (the temperature is T = 100 K). Considering Figs. 9(a) and 9(c), we note that prior to the pulse, at t = −1 ps, for example, clearly four frequencies are seen which we attribute to the vibrations listed in Table II. The centers of the prepulse vibrational frequency distributions arising from the MDEF-GLO calculations are also included in the table. The fact that these frequencies appear already prior to the pulse is a consequence of the fact that at T = 100 K (the prepulse temperature), these modes are excited in classical mechanics (but not necessarily in quantum mechanics). The slight shift with respect to other SAP-derived values in Table II is due to (classical) dynamics and the (different) treatment of surface motion.
When the pulse is switched on, the same features as for the single AIMDEF trajectory are observed: Vibrational signals shift and change their intensity with time. For the high fluence [ Figs. 9(a) and 9(b)], we note that (i) the Z (or S) mode (at 400 cm −1 ) vanishes; (ii) the θ, φ (or R) modes broaden and redshift; and (iii) the r (CO stretch) mode shifts to the red, becomes weaker, and finally shifts back to closer to its original frequency. Concerning the C-O stretch mode, we note that its redshift is in qualitative agreement with the pump-probe experiment [9]. Quantitatively, from the blow-up of the C-O stretch region in Fig. 9(b), we see that the shift is too large (∼100 cm −1 in theory vs ∼15 cm −1 in experiment) and occurs at later times (∼2 ps in theory vs less than 0.5 ps in experiment).
For the lower fluence (F = 50 J/m 2 ), similar observations are made, but now the influence of the pulse on vibrations is generally less important. For the C-O stretch mode, the redshift is now ∼15 cm −1 , maximal at around 1.5 ps after the pump pulse maximum.
The comparatively large redshift of the C-O stretch mode is consistent with the lateral probability distributions P(X, Y ) reported in Fig. 6, with areas close to bridge and even hollow sites being populated by the action of the laser pulse. Since the latter show strongly softened C-O vibrations according to Fig. 3, this leads to the observed redshift. At longer times (e.g., at t = 4.5 ps in Fig. 6), one finds that bridge and in particular hollow sites start to be depopulated again, reducing the redshift.
While our MDEF-GLO simulations are qualitatively consistent with experimental findings, at least for the redshift of the C-O vibration quantitative differences remain, i.e., for the magnitude and the timescale. A few more observations should be mentioned. First, in Fig. 9(b) we show the (averaged) results of a dynamical, quantum mechanical INM analysis as an alternative to the time-resolved Fourier analysis of classical correlation functions. The light-blue and green curves shown refer to a case where the INM has been done for all nondesorbing trajectories as described in Appendix G (r unopt.), and to a situation where the r mode has been reoptimized (at fixed other five coordinates) before INM analysis (r opt.), respectively. For the two curves, a small frequency offset compared to the Fourier analysis is found, but, more importantly, a still too large redshift after photoexcitation when compared to experiment, with and without reoptimization. Second, however, when excluding configurations moving away more than 0.7 Å from the top site in our INM analysis with reoptimized r values, we find a redshift in C-O vibrations of ∼10 cm −1 at 0.4 ps, and ∼15 cm − 1 at 1 ps after the pulse maximum, for the 170 J/m 2 pulse-see the yellow curve top, r opt. That is, at least the magnitude of the shift appears now consistent with experiment [9].
Which conclusions can be drawn from the latter observation? The constrained INM analysis (yellow curve) gives a frequency shift ω r close to the one found in the 2D (Z, θ ) Langevin simulation of Ref. [9], for which the same SAP potential as here was used. Based on our 6D Langevin dynamics we conclude that hot-electron driven lateral motion of CO cannot be neglected (when using the SAP potential). There is the possibility, however, that either the motion of the CO molecules away from the top site is hindered in reality, and/or that CO away from on top is not VSF detectable. Indeed, concerning this second possibility, it is interesting to note that in Ref. [9] it was argued that the VSF signal came mainly from molecules close to the adsorption well. Concerning the first possibility, we stated above that diffusion barriers are still controversial and so it is also not clear how accurate the SAP diffusion barrier is. Furthermore, as also stated our periodic setup with a single CO per unit cell implies an effective, low coverage (→ zero, see above), whereas independently moving CO molecules could block adsorption sites to be visited by other molecules. In that sense, our current model can serve as a prediction for the low-coverage regime.
We close this section by relating our MDEF approach to the theoretical work done in Ref. [20]. There, density functional perturbation theory was used to analyze the harmonically idealized, C-O stretch mode for a periodic array of CO molecules, adsorbed on top. The formalism accounted for electron-phonon coupling up to second order, including electron-hole pair excitations (first order nonadiabatic coupling) and electron-mediated phonon-phonon coupling [18,20] (second order nonadiabatic coupling). It was shown that, whereas the linewidth changes of the internal stretch mode are dominated by electron-mediated phonon-phonon coupling, the frequency shift is governed by electron-hole pair excitations. As a consequence, the redshift and blueshift follow closely the increase and decrease, respectively, of the electronic temperature, in agreement with the experimental observations. The model of Ref. [20], by construction, describes nonadiabatic effects more explicitly than the friction model, thus better accounting for fine details of the frequency shift at subpicosecond timescales, but lacks lateral displacements and desorption-which may not show up in VSF spectra as argued in Ref. [9]. Therefore, both models nicely complement each other and, when taken together, give a more complete picture of the ongoing dynamics.

V. SUMMARY AND CONCLUSIONS
In the present work, time-resolved, hot-electron driven dynamics were studied by a molecular dynamics with electronic friction model accounting for all six molecular degrees of freedom, a potential energy surface based on DFT, electronic friction treated in the local density friction approximation, coupling to phonons by a generalized Langevin oscillator model, and hot-electron creation by an improved two temperature model which takes band structure effects into account. Our simulations on precomputed potentials allow for excellent statistics for timescales of interest in recent experiments [9], and beyond. By computing time-resolved frequency maps, instantaneous frequencies and probability distribution functions, we find that the CO molecules are quite mobile on the surface, probing also non-top sites on subpicosecond timescales. We see that all modes of CO on the surface play a role in laser-induced dynamics. In particular, lateral large-amplitude motion leads to diffusive behavior, with enhanced rates compared to diffusion under ambient conditions. Desorption takes also place, however, only in rare cases. Based on computed frequencies alone, a strong redshift of the C-O vibration is predicted, of ∼100 cm −1 for the larger fluence, maximal at around 2 ps after the pulse, and then followed by slow recoverage phase. This shift is due to the hot-electron driven population of sites close to the bridge and even hollow positions, which are characterized by weakened C-O bonds with "softened" C-O frequencies. While the redshift is in qualitative agreement with a 2D (Z, θ ) Langevin model [9], in our case the shift is larger and "slower," and clearly dominated by lateral motion (along X and Y ). While it is not fully clear at the moment if diffusion barriers are accurate enough with the chosen, SAP potential, we also argue that, for better quantitative agreement with experiment, it could be that molecules far outside the top region are largely "invisible" in VSF probe experiments. In order to test both of these possibilities, we suggest to use AIMDEF simulations to overcome certain restrictions of the current model (only one CO per unit cell, restriction to the SAP potential, quite approximate treatment of phonons), and to calculate VSF signals, frequencies, and intensities, explicitly, as done, for example, in Ref. [47]. From the current MDEF trajectories we can already predict that, if the molecules are not blocked in their lateral motion (e.g., at low coverage), and if they can be probed outside the potential well (e.g., by another method than VSF), then the redshift should be larger and more long lived than found in Ref. [9]. We also predict that low-frequency modes S, R, and T, not probed so far in experiments, should show an interesting time-energy behavior in response to hot electrons.
Apart from these findings, other aspects of the present work are worth mentioning. First, we found that while electronic excitation and relaxation channels dominate, the surface atom motion cannot be neglected and in fact, both act in concert. Here the GLO model serves well, however, it is clearly less accurate than explicit inclusion of phonons. As mentioned and demonstrated, explicit inclusion of surface motion can be done by AIMDEF calculations "on the fly". Second, we have seen that the TTM for copper has some shortcomings, when linear, free-electron or Sommerfeld-like, approximations are used. Third, we have suggested several ways of how to compute vibrational frequencies in a time-resolved fashion, for interpretation of pump-probe vibrational signals.
In final conclusion, despite being improvable in the analysis and fine details of the VSF signal, the present investigation gave microscopic insight in the femtosecond-laser driven, hot-electron mediated dynamics of CO adsorbed on Cu(100). We think that both the conclusions and the methodological framework can be transferred to a broader class of systems and materials, relevant for femtosecond-laser surface chemistry and spectroscopy.

APPENDIX A: ELECTRONIC FRICTION COEFFICIENTS AND RANDOM FORCES
Atomic electronic friction coefficients are calculated within the local density friction approximation (LDFA) [8,31]. There, friction is derived from scattering phase shifts of an atom in a free electron gas of density ρ, closely related to the Wigner-Seitz radius r s (r) = (3/(4πρ(r))) 1/3 . In this work we fit the LDFA friction coefficients by double exponentials η el,k (r) = At lower densities (corresponding to large r s values), we set the friction coefficients to zero. Note that in Ref. [8] slightly less accurate single-exponential fits have been used instead. Note further that the friction coefficients generally also depend on T el [14,49]. Here we neglect this effect, being aware that it may play a role at very high electronic temperatures.
Both ρ and r s are coordinate (r) dependent. The density is calculated for each point in space possibly visited by a moving atom, using gradient-corrected density functional theory. Specifically, ρ(r) was calculated for the naked surface, using a five-layer slab model (lateral lattice constant 2.556 Å, relaxed along the surface normal), the Perdew-Burke-Ernzerhof, and a (15,15,1) k-point mesh, using the Vienna atomistic simulation package (VASP) [45,46].
Given the atomic friction coefficients and the timedependent electron temperature calculated from the twotemperature model, the random forces R el,k (t ) entering the Langevin equation are calculated as Gaussian white noise with properties R el,k = 0 and according to the second fluctuation-dissipation theorem [50]. Details of the numerical realization of Eq. (A2) are described in Ref. [25].

APPENDIX B: PARAMETERS USED FOR THE TTM
1. The source term in the TTM. We assume for the source term in Eq. (2), Here τ = FWHM/(2 √ 2ln2) is related to the full width at half maximum (FWHM) of the pulse, t 0 is its peak time, ξ is the penetration depth (inverse of absorption coefficient), and F is the laser fluence. The penetration depth is ξ = 14.9 nm for copper for a laser wavelength λ = 400 nm [51], which was chosen in the experiments of Ref. [9] and here. In further agreement with Ref. [9], we also use as the FWHM of the pulse 150 fs, and we consider two fluences, F = 50 and 170 J/m 2 , respectively. The temperatures T el,ph (t, z) and also S(z, t ) depend on time t and also on the position z within the metal, relative to the surface. The temperatures (and source term) used in this work are for z = 0, i.e., on the surface.
2. The electron-phonon coupling constant g. The coupling constant g was approximated to be the experimental value at room temperature value (g = 1 × 10 17 W m −3 K −1 [52]). We take this value for all temperatures, also because the precise temperature dependence of g for Cu is still under dispute [37,53,54].
3. Phonon heat capacity C ph . As is usual practice, we calculate C ph from the Debye model [36] of the Cu crystal, as with ρ Cu being the density of the copper metal (ρ Cu = 8.96 g/cm 3 [55]), M Cu is the atomic mass of Cu (M Cu = 63.546 atomic mass units [56]), N A is Avogadro's constant, and T D is the Debye temperature of copper. For the latter, we take T D = 343 K, the value for Cu at 0 K [36]. The quantity x in Eq. (B2) ishc s k/(k B T ph ) in the Debye model, with c s the speed of sound in Cu and k a phonon wave vector [36].
4. Nonlinear electron heat capacity C el . The T el dependence of C el is usually linearly approximated as C el (T el ) = γ el T el , where γ el is the specific heat capacity, γ el = 98 J m −3 K −2 for copper [37]. This equation is obtained from the Sommerfeld expansion of the electronic free energy and only accurate for single-band metals like aluminum and/or at low (electronic) temperatures. For copper with its localized 3d band and a delocalized 4s band, nonlinear behavior may be observed. A more general expression for C el (T el ) is [37] where f (ε, μ, } is the Fermi-Dirac distribution with chemical potential μ and D(ε) is the electronic density of states at energy ε. In fact, as shown in Fig. 3(c) of Ref. [37], when calculating μ and D(ε) using periodic DFT-GGA and VASP, deviations from linearity are found at temperatures above T el ∼ 2000 K. Using the nonlinear C el (T el ) curve determined in Ref. [37], and splining, we obtain indeed lower electronic temperatures from the TTM as demonstrated in Fig. 2 (a). To obtain this curve, the tabulated values of the C el (T el ) curve of Ref. [37] (on an extended temperature range) were used, which can be found under: http://www.faculty.virginia.edu/CompMat/ electron-phonon-coupling/. There, values are given for temperatures above 300 K only. For temperatures below 300 K, the Sommerfeld model Eq. (4) was used instead in the present work. 5. Nonlinear electron heat conductivity κ el . In the linear approximation of κ el according to Eq. (5), we used for κ RT the heat conductivity at room temperature (300 K), κ RT = 401 W m −1 K −1 for Cu [55]. Also this approximation becomes questionable at high electron temperatures. For improving on that, Petrov et al. developed a first-principle version of Eq. (5), with [38] κ el (T el , T ph ) = 1 κ s,el (T el ) Here the dimensionless t el is t el = 6k B T el /ε F with ε F being the Fermi energy, and

APPENDIX C: THE GENERALIZED LANGEVIN OSCILLATOR (GLO) MODEL
In the GLO, surface motion is described in terms of a threedimensional (3D) harmonic oscillator (mass m s , coordinate vector r s ), and associated, diagonal 3 × 3 frequency matrix . Dissipation and thermal fluctuations are modeled via a ghost 3D oscillator (coordinate vector r g ), which feels friction and random fluctuation forces that model 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 (C2) respectively. Here η ph is a phononic friction coefficient, and R ph is Gaussian white noise calculated similar to above, there for the electronic case. In particular, Eq. (A2) holds here also, however, with η el and T el (t ) replaced by η ph and T ph (t ), respectively. T ph is the phonon temperature as obtained from the TTM model. To adopt the GLO to CO:Cu(100), the following parameters (frequencies based on experiments [57]) have been used here (all in atomic units): xx = gs,xx = ω x = 10.6 × 10 −4 , yy = gs,yy = ω y = 10.6 × 10 −4 , zz = gs,zz = ω z = 6.2 × 10 −4 . For the friction coefficient, η ph = 5.75 × 10 −4 was taken, leading to an equilibration (at fixed T ph ) on the timescale of several picoseconds.

APPENDIX D: DETAILS OF THE SAP POTENTIAL ENERGY SURFACE
For V (r 1 , r 2 ) in Eq. (1), the SAP-PES of Ref. [11] for CO:Cu(100)-( √ 2 × √ 2) (coverage 1/2) was employed. The PES is an analytical fit to periodic DFT-GGA (Perdew-Wang) calculations using a five-layer slab, an atomic orbital basis of triple-zeta with polarization functions quality, and the ADF-BAND code. Dispersive interactions were included in the fit via a special form of three-body terms. The potential has been shown to nicely reproduce experimental vibrational frequencies of CO:Cu(100)-( √ 2 × √ 2) and desorption energies. For example, the SAP desorption energy is 0.76 eV when CO adsorbs perpendicular on top of a Cu atom, with C pointing downward, the most stable adsorption site (experiment: 0.79 eV at low coverages). Furthermore, the anharmonically calculated CO stretching (2053 cm −1 ) and frustrated rotation (289 cm −1 ) fundamental vibrational wave numbers deviate from experiment by only −26 and 4 cm −1 , respectively [11].
The classical barrier for diffusion from top to top via bridge is 57 meV according to SAP. Including zero-point corrections, this value is reduced to 33 meV [11]. Experiments on measured diffusion barriers for CO on Cu(100) are controversial: While in Ref. [41], a barrier of 30 meV was found, in more recent papers based on 3 He spin-echo spectroscopy, a barrier of around 130 meV was suggested [42,43] instead (for diffusion via the bridge site).
To gain some further insight into the SAP-PES, we show certain features of it in the form of two-dimensional functions f (X, Y ) in Fig. 3, and in Fig. 3(a) of the main text. There and here, instead of Cartesian atomic coordinates x k , y k , and z k (k = C, O), 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π ] are used. These coordinates are illustrated in Fig. 1. The relations between curvilinear and Cartesian atomic In Fig. 3(a) in the main text we already showed the potential V (X, Y ), with all other coordinates (r, θ , φ, Z) being optimized. Figure 10(a) indicates that, not surprisingly, for adsorption atop the distance of the molecule's center of mass from the surface (Z) is largest with Z around 2.5 Å, while the molecule moves closer to the surface when approaching bridge and, in particular, hollow sites. At the same time, the CO distance r slightly increases from about 1.15 Å to about 1.18 Å, for adsorption in the hollow site [ Fig. 10(b)]. Furthermore, Fig. 10(c) shows that the molecule is perfectly upright (θ = 0 • ) when on top, while tilting sideways to θ ∼ 15 • on the bridge, and θ close to 30 • at hollow sites.

APPENDIX E: PROBABILITY DISTRIBUTIONS
In the main text, lateral probability distributions P(X, Y ; t ) of CO molecules on the surface are shown in Fig. 6. These 2D distributions were obtained with the fast kernel density estimation (fastKDE) method [58]. Prior to the estimation procedure the (X, Y ) data were projected onto the irreducible wedge (spanned by top, bridge, and hollow), and then mirrored to fill the whole unit cell. This was done to enforce the symmetry of the surface within the distributions P(X, Y ; t ).
In addition, we also computed distribution functions P(Z, θ; t ) for different times relative to the laser pulse, with a fluence of F = 170 J/m 2 . These are shown in Fig. 11.

APPENDIX F: COMPUTATION OF DIFFUSION COEFFICIENTS
To compute lateral diffusion coefficients, as in previous work [8], first mean squared displacements in the surface plane were calculated as where R(t ) = X (t ) 2 + Y (t ) 2 and R 0 = X (0) 2 + Y (0) 2 = 0 is the starting point of all trajectories. The squared displacement can be related to a (2D) surface diffusion coefficient D(t ), which is time dependent due to time-dependent temperatures:

APPENDIX G: INSTANTANEOUS NORMAL MODE ANALYSIS FOR THE CO STRETCH
In this Appendix we first demonstrate how static CO stretch frequencies which parametrically depend on the lateral position (X , Y ) of CO were determined. Accordingly, we (i) fix the two lateral coordinates X and Y , (ii) optimize all other coordinates θ , φ, r, and Z at a selected (X, Y ) point, (iii) set up a (2 × 2) Hessian matrix (at each point) for the r and Z modes, H i j = ∂ 2 V/∂Q i ∂Q j (with V being the 6D SAP-PES, Q i and Q j being r and/or Z coordinates), and (iv) diagonalize the Hessian. This way we obtain eigenvalues λ r and λ Z for both modes, of which we are for now interested in the r-mode only, to calculate vibrational frequencies according to ω r = √ λ r m r . Here m r is the reduced mass of C-O vibration in the 2D (r, Z) Hessian model, given as m r = √ cos 2 (α)μ 2 CO + sin 2 (α)M 2 CO , where α is the "mixing angle" appearing in the (2 × 2) rotation matrix S which diagonalizes the Hessian according to S T H S = diag(λ r , λ Z ). Furthermore, μ CO and M CO are the reduced mass of the (pure) C-O vibration and the mass of the CO molecule, respectively.
To perform dynamic, time-resolved INM analyses, the instantaneous frequency of the r mode is computed as a weighted, moving average separately for each MDEF trajectory i, ω r,av,i (t ). The moving weights are determined from a cos 2 function in time, with a FWHM of 50 fs. The moving averages of each trajectory are further averaged over all FIG. 11. Probability distribution functions P(Z, θ ) at different times, averaged over 20 000 trajectories, for MDEF-GLO and F = 170 J/m 2 . The contours enclose regions to find CO with 25% probability (1st Q), 50% probability (2nd Q), and 75% probability (3rd Q). The abscissa Z is the change in the center of mass height respect to its equilibrium value at the atop adsorption site.
(or a subset of) trajectories to give an overall, averaged, timedependent frequency ω r,av (t ): Here N (t ) is the number of trajectories being averaged over. This number depends on time, because, again, we only count trajectories which are not desorbed at time t.

APPENDIX H: KINETIC ENERGY DISTRIBUTIONS
In the main text, Sec. IV A, we refer to the time evolution of the total kinetic energy of CO, as a result of energy uptake from the hot-electron bath, for both fluences and both for MDEF-GLO and MDEF models. These are shown in Fig. 12 in this Appendix. The figures give energy-and time-resolved distributions, and the time-resolved averaged (over 20 000 trajectories) kinetic energy E kin (t ) (as lines).

APPENDIX I: FRACTION OF MOLECULES NEAR ON TOP AND TIMING OF DIFFUSION
In this Appendix, we show in Fig. 13 the fraction P top of trajectories which are near the top position (within a circle of 0.7 Å), for MDEF-GLO and F = 170 J/m 2 . Also, curves are shown which give further details on the diffusion behavior.
In the main text it was hypothesized that molecules far away from the on-top site, contribute only weakly to a SFG signal. In fact, it seems this conclusion is supported by a different aspect of the experiment: the experimental relative amplitude R of the SFG signal (cf. Fig. 3(b) of Ref. [9]), i.e., the difference between the signal strength for the pumped system vs the unpumped one, contains a feature at t = 0.2 ps which the 2D model of Ref. [9] was not able to reproduce.
Interestingly, this feature and the curve shape of the relative amplitude in general correlate well with the proportion P top FIG. 13. Fraction of trajectories which are: (P top ) around the top (i.e., within a radius of 0.7 Å); (P top,notdiffused ) around the top and did not cross a diffusion barrier yet; (P top,diffused ) did so at least once but are around the top again; (P diffused ) diffused at least once (without checking if on top or not). of trajectories that are near the on-top position (i.e., within a radius of 0.7 Å), also shown in Fig. 13.

APPENDIX J: AIMDEF CALCULATIONS AND THEIR ANALYSIS
The AIMDEF calculations including electronic friction were done for CO:Cu(100)-( √ 2 × √ 2) with coverage one half. The unit cell consisted of one CO molecule and six layers of Cu atoms, each containing two atoms per cell. We further adopted the revPBE exchange-correlation functional [44], an energy cutoff of 400 eV, and a 5 × 5 × 1 -centered, Monkhorst-Pack k-point mesh [59] to sample the Brillouin zone. The electron-core interaction was treated with the projector augmented-wave (PAW) method [60] implemented in VASP [61]. Periodic supercells were used, with fixed lateral lattice constants of 3.611 Å (optimized for bulk copper with revPBE), and a vertical lattice constant (perpendicular to the surface) of 25.28 Å. To get starting geometries for AIMDEF runs, the lowest Cu layer and all lateral positions of the Cu atoms were fixed, and vertical positions and all C and O coordinates (with CO on top of one of the Cu atoms in the first surface layer) were optimized. The same, nonlinear two-temperature model with the 170 J/m 2 laser pulse as before [cf. Fig. 2(b)] was adopted. All calculations were done with a modified version of the Vienna atomistic simulation package (VASP) [33,45,46].
In a first set of calculations (with one example shown in the main text), only the electron temperature was included and the surface atoms were kept fixed. In a second set of calculations (with one example shown in the main text also), both temperatures T el and T ph were included and, besides the CO molecule, also the uppermost four Cu layers were allowed to move. In the T el + T ph calculations, C and O only feel T el , while T ph acts on the upper three Cu layers (the lowest two Cu layers were fixed, and the remaining third-lowest layer was allowed to move, but not exposed to temperature). The initial phonon and electron temperatures were 100 K. Algorithms for the T el -only calculations are published in Ref. [33], while details for T el + T ph calculations are published in Ref. [34].
Since the standard LDFA approach with friction coefficients calculated from electron densities of a rigid, naked surface is expected to be inaccurate when the surface moves, the electronic density used was taken from a Hirshfeld partitioning analysis [62,63].
In total, we ran five trajectories for the T el only, and four for the T el + T ph case. Each trajectory was 5.3 ps long (starting at −0.3 ps, ending at +5.0 ps, the laser pulse maximum being at t 0 = 0), and time steps of 1 fs were used.
In the AIMDEF calculations, time-frequency maps (like the one in Fig. 8) were obtained from (squared) fast Fourier transformations (FFTs) of velocity components v A i (t ) for component i (= x, y, z) of atom A = C, O, after multiplication with Hann window functions. In the figure, only signals related to transforms of v C z (t ) are shown. In the FFT/floating Hann window method, the function to be transformed is multiplied with a "sliding" Hann window function H (t − t H,n ) with width w (=1000 fs in our case), centered at times t H,n = t H,0 + (n − 1) t H (with t H = 20 fs in our case), and covering the whole function to be transformed. Each resulting function is Fourier transformed. Note that, since we work with unnormalized functions v A i (t ), and since in the figure we use squared signals |S(ω, t )| 2 , the scale obtained in Fig. 8 is "arbitrary" and cannot be compared to the scales used in