Photoinduced Desorption Dynamics of CO from Pd(111): A Neural Network Approach

Modeling the ultrafast photoinduced dynamics and reactivity of adsorbates on metals requires including the effect of the laser-excited electrons and, in many cases, also the effect of the highly excited surface lattice. Although the recent ab initio molecular dynamics with electronic friction and thermostats, (Te,Tl)-AIMDEF [AlducinM.;Phys. Rev. Lett.2019, 123, 246802]31922860, enables such complex modeling, its computational cost may limit its applicability. Here, we use the new embedded atom neural network (EANN) method [ZhangY.;J. Phys. Chem. Lett.2019, 10, 496231397157] to develop an accurate and extremely complex potential energy surface (PES) that allows us a detailed and reliable description of the photoinduced desorption of CO from the Pd(111) surface with a coverage of 0.75 monolayer. Molecular dynamics simulations performed on this EANN-PES reproduce the (Te,Tl)-AIMDEF results with a remarkable level of accuracy. This demonstrates the outstanding performance of the obtained EANN-PES that is able to reproduce available density functional theory (DFT) data for an extensive range of surface temperatures (90–1000 K); a large number of degrees of freedom, those corresponding to six CO adsorbates and 24 moving surface atoms; and the varying CO coverage caused by the abundant desorption events.


