Femtosecond laser induced desorption of H 2 , D 2 , and HD from Ru ( 0001 ) : Dynamical promotion and suppression studied with ab initio molecular dynamics with electronic friction

J. I. Juaristi,1,2,3,* M. Alducin,2,3,† and P. Saalfrank4,‡ 1Departamento de Física de Materiales, Facultad de Químicas (UPV/EHU), Apartado 1072, 20080 Donostia-San Sebastián, Spain 2Centro de Física de Materiales CFM/MPC (CSIC-UPV/EHU), Paseo Manuel de Lardizabal 5, 20018 Donostia-San Sebastián, Spain 3Donostia International Physics Center (DIPC), Paseo Manuel de Lardizabal 4, 20018 Donostia-San Sebastián, Spain 4Institut für Chemie, Universität Potsdam, Karl-Liebknecht-Strasse 24-25, D-14476 Potsdam, Germany (Received 11 January 2017; revised manuscript received 10 March 2017; published 29 March 2017)


I. INTRODUCTION
The study of femtosecond laser induced chemical reactions of adsorbates at surfaces constitutes an active field of research [1,2].In metals, the reactions are typically substrate mediated and initiated by the excitation of the metal electrons by the laser pulse.The energy is subsequently transferred from the electronic system to the adsorbates leading to reactions among them.The important role played by hot electrons in these processes is manifested in several characteristics such as the opening of new reaction channels that cannot be accessed by thermal activation, a nonequal energy partitioning among the different degrees of freedom of the eventually desorbing molecules, and large isotope effects [1][2][3][4][5][6].
One of the reactions that has been extensively studied is the recombinative desorption of H 2 and D 2 from a saturated Ru surface [1,[3][4][5][7][8][9][10][11][12][13].Experimentally, a large isotope effect was observed, with H 2 desorbing from a fully H-covered surface much more readily than D 2 from a fully D-covered surface.Moreover, also experiments with varying mixtures of the two isotopes H and D were performed.In this way, it was found that, whereas desorption of the heavier isotope D 2 is facilitated by the presence of the lighter one (H) at the surface, on the contrary, desorption of the lighter isotope H 2 was hindered by the presence of the heavier one (D) at the surface [3].The enhancement of D 2 yields in the presence of H was attributed to a phenomenon, called "dynamical promotion" of a surface reaction in Ref. [3].Similarly, one is tempted to call the suppression of H 2 desorption in the presence of adsorbed D, "dynamical poisoning" of a surface reaction.These notions generalize the usual, well-known effects of static promotion and poisoning of reactions at catalytic surfaces, by dopants that affect their electronic properties and reactivity.In the experiments of Refs.[3,4], intense ultrashort femtosecond laser pulses are used and the desorption process is induced by multiple electronic transitions (DIMET) via "hot electrons".In Ref. [3], dynamical promotion was suggested to arise from the fact that coadsorbed H is efficiently heated and transfers energy to nearby D atoms to trigger the desorption of D 2 .However, neither clear evidence nor a detailed dynamical picture of the proposed dynamical promotion (nor the dynamical suppression in case of H 2 desorption) could be given.It is the purpose of the present contribution to shed light on the underlying nonequilibrium, many-atom dynamics for this reaction.
An approach that has shown to be advantageous to model DIMET is so called molecular dynamics with electronic friction (MDEF) [11,[14][15][16][17]. Briefly, the motion of the desorbing species is described by solving a Langevin equation in the ground state potential energy surface, which can be efficiently calculated within density functional theory (DFT).The nonadiabatic coupling to the electrons excited by the laser pulse is included via an electronic friction force and corresponding stochastic forces that arise from a time-dependent electronic temperature T e (t).An efficient way to incorporate the nonadiabatic forces in multidimensional dynamics is the so called local density friction approximation (LDFA) [18].Within this method, the coupling to the electrons is obtained in terms of the value of the bare-surface electronic density at the position of the adsorbate atom.Finally, in order to obtain T e (t) in terms of the characteristics of the laser pulse and the properties of the substrate the two temperature model (2TM) is customarily used [19].Here we use the 2TM also, as in Ref. [11], being aware that the concept of an electron temperature in laser-induced processes is an approximation in particular on very short time scales [20].
In Ref. [11], this formalism was applied to study the dynamics of the laser induced desorption of H 2 and D 2 from either the fully H-or D-covered Ru(0001) surface, respectively.In that work, MDEF simulations in a precalculated six-dimensional DFT potential energy surface (PES) that included all the degrees of freedom of the desorbing molecule were used.In this way, reasonable agreement was obtained with several experimental observations such as the isotope effect, the dependence of the desorption yield on the laser fluence, and the energy partitioning among the different degrees of freedom of the desorbing molecules.However, in order to study the effect of having different mixtures of H and D, the number of adsorbates that one must include in the simulations is much larger than two.This means that the number of degrees of freedom involved (three per adsorbate atom) is too large to allow a calculation based on a precalculated PES.This problem can be tackled with the recently developed ab initio molecular dynamics with electronic friction (AIMDEF) method that is based on the LDFA scheme [21][22][23][24][25].
In this work we have extended the AIMDEF method in order to include random forces dependent on a timedependent electronic temperature.This allows us to integrate the Langevin equation calculating on the fly the adiabatic, friction, and stochastic forces with, in principle, an arbitrary number of adsorbates and with different H:D ratios.Our results show that, in agreement with the experiments, the presence of H on the surface favors the desorption of D 2 , whereas, on the contrary, the presence of D on the surface hinders the desorption of H 2 .Looking at the energy gained by the adsorbates as a function of time for different relative concentration of isotopes we shed light on the mechanism behind this effect.
The paper is organized as follows: In Sec.II we describe the theoretical model employed and the details of the simulations.In Sec.III we show and discuss the main results of our work.Finally, Sec.IV is devoted to summarize the main conclusions of the work.

II. METHODS: AIMDEF APPLIED TO DIMET
The femtosecond laser induced desorption dynamics is simulated by combining, on the one hand, the 2TM that defines as a first step the electronic temperature T e (t) describing the laser-induced hot electrons and, on the other hand, Langevin equations that describe as a second step the effect of the hot metal electrons on each adsorbate i, that is, where m i , r i , and η e,i are the mass, position vector, and electronic friction coefficient of the ith adsorbate.The first term on the right-hand side of the equation represents the adiabatic force that depends on the position of all (adsorbates and surface) atoms.The third term R e,i is the random fluctuating force that mimics the effect of the hot metal electrons on the ith adsorbate.This force is related to the electronic friction force (second term) through the fluctuation-dissipation theorem.Following previous works [11,[14][15][16], R e,i is modeled by a Gaussian white noise with variance where k B and t are the Boltzmann constant and the timeintegration step, respectively.The AIMDEF package [21][22][23][24][25], which relies on the LDFA [18] and is implemented in the DFT-based VASP [26] code, is here adapted to include and calculate at each integration step the time-dependent electronic temperature T e (t) and all the forces, i.e., also the Langevin equations ruling the dynamics of the adsorbates are solved ab initio and "on the fly".
All calculations are performed with VASP (version 5.2.12) using DFT and the generalized gradient approximation to the exchange-correlation functional by Perdew-Burke-Ernzerhof (PBE) [27].The electron-core interaction is treated with the projector augmented-wave (PAW) method [28,29].The PAW potential for Ru has eight valence electrons.Electron wave functions are expanded in a plane-wave basis set with an energy cutoff of 350 eV.Fractional occupancies are determined through the Methfessel and Paxton broadening scheme of first order using a width of 0.1 eV [30].The energy criteria for total energy self-consistency is 10 −6 eV.The initial clean and H(1 × 1)-covered Ru(0001) surface optimization as well as all AIMDEF simulations are performed using a -centered 3 × 3 × 1 Monkhorst-Pack grid [31] for the Brillouin zone integration.
The Ru(0001) surface is modeled with a periodic five-layer slab defined by a (4 × 4) surface unit cell and a supercell vector along the surface normal of 30.129Å.The Ru bulk lattice constant a = 2.707 Å is previously calculated for the stacking c/a = 1.59 Å [32] using a -centered 11 × 11 × 11 Monkhorst-Pack grid.After surface relaxation by either keeping the third layer fixed or by only relaxing the two topmost layers, the first and second interlayer distances are, respectively, reduced to d 1 = 2.07 Å and d 2 = 2.16 Å, to be compared to the experimental d 1 = 2.06 Å and d 2 = 2.17 Å [33].
Next, the H(1 × 1)-covered Ru(0001) surface is modeled using the same periodic slab and placing one H atom in each of the 16 fcc sites.The equilibrium configuration is obtained by allowing relaxation of the H atoms and the two first Ru layers.In agreement with other studies [34,35], the first interlayer distance increases upon H adsorption.More precisely, the interlayer distances in the H-saturated Ru(0001) surface are d 1 = 2.11 Å and d 2 = 2.15 Å, while the height of the adsorbed H atoms with respect to the Ru topmost layer is 1.05 Å.For the sake of comparison to the 6D PES calculated in Refs.[8,11] using the RPBE exchange-correlation functional, we have calculated the elbow plot V (r,Z) for two adjacent H atoms desorbing with its molecular axis oriented parallel to the surface.Taking the relaxed H(1 × 1)/Ru(0001) surface as the reference zero energy, the corresponding PBE energy barrier for desorption is E des ∼ 1.2 eV and it is located at r ∼ 0.79 Å and Z ∼ 2.1 Å, while the RPBE-DFT calculations of Ref. [11] yield E des = 1042 meV, r ∼ 0.77 Å, and Z ∼ 2.25 Å.Note that the PES of Ref. [11] is a modified version of the PES of Ref. [36], using a three-layer slab model (originally for a lower coverage of 1/4).In Ref. [7], Luntz and co-workers studied the same system at full coverage using PW91-DFT and a six-layer slab model.They obtained a minimum desorption barrier of about 1.15 eV that corresponds to a reaction path slightly different from the aforementioned 2D desorption configuration, for which their calculated desorption barrier is about 1.2 eV.Experimentally, a barrier of about 0.9 eV was given for a coverage close to one [37].
All AIMDEF calculations are performed keeping the Ru surface atoms fixed at their equilibrium position (frozen surface approximation).This approximation is reasonable in view of the small mass ratio m H /m Ru = 9.97 × 10 −3 as demonstrated in Refs.[23][24][25].Furthermore, previous MDEF simulations performed on the H(1 × 1)-and D(1 × 1)-Ru(0001) surfaces, which neglect the laser-induced phonon excitations, reproduce rather well the experimental data [11].Such a good agreement suggests that the dynamics of the adsorbates is mainly ruled by the excited electrons rather than by the excited phonons.Actually, the same conclusion was achieved by means of two-pulse correlation experiments in Ref. [3].In order to reduce the computational cost, only the first three (frozen) Ru layers are included in the AIMDEF simulations, while keeping the same supercell dimensions described above [the (4 × 4) surface unit cell and perpendicular vector of 30.129Å].We note that despite the H adsorption height being the same, the desorption energy barriers become smaller using the three-layer slab and PBE.In particular, the energy barrier to desorption for the aforementioned configuration of two adjacent H atoms desorbing parallel to the surface reduces to E des ∼ 0.88 eV.Assuming that, as it was the case in Ref. [7], the minimum energy desorption barrier is around 50 meV lower than the one found for this configuration, we can estimate a minimum energy desorption barrier in our simulation system of the order of 0.83 eV.Though accidentally, we note that these values for the desorption barrier are rather close to the experimental one (see above), at least if zero-point corrections are neglected, as done here.
In order to investigate the experimentally observed dependence of the H 2 , D 2 , and HD desorption yields on the D fraction of saturation coverage D , we performed AIMDEF simulations for the following H:D mixtures: (i) 16:0, (ii) 12:4, (iii) 8:8, (iv) 4:12, and (v) 0:16.The selected H:D mixtures represent, respectively, the following D fraction of saturation coverage D = 0.0, 0.25, 0.5, 0.75, and 1.0.For each D we run a total of 40 (classical) trajectories with the 16 adsorbates initially located at the bottom of each fcc adsorption well, i.e., the zero point energy of the adsorbates is neglected.Note that the initial conditions for each trajectory differ on the random force value.Additionally, in cases (ii)-(iv), the adsorbed H and D atoms are initially randomly distributed in each of the AIMDEF trajectories.
Simulations are performed for a laser fluence F = 60 J/m 2 and the laser pulse properties of Ref. [3], namely, a Gaussian shape with 800 nm wavelength and 130 fs of full width at half maximum (FWHM).The corresponding electronic temperature obtained from the 2TM for the Ru(0001) surface is plotted in Fig. 1.The Ru parameters entering the 2TM are those of Refs.[9,11,38] and are summarized in the caption of Fig. 1.The initial electronic temperature previous to the arrival of the laser pulse is 170 K as in the experiments of Ref. [3].However, by performing a few preliminary simulations to optimize the FIG. 1. Electronic temperature calculated from the 2TM using the implementation and the same material parameters of Refs.[9,11,38] for Ru(0001) and a Gaussian laser pulse of 800 nm wavelength, 130 fs FWHM, and fluence F = 60 J/m 2 .The material parameters are: g = 1.85 × 10 18 W/(m 3 K) (electron-phonon coupling constant), γ = 400 J/(m 3 K 2 ) (specific electronic heat capacity), κ 0 = 117 W/(m K) (thermal conductivity), T D = 600 K (Debye temperature), ρ = 12 370 kg/m 3 (material density), d = 15.6 nm (optical penetration depth at 800 nm).The laser pulse hits the surface at t = −100 fs when T e is equal to the experimental surface temperature of 170 K.The AIMDEF simulations start at t = 0. computational parameters we observed that the adsorbates hardly move during the first 100 fs after the application of the laser.For this reason, and in order to shorten the required computational time, our AIMDEF simulations start after these first 100 fs.In the following we define as t = 0 the instant at which our dynamics simulation start.In our frozen-surface approach, the phonon temperature T ph (t), which we also obtain from the 2TM, has no influence on the dynamics.The yields plotted in the left panel are normalized such that for pure coverages ( D = 0 and 1) they correspond to the absolute yield per shot.In other words, the values are obtained by dividing the total number of desorbing H 2 , D 2 , and HD molecules by the total number of trajectories run for each D (40) and by the number of molecules of one kind available for pure coverages (8).In order to show more clearly the effect of the relative concentration of the adsorbates, the middle panel shows the H 2 and D 2 desorption yields normalized to the values obtained for the pure H and D adlayers, respectively.

III. RESULTS OF THE AIMDEF SIMULATIONS
Let us start by discussing the isotope effect between isotopically pure coverages.For the value of the laser fluence that we consider in this work (F = 60 J/m 2 ) the MDEF simulations performed by Füchsel et al. [11] on a precalculated six-dimensional PES give a ratio between the H 2 and D 2 desorption yields of around 8 (see Fig. 4 in that reference).Since the isotope ratio is a very sensitive quantity, this result was considered to be in a quite reasonable agreement with The reason for the less pronounced isotope effect found here compared to Ref. [11] can be linked to the above discussed lower energy barriers towards desorption in our PBE three-layer model compared to Ref. [11].A lowering of the energy barrier to desorption is more critical for the less mobile isotope D and it favors more efficiently its desorption as a D 2 molecule.For this reason it is not surprising to obtain a smaller isotope ratio in our simulations compared to Ref. [11].More importantly, compared to the experimental yield ratio of ∼10, we see that our multidimensional AIMDEF calculations neither can fully quantitatively account for the experimental desorption yields.Closer analysis of our data suggests that the absolute yields per laser-shot Y H 2 are by a factor of ∼2 too high when compared to experiment (which is ∼0.16 [4]), while Y D 2 is by a factor of ∼7 too high.This explains the too small isotope effect in our AIMDEF simulations compared to experiment.In contrast, the isotope effect found in Ref. [11] is similar to experiment as mentioned above.However, absolute yields are too low there by more than one order of magnitude compared to experiment for both isotopes.The latter might be due to an overestimated desorption barrier, but also to the use of a smaller surface cell that might limit the interadsorbate energy exchange effect discussed below.Thus, considering that, on the one hand, our absolute values of the desorption yields are closer to the experimental ones than those obtained in Ref. [11] and, on the other hand, the isotope effect is better captured by the results of Ref. [11], it is very difficult to decide which value of the desorption barrier is the most realistic one.All in all, this shows that absolute theoretical yields are sensitive quantities that significantly depend on the accuracy of the DFT PES, what limits somehow the quantitative comparison between theory and experiments.In spite of it, our present AIMDEF simulations can perfectly capture in a qualitative manner how the mobility of the adsorbates depends on the relative H:D concentration at the surface and disentangle the mechanism behind such a dependence.
Before analyzing the results obtained for mixed coverages, let us first discuss the translational energy of desorbing molecules in the case of pure coverages, which is a quantity that is also measured from time-of-flight spectra [5].In our dynamics the statistics is not good enough to compute timeof-flight spectra.However, the average translational energy of the desorbing H 2 and D 2 can be determined as 474 and 390 meV, respectively.These results compare well with 490 and 415 meV that are the theoretical results obtained in Ref. [11].From time-of-flight measurements, a translational energy of 535 meV has been reported for H 2 [5].No result for D 2 for a 60 J/m 2 fluence was given in that reference, but the result for a lower fluence (50 J/m 2 ) was 336 meV.Therefore, it can be concluded that the agreement between our results and the experiments for the translational energy of the desorbing molecules is satisfactory.
The fact that our simulations account for the enhanced (reduced) desorption yield in the presence of H (D) adsorbates at the surface can be seen for instance in the normalized desorption yields shown in the middle panel of Fig. 2. In case the relative concentration of adsorbates at the surface had no influence in the desorption rate, the normalized yields for the 50% concentration should be the same for both H 2 and D 2 desorption.However, the normalized yield is significantly larger for D 2 at this concentration, showing the promotion effect due to the presence of the lighter H isotope.
The same conclusion can be extracted by looking at the solid lines shown in Fig. 2 (left and middle panels) and comparing them to the results of our simulations.The solid lines are obtained assuming that, as in the case of thermal desorption experiments [3], the desorption rates follow second order rate equations of the form In this approach, it is assumed that the rate constants k H and k D are independent of the concentration of H and D on the surface.We also assume that the rate constant for HD formation is the geometric mean of the "pure" rate constants.In order to obtain the solid lines, these equations are integrated for different initial values of [H] and [D] under the constraint [H(t = 0)] + [D(t = 0)] = 16, since 16 is the total number of adsorbates in our simulations.The values k H and k D are obtained by performing the integration for full H and D coverage, respectively, and fitting the results to the desorption yield values from the simulations.When solving Eq. ( 3) the final integration time is arbitrary, but it must be the same independently of the relative concentrations and equal to the one used in the calculations of k H and k D .Test calculations have shown that considering 16 adsorbates at the surface is indeed a very good approximation to represent a large surface area covered with adsorbates [39].
FIG. 3. Snapshots of one of the simulations for a 8:8 isotope mixture ( D = 0.5) that ends with two HD desorbing molecules.White and black spheres are, respectively, the nondesorbing H and D atoms, (labeled) yellow and blue spheres are the desorbing H and D atoms, and golden spheres are Ru atoms.The AIMDEF simulation time is indicated in each panel (t = 0 as in Fig. 1).
We observe in Fig. 2 that, with few exceptions, in general the curves obtained from the kinetic model lie above the data points in the case of H 2 desorption and below the data points in the case of D 2 desorption.These results are consistent with the experimental observation that the presence of D on the surface hampers the mobility of H, whereas, on the contrary, the presence of H favors the mobility of D. In other words, the observation supports the concepts of "dynamical promotion" [3] and "dynamical suppression" of surface reactions, respectively.
In order to understand why the relative concentration of isotopes at the surface influences the desorption yields, it is illustrative to follow the dynamics of the adsorbates.Figure 3 shows as an example the snapshots of one of our simulations for D = 0.5.The selected trajectory corresponds to a case in which two HD molecules are desorbed at the end of the simulation.Interestingly, we observe that all 16 atoms in the simulation cell abandon their adsorption sites and move at the metal surface upon their interaction with the heated electrons well before the first HD desorbs.The snapshots collected up to t = 550 fs show how the adsorbates collide frequently which each other until a molecule can eventually be formed and desorbs.Clearly the desorption of a molecule is not just caused by the excitation of the two adsorbates forming the molecule in an otherwise static environment.On the contrary, the picture that emerges is that the desorption of the molecule is the result of a collective excitation of adsorbates that form a heated environment in which the desorbing molecules are created.Since the heating of the adsorbate system comes from the coupling of the adsorbates to the electronic system, this process is adsorbate mass dependent.As a consequence, the degree of heating of the environment and, ultimately, the desorption rate depend on the relative concentration of isotopes at the surface.Another manifestation of this collective character of the desorption process is that the desorbing molecule is not always formed by two first nearest neighbors.For instance, in the simulation shown in Fig. 3 the atomic constituents of the two desorbing HD come from second and fourth neighboring sites.The right panel of Fig. 2 shows the distance distribution of the adsorbates forming the desorbed species for each D .Even if in all cases the distribution is dominated by the first nearest neighbors, the contribution of the second nearest adsorbates to desorption is also important, in particular for coverages D < 0.75, which are characterized by an important presence of extremely mobile H atoms.
The time evolution of the kinetic energy of the adsorbates provides direct information on the heating of the adsorbate system.Figure 4 shows the mean kinetic energy E kin of desorbing and nondesorbing adsorbates as a function of time.The mean value is calculated as an average over all trajectories contributing to a specific event [40].The results for different H:D mixtures are plotted in different panels, except for the top-left panel that shows E kin for the two isotopically pure coverages ( D = 0 and 1).Precisely this panel shows that H gains energy more rapidly than D, as expected.Also for mixed H:D coverages we observe that the H adsorbates gain energy more rapidly than D. However, at around t = 400-500 fs (500-600 fs after the arrival of the laser pulse) the mean kinetic energy of the nondesorbing adsorbates for a given D goes to the same value independently of the kind of isotope.In other words, the energy gain of the H atoms is reduced in the presence of D adsorbates due to frequent collisions among the adsorbates.The opposite happens for D atoms in the presence of H adsorbates.As speculated in Ref. [3], the more efficient energy uptake of the H adsorbate may promote a nearby surface reaction, D 2 desorption in this case.Now that the role of the adsorbates mobility on the desorption process is well established, the remaining point is to clarify how the efficiency of this heating process depends on the relative concentration of isotopes on the surface.In this respect, the analysis of the adsorbates that remain on the surface still in contact with the rest of heated adsorbates and metal electrons is more meaningful.Thus, the time dependence of E kin for the nondesorbing D (left panel) and H (right panel) is plotted in Fig. 5 for different D .The enhancement (reduction) of the mean kinetic energy of D (H) in the presence of H (D) at the surface is observed in the left (right) panel of Fig. 5.At around t = 300 fs (400 fs after the laser pulse hits the surface) the E kin curves of the adsorbed D atoms separate from each other so that the kinetic energy increases as D decreases (see left panel).The opposite behavior is observed for H in the right panel, although the effect is less apparent and somehow masked by the poor statistics.All in all, since the desorption rates are expected to increase the larger is the energy gained by the adsorbates, this explains why in our simulations the desorption of D 2 (H 2 ) is facilitated (hampered) in the presence of H (D) at the surface, and ultimately, the experimental observations regarding the dependence of the desorption yields on the relative isotope concentration.In other words, the energy gain of the H atoms is reduced in the presence of D adsorbates due to frequent collisions among the adsorbates.The opposite happens for D atoms in the presence of H adsorbates.It is worth mentioning that a similar interadsorbate energy transfer has recently been proposed in Ref. [41] in order to understand the (laser-induced) ultrafast desorption dynamics of CO adsorbed on Pd(111) at saturation coverage.

IV. CONCLUSIONS
In this work we have extended the AIMDEF method, already implemented in the VASP code, in order to incorporate a random force dependent on a electronic temperature that varies with time.This has allowed us to solve the Langevin equation for adsorbates at a surface by calculating on the fly all the forces (adiabatic, frictional, and stochastic).Making use of the two temperature model to obtain the time-dependent electronic temperature, we have applied the model to simulate the femtosecond laser induced desorption of H 2 , D 2 , and HD from the Ru(0001) surface at saturation coverage and different isotope concentration ratios.
The main result of our simulations is that upon laser excitation, and before the desorption of the first molecule takes place, the adsorbate system is collectively excited.As a result, formation and desorption of a molecule is not the result of the excitation of two neighbor atoms that takes place independent of the dynamics of the other adsorbates in the system.On the contrary, desorption of the molecule takes place in the presence of a heated adsorbate environment that influences the desorption dynamics.For this reason, the degree of heating of the complete adsorbate system governs the desorption rate and yield.Our simulations show that the energy gain of the adsorbates system coupled to the electronic system excited by the laser depends on the isotope concentration ratio.In mixed H:D coverage situations, though at the first stages H atoms gain energy more rapidly than D atoms, a state is reached at about 500-600 fs after the arrival of the laser pulse in which both kinds of isotopes have the same energy.This energy is higher (lower) the larger (smaller) is the relative number of H (D) atoms at the surface, i.e., the smaller (larger) is the average mass of the adsorbate system.This picture explains another important result of our simulations.We have found that the desorption yield of lighter (heavier) adsorbates is reduced (enhanced) in the present of the heavier (lighter) adsorbate.Note that this effect, dynamical promotion (suppression) of a surface reaction was observed experimentally in Ref. [3].However, up to now no theoretical modeling of the dynamics had been performed that could account for this effect.In our case, this has been possible due to the extension of the AIMDEF method that we present in this work that allows us to integrate the Langevin equation for varying electronic temperatures with all forces calculated on the fly.The importance of the method, in our case, is that it allows us to treat the dynamics of a large number of adsorbates, which is necessary to describe the "collective excitation" of the adsorbate system.Last but not least, we would like to stress that the methodology presented here is completely general and it can be used for different adsorbate/surface systems, different coverages, laser fluences, number of adsorbates in the simulation cell, and reactions.Some of these applications will require us to go beyond the frozen surface approximation, which can be readily done within the AIMDEF scheme.

Figure 2
Figure2shows the results of the AIMDEF simulations for the H 2 , D 2 , and HD desorption yields (Y H 2 , Y D 2 , and Y HD ) as a function of the D fraction of saturation coverage D .The yields plotted in the left panel are normalized such that for pure coverages ( D = 0 and 1) they correspond to the absolute yield per shot.In other words, the values are obtained by dividing the total number of desorbing H 2 , D 2 , and HD molecules by the total number of trajectories run for each D(40) and by the number of molecules of one kind available for pure coverages(8).In order to show more clearly the effect of the relative concentration of the adsorbates, the middle panel shows the H 2 and D 2 desorption yields normalized to the values obtained for the pure H and D adlayers, respectively.Let us start by discussing the isotope effect between isotopically pure coverages.For the value of the laser fluence that we consider in this work (F = 60 J/m 2 ) the MDEF simulations performed by Füchsel et al.[11] on a precalculated six-dimensional PES give a ratio between the H 2 and D 2 desorption yields of around 8 (see Fig.4in that reference).Since the isotope ratio is a very sensitive quantity, this result was considered to be in a quite reasonable agreement with

FIG. 2 .
FIG. 2. Left panel: Femtosecond laser-induced desorption yields of H 2 (green triangles and line), HD (blue squares and line), and D 2 (red circles and line) from a Ru(0001) surface against the D fraction of saturation coverage D .Middle panel: Normalized H 2 and D 2 desorption yields at the same conditions as in the left panel.In both panels, symbols denote the AIMDEF results, while solid lines denote the results from the second order rate equations (3).Right panel: Nearest neighbors (NN) distribution of the adsorbates forming the desorbing molecules.Different colored symbols correspond to different D .

FIG. 4 .
FIG. 4. Mean kinetic energy of the desorbing D and H (thick black and red lines, respectively) and nondesorbing D and H (thin black and red lines, respectively) as a function of time.In each case, the average is calculated over the corresponding total number of atoms that experience a specific event.Each panel corresponds to different D fraction of saturation coverage: Top-left, 0 and 1, top-right, 0.75, bottom-left, 0.5, and bottom-right, 0.25.