Understanding the Photoinduced Desorption and Oxidation of CO on Ru(0001) Using a Neural Network Potential Energy Surface

The study of ultrafast photoinduced dynamics of adsorbates on metal surfaces requires thorough investigation of laser-excited electrons and, in many cases, the highly excited surface lattice. While ab initio molecular dynamics with electronic friction and thermostats (Te, Tl)-AIMDEF addresses such complex modeling, it imposes severe computational costs, hindering quantitative comparison with experimental desorption probabilities. In order to bypass this limitation, we utilize the embedded atom neural network method to construct a potential energy surface (PES) for the coadsorption of CO and O on Ru(0001). Our results demonstrate that this PES not only reproduces the short-time ab initio dynamics but is also able to yield statistically significant data for long lasting trajectories that correlate well with experimental findings. Furthermore, the analysis of the laser-induced dynamics reveals the existence of a dynamic trapping state that acts as a precursor for CO desorption, and it is not observed under thermal conditions. Altogether, our results validate the underlying theoretical framework, providing robust support for the description of not only the photoinduced desorption but also the oxidation of CO in terms of nonequilibrated but thermal hot electrons and phonons.

The (T e ,T l )-AIMDEF simulations performed by Tetenoire et al.S1 were carried out with vasp S2,S3 (version 5.4) and the AIMDEF module S4-S8 using a periodic supercell that consisted of a vector along the surface normal of 30.2225Å that included about 19 Å of vacuum and the same (4×2) surface unit cell that we use here in (T e ,T l )-MDEF.Within this supercell, the (0.25ML CO+0.5MLO)/Ru(0001) surface was described by five layers of Ru atoms and one (CO+2O) layer adsorbed on the topmost Ru surface layer.In the adlayer, each CO adsorbs atop a Ru atom, while the O atoms occupy the second nearest hcp and fcc sites forming a honeycomb arrangement around the CO.
During the AIMDEF simulations, the adiabatic forces were calculated with non spinpolarized DFT using the van der Waals exchange-correlation functional proposed by Dion et al.S9 The equations of motion were integrated using a time step of 1 fs and the Beeman integrator implemented in both the AIMDEF module S4 and the (T e , T l )-MDEF dynamics code.S10 On each integration step, the electronic ground state was determined by minimizing the system total energy up to a precision of 10 −6 eV.Integration in the Brillouin zone was performed using a Γ-centered 3×6×1 Monkhorst-Pack grid of special k points S11 and the Methfessel and Paxton scheme of first order with a broadening of 0.1 eV for occupations.S12 An energy cut-off of 400 eV was used in the plane-wave expansion of the valence electrons, whereas the core electrons were treated with the projected augmented wave (PAW) method S13 using the VASP implementation.S14 Using all the above parameters, 200 (T e , T l )-AIMDEF trajectories with total simulation time of 4 ps each were run for the CO/2O/Ru(0001) surface.The latter was initially thermalized at 100 K and next irradiated with an 800 nm Gaussian pulse of 110 fs duration and an absorbed fluence F = 200 J/m 2 .The T e (t) and T l (t) curves were calculated with the 2TM using the input parameters summarized in Table S2 of

S6
We have calculated using our EANN-PES the CO desorption energy and activation energy for CO oxidation in the equilibrium surface.These values compare rather well with the corresponding DFT data considering that no information on desorption or oxidation under equilibrium conditions was used for training the NNPES.Table S1: CO desorption energy (E des (CO)) and activation energy for CO oxidation (E a (CO 2 )) in the equlibrium surface (in eV) calculated with the EANN PES are compared to their corresponding values calculated with DFT+vdW-DF in S15.

