Surface electron density models for accurate ab-initio molecular dynamics with electronic friction

D. Novko, M. Blanco-Rey, 1 M. Alducin, 1 and J.I. Juaristi 3, 1 1 Donostia International Physics Center (DIPC), Paseo Manuel de Lardizabal 4, 20018 Donostia-San Sebastián, Spain 2 Departamento de F́ısica de Materiales, Facultad de Qúımicas UPV/EHU, Apartado 1072, 20080 Donostia-San Sebastián, Spain 3 Centro de F́ısica de Materiales CFM/MPC (CSIC-UPV/EHU), Paseo Manuel de Lardizabal 5, 20018 Donostia-San Sebastián, Spain (Dated: June 14, 2016)


INTRODUCTION
Many are the theoretical studies confirming that the fundamental properties in most elementary gassurface processes are satisfactorily described by the Born-Oppenheimer approximation 1,2 . Common to all these studies is the use of an accurate adiabatic potential energy surface (PES) calculated with density functional theory (DFT) that at least accounts for the degrees of freedom of the gas species involved in the event.
Still, the challenge in present thermal and hyperthermal gas-surface simulations is to provide a reliable description of the two main energy exchange channels that may affect the dynamics and reactivity of gas-phase species on solid surfaces, namely, phonon excitations and electron-hole (e-h) pair excitations. In the end, these are the mechanisms that dictate the thermalization rate and, hence, the mean traveled length of the nascent adsorbates. Even more generally, these mechanisms are expected to contribute actively in any gas-surface process that involves strong and long-lasting interactions. The promising femtosecond laser induced chemistry [33][34][35][36][37][38] is another good example of it. There are various theoretical studies showing that phonons and, particularly, e-h pairs are the driving ingredients in many cases [39][40][41][42][43][44][45] .
In the past, phonons or, more precisely, the effect of energy exchange with the lattice in gas-surface dy-namics have been reasonably described by using thermostats coupled to an adiabatic PES that paradoxically neglects the individual surface atoms degrees of freedom 8,18,31,32,[46][47][48][49][50][51] . The latter can become a serious limitation when large energy exchange and long interaction times are at work, since the distortions created locally on the surface can continuously modify the PES. In this respect, ab-initio molecular dynamics (AIMD), which based on DFT and the Hellmann-Feynman theorem allows incorporating the surface atoms movement, is nowadays the state-of-the-art methodology to account for the aforementioned phonon excitations effects [52][53][54][55][56][57][58] . The QM/Me model developed by Meyer et al. 59 has been recently proposed as an improvement over the usual AIMD method because it avoids the spurious periodic distortions that may appear in AIMD in case of using too small surface unit cells. Notably, there are also theoretical studies that include a quantum treatment of the phonon excitations 60 , but in those cases the gas-surface interaction is described through simplified model potentials 61,62 .
While searching for an accurate and joint description of the electronic and phononic energy dissipation channels, the recently developed AIMD with electronic friction (AIMDEF) method 63 that is based on the local density friction approximation (LDFA) 64 , constitutes a promising tool to meet this goal. However, the original AIMDEF of Ref. 63, which is based on the rigid surface electron density, may fail in cases of large surface atoms displacements that cause non-negligible changes in the surface electron density.
In this paper, we propose and analyze different meth-ods to successfully overcome this limitation. Having these powerful AIMDEF approaches, we apply them to investigate a central issue in gas-surface interactions, that is, the adsorption and relaxation of hot gas species on metal surfaces. More precisely, we investigate three different adsorption scenarios that cover a representative range of adsorption energies E ads and adsorbateto-surface mass ratios. Specifically, dissociated H 2 on Pd(100), N on Ag(111), and N 2 on Fe(110). These systems were recently analyzed in Refs. 65 and 66, using one of our proposed methodologies. Here, we extend and present a more detailed analysis on how phonons and e-h pair excitations depend on each other. We find that the electronic mechanism can be noticeably sensitive to the phononic one, in particular, in those cases where large surface atoms displacements are likely to occur. Nevertheless, we also find that the adsorption and relaxation processes themselves are not much affected by the details of these excitations, at least for the typical energies of few eVs that matter in practical gas-surface reactions. Therefore, our new results also confirm the robustness of the conclusions stated in Ref. 65, namely, (i) that the adsorption of light and heavy gas species are dominated by the electronic and the phononic excitations, respectively, and (ii) that independently of the gas species considered, the electronic mechanism is crucial during the final accommodation of the adsorbate on the adsorption well.
The outline of the article is as follows. Section II starts with the basics of the LDFA of Ref. 64 and continues with a detailed description of the surface electron density models we propose to use in AIMDEF simulations. The section ends by analyzing the performance of each density model under extreme reliable conditions of large surface atoms displacements as those occurring upon equilibration of the adsorbates on the adsorption wells. Their behaviour in representative AIMDEF simulations is discussed in Sec. III. In particular, we analyze how the adsorption probabilities and energy dissipation into e-h pair and phonon excitations depend on the surface density description for the above mentioned three adsorption scenarios. The summary and conclusions are given in Sec IV.

