Ab initio molecular dynamics with simultaneous electron and phonon excitations : Application to the relaxation of hot atoms and molecules on metal surfaces

D. Novko,1 M. Blanco-Rey,2,1 J. I. Juaristi,2,3,1 and M. Alducin3,1 1Donostia International Physics Center (DIPC), Paseo Manuel de Lardizabal 4, 20018 Donostia-San Sebastián, Spain 2Departamento de Fı́sica de Materiales, Facultad de Quı́micas UPV/EHU, Apartado 1072, 20080 Donostia-San Sebastián, Spain 3Centro de Fı́sica de Materiales CFM/MPC (CSIC-UPV/EHU), Paseo Manuel de Lardizabal 5, 20018 Donostia-San Sebastián, Spain (Received 30 July 2015; revised manuscript received 16 September 2015; published 24 November 2015)

In dynamic gas-surface environments, where gas-phase atomic and molecular species impinge on the surface at energies of the order of up to a few eV, energy dissipation occurs by the excitation of electron-hole (e-h) pairs and the excitation of lattice vibrations, i.e., phonons [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17].In the adsorption processes of atomic and molecular species, dissociative as well as nondissociative, the species trapped by the surface gradually lose their energy until they become thermalized on the surface.The competition between the e-h pairs and phonon channels governs the relaxation dynamics of the transient hot species, and thus it plays a decisive role in the system reactivity properties.The reason is that it rules the traveled length and relaxation time of a hot atom or molecule on the surface and, consequently, the probability to undergo a recombination reaction with another adsorbate [18][19][20][21][22][23].
Recent ab initio molecular dynamics (AIMD) simulations with electronic friction (AIMDEF) have shown that e-h pair excitations are the dominant relaxation mechanism for hot H atoms on Pd(100) that originate from the dissociative adsorption of H 2 [16].More particularly, this channel dissipates energy at a five times faster rate than the phonon channel [10].The two main reasons behind this behavior are the long H-Pd interaction time, of hundreds of fs, and the low adsorbate-tosurface atom mass ratio, γ = m H /m Pd = 0.0094.The case of H on Pd(100) represents a limiting case.For heavier adsorbates, the relative weight of e-h pairs and phonons in the energy loss is expected to vary.The energy transfer to the substrate will be determined not only by kinetic factors, such as the value of γ and the incidence conditions, but also by the topography of the multidimensional potential energy surface (PES) and the electronic structure details of the configurations probed along the relaxation trajectory.
In this Rapid Communication, we investigate the relaxation dynamics of hot species in three adsorption scenarios that are representative of different energy loss regimes.We have chosen atomic N on Ag(111) (γ = 0.13), N 2 on Fe(110) (γ = 0.5), and the aforementioned H on Pd(100) as case studies.Our choice is also motivated by the results reported from chemicurrent experiments showing that the number of low-energy metal electrons excited during the adsorption of different gas species, i.e., the chemicurrent intensity, scales with the adsorption energy [4,5].The wide range of E ads values covered by the present case studies, between 0.2 and 2.3 eV, allows one to elucidate the observed scaling law.Modeling on an equal footing the e-h pair and phonon dissipation mechanisms is crucial to properly account for the hot species relaxation dynamics, irrespective of the kinetic factors considered.Thus, our challenge here has been to adapt the efficient AIMDEF methodology in order to accurately treat e-h pair excitations in those situations of non-negligible electron density changes produced by the movement of surface atoms.
In AIMDEF [16,24], the effect of low-energy e-h pair excitations on the gas-species dynamics is included by a friction force that is calculated on the fly within the local density friction approximation (LDFA) [25].This approximation assumes that the friction coefficient η(r A ) acting on each atom of the hot species at its position r A is the same friction coefficient the atom would experience when embedded in a homogeneous free electron gas with electron density n(r A ).The latter is the density of the bare surface at r A , which is calculated with density functional theory (DFT).Obviously, obtaining n(r A ) on the fly becomes a complicated task in simulations where surface atoms are allowed to move.The problem is that at each time step of the trajectory n(r A ) should be that of the actual configuration of the (moving) surface atoms at this precise instant.Here, we propose a scheme to evaluate the embedding electronic density that allows one to incorporate density variations due to the movement of the surface atoms.The idea consists in evaluating the gas-species electron density within a Hirshfeld partitioning analysis [26] and using it to remove its contribution from the self-consistent gas-surface electron density n SCF (r A ). Therefore, the embedding bare-surface density for each atom can be calculated on the fly at each integration step as where the summations on i and j run over all the N atoms of the system and over the N atoms conforming the gas species, respectively, and n atom i,j is the electron density associated with the isolated atoms.The Hirshfeld partitioning scheme has been successfully applied to study the vibrational lifetimes of molecular adsorbates within the LDFA framework [27].Here, we have checked that the energy losses obtained in frozen-surface simulations by using either the self-consistent bare-surface electron density or the present n sur are basically undistinguishable.
We start by revisiting the relaxation of hot H atoms formed from H 2 dissociation on Pd(100) at normal incidence and with an initial kinetic energy E i (H 2 ) = 0.5 eV.Our simulations (178 individual H trajectories) start at the instant (t = 0) in which H 2 is already dissociated (see Ref. [16] for further technical details).The adsorption cases of N atoms on Ag(111) and N 2 on Fe(110) are studied for normal incidence and initial energies E i = 0.1 and 0.75 eV, respectively.These values are chosen to ensure adsorption probabilities S 0 > 0.85 [28] and >0.7 [29], respectively, that will yield reliable statistical averages.A total of 20(80) trajectories departing from height z i = 4 Å are simulated using as the time step 0.5 fs (0.7 fs) for N atoms (N 2 molecules).Special care is taken in describing properly the N spin state at different heights.The supercell size used in each system ensures negligible lateral interactions [30].In order to disentangle the contributions of the e-h pair and phonon excitations, we perform three types of AIMD simulations on each system: (i) frozen surface with electronic friction (FS+EF), (ii) nonfrozen surface without electronic friction (NFS), where the two outer metal layers are allowed to move, and (iii) nonfrozen surface with electronic friction (NFS+EF).
Figure 1 shows, for each simulation type and system, the adsorbate kinetic energy as a function of time averaged over the available number of adsorption trajectories E A K .Starting with the H/Pd(100) case, we first note that the NFS+EF result is not significantly affected by the use of n sur instead of the bare frozen-surface density approach taken in Ref. [16].The reason is that Pd atoms are only slightly displaced by collisions with H.The figure shows the two typical features of a system in which relaxation is dominated by e-h pair excitations: (i) E A K (t) decays more rapidly in the FS+EF than in the NFS simulations and, as a consequence, (ii) the NFS+EF results, which include both dissipation channels, lie close to the FS+EF curve, which only includes the dominant electronic channel.
In contrast, the N/Ag(111) results are those representative of a system in which relaxation is dominated by phonon excitations, i.e., a much faster decay of E A K (t) when surface atom movement is allowed, which is very similar irrespective of whether electronic friction is included (NFS+EF) or is not included (NFS).Regarding the details of the dynamics, at the beginning of the N/Ag(111) simulation (t < 0.2 ps), we observe a rapid kinetic energy gain, which is the result of the barrierless strong attraction felt by the N atom towards the Ag(111) surface.This increase is followed by a rapid decrease that is a consequence of N probing the repulsive part of the PES.In the NFS and NFS+EF cases there is an additional energy transfer to the surface atoms acting as soon as N approaches the surface, as we will discuss below.Therefore, after the first collision, the N atom is left with about 0.5 eV of kinetic energy (0.8 eV in the FS+EF case) to be lost as a regular hot atom.
The case of nondissociative adsorption of N 2 on Fe(110) is even more extreme regarding the dominant role played by phonons over e-h pair excitations.In fact, there are no molecular adsorption events in the absence of surface atom movement.Furthermore, the NFS and NFS+EF results are almost undistinguishable.In contrast to N/Ag(111), we do not observe a fast kinetic energy gain upon the approach of the projectile to the surface.The reason is the presence of energy barriers at z ∼ 2.5 Å [29,31] that the molecules have to overcome before accessing the adsorption wells.As a result, we observe a drop of E A K starting at t = 0 that is the combined effect of an increase of the potential energy and of energy transfer to the surface lattice.Moreover, unlike in the two previously studied atomic adsorption cases, N 2 molecules undergo little lateral displacement as hot species on the surface.
All in all, Fig. 1 highlights that, under the usual conditions relevant for gas-surface reactions, e-h pair excitations will dominate energy dissipation only for light atoms (γ 1).In principle, the kinetic energy of a free atom in a free electron gas (FEG) decays at a rate 2η/m A , where m A is the atom mass and its friction coefficient η depends nontrivially on the atomic number Z and the FEG density (Z oscillations) [32].Nevertheless, for typical electron densities probed by hot atoms and molecules on metal surfaces, η varies slowly with Z [32][33][34] and thus the electronic decay rate is dominated by m A .As a general trend, this causes the behavior observed in Fig. 1, although the actual decay rates deviate from the freeatom values due to the PES topography of each system [16].
Next, we analyze the amount of energy transferred from the adsorbates to e-h pairs and to the lattice.For a single atomic trajectory, the e-h pair energy contribution as a function of time can be evaluated by integration of the friction force on each gas atom, where v A (t ) is the instantaneous atom velocity.The average value of this quantity over the available number of trajectories E ehp (t) is shown in Fig. 2 for the FS+EF and NFS+EF simulations performed on the three systems (for N 2 the contributions of both atoms are added).In order to extract the energy transferred to the lattice atoms, we need both the instantaneous kinetic energy of the surface atoms and the instantaneous variation of the potential energy due to lattice distortions.The former is provided on the fly in the NFS+EF simulations, and its statistical average over all the trajectories E S K (t) is also shown in Fig. 2. The energy stored as potential energy of the surface atoms E S P can be obtained from an analysis a posteriori of the simulated trajectories.For a few individual time steps {t n }, we remove the adsorbate and consider only the surface atomic coordinates.In these distorted surface configurations, we evaluate the increase in the potential energy with respect to the equilibrium bare-surface configuration.This calculation is carried out for the available trajectories and then averaged to obtain estimates of E S P (t n ) . Figure 2 shows that, for the three systems, E S P (t n ) and E S K (t) lie close to each other.This can be interpreted as a fingerprint of a harmonic oscillatory regime for the surface atoms.The average total energy transferred to the lattice, i.e., the phononic part of the energy loss, is In the e-h pair dominated H/Pd(110) system, Fig. 2 shows that after just 0.5 ps the energy transferred to e-h pairs is around five times larger than the total energy transferred to the lattice ( E ph 0.1 eV).We also observe that the energy dissipation to each of the channels is not additive, since the inclusion of phonon excitations reduces the amount of energy going to e-h pairs.This effect is more pronounced in the N/Ag(111) system, in which the energy transfer to e-h pair excitations is around a factor of 2 larger in the FS+EF than in the NFS+EF simulations after 3.5 ps.The main reason is that the average v A values entering Eq. ( 2) are lower in the NFS+EF than in the FS+EF simulations due to the additional energy release into the competing phononic channel.
There are some common features in N/Ag(111) and N 2 /Fe(110), i.e., the systems in which phonons dominate the adsorption relaxation.In the first stages of the dynamics, at t < 0.5 ps in N 2 /Fe(110) and t < 1.5 ps in N/Ag(111), a rapid energy transfer takes place from the molecule (atom) to the surface lattice.In fact, this fast energy transfer fully accounts for adsorption in the N 2 /Fe(110) system and is its main cause in N/Ag(111).After this short period of time, while E ph reaches a plateau, E ehp increases monotonically.Interestingly, we observe that, despite the fact that phononic channel dominates N and N 2 relaxation (see Fig. 1), the energy lost to e-h pairs is by no means negligible.Note that e-h pair excitation is much more important in N/Ag(111) than in N 2 /Fe(110), though the initial kinetic energy of the projectile is larger in the latter.This is, in fact, related to the larger adsorption energies of N on Ag(111) than of N 2 on Fe(110).The picture that emerges for these two systems is the following.At the first stages of the interaction, the projectile transfers energy mainly to the lattice atoms in large momentum transfer collisions and gets trapped on the surface.Subsequently, it may travel as a hot species before being accommodated in the adsorption well by exciting, mainly, e-h pairs.Therefore, the amount of electronic excitation is closely related to the adsorption energy, i.e., to the energy that the adsorbate must dissipate to get fully relaxed.This is consistent with the reported scaling of the intensity of chemicurrents measured upon adsorption of different gas species with their adsorption energy [4].
The richness of the N 2 /Fe(110) system regarding its adsorption properties allows us to better characterize this effect.In this system, the six-dimensional (6D) PES calculated within the frozen-surface approximation features two molecular adsorption minima at the top (local minimum) and hollow (global minimum) sites, where the N-N bond lies perpendicular and parallel to the surface, respectively [29,31].Remarkably, if relaxation of Fe atoms is allowed, a new molecular adsorption local minimum appears at the surface bridge site with parallel orientation (see Fig. 3 inset and Table I).The existence of such distinct adsorption wells permits one to isolate the effect of E ads on the e-h pair excitations from other factors.Figure 3 shows that the lowest E ehp occurs for adsorption on the top site, despite its E ads being about 100 meV larger than on the bridge well.The reason is that the molecules adsorbed on top are farther away from the surface, where the electron density and, therefore, the probability to excite e-h pairs is considerably smaller.The scaling of E ehp with E ads is recovered when comparing the results for adsorption on the hollow and bridge wells, since in both cases the molecule is 201411-3 close to the surface in regions of relatively high and similar density.
In summary, we have adapted the AIMDEF methodology to accurately account for the surface electron density variations caused by the motion of the surface atoms.With this feature, the AIMDEF scheme is now perfectly suited to treat from first principles all possible gas-surface elementary processes where the surface temperature and e-h pair excitations are determinant, such as adsorption, scattering, and desorption of both light and heavy gas species.This AIMDEF method is used here to study the role of e-h pair and phonon excitations during the thermalization of hot species on metal surfaces.Different representative systems allow us to extract general conclusions for reactions on surfaces.The thermalization of light hot reactants and intermediate products on the surface, e.g., H, will happen at a faster rate than for heavier species that involve C, N, and O, because e-h pairs are more efficiently excited by the lighter atoms.For heavier atoms, we find that dissipation is dominated by lattice vibrations mainly at the initial stages of the hot species interaction with the surface, along with a substantial excitation of e-h pairs that is active during long time scales.We find that more energy is diverted into the e-h pair channel for higher adsorption energies, in consistency with the experimental observations on chemicurrents, although deviations from that behavior can be induced by the particular details of the surface electron density distribution.These conclusions can only be drawn within the theory level used here.Ultimately, the weight of each channel in the energy loss is the result of the nontrivial interplay of atomic masses, PES topography, and surface electronic structure.

FIG. 3 .
FIG. 3. (Color online) Energy lost to e-h pair excitation E ehp during NFS+EF simulations for N 2 trajectories ending up at each adsorption site on Fe(110).Inset: Molecular adsorption configurations.

TABLE I .
Molecular adsorption energies E ads , center of mass height Z c.m. , and molecule bond length r at the three energy minima configurations of N 2 on Fe(110) (see Fig.3 inset).The corresponding adsorption probabilities S 0 obtained from NFS+EF simulations are also shown.