Photo-Induced CO Desorption Dominates over Oxidation on Different O+CO Covered Ru(0001) Surfaces

The photo-induced desorption and oxidation of CO on Ru(0001) is simulated using ab initio molecular dynamics with electronic friction that accounts for the non-equilibrated excited electrons and phonons. Different (O,CO) coverages are considered, the experimental room temperature coverage consisting in 0.5ML-O+0.25ML- CO (low coverage), the saturation coverage achieved experimentally at low temperatures (0.5ML-O+0.375ML-CO, intermediate coverage), and the equally mixed monolayer that is stable according to our calculations but not experimentally observed yet (0.5ML-O+0.5ML-CO, high coverage). The results of our simulations for the three coverages are consistent with femtosecond laser experiments showing that the CO photo-desorption largely dominates over CO photo-oxidation. These results cannot be explained in terms of the distinct activation energies calculated for the relaxed surfaces. Different (dynamical) factors such as the coupling to the laser-excited electrons and, more importantly, the interadsorbate energy exchange and the strong surface distortions induced in the more crowded surfaces are fundamental to understand the competition between these two processes under the extremely non-equilibrated conditions created by the laser.


Introduction
Understanding the oxidation of CO on transition metal surfaces is one of the long-standing important problems in surface science that continues motivating extensive fundamental research. [1][2][3][4][5][6][7] A particularly interesting surface is Ru(0001) with coadsorbed CO and atomic O.
It has been demonstrated that under ultra high vacuum (UHV) conditions the oxidation process cannot be thermally activated. 8,9 Interestingly, the reaction can be propelled by exciting the system with electromagnetic radiation. [9][10][11][12] Specifically, this has been achieved by using near-infrared femtosecond laser pulses that excite the electrons of the system. Subsequently, those laser-excited electrons transfer energy directly to the adsorbates, but also indirectly, via excited phonons that result from electron-phonon coupling. Under these conditions, both CO and CO 2 are desorbed from the surface.
A proper characterization of this kind of experiments requires an accurate simulation of the reaction dynamics. In a successful approach to this problem, the response of the metal to the laser pulse is described in terms of coupled electronic and phononic excitations characterized, respectively, by time-dependent electronic (T e ) and phononic (T l ) temperatures. 13 Subsequently, simulations of the adsorbates dynamics in this highly excited system are performed, accurately describing the adsorbate-surface interaction with potential energy surfaces based on density functional theory (DFT). [14][15][16][17][18] In this respect, the development of the ab initio molecular dynamics with electronic friction method (AIMDEF), 19-21 and in particular its version that incorporates electronic and phononic temperatures ((T e , T l )-AIMDEF), has allowed to treat the multidimensional nonadiabatic dynamics of this kind of photo-induced reactions at the DFT level taking into account the coupling of the adsorbates to only the excited electrons in some cases 22 and to both the excited electronic and phononic systems in others. [23][24][25] Also the coadsorbate energy exchange that may become important as surface coverage increases 26,27 can be readily included in this kind of simulations. 22,23 The (T e , T l )-AIMDEF scheme is of special relevance for CO oxidation at photo-excited metal surfaces due to the high dimensionality involved in the reaction. In addition to the minimal 9 dimensions of the reacting O and CO, the mentioned coadsorbate energy exchange, energy exchange of the adsorbates (particularly CO) with the surface atoms, and surface distortions caused not only by the excited phonons, but also by the excited adsorbates, including possible formation of subsurface O absorbates at high surface temperatures, 28,29 are important factors that must be included in a realistic simulation.
The preparation of the adsorbate layer in the experiments of refs. 9-11 is performed in the following way. First, the Ru(0001) surface is dosed with oxygen up to saturation which consists in a 0.5 ML O coverage. 8 In a second step, the dosing of the surface with CO is produced also up to saturation. Theoretical studies motivated by these experiments, devoted to the description of the energetics of the reaction path 11 and to the use of kinetic models 10 to understand the measured desorption yields, assumed a 0.25 ML CO coverage.
Also, in a recent work, the (T e , T l )-AIMDEF method has been used for the first time to simulate the femtosecond laser pulse induced oxidation of CO in the Ru(0001) surface for this coverage. 25 In that work, the simulations reproduced and allowed to understand the large experimental branching ratio between desorption and oxidation that exceeds one order of magnitude, in spite of the fact the CO oxidation is energetically favored over CO desorption at that 0.5 ML O + 0.25 ML CO surface coverage. Nevertheless, though this is the reported saturation coverage at room temperature, 8 in the experiments the temperature was lowered before the laser excitation to 100 K, [9][10][11] for which a higher CO saturation coverage cannot be excluded. 8 The above discussed facts constitute our motivation to perform (T e , T l )-AIMDEF for different surface coverages. More precisely, we have simulated the photo-induced CO desorption and oxidation on three different O+CO covered Ru(0001) surfaces, for which the equilibrium structures and minimum energy reaction paths for CO desorption and oxidation were studied and characterized using DFT in a previous work. 30 In the three cases, the O coverage is fixed to 0.5 ML and only the CO coverage, which is not clearly determined in the femtosecond laser experiments, is different. We will denote these CO coverages as low, intermediate, and high coverages following the nomenclature of ref. 30 30 and hence useful to understand the coverage dependence of the photo-induced desorption and oxidation of CO.
The paper is organized as follows. The theoretical methods and computational settings used to model the different (O,CO)-covered Ru(0001) surfaces are described in the Methods section. The results of the photo-induced oxidation and desorption dynamics, obtained from our (T e , T l )-AIMDEF simulations for each of the three coverages considered, are presented and comparatively discussed in the Results and Discussion section. The main conclusions extracted from our dynamics simulations are summarized in the Conclusions section.