II. LOCAL DENSITY FRICTION APPROXIMATION FOR MOVING SURFACE ATOMS
Non-adiabatic effects that come from the energy exchange between the nuclear and electronic degrees of freedom can be effectively included in the nuclei classical equations of motion in terms of a friction force 67 . The crucial point is to determine a realistic value of the friction coefficient η(r i ) acting on the gas-atom i at each point r i along its trajectory. In this respect, different theoretical studies have been performed during the last years treating this issue 64,68,69 . The local density friction approximation (LDFA) 64 is one of the formalisms that, in spite of its simplicity, captures the relevant physical aspects of the low energy e-h pair excitations [70][71][72][73] as those created by slowly moving gas species. This is one of the reasons of being widely applied to study the effect of electronic excitations in the dynamics of atoms and molecules on metal surfaces 14,28,44,45,47,48,57,[63][64][65][66][74][75][76][77][78] . More recently, the LDFA has been shown to accurately describe the electronic energy loss in the scattering of H from Au(111) 79 . The LDFA assumes that η(r i ) is equal to the friction coefficient that the same atom i would have in case of being moving within a homogeneous free electron gas (FEG) of density n 0 = n sur (r i ), where n sur (r i ) is the electron density of the bare metal surface at the position r i .
The friction coefficient of a slowly-moving atom inside a FEG was derived, using quantum scattering theory, from the energy that loses per unit path length. Its expression in atomic units (a.u.) reads 70,71,80 , where k F is the Fermi momentum and δ l (k F ) are the scattering phase-shifts at the Fermi level that are calculated from the DFT scattering potential of an atom within the FEG. The latter turned out to be a crucial step to reproduce available experimental data on the stopping power of atoms and ions in metal solids and surfaces [70][71][72][73] .
For the case of a molecular projectile the original LDFA of Ref. 64 calculates the friction coefficient on each atom in the molecule as if they were non-interacting atoms (independent atom approximation, IAA). The latter has been shown to be a reasonable approximation for the translational degrees of freedom 81 , but it may introduce errors when treating the coupling of the molecular vibrational movement with the metal electrons in situations of strong and long-time molecule-surface interactions. The extreme realization of the latter conditions is undoubtedly the electron-vibration coupling that affects the lifetime of the molecular stretching mode when adsorbed on metals 82,83 . In the recent study of Ref. 83, the authors use the Fermi golden rule formulation 84 and also analyze the anisotropies of the friction tensor in nonuniform systems that cannot be captured by the isotropic LDFA. We note, however, that a direct quantitative comparison between these results and those of the LDFA is not straightforward because in the calculation of the friction tensor only electronic transitions that conserved the crystal momentum were included (i.e. initial and final electronic wavevectors are equal, k = k ), whereas the calculation of the friction tensor for a single adsorbate requires the Fermi surface integration over transitions not conserving the crystal momentum.
The LDFA was first used in AIMDEF simulations to study the relaxation of the hot H atoms formed upon dissociation of H 2 on Pd(100) 63 and, more recently, to investigate quantum-size effects in the vibrational lifetime of H adsorbed on Pb films 76 . In both cases, the bare surface electron density, which determines the friction coefficients η(r i ) within the LDFA, was approximated by the electron density of the bare frozen surface (FS) calculated self-consistently with DFT, which will be denoted as n FS sur in the following. Although the use of n FS sur is only justified in simulations where the surface atoms are fixed at the equilibrium positions, it is still a reasonable approximation in those cases where the surface atoms are barely moving 63,76 . However, in most cases the surface atoms displacements are expected to cause appreciable changes in the surface electron density n sur . Thus, the latter needs to be known at any instant t and it complicates the use of the LDFA in usual AIMDEF simulations because only the electron density of the whole system, i. e., gas species and surface atoms, is calculated selfconsistently at each integration step.
Here, we introduce two methods that facilitate the applicability of the LDFA for moving surface atoms. In the first one, the surface electron density is calculated at each time step t as the superposition of the ground state electron densities of the isolated individual surface atoms n atom j , i.e., where the summation index j runs over all surface atoms N sur . This method successfully accounts for the movement of the surface atoms at each time step, but it obviously misses the charge redistribution upon formation of bonds between the surface atoms. The second method corrects this misbehavior by making use of the Hirshfeld partitioning scheme 85 , which has been successfully applied to study the vibrational lifetimes of molecular adsorbates within the LDFA framework 82 . Here, we use it in order to subtract the contribution of the gas-phase atoms from the self-consistent density of the whole system n SCF (r i , t). More precisely, the bare surface electron density is approximated at each t by where the indexes m and n run, respectively, over the total numbers of atoms in the system N and in the adsorbate N A . In this equation, the Hirshfeld weighting factor w n (r i , t) represents the contribution of the n-th atom to the electron density of the whole system at r i . Thus, the factor 1 − NA n=1 w n (r i , t) defines the weight corresponding to the system without the contribution of the adsorbate.
The described electron density methods (n FS sur , n AS sur , and n H sur ) have been implemented in the Vienna Ab-initio Simulation Package (vasp) 86 to perform AIMDEF calculations 63,65 . A. Performance of the proposed surface density models We start by examining the adequacy of n FS sur , n AS sur , and n H sur in describing the self-consistent bare surface electron density n sur once the adsorbates are fully relaxed and accommodated on the surface. The latter constitutes one of the possible real extreme conditions under which the surface density can be significantly altered as a consequence of the charge redistribution between the adsorbate and the surface. In fact, this redistribution can even distort locally the surface lattice.
All DFT calculations presented in this work were done using the vasp package, which uses a plane wave basis set. The core-electron interaction was approximated by projector-augmented wave potentials 87 . Following previous works 48,49,63 , a different generalized gradient approximation (GGA) for the exchange and correlation functional was used in each of the studied systems, namely, the Perdew-Wang 1991 (PW91) functional 88 for H/Pd(100) and N/Ag(111) and the Revised Perdew-Burke-Ernzerhof functional (RPBE) 89 for N 2 /Fe(110). The surfaces were modelled by five-layer (2 × 2) periodic cells for H/Pd(100) and N/Ag(111) and by a four-layer (3 × 3) periodic cell for N 2 /Fe(110). The Brillouin zone was sampled with n × n × 1 Monkhorst-Pack meshes 90 , with n = 6 for H/Pd(100), n = 5 for N/Ag(111) and n = 3 for N 2 /Fe(110). The energy cut-offs for plane wave basis sets were 350 eV for H/Pd(100), and 400 eV for the other systems.
In Fig. 1 we compare the spatial distributions of the three bare surface electron density models with n sur for  (100) surface. In each case, n sur is calculated using the surface atom positions of the relaxed adsorbate-surface structure. For completeness, the selfconsistent electron density of the whole system n SCF is also shown. The curves display the one dimensional (1D) cut of the electron densities along the line normal to the surface that contains the adsorbate. The n FS sur data practically coincide with n sur due to the negligible displacements that the adsorbed H causes on the surrounding Pd (the displacements of the adsorbates nearest neighbors d nn are shown in Table I). Obviously, dissimilarities appear in the region of large density gradients. The other two methods also reproduce well the values of n sur . The agreement is particularly good at the position of the adsorbate z A (dashed vertical line), which is the value of interest for applying the LDFA. It is also worthy to remark the goodness of the Hirshfeld partitioning scheme in removing the contribution of the adsorbate from the electron density of the whole system n SCF . This becomes apparent when noticing how different n SCF is from n sur , but how finely the latter is reproduced by n H sur . Deviations of n H sur and n AS sur from the correct n sur are more apparent at distances z < z A , i.e., in the region where the adsorbate-surface bond and the surface metal bonds are formed, especially for the top and bridge sites. On the one hand, n H sur reproduces rather well the surface density for the bridge site (notably, around z = 0, where the metallic character of the bonds is manifest), but it overestimates it for the top case. On the other hand, n AS sur behaves well for the top site, but it underestimates n sur around the topmost Pd layer (z = 0) in the bridge case. As argued below, these results are a consequence of the different charge redistribution occurring at each site.
The electron densities for N adsorbed at the top, hcp hollow, f cc hollow and bridge sites of the Ag(111) surface are compared in Fig. 2. As in the H/Pd(100) system, the three models are good in reproducing n sur at the adsorbate position and above (z ≥ z A ). The n H sur values are impressively good along the whole z-range, while both n FS sur and n AS sur deviate from n sur in the region around the Ag topmost layer. In particular, n FS sur and n AS sur are sys-  tematically overestimating and underestimating, respectively, the surface density. The somewhat large errors obtained with n FS sur when N is adsorbed on bridge are a consequence of the large displacements of about 0.18Å that N is causing around it (see Table I for details on the other sites).
As depicted in the insets of Fig. 3, N 2 can adsorb on the Fe(110) surface with upright orientation on top of the Fe atoms and parallel to the surface with its center of mass over either the hollow or the bridge site 49,65 . It is worthy to note that the latter adsorption configuration only appears on the relaxed Fe(110) surface 65 . Figure 3 shows for the three N 2 adsorption sites that any of the proposed surface density models succeeds in reproducing n sur at the position of each N atom conforming the molecule. Regarding the performance of the three models along z, the smallest to the largest errors are obtained by n FS sur , n H sur , and n AS sur , following this order. As noted previously, the errors are expectably larger in the regions where the density changes rapidly with z.
All in all, the analysis of the different cases studied in Figs. 1, 2, and 3 allows us to extract the following conclusions. The three surface density models provide a good description of n sur at distances from the surface close or larger than the equilibrium adsorption heights. In the bonding region between the adsorbate and the surface, the possible errors introduced by n H sur and n AS sur are reasonably small, while the adequacy of n FS sur depends strongly on the size of the lattice distortions, particularly, in those areas of large density gradients. By construction, n H sur is expected to overestimate (underestimate) the density whenever the adsorbate-surface interaction causes a negative (positive) induced density, i.e., a removal (piling up) of electrons, whereas n AS sur underestimates the electronic density in the interstitial region where the metal character of the surface atoms bonding is manifested.
As additional stringent tests, we have also compared the performance of the proposed density models for surface lattice distortions we encounter in real AIMD simulations as those presented in the next section. Figure 4 shows the results for one trajectory that is characterized by large lattice distortions, with averaged displacements with respect to the equilibrium positions that vary between 0.2Å and 0.4Å in the topmost surface layer. This trajectory corresponds to the adsorption of a 0.1 eV incident N on Ag(111). The self-consistent bare surface electron density n sur and the model densities n FS sur , n AS sur , and n H sur are calculated using the surface atom positions at different instants along the trajectory. The bottom panel shows the values of the densities at the position where the N atom is located at each instant during the simulation. Clearly, both n H sur and n AS sur in this order are the best approximations to n sur . In contrast, n FS sur while valid for small distortions, fails otherwise. The upper panels show a 1D cut of the surface densities along the same line used in Fig. 2 at two distinct instants. In both cases, n H sur is the model density that gives an overall better description of n sur .