INTRODUCTION
The use of intense (∼1 mJ/cm 2 ) femtosecond (fs) laser pulses in the ultraviolet, visible, and near-infrared regime has been shown to be a very efficient way to promote reactions at adsorbate-covered metal surfaces. 1−3 At these wavelengths, a large fraction of the light is absorbed by the metal giving rise to electronic excitations. Subsequently, energy transfer to the lattice atoms takes place via electron−phonon coupling. As a result, the adsorbates encounter a combined electronic and phononic excited system from which they can gain energy and experience different reactions, diffusion, and even desorption from the surface. 4−9 Interestingly, this kind of excitation mechanism can increase significantly the cross section of reactions with respect to what is observed under ordinary thermal excitation conditions and even open new reaction channels. Two pulse correlation experiments are customarily employed to obtain the time scale of the energy transfer between the adsorbate and the substrate. In this way, in principle, whether a specific reaction is mainly governed by the excited electrons or phonons can be disentangled experimentally. However, in several cases, this information is not unequivocally obtained from the experiments and theoretical modeling is necessary.
From a theoretical aspect, the modeling of these experiments requires performing molecular dynamics simulations in an excited environment. 10−17 First, the excitation generated by the laser pulse in the substrate is described in terms of timedependent electronic (T e ) and phononic (T l ) temperatures that are obtained using the two-temperature model (2TM). 18 Subsequently, the motion of the adsorbates is determined by solving Langevin equations of motion in the ground-state potential energy surface (PES). In this way, the coupling of the adsorbates to the electronic system is modeled in terms of electronic friction forces and associated stochastic forces that depend on T e . A nonempirical and accurate potential energy surface is typically obtained by characterizing the adsorbate− metal surface interaction at the level of density functional theory (DFT). Until very recently, due to the large computational cost involved in the DFT calculations, potential energy surfaces of reduced dimensionality, only involving the degrees of freedom (DOFs) of the adsorbate, were used. As a result, the effect of the heated phonon system was either not included at all in the dynamics or in a rather approximate way using the generalized Langevin oscillator (GLO) model that does not account for independent surface atom move-ment. 13−17 Moreover, a maximum of two atomic adsorbates (or a single diatomic molecule) were described with sixdimensional potential energy surfaces, which did not allow one to model interadsorbate energy exchange and study effects related to the local reduction of the coverage that is caused by sequential desorption events.
In this respect, only very recently have these limitations been overcome using ab initio molecular dynamics with electronic friction (AIMDEF). 19−27 The model, hereafter denoted as (T e ,T l )-AIMDEF, 27 incorporates both the electronic and phononic excitation channels: the former by solving the Langevin equation for the adsorbates using T e and the latter by coupling the surface atoms to a thermostat at a temperature T l . This methodological approach naturally includes all of the system's degrees of freedom as required. Specifically, in principle, any number of surface atoms are allowed to move independently and multiple adsorbates can be treated. Using (T e ,T l )-AIMDEF, the importance of including not only T e but also T l (particularly for those surfaces that may reach a high T l ), as well as the interadsorbate interactions, has been demonstrated. 27 The main shortcoming of the approach is that AIMDEF is extremely computationally demanding. This means that, in practice, a reduced number of the order of few hundreds of trajectories can be realistically computed for a given set of experimental conditions. This results in limited statistics. For the same reason, with the required time steps of the order of femtoseconds, the integration time is limited to around 2−4 ps.
In the last few years, the use of neural network (NN)generated multidimensional PESs has become an accurate alternative to ab initio molecular dynamics (AIMD) to describe the dynamics of diverse gas−surface processes 28−39 and also the dynamics at solid−liquid water interfaces. 40−43 In particular, for these studies, the development of the atomistic neural network (AtNN) approach has constituted a major advancement. 44−46 Within AtNN, the PES is constructed in terms of atomistic contributions, which allows for obtaining NN-PESs that are a function of all the atomic positions in systems of arbitrary size. Application of this methodology to gas−surface dynamics studies has allowed for constructing PESs not only for diatomic 31−33,35,39 but also for polyatomic molecules 29,30,34,36,37 interacting with surfaces. Moreover, since the NN-PES can also be made a function of the surface atom coordinates, the treatment of both the independent surface atom movement and the surface temperature effects has also been performed, which has allowed for accounting for energy exchange between the molecule and the surface along the dynamics. 31 −33,35−39 However, it must be emphasized that the requirements imposed on a NN-PES capable of describing femtosecond laser-induced reactions are extremely demanding as compared to those required in usual elementary gas−surface processes. Since it is necessary to model the movement of multiple adsorbates and surface atoms, the number of degrees of freedom to be accounted for is huge. In this respect, it is worth mentioning that, to our knowledge, all the NN-PESs generated up to now for gas−surface studies are restricted to a single gas molecule. Moreover, the PES must be able to describe accurately the very distinct and changing adsorbate coverages that exist during the photoinduced dynamics. This means that it is necessary to ensure a precise description of adsorbate− substrate and interadsorbate interactions under very different conditions involving local changes of the configuration space of neighbor adsorbates and strong lattice distortions. In this respect, note that in these experiments, the lattice temperatures (T l ) may vary rapidly in the range of 90−1000 K, which implies that the configuration space corresponding to the surface atoms is very large. Therefore, altogether the requirements for the atomistic NN-PES are unprecedentedly extreme and demanding for these kinds of processes. Here, we show that the recently developed embedded atom neural network (EANN) method, 47 which has already been successfully applied to construct the PESs for diatomic and polyatomic molecules interacting with multiple metal facets, 48,49 is indeed impressively accurate and flexible to account for all these necessities.
For this purpose, we study the femtosecond laser-induced desorption of CO from Pd(111). 9 In particular, we concentrate on the CO saturation coverage of 0.75 monolayer (ML), in which the CO molecules adsorb in atop, face-centered cubic (fcc), and hexagonal close-packed (hcp) sites. (T e ,T l )-AIMDEF results for this system were recently presented in ref 27. In the present work, we use the configurations encountered along these dynamics as the input data to generate our EANN-PES.
The work is organized as follows. In Section 2, the procedure used to construct EANN-PES is detailed and accuracy tests are presented. Next, Section 3 is devoted to the description of the theoretical framework used to perform molecular dynamics simulations of the laser-induced reactions at metal surfaces. Also, in this section, we present the method used to obtain the background electronic density at the position of the moving adsorbates at each time step of the dynamics, which is required to obtain the friction coefficients that describe the coupling of the adsorbates to the heated electronic system in the Langevin equation of motion. Subsequently, in Section 4, the molecular dynamics simulations performed in the precalculated 0.75 ML-CO/Pd(111) EANN-PES to model the laser-induced desorption of CO from the Pd(111) surface are presented. The results of the simulations are compared to the (T e ,T l )-AIMDEF results of ref 27 showing, conclusively, the validity of our EANN-PES to perform molecular dynamics in these unprecedentedly exigent conditions. Finally, in Section 5, the main conclusions of the work are summarized.