Photo-Induced Desorption Model
The photo-induced desorption and oxidation of CO from the (O,CO)-covered Ru(0001) surface is simulated with ab initio classical molecular dynamics using the (T e , T l )-AIMDEF methodology. 23 These kinds of simulations incorporate the effect of the laser-induced excited electrons in both the adsorbates dynamics and the surface atoms dynamics. An extensive description of this method and its applicability to model femtosecond laser-induced reactions on metals can be found in refs. 23,25. Hence, only its main ingredients will be reviewed next.
As a first step, the response of the metal surface to the near-infrared laser pulse is described within the two-temperature model (2TM). 13 In this model, the laser-induced electronic and ensuing (electron-induced) phononic excitations are represented by two coupled heat thermal baths, being T e (t) and T l (t) their associated time-dependent temperatures obtained from the following equations: where C e is the the electron heat capacity, C l is the phonon heat capacity, κ is the electron thermal conductivity, g is the electron-phonon coupling constant, z the perpendicular position relative to the surface and S(z, t) is the absorbed laser power per unit volume that depends on the shape, wavelength, and fluence of the applied pulse.
Once T e (t) and T l (t) are known, the effect of the laser-excited electrons on each adsorbate is described through the following Langevin equation: where m i , r i , and η e,i are the mass, position vector, and electronic friction coefficient of the i th atom conforming the set of adsorbates. The first term in the right hand side of the equation Furthermore, the heating of the surface lattice due to the laser-induced electronic excitations, which is described through T l (t), is assured by coupling the surface atoms in the two upper layers to the Nosé-Hoover thermostat. 34,35 Hence, these surface atoms with mass m j and position vector r j are subjected to the following equations of motion, where N is the number of atoms coupled to the thermostat (in our case the first two layers), Q is a parameter with dimensions of energy×time 2 that acts as the mass of the dynamical variable s with associated conjugated momentum p s , and ξ = Q −1 sp s is the thermodynamic friction coefficient. 36 Finally, the fourth and fifth Ru layers are kept frozen, while the third layer is used as a transition layer to simulate that the lattice heating flows from the surface to the bulk. Hence, the movement of the third layer atoms is described by the classical Newton equations of motion and the adiabatic approximation. The (T e , T l )-AIMDEF method and similar methods that combine the 2TM and the Langevin dynamics for the adsorbates have been widely used to simulate the femtosecond laser induced dynamics and reactions of adsorbates at metal surfaces. [14][15][16][17][18][22][23][24][37][38][39] General DFT Computational Settings  During the (T e , T l )-AIMDEF simulations, the adiabatic forces are calculated with non spin-polarized DFT using the van der Waals exchange-correlation functional proposed by Dion et al. 48 and the same computational parameters that were used in our previous study on the energetics of CO desorption and oxidation at different coverages. 30 Specifically, the electronic ground state is determined at each integration step by minimizing the system total energy up to a precision of 10 −6 eV. In doing so, integration in the Brillouin zone is performed using a Γ-centered 3×6×1 Monkhorst-Pack grid of special k points 49 and the Methfessel and Paxton scheme of first order with a broadening of 0.1 eV to describe partial occupancies of each state. 50 The latter are expanded in a plane-wave basis set with an energy cut-off of 400 eV, whereas the electron-core interaction is treated with the projected augmented wave (PAW) method 51 that is implemented in VASP. 52 All dynamics simulations are performed with the Beeman integrator implemented in our AIMDEF module 42 using a time step of 1 fs.
Each trajectory is propagated up to 4 ps. As possible outcomes, we observe CO desorption and CO oxidation (i.e., the O+CO recombinative desorption). Atomic O is strongly adsorbed and no desorption events are obtained. Oxygen subsurface migration is however temporarily possible, being more probable the higher the coverage.