III. PERFORMANCE OF THE SURFACE DENSITY MODELS IN AIMDEF SIMULATIONS
In this section we study the performance of the three density models in a gas-surface dynamics problem, namely, the adsorption and relaxation of hot species on metal surfaces. Although we have shown that the differences in the densities are small, it is not clear that they will be manifested also as small differences in the dynamical magnitudes for, at least, three reasons: (i) the surface atom displacements vary in magnitude and are in persistent change along the trajectory, resulting in configurations where the different models can provide a fluctuactingly faithful description of the bare surface density (as shown in Fig. 4, where the density can be overestimated as well as underestimated); (ii) the friction coefficient η is not linearly dependent on n sur ; and (iii) since the friction force is also proportional to the projectile velocity, the electron density alone gives incomplete information about the e-h pairs excitation. Therefore, a detailed dynamical analysis is revealed as a necessary complement to the static one.
In this respect, H/Pd(100), N/Ag(111), and N 2 /Fe(110) are well suited for the present analysis because they cover the limiting cases in which the energy exchange with the surface is dominated by either e-h pairs or phonons excitations 65 . For each system, we will examine how the differences in the densities originated by each density model affect: (i) the adsorption probability (Sec. III A), (ii) the surface atom displacements (Sec. III B) and the friction coefficients experienced by the hot species (Sec. III C), which are the factors determining the energy dissipation mechanisms, and, importantly, (iii) the kinetic energy loss of the hot species (Sec. III D), which is the central quantity of the problem.
The results that follow for each system are statistical averages obtained from a suitable number of AIMDEF trajectories, in which the two outer layers of the surface are allowed to move (unless otherwise stated). In the simulations, the Beeman predictor-corrector algorithm is used to integrate the classical equations of motion 91 , where 0.1, 0.5, and 0.7 fs are the time steps for H/Pd(100), N/Ag(111), and N 2 /Fe(110), respectively. In the following, the three sets of AIMDEF simulations carried out using n FS sur , n AS sur , and n H sur are correspondingly denoted as FSM, ASM, and HM. We note that the friction coefficient is in all cases neglected when the hot species move in very-low-density regions. Specifically, η = 0 for surface densities smaller than 7.46×10 −3 e/Å 3 (r s > 6 a.u.).
For each density model, 50 hot H atom trajectories are simulated on Pd(100) that result from the dissociation of 25 H 2 molecules on the surface, where they impinge at normal incidence with initial kinetic energy E i = 0.5 eV. As in previous works 63, 65 , the initial coordinates (x i , y i , z i ) and velocities of the individual H atoms are taken from adiabatic frozen surface molecular dynamics simulations on a six-dimensional PES 92 that describes H 2 dissociation on Pd(100). When the H-H distance in those simulations reaches three times that of the H 2 equilibrium bond length in the gas-phase, we set the time t = 0 for the present AIMDEF simulations.
The AIMDEF simulations with N atoms and N 2 molecules account, instead, for the complete adsorption process at normal incidence, where the used E i values ensure large adsorption probabilities. For N atoms on Ag(111), 20 trajectories are simulated for each density model with E i = 0.1 eV, z i = 4Å, and random (x i , y i ) values. The open-shell character of N requires the use of spin-polarized DFT. However, since the spin is quenched upon N-Ag interaction, computational effort can be saved by doing non-spin-polarized calculations when N lies close to the surface (see Appendix for a practical description of how this effect is considered in the AIMDEF simulations).
In the simulations of non-dissociative adsorption of N 2 on Fe(110), the molecules impinge normal to the surface with initial translational energy E i = 0.75 eV and zero rotational and vibrational energies (i.e., the zero point energy is neglected). The initial coordinates of the N 2 center of mass are, as in the previous case, z i = 4Å and random (x i , y i ). For each density model, 80 trajectories are calculated.

