Energy loss and surface temperature effects in ab initio molecular dynamics simulations: N adsorption on Ag(111) as a case study

Dino Novko,1,* Ivor Lončarić,2,† María Blanco-Rey,1,3,‡ J. Iñaki Juaristi,1,2,3,§ and Maite Alducin1,2,‖ 1Donostia International Physics Center DIPC, Paseo Manuel de Lardizabal 4, 20018 Donostia-San Sebastián, Spain 2Centro de Física de Materiales CFM/MPC (CSIC-UPV/EHU), Paseo Manuel de Lardizabal 5, 20018 Donostia-San Sebastián, Spain 3Departamento de Física de Materiales, Facultad de Químicas, Universidad del País Vasco (UPV/EHU), Apartado 1072, 20080 Donostia-San Sebastián, Spain (Received 7 June 2017; published 28 August 2017)

We study surface temperature effects on the adsorption and relaxation of N atoms on Ag(111). To this aim, we perform ab initio molecular dynamics simulations with electronic friction, in which the surface is coupled to a thermostat that fixes the desired surface temperature. Simulations performed at 80 and 700 K show that the surface temperature has minor effects on magnitudes such as the initial adsorption probability, the relaxation rate of the adsorbing N, and the energy lost in electronic excitations. Slight differences are identified in the adsorption paths with the appearance of subsurface absorption events at 700 K that are not observed at 80 K. Furthermore, we perform additional simulations without a thermostat in order to examine the validity of commonly used ab initio molecular dynamics simulations in which no heat dissipation from the simulation cell is allowed. Our results show that such a methodology may not suffice to simulate the low-temperature regime since the surface becomes unphysically heated within a few picoseconds upon adsorption of the N atom. However, neither in this unfavorable case are the magnitudes defining the dynamics of the adsorbates at the same time scale significantly modified from those obtained at constant surface temperature. DOI: 10.1103/PhysRevB.96.085437

