Femtosecond laser driven molecular dynamics on surfaces: Photodesorption of molecular oxygen from Ag(110)

We simulate the femtosecond laser induced desorption dynamics of a diatomic molecule from a metal surface by including the effect of the electron and phonon excitations created by the laser pulse. Following previous models, the laser induced surface excitation is treated through the two temperature model, while the multidimensional dynamics of the molecule is described by a classical Langevin equation, in which the friction and random forces account for the action of the heated electrons. In this work, we propose the additional use of the generalized Langevin oscillator model to also include the effect of the energy exchange between the molecule and the heated surface lattice in the desorption dynamics. The model is applied to study the laser induced desorption of O$_2$ from the Ag(110) surface, making use of a six-dimensional potential energy surface calculated within density functional theory. Our results reveal the importance of the phonon mediated process and show that, depending on the value of the electronic density in the surroundings of the molecule adsorption site, its inclusion can significantly enhance or reduce the desorption probabilities.


I. INTRODUCTION
Laser driven photochemistry has proven to be a useful tool for promoting reactions at surfaces or even as a way to open new reaction channels not accessible by thermal activation [1][2][3][4][5].In particular, one important reaction is the photodesorption of a molecule from a metal surface.Generally, desorption on metals can be induced either by directly exciting the molecule (IR photons) or it can be substrate mediated (UV/Vis/Near IR photons).Among substrate mediated processes, one usually distinguishes between desorption induced by electronic transitions (DIET) and desorption induced by multiple electronic transitions (DIMET) [6].DIET is practically realized by using continuous wave or nanosecond-pulse lasers with low intensity, resulting in small desorption yields that increase linearly with laser fluence.In DIET on metals, the adsorbate captures a hot electron and forms a short lived excited state (negative ion resonance).After decaying to the electronic ground state, the adsorbate may gain enough energy and desorb.On the contrary, DIMET, which is the subject of the present study, is realized by intense femtosecond laser pulses.Such pulses are short in comparison to typical relaxation times of adsorbate excited states and, consequently, they can produce multiple excitations of an adsorbate that lead to desorption.DIMET results in relatively large desorption yields that increase superlinearly with laser fluence [1].
Different methods have been used to model DIMET [3][4][5].Several of them are the so-called excitation-deexcitation models, in which the system jumps between two or more electronic states (see review [3] for a complete list).However, these methods due to their complexity have only been applied to a reduced number of degrees of freedom.In this work, we use an alternative model [7][8][9][10] that permits treating all the molecular degrees of freedom.Instead of treating excited states explicitly, in this model the nuclear motion is classical in the ground state potential and all the electronic degrees of freedom are included via friction and associated fluctuation forces.The friction force accounts for the dissipation of the adsorbate energy on the surface by creation of low energy electron-hole pairs, while fluctuation forces represent the inelastic scattering of hot electrons on the adsorbate nuclei.The magnitude of the fluctuation force is obtained in terms of the temperature of the laser-induced hot electrons.This electronic temperature can be estimated from the properties of the laser pulse and the metal substrate.The first important ingredient of this model is an accurate ground state potential, which can be modelled with a range of methods with increasing accuracy and theoretical, as well as, computational complexity, starting with simple two body potentials up to accurate quantum chemistry methods.Early works that used the molecular dynamics with electronic friction model to simulate the laser-induced desorption were based on empirical potentials [7,8].Nowadays, one can obtain better accuracy and predictability by state of the art non-empirical theoretical methods.Particularly, a good balance between accuracy and computational com-plexity is achieved by density functional theory (DFT).This method, already at its semi-local level, is able to capture reasonably well both metallic delocalized states and molecular localized states and their interaction.
Ab initio molecular dynamics, in which DFT is used at each integration step to calculate the forces, keeps both the DFT accuracy and the full dimensionality of the problem.However, it is still computationally too demanding to treat low probability processes or even to run long time (more than few ps) dynamics.In this work we are interested in phenomena that typically demand both large statistics and long time dynamics.Therefore, it is more advantageous to follow an alternative scheme that consists in constructing the adiabatic potential energy surface (PES) from a large set of DFT energies.Several interpolation techniques have been developed to obtain PESs that preserve the DFT accuracy [11][12][13].The main prerequisite is to reduce the dimensionality of the problem to the molecular degrees of freedom, to decrease the computational demand, thus keeping the surface frozen.This means that in the case of diatomic molecules, the interaction of the molecule and the surface is described by a six dimensional (6D) PES.This approach has been successfully applied to study the dynamics of different molecules on different metal surfaces [14][15][16][17][18].
Recently, the laser-induced associative desorption of H 2 on Ru(0001) has been successfully modelled by using such DFT-based 6D PES [10].In that work, the metal surface is kept frozen and the laser excitation is only modelled by random scattering of hot electrons with the nuclei of the molecule.Here, we extend this model by allowing for lattice movement that enables us to incorporate laser-induced phonon excitations.The study of the effect of phonons in photodesorption, compared to that of electronic excitations, is one of the main objectives of this work.
As an example, we will employ this methodology to study the laser-induced desorption of O 2 on Ag(110).Due to the importance of oxygen adsorption on silver surfaces in the process of ethylene epoxidation, this system has been a subject of numerous studies [19][20][21][22][23][24].Using the corrugation reducing procedure (CRP) [11], we have recently constructed the first DFT-based 6D PES for O 2 on Ag(110) and used it to study the dissociative adsorption process [25].The O 2 molecule can adsorb on Ag(110) on several adsorption sites that are characterized by different adsorption energies and electronic densities and, as such, it is an interesting model system.It gives us the possibility to investigate the importance of including the phonon excitations in the model for desorption from adsorption wells of different characteristics.
Photochemistry of O 2 on Ag(110) after substrate mediated photoexcitation under DIET conditions has been studied experimentally in Refs.[26][27][28].Photodesorption, photodissociation, and also CO 2 formation were observed there.To our knowledge, no experimental studies under DIMET (femtosecond laser) conditions have been carried out so far.As such, our investigation has a strong predictive character.
The paper is organized as follows.The theoretical model and its implementation are described in Sec.II.Application of this model to the desorption of O 2 from Ag(110) is examined in Sec.III.The main conclusions of the paper are summarized in Sec.IV.