A. Adsorption probabilities
Previous MD calculations performed on a precalculated (frozen surface) three-dimensional N/Ag(111) PES show that the initial adsorption probability for N imping-ing at off-normal incidence with E i = 0.1 eV is S 0 0.98 for an ideal surface temperature T s = 0 K 48 . In those simulations, the effect of e-h pair excitations is described through the LDFA, while energy exchange with the surface lattice (phonon excitations) is included by means of the generalized Langevin oscillator (GLO) model 46,93-95 . The same value was obtained in pure GLO calculations that only included phonons excitations, while the authors found a slightly smaller value S 0 0.87 when only e-h pair excitations were considered.
In the present AIMDEF simulations, we obtain S 0 = 1 irrespective of the adopted surface density model (FSM, ASM, HM). For comparison purposes, we have performed two additional types of simulations: (i) AIMDEF calculations in which the surface atoms are fixed at their equilibrium positions (FS+EF) and (ii) usual AIMD calculations without the electronic friction force, in which the two outermost surface layers are allowed to move (NFS). Overall, our S 0 values are consistent with the previous MD results: on the one hand, NFS simulations yield S 0 = 1 and, on the other hand, FS+EF simulations yield a slightly lower S 0 = 0.85.
In the case of N 2 /Fe(110), the global adsorption probability is S 0 = 0.75 and the site-specific adsorption probabilities are 0.31, 0.13, and 0.31 for top, hollow, and bridge configurations, respectively. Here, too, the latter values remain barely unchanged when using any of the proposed surface electron density models. Also in this system we have performed the two additional types of simulations described above. NFS simulations yield S 0 = 0.71, which is in good agreement with GLO simulations carried out with a six-dimensional N 2 /Fe(110) PES for the same incidence conditions and low T s 49 . Interestingly, no adsorption event is observed with the FS+EF calculations that neglect energy exchange with the surface lattice.
In the following, we will focus on the adsorption process and restrict all the analysis to the results obtained from the adsorption trajectories exclusively.