I. INTRODUCTION
Ab initio molecular dynamics (AIMD) [1,2], based on density functional theory (DFT), and the recently developed ab initio molecular dynamics with electronic friction (AIMDEF) [3][4][5] that also incorporates low-energy electronic excitations are advanced and powerful methods for the theoretical treatment of gas-surface dynamics. AIMD and particularly AIMDEF provide a joint and accurate description of the adiabatic forces and the relevant energy dissipation/exchange channels, namely, phonons and (low-energy) electron-hole (e-h) pairs. Still, surface temperature, which can be a determinant factor in the long-term gas-surface reactivity, is usually overviewed in most AIMD(EF) studies. Usually, the reason is that these studies are aimed to understand the evolution of elementary gas-surface processes, such as dissociative and nondissociative adsorption, scattering, and atom-atom abstraction, among others, in the short-time femtosecond and picosecond scale for which surface temperature effects are expected to be less important.
Even if the outcome of most of such processes is indeed decided within this short time scale, this argument alone may be questionable for the common room and higher surface temperatures (T s 300 K) under which most molecular beam experiments aimed to explore the basic properties of the gas-surface dynamics and reactivity are conducted. Properties such as the energy loss of the scattered beam, the rovibrational state of the scattered molecules, the thermalization rate and paths of the adsorbing species, and even the probability of the elementary processes, are presumably affected by T s [6][7][8]. Therefore, one in principle should be aware of possible surface temperature effects when comparing theory and experiments, particularly for T s 300 K.
Molecular beam experiments are typically performed in ultrahigh vacuum chambers keeping the surface at a constant temperature, while the incidence energy of the gas species (molecules/atoms) leaving the nozzles toward the surface ranges from tenths of meV to few eVs [9][10][11][12]. These conditions have been simulated with the generalized Langevin oscillator model (GLO) [13][14][15] in molecular dynamics calculations based on a precalculated adiabatic potential energy surface (PES) that usually neglects the degrees of freedom of surface atoms. Surface temperature effects and energy exchange with the surface lattice have been investigated using this methodology in different elementary processes and systems: adsorption [15][16][17][18][19][20][21][22][23][24][25], scattering [26][27][28][29], Eley-Rideal and hot atom abstraction [30][31][32][33][34]. However, the GLO simulations obviously miss effects caused by the surface lattice distortions that, in contrast, are naturally included in AIMD(EF) [35]. In this respect, one procedure to mimic the experimental constant temperature conditions while using constant-energy (microcanonical) AIMD simulations is the one followed in Refs. [6,[35][36][37][38]. Briefly, a first calculation is performed in which only the surface is thermalized to the desired temperature. Afterwards, customary microcanonical AIMD simulations of the gas-surface system are performed using as initial conditions for the gas species those of the experimental molecular beam and as initial conditions of the (moving) surface atoms the positions and velocities of the prethermalized surface. If the kinetic energy of the molecular beam E i is small compared to the sum of the kinetic energy of the surface atoms, this approach is expected to be accurate. However, in the case of a long-lasting gas-surface interaction process (e.g., adsorption), if E i is comparable to the total kinetic energy of the surface, the system will end with a temperature much higher than the desired initial surface temperature because there is no heat dissipation. This is not an uncommon situation in standard AIMD(EF) studies as the surface is typically represented by only a few tens of (mobile) atoms. At the common room temperatures, the total kinetic energy of the model surface would be hence in the range of a few eVs, i.e., comparable to the total (kinetic and adsorption) energy of the atoms and molecules of interest in surface reactions.
There are very few studies that have recently investigated the possible limitations of microcanonical AIMD simulations associated to the inherent use of a finite-size cell when studying the adsorption and scattering of gas species on surfaces. The QM/Me model developed by Meyer et al. [39] suggests that such size effects can be important in the adsorption and subsequent relaxation of O 2 on Pd(100). In contrast, various adsorption and scattering properties of O 2 on Pt(111) predicted from microcanonical and constant surface temperature AIMD simulations appear to be similar in the short-time scale [40].
Within this context, the aim of this work is twofold. On one side, we explore the size-cell limitations of customary AIMD and also AIMDEF simulations using an extreme example of large energy exchange between the gas species and the surface. Specifically, we investigate the adsorption on Ag(111) of 0.1-eV N atoms that, being rapidly accelerated by the largely attractive atom-surface interaction, cause an important heating of the surface [4,5,18,41]. Thus, we find it meaningful to test the reliability of the AIMD and AIMDEF methods in describing the gas-surface dynamics under such extreme conditions. As a second objective, we also investigate the dependence on T s of different properties of the adsorption process itself and of the energy dissipation channels (phonons and e-h pairs). To this aim, we have performed standard AIMD and AIMDEF simulations as well as AIMD and AIMDEF simulations with a thermostat coupled to the surface atoms that assures the desired constant T s . Among the various thermostats used in molecular dynamics simulations (see Ref. [42], for instance), here we have adapted the Nosé-Hoover thermostat [43,44] such that only the temperature of the surface atoms is fixed when performing either AIMD or AIMDEF calculations. The practical advantages of using this thermostat are that it provides relatively realistic trajectories, its dynamics is time reversible, and it conserves the energy of the extended system, which permits to monitor the accuracy of the numerical integration.
The paper is organized as follows. Section II describes the theoretical method and the practical implementation of the Nosé-Hoover thermostat that has been adapted to be applicable in general gas-surface dynamics simulations based on AIMD and AIMDEF. The methodology is applied in Sec. III to study surface temperature effects on the adsorption of N on Ag(111). The applicability of standard AIMD and AIMDEF in this system is also discussed in this section. A summary of the results and the main conclusions of this work are provided in Sec. IV.

II. THEORETICAL METHODS
In the present general formulation of gas-surface dynamics, the system formed by all allowed moving atoms is partitioned in two or three coupled subsystems to fulfill the most general experimental conditions in which the surface is kept at a constant temperature T s and the gas species, instead of being thermalized with the surface, has a distinct kinetic energy E i . Each subsystem is defined such that its corresponding equations of motion account for the main factors contributing to the dynamics of this particular subsystem as detailed next.
The first subsystem consists of the gas species that, in general, is not initially thermalized with the surface. The gas species dynamics is assumed to be ruled by the adiabatic multidimensional potential energy surface. Furthermore, the effect of the low-energy electronic excitations and deexcitations can eventually be introduced in terms of additional frictional and stochastic forces. Thus, the equation of motion for each atom i in the first subsystem is where m i , r i , and η i are its corresponding mass, position vector, and electronic friction coefficient. The first term in the right-hand side of this equation is the adiabatic force that depends on the position of each atom in the gas-surface system. This force is calculated in usual AIMD/AIMDEF simulations within DFT, applying the Hellmann-Feynman theorem. The friction force (second term) is calculated following the AIMDEF scheme of Refs. [4,5,41], which is based on the local density friction approximation (LDFA) [45]. The latter approximates the electronic friction coefficient η i (r i ) as the one that this very same atom would have when moving within a free electron gas of electronic density equal to the electron density of the bare surface at the atom position n sur (r i ). Among the different models proposed in Ref. [5] to evaluate on-the-fly n sur (r i ), the present AIMDEF simulations are performed following the one based on the Hirshfeld partitioning scheme [46]. The random force R[T e ,η i (r i )] (third term) describes the effect of e-h pair excitations and deexcitations for nonzero electronic temperatures T e and it is commonly approximated by a Gaussian white noise of variance Var[R(T e ,η i )] = (2k B T e η i )/ t, where k B is the Boltzmann constant and t is the integration time step. The R(T e ,η i ) term is crucial when describing desorption processes induced by femtosecond laser pulses [47][48][49][50][51][52][53][54][55], but it can safely be neglected in simulations of common molecular beam experiments, which are usually performed at surface temperatures of a few hundred Kelvin for which the e-h pair distribution hardly deviates from its ground state at 0 K. Regarding the surface atoms, the constant surface temperature conditions are simulated by thermalizing all or part of the mobile surface atoms (see below) with the Nosé-Hoover thermostat [43,44]. In the Nosé-Hoover thermostat, the real system is coupled to a dimensionless dynamical variable s that together with its conjugate momentum p s act as a thermal bath, i.e., ensure a canonical distribution in the real system. Thus, the equations of motion of these surface atoms j with mass m j and position vector r j that constitute the second subsystem in the present formulation are where the thermodynamic friction coefficient ξ is related to the dynamical variable and momentum by ξ = Q −1 sp s [42]. In these expressions, Q is a parameter with dimensions of energy×time 2 that acts as the mass of the dynamical variable s and N 2 is the number of atoms in the second subsystem. It can be proven that any positive and finite Q leads to a canonical distribution (for Q → ∞ a microcanonical distribution is recovered). However, Q should be chosen rather carefully as too large values of Q ensure a canonical distribution only after very long simulation times, while too small Q values lead to large high-frequency oscillations of the instantaneous temperature. In a canonical ensemble, the fluctuations of the system total energy E around its average value are related to the system heat capacity C V and temperature T through the The latter together with the frequency of the instantaneous-temperature oscillations can guide the selection of a reasonable Q value as specified in Sec. III A. Even though the Nosé-Hoover thermostat gives relatively realistic trajectories of the canonical system (i.e., the dynamics is not as disturbed as in other popular thermostats) [42], it sounds reasonable to consider that the temperature of the atoms that directly suffers the possibly strong impacts with the gas species must considerably deviate from the macroscopic T s during a time scale comparable to the thermalization of the gas species. Thus, following the spirit of the primary and secondary regions defined in the generalized Langevin oscillator model [13,14], one may also introduce in the present formulation a third subsystem that consists of the surface atoms that directly interact with the gas species. The dynamics of each atom k in this third subsystem is simply described by the classical Newton equations and the adiabatic approximation where m k and r k are its corresponding mass and position vector. As done in Eqs. (1) and (2), the adiabatic force is calculated with DFT and the Hellmann-Feynman theorem. This subdivision of the simulation cell in three subsystems was already performed in Ref. [40]. Next, we apply this formulation to study temperature effects in the adsorption dynamics of N on Ag(111).

III. ADSORPTION AND RELAXATION OF N ON Ag(111)
In this section we use the adsorption dynamics of N on the Ag(111) surface as a case study to both analyze methodological and physical aspects relevant to present gassurface simulations. Previous calculations performed on this system for different hyperthermal incidence energies showed already that the adsorption process is initially dominated by the energy exchange with the Ag atoms [4,5,18,41]. Actually, such an energy exchange may cause significant surface lattice distortions according to the AIMDEF simulations of Refs. [4,5,41]. The latter makes the system particularly appealing to investigate the possible limitations of usual microcanonical AIMD simulations in describing the energy flow from the atom-surface collision region into the bulk using the standard supercell approach. Furthermore, the theoretical adsorption energy of N on either the hollow fcc (2.26 eV) or hcp (2.38 eV) sites [41] is large enough to study surface temperature effects using quite extreme, though reliable, T s values.

A. Computational details
The DFT calculations presented in this work are done with the VASP package [56,57], which uses a plane-wave basis set. The energy cutoff for the latter is set to 400 eV. The coreelectron interaction is approximated by projector augmentedwave (PAW) potentials [58]. In order to be consistent with the previous work [4,5,18,41], the PW91 exchange and correlation functional is used [59]. The surface is modeled by a five-layer (2 × 2) periodic slab in the supercell approach. The Brillouin zone is sampled with a 5 × 5 × 1 Monkhorst-Pack mesh [60].
All the AIMD and AIMDEF results presented below are obtained from a statistical average of 20 trajectories in which the impinging N atom starts at z i = 4Å from the surface with a vertical kinetic energy of E i = 0.1 eV and random (x i ,y i ) values [4,5,41]. 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 in Ref. [5] for more details). In the simulations, the Beeman predictor-corrector algorithm is used to integrate all the classical equations of motion [i.e., Eqs. (1)-(4)], with the time step of 0.5 fs. The simulations are run for 3.5 ps, which is sufficient for the purpose of this study because the N atoms transfer most of their energy to the surface during this period of time (see below).
In order to comprehend the effects of using a constant surface temperature, we perform two well-differentiated types of AIMD and AIMDEF calculations, namely, the customary AIMDEF and microcanonical AIMD description as in Refs. [4,5,41], in which the surface atoms are initially at rest in their equilibrium positions (i.e., with an initial surface temperature T s0 = 0 K) and only the two topmost surface layers are allowed to move, and another with constant surface temperature T s using the Nosé-Hoover thermostat as described in Sec. II. In this latter type of calculation, we distinguish the low-and high-temperature regimes with T s = 80 and 700 K, respectively. For simplicity, the electronic temperature in the corresponding AIMDEF simulations is, however, kept at 0 K in all cases, as argued below.
The procedure for the constant-temperature calculations goes as follows. First, the three topmost layers of the bare Ag(111) surface are thermalized to a desired constant temperature by using Eqs. (2) and (3) in AIMD simulations of the bare Ag(111) surface. As already pointed out in Sec. II, the parameter Q that enters these equations needs to be chosen with great care, and one reasonable option would be to choose a Q value that satisfies the relation between the system total-energy fluctuation and the corresponding heat capacity, i.e., E 2 − E 2 = k B T 2 s C V . By making a simple approximation that surface atoms constitute a collection of N threedimensional harmonic oscillators, we can cast the total-energy fluctuations into (instantaneous) temperature fluctuations, so that the above expression now reads as (3N ). However, in our case the frequencies of the instantaneous-temperature fluctuations obtained from the latter equation deviate strongly from the usual phonon frequencies of the Ag(111) surface [61]. Therefore, to make the T frequencies closer to those of the Ag(111) phonons, but still in connection to C V , we allow for twice as large temperature fluctuations, i.e., (3N ). Now, the Q values that satisfy this equation for T s = 80 K (C V = 17 Jmol −1 K −1 ) and T s = 700 K (C V = 25 Jmol −1 K −1 ) [62] turn out to be 0.53 × 10 −29 a.u. (atomic units) and 0.53 × 10 −28 a.u., respectively. After these values are set, we perform thermalization simulations up to 5 ps for each of the surface temperatures. Following a 2.5-ps equilibration period, we pick 20 uncorrelated surface configurations from the remaining time period and we use them as the initial surface condition for the constant-T s AIMD and AIMDEF simulations of the N/Ag(111) system. Finally, we divide the N/Ag(111) system as explained in Sec. II: the N atom dynamics is described by means of Eq. (1) (first subsystem), the third topmost Ag(111) layer by means of Eqs. (2) and (3) (second subsystem), and the first two topmost layers by means of Eq. (4) (third subsystem). Recall that the separation of the mobile surface atoms into the second and third subsystem is particularly recommended in those cases for which a large gas-surface energy exchange is expected since the thermostat would otherwise cause an unrealistically fast dissipation of the excess energy from the surface atoms in close contact with the adsorbing species and from the adsorbing N. This is what we observe in Fig. 1 by comparing the constant-T s AIMD results in which the two moving surface layers are coupled to the thermostat (blue curves) with those obtained from a two-subsystem description of the surface (red and green curves). However, assuming that the surface is already separated into two subsystems, the specific definition of each of them seems to be less relevant. Figure 1 shows that the mean kinetic energies of the adsorbing N and of the Ag atoms in the topmost layer as obtained from constant-T s AIMD simulations in which either the second or the third layer is coupled to the thermostat are quite similar.
In addition to the aforementioned two main types of calculations, we perform a few sets of selected calculations to explore other related issues. First, the possible effect of the electronic temperature is analyzed by constant temperature AIMDEF simulations in the (extreme) high-temperature regime T s = T e = 700 K, using 10 trajectories. No significant changes in the N adsorption dynamics and relaxation rates are found, which justifies the use of T e = 0 K in the rest of AIMDEF simulations that are discussed in Sec. III B. Finally, the limitations of the simulations without thermostat for less unfavorable conditions that would correspond to high T s are examined with microcanonical AIMD simulations using three mobile surface layers with an initial surface temperature T s0 = 700 K.