EANN PES
where C e is the the electronic heat capacity that is assumed to depend linearly on T e (i.e., C e = γ e T e ), C l is the phonon heat capacity described within the Debye model, S16 κ e is the electronic thermal conductivity described as κ e (T e , T l ) = κ 0 T e /T l , S16 g is the electron-phonon coupling parameter, z the perpendicular position relative to the surface, and S(z, t) is the absorbed laser power per unit volume that for a Gaussian pulse of wavelength λ and absorbed fluence F is given by In this equation, τ = FWHM/(2 √ 2 ln 2), t 0 is the instant at which the pulse intensity is maximum, and ξ is the light penetration depth of the metal surface for wavelength λ, which is calculated from the imaginary part of the reflective index k λ as ξ = λ/(4πk λ ).Recall that the above equations are implicitly assuming a constant absorbed laser fluence in the plane parallel to the metal surface and neglecting lattice thermal diffusion into the bulk.
The diameter of the experimental laser beams and the time scale of interest justify each of those assumptions.S17,S18 The corresponding T e (t) and T l (t) curves are shown in Figure S7 for a representative selection of absorbed laser fluences.

S4 Surface electron density
Within the Local Density Friction Approximation (LDFA) S21,S22 used in our simulations, the electronic friction coefficient applied at each adsorbate depends on the value of the electronic density of the bare Ru(0001) surface at the position of the adsorbate, ρ 0 (R i ).In the (T e , T l )-AIMDEF approach, ρ 0 (R i ) is calculated using the Hirshfeld partitioning scheme S23 to remove the adsorbate density contribution from the (DFT-calculated) self-consistent electron density of the complete system (adsorbates and Ru surface) that is calculated at each integration step (see refs S5,S6,S22 for details).Conversely, in (T e , T l )-MDEF the electron density of the bare Ru(0001) surface at the position of each adsorbate is modeled as a sum of the electronic densities contributed by individual Ru atoms at this position.The contribution of each Ru atom to the density is described by two decaying exponential functions.As a result, the expression for ρ 0 (R i ) is the following:

S5 Statistics of the (T e ,T l )-MDEF simulations
The specific number of trajectories run for each absorbed laser fluence as well as the total simulation time used in each case are given in Table S3.It is worth to remark that a single (T e , T l )-AIMDEF trajectory with a simulation time of 4 ps and a time step of 1 fs (i.e., 4000 integration steps) took between 5-7 days on 24 cores (Intel(R) Xeon(R) Platinum 8168 CPU @ 2.70GHz), whereas a single (T e , T l )-MDEF trajectory with a simulation time of 30 S10 ps and a time step of 0.5 fs (i.e., 60 000 integration steps) takes between 2 to 3 minutes on 36 cores (mixed cluster with Intel(R) Xeon(R) Platinum and Intel(R) Xeon(R) Gold, CPUs @ 2.70GHz-3.00GHz).
Table S3: Details of the (T e , T l )-MDEF simulations carried out for each absorbed laser fluence: number of trajectories N traj and simulation time t end .

S6 Time-convergence of the Desorption Probability
As defined in the Results and Discussion section of the main text, in the (T e ,T l )-MDEF simulations presented there, a molecule is considered to be desorbed if its center of mass height from the topmost Ru surface layer is ≥ 10.5 Å and its center of mass velocity along the surface normal is positive.Following this definition, the time evolution of the desorption probability is calculated at each integration step t as where N traj and N CO are the total number of trajectories for a given fluence and total number of CO molecules present in the unit cell, respectively.N des,i (t) is the number of desorbed molecules at the time-step t for a given trajectory i.
where the variance is obtained as