B. Surface atoms displacements
The top and bottom panels of Fig. 5 show the mean displacements of the surface atoms within the first d 1 and second d 2 layers, respectively, as a function of time for each electron density model and for each system. The displacements are evaluated with respect to the equilibrium position for each trajectory and time step as, where l = 1, 2 indicates the topmost and second layers, respectively, r n are the surface atom positions, and the sum runs over the N l atoms within the l-th layer of the unit cell. A common trend in the three systems is that d 1 > d 2 . This is a reasonable result, since the projectile exchanges energy and momentum directly with the For H/Pd(100) and N 2 /Fe(110) the differences between the displacements calculated with the three surface density models are almost negligible. For N/Ag(111) small differences appear after 1.3 ps. It is only after this time that the displacements of the HM model are visually lower than the ones of the FSM and ASM models.
Comparing the three systems, the Ag(111) displacements are much larger and increase faster with time than those of Pd(100) and Fe(110), which reach a plateau at earlier times. Notice that the values for Fe arrive close to their maximum less than 1 ps after the moleculesurface collision, while Ag movements continue to increase in amplitude even 3 ps after the collision. Interestingly, this different behaviour is not correlated with the different projectile-to-surface atom mass ratios γ we have for H/Pd(100) (γ=0.0095), N/Ag(111) (γ=0. 13), and N 2 /Fe(110) (γ=0.5). Within simple binary collision models 96,97 , this parameter relates (albeit not exclusively) to the projectile-to-surface momentum transfer in the successive collisions with the metal atoms. The low γ value of H/Pd(100) is in line with the small Pd displacements, but this argument alone cannot explain the N/Ag(111) and N 2 /Fe(110) results of Fig. 5. These results are neither explained by the E i value, which is smaller for N than for N 2 . Instead, Fig. 5 is to be understood by considering the different PES topographies of these systems. The N 2 /Fe(110) system has an entrance energy barrier above the surface, which diminishes the kinetic energy of the impinging molecules 49 , whereas the N/Ag(111) PES is barrierless and strongly attractive at long range, which accelerates the N atoms 30 .