EANN-PES GENERATION AND QUALITY CHECK
The analytical representation of the adiabatic PES E({r k }) ruling the desorption of CO from Pd(111) with 0.75 ML coverage is calculated with the recently developed embedded atom neural network (EANN) method. 47 Similar to AtNN by Behler and Parrinello, 44 the total energy of an N atom system is expressed as the sum of the energy of each atom that conforms it, E i ({r k }). In the EANN framework, E i ({r k }) is described in terms of the electronic embedding density, i.e.
where NN i is the species-dependent atomic neural network of the ith atom in the system that depends on the embedding density vector ρ i whose components represent the local electron density provided by the surrounding near atoms. Specifically, the set of local density components defining ρ i are given by are Gaussian-type orbitals (GTOs) centered at each of the n a atoms j that are located within a radius r c from the embedded atom i. In these equations, r i = (x i , y i , z i ) and r j = (x j , y j , z j ) are the Cartesian position vectors of atoms i and j, respectively, with r ij = |r i − r j |; α and r s determine the width and the center of the Gaussian-like term in eq 3, respectively (and thus control the shape of the radial distribution related to each GTO), while l x , l y , and l z are the values of the orbital angular momentum in each coordinate, whose sum equals the total angular momentum L, i.e., L = l x + l y + l z . In eq 2, f c (r ij ) is the commonly used cosine type cutoff function that makes the interaction to smoothly decay to zero as r ij approaches the cutoff radius r c , 44 and c j is the element-and orbital-dependent weight 50 53 and the AIMDEF module 20−26 that was extended to include excited electrons and excited phonons effects through time-dependent electronic and lattice temperatures (see Section 3). 27 A total of 100 trajectories were calculated in a four-layer (4 × 2) supercell that contained six CO adsorbates equally distributed among atop, hcp, and fcc sites and 32 Pd atoms describing the Pd(111) surface. Note in passing that the use of such a large cell (twice the unit cell for 0.75 ML) was aimed to account for out-of-phase movements of the coadsorbed CO molecules, providing a more realistic description of the interadsorbate interactions. Each trajectory was integrated up to 3.5−4 ps using a time step of 1 fs. Altogether, the whole data set consists of 352 505 configurations, each defined by the position vectors of the 44 atoms in cell {r i }, for which the corresponding DFT potential energies E DFT and DFT forces on each atom F i DFT are well characterized. The training process can be greatly simplified if, instead of using such a huge amount of data, we select a smaller subset that correctly represents the relevant configurational space of the system. There are different properties that can guide this selection, e.g., distance of the gas species to the surface, 31 forces, 37 etc. In our case, the photoinduced desorption dynamics is characterized by trajectories that yield none, one, or two desorbing CO, each conceivably providing information on the interaction at variable coverages. Thus, an initial subset is constructed with 4500, 6000, and 4500 configurations that are randomly selected from the set of trajectories with zero, one, and two desorption events, respectively. The ulterior analysis of the covered potential energy range served us to validate that the configurational space probed in the AIMDEF simulations is in principle well sampled. Figure 1 shows that the E DFT distributions of the whole and of each of the three types of desorbing trajectories (top panel) are accurately reproduced by the subset of 15 000 trajectories (bottom panel). It is worth mentioning that the values E DFT in Figure 1 are the output of the VASP program. As shown in this figure, the relevant potential energy variation range in the dynamical configurational space covers around 12 eV, which corresponds to around 0.3 eV per moving atom (see below).
The analytical EANN-PES for the 0.75 ML-CO/Pd(111) system uses 60 density descriptors for each atomic species that correspond to take L = 0−3 combined with 15 Gaussian functions with α = 0.93 Å −2 and r s varying within the interval [0,r c ] in increments Δr s = 0.46 Å for our chosen value r c = 6.5 Å. Different architectures were tried, but an optimal balance between small errors in energy and the required computation time is achieved using two hidden layers with 60 neurons each for every atomic NN i . The EANN code takes advantage of the usual random separation of the reference data set into training and validation subsets to produce different NNs during the same run. In our case, five EANN-PESs have been trained using 90 and 10% of the configurations as training and validation subsets, respectively. An efficient extreme machine learning Levenberg−Marquardt algorithm is used in the optimization of the fitting parameters (c j ). 54 The cost function used to evaluate at each iteration the quality of an EANN-PES involves the energies and atomic forces. As discussed in ref 47, convergence is very efficient in the EANN method. For our settings, the required accuracy is achieved in the five training PESs in less than 50 iterations, being the root-mean-squareerrors (RMSEs) in the energy per moving atom of 0.43−0.58 Since the photoinduced desorption dynamics of interest implies that the system is exposed to extreme conditions characterized by high surface temperatures (1000 K) and highly excited adsorbates, we consider it important to evaluate the accuracy of the PES in predicting not only the energy but also, especially, the atomic forces. A new data set formed by 87 382 configurations randomly taken from the AIMDEF data set not used in the fitting process and representative of the three types of desorbing trajectories is used for this purpose (predict data set). The results show that the obtained PESs are excellent with the RMSE in energy of only 0.86−0.95 meV per moving atom and the RMSE in the Cartesian components of the forces of 0.05−0.06 eV/Å. Nevertheless, we observe that the maximum errors in the forces are in some cases large (between 1.1 and 5.9 eV/Å depending on the coordinate and trained EANN-PES considered). Thus, new configurations are added to the initial set of 15 000 structures to improve the quality of our EANN-PES. The criterion to select these new configurations is as follows. First, among the five trained EANN-PESs, we select the one with the smallest maximum absolute errors in the atomic forces, i.e., where β = x, y, z refers to the force component. From this analysis, we select the 10 largest errors for each of the 36 moving atoms (i.e., 12 atoms forming the six CO adsorbates and 24 Pd atoms that correspond to the moving three of the four layers describing the surface). Figure 2 shows the distribution of these maximum errors for each of the force components as obtained with the EANN-PES with the smallest errors in the forces (black histogram bars). The corresponding configurations (883 because errors in different force coordinates can occur for the same configuration) are added to the original 15 000 data input to develop five new EANN-PESs, using the same NN settings (i.e., basis set and architecture).
The quality of the five newly developed EANN-PESs is impressively good. Since the accuracy in the forces is also crucial for reproducing the system dynamics and thus its physical behavior properly, we select the EANN-PES with the smallest errors in the atomic forces. In this PES, the RMSE in the energies per moving atom is 0.41 and 0.87 meV for the training and validation sets, respectively. Its accuracy is further confirmed by the energy results obtained with the predict data set. As shown in Figure 3 (left panel), the maximum error per moving atom is 8.27 meV but the small RMSE of 0.85 meV marks the minor error introduced in most of the configurations. The error distribution plotted in Figure 3 (right panel) clearly shows that this is the case.
Regarding its accuracy in predicting the atomic forces, the RMSE is not greatly improved in the retraining process and it is still around 0.05 eV/Å, whereas the mean absolute errors are around 0.04 eV/Å. However, a comparison of the 10 maximum absolute errors for each atomic force component |ΔF β | obtained with our final EANN-PES (red bars in Figure 2) with the ones of the first training process (black bars) shows a remarkable reduction in the errors. The reduction in the maximum errors is as follows: |ΔF x | ≈ 1.13 → 0.66 eV/Å, |ΔF y | ≈ 1.04 → 0.58 eV/Å, and |ΔF z | ≈ 2.56 → 0.62 eV/Å. This improvement is also noticeable when comparing the mean absolute error (taken over the 10 maximum errors for each of the 36 moving atoms) of F z , plotted in black and red histograms, that decreases from |Δ | ≈ As an additional stringent test, we also show the maximum error predictions of the final EANN-PES on the 352 505 configurations forming the full AIMDEF data set (green bars of Figure 2). In this case, we obtain |ΔF x | = 0.68 eV/Å, |ΔF y | = 0.58 eV/Å, and |ΔF z | = 0.63 eV/Å so that the maximum absolute-valued errors show none to a very little increase with respect to the predictions of the red histogram. The corresponding mean absolute error values  The results from all these validation tests imply that we should also observe good convergence regarding any atom of the CO/Pd(111) cell. Figures 4 and 5 show two examples of these good agreements between F β EANN and F β DFT , β = x, y, z, for two representative atoms of the system: one C atom from a CO at a fcc site and one O from a hcp site-placed CO. As a result, after the full static analysis, we are confident about this final EANN-PES's robustness, which we further determine by trying to reproduce the dynamics of the CO/Pd(111) system as described in ref 27.