B. Results from AIMDEF simulations: Surface temperature effects
The evolution of the surface kinetic energy during the N adsorption process is the quantity that presumably will be most affected by the initial T s and by the use or not of thermostats assuring that the temperature of the deeper surface region is kept at the desired T s . The results from the AIMDEF and AIMD simulations for the surface kinetic energy are shown in Figs Fig. 2 for AIMD simulations without electronic friction: Results obtained from simulations without thermostat at initial surface temperatures T s0 = 0 K (black curves) and 700 K (blue curves) and from simulations with thermostats at T s = 80 K (red curves) and 700 K (green curves).
was identified as the main energy dissipation mechanism during the initial stages of the adsorption process [4,5,18,41]. Thus, the very similar results obtained when the effect of e-h pair excitations is (Fig. 2, AIMDEF results) or not (Fig. 3, AIMD results) included are not surprising. In both cases, the third layer, which is the one coupled to the thermostat, exhibits the expected quasioscillatory behavior around the value E (3 rd L) kin = 3N (3 rd L) /2 (k B T s ) associated to the nominal T s (i.e., 0.04 eV at T s = 80 K and 0.36 eV at 700 K, with N (3 rd L) = 4 in our simulations). Obviously, the first layer, being directly involved in the collision with the incoming N, is the one from the two freely moving layers that experiences the largest E kin variations. Furthermore, as expected from classical dynamics, the kinetic energy gain during the very first collisions (t 0.5 ps) is larger when the surface is at 80 K rather than at 700 K and, hence, the cooling down of the two topmost layers also occurs at a faster rate for T s = 80 K. 1 Focusing now on the standard AIMDEF and AIMD simulations performed without thermostat for an initial surface temperature T s0 = 0 K (black curves), we observe that the energy exchange with the moving layers is correctly described at the beginning of the adsorption process since E (1 st L) kin and E (2 nd L) kin are similar to the results obtained with the thermostat at T s = 80 K. However, after about 1.2 and 1.8 ps, the kinetic energies of the second and first layers, respectively, continue to rise and are about to reach the corresponding curves at 700 K at the end of the simulations (3.5 ps). This faulty behavior appearing at long integration times is precisely caused by the use of a finite supercell, i.e., a limited number of atoms, to represent the surface. In absence of a thermostat allowing for energy dissipation, the few surface atoms used in the simulations will accumulate most of the energy exchanged with the adsorbing atom, what may give rise to unrealistic surface dynamics in the long term. However, it is worthy to note that this error is expected to diminish as the number of mobile surface atoms increases and the gas-surface energy exchange decreases. This is precisely what our additional results from microcanonical AIMD calculations performed with the surface atoms initially at 700 K and three mobile layers show in Fig. 3. Despite the fact that the mean kinetic energy of the third layer separates from the constant-T s AIMD results at T s = 700 K, both types of simulations provide a rather similar description of the first-and second-layer dynamics. Recall also that even for the less favorable conditions of low T s , one important point is to determine if the finite-cell problem detected already in the surface atoms after the first 1-2 ps affects the N relaxation dynamics and particularly any of the magnitudes of interest, such as the adsorption probability as well as the mean relaxation rate, time, and paths of the adsorbing N atoms.
Regarding the initial adsorption probability, 2 trajectories out of 20 are scattered for AIMDEF at T s = 700 K, while only one for AIMD at T s = 700 K and AIMDEF at T s = 80 K. In all the other cases, the results at the end of the simulation (3.5 ps) point to 100% adsorption. These values suggest a decrease in the number of adsorption events when increasing T s . However, the number of trajectories is statistically insufficient to be conclusive in regards to the temperature dependence of the initial adsorption probability. Information on the relaxation rate and its dependence on T s can be extracted by analyzing the time evolution of the mean kinetic energy E (N) kin of the adsorbing N. Figure 4(a) shows that the values and dissipation rate of E (N) kin obtained from the different AIMDEF simulations are in all cases very similar and, importantly, when comparing the AIMDEF results without and with the 085437-5 thermostat at T s = 80 K. Similar values of E (N) kin (not shown) are also obtained from the different AIMD simulations that confirm the conclusions of Refs. [4,5,18,41] stating that the adsorption and relaxation of N on Ag(111) is dominated by phonon rather than by e-h pair excitations. As explained in Refs. [4,5,41], the sharp initial increase of E (N) kin of about 1 eV is a consequence of the strong attraction that the N atom experiences along its approach to the surface. After this point, the kinetic energy is reduced to a half in less than 0.5 ps. Interestingly, at the end of our simulations, the adsorbing N atoms are practically thermalized with the (hot) atoms of the topmost layer for T s = 700 K, but not for T s = 80 K.
Surface temperature may affect other characteristics of the adsorption dynamics directly related to the adsorption paths. Actually, a remarkable feature is the presence of subsurface N between the first and second layers at T s = 700 K that is not observed at T s = 80 K. Such absorption events seem to be a characteristic associated to the surface temperature since they occur whenever T s = 700 K is used as initial condition for the surface atoms in any of the simulations discussed so far (i.e., also in the microcanonical AIMD case with T s0 = 700 K). Only the number of events varies from one simulation to another and, as a general trend, e-h pair excitations contribute to reduce the number of absorption events in about a factor 2, at least within our limited statistics. A remark on the role of e-h pairs is necessary at this point. Even if the adsorption of N on Ag(111) itself is decided by the energy exchange with the surface lattice that also dominates the adsorbate dynamics at the early stage, one cannot overview the crucial contribution of e-h pairs in the subsequent stage that involves the relaxation and accommodation of the N atoms on the surface [4,5,18,41]. In practice, both energy dissipation/exchange channels affect directly the properties of the relaxation process. It is for this reason that in the rest of this section, which focuses on the details of the N adsorption dynamics, we will only discuss the results obtained from the AIMDEF calculations. Qualitatively, the conclusions from AIMD are the same.
The total traveled length at the end of the simulations is very similar irrespectively of including or not the thermostat and the value of T s . Specifically, it varies from about 63 to 66Å when changing T s from 700 to 80 K. Notwithstanding, such large mean traveled lengths are associated to comparatively small lateral displacements of 4-5.5Å measured from the initial (x i ,y i ) position of the incident N. This is the result of intricate adsorption paths that are characterized by multiple rebounds along the perpendicular and parallel directions. The time evolution of the mean lateral displacement d is shown in Fig. 4(b). We observe that after 1 ps the lateral displacement is basically constant at T s = 80 K, but it still increases slowly until the end of the simulations at T s = 700 K. In the latter case, the differential analysis of the dynamics followed by the on-surface and subsurface atoms reveals that the increase in d is exclusively due to the former group. All these results together suggest that at low T s , the adsorbing N is more rapidly trapped within the neighborhood of its final adsorption well, where it will remain until it thermalizes with the surface. In contrast, even if the adsorbing N thermalizes faster on the surface at high T s , the eventual excitation/deexcitation of the surface atoms impedes the stabilization of the atom in a close region. Finally, note that the unrealistic hot lattice created by the AIMDEF results without thermostat only causes small differences of about 0.5Å in d (compare the black and red curves). All in all, the detailed analysis of the adsorption paths shows that only the absorption process at high temperatures is a distinct characteristic associated to T s in our system. Apart from it, there are no significant differences between the different type of simulations regarding the total traveled length, number of rebounds, and lateral displacements at this stage of the adsorption process. Differences due to surface temperature or finite-size cells may start at longer times. Thus, our present analysis suggests that the commonly used microcanonical description of gas-surface elementary reaction dynamics is well justified. Similar conclusions were obtained for the adsorption of O 2 on Pt(111) [40].
Another issue of interest is to analyze possible surface temperature effects on the competition between the two main mechanisms contributing to energy dissipation, that is, phonons and e-h pairs. The theoretical analysis based on AIMDEF simulations that allow (nonfrozen surface) or not (frozen surface) the surface atoms movement showed that the energy dissipated into e-h pairs E ehp is typically overestimated when the energy exchange to the surface lattice is neglected (frozen surface) [4]. Such an effect was observed in both the adsorption of H on Pd(100), which is an example of energy dissipation dominated by the electronic channel, and the adsorption of N on Ag(111), which is dominated by energy exchange with the lattice. The reason of the E ehp overestimation is that the adsorbate velocity is larger when the competing energy dissipation channel was neglected. Similar findings and conclusions were also reported when studying the competition between phonons and e-h pairs in the H-H and N-N Eley-Rideal recombinative abstraction from W(100) and W(110) [34]. In view of all these previous results, our interest here is to determine a possible dependence of E ehp on T s .
The energy lost into e-h pair excitations for each individual trajectory can readily be calculated on-the-fly as where r i and η(r i ) are, respectively, the position vector and friction coefficient of the atom at the time-integration step n and t is the constant time step used in the numerical integration. Figure 5(a) shows the energy lost into e-h pairs averaged over the adsorbing (and absorbing) trajectories E ehp (t) as obtained from the different AIMDEF simulations. Similar results are obtained at T s = 700 K if only the (onsurface) adsorbed trajectories are included in the average. Although the E ehp (t) values are slightly larger in this case, they are still smaller than than those at T s = 80 K by about 20 meV. Precisely this dependence of E ehp (t) on T s after 1-2 ps is in principle surprising because the energy transferred into the surface lattice is expected to decrease the larger is T s in comparison to the incidence energy of the incoming atom. The latter is actually consistent with the larger increase in E (i th L) kin obtained at T s = 80 K than at 700 K (see Fig. 2). Figure 5(b) shows that the behavior of E ehp (t) cannot be explained in terms of the different friction coefficients experienced by the adsorbing N since η(t) is on average similar or even slightly 085437-6 larger at 700 K than at 80 K. Therefore, the reason should be the somewhat smaller values of E (N) kin (t) at 700 K that we can observe in Fig. 4(a) within the time interval t 1.25-2 ps. Nevertheless, we have verified that the differences in both E ehp (t) and E (N) kin (t) are in the order of the numerical energy drift occurring during the integration of the equation of motion. Therefore, we must conclude that within the accuracy of the present simulations, no sizable surface temperature effects are observed in the energy dissipated into e-h pairs for this system.