C. Friction coefficients
In this subsection, we focus in the friction coefficients η experienced by the hot atoms during the simulations. Figure 6 compares their statistical averages η as a function of time for the three systems and for the different surface electron density models. In the case of N 2 /Fe(110), we find meaningful to show the results for each of the adsorption configurations because the probed surface regions and hence the friction coefficients experienced by the molecules depend strongly on whether N 2 is adsorbed on the top-vertical or on the bridge-and hollow-parallel wells.
In accordance with the small surface atoms displacements occurring in H/Pd(100), the η values in the three density models are very similar. However, for N and N 2 adsorbed on any of the parallel configurations, Fig. 6 shows that FSM yields an overall overestimation of η . This behaviour becomes very clear at times t > 800 fs in N/Ag(111) and t > 500 fs in N 2 /Fe(110), namely when the energy lost into the phonons channel starts to saturate 65 . Comparing the time-averaged values of the friction coefficients η t computed for t > 300 fs, the relative differences between the FSM and HM simulations, calculated as ∆ = ( η FSM t − η HM t )/ η HM t amounts to 25% for N/Ag(111) and 16-20% for parallel N 2 on Fe(110). The FSM overestimation is a consequence of using the undistorted bare surface density n FS sur for moving surface atoms in cases in which the adsorption dynamics is dominated by on-surface (z < 2Å) movements. Along those trajectories, the large density regions existing within very short distances of the surface atoms are accordingly very repulsive and, hence, inaccessible for the typical hot species energies. However, if the surface atoms move from their equilibrium positions, the hot species may access those, otherwise, forbidden regions, where the n FS sur values are large because they correspond to the undistorted surface. Since the electron density gradient increases rapidly as the distance to the surface atomic cores decreases, it is understandable that the time spent by the hot atoms nearby these regions, though short, will have more weight in the statistical average and thus result in the overall η overestimation observed in the FSM curves of Fig. 6. On the contrary, with the HM and ASM models, the probed densities, and consequently the friction coefficients that enter the average, are always similar to the actual distorted surface density values, as shown in Fig. 4.
The performance of the FSM model for N 2 adsorbed on the top-vertical configuration is completely different. Figure 6 shows that FSM largely underestimates η as compared to ASM and HM. The discrepancies start at t > 250 fs, while the relative difference in the timeaveraged η between the FSM and HM simulations is around −37%. In this case, the molecule is mainly moving along the surface normal at 2-3Å above the surface in a concerted N 2 -Fe motion that brings the Fe atom inwards and also outwards the topmost layer. In contrast to the parallel-N 2 adsorption dynamics, the large density-gradients along the surface normal appear in the low-density regions of the undistorted surface that are probed by the top-vertical N 2 during the outwards motion. Therefore, the same large-density-gradient argument explains that during the concerted N 2 -Fe movement n FS sur is now predominantly underestimating the density. When FSM is used, a modulation in η is clearly visible for N/Ag(111) and N 2 /Fe(110), which consists in large-amplitude low-frequency oscillations with periods ∼ 0.8 and ∼ 0.3 ps, respectively. The trajectories that enter the statistical averages are not correlated and thus this modulation is to be interpreted as a mere statistical artifact. As a matter of fact, it is observed that the projectiles impact on different positions within the surface unit cell and that the paths followed by the hot species on the surface are very different. Nonetheless, there are cases in which an overlying low-frequency modulation that we tend to ascribe to the surface atoms movement seems to be also present and could explain the modulation in the FSM η . A considerably much larger statistics would be needed to confirm that surface phonons and not statistical errors are at the origin of these FSM oscillations. However, this point is completely out of the scope of the present study.

