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

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 Wodtke et al. (2004); Nienhaus (2002). 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 nuclear-electron wave function are restricted to systems of the complexity level of gas-phase ${\text{H}_{2}}^{+}$ Abedi et al. (2010). 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 quantum-classical dynamics. The imposed computational burden nevertheless still restricts their practical use to simple metals and sub-picosecond time scales Lindenblatt and Pehlke (2006); Grotemeyer and Pehlke (2014), to symmetric adsorbate trajectories Timmer and Kratzer (2009); Meyer and Reuter (2011), or to only qualitative accounts of the metal electronic structure Shenvi et al. (2009a, b). Presently, it is thus only the concept of electronic friction Echenique et al. (1981); Hellsing and Persson (1984); Echenique et al. (1986); Head-Gordon and Tully (1995) and its incorporation into classical molecular dynamics (MD) simulations on the Born-Oppenheimer potential energy surface (PES) $V_{\text{PES}}$ Trail et al. (2003); Luntz and Persson (2005); Juaristi et al. (2008); Füchsel et al. (2011); Martin-Gondre et al. (2012); Füchsel et al. (2013); Blanco-Rey et al. (2014); Saalfrank et al. (2014) that promises predictive-quality and material-specific trajectory calculations over an extended period of time.

Particularly the local-density friction approximation (LDFA) Li and Wahnström (1992); Juaristi et al. (2008) and for molecular adsorbates an additional independent atom approximation (IAA) Juaristi et al. (2008); Füchsel et al. (2011); Martin-Gondre et al. (2012) provide a further decrease in computational cost. This has allowed for first accounts of electronic non-adiabaticity in large-scale MD simulations based on a first-principles and high-dimensional description of the underlying PES – either interpolated Juaristi et al. (2008); Füchsel et al. (2011); Martin-Gondre et al. (2012); Füchsel et al. (2013) or even evaluated on-the-fly within ab initio MD simulations Blanco-Rey et al. (2014); Saalfrank et al. (2014). 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 Juaristi et al. (2008); Luntz et al. (2009); Juaristi et al. (2009). 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 Luntz et al. (2009). 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 Juaristi et al. (2008).

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 Tully et al. (1993); Arnolds (2011); Saalfrank (2006). 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 have been studied most extensively and conclusively by experiments: CO adsorbed on Cu(100) and Pt(111) Arnolds (2011). 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 Beckerle et al. (1991); Morin et al. (1992) and theoretical Forsblom and Persson (2007); Krishna and Tully (2006) 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 $\mathbf{R}_{i}$ in the system. The fluctuating white noise force $\mathcal{F}_{i\alpha}(T)$ becomes negligible at very low temperatures and vanishes exactly at $0\,\text{K}$. Every element of the friction tensor $\bm{\eta}$ 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 Head-Gordon and Tully (1995); Saalfrank (2006).