PHOTOINDUCED MOLECULAR DYNAMICS ON SURFACES
The desorption of CO from the Pd(111) surface induced by femtosecond laser pulses is simulated with molecular dynamics with electronic friction and thermostat [(T e ,T l )-MDEF] calculations performed in our developed EANN-PES. To do so, we have modified our implementation of the (T e ,T l )-AIMDEF methodology 27 to compute trajectories in which forces and electronic friction coefficients are evaluated at each visited configuration from an arbitrary PES and an arbitrary electronic density generator function (DGF), respectively. In Figure 6, we show a concise scheme of the (T e ,T l )-(AI)MDEF model. Macroscopically, a laser pulse heats the electrons of the metal surface, which subsequently excite surface phonons. This response is modeled with the same 2TM as in the original (T e ,T l )-AIMDEF calculations. The macroscopic response of the system is thus followed microscopically by CO and Pd degrees of freedom (DOFs) with different equations of motion.
CO DOFs are subjected to the following Langevin equations of motion where m i , r i , and η e,i are the mass, position vector, and the electronic friction coefficient of the ith atom conforming the set of adsorbates, respectively. 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. In the (T e ,T l )-AIMDEF simulations of ref 27, this force was calculated on the fly using the Hellmann−Feynman theorem. In the (T e ,T l )-MDEF calculations that we present here, this force is calculated using the precalculated EANN-PES. The third term R e,i is the random fluctuating force that mimics the effect of the hot metal electrons on the adsorbates, and it is responsible for the excitation of the latter by the excited electronic system. This force is related to the electronic friction force (second term) through the fluctuation−dissipation theorem. Specifically, R e,i is modeled by Gaussian white noise with variance where k B and Δt are the Boltzmann constant and the timeintegration step, respectively.  The equations of motion followed by the Pd surface atoms incorporate two kinds of force terms. The first term consists of the adiabatic force that, as in the case of adsorbates, depends on the position of all atoms. Again, in the (T e ,T l )-AIMDEF simulations of ref 27, this force is calculated on the fly using the Hellmann−Feynman theorem, whereas in the (T e ,T l )-MDEF calculations, it is computed using our EANN-PES. The second term accounts for the heating of the surface layers due to the electronic excitation generated by the laser pulse. As done in ref 27, this effect is simulated by coupling the atoms of the two topmost surface layers to a thermal bath described by a Nose−Hoover thermostat, 55,56 in which the temperature is the time-dependent temperature T l (t) obtained from the 2TM. Thus, the equations of motion of these surface atoms j with mass m j and position vector r j are as follows where N is the number of atoms of the first two layers in our simulation cell, Q is a parameter with dimensions of energy × time 2 that acts as the mass of the dynamical variable s, and ξ = Q −1 sp s is the thermodynamic friction coefficient. 57 Finally, the third layer is used as a transition region between the hot surface and the inner (not heated) bulk. Therefore, the movement of each atom k in the third surface layer is described by the classical Newton equations of motion and the adiabatic approximation where m k and r k are its corresponding mass and position vector, respectively. The fourth-layer atoms are kept frozen in our simulations.
Regarding the coupling of the adsorbates to the electronic system via the second and third terms in the right-hand side of eq 4, the electronic friction coefficient for each atom in each adsorbate η e,i is calculated within the original local density friction approximation (LDFA) 19,58 that in the case of molecular adsorbates treats the molecule as formed by independent atoms (independent atom approximation, IAA). This means that η e,i depends on the value of the bare surface electronic density at the position of the atom i forming the adsorbate, i.e., n sur (r i ). In the (T e ,T l )-AIMDEF simulations, 27 this is achieved using the Hirshfeld partitioning scheme 59 to subtract the contribution of the adsorbates from the selfconsistent electronic density that is calculated at each time step. 22,23 In the case of (T e ,T l )-MDEF, apart from an EANN-PES, the friction coefficients must also be obtained from precalculated results. Moreover, to use an EANN-PES efficiently, it is clear that the friction coefficients have to be calculated at least as fast as the potential. One could use a similar NN scheme to interpolate the density or the electronic friction coefficients based on the AIMDEF data set. 60 Actually, even a more complicated object such as the electronic friction tensor can be interpolated by NNs. 61−63 However, here we choose to employ a more simple approach that we show to be accurate and compared to a NN interpolation, faster and more stable.
We start by writing the Pd(111) electron density at the position of each atom of each CO adsorbate n(r C,O ) as a sum of the electronic densities contributed by individual Pd atoms at this position We then assume that ρ(|r C,O − r Pd |) can be described by two exponentially decaying functions Using AIMDEF data points of n(r C,O ), one can fit the four parameters: a, b, c, and d. We did also try other functions such as Gaussian functions and their combination with the exponential function, but two decaying exponential functions gave the best results retaining the simplicity and a small number of parameters. Since it is easy to evaluate such a function, we decided to use all available (T e ,T l )-AIMDEF surface electronic densities from which the corresponding friction coefficients were calculated. The fitting procedure results in the following parameters: a = 3.15975 au, b = 4.25214 Å −1 , c = 0.29080 au, and d = 2.52252 Å −1 (au stands for atomic units).
Once the density n(r C,O ) is known, the friction coefficients are calculated within the LDFA, 19,58 as it was also done in the (T e ,T l )-AIMDEF simulations. In particular, the C and O friction coefficients are individually fitted to the following analytical function In Figure 7, we compare the LDFA friction coefficients obtained with the fitted electron density n(r C,O ) to their corresponding (T e ,T l )-AIMDEF friction coefficients. It can be seen that both the maximum error of 0.0086 au and RMSE of 0.0017 au are small, especially compared to the errors associated with different calculations of embedding densities. 23 In addition to accuracy, the advantage of this approach is that it is extremely fast in evaluating friction coefficients compared to the potential evaluation. Lastly, by the construction of two exponentially decaying functions, it is also ensured that there is no problem of overfitting or large errors in the extrapolation range that could occur with approaches such as NN. All this guarantees that eq 9 together with eq 10 constitutes an adequate DGF to be used in the (T e ,T l )-MDEF simulations to calculate, at each time step, the friction and stochastic forces in eq 4.

