Vibrational lifetimes of hydrogen on lead films : An ab initio molecular dynamics with electronic friction ( AIMDEF ) study

Using density functional theory and Ab Initio Molecular Dynamics with Electronic Friction (AIMDEF), we study the adsorption and dissipative vibrational dynamics of hydrogen atoms chemisorbed on free-standing lead films of increasing thickness. Lead films are known for their oscillatory behaviour of certain properties with increasing thickness, e.g., energy and electron spillout change in discontinuous manner, due to quantum size effects [G. Materzanini, P. Saalfrank, and P. J. D. Lindan, Phys. Rev. B 63, 235405 (2001)]. Here, we demonstrate that oscillatory features arise also for hydrogen when chemisorbed on lead films. Besides stationary properties of the adsorbate, we concentrate on finite vibrational lifetimes of H-surface vibrations. As shown by AIMDEF, the damping via vibration-electron hole pair coupling dominates clearly over the vibration-phonon channel, in particular for high-frequency modes. Vibrational relaxation times are a characteristic function of layer thickness due to the oscillating behaviour of the embedding surface electronic density. Implications derived from AIMDEF for frictional many-atom dynamics, and physisorbed species will also be given.


I. INTRODUCTION
Molecular dynamics at metal surfaces are often influenced by non-adiabatic effects, i.e., by the coupling of molecular degrees of freedom to a continuum of electronic excitations of the metal substrate.This becomes not only apparent for surface photochemistry and electron-induced processes, for which electronic excitations are the rule rather than the exception, 1 but also for dynamics which are usually considered to proceed in the ground electronic state.Examples are line-broadening in vibrational spectroscopy, 2 deexcitation of vibrationally excited molecules during inelastic moleculesurface scattering, [3][4][5][6] molecule-surface reactions, [7][8][9][10][11] atom-surface scattering, [12][13][14] or the cooling of "hot atoms" at metal substrates. 15These processes, when dominated by vibration-electron hole pair (EHP) coupling, often take place on the timescale of picoseconds and below.Perhaps the simplest manifestation of vibration-EHP coupling is the fast vibrational relaxation of adsorbed atoms and molecules. 16n this work, we study the vibrational relaxation of the simplest adsorbate, a hydrogen atom, on overlayers of lead of increasing thickness.It will be argued that also for this system (not yet studied by experiment so far according to our knowledge), coupling of the H-surface vibration(s) to EHPs is the most important channel for energy dissipation.0][21][22] QSEs arise from the fact that upon layer-by-layer growth of a film, the resulting film thicknesses are usually not commensurate with the wavelength of the particle-in-a-box like electronic wavefunctions along the zdirection (perpendicular to the surface).This leads to nonmonotonic behaviour of certain features with film thickness.For example, it was found experimentally that Pb overlayers grown on solid substrates, show characteristic variations of the electron density above the film as it becomes thicker. 17n these experiments, constant electron densities ρ 0 were probed by low-energy He atom scattering (HAS) from Pb films on copper substrates.The HAS signals reveal characteristic oscillations of the "apparent step heights" with increasing number of Pb layers, N, 17 i.e., different heights above the film which correspond to the particular electron density value probed by the He atom.In another key experiment, the nature and the decay of excited "quantum well states" in Pb(111) films grown on Cu(111) were studied by twophoton-photoemission spectroscopy. 18Characteristic oscillations with film thickness were observed for the lifetime of the lowest unoccupied quantum well states, 18 and traced back theoretically to discontinuous contributions of electron-electron scattering to spectroscopic linewidths. 21ere, we investigate possible QSEs on Pb(111) films when covered with an adsorbate, hydrogen in this case.We expect similar characteristic variations for properties associated with the varying embedding density for H in front of the surface.In particular, we are interested in vibrational relaxation times as a function of layer thickness.In this context we also mention the effect of size quantization on electronic level broadening of H atoms near thin Al(111) films, as predicted theoretically in Ref. 23.
In order to tackle this problem, we use density functional theory (DFT) in a supercell approach, together with Ab Initio Molecular Dynamics with Electronic Friction (AIMDEF) for vibrational dynamics.This approach, recently developed to treat the coupling of nuclear motion to EHPs of a metal surface as "electronic friction," 15 is a parameter-free "on the fly" classical dynamics method which avoids the a priori determination of a global potential surface, and also the explicit inclusion of electronically excited states.The latter are implicitly included via atomic friction coefficients, obtained within the Local Density Friction Approximation (LDFA). 8In the LDFA, the friction coefficient is calculated from the embedding density of the clean, adsorbate-free surface.The method is based on well-defined approximations which make it applicable for multidimensional molecular, quasi-nonadiabatic dynamics including anharmonicity, mode coupling, and surface atom motion.
The paper is organized as follows.In Sec.II, the models and numerical methods will be outlined.Section III presents and discusses results obtained for adsorption and vibrational relaxation of hydrogen on Pb(111) overlayers of increasing thickness.Finally, Sec.IV summarizes and concludes our work.

II. MODELS AND METHODS
We use DFT to study the adsorption of hydrogen atoms on lead films of (111) orientation and increasing thickness.For systems with N layers of lead, we used hexagonal supercells with in-plane lattice parameters corresponding to the experimental lattice constant of a = 4.95 Å of the bulk metal, i.e., a 1 = (a/ √ 2, 0, 0) and a 2 = (a/(2 √ 2), a √ 3/ √ 8, 0), respectively.Each cell contains one Pb atom per layer, corresponding to a Pb-Pb distance of a/ √ 2 = 3.50 Å within a layer.In most examples below, unrelaxed Pb N films were considered with fixed interlayer separations of d = a/ √ 3 = 2.86 Å.Effects of lattice relaxations were also studied.The third lattice constant of the supercell was chosen as a 3 = (0, 0, (N + 5)d), i.e., with a vacuum gap corresponding to 6d = 17.15Å. Test calculations with constant cell size, namely, a 3 = (0, 0, 17d) for all N, showed only marginal differences to the "constant gap" approach.The stacking sequence of layers is ABC-ABC. . .as required for a facecentered cubic (fcc) metal.The same unit cell was adopted for Pb surfaces when covered with atomic hydrogen.For H adsorption on Pb films consisting of N layers (often denoted as H-Pb N below), a single H atom is placed on the topmost layer at various positions.In particular, the fcc and hcp (hexagonal closed packed) sites, as well as atop and br (bridge) sites were considered, cf.Fig. 1 for an illustration.This corresponds to a high coverage situation (coverage = 1), i.e., H:Pb N (111)(1 × 1).
All calculations were performed with the Vienna Atomistic Simulation Package (VASP). 24For Pb and H, we applied the Projector Augmented-Wave (PAW) method 25 and the generalized gradient approximation (GGA) by Perdew-Burke- Ernzerhof (PBE) for the exchange-correlation functional. 26or Pb, only the four valence electrons were explicitly considered.A plane-wave basis with an energy cutoff of 250 eV was adopted, and a Monkhorst-Pack k-point mesh 27 of size 9 × 9 × 1. Convergence tests have been carried out for bulk Pb.Geometry-optimization of the lattice parameter for bulk Pb gave a the = 5.03 Å, i.e., the same slightly elongated value (compared to experiment), as in Ref. 20.
Besides stationary calculations, AIMDEFs were performed to study vibrational dynamics of a single H atom, using a new version of VASP as developed in Ref. 15. H moves under the influence of frictional forces, the latter due to coupling to EHPs in the metal.Specifically, for H we solve a Langevin-type equation of motion, Here, m H is the hydrogen mass, −∇ H V the force acting on H due to the conservative potential, and −η H dr H dt is the frictional force.Both forces are calculated "on the fly."The friction coefficient η H is determined within the LDFA. 8Here as in Refs.8 and 15, η H is calculated at every position of the H atom visited along the trajectory r H (t), which probes different embedding densities ρ(r H (t)).The friction coefficient was extracted from the scattering properties of the Kohn-Sham wavefunctions in a static DFT calculation for a H atom embedded in a homogeneous electron gas with density ρ 0 = ρ(r H ). 28 The numerical results η H = f(ρ 0 ), formally given as η H = A r 2 s l (l + 1)sin 2 [δ l − δ l+1 ] (with A = (12/π ) 1/3 , δ l the Fermi phase shift of the lth partial wave and r s explained shortly, see Refs. 8, 29, and 30 for details), were then fitted to the analytical form ( The dependence of η H is now expressed as a function of the mean electron radius r s of the free electron gas, where r s = [3/(4πρ 0 )] 1/3 .In order to faithfully reproduce phase shifts and friction coefficients in an interval of r s values between r s = 1 to r s = 10 a 0 (bohr), i.e., from rather high (low r s ) to rather low densities (high r s ), we use a new fit (relative to the one that was used in Ref. 15) with M = 3 and nine fit parameters a 1 to c 3 as given below (in the caption of Fig. 6).This new fit extends the fit used in Ref. 15, which was obtained for r s = 1 − 6 a 0 , to lower densities.The use of the LDFA, based on the independent atom approximation, can be controversial if, unlike here, the coupling of molecular (internal) vibrational modes to EHPs were of interest.However, in Ref. 30 it was demonstrated (for Hcovered Ru (0001)) that LDFA gives quite reasonable results even for vibrational modes when compared to methods such as tensorial friction 7 and Golden-Rule type expressions, 31 both of them based on gradient-corrected DFT.For single H atoms as studied here the LDFA method is expected to be very accurate.In fact, in experiments in which the energy loss of different atomic particles scattered from metal surfaces was measured, it has been proven that the method provides friction coefficients that reproduce accurately the experimentally obtained friction force close to the surface. 32quation (1) has been propagated by Beeman's predictorcorrector scheme as outlined elsewhere. 15We will monitor the vibrational relaxation of selected modes of H adsorbed on the surface, due to coupling to electron-hole pair excitations in the metal.If not stated otherwise, we will choose as initial conditions for Eq. ( 1) the equilibrium position of the H atom in a particular adsorption site, r H (t = 0) = r eq H , and an initial momentum of H corresponding to |p Here, E vib is the vibrational energy of (one quantum of) the chosen vibration of H relative to the surface, calculated from normal mode analysis (diagonalization of the dynamical matrix), NMA.Typically, we fix the positions of the Pb lattice during normal mode analysis, resulting in three modes for the vibrating H: One "high-frequency" mode perpendicular to the surface, often called the "z-mode" in the following, and two (ideally degenerate) low-frequency modes parallel to the film, the "(x,y)-modes."The direction of the initial momentum is chosen along the computed normal modes and its sign to be negative, i.e., in the case of the z-mode this corresponds to initial motion towards the surface.(Choosing the other initial direction will merely phase-shift the vibration, causing a delayed onset of dissipation.) In most cases, the Pb lattice was held rigid also during time-evolution, however, in selected cases the lattice was allowed to move.This way we can study vibration-phonon coupling and effects of lattice relaxation on the timescales considered here (∼1 ps).Even for moving lattices, friction was restricted to the H atom only, with friction coefficients taken for a rigid lattice.(Test calculations have shown that for the small lattice displacements present in this problem this is an excellent approximation.)

A. Adsorption of H on Pb films of increasing thickness
In a first step we studied free-standing lead overlayers, without an adsorbate, similar to Ref. 20.For Pb films with (111) orientation and increasing number, N, of layers (from N = 1 to N = 12), we obtain total energies E(N)/N per layer as shown in Fig. 2(a) as a function of N.An oscillatory behaviour for the quantity becomes clearly visible in Fig. 2(b).Here, "total energy" refers to cohesive energies where the energies of the reference atoms have been subtracted.The energy (per Pb atom) obtained for the bulk metal (with a = a exp = 4.95 Å and with a 6 × 6 × 6 k-point mesh) is shown for comparison in Figs.2(a) and 2(b).These oscillations are in almost quantitative agreement with those found in Ref. 20, where 14 rather than four valence electrons were considered per Pb atom, and the Perdew-Wang exchange-correlation functional 33 instead of PBE.In Ref. 20, it has been shown that these oscillations are due to quantum size effects perpendicular to the surface.They lead to nonmonotonic behaviour for certain features such as E(N), or the "spillout" of electron density as outlined above.From Fig. 2(b) we note that the oscillations in E(N) have a period of two, up to N = 6, becoming more irregular thereafter.Also, the amplitude of the oscillations is damped towards medium N, then slightly increasing again, in agreement with Ref. 20.In that reference it has also been demonstrated that lattice relaxation and the inclusion of geometric constraints to mimic a support on which the overlayer is grown has quantitative, but no qualitative, impact on energetic QSE in thin lead films.
As a next step, the chemisorption of H atoms at high coverages, H:Pb N (111)(1 × 1) was considered.Different minima of a single H atom adsorbed on Pb N (111) were found by geometric relaxation, starting from different initial geometries.The Pb lattice was initially fixed at the geometries of the adsorbate-free films.Effects of reoptimizations of the lattice will be considered shortly.The adsorption energies are , when H is adsorbed in fcc sites (a) or hcp sites (b), respectively.The upper panels show adsorption energies E ads , the middle panels the optimized distance of H atom above the top Pb layer, z H , and the lower panels the energy of the highest-frequency vibration (z-mode) of H. Crosses in the upper panels mark adsorption energies in bridge sites (for N = 1) and atop sites (for N = 5), respectively.calculated as differences between free and adsorbed state energies, Here, E H−Pb N is the PBE energy of the geometry-optimized adsorbed species, and is approximated by a single-point calculation with the adsorbed H atom being displaced from its equilibrium geometry along z by three Pb interlayer distances (i.e., by 3d = 8.57Å). (An approximation, which has been checked to introduce only a marginal error ≤1 meV).Again, up to N = 12 layers were considered.
As expected and known for other metal surfaces, hydrogen prefers high-symmetry sites with high metal coordination, also for H:Pb N (111) (1 × 1).All adsorption energies found here are around and slightly above 2 eV per H atom.In particular, the fcc and hcp sites in the center of a triangle of three top-layer Pb atoms are typically most stable.In the hcp site, the H atom adopts the same lateral coordinates as the next Pb layer below the top Pb layer, whereas fcc has lateral coordinates equal to those of the second-next layer below the top layer.The adsorption energies for H in fcc and hcp sites on Pb N films are shown in the upper panels of Figs.3(a) and 3(b), respectively.For most overlayers, the fcc site is slightly more stable than hcp, typically by about 0.04 eV.The lowercoordination high-symmetry sites such as "bridge" and "atop" are usually less stable.For example, for H-Pb 5 , the "atop" site is destabilized with respect to fcc by 0.18 eV and 0.14 eV less stable than hcp (crosses for N = 5 in the upper panels of Figs.3(a) and 3(b)).
An interesting exception to the high-coordination rule is H-Pb 1 , i.e., H adsorbed on a single Pb(111) layer, where the "bridge" adsorption energy (with H coordinated by only two nearest-neighbour Pb atoms), is 2.40 eV and therefore by 0.17 eV more stable than the fcc = hcp site with E ads = 2.23 eV (crosses for N = 1 in the upper panels of Figs.3(a) and 3(b)).For H-Pb 1 , fcc and hcp sites are identical, since there are no further, deeper-lying Pb layers.Obviously due to the lack of lower-lying layers, H tries to maximize its interaction with Pb atoms by adopting a position with lower coordination, but shorter H-Pb distance.In fact, the adsorption energy of 2.40 eV is the highest found for all systems studied in this work.Another, related exception is N = 2, where hcp becomes more stable than fcc, by about 0.05 eV, while for N ≥ 3 the fcc site is always most stable.For H-Pb 2 the enhanced stability of hcp is due to the absence of a third layer of Pb (making the definition of hcp and fcc less straightforward), and the hcp site can more favourably interact with the "in-phase" Pb layer below the top Pb layer than the fcc site, which is "out of phase." Before discussing more details, a few general remarks are in order.First, the approximations made in our approach which are, in addition to the "standard" GGA-DFT, finite cutoff, pseudopotential, and finite vacuum gap approximations, the neglect of lattice relaxation and film support will not allow one to calculate adsorption energies with accuracies of (a few tens of) meV.Also, the calculations were done in spinunpolarized fashion.However, we expect that within the same computational model and settings, trends obtained for different sites and layer thicknesses are meaningful.
Given the classical binding energy of about 4.7 eV for molecular hydrogen, also the question arises if atomically adsorbed hydrogen at Pb(111) surfaces is thermodynamically stable against H 2 formation.In fact, the studied coverage may indeed be only metastable at finite temperatures.While this would not affect any of our findings below, let us elaborate on this point a bit further.For this purpose we did test calculations for a Pb(111)(2 × 2) unit cell with five layers, and four Pb atoms per layer.Covering with four H atoms in fcc sites, corresponding to coverage 1, and optimizing, very similar geometries and total energies were obtained as for the (1 × 1) unit cell.From spin-unpolarized calculations we get an adsorption energy for a single H of 2.28 eV, obtained from displacing one of the four H atoms upward by 3d.This is higher than the 2.16 eV for the (1 × 1) cell (cf.Fig. 2(a)), where formally every H atom is "desorbed."For H:Pb(111)(2 × 2), coverage 1, normal mode analysis shows real frequencies only, which proves the minimum character of the found geometry.
By removing two of the H atoms from nearest-neighbour fcc sites, we have coverage 1/2, which is also a true minimum, with a binding energy of 2.31 eV per H.Further reducing the coverage to 1/4 (one H per (2 × 2) unit cell) gives a binding energy of 2.33 eV for the fcc site, and a stable minimum again with real frequencies only.Thus, the binding energy increases slightly at lower coverages.
We also found that H atoms moving on the surface encounter diffusion barriers.For example, at coverage 1, a H atom diffusing from a hcp site towards a neighbouring H atom in a fcc site, feels a barrier of about 0.22 eV.Once these barriers have been overcome, H 2 molecules may form: For coverage 1, the formation of H 2 is exoenergetic by 1.59 eV, according to spin-unpolarized calculations for (2 × 2) where two neighbouring H atoms (out of four) have been moved well into the vacuum and the H-H distance was optimized.From a naive picture one would expect a reaction energy near zero, namely, 4.7 eV (the classical binding energy of H 2 ) minus 2 × 2.3 eV (twice the binding energy of H on Pb 5 (111)).Indeed, for spin-polarized calculations for a single, desorbed H and enforcing one more up than down spin, the adsorption energy per H atom becomes 1.47 eV for (1 × 1) and 1.44 eV for (2 × 2), both at coverage 1. Spin-polarized calculations for the adsorbed states, in this case without enforcing N ↑ − N ↓ = 1, converged to the spin-unpolarized solutions.With the modified adsorption energies per H atom, the energy balance is much closer to the expected value.
To summarize, we conclude that a H-covered Pb(111) film is not particularly stable at high temperatures and high coverages, but metastable.Stability increases at low temperatures and if coverages are not too high.Note that the numerical values for adsorption energies per H atom in Fig. 3 will be (constantly) shifted down by several tenths of an eV in spin-polarized calculations.Since we focus on adsorbed states below, we continue with spin-unpolarized calculations.
Returning now to the trends for (spin-unpolarized) adsorption energies of H at Pb(111)(1 × 1) again, from both upper panels of Fig. 3 we note an oscillatory behaviour of E ads with increasing N. The oscillations are most pronounced for low N values, dying off at larger N.They follow, up to N = 7, nicely those of E(N) in Fig. 2. We find that for up to N = 6, even N corresponds to extra stable Pb N layers (cf.Fig. 2) whose adsorption energies are comparatively small.Differences between the adsorption energies at fcc sites, for example, are up to more than 0.4 eV (or about 20% of a typical adsorption energy), which should be measurable, perhaps by thermal desorption spectroscopy, TDS.Further, Figs.3(a) and 3(b), middle panels, show the optimized height z H (N ) = z eq H (N ) of adsorbed H in fcc and hcp sites, as a function of N. Also here, characteristic oscillations are observed.At least in the low-N regime up to N ∼ 6 where oscillations are most pronounced, a clear correlation between z H (N) and E ads (N) exists, with larger adsorption energies corresponding to longer bond lengths.The bond lengths oscillate for the fcc series, with a maximal "amplitude" of about 0.08 Å around a mean of ∼1.07 Å.For adsorption in hcp sites slightly longer bond lengths are predicted by theory, and oscillations seem less pronounced.Overall, differences in bond lengths between fcc and hcp and between different layer thicknesses are small but non-negligible.A possible explanation for the observed correlations is as follows: For the smaller N at least, layers with even N are extra stable and have a smaller electron spillout compared to N odd. 20An adsorbing H atom needs thus to get closer to the substrate to find its optimal embedding density, causing shorter z eq H values. Finally, Figs.3(a) and 3(b), lower panels, show the vibrational energies E vib = ¯ωz for the highest-frequency mode of H on rigid Pb films, obtained from NMA.As stated, these high-frequency modes are pure z-vibrations normal to the surface.Here, the plausible trend emerges of most tightly bound species having the highest frequencies (with N = 1 being exceptional again), with the consequence that oscillations with film thickness also arise for them.For fcc sites, the harmonic frequencies of z-modes are centered around 88 meV (∼700 cm −1 ), with variations at most of about 10 meV or 80 cm −1 .Again, this is not out of experimental reach and could be accessible, e.g., by optical methods, by electron energy loss spectroscopy, EELS, or by HAS.
As mentioned above, the other two (x,y)-modes are degenerate (not quite in computational practice), describing frustrated translations parallel to the substrate.For the fcc sites, the corresponding vibrational energies vary between about 15 meV (∼130 cm −1 ) for N = 3, and close to 40 meV (∼310 cm −1 ) for N = 8.We refrain here from listing the exact values, since these are not easy to converge numerically and also because the harmonic approximation becomes questionable for the low-frequency modes.Sufficient to note that also here, an oscillatory behaviour with N is found, typically with higher frequencies for even N.We mention that the normal-mode analysis gives only real frequencies for the fcc sites, while lateral modes are sometimes characterized by imaginary frequencies for the hcp sites, questioning the true minimum-character for hcp in these cases, at least within the rigid-surface approximation.We further note that the Debye frequency of bulk lead is around 10 meV.For Pb films (on a Ni substrate), Debye frequencies nearly twice as high have been found. 34We therefore expect that the low-frequency (x,y)vibrations can more efficiently couple to substrate phonons in comparison to the z-modes.
In a further step, the rigid-surface approximation has been tested.As an example, the adsorption of H on a Pb 5 layer was studied by allowing the upper four Pb layers to relax, while keeping the bottom one fixed.In this case the optimized H-Pb 5 structure (fcc adsorption) is characterized by alternating distances between Pb layers, but the total and adsorption energies are hardly affected.In the rigid-surface case, the (constant) interlayer distance is d = 2.86 Å, changing to the alternating series d 12 = 2.75, d 23 = 3.14, d 34 = 2.91, and d 45 = 2.93 Å, respectively, if the lattice was allowed to relax.(d ij denotes the distance between layers i and j, with i = 1 being the top layer.)At the same time, the total energy is lowered by merely 45 meV.The optimized height of H above the top layer, stays 1.08 Å as before.Displacing the adsorbed H atom by 3d upward (as was done for the rigidsurface case to obtain adsorption energies), and reoptimizing the lattice, an adsorption energy of E ads = 2.18 eV is found, instead of 2.16 eV for the rigid surface.(Without reoptimization of surface degrees, we obtain a value of 2.20 eV.)The effects of surface relaxation on vibrations of H-covered films are also small: The z-mode of H-Pb 5 (fcc) changes its vibrational energy from 89.7 meV (723 cm −1 ) for the unrelaxed, rigid surface, to E vib = 90.3meV (728 cm −1 ) for the relaxed surface.In summary, the effects of surface relaxation on stationary properties of H-Pb N are present, but small.Nevertheless, for selected dynamical calculations surface atom motion will be included.
In this context it is worthwhile to discuss the issue of substrate atom motion slightly deeper.When including substrate atoms in the NMA we get substrate "phonon" modes in addition to adsorbate vibrations, with the highest (nearly pure) Pb mode being 8.3 meV in energy.When using a naked Pb 5 film, optimizing the interlayer distances, and doing a normal mode analysis (with the lowest layer fixed), we obtain "phonon" modes ranging from 0.6 to 8.9 meV.Note that these frequencies are -point frequencies only, and not arising from true phonon calculations with k-dispersion.Experimentally, Braun et al. 35 found -point phonon energies for a Pb 5-layer slab (on Cu(111)), ranging from 0 meV (for the so-called α 1mode) to 9.4 meV (for the ε 2 -mode).Thus, given our simplified slab model the agreement between measured and computed "lattice phonon" energies is quite good.

B. AIMDEF and vibrational lifetimes of the z-mode
We now study vibrational relaxation of excited H atoms, due to EHP-coupling.We begin again with rigid, unrelaxed Pb films.As outlined in Sec.II, relaxation is studied by solving the Langevin-type equation ( 1) with friction coefficients obtained from the LDFA, and H starting in a particular adsorption site with momentum directed along a selected normal mode corresponding to a kinetic energy E kin = E vib .
In Fig. 4, we monitor the time evolution of H adsorbed in an fcc site of Pb 5 , when excited along the z-mode initially.The zero of total energy refers to the rigid surface with H in its equilibrium position, unexcited.The upper panel shows the total energy without (black curve, i.e., η H set to zero) and with friction included (blue curve, i.e., η H = 0), respectively.We note a decrease of E tot towards the expected final excess energy of zero, for a H atom which came to rest in the fcc site.The kinetic energy of the H atom is shown in the lower panel of Fig. 4, again for frictional and non-frictional cases.The kinetic energy for frictional dynamics reflects a damped, quasi-harmonic oscillation, with an oscillation period of t vib (AIMDEF) = 47 fs, which is close to the harmonic , when H was in the fcc site and excited with one vibrational "quantum" of the z-mode with E vib = 89.65 meV, initially moving towards the surface.Upper panel: Total energies for the non-dissipative and dissipative cases, with dashed red curves indicating energies expected after t = τ vib and t → ∞, respectively (see text).The zero of total energy refers to the rigid surface with H in its equilibrium position, unexcited.The vibrational lifetime obtained from AIMDEF is indicated.Lower panel: Kinetic energy of the H atom for both cases.vibrational period of t vib (HO) = 46 fs as obtained from normal mode analysis.As a small detail (hardly visible on the scale of the figure), a slightly longer oscillation period of the undamped motion compared to the damped one is found.
At the classical turning points of the vibrating H, its velocity is zero and therefore the frictional force vanishes.As a consequence, one finds narrow plateau regions in the E tot curve around times corresponding to the classical turning points.The largest energy loss rate is observed when the kinetic energy of the H atom is high, i.e., in between two turning points.For this reason the decay of E tot in the upper panel of Fig. 4 is not strictly exponential.Nevertheless, a vibrational lifetime τ vib (AIMDEF) can be defined as the period after which the total energy has dropped to a value E tot = E vib /e, which is shown as the upper, dashed red curve in the figure .Here, we find τ vib (AIMDEF) = 343 fs, indicating a sub-picosecond relaxation of vibrating H due to electronic friction.
It should be noted that in practice a precise determination of τ vib by AIMDEF is hampered by the non-exponential decay, but also due to numerics.Small timesteps are needed to prevent a "drift" of the total energy even for the nonfrictional case (black curves in the figure), which can have an influence on the numerical value of τ vib .In the example of Fig. 4, a timestep of t = 0.1 fs was chosen for solving Eq. ( 1), leading to a slight upward drift of E tot of about 0.6 meV at t = τ vib (AIMDEF), or about 0.7% of the vibrational quantum E vib = ¯ωz of the z-mode of H-Pb 5 (fcc).
We have also determined τ vib (AIMDEF) for the z-modes of H adsorbed in fcc sites for N = 1, 2, 3, 4, and 6, using the  5), for H adsorbed in fcc sites of Pb N (N = 1 to 6, see numbers attached to the data points).For the AIMDEF calculations, short timesteps for which the drift in the dissipation-free case was small or negligible, were used (see text).
same methodology.In most cases an even shorter timestep ( t = 0.025 fs) had to be adopted to keep the non-frictional drift of E tot to values below 1% of ¯ωz at the resulting lifetime, τ vib (AIMDEF).All of these lifetimes are of the same order of magnitude as for N = 5, however, with characteristic variations.Before discussing these variations, we note that a reliable estimate of τ vib (AIMDEF) arises from the fact that the vibrational amplitudes are small, hence only very similar r s values are probed by a given trajectory.Therefore, without doing any dynamics, we can approximate τ vib by the simplifying assumption that during a vibration, r s does not vary at all and is equal the value at the computed equilibrium position of H, r eq H , with corresponding r s value, r eq s , i.e., . ( Note that a similar expression has also been used by Persson and Hellsing to estimate electronic damping of adsorbate vibrations at metal surfaces. 36Figure 5 shows indeed a good correlation between the model vibrational lifetimes τ vib (model) and the numerical lifetimes obtained from AIMDEF, τ vib (AIMDEF), for (converged) AIMDEF values for N = 1 to 6.With this in mind, we can estimate the vibrational lifetimes easily from Eq. ( 5), for all H-Pb N systems with H adsorbed either in fcc, or hcp sites.The result is shown in Fig. 6(a) for fcc sites.
From there we note qualitatively similar behaviour as before for other quantities, i.e., a damped oscillation.It is seen that the vibrational lifetimes cover a range between 315 and 345 fs.For the hcp sites (not shown), lifetimes exhibit a similar N-dependence with slightly longer values scattering between 330 and 360 fs.The observed trend is easy to understand: The τ vib curve follows the equilibrium distance z H of Fig. 3(a), middle panel.The closer the atom to the surface, the higher the embedding density ρ, the lower r s and the higher the friction coefficient η H .For fcc chemisorption, the r eq s values vary between 3.44 a 0 for N = 1 and 3.29 a 0 for N = 4. (For hcp, the variation is between 3.54 for N = 2 and 3.36 a 0 for N = 4.) For fcc, this results in friction coefficients η H between 0.129 (N = 1) and 0.140 (N = 4) atomic units.The variation of η H with r s is plotted in Fig. 6(c), with the range of r s values relevant for chemisorption of H in fcc sites indicated as a shaded area.
Because of the strong correlation between z eq H distances and vibrational lifetimes, the variation of the "electron spillout" caused by QSEs seems not to be the major determinant for the variations of τ vib .In fact, when we replace in Eq. ( 5) the r s value r eq s for the optimized adsorption geometry (at z eq H ), by an r s value obtained from fixing H at a constant height above the fcc site instead, one can eliminate effects due to different bond lengths.We then obtain a "lifetime" within this constant-height scenario, which is shown, for a distance of 1.1 Å above fcc (which is close to most optimized z eq H values), as the dashed curve with stars in Fig. 6(a).
Characteristic variations of the "lifetime" with film thickness also emerge in the constant-height modus, which are now only due to variations of the electron density spillout (and hence r s ).Note, however, that now the computed variations are, at least for the lowest N-values (N = 2 in particular), not in phase with those observed for fully optimized geometries, and also the range of variations of τ vib is smaller.Both indicate that indeed, the different spillout of electron density is not dominant at distances typical for bond lengths of chemisorption.Instead, the variation of equilibrium z eq H values is more relevant.
From an experimental viewpoint the predicted variations of lifetimes may be difficult to measure: The effect is in the order of 10% of a signal which is difficult to measure by itself.In a time-resolved experiment a resolution of a few tens of fs would be needed, and in a linewidth measurement (e.g., from HAS) one must be able to detect linewidth changes of roughly 0.2 meV of a signal of width of about 2 meV.(Assuming that electronic friction is the only channel which causes line broadening.) We expect absolute smaller effects on linewidths, but larger contrasts, i.e., relative differences between different layer thicknesses, for (H) distances farther away from the surface.To this end, we simply place and fix, the H atom 2.5 Å above an fcc site of H-Pb N (111), and use Eq. ( 5) to estimate relaxation times for these hypothetical arrangements.The relevance of this approach for physisorbed species will be outlined shortly.Let us first discuss the result as shown in Fig. 6(b), where lifetimes in case of hypothetical "physisorption" are plotted as a function of N. One finds relaxation times which are now in the order of several picoseconds.Their variation with N, i.e., the contrasts, are indeed large, particularly at low N: Oscillations are clearly visible, with a maximal amplitude of close to one ps (about 40% of the mean value).The oscillations extend to larger N-values than previously, but again in damped fashion.Now the larger differences between various film thicknesses are due to larger variations of r s values in the low-density regime farther away from the surface, compared to the high-density chemisorption regime.The range of r s values at constant z H = 2.5 Å is between 7.56 and 8.56 a 0 (see the second shaded area in Fig. 6(c)).In the end, this leads to larger, relative variations of vibrational relaxation times.Note also from Fig. 6(b) that up to N = 6, the phase of the oscillations for τ vib is shifted relative to the "chemisorbed" species of Fig. 6(a).For example, the lifetime for N = 3 is a maximum in Fig. 6(a), but a minimum in Fig. 6(b).
The low-density, "spillout" region, QSE associated with it, and its relevance for low-energy HAS has been emphasized above and elsewhere. 17,20 ere, we argue that this regime will in fact be most relevant for physisorbed species-which is why the terminus "physisorbed" was used in Fig. 6(b).Namely, some time ago Wöll and co-workers have demonstrated by HAS that the linewidth of the so-called FT z -mode (Frustrated Translation along z) of hydrocarbons on metal surfaces can have a strong, possibly dominating contribution from vibration-EHP coupling. 37They were further able to show that hydrocarbons such as ethylene or ethane can be physisorbed also on Pb(111), exhibiting the FT z -mode there as well. 38Typical physisorption wells for hydrocarbons on metal surfaces are around 2.5 Å above the surface. 38We suggest as a possible approach to obtaining vibrational lifetimes of polyatomic, physisorbed molecules to perform a normal mode analysis for the adsorbed species, followed by AIMDEF along the lines given above.We are fully aware of the fact that the current theoretical approach will be insufficient to do so, since dispersion interactions are not included so far.Also, standard (semi-)local exchange-correlation functionals may fail in describing accurately enough the low-density region several Å above a metal surface. 39(See, however, also Ref. 20, where this last point was found to be less severe for the interpretation of HAS experiments.)Nonetheless, based on Fig. 6, we can already now predict measurable oscillations in (picosecond) lifetimes of organic molecules, when adsorbed on Pb films of increasing thickness.

C. The role of phonons
So far we neglected surface atom motion in AIMDEF, i.e., vibration-phonon coupling.Indeed, for the z-modes and not too large initial energies at least, this is an excellent approximation on the sub-ps timescales relevant for vibrational relaxation of H on Pb N .
To demonstrate this, we plot in Fig. 7(a), like in Fig. 4(a) the total energy of H-Pb 5 within the first ps, when H was adsorbed in the fcc site on the rigid surface and initially excited along the z-mode, with a kinetic energy of ¯ωz (blue curve).Similarly, Fig. 7(b) shows once more the corresponding (total) kinetic energy in blue.In addition, we provide in the same figures as black curves, the analogous quantities when the top four layers of the five-layer slab were allowed to move, while keeping the bottom layer fixed.From Fig. 7(a) we note that up to the previous vibrational lifetime of 343 fs (cf.Fig. 2), there is almost no difference between the E tot (t) curves, when either surface relaxation was accounted for or not.After about 500 fs or so, we see a small difference with the energy for the nonrigid surface decreasing slightly more slowly than for the rigid surface.However, even at 1 ps the difference is merely a few meV.In Fig. 7(b), we see that the kinetic energy (apart from its oscillating nature) continuously drops for the rigid surface, while this is not so in the non-rigid surface case.We recognize after a few hundred femtoseconds the excitation of phonons, with a slow beating pattern on the timescale of several hundred fs, corresponding to a "mean phonon energy" in the order of 10 meV or a frequency of several THz.This is consistent with the lattice vibrations of H-Pb 5 as described above.
The observation of a slightly slower loss of total energy in Fig. 7(a) for the non-rigid surface may seem counterintuitive at first glance, because vibration-phonon coupling provides an additional energy loss channel.We wish to discuss this effect, albeit small, a bit more.First we note that Figs.7(a) and 7(b) refer to the total system, i.e., H plus substrate, and not to the H subsystem alone.Further, in our model true dissipative dynamics is only due to H-vibration electron hole pair coupling while vibration-lattice phonon coupling conserves total energy, i.e., a faster decay of excess energy upon including the latter cannot be expected.Rather, if part of the energy of vibrating H is transferred to the lattice, it cannot be dissipated anymore.In order to do so, also the lattice phonons would have to be coupled to a thermal bath.
To gain insight into subsystem and long-time dynamics, we also studied features of the vibrating H alone, and we considered longer propagation times, up to 10 ps.In Fig. 7(c), the velocity of the H atom (only the z-component) is shown, during the first ps, for rigid and non-rigid surfaces.Differences between both cases appear to be small again.However, closer inspection reveals that at early times (during the first two vibrational periods), the deceleration of H is a tiny bit smaller in the rigid-surface case, but becomes larger towards longer times.(Also, the oscillation period becomes slightly longer.)This long-time effect-faster energy loss in the adsorbate subsystem for the rigid surface-is more clearly seen when we go beyond 1 ps, as in Fig. 7(f), where v z (H) is plotted on a narrow ordinate interval between −0.005 and +0.005 Å/fs, but now up to 10 ps.There we see that in fact the loss of kinetic energy for H is never complete if the surface atoms are allowed to move, in contrast to the rigid surface where the H atom comes to a full rest after about 4 ps or so.If lattice phonons are accounted for, those are excited along the trajectory, establishing some sort of thermal, phonon "heatbath" with "thermal fluctuations" becoming clearly visible in Fig. 7(f), up to the final propagation time.Our system size is presumably too small and propagation times too short to really speak of a temperature, which is why we use the expression "thermal" with reservation.Nonetheless, as a consequence of the substrate motion and associated impulsive "kicks" experienced by the H atom, which then suffers friction, the total energy decreases more slowly beyond 1 ps when substrate motion is accounted for, compared to the rigid surface.This is demonstrated in Fig. 7(d), which shows E tot of the total system up to 10 ps.Note that at around 2 ps, the total energy for the non-rigid surface becomes lower than for the rigid surface.In fact, two different dynamical periods can be observed: A first period up to about 2 ps, with a fast dissipation of energy due to frictional motion of the vibrating H atom. Thereafter in a second period, the H atom came basically to rest and only the occasional hits by surface atoms lead to dissipative motion of H and an associated, slow energy loss.Eventually, at very long timescales, the total energy will be the one of the fully relaxed system (H and substrate).From our relaxed-surface optimizations as mentioned above we found this value to be about −45 meV on the energy scale given, which refers to the rigid surface.After 10 ps, the total energy in our AIMDEF calculation is −9 meV, i.e., still far away from the lowest-energy state of the total, non-rigid system.The two different relaxation phases can also nicely be noticed from Fig. 7(b).Accordingly, the oscillations of H die out quickly, and then the slower oscillations of the substrate "phonons" surface dominate, which decay only slowly.Of course, in a more realistic treatment direct phonon damping should be accounted for, presumably leading to a somewhat faster total energy loss also in the longtime regime.However, the emergence of two clearly different decay phases, a fast and a slow one, will survive also in this case.
In summary, for the present low-energy scenario the effects of vibration-phonon couplings are small on the timescale of vibrational relaxation (below 1 ps).Still small but subtle effects are found on the longer timescales of phonon motion (1 ps and later).When more energy is pumped into the system, e.g., during the relaxation of "hot atoms" on metals, the effects of surface phonons can be larger and the loss of excess energy in the H subsystem can be accelerated as argued elsewhere. 13,15

D. Other vibrations
A more efficient coupling of vibrations to phonons is also expected at the other end of the energy scale, namely, for low-energy vibrations which are comparable to phonon energies.In Fig. 8, we study the relaxation of one of the (x,y)-modes obtained from normal mode analysis of H-Pb 5 , when hydrogen is adsorbed in an fcc site.As mentioned earlier, in numerical practice NMA does not give exactly degenerate lateral modes, for systems with threefold symmetry and degenerate irreducible representations of the respective point group.Keeping the Pb lattice fixed, we get two slightly  non-degenerate (and localized) modes around 200 cm −1 , the lower one being at 195 cm −1 (24.2 meV).This mode is not a true (but close to an) eigenstate of the dynamical matrix, with a "normal mode vector" oriented in x-direction (with components (0.999, −0.050, 0.001).)This mode was selected to illustrate the coupling of a low-frequency adsorbate mode to EHPs, and to substrate phonons.For this purpose we started trajectories in the fcc minimum with a momentum corresponding to E vib and pointing in negative x-direction.In the calculation, timesteps of 0.1 fs were chosen.Both rigid and non-rigid surface models were considered, the latter by allowing the top four Pb layers to move.From Fig. 8(a) we see that the vibrational decay proceeds on the timescale of a few hundred fs also for the xmode.In fact, according to the simplified formula (5), lifetimes should be the same for all hydrogen-only modes but this is not quite true in an actual AIMDEF calculation.In the present example, for the rigid surface, we find a lifetime of ∼367 fs.This is a larger difference from the value of 342 fs suggested by the model expression (5) than for the z-mode (where AIMDEF gave 343 fs).This is presumably due to the larger anharmonicity of the x-mode and its strong coupling to other modes for which AIMDEF accounts, but we should also say that the numerical error is larger as indicated by a larger drift of the total energy in the frictionless case (by 2.7 meV up to 367 fs, or about 10% of the vibrational quantum).
Figure 8(a) also demonstrates that surface atom motion has a larger impact on the dynamics compared to the z-mode, for the reasons outlined above: While on the timescale of a ps inclusion of substrate motion leads to only marginal differences in E tot for z-mode excitation, the difference is more substantial for the x-mode.This is corroborated by the analysis of the total kinetic energy in Fig. 8(b), showing an appreciable amount of phonon excitation already on the sub-ps timescale.In fact, after about 100-200 fs, the kinetic energy is dominated by phonon contributions rather than H atom motion alone, when substrate atom motion is included.
In Fig. 8(c), we demonstrate the impact of phonons on the velocity of the vibrating H atom, for its x-component.[This motion is also strongly coupled to the y-component of the moving H atom, not much to the z-component (both not shown).]For this quantity, up to about 100 fs the influence of substrate phonons is basically negligible, then becoming more and more important at longer times.
Finally, in Fig. 8(d) we show the x-, y-, and z-components of the velocity of the top-layer Pb atom to which the H atom is bound, for the non-rigid surface.Note that the v xcomponents of H and Pb(top) are strongly correlated, showing both the same period of about 180 fs.This is quite consistent with the harmonic vibrational period of the x-mode, of about 170 fs.To a somewhat smaller extent also the y-and z-motions of the Pb atom are excited by the vibrating H atom. Their periods of vibration are higher and smaller, respectively, than for x, for reasons which we do not analyze here any further.
Sufficient to say that in summary, the frictional dynamics of low-frequency adsorbate modes is more complicated than for the near-harmonic z-mode, and also more strongly affected by surface phonons.The latter is due to the greater similarity of frequencies of lateral H-vibrations and substrate phonons.All of this is in line with our expectations, but it also proves that a detailed understanding of low-frequency adsorbate modes requires a sophisticated dynamical treatment.

IV. SUMMARY AND CONCLUSIONS
Using density functional theory and AIMDEF, we have studied the adsorption and vibrational relaxation of H atoms on films of Pb(111).It is suggested that AIMDEF in the LDFA is an approximate, but well-defined, parameter-free, and powerful, method to treat vibrational relaxation when dominated by electronic friction.The AIMDEF approach also allows one to include surface phonons, anharmonicity, and modecouplings.It is a computationally expensive method, however, probably the only practical way at present when many-atom dynamics under the influence of electronic friction is of interest.Then, no potential energy surfaces can be easily constructed, and the handling of tensorial friction 40 (not speaking of explicit inclusion of excited potential energy surfaces and couplings between them), is impractical.
As an example of the application of AIMDEF, we focused on vibrational relaxation of an adsorbate (hydrogen) at metal surfaces.Specifically, we investigated if the dissipative vibrational dynamics are influenced by quantum size effects, for Pb films of increasing thickness.Since the latter are known to exhibit oscillating (with film thickness) electron densities in front of the surface, different vibrational lifetimes caused by different embedding densities for the adsorbed atom were expected.In fact, efficient vibrational quenching of H-Pb vibrations by vibration-electron hole pair coupling on the timescale of a few hundred fs is predicted by our theory.Also, oscillations in lifetimes and other properties of the adsorbate emerge from our work.Some of these may be not easily detectable in practice, however, larger contrasts are predicted for organic molecules when physisorbed on lead films, and presumably other surfaces as well.Work along these lines are in progress in our laboratories.

FIG. 2 .
FIG. 2. (a) DFT (PBE) energy E/N per Pb layer as a function of N, for films consisting of N Pb layers with (111) orientation (see inset, d is the interlayer distance).(b) The same for the quantity E(N) as defined in Eq. (3).The energy per Pb atom obtained for the bulk (see text) is shown as a horizontal line in both figures.

FIG. 4 .
FIG.4.AIMD (black curves) and AIMDEF calculations (blue curves) for H-Pb 5 , when H was in the fcc site and excited with one vibrational "quantum" of the z-mode with E vib = 89.65 meV, initially moving towards the surface.Upper panel: Total energies for the non-dissipative and dissipative cases, with dashed red curves indicating energies expected after t = τ vib and t → ∞, respectively (see text).The zero of total energy refers to the rigid surface with H in its equilibrium position, unexcited.The vibrational lifetime obtained from AIMDEF is indicated.Lower panel: Kinetic energy of the H atom for both cases.

1 FIG. 5 .
FIG.5.Correlation between vibrational lifetimes obtained either from AIMDEF calculations, or the simple "model" lifetimes as defined in Eq. (5), for H adsorbed in fcc sites of Pb N (N = 1 to 6, see numbers attached to the data points).For the AIMDEF calculations, short timesteps for which the drift in the dissipation-free case was small or negligible, were used (see text).

FIG. 7 .
FIG.7.AIMDEF calculations for H-Pb 5 , when H was in the fcc site and excited with one vibrational "quantum" of the z-mode with E vib = 89.65 meV, initially moving towards the surface.(a) Total energies for the rigid surface (blue, same as in Fig.4(a)) and the non-rigid surface (black), within the first ps.(b) (Total) Kinetic energies for the rigid surface (blue, as in Fig.4(a)) and the non-rigid surface (black), within the first ps.(c) Velocity (z-component only) of the H atom for the rigid (blue) and non-rigid surfaces (black), within the first ps.(d) Same as (a), but for a longer time interval, up to 10 ps.(e) Same as (b), for up to 10 ps.(f) Same as (c), for up to 10 ps and a smaller ordinate interval.In (a) and (d), the zero of total energy refers to the rigid surface with H in its equilibrium position, unexcited.
Pb 5 , when H was initially in the fcc site and excited with one vibrational "quantum" of the x-mode with E vib = 24.18meV, initially moving in x-direction.(a) Total energies for the rigid surface (blue) and the non-rigid surface (black).Dashed red curves indicate energies expected after t = τ vib and t → ∞, respectively, for the rigid surface.The zero of total energy refers to the rigid surface with H in its equilibrium position, unexcited.The vibrational lifetime obtained for the rigid-surface case is indicated as an arrow.(b) (Total) Kinetic energies for the rigid surface (blue) and the non-rigid surface (black).(c) Velocity (x-component only) of the H atom for the rigid (blue) and non-rigid surfaces (black).(d) Velocity components (x, y, and z, see attached letters) of the top layer Pb atom, for the non-rigid surface.