Initial Conditions
For each coverage, a preliminary simulation is performed in which the three outmost Ru layers are equilibrated at 100 K during 20 ps using the Nosé-Hoover thermostat 34,35 [Eqs. (3) and (4)

Calculation of observables
A molecule (CO or CO 2 ) is counted as desorbed, if its center of mass crosses the plane z=6.5Å, measured respect to the topmost Ru surface layer, with positive velocity along the surface normal. Since the Ru atoms are also moving, the reference z plane is calculated at each instant t as the mean height of the Ru atoms in the topmost layer. The CO desorption and oxidation probabilities per CO molecule are calculated for each coverage as where N des (A) is the number of desorbing molecules (CO or CO 2 ), N t is the total number of trajectories, and N CO is the number of CO molecules in the simulation cell (2, 3, and 4 for low, intermediate, and high coverages, respectively). In the intermediate coverage, the site-resolved desorption probabilities refers to the desorption probabilities for each of the three distinct CO molecules in the cell, i.e., P des (CO1), P des (CO2), and P des (CO3).
The mean total kinetic energy E kin (t) and mean center-of-mass kinetic energy E cm (t) per adsorbate type (i.e., O or CO) are calculated at each instant t as an average over all the trajectories and all the atomic or molecular species under consideration.

Results and Discussion
All the simulations correspond to irradiating the surface with a 800 nm Gaussian pulse of 110 fs duration and absorbed fluence F=200 J/m 2 as in experiments. 9 The corresponding T e (t) and T l (t) curves calculated with Eq. (1) (input parameters for the Ru slab as in refs. 14, 18,22,25) are shown in Figure 2. The CO desorption and CO oxidation probabilities obtained from the (T e , T l )-AIMDEF simulations for each coverage are summarized in Table 1.
Starting with the CO desorption process, the largest to the lowest probabilities correspond to the intermediate, high, and low coverages. As also shown in the table, this ordering agrees with the values of the CO desorption barriers obtained in our previous DFT analysis of the coverage-dependent energetics of the two reactions. 30 The photo-oxidation process also occurs on the three coverages. The highest probability for this reaction is unexpectedly obtained for the high coverage, i.e., the coverage for which the activation barrier for CO oxidation is the largest. The analysis of the oxidizing trajectories will allow us to rationalize this result. Compared to the CO desorption process, the CO oxidation probabilities are considerably lower in the three cases. In spite of having a limited statistics for CO oxidation, let us remark that the resulting large branching ratios between desorption and oxidation are in line with the experimental observations, P des (CO)/P des (CO 2 )= 35 9 and 31. 10 As shown in ref. 25 for the low coverage, the extremely small P des (CO 2 ) values are due to the very reduced configurational space leading to CO oxidation as compared to that for CO desorption. More precisely, at low coverage, we found that oxidation is unlikely due to the difficult access to the transition state region that requires the O atom crossing the bridge site and finding the CO conveniently close and tilted to form the CO 2 molecule and due also to the fact that this access does not guarantee a successful recombination. In this work, we will hence focus in understanding how the CO desorption and CO oxidation dynamics depend on the Ru (0001) coverage. Table 1: (T e , T l )-AIMDEF CO desorption probability P des (CO), CO oxidation probability P des (CO 2 ), and CO to CO 2 branching ratio calculated for an absorbed fluence F=200 J/m 2 and three different (O,CO)-Ru(0001) coverages. N t is the number of trajectories run for each coverage. Activation energies for CO desorption E TS (CO) and CO oxidation E TS (CO 2 ) are from ref 30.

CO Desorption Dynamics
The time evolution of the CO desorption probabilities P des (t) for the three coverages are compared in Figure 2  for the high coverage. During this initial stage, the main mechanism contributing to the adsorbates excitation is the coupling to the hot electrons (compare T e to T l in Figure 2).
Since CO adsorbs on hollow sites in the high coverage case (i.e., very close to the surface) and at near-top and top sites in the intermediate and low coverage cases, respectively, this coupling is stronger in the former than in the two latter, explaining the slightly larger initial energy gain in the high coverage. This larger energy gain may also explain that P des (CO) is initially larger for this coverage and not for the intermediate one, which has a lower activation energy for desorption (see Table 1). There are various interesting features when comparing the mean kinetic energy among the three coverages. In the case of the desorbing CO, E kin is the largest during the whole integration interval for the high coverage.     coverages and, in the intermediate coverage, we verified that the desorption probability is unaffected by the slight differences between the three initial adsorption sites (CO1, CO2, and CO3 in Figure 1), since the site-resolved desorption probabilities are basically the same for the three distinct CO. We do observe differences when comparing the mean kinetic energies E kin (t). Figure 5 shows that at a certain instant (different for each coverage) the desorbing CO molecules (dashed lines) start to have more kinetic energy than the nondesorbing (solid lines) molecules (compare for each coverage the black dashed to the black solid lines). Note also that the extra energy gained by the desorbing CO mainly goes to the center-of-mass degrees of freedom. This can be inferred from the figure because the difference between the mean total kinetic energy of the desorbing and non-desorbing CO (i.e., difference between black dashed and black solid lines) is quite similar to the difference between the mean center-of-mass kinetic energy of the desorbing and non-desorbing CO (i.e., difference between red dashed and red solid lines). By decomposing the mean center-of-mass kinetic energy into parallel ( E par cm (t), in blue) and perpendicular to the surface ( E z cm (t), in green) contributions, it is clear that the desorbing molecules are characterized by having the largest E z cm (t) (green dashed lines), as one would naively expect. Interestingly, the parallel contribution of the desorbing and non-desorbing molecules is basically identical for each coverage. Altogether, we conclude that the desorbing molecules are characterized by a larger net energy gain that basically goes into the center-of-mass perpendicular motion. Note in passing that at the end of the simulations, the relationships between the kinetic energies for the non-desorbing CO verify roughly equipartition among the different degrees of freedom, i.e., E kin : E par cm : E z cm ≡ 6:2:1, that suggest a system arriving to a quasithermalized state at the end of the simulations.

CO Oxidation Dynamics
The number of oxidation events in our simulations is very scarce which, in principle, represents a question mark about the statistical validity of the results. However, taking this limitation in mind, the following observations on the differences found on the oxidation dynamics at the different coverages look very general and robust. The instants at which CO 2 desorbs from each covered surface are compared in Figure 2 (top panel). The two oxidation events obtained in the low coverage occur at well separated instants, t ∼2 and 4 ps, respectively. In contrast, in the high and intermediate coverages, all the CO 2 desorption events occur within a shorter time interval of about 1 ps. The important difference between these two coverages is that the first event takes place at t ∼1.4 ps in the former and at t ∼2.6 ps in the latter. Considering that the high coverage has the largest activation energy for CO oxidation, not only is it surprising to find that the oxidation probability is the highest for this coverage (see Table 1) but also that the process initiates more rapidly on this surface.
There are various reasons that can explain this unexpected result. As above discussed, the adsorbates gain more energy and at faster rate in the high coverage because the coupling to the excited electrons is initially the strongest and, importantly, because the interadsorbates energy exchange is more efficient in this crowded surface. Furthermore, we cannot discard that the energy barriers on the highly excited surface are different from those calculated on the relaxed surface. In this respect, the analysis of the CO oxidation trajectories below allows us to determine up to what extent the followed oxidation paths coincide with the minimum reaction paths we identified for each coverage in our previous static calculations. 30 As already discussed in ref 25, the two CO oxidation trajectories in the low coverage follow the corresponding DFT+vdW-DF minimum energy path (MEP), 30 in which the less bound O fcc surmounts its adsorption well, crosses the bridge site between two Ru atoms, and recombines with the nearby CO that tilts to form the chemisorbed bent CO 2 (bCO 2 ).
However, compared to the low coverage, not only the adsorbates become more rapidly very mobile in the intermediate and high coverages (see Figure 4)  We observe that the configurational space leading to CO desorption is quite ample. As a consequence, the CO desorption is a rather simple and direct process only limited by the energy the CO molecules need to gain in order to overcome the energy barrier to desorption. This is reflected in the dependence of the CO desorption probability on coverage. This probability is the largest for the intermediate coverage, which has the lowest energy barrier to desorption (around 0.6 eV), and it is the smallest for the low coverage, which presents the highest energy barrier to desorption (1.57 eV). By comparing the time evolution of the desorption probabilities, we identify other factors that also contribute to the dependence of the desorption dynamics on coverage, the coupling to the electronic system and the interadsorbate energy exchange. The former, being quite similar for the three coverages, is only responsible of the slightly larger energy that the CO molecules gain initially in the high coverage. The efficient interadsorbate energy exchange that is taking place in the high and intermediate coverages explains the faster energy uptake and faster desorption dynamics observed in these coverages as compared to the low coverage. In the three cases, the desorbing CO are characterized by gaining more energy than the non-desorbing CO. The extra energy is mainly deposited in the translational degree of freedom normal to the surface.
In contrasts to CO desorption, the CO oxidation probabilities do not correlate with the values of the activation energies for CO oxidation that were calculated in the relaxed surface. For instance, the CO oxidation probability is the highest for the high coverage, i.e., the coverage with the largest activation energy. Our dynamics simulations show that CO oxidation in these nonequilibrium conditions presents a much more complex scenario than CO desorption. Except for the low coverage case, the oxidation events taking place in our dynamics do not follow the minimum energy paths calculated under equilibrium conditions.
Moreover, oxidation is sometimes preceded in the intermediate and high coverages by CO desorption events that modify the local environment. The dynamics simulations also show that the excited adsorbates, that become particularly mobile in the intermediate and high coverages, can cause profound surface distortions. In particular, a common event occurring in the high coverage that is favoring the formation of the chemisorbed bent CO 2 , and subsequent CO 2 desorption, is the pronounce outwards displacement of the Ru atom that binds the reactive CO. As a consequence of all it, the relevant energy barriers for the oxidation process in the highly excited surface are different from those calculated on the relaxed surface. All in all, these results show the importance of accounting fot the adsorbate and surface dynamics under extremely non-equilibrium conditions in order to understand the laser induced photochemistry at metal surfaces.