RESULTS OF THE DYNAMICS SIMULATIONS
We have hitherto shown how our best developed EANN-PES (Section 2) and our atomwise additive electronic density function approach (Section 3) yield results in being close to the original 0.75 ML-CO/Pd(111) AIMDEF data set. Despite this being a reasonable quality check, it is not sufficient to prove that using both functions in MDEF calculations produces exactly the same results as in AIMDEF. This is due to two possible undesirable situations. 46 First, MDEF trajectories may sample regions of the configurational space that, even when close to the EANN training set or the DGF data set, are not yet correctly described by the optimized parameters, and second, these trajectories may enter regions of the configurational space far from the AIMDEF data set where both fitted functions have to extrapolate information. In addition, the probability of encountering any of the described problems is increased by the high number of DOFs, in our case 108, that play a role during this type of laser-induced desorption dynamics simulations.
To rule out these sources of error and further address the accuracy of our developed EANN-PES and additive DGF, we have simulated the desorption of CO from 0.75 ML-CO/ Pd(111) with (T e ,T l )-MDEF, assuming the same experimental conditions as in the original (T e ,T l )-AIMDEF calculations. 27 Specifically, the Pd(111) being initially at 90 K is heated by a laser pulse of 780 nm wavelength, 100 fs of full width at halfmaximum (FWHM), and absorbed fluence F = 13 mJ/cm 2 . In practice, the heated surface is modeled with the same timedependent temperatures T e (t) and T l (t) used in the (T e ,T l )-AIMDEF simulations, 27 in which the maximum of the laser pulse arrives at the instant t = 410 fs.
Using these conditions, two sets of dynamics calculations have been run. In the first one, hereinafter called (T e ,T l )-MDEF-1, we have used the same set of 100 initial configurations as utilized in ref 27. This allows us to compare step by step the extent to which (T e ,T l )-MDEF and (T e ,T l )-AIMDEF trajectories are similar and detect if there are artifacts in the potential or electronic friction coefficients that could push trajectories away from the EANN-PES confidence zone, even when starting from (T e ,T l )-AIMDEF initial configurations. In the second set of calculations, hereinafter called (T e ,T l )-MDEF-2, we have used 2000 configurations selected randomly from a set of 10 000 structures generated by letting the 100 initial configurations of ref 27 evolve in 1 ps with a constant temperature of 90 K. This enables us to test how robust (T e ,T l )-MDEF dynamics is when trajectories start from configurations not included in any of the data sets used to fit the EANN-PES or DGF.
We find that (T e ,T l )-MDEF-1 trajectories lie close to their (T e ,T l )-AIMDEF counterparts. In fact, configurations visited during the first 400−500 fs are practically identical to those visited by the original trajectories. To illustrate this finding, in Figure 8 we show how a typical (T e ,T l )-MDEF-1 trajectory (blue line) compares with its corresponding (T e ,T l )-AIMDEF trajectory (green line) in terms of time evolution of C atom friction coefficients during 1 ps. The agreement for times below 500 fs is very high, independent of the initial position of the C atoms. After 500 fs, friction coefficients start to diverge as trajectories commence to follow different pathways due to cumulative differences in the dynamics. This behavior is shared among all (T e ,T l )-MDEF-1 simulations and reflects that the quality of our fitted EANN-PES and DGF is good enough to keep (T e ,T l )-MDEF trajectories close to their (T e ,T l )-AIMDEF analogues in a very detailed way for several hundred femtoseconds. This is a remarkable quality achievement as the agreement was attained despite the possible dynamical instability, the stochastic nature of the Langevin equations of motions, and the different integration time steps used in the simulations (Δt = 0.2 fs in MDEF, Δt = 1 fs in AIMDEF) for an amount of time in which the electronic temperature of the 2TM changes from 90 K (t = 0 fs) to 5500 K (t = 500 fs).
In Figure 9, we show (T e ,T l )-MDEF-1 (blue line), (T e ,T l )-MDEF-2 (red line), and (T e ,T l )-AIMDEF (black line) results for total (top panel) and initial site-dependent (bottom panels) CO desorption probabilities as a function of time. In each case, the desorption probability p CO is defined as the (cumulative) number of CO molecules desorbed at a given time divided by the total number of trajectories (N traj ) and the total number of Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article CO molecules in our simulation cell (six; see Figure 6). Since this definition is equivalent to consider that the number of CO desorption events distribute as a binomial distribution with six N traj independent trials and success probability p CO , we associate to each evaluation of p CO a confidence interval calculated with Wilson asymmetric score intervals 64 for the confidence of 99%. It is apparent that all desorption probabilities calculated with (T e ,T l )-MDEF-1 and (T e ,T l )-MDEF-2 are in good agreement with the (T e ,T l )-AIMDEF results, especially in the case of total and fcc-CO desorption probabilities. On a closer look, we can see that site-dependent (T e ,T l )-MDEF-1 desorption probabilities present more discrepancies with (T e ,T l )-AIMDEF than (T e ,T l )-MDEF-2. In particular, (T e ,T l )-MDEF-1 top-CO desorption ratios tend to be higher than those of (T e ,T l )-AIMDEF, whereas hcp-CO desorption ratios tend to be lower. These deviations are mostly from the statistical variability associated with their respective 100 trajectories' ensembles, since all (T e ,T l )-MDEF-1 and (T e ,T l )-AIMDEF desorption probabilities are within their 99% confidence intervals reciprocally. In the case of (T e ,T l )-MDEF-2, all desorption probabilities lie closer to those of (T e ,T l )-AIMDEF despite not sharing exactly the same initial configurations. This supports that our EANN-PES and DGF are accurate enough to describe the dynamic energy barriers that CO molecules encounter during the laser-induced desorption process on the same footing as the original AIMDEF calculations. Compared to (T e ,T l )-AIMDEF, the improved statistics of the (T e ,T l )-MDEF-2 simulations allow us to establish with more precision the time scale of the desorption process. It takes more than 0.6 ps since the arrival of the laser-pulse maximum (occurring at t = 410 fs) to observe the first desorption events. The site-resolved desorption probabilities show that these initial events correspond to atop-CO adsorbates. Around an additional 500 fs are required for desorption from either the hcp or fcc sites.
Having demonstrated the accuracy of our methodology to reproduce the (T e ,T l )-AIMDEF trajectories for the first few hundred femtoseconds and the final outcome in terms of CO desorption probabilities, we now focus on the time evolution of the kinetic energy of CO molecules, which is more sensitive to the fine details of the paths followed on the PES. In Figure  10 (left panel), we plot the mean translational (E kin trans , reddish, full thick lines) and mean rovibrational (E kin rovibr , bluish, full thick lines) kinetic energy of adsorbed CO molecules averaged over trajectories as a function of time. We define E kin trans as only the mean translational energy of the center of mass of CO molecules, while E kin rovibr is calculated as the mean of C and O total kinetic energy minus E kin trans . Starting with E kin trans , we observe that the (T e ,T l )-MDEF-1 (red) and (T e ,T l )-MDEF-2 (light red) results lie very close to the (T e ,T l )-AIMDEF (dark red) remarkably well with the (T e ,T l )-AIMDEF E kin rovibr values (dark blue). But together with the mean energy values, it is important to confirm whether the instantaneous energy distributions of the adsorbates are also well reproduced by our EANN-PES. The latter can be estimated by comparing, for instance, the associated standard deviations of the translational ΔE kin trans and rovibrational ΔE kin trans kinetic energy distributions. Thus, there are two additional (thin dotted) curves representing the values E kin trans(rovibr) ± ΔE kin trans(rovibr) associated to each mean kinetic energy curve in Figure 10. As shown in the figure, there is also a good agreement between the (T e ,T l )-AIMDEF instantaneous standard deviations and those obtained from the (T e ,T l )-MDEF-1 and (T e ,T l )-MDEF-2 simulations, which reinforces the high quality of the EANN-PES.
In Figure 10 (right panels), we also show the empirical cumulative distribution function (ECDF) of ensembleaveraged mean translational (top-right panel) and mean rovibrational (bottom-right panel) kinetic energies of desorbed CO molecules. These energy distributions are defined in the same way as in the previous panel, but they are only calculated when desorbed CO molecules are far away from the surface (center-of-mass-to-surface-plane distance greater than 6 Å) at the end of the dynamics (3500 fs). We have marked with shaded areas around each ECDF their associated 99% Dvoretzky−Kiefer−Wolfowitz confidence intervals. 65 These areas define for each empirical distribution the region inside which their exact cumulative distribution functions lie with a 99% confidence. From these panels, we can see that (T e ,T l )-MDEF-1 (blue line) and (T e ,T l )-MDEF-2 (red line) evaluated mean energies agree considerably with (T e ,T l )-AIMDEF (black) results, as all MDEF curves lie within the confidence interval of the original calculation. It is also apparent that the (T e ,T l )-MDEF-2 confidence intervals are much smaller than the others (and therefore closer to the exact distribution), due to the higher number of trajectories sampled. These results further show that both, rovibration and translation of CO molecules on top of the surface and after desorption, are in good consonance with the findings in ref 27.
Taking advantage of the improved statistics provided by (T e ,T l )-MDEF-2, we can reliably determine the rovibrational state of the desorbed CO molecules (ν f , j f ). In performing such an analysis, the rotational quantum number j f is computed from the classical angular momentum L cl as the closest integer , while the vibrational state ν f is determined from the vibrational action α f as the nearest integer that verifies ν f = α f /h − 1/2, where h is the Planck constant. 66 The results from this quasi-classical analysis show that most of the desorbed molecules (83.3%) are in the vibrational ground state and slightly rotationally excited (j f ≤ 20).