D. Energy loss
In the previous subsections we have demonstrated that the use of different models to evaluate the bare surface electron density during AIMDEF simulations of hot species on surfaces results in substantially different mean friction coefficients and, in some cases, also lattice distortions. These quantities determine the energy loss rate of the hot species, which after all, is the key quantity in the modeling of reactive processes on surfaces. The hot species kinetic energy is directly linked to several experimentally observable magnitudes, such as the maximum distance traveled on the surface, the relaxation process time-scale, and the amount of energy transferred to the substrate. Last but not least, there is the question of how this energy is partitioned into e-h pair and phonon excitation contributions, and to what extent the accuracy in the description of the bare surface charge density influences the partitioning. Since the friction force is also velocity-dependent, the variation of the relaxation rate with the density model cannot be predicted. Figures 7, 8 and 9 show, for the three systems and for the different density models, the kinetic energy of the hot species averaged over the adsorbed trajectories as a function of time, denoted E A kin for the atom and E M kin for the molecule in the following. Importantly, the general observation is that this quantity is not sensitive to adopting FSM, HM or ASM to describe the electron densities (there may be subtle differences that are, nonetheless, of similar magnitude as the oscillations in the curves). Considering that for N/Ag(111) and N 2 /Fe(110) the η values obtained with FSM deviate from those obtained with ASM and HM (see Fig. 6), it is unexpected to find hardly any difference between the corresponding E A,M kin curves (see Figs. 8 and 9). This is indeed a remarkable observation, since it stresses that uncertainties in the friction coefficients do not necessarily translate into the final dynamics and the measurable magnitudes of interest.
The PES corrugation. In the trajectories that enter in the averages, the hot species follow disparate routes on the surface, sampling thus PES regions of very different energies. Contrary to the rapidly oscillating behaviour of E A,M kin , the average energy dissipated into e-h pair excitations E eh (t) shows a smooth time dependence (see upper panels of Figs. 7, 8, and 9). Therefore, the observation of E eh (t) allows a sound comparison of the performance of the three density models on the relaxation rates. For each individual trajectory, this energy is evaluated as where the summation runs over the atoms that constitute the hot species, and where v n (t ) is the atom instantaneous velocity. At the end of the simulation time, when the H atoms are close to being thermalised on Pd(100), they have lost 0.53 eV into e-h pair excitations and the energy differences between models are minimal, below 0.01 eV. For N 2 /Fe(110), where we separately consider the energy loss for each adsorption site (see Fig. 9 upper panels), the differences between FSM and the other two models are more noticeable, and they are also in line with the underestimation or overestimation behaviors expected from the η values. Moreover, it must be taken into account that the molecule is far from being relaxed, and thus we can anticipate that the incipient deviations observed in Fig. 9 top panels will grow at longer times. This effect is manifested in N/Ag(111) too, where a very clear monotonously increasing deviation of FSM with respect to HM exists at the end of the simulation time.
Here, the amount of energy used to excite e-h pairs is 0.499 eV with HM, and 0.054 (0.018) eV more (less) than that with FSM (ASM). Again, the larger E eh values provided by FSM in the latter system are consistent with the η overestimation shown in Fig. 6. The general conclusion we extract from the behaviors of E A,M kin and E eh is that the calculated energy loss rates and the relaxation times are, for practical purposes, density-model-independent. In other words, the three models under scrutiny are able to provide similar descriptions of the hot species adsorption dynamics of diversely behaving systems. Nevertheless, among the studied models, HM is the one that provides the best description of the bare surface electron density, and therefore its use should be recommended in simulations when there is no a priori knowledge of the dependence of the energy loss on the friction coefficient values and on the surface atoms motion. HM overcomes the limitations of FSM and ASM to describe, respectively, the instantaneous density when the surface atoms are free to move and the bonds between surface atoms.
For completeness, additional curves are shown in Figs. 7, 8, and 9 that correspond to the two types of calculations described in Sec. III A, namely, NFS simulations without electronic friction and FS+EF simulations with electronic friction but without surface atoms movement. Their comparison to the AIMDEF results provide a reference picture of the e-h pairs excitation importance during adsorbate relaxation in these systems. The conclusion obtained is that e-h pairs excitation is the dominant dissipation channel for the hot H atoms, while phonons dominate N relaxation. Actually, as mentioned earlier, no N 2 adsorption is found when only the electronic energy dissipation channel is considered for N 2 /Fe(110). A detailed explanation of these results can be found in Refs. 65 and 66.
One remaining issue when accounting for the energy dissipation mechanisms will be to incorporate the e-h pair excitations created by the kinetically excited surface atoms. This effect can easily be incorporated with our present AIMDEF method. However, for the systems and incidence conditions considered here, we find that the latter can be neglected within the time scale of the considered adsorption processes. We have made some estimations for the case of N/Ag(111), which is the system where the substrate atoms acquire the largest kinetic energies upon interaction with the impinging atom. The electron density for a surface Ag atom can be estimated by taking the electron density at a surface vacancy position. The corresponding mean electron radius would be r s ∼ 5 a.u., which results in a friction coefficient η = 0.044 a.u. for Ag. Using this value to estimate the energy relaxation rate for an ideal damped oscillator, we obtain η/m(Ag) = 0.009 ps −1 . Interestingly, taking the value of the electron-phonon coupling factor for Ag at 300 K 2.5 × 10 16 W/(m 3 K) 98 and dividing it by the specific heat of silver at the same temperature 2.52 × 10 6 J/(m 3 K) , one obtains a similar relaxation rate of 0.010 ps −1 . Both estimations show that the rate at which the mobile Ag atoms dissipate energy to the electronic system is two orders of magnitude slower than the time scales of the hot atom relaxation processes studied in the present work.