II. THEORETICAL MODEL
A. Langevin dynamics for laser-driven nuclear dynamics at surfaces The response of a metal surface to the excitation generated by an ultra short laser pulse can be described by the so called two temperature model (2TM) [29].In this model, the equilibration between the electron and lattice heat baths with temperatures T el and T ph , respectively, is described by the following coupled diffusion equations, where C el is the the electron heat capacity, C ph is the phonon heat capacity, κ is the electron thermal conductivity, g is the electron-phonon coupling constant, and S(z, t) is the laser source term.In the regime of intense laser pulses that are studied here, metal electrons are rapidly heated to several thousands K due to the low electron heat capacity C el of metals.The formed hot electrons can either diffuse to the bulk [first term in the right hand side of Eq. (1)] or transfer heat to the lattice phonons [term g(T el − T ph ) in Eqs.
(1) and ( 2)].The heat source term S(z, t) is calculated for a metal film of thickness d by where I(t) is the adsorbed fraction of a laser pulse intensity and α −1 is the optical penetration depth.The latter is calculated from the laser wavelength λ and the imaginary part of the refractive index of the surface k as α −1 = λ/(4πk).Following Ref. [7], the laser-induced dynamics of the adsorbed molecule is modelled using a Langevin equation for each atom i in the molecule, where m i , r i , and η el,i are their corresponding mass, position vector, and electronic friction coefficient (for diatomic molecules as studied here i, j = 1, 2).The first term on the right hand side of the equation is the adiabatic force exerted on each atom of the molecule that originates from the interaction between the molecule and the surface.Here, this force is calculated as the gradient over the atomic coordinates of the 6D DFT PES.Note that the positions of the atoms are set relative to the position of the surface topmost layer r s .Surface movement r s (t) and the subsequent energy transfer between the molecule and the lattice atoms are calculated within the generalized Langevin oscillator (GLO) model [30][31][32] as explained below in Sec.II B. Non-adiabatic effects due to the coupling of the atoms in the molecule with the surface electronic excitations enter Eq. ( 4) through the friction force −η el,i dri dt and the random fluctuating force R el i , as described in Sec.II C.