CONCLUSIONS
Using the recently developed embedded atomic neural network (EANN) 47 framework, we generate an accurate and complex potential energy surface that is able to describe the dynamics of the femtosecond laser-induced desorption of CO from Pd(111) with a coverage of 0.75 ML. As the training data set, we use around 16 000 configurations taken from ab initio molecular dynamics simulations [(T e ,T l )-AIMDEF] 27 that incorporate both the effect of the laser-excited electrons and the concomitant excitation of the surface phonons. The   65 For all panels, (T e ,T l )-AIMDEF represents results extracted from ref 27 calculations, (T e ,T l )-MDEF-1 represents dynamics results obtained from our best trained potential using the same 100 (T e ,T l )-AIMDEF initial conditions, and (T e ,T l )-MDEF-2 represents dynamics results obtained with the same EANN potential using 2000 random initial conditions. Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article otherwise random selection of AIMDEF configurations has been biased so that the final training set incorporates a defined proportion of configurations from trajectories characterized by a different number of CO desorption events (in practice, none, one, and two in the proportion 3:4:3, respectively). This procedure enforces the training set to contain an equilibrated amount of information about different desorption channels and different instantaneous surface CO coverages. The quality of our EANN-PES in predicting not only the system energy but also the atomic forces is validated against a huge set of almost 90 000 AIMDEF configurations not included in the training process. An impressively small RMSE of 0.85 meV is achieved in the energy per moving atom and of around 0.05 eV/Å for the forces. The EANN-PES robustness is further checked by evaluating its performance in simulating the desorption dynamics of CO from the Pd(111) surface covered with 0.75 ML. Following the (T e ,T l )-AIMDEF method, the excitation created by the laser in the electronic system and the subsequent excitation of the surface phonons is represented within the two-temperature model as two coupled thermal baths with temperatures T e (t) and T l (t). Next, the coupling of each adsorbate to the highly excited electrons is described by a Langevin equation that depends on the time-dependent electronic temperature T e (t), and the Nose−Hoover thermostat is used to describe the timedependent temperature of the surface atoms T l (t). In all these equations, the adiabatic forces over each adsorbate and surface atom, which are obtained from the Hellmann−Feynman theorem in AIMDEF, are here calculated by means of the 0.75 ML-CO/Pd(111) EANN-PES. The value of the density at the position of the adsorbates at each time step, necessary to evaluate the friction coefficients entering the nonadiabatic forces of the Langevin equation within the LDFA scheme, is obtained with an efficient DGF based on the fitting of the (T e ,T l )-AIMDEF friction coefficients. In this way, the (T e ,T l )-AIMDEF results are reproduced with a remarkable level of accuracy. This demonstrates the outstanding performance of the obtained EANN-PES that can cover an extensive range of surface temperatures (90−1000 K); a large number of degrees of freedom, those corresponding to multiple adsorbates and surface atoms, i.e., 108 in our simulation cell; and an extremely complex configurational space, characterized by very mobile adsorbates and a changing CO coverage caused by the sizable desorption.
Application of this NN-PES for future computational tests of system dynamics under different initial conditions should be straightforward. Up to now, molecular dynamics simulations at the ab initio-DFT level of femtosecond laser pulse-induced reactions at surfaces have been limited either by the impossibility to account for all the relevant degrees of freedom of the system or, in the case of (T e ,T l )-AIMDEF calculations, by the huge computational cost of AIMD that restricted the achievable statistics. In this respect, the usefulness of NN-PES in the present work cannot be underestimated. For instance, it will allow us to perform simulations, at the level of (T e ,T l )-AIMDEF with low computational costs and higher statistics, for different laser fluences and adsorbate coverages. Also, the theoretical study of the time delay dependence of two pulse correlation experiments will be tractable. It will be also possible to extend the simulation times to values much larger than the usual 2−4 ps in (T e ,T l )-AIMDEF to guarantee that the computed reaction/desorption probabilities are saturated. The increase of the simulation cell to avoid finite-size effects and to account more accurately for interadsorbate energy exchange will also be easily accessible. Indeed, it opens the path to, from a theoretical point of view, fully characterize and understand these kinds of experiments. Last but not least, it must be stressed that our work constitutes a strong support for utilization of the EANN methodological framework for the development of accurate NN-PESs for other complex gas− solid interfaces.

■ ACKNOWLEDGMENTS
The authors acknowledge financial support by the Gobierno Vasco-UPV/EHU Project no. IT1246-19 and the Spanish Ministerio de Ciencia e Innovación [Grant no. PID2019-107396GB-I00/AEI/10.13039/501100011033]. This work has been supported in part by the Croatian Science Foundation under project UIP-2020-02-5675. This research was conducted in the scope of the Transnational Common Laboratory (LTC) "QuantumChemPhysTheoretical Chemistry and Physics at the Quantum Scale". Computational resources were provided by the DIPC computing center.