IV. CONCLUSIONS
We have examined the performance of different models of the bare surface electron density n sur in AIMDEF simulations of the adsorption dynamics of atoms and molecules on metals, using H, N, and N 2 on Pd(100), Ag(111), and Fe(110), respectively, as case studies. In the original formulation of AIMDEF the surface electron density n sur , which is used to calculate the electronic friction force acting on the adsorbing species, was approximated by that of the frozen surface n FS sur (FSM model). Here, we improve the methodology by using models that account for the n sur changes brought by the displacements of the surface atoms during the simulations, which can reach considerably large values of up to ∼ 0.4Å in some of the studied surfaces. The proposed n sur models are constructed on-the-fly at each simulation time step from either superposition of atomic electron densities (ASM model) or a Hirshfeld partitioning scheme (HM model) of the total self-consistent density.
From static analyses for a few fixed geometries, we deduce that all the models accurately reproduce n sur at the hot atom positions, as required by the simulations, and also that they provide good estimates at other positions. In a subsequent dynamical analysis, we find that the three of them yield similar energy loss rates, despite the limitations of FSM to model the distorted surface density as compared to ASM or HM. An in-depth exam-ination of the trajectories reveals that FSM can produce significant deviations in the friction coefficients that depend closely on the surface density regions visited by the adsorbates.
Although the results presented in this work apply to a particular class of dynamical processes, they allow us to establish some guidelines for the applicability of each model in a broader context. First, we have shown that, when the dynamics involves large displacements of the surface atoms, n FS sur clearly deviates from the average electron densities experienced by the hot species. Therefore, when modeling surface processes of similar characteristics, such as temperature effects in gas-surface interactions, ASM and HM will prove more reliable. Secondly, if the electronic structure of the surface under study is sensitive to changes in the interatomic distances, then HM is to be preferred over ASM, because it accounts more realistically for the charge distribution at the crystal bonds. This is particularly relevant in the description, for example, of surface penetration dynamics, where the projectile travels across both surface and bulk environments, of distinct electronic structure. Such penetration processes are more likely to occur, for instance, for faster impinging atoms. In practical terms, the evaluation of the Hirshfeld partitioning of charge in HM simulations does not imply a major computational cost increase with respect to the other methods. For all the reasons stated here, we can conclude that the use of a surface electron density model based in a Hirshfeld partitioning scheme can be considered a highly accurate and efficient strategy to describe e-h pair excitations in AIMDEF simulations. Finally, note that our new AIMDEF methodology is also well-suited to incorporate, when necessary, the effect of the e-h pair excitations created by the moving surface atoms.

ACKNOWLEDGMENTS
This work has been supported in part by the Basque Departamento de Educación, Universidades e Investigación, the University of the Basque Country UPV/EHU (Grant No. IT-756-13) and the Spanish Ministerio de Economía y Competitividad (Grant No. FIS2013-48286-C2-2-P). The authors thankfully acknowledge the computer resources, technical expertise and assistance provided by the Red Española de Supercomputación and by the DIPC computing center.

Appendix: Molecular dynamics zones
In usual AIMD (also AIMDEF) simulations, the converged wave functions at each integration step t 0 are used to extrapolate the wave functions at the next integration step. This scheme, which certainly facilitates the AIMD calculation, might be problematic when dealing with open-shell atoms or molecules for which the spin state changes with their distance to the surface, more specifically, when it changes from zero to a finite value because the wave functions used in the extrapolation are spin-degenerated. Such is the case of N incident on Ag(111). Figure 10 shows the spin magnetization of the system as a function of the distance from the surface z, for N located above a Ag surface atom. It is zero close to the surface and increases rapidly to the gasphase value of 3 µ B far from the surface. Similar variations are obtained for N located at other positions over the Ag(111) surface. The main numerical difficulty consists in converging to the correct non-zero spin-polarized ground state as the N-Ag(111) distance grows. Thus, in order to break the spin-degeneracy of the wave functions during the AIMDEF simulations, we define the following three different zones within the supercell (see inset of Fig. 10) and adopt a different strategy within each of them, • zone 1 (z < 1.2Å), where the N spin is completely quenched. Therefore, as soon as N enters and stays in this zone a standard non-spin polarized AIMD calculation is performed.
In this zone, the mentioned numerical problems in breaking the spin-degeneracy may appear when the N atom enters this zone from the surface (zone 1) with zero spin magnetic moment. In this case, the simulation is stopped at each integration step once the electronic wave functions and the forces on the atoms are converged. For the next integration step, a new calculation is launched from scratch, but using an initial magnetic moment of 3 µ B for N. There is no need to stop the simulation when N comes from the region z > 2.6Å.
• zone 3 (z > 2.6Å), where the spin magnetic moment is that of the gas-phase N. In this zone, a standard spin-polarized AIMD calculation is performed.