S7 Dynamic trapping
An interesting finding visible in Figure 3 of the main text is the dynamic trapping of CO molecules.In order to understand this behaviour, analysis was conducted by displacing the center of mass of one CO molecule in the x − y plane over a range of CO center-of-mass heights from the ruthenium surface, Z cm = 3 − 6.5 Å.For each CO height and position over the surface plane, the polar angle of the molecular axis was varied within the range [0 -π]. Figure S9: 2D cuts of the potential energy surface calculated by the NNPES for the thermalized surface at 100 K (left), and for a distorted surface taken from one of the MDEF trajectories exhibiting CO dynamics trapping (right).The cuts were calculated at distances of the CO center of mass from the ruthenium surface, Z cm = 3.1 Å (top panels), Z cm = 4.5 Å (middle panels), and Z cm = 6 Å (bottom panels), respectively.In all cases, the polar angle of the molecular axis for displaced CO molecule is π/4.Top views of the surface unit cell for the thermalized and distorted surfaces are shown atop the left and right panels, respectively.In both cases, the CO that is being displaced to create the 2D plots is the rightmost one.
Figures S1 and S2 show the distribution of the potential energies and atomic force components (F x , F y , F z ) that form the final training data set.As explained in the main manuscript, the data correspond to configurations that were probed during the 200 (T e , T l )-AIMDEF trajectories with total simulation time of 4 ps each that were calculated in ref S1 for the CO/2O/Ru(0001) surface and an absorbed laser fluence of 200 J/m 2 .The system can probe an ample configurational space that is characterized by a wide potential energy range of about 15.6 eV.

Figure S1 :
Figure S1: Distribution of the system potential energy probed during the 200 (T e , T l )-AIMDEF trajectories calculated in ref S1 for the CO/2O/Ru(0001) surface and F = 200 J/m 2 .

Figure S2 :
Figure S2: Same as Figure S1 but for each atomic force component.

Figures S3 ,
Figures S3, S4, S5, and S6 show the comparison of the EANN-PES predicted energies and the three cartesian components of the atomic forces with the corresponding DFT values for the configurations of the test data set.

Figure S3 :
Figure S3: Comparison of the potential energies computed with the EANN-PES for the configurations in test data set and the corresponding DFT values.

Figure S4 :
Figure S4: Comparison of the z component of the atomic forces computed with the EANN-PES for the configurations in test data set and the corresponding DFT values.

Figure S5 :
Figure S5: Comparison of the x component of the atomic forces computed with the EANN-PES for the configurations in test data set and the corresponding DFT values.

Figure S6 :
Figure S6: Comparison of the y component of the atomic forces computed with the EANN-PES for the configurations in test data set and the corresponding DFT values.

Figure S7 :
Figure S7: Electronic (solid lines) and lattice (dashed lines) temperatures calculated with the two temperature model for a Ru surface irradiated by an 800 nm Gaussian laser pulse, FWHM 110 fs, and various absorbed fluences F .The maximum of the laser pulse intensity is at t = 236 fs.
S3)   where r Ru is the position of the Ru atom.After fitting this function to the electronic density data extracted from the (T e , T l )-AIMDEF simulations, the values of the parameters are: a = 4.178222 a.u., b = 4.022414 Å−1 , c = 0.061987 a.u., and d = 1.876455Å−1 (a.u.stands for atomic units).

Figure S8 shows
FigureS8shows the time evolution of the CO and CO 2 desorption probabilities for the representative absorbed fluences, 125 J/m 2 , 160 J/m 2 , and 200 J/m 2 .The values of the CO and CO 2 desorption probabilities that are compared to experimental data correspond to the values obtained at the end of simulation time for each fluence, P des (t = t end ).As shown in this figure, all the P des (t) curves are well saturated at the end of the simulation time, demonstrating that the desorption probabilities reported in the Results and Discussion

Figure S8 :
Figure S8: Time dependent CO desorption probability (blue) and CO 2 desorption probability (orange) for three representative absorbed laser fluences that cover the whole F -range used in our calculations: 125 J/m 2 (left), 160 J/m 2 (middle), and 200 J/m 2 (right).

Figure S9
Figure S9 illustrates an example of an energy landscape calculated by the NNPES for the thermalized, slightly perturbed surface at 100 K, and for a distorted surface taken from one

Table S2 :
Constants and parameters used in the 2TM model for the Ru surface, taken from refs S19,S20.