IV. CONCLUSIONS
Surface temperature effects on the adsorption of N atoms on the Ag(111) surface have been studied within the AIMD and AIMDEF schemes, by coupling the third surface layer to a Nosé-Hoover thermostat. The low and high surface temperature regimes are covered by simulations performed at T s = 80 and 700 K. We find that for this system the adsorption probabilities, the relaxation rate of the adsorbates at the surface, and their total traveled length l ads are not much affected by the surface temperature. No large differences are either obtained for the lateral displacements of the adsorbates from their initial lateral position, which are of the order of 4.5-5Å and considerably smaller than the l ads due to intricate adsorption paths with multiple rebounds. We observe, however, that whereas at T s = 80 K this lateral distance is roughly constant after 1 ps, at T s = 700 K is somewhat larger and it is still increasing slowly at the end of the simulation time (3.5 ps). It is also noteworthy that at T s = 700 K, a number of simulations lead to absorption events with N atoms remaining between the first and second surface topmost layers, but not in the case of T s = 80 K. Finally, no sizable surface temperature effect is observed in the amount of energy lost to electronic excitations evaluated by means of the AIMDEF simulations.
Another objective of this work was to test the validity of the commonly used AIMD and AIMDEF simulations performed without the inclusion of thermostats in, inevitably, finite-size simulation cells. As expected, we find that since there is no heat dissipation in this kind of simulations, the first-and second-layer atoms are unphysically heated at long simulation times, particularly in simulations aimed to describe the adsorption process at low T s . However, magnitudes related to the dynamics of the adsorbates such as the adsorption probabilities, relaxation rates, traveled lengths, and the energy lost to excitation of electron-hole pairs, are not significantly different from those obtained when the surface temperature is correctly accounted for, at least within the time scale of our simulations (3.5 ps). Therefore, considering that the adsorption of N in Ag(111) is a highly exothermic process, we can conclude that the customary AIMD (or AIMDEF) simulations that do not include coupling to thermostats are well justified to model adsorption and relaxation processes taking place at the picosecond time scale.
Still, there are processes also occurring at this time scale for which the use of thermostats can be relevant as, for example, reactions induced by femtosecond laser pulses. In this respect, it has recently been shown that the desorption of O 2 from Ag(110) [52,53] and, in particular, the desorption of CO from Ru(0001) [54] is induced to a great extent by the hot surface lattice created by the laser. Thus, we anticipate that the use of thermostats together with the AIMDEF scheme, which is also adapted for describing time-dependent electronic temperatures [55], can be a powerful tool to simulate and understand surface femtochemistry.