Electronic friction-based vibrational lifetimes of molecular adsorbates: Beyond the independent atom approximation

We assess the accuracy of vibrational damping rates of diatomic adsorbates on metal surfaces as calculated within the local-density friction approximation (LDFA). An atoms-in-molecules (AIM) type charge partitioning scheme accounts for intra-molecular contributions and overcomes the systematic underestimation of the non-adiabatic losses obtained within the prevalent independent atom approximation. The quantitative agreement obtained with theoretical and experimental benchmark data suggests the LDFA-AIM as an efficient and reliable approach to account for electronic dissipation in ab initio molecular dynamics simulations of surface chemical reactions.

We assess the accuracy of vibrational damping rates of diatomic adsorbates on metal surfaces as calculated within the local-density friction approximation (LDFA). An atoms-in-molecules (AIM) type charge partitioning scheme accounts for intra-molecular contributions and overcomes the systematic underestimation of the non-adiabatic losses obtained within the prevalent independent atom approximation. The quantitative agreement obtained with theoretical and experimental benchmark data suggests the LDFA-AIM as an efficient and reliable approach to account for electronic dissipation in ab initio molecular dynamics simulations of surface chemical reactions. A central challenge in energy and catalysis applications is to transfer energy specifically into those degrees of freedom that actually drive a desired surface chemical reaction -and to keep this energy in these degrees of freedom for a sufficiently long time. In this transfer of energy, losses due to electronic non-adiabaticity can be an important dissipation channel [1,2]. In aiming to assess this channel for systems of technological interest, predictive-quality calculations would be a valuable addition to experimental endeavors. Especially for chemical reactions at frequently employed metal substrates, however, a corresponding methodology has not yet been well established.
To date, most accurate solutions of the full nuclearelectron wave function are restricted to systems of the complexity level of gas-phase H 2 + [3]. In the limit of weak non-adiabaticity as pertinent to electron-hole (eh) pair excitations during adsorbate dynamics on metal surfaces, less rigorous approaches rely on mixed quantumclassical dynamics. The imposed computational burden nevertheless still restricts their practical use to simple metals and sub-picosecond time scales [4,5], to symmetric adsorbate trajectories [6,7], or to only qualitative accounts of the metal electronic structure [8,9]. Presently, it is thus only the concept of electronic friction [10][11][12][13] and its incorporation into classical molecular dynamics (MD) simulations on the Born-Oppenheimer potential energy surface (PES) V PES [14][15][16][17][18][19][20][21] that promises predictive-quality and material-specific trajectory calculations over an extended period of time.
Particularly the local-density friction approximation (LDFA) [16,22] and for molecular adsorbates an additional independent atom approximation (IAA) [16][17][18] provide a further decrease in computational cost. This has allowed for first accounts of electronic nonadiabaticity in large-scale MD simulations based on a first-principles and high-dimensional description of the underlying PES -either interpolated [16][17][18][19] or even evaluated on-the-fly within ab initio MD simulations [20,21]. However, due to the drastic simplifications introduced with the IAA, the validity of the LDFA formalism for molecular adsorbates per se has been controversially discussed [16,23,24]. By construction the IAA does not resolve the electronic structure of the interacting molecule-surface system and in particular the location of the molecular frontier orbitals in the surface band structure. It can thus, for instance, not reproduce the enhancement of friction coefficients close to the transition state of a molecular dissociation event on metal surfaces [23]. On the other hand, the necessity of an accurate description of such regions of enhanced friction for the overall non-adiabatic energy dissipation has been questioned, as the typically low velocities in these regions effectively suppress the contribution of the friction term within the dynamics [16].
Despite the success in recent applications, it thus remains elusive to which extent the limitations of the prevalent IAA carry over to actual observables. In this situation, the vibrational lifetimes of high-frequency adsorbate modes can provide a sensitive measure, as they are largely governed by energy dissipation in the electronic non-adiabatic channel [25][26][27]. Accurate experimental reference data then allows for a substantiated assessment of the quality of the non-adiabatic description. In this study we perform such an assessment, primarily focusing on the internal stretch mode of two systems which arXiv:1506.06877v1 [cond-mat.mtrl-sci] 23 Jun 2015 have been studied most extensively and conclusively by experiments: CO adsorbed on Cu(100) and Pt(111) [26]. Despite the largely different surface frontier orbital locations and concomitant hybridizations at the transition and noble metal surface, we find the LDFA-IAA to already exhibit a good qualitative performance with respect to experimental [28,29] and theoretical [30,31] benchmark data. Rather than an explicit account of the surface band structure, our analysis suggests missing intra-molecular contributions as reason for the remaining differences. Approximately incorporating such contributions through a numerically efficient atoms-in-molecules charge partitioning indeed yields consistent lifetimes for a range of diatomic adsorbate systems.
In friction theory, all non-adiabatic effects due to the excitation of eh pairs in the metal substrate are condensed into a single velocity-dependent dissipative force that augments the classical equations of motion Here, small latin and greek subscripts denote atoms and Cartesian degrees of freedom, respectively, and N is the total number of atoms of mass m i and position R i in the system. The fluctuating white noise force F iα (T ) becomes negligible at very low temperatures and vanishes exactly at 0 K. Every element of the friction tensor η is, in principle, a function of all nuclear coordinates. Within the spirit of weak coupling, the focus is usually on the diagonal contributions describing the electronic friction felt by each atom [13,27]. These atomic friction coefficients η iαjβ = η iα δ ij δ αβ can e.g. be calculated within the quasi-static regime building on time-dependent density-functional theory (DFT) as suggested by Persson and Hellsing [11,32]. While insightful and generally in good agreement with experiments, the accurate numerical evaluation of this approach is challenging in practice and has hitherto been limited to low-dimensional potentials describing the adsorbatemetal interaction [14,15,31,33,34]. This shifts interest to more effective schemes and there in particular to the local-density friction approximation. The LDFA introduces isotropic scalar atomic friction coefficients η LDFA These coefficients can be calculated very efficiently from the scattering properties of an atomic impurity embedded in a free electron gas (FEG) [10,12,35]. In order to relate this model to the actual motion of adsorbates, the embedding density of the FEG ρ emb is then chosen as the electron density value of the clean metal surface ρ surf at the position of the adsorbate atoms. In its application to molecular adsorbates this implies an independent atom approximation -which has been regarded and discussed as being inherently included in the LDFA formalism in this context [16,23,24]. As a result of the IAA, the employed atomic friction coefficients are insensitive to the molecular nature of the adsorbate, yet can be evaluated very efficiently. The only input variable is the clean surface electronic density, which, assuming a frozen surface, has to be calculated only once in advance.
We obtain the two-dimensional PES V PES (d CO , Z COM ) of an adsorbed CO molecule as a function of its bond length d CO and center-of-mass (COM) height Z COM above the frozen surface through DFT calculations. Each PES is supported by 442 data points that are calculated with the CASTEP plane-wave pseudopotential code [36] and subsequently interpolated using bivariate cubic splines. Electronic exchange and correlation (xc) is treated on the GGA level in terms of the PBE functional [37]. The metal surfaces are modeled by five layer slabs with a separating 20Å vacuum distance. We consider top-site adsorbed CO molecules on one side of the slab within a c(2 × 2) and ( • surface unit-cell on Cu(100) and Pt(111), respectively. In both cases, the molecular axis is perpendicular to the surface, with the C atom coordinated to the metal atom. At the employed computational settings (600 eV cut-off energy, ultrasoft pseudopotentials [38], (10 × 10 × 1) and (11 × 11 × 1) Monkhorst-Pack k-point grids [39] for Cu(100) and Pt(111), respectively) the PES data points are converged to < 5 meV. Our investigation is not affected by the well-known CO adsorption puzzle and the concomitant wrong absolute depth of the adsorption well [40]. This is confirmed by essentially identical lifetimes we obtain when using PESs generated with a van der Waals-xc-functional [41] that leads to a stabilization of the top site [42].
Vibrational lifetimes τ are extracted from classical MD simulations on the interpolated PESs by numerically solving Eq. (1). Within the LDFA the atomic friction coefficients η LDFA i are calculated from the scattering phase shifts of the Kohn-Sham orbitals at the Fermi momentum δ F l = δ l (k F ) for an atomic impurity embedded in a FEG of density ρ emb,i [10,12,35] where ρ IAA emb,i = ρ surf (R i ) within the IAA as described before. Assuming a constant energy dissipation rate and thus an exponential decay of the vibrational energy E vib , the lifetime τ can be extracted from the simulations by a logarithmic fit of E vib versus time t. To initialize our simulations, we assign the adsorbate stretch-mode a projected kinetic energy of ω, where ω is the normal mode frequency. Higher initial kinetic energies up to 5 ω result in minute lifetime changes of less than 0.1 ps. Figure 1 shows the vibrational lifetimes that result from our simulations. We compare them to experimental values obtained from pump-probe spectroscopy [28,29] and theoretical values as published by Forsblom and Persson (FP) [30], and Krishna and Tully (KT) [31]. The latter two are both based on the orbital-dependent Persson/Hellsing expression mentioned above [11,32], albeit obtained from different derivations and relying on slightly different numerical treatment. For both systems the LDFA-IAA lifetimes agree fairly well with the theoretical values from FP and KT, as well as with experiment. With all numbers falling within one order of magnitude, the current data thus does not further support the harsh criticism the LDFA-IAA was faced with before [23]. Generally, they are instead consistent with the good LDFA-IAA performance reported in earlier studies on non-adiabatic energy losses of various ions scattered off metal surfaces [43] and on the vibrational damping of atoms on metal surfaces [11,32,44]. Conspicuously, however, in these studies on adsorbate atoms the agreement was even more quantitative and lacked the systematic underestimation (overestimation) of LDFA-IAA nonadiabatic energy losses (lifetimes) apparent in Fig. 1.
Rather than from a generally insufficient account of the electronic structure of the interacting adsorbate-surface system, this suggests that LDFA-IAA deficiencies arise particularly in the treatment of molecular adsorbates as isolated adatoms. In this respect -within the underlying atomic embedding model -systematic shortcomings lie in the complete neglect of both adsorbate-substrate as well as intra-molecular contributions to ρ emb . Considering only the density of the clean surface, ρ emb is systematically underestimated; consistent with the underestimated energy losses observed in Fig. 1. Even more, by construction a thus defined ρ emb is also blind to dynamical changes of intra-molecular bond distances and characters. Figure 2 illustrates this for the vibrational motion of CO at Cu(100). The LDFA-IAA friction coefficient of the carbon atom decreases with smaller C-O bond lengths d CO and thus larger distances Z C of the C atom from the metal surface. This is due to the fact that the IAA only accounts for changes of ρ surf along the vibrational coordinate. The LDFA-IAA thus effectively decomposes the molecular vibration into two independent adatom vibrations. Yet, even without any intramolecular bonding there will be an increasing overlap of the O and C atomic densities with decreasing d CO . This intra-molecular contribution to ρ emb at the position of every constituent atom is thus completely missed in the LDFA-IAA approximation.
A straightforward way to account for such contributions is to perform a charge decomposition through a projection scheme like Hirshfeld's analysis [45]. For any given adsorbate configuration this provides the projected density of any adsorbate atom ρ H i , and corresponding sharing function w H i at any position R [46]. Within an atoms-in-molecules (AIM) picture we can then define as embedding density for atom i at position R i We thus consider as embedding density the full selfconsistently calculated density ρ SCF of the entire interacting adsorbate-surface system just without the contri-bution associated with atom i. This naturally contains density contributions from all other atoms in the system (both substrate and adsorbate) and consequently carries an implicit dependence ρ AIM emb,i ({R i =j }). Figure 2(b) contrasts the correspondingly obtained ρ AIM emb,C to ρ IAA emb,C during the CO vibrational motion. Intriguingly, the two quantities differ not only quantitatively, but even exhibit a reversed phase behavior. The account of the density contribution of the O atom in ρ AIM emb,C thus outweighs the influence of the clean surface density seen by the IAA: At the largest vertical heights of the carbon atom Z C , where ρ IAA emb,C is smallest, also d CO is smallest. Already in the simple diatomic this leads to such a large intra-molecular density contribution that the overall ρ AIM emb,C is largest. These new features also carry over to the corresponding friction coefficients, and in turn to the calculated lifetimes. As apparent from Fig. 1 the LDFA-AIM approach cures the systematic overestimation of lifetimes and yields values that are now en par with the FP calculations and experiments. Similar findings have been obtained when applying the LDFA-AIM to two further systems [47] for which reference lifetimes have been reported from orbital-dependent theories, CN on Pt(111) [30] and H 2 on Ru(0001) [17]. With τ AIM = 0.9 ps vs. τ FP = 2.4 ps and τ AIM = 210 fs vs. τ [17] = 215 fs, respectively, LDFA-AIM yields in both cases lifetimes that are fully consistent with the reference numbers, yet at a fraction of the numerical cost of the orbital-dependent theories.
Adjusting ρ emb to also take into account influences from atoms other than only the clean metal surface hence seems to provide a simple, but effective correction for a molecular treatment within the LDFA. The idea underlying Eq. (3) is thereby similar to the subtraction of a free atom density as suggested in the context of vibrational damping of adatoms [44]. However, employing Hirshfeld sharing functions [45] offers the advantage that the embedding densities are guaranteed to be within physically well defined boundaries 0 ≤ ρ AIM emb,i ≤ ρ SCF (R i ). This way, our proposed scheme also preserves the molecular dissociation limit by construction: At large bond distances, the respective friction coefficients smoothly go over into friction coefficients virtually identical to the ones obtained for independent atoms. This is nicely illustrated by essentially identical vibrational lifetimes obtained within the IAA and AIM for the H 2 on Ru(0001) system, where the individual atoms are separated by about 2.7Å on the surface. Furthermore, a Hirshfeld analysis typically requires only a minute computational effort compared to achieving self-consistency for electronic energies and forces (or compared to the multiple self-consistent calculations required to obtain derivatives from finite differences in the FP or KT approaches). The LDFA-AIM scheme proposed here can thus be easily carried out at every time step of ab initio MD simulations, allowing surface motion to be explicitly taken into ac-count. In this respect, Eq. (3) defines an embedding density not only for adsorbate but also for bulk atoms. Friction coefficients derived therefrom could thus, in principle, as well be invoked to evaluate energy losses due to electron-phonon coupling in the bulk.
In conclusion, we have shown that the vibrational damping of high-frequency adsorbate modes on metal surfaces can be added to the list of non-adiabatic phenomena that are reasonably well described by means of electronic friction. The complete neglect of intramolecular effects in the prevalent LDFA-IAA approximation is thereby likely to underestimate the electronic energy dissipation. The here presented atoms-inmolecules alternative instead accounts for them approximately through a charge partitioning scheme. As it thus effectively treats the molecular electrons as part of the metallic substrate, we expect the AIM friction concept to generally rather overestimate non-adiabatic energy losses and to perform best for chemisorbed adsorbates at close distances to the surface. Being a direct descendant of the LDFA, our scheme is, of course, also unlikely to overcome fundamental limitations that come with its heritage. As such it is unlikely to properly capture the strong enhancement of friction coefficients directly at the transition state leading to molecular dissociation. However, as important as it may be, this actual dissociation event only constitutes a fraction of the dynamics that is relevant to chemical surface reactions. Other important aspects like a vibrational pre-excitation or an ensuing hot adatom motion may dominate the overall (non-adiabatic) energy dissipation [20,48], yet take place over much longer time-scales. They can thus only be assessed with a numerically highly efficient method like the LDFA. In this respect, our results further consolidate the trust in LDFA-based results for these important long-term events.