B. Phonon excitations: The generalized Langevin oscillator model
As discussed in Sec.I, the use of a 6D PES in the dynamics equations limits the possibility of including the dynamical energy exchange between the molecule and the surface lattice, i.e., of allowing phonon excitations/deexcitations. One successful model that is able to keep the accuracy of a DFT based PES and at the same time provides a reasonable description of the surface movement is the generalized Langevin oscillator (GLO) model [30][31][32][33].In the GLO model, surface motion is described in terms of a three dimensional (3D) harmonic oscillator of mass m s with position vector r s and associated diagonal 3 × 3 frequency matrix Ω 2 .Energy dissipation and thermal fluctuations are modelled with the help of a ghost 3D oscillator with position vector r g .The mass and the associated frequency matrix for the ghost oscillator are also m s and Ω 2 .The equations of motion for the surface and ghost oscillators, which are coupled by the coupling matrix Λ gs , are the following, (6) The friction force −η ph drg dt models energy dissipation from the interacting surface atoms to the bulk thermal bath.Following Ref. [31], the friction coefficient η ph is calculated from the Debye frequency ω D as η ph = m s πω D /6.The random fluctuation force R ph , which models the heating of the surface atoms due to the thermal motion of the bulk atoms, is a Gaussian white noise with variance where ∆t is the time integration step, k B is the Boltzmann constant, and T ph is the time-dependent phonon (surface) temperature that is calculated in the 2TM model.The friction and random fluctuation forces are linked by the fluctuation-dissipation theorem to ensure that the surface atoms are coupled to a thermal bath of T ph .Oscillator frequencies (Ω 2 ) ii = 2ω 2 i and coupling matrix elements (Λ gs ) ii = ω 2 i are obtained from the surface phonon frequencies ω i (i = x, y, z) at the edges of the surface Brillouin zone, as implemented in Refs.[32,33].

C. Electronic non-adiabatic effects. Local density friction approximation
The Born-Oppenheimer (adiabatic) approximation, in which the electrons react instantaneously to the nuclear motion, is a cornerstone in gas-surface dynamics.Nevertheless, the existence of a nonadiabatic energy dissipation upon adsorption of gas species (atomic or molecular) on metal surfaces through electron-hole pair excitations is well established [34,35].Several methods have been used to model this dissipation mechanism [36][37][38][39][40].Among them, a method that has proven to be both accurate and suitable to perform multidimensional molecular dynamics is the local density friction approximation (LDFA) [38].In this model, electronic non-adiabatic dissipative effects are introduced in the dynamics via a friction force proportional to the velocity of the atom, as in Eq. ( 4).The friction coefficient η el is obtained in terms of the scattering of electrons by an atom inside a homogeneous free electron gas (FEG) as In this equation, n is the FEG density and k F is the Fermi momentum.The δ l (k F ) are the scattering phase shifts evaluated at the Fermi level corresponding to the potential induced by the atom in the FEG, which is calculated within DFT.The friction coefficient of Eq. ( 8) has successfully been used to calculate the stopping power of atoms and ions in metal solids and surfaces [41][42][43][44].Within the LDFA, the electronic density entering Eq. ( 8) is chosen at each point of the trajectory as that of the bare surface at the position of the atomic nuclei n(r i −r s ).The latter can be easily obtained from a DFT calculation as described in Sec.II D below.The LDFA has been applied to study the effect of electronic excitations in the dynamics of atoms and molecules on metal surfaces [38,[45][46][47][48][49][50].
The scattering of heated electrons with the adsorbate results also in a fluctuating force R el i [see Eq. ( 4)] that is connected through the fluctuation-dissipation theorem to the electronic friction force via the electronic temperature T el [7].Here, R el i is modelled by a Gaussian white noise with variance Note that R el i is usually negligible for the typical thermal surface temperatures used in gas/surface experiments and it can safely be neglected.However, this term gives a large contribution in case of the high T el that are obtained in the wake of the laser excitation, particularly, for adsorbates embedded in high electron density regions of the surface.

D. Implementation of the method
We start by solving the 2TM differential Eqs.(1)-( 3) to obtain T el and T ph as a function of time for the specific surface and laser pulse properties of interest.The calculated time dependent electronic and phonon temperatures are saved on a grid (in practice in steps of 0.05 ps) and used as inputs in the molecular dynamics calculations [Eqs.( 4)-( 9)].Another required input that is needed to obtain η el is the electronic density of the bare surface n(r).Here, n(r) is calculated with DFT and saved on a real space grid.
We perform classical dynamics calculations that neglect the zero point energy of the adsorbate.Each trajectory starts with the molecule resting in one of the adsorption wells.The initial position of the surface r s (and the corresponding momenta) are sampled by a conventional Monte Carlo procedure, such that they correspond to the initial surface temperature.The dynamics equations ( 4)-( 6) are integrated with a Beeman algorithm [51] as implemented in Refs.[32,52].At each integration step, the corresponding T el and T ph are obtained by a cubic spline interpolation.The electronic density at the position of each atom in the molecule n(r i ) is obtained with a 3D cubic spline interpolation of the DFT calculated bare surface density.
Using the same implementation that solves Eqs. ( 4)-( 6), one can also perform dynamics simulations that only include the electronic or the phonon contribution by setting, respectively, r s = 0 or η el,i = 0 in Eq. ( 4).In the following, the three types of calculations will be denoted as, LDFA+GLO, when including both the electronic and phonon contributions, LDFA, when including only the electronic channel, and GLO, when only phonons are included.

A. System properties: Results from DFT calculations
The ground state properties of O 2 on Ag(110) are described by a recently constructed 6D PES [25] that is obtained from the CRP interpolation of ∼ 25000 spinpolarised DFT energies.The latter were calculated with the Perdew-Burke-Ernzerhof (PBE) exchange correlation functional [53] as implemented in the VASP code [54,55] and with a plane-wave basis set energy cut-off of 400 eV.The surface was modelled by a supercell consisting of a (2×3) surface unit cell, a five-layer thick slab, and 14 layers of vacuum.Additional details are given in Ref. [25].The oxygen molecule on Ag(110), as predicted by DFT with the PBE exchange-correlation functional, features the four adsorption wells sketched in Fig. 1: In the SB well, the molecular center of mass (CM) is at the short bridge site with the molecular axis oriented along the [1 10] direction.In the LB well, the molecular CM is at the long bridge site with the molecular axis oriented along the [001] direction.In the H001 well, the molecular CM is at the hollow site with the molecular axis oriented along the [001] direction.In the H110 well, the molecular CM is at the hollow site with the molecular axis oriented along the [1 10] direction.In all these adsorption sites the molecule is parallel to the surface.Table I summaries the main features of each adsorption position, namely, the adsorption energy E a , the distance Z from the surface of the molecular CM, the O 2 internuclear distance r, and the value of the bare surface electron density at the position of each O atom, which is given in terms of the mean free electron radius r s .The H110 well is the closest to the surface and its adsorption energy is E a = −0.21eV.The H001 well is somewhat further from the surface, but with a larger adsorption energy of −0.24 eV.It is important to remark that all these values correspond to the results obtained with the frozen surface 6D PES that is used in our simulations of Sec.III C.However, we have checked that if the surface is allowed to relax, the H001 and H110 adsorption energies increase (in absolute value) to a similar value of around −0.36 eV.Note that the DFT-PBE description of the H001 and H110 adsorption wells seems to be in good agreement with experimental observations [22,23,56,57].
Adsorption wells at the bridge sites (LB and SB) are further away from the surface compared to the hollow wells.Considering the Z values of the bridge wells, one could be tempted to assign them to the measured physisorption wells [21,22,58].However, Table I shows that DFT-PBE predicts too large adsorption energies in these wells to be considered as physisorbed states [24].In spite of this, our study is meaningful since it uses a state of the art PES for O 2 on Ag(110).Additionally, the use of a PES for a system that presents several adsorption sites with different characteristics is advantageous for a theoretical study over systems with just one adsorption site.Having one system with several wells gives us the opportunity to more clearly study the dependence of the results on the properties of the wells, such as the adsorption energy, the distance from the surface, and thus, the electron density in which the molecule is embedded.In this respect, our results can be predictive for systems in which adsorption wells with similar characteristics exist.

B. Computational details
Our simulations are performed for laser pulses of Gaussian shape with 800 nm wavelength, 130 fs of full width at half maximum (FWHM), and absorbed fluences in the range F = 50 − 200 J/m 2 .Laser pulses with these properties were used in desorption experiments performed on other systems [2,59].The laser-induced T el and T ph are calculated using the following material constants for Ag: C el = 63.3J/m 3 K, κ = 429 W/mK, γ = 63.3J/m 3 K, and k = 5.29 [60][61][62].The phonon heat capacity C ph is calculated in the Debye model, with Debye temperature T D (Ag) = 225 K.The metal slab thickness d in Eq. ( 3) is set to 0.5 µm.We have checked that with this d-value the calculated T el and T ph are well converged.
The electronic friction coefficient entering Eqs. ( 4) and ( 9) as a function of the embedding density is given by η el (r s ) = 1.365 r −1.828 s e −0.082rs + 50.342 r 0.491 s e −2.704rs , (10) where both r s and η el are in a.u.This function fits the friction coefficients of an oxygen atom calculated for embedding FEG densities varying in the range r s = 1−6 a.u.This range covers all the electronic density values that are relevant in our dynamics.
In all the calculations presented below, the initial surface temperature is set to T ini = 100 K. To enable the thermalization of the molecule prior to the laser excitation, the laser pulse is turned on after 1.5 ps, thus keeping the initial temperature constant during this time interval.However, we have checked that the results of the dynamics do not depend on this thermalization time and that the laser pulse could be turned on at the beginning of the dynamics calculation without altering the final results.We have also checked that the largest integration step that can be used keeping the results of the dynamics stable is 1 fs.In all cases the integration time is 50 ps and the instant t = 0 corresponds to the start of the trajectory calculation.
As an outcome of our dynamics we consider that a molecule has been desorbed when its center of mass arrives at 6 Å from the surface and its velocity direction points away from the surface.We also distinguish another possible exit channel, dissociation, if the interatomic distance r is larger than 2.5 Å with positive radial velocity.

C. Results and discussion
Calculated desorption yields Y as a function of the laser fluence F are shown in Fig. 2 for the four different adsorption wells.These values have been obtained from the number of desorption events out of 30 000 trajectory calculations performed for each laser fluence and initial adsorption position.Characteristic super-linear desorption yields, which follow a power law Y = aF n with n > 1, are observed for the four wells.The values of the exponent n are in the range 2.6 − 5.8.These values are similar to those obtained for other systems [3], such as CO/Cu(100) with n = 5 − 8 [8,65], CO/Pd(111) with n = 7 − 9 [66], NO/Pt(111) with n = 6 [1, 67], O 2 /Pt(111) with n = 6 [68,69], O 2 /Pd(111) with n = 6−9 [70], and associative desorption of H 2 /Ru(0001) with n = 3 [10,59].Independent of the considered laser fluence, the highest desorption yields are obtained for H110, followed by H001, LB, and SB.The exponent n of the power law is also different for each well, its value decreasing from n = 5.8 for desorption from the SB well down to n = 2.6 for desorption from the H110 well.Both results can be mostly related to the differences in the adsorption energies of the different wells (see E a in Table I).The highest desorption yield and lowest exponent correspond to the well with the lowest adsorption energy and vice versa.However, the adsorption energy itself is not the only property ruling the desorption process.The LB and H001 wells have the same adsorption energy (−0.24 eV), but the yields are consistently larger for desorption from the H001 well than from the LB well.As shown below, this effect is related to the different mechanisms that rule desorption from the hollow and the bridge sites.
Figure 3 shows the time dependence of T el and T ph as obtained from the 2TM for F = 200 J/m 2 in comparison with the time evolution of the desorption rate from each of the adsorption wells.There are remarkable differences between the bridge wells (LB and SB) and the hollow wells (H001 and H110) observed not only in the magnitude of the desorption rates, but also in their time evolution.The desorption rates for the hollow wells seem to follow the time evolution of T el , but with a delay of around 3.5 ps.In contrast, the desorption rates from the bridge sites do not seem to be very much affected by the high increase of T el at short times.In these cases, the highest values of the desorption rates occur at longer times, once T el and T ph are equilibrated.It is worth to mention that the desorption rate from SB seems to follow the time evolution of T ph , but also with a certain delay.On the one hand, these observations suggest that desorption from the hollow sites is mainly an electron mediated effect, where the energy transfer from the electrons excited by the laser pulse to the adsorbed molecule plays a dominant role.On the other hand, these results also suggest that the heating of electrons is not that important for desorption from the bridge sites and that the laser mediated phonon excitation is the relevant mechanism in this case.In order to confirm these ideas and gain further insight in the relative importance of the .Note the different scales used for the bridge sites (upper panels) and the hollow sites (lower panels).Electron (orange line) and phonon (blue line) temperatures calculated from the 2TM are also shown (left ordinate).The electronic temperature peaks at values T el > 6000 K (see Fig. 4).The histograms are obtained by counting desorption events in intervals of 1 ps.
electron and phonon mediated mechanisms, we have performed the two additional types of calculations described in Sec.II D above, in which only the effect of either the heated electrons (LDFA) or heated phonons (GLO) is included in the desorption dynamics.The desorption yield obtained from the four adsorption wells for two different laser fluences, using the three different models (LDFA, GLO, and LDFA+GLO), are given in Table II.Additionally, the desorption rates for laser fluence F = 200 J/m 2 calculated with the LDFA and GLO models are shown in Figs. 4 and 5, respectively.The new LDFA and GLO results confirm the ideas inferred above.In the case of desorption from the hollow sites, the LDFA yields and rates are significantly larger than the GLO ones, while the opposite behavior is observed for desorption from the bridge wells although the differences between GLO and LDFA are smaller in these cases.Focusing on the LDFA calculations, it is clear that the desorption yields (Table II) and rates for the bridge wells (Fig. 4) are reduced to marginal levels as compared to the ones obtained for the hollow sites.However, Fig. 5 and the GLO values in Table II show that the phonon-mediated contribution to desorption is rather similar among the four wells.In fact, the small differences we observe seem to be correlated with the differences in adsorption energy.Thus, the lowest yield  Note the different scales used for the bridge sites (upper panels) and the hollow sites (lower panels).
corresponds to the SB site, the one with the largest E a , and the largest yield to the H110, the one with the lowest E a .The intermediate cases represented by the LB and H001 sites, which have the same E a , show very similar desorption yields.The absence of a similar one to one correspondence between E a and the LDFA yields points to the electronic-mediated mechanism as the one responsible for removing that correlation in the LDFA+GLO yields, since in both cases the largest to lowest values for desorption follow the order H110, H001, LB, and SB.Yet, it remains to be understood what property (together with E a ) rules the efficiency of the electronic mechanism.
As explained in Sec.II C, the electronic contribution to desorption is determined within the LDFA description by the value of the bare surface electron density at the position of each adsorbate (in our case the O atoms).The density profile along the plane normal to the surface that contains the molecule is shown in Fig. 6 for each of the adsorption configurations, together with the corresponding O atom positions.The inset shows the friction coefficient of one O atom as a function of the electronic density.Clearly, the embedding electron density is higher when the molecule is adsorbed on the hollow than on the bridge wells (see also Table I).This nicely fits with the results we have obtained.When T el is high, the fluctuation forces acting on O 2 are correspondingly larger if adsorbed on the hollow wells than if adsorbed on the bridge wells.Therefore, despite the similar adsorption energies of H001 and LB, desorption is more efficient from the former because of the larger embedding density.Further insight regarding the competition between the electron-and phonon-mediated mechanisms can be gained by comparing the LDFA and GLO results to those obtained with the LDFA+GLO simulations.Thus, the LDFA model predicts larger yields for the electrondominated desorption cases (H110, H001) than the LDFA+GLO.The reason is that the adsorbed O 2 , being efficiently heated during the initial time interval in which T el is high, reaches temperatures larger than T ph at least during this period.Therefore, when surface motion is also included in the dynamics, the surface takes energy from the electronically heated molecule and the desorption probability is reduced in respect to the ideal case in which no surface motion is allowed.In the case of the GLO simulations, the yields are slightly larger than the LDFA+GLO ones for the LB site, suggesting energy uptake by the electronic system, while for the SB well the LDFA+GLO yields almost coincide with the sum of the GLO and LDFA values.
In the following, we analyze the characteristics of the molecules desorbed from the different wells.The corresponding angular distributions of the desorbed molecules are shown in Fig. 7.These distributions are rather symmetrical around a desorbing angle of 45 • relative to the surface normal for all the adsorption wells.Nevertheless, in the case of molecules desorbed from the H110 site a slight tendency to desorb into directions closer to the surface normal is observed.We fit the obtained angular distribution to the velocity integrated flux-weighted Maxwell-Boltzmann distribution [10,71] that gives P (α) = (β + 1) sin α cos β α.The parameter β is a measure of the alignment of the desorption flux.For large values of β the flux is aligned normal to the surface and the distribution is narrow, while β = 1 corresponds to a cosine distribution.As seen in Fig. 7, β is practically one for desorption from the H100, LB, and SB wells and it is somewhat larger than one for the H110 well.These (small) values contrast with the values β > ∼ 3 obtained for the associative desorption of H 2 from Ru(0001) [10].In that case, the deviation from the cosine distribution was explained by the presence of a late barrier towards desorption, causing a channeling effect and narrow angular distribution.However, in our case, the potential energy defining the molecule-surface interaction is monotonically increasing from the wells to the vacuum region (see Fig. 5 of Ref. [25]), which results in β ∼ 1.
Next, we analyze how the energy of the desorbing molecules are partitioned in translational and internal (vibrational and rotational) degrees of freedom.The dependence of the translational, vibrational, and rotational energies of the molecules desorbed from the four adsorption wells as a function of the laser fluence is shown in Fig. 8. Equipartition of energy between the different degrees of freedom of a free diatomic molecule means that the values of its translational, vibrational, and rotational energies are ordered according to the ratio 3 : 2 : 2 [72].Figure 8 shows that only in the case of LB, where desorption is dominantly phonon-mediated, the ideal thermal desorption is approximately fulfilled.Deviations from the ideal ratio are already observed for desorption from the SB well, which is also a phonon-dominated process, and more clearly for H001 and H110 (both electron dominated).In these three cases, and for all the laser fluences, the translational energy is the largest and the rotational energy is the lowest.Within a good approximation, a linear increase with the laser fluence of the vibrational and translational energies of the desorbed molecules is observed.This is considered to be one of the hallmarks of DIMET [5].Finally, it is worth noting that we also observe few dissociation events for molecules initially adsorbed in the H110 well.This is a very unlikely process that has only been observed at the highest fluences (F > 175 J/m 2 ) and with probabilities lower than 10 −4 .In these conditions it is not possible to perform a more detailed analysis of the process.Still, the occurrence of dissociation events is an interesting result considering that the energy barrier to dissociation is 0.57 eV [25], significantly larger than the well depth of −0.21 eV.

IV. SUMMARY
In summary, we have extended the approach of Ref. [10] to simulate the multidimensional dynamics of a molecule adsorbed on a metal surface excited by an ultrashort laser pulse by including surface movement (phonons) via the GLO model.This allows us to treat simultaneously the laser induced electron and phonon excitations and their effect on the dynamics of the eventually desorbing molecule.
Using this new approach we have studied the laser induced desorption of O 2 from Ag(110).An interesting feature of this system is that it possesses four distinct molecular adsorption wells.This enables us to study how the desorption mechanisms are connected to the properties of the adsorption configuration.In general, we find that the effect of the laser-heated phonons in this system cannot be disregarded.Importantly, the phonon contribution to the desorption yield can be either positive or negative depending on the adsorption site.More precisely, when the molecule is initially adsorbed on the bridge sites inclusion of phonons increases the desorption probability.In fact, for these sites, coupling of the molecule to the phonon excitations constitutes the main desorption mechanism.However, for molecules adsorbed on the hollow sites not only the electronic channel is the dominant mechanism, but inclusion of phonons reduces the desorption probabilities because they take energy from the excited molecule.The subsequent reductions of the desorption yields can be rather high, in the range of a factor 2 − 7, depending on the laser fluence.These observations are rationalized in terms of the distances from the surface at which the adsorption sites are located and the subsequent values of the electronic density in their surroundings.Hollow sites are closer to the surface than bridge sites and, consequently, in regions of higher electronic density.For this reason the electron channel dominates desorption in the former and the phonon channel in the latter.
Our results also suggest which desorption mechanism will be dominant in systems that present both physisorbed and chemisorbed species.Since physisorbed molecules are located in low electronic density regions their desorption behavior is expected to be similar to the one we obtain for the bridge sites, whereas for chemisorbed states our findings for hollow sites apply.
FIG. 1. (Color online) Sketch of the position of the molecule in the four adsorption wells as predicted by DFT: long bridge (denoted LB), short bridge (denoted SB), hollow along the [1 10] direction (denoted H110), and hollow along the [001] direction (denoted H001).First layer surface atoms are shown in light grey, while second layer atoms are shown in dark grey.

FIG. 2 .
FIG. 2. (Color online) Desorption yields Y from the four adsorption wells (shown with different symbols and colors) as a function of the laser fluence F .For every well the coefficient n is calculated by fitting the data to the equation Y = aF n .

FIG. 3 .
FIG. 3. (Color online)Desorption rates as a function of time for F = 200 J/m 2 from the four adsorption sites (right ordinate).Note the different scales used for the bridge sites (upper panels) and the hollow sites (lower panels).Electron (orange line) and phonon (blue line) temperatures calculated from the 2TM are also shown (left ordinate).The electronic temperature peaks at values T el > 6000 K (see Fig.4).The histograms are obtained by counting desorption events in intervals of 1 ps.

FIG. 4 .
FIG. 4. (Color online)Desorption rates as a function of time for F = 200 J/m 2 from the four adsorption sites calculated with the surface frozen (LDFA model) (right ordinate).Laser excitation of the surface is modelled only by the electronic temperature (orange line) given by the 2TM (left ordinate).Note the different scales used for the bridge sites (upper panels) and the hollow sites (lower panels).

FIG. 5 .
FIG. 5. (Color online) Desorption rates as a function of time for F = 200 J/m 2 from the four adsorption sites calculated neglecting electronic excitations (GLO model) (right ordinate).Laser excitation of the surface is modelled only by phonon temperature (blue line) as given by the 2TM (left ordinate).

FIG. 6 .
FIG. 6. (Color online) Contour plot of electronic densities, expressed in terms of the mean free electron radius rs, for the configurations of the four adsorption sites.Contour lines are separated by 0.25 a.u., in the range from 2 a.u. to 6 a.u. as shown by the color map.The positions of the oxygen atoms in the adsorption sites are shown with red dots.The inset shows the friction coefficient as a function of electronic density given in terms of rs.

66 ❍✶✶✵FIG. 7 .
FIG. 7. (Color online) Angular distribution of the molecules desorbed from the four adsorption wells for a laser fluence of F = 200 J/m 2 .Red squares show the results of the dynamics.The blue line is obtained by fitting the data to the following function: P (α) = (β + 1) sin α cos β α.

F
FIG. 8. (Color online) Partition of the energy into translational (black squares), vibrational (blue circles), and rotational (red triangles) degrees of freedom of the molecules desorbed from the four wells as a the function of laser fluence F .Values for low fluences and SB well are not shown due to the poor statistics.

TABLE I .
Properties of the adsorption states of O2 on Ag(110): Adsorption energy Ea, O2-surface distance Z, interatomic O-O distance r, and electronic density in which oxygen atoms are embedded (expressed in terms of the mean free electron radius rs given in atomic units, a.u.).

TABLE II .
Desorption yields from the four sites and for two different laser fluences calculated with the full model (YLDFA+GLO), the model in which the surface is frozen (YLDFA), and the model in which electronic excitations are neglected (YGLO).