These atomic friction coefficients $\eta_{i\alpha j\beta}=\eta_{i\alpha}\delta_{ij}\delta_{\alpha\beta}$ can e.g. be calculated within the quasi-static regime building on time-dependent density-functional theory (DFT) as suggested by Persson and Hellsing Hellsing and Persson (1984); Persson and Hellsing (1982). 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 adsorbate-metal interaction Lorente and Persson (2000); Trail et al. (2003); Luntz and Persson (2005); Krishna and Tully (2006); Luntz et al. (2006). 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 $\eta_{i}^{\text{LDFA}}$, such that $\eta_{i\alpha j\beta}^{\text{LDFA}}=\eta_{i}^{\text{LDFA}}\delta_{ij}\delta_{\alpha\beta}$. These coefficients can be calculated very efficiently from the scattering properties of an atomic impurity embedded in a free electron gas (FEG) Echenique et al. (1981); Puska and Nieminen (1983); Echenique et al. (1986). In order to relate this model to the actual motion of adsorbates, the embedding density of the FEG $\rho_{\text{emb}}$ is then chosen as the electron density value of the clean metal surface $\rho_{\text{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 Juaristi et al. (2008); Luntz et al. (2009); Juaristi et al. (2009). 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_{\text{PES}}(d_{\text{CO}},Z_{\text{COM}})$ of an adsorbed CO molecule as a function of its bond length $d_{\text{CO}}$ and center-of-mass (COM) height $Z_{\text{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 Clark et al. (2005) and subsequently interpolated using bivariate cubic splines. Electronic exchange and correlation (xc) is treated on the GGA level in terms of the PBE functional Perdew et al. (1996); *Perdew_PhysRevLett_1997. The metal surfaces are modeled by five layer slabs with a separating $20\,\text{\AA }$ vacuum distance. We consider top-site adsorbed CO molecules on one side of the slab within a $c(2\times 2)$ and $(\sqrt{3}\times\sqrt{3})\text{R}30^{\circ}$ 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\,\text{eV}$ cut-off energy, ultrasoft pseudopotentials Vanderbilt (1990), $(10\times 10\times 1)$ and $(11\times 11\times 1)$ Monkhorst-Pack $\mathbf{k}$-point grids Monkhorst and Pack (1976) for Cu(100) and Pt(111), respectively) the PES data points are converged to $<5\,\text{meV}$. Our investigation is not affected by the well-known CO adsorption puzzle and the concomitant wrong absolute depth of the adsorption well Feibelman et al. (2001). This is confirmed by essentially identical lifetimes we obtain when using PESs generated with a van der Waals-xc-functional Dion et al. (2004) that leads to a stabilization of the top site Lazić et al. (2010).

Vibrational lifetimes $\tau$ are extracted from classical MD simulations on the interpolated PESs by numerically solving Eq. (1). Within the LDFA the atomic friction coefficients $\eta^{\text{LDFA}}_{i}$ are calculated from the scattering phase shifts of the Kohn-Sham orbitals at the Fermi momentum $\delta_{l}^{\text{F}}=\delta_{l}(k_{\text{F}})$ for an atomic impurity embedded in a FEG of density $\rho_{\text{emb},i}$ Echenique et al. (1981); Puska and Nieminen (1983); Echenique et al. (1986)

where $\rho_{\text{emb},i}^{\text{IAA}}=\rho_{\text{surf}}(\mathbf{R}_{i})$ within the IAA as described before. Assuming a constant energy dissipation rate and thus an exponential decay of the vibrational energy $E_{\text{vib}}$, the lifetime $\tau$ can be extracted from the simulations by a logarithmic fit of $E_{\text{vib}}$ versus time $t$. To initialize our simulations, we assign the adsorbate stretch-mode a projected kinetic energy of $\hbar\omega$, where $\omega$ is the normal mode frequency. Higher initial kinetic energies up to $5\hbar\omega$ 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 Morin et al. (1992); Beckerle et al. (1991) and theoretical values as published by Forsblom and Persson (FP) Forsblom and Persson (2007), and Krishna and Tully (KT) Krishna and Tully (2006). The latter two are both based on the orbital-dependent Persson/Hellsing expression mentioned above Persson and Hellsing (1982); Hellsing and Persson (1984), 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 Luntz et al. (2009). 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 Winter et al. (2003) and on the vibrational damping of atoms on metal surfaces Persson and Hellsing (1982); Hellsing and Persson (1984); Tremblay et al. (2010). Conspicuously, however, in these studies on adsorbate atoms the agreement was even more quantitative and lacked the systematic underestimation (overestimation) of LDFA-IAA non-adiabatic 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 $\rho_{\text{emb}}$. Considering only the density of the clean surface, $\rho_{\text{emb}}$ is systematically underestimated; consistent with the underestimated energy losses observed in Fig. 1. Even more, by construction a thus defined $\rho_{\text{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_{\text{CO}}$ and thus larger distances $Z_{\text{C}}$ of the C atom from the metal surface. This is due to the fact that the IAA only accounts for changes of $\rho_{\text{surf}}$ along the vibrational coordinate. The LDFA-IAA thus effectively decomposes the molecular vibration into two independent adatom vibrations. Yet, even without any intra-molecular bonding there will be an increasing overlap of the O and C atomic densities with decreasing $d_{\text{CO}}$. This intra-molecular contribution to $\rho_{\text{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 Hirshfeld (1977). For any given adsorbate configuration this provides the projected density of any adsorbate atom $\rho_{i}^{\text{H}}$, and corresponding sharing function $w_{i}^{\text{H}}$ at any position $\mathbf{R}$ McNellis et al. (2009). Within an atoms-in-molecules (AIM) picture we can then define as embedding density for atom $i$ at position $\mathbf{R}_{i}$

We thus consider as embedding density the full self-consistently calculated density $\rho_{\text{SCF}}$ of the entire interacting adsorbate-surface system just without the contribution 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 $\rho^{\text{AIM}}_{\text{emb},i}(\{\mathbf{R}_{i\neq j}\})$. Figure 2(b) contrasts the correspondingly obtained $\rho^{\text{AIM}}_{\text{emb},\text{C}}$ to $\rho^{\text{IAA}}_{\text{emb},\text{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 $\rho^{\text{AIM}}_{\text{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_{\text{C}}$, where $\rho^{\text{IAA}}_{\text{emb,C}}$ is smallest, also $d_{\text{CO}}$ is smallest. Already in the simple diatomic this leads to such a large intra-molecular density contribution that the overall $\rho^{\text{AIM}}_{\text{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 ¹¹1Consistent with the theoretical work by FP Forsblom and Persson (2007), we consider perpendicular CN adsorption with the C atom coordinated to the metal on a top site within a $(2\times 2)$ surface unit cell. The H₂ on Ru(0001) system is modeled by analogy with the equilibrium configuration mentioned by Füchsel et al.. Füchsel et al. (2013), i.e., using a $(2\times 2)$ surface unit cell with two H atoms adsorbed in adjacent fcc hollow sites. For both systems, we use a $(9\times 9\times 1)$ Monkhorst-Pack $\mathbf{k}$-point grid Monkhorst and Pack (1976). All remaining computational details are the same as outlined for CO on Cu(100) and Pt(111), respectively. for which reference lifetimes have been reported from orbital-dependent theories, CN on Pt(111) Forsblom and Persson (2007) and H₂ on Ru(0001) Füchsel et al. (2011). With $\tau^{\text{AIM}}=0.9\,\text{ps}$ vs. $\tau^{\text{FP}}=2.4\,\text{ps}$ and $\tau^{\text{AIM}}=210\,\text{fs}$ vs. $\tau^{\text{\cite[cite]{\@@bibref{Authors Phrase1YearPhrase2}{Fuechsel_PhysChemChemPhys_2011}{\@@citephrase{(}}{\@@citephrase{)}}}}}=215\,\text{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 $\rho_{\text{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 Tremblay et al. (2010). However, employing Hirshfeld sharing functions Hirshfeld (1977) offers the advantage that the embedding densities are guaranteed to be within physically well defined boundaries $0\leq\rho^{\text{AIM}}_{\text{emb},i}\leq\rho_{\text{SCF}}(\mathbf{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₂ on Ru(0001) system, where the individual atoms are separated by about $2.7\,\text{\AA }$ 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 account. 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 intra-molecular effects in the prevalent LDFA-IAA approximation is thereby likely to underestimate the electronic energy dissipation. The here presented atoms-in-molecules 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 Blanco-Rey et al. (2014); Meyer and Reuter (2014), 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.