Testing Electronic Friction Models: Vibrational De-excitation in Scattering of H2 and D2 from Cu(111).

At present, molecular dynamics with electronic friction (MDEF) is the workhorse model to go beyond the Born-Oppenheimer approximation in modeling dynamics of molecules at metal surfaces. Concomitant friction coefficients can be calculated with either the local density friction approximation (LDFA) or orbital-dependent friction (ODF), which, unlike LDFA, accounts for anisotropy while relying on other approximations. Due to the computational cost of ODF, extensive high-dimensional MDEF trajectory calculations of experimentally measurable observables have hitherto only been performed based on LDFA. We overcome this limitation with a continuous neural-network-based representation. In our first application to the scattering of vibrationally excited H2 and D2 from Cu(111), we predict up to three times higher vibrational de-excitation probabilities with ODF than with LDFA. These results indicate that anisotropic electronic friction can be important for specific molecular observables. Future experiments can test for this "fingerprint" of different approximations underlying state-of-the-art MDEF.

The motion of atomic and molecular adsorbates on metal surfaces underlies every elementary reaction step in heterogeneous catalysis. Due to the absence of an energy gap between valence and conduction band electrons, these motions can result in the excitation of electron–hole pairs (EHPs) and thus violate the Born–Oppenheimer approximation.[1−3] A growing number of experiments points to the importance of this nonadiabatic energy loss channel.[4] On the other hand, the development of suitable theoretical models to account for these nonadiabatic effects is still an ongoing process.[5−9] For systems with weak nonadiabatic coupling, molecular dynamics with electronic friction (MDEF)[10] is currently the most popular approach.[11] MDEF relies on a potential energy surface (PES), nowadays typically obtained from density functional theory (DFT),[3] and accounts for the effects of the EHPs on the motion of the nuclei by electronic friction coefficients.[10] One state-of-the-art technique for calculating these coefficients as functions of the adsorbate positions relies on mapping to an atom-in-jellium model for which only the surface electron density is considered (local density friction approximation, LDFA[12,13]). Alternatively, the electronic states of the molecule–surface system can be taken into account (orbital-dependent friction, ODF[10,14]). For the inelastic scattering of H atoms from Au(111), millions of MDEF trajectories based on a high-dimensional PES[15] and LDFA have recently been demonstrated to yield accurate scattering probabilities in excellent agreement with experimental data.[16]

The situation is quite different for molecules. Due to its combination with the independent atom approximation, the LDFA completely neglects any molecular effects.[12] ODF on the other hand accounts for the anisotropic tensorial character of friction coefficients on corrugated metal surfaces and along adsorbate-internal bonds,[17−19] which is why ODF has been argued to be “theoretically” more accurate for (diatomic) molecules.[20] However, this discussion[12,20,21] has still remained inconclusive because an evaluation of ODF comes at very high computational costs. Consequently, extensive MDEF trajectory calculations for molecules including all relevant degrees of freedom (DOF) can be easily performed with LDFA,[11,12] whereas only two molecular DOF have so far been included for ODF.[22] The very recent on-the-fly evaluation of ODF within ab initio molecular dynamics by Maurer et al.[19] is an important step, but the less than 20 calculated trajectories make direct validation via molecular beam experiments impossible. Modeling the nonadiabatic contribution to vibrational lifetimes of molecules adsorbed on metal surfaces on the other hand does not require such extensive statistical averaging.[23,24] The most recent implementations of LDFA and ODF both yield results that agree with experimental data within the error bars.[25,26] Furthermore, Novko et al. have shown recently in this context[27] that the numerical evaluation of friction tensors within ODF[19,26] effectively includes potentially spurious electronic memory effects with unclear consequences for MDEF.[28] Given this situation, theoretical understanding and modeling relying on MDEF faces an important question: Is the molecular anisotropy as described by ODF important for any observables that can be validated by high-precision molecular beam experiments like for atoms?[16]

In this work, we provide an answer to this question using H2 and D2 on Cu(111). For this system, weak nonadiabatic coupling and static surface approximation are well justified,[29−31] and an accurate DFT-based PES relying on the semiempirically constructed “SRP” exchange–correlation functional is available.[32,33] We construct a six-dimensional neural-network-based continuous representation of ODF that allows us to perform extensive MDEF trajectory calculations on equal footing with LDFA. While dissociative sticking probabilities are hardly affected in general and by the type of electronic friction coefficients used, we find vibrational de-excitation probabilities to be a “fingerprint” that can be used to distinguish and validate LDFA and ODF in future experiments.

Quasi-classical trajectory calculations[32] within MDEF rely on a generalized Langevin equation[10]where i,j indicate atoms and α,β Cartesian coordinates. Atomic masses and positions are denoted by m and R, respectively. For a H2 or D2 molecule on a static surface, the total number of moving atoms N is two, resulting in six DOF. In addition to the forces from the PES , which yield the adiabatic dynamics, nonadiabatic effects on the nuclear dynamics originate from electronic friction forces Fiαfric(R) and thermal white noise , respectively. In this work, V(R) is mainly taken to be the static surface PES based on the SRP48 exchange–correlation functional from ref (33), but we also compare with the PW91-PES from earlier work.[34] The friction forces are linear in nuclear velocities Ṙ and are in general given by a symmetric 6 × 6 friction tensor η(R), which consists of 21 independent elements each depending on six nuclear coordinates. These coordinates can be Cartesian R = (R1,R2) or expressed in the center-of-mass centered spherical coordinate system R = (X,Y,Zd,θ,ϕ), which is commonly used for diatomics and described by Figure A.

Figure 1

(A) Molecular coordinate system denoting the center of mass positions (X,Y,Z), bond length d, as well as spherical orientation (θ,ϕ). (B) Top view of a reference configuration with , Y = 0, and θ0 = ϕ0 = 90° from the minimum-energy reaction path for dissociative chemisorption over the bridge site,[34] where a denotes the surface lattice constant. Cu atoms in the first, second, and third layer are depicted by increasing transparence. Note that X,Y,Z = 0 corresponds to the position of a Cu atom in the surface plane (top site).

Within ODF, these 21 friction coefficients are obtained according to a Fermi golden rule-like expression resulting from time-dependent perturbation theory, which can be written in the quasi-static limit as[10,24,26,35]The electron–phonon matrix elements g describe the nonadiabatic coupling between two electronic states of the molecule at the metal surface with band indices a and b at wave vector k due to the motion of (adsorbate) atom i along direction α. In general, the ODF tensor can have different diagonal elements even for the same atom. This anisotropy yields very different friction forces when the atoms move (with the same velocity) in different directions. Its generally nonzero off-diagonal elements are responsible for coupling the motion in different directions and between both atoms in a way that is not accounted for by the PES. In particular, this can lead to a strong damping of the molecular stretch vibration of a diatomic molecule and thus a pronounced molecular anisotropy.[19,20,22]

In order to use the so-calculated ηiODF(R) in MDEF trajectory calculations of generic experimentally measurable observables, a continuous representation of this 6 × 6 tensor is required that can be evaluated at low computational cost. We have designed such a representation based on a symmetry-adapted neural network fit that is briefly described in the Supporting Information and will be discussed extensively in a forthcoming publication.

Within LDFA, friction coefficients for hydrogen atoms ηH(ρ) are obtained from a spherical atom-in-jellium model with background density ρ, which is solved via DFT at the level of the local density[36] or generalized gradient approximation.[37] Mapping of the actual surface problem is accomplished by taking the electron density of the bare surface (without the molecule) at each atom’s position ρ(R) as the background density of the jellium.[12] This independent-atom approximation (IAA) results in electronic friction coefficients that are isotropic for each atom and depends on its own three coordinates alone. In Cartesian coordinates, only diagonal elements of the friction tensor are nonzeroA continuous representation of ηiLDFA(R) for extensive MDEF trajectory calculations can be easily constructed.[12,38]

Going beyond the IAA within LDFA is possible, for example, by determining the background electron density using an atoms-in-molecules technique (LDFA-AIM).[25] However, this approach does not lift the isotropy and, as detailed in the Supporting Information, cannot be applied to H2 and D2 molecules. The other way round, isotropic friction can be constructed from ODF by neglecting the coupling between different directions and atoms plus averaging the remaining (generally anisotropic) friction over different directionsThis ODF-iso allows one to disentangle the influence of anisotropy from the very different electronic structure inherent to ODF and LDFA.

Figure A–C shows η, η, and η, respectively, as obtained from eqs –4 along the minimum-energy reaction path for dissociative chemisorption over the bridge site, as depicted in Figures D and 1B. We focus here on these three particular friction coefficients in order to compare with the earlier two-dimensional ODF calculations.[20,22] The agreement is quite good except for some differences close to the transition state for ηODF. As the molecule approaches the surface, each model yields increasing friction for the six diagonal elements of the friction tensor, and the absolute values of the off-diagonal elements increase likewise in the case of ODF. Furthermore, ODF directly reflects the strong rearrangement of (Kohn–Sham) orbitals when approaching the dissociation barrier by significantly higher friction along the molecular bond (and thus reaction) coordinate, resulting in ηODF ≈ 3ηLDFA at the transition state, in agreement with earlier work.[20] For the observables calculated below, friction beyond the dissociation barrier is not relevant. Quite remarkably, ηODF-iso (ηODF-iso) and ηLDFA (ηLDFA) are almost identical up to the transition state and thus much more alike than what has originally been found for the diffusion of H atoms on Pd(100).[17,39]

Figure 2

(A–C) η, η, and η in the molecular coordinate system (see Figure A), respectively, along the minimum-energy reaction path for dissociative chemisorption over the bridge site (see Figure B) as depicted in (D) together with the corresponding two-dimensional PES cut. The blue, red, and purple lines indicate the continuous representation from this work for ODF, LDFA, and ODF-iso, respectively, as obtained from eqs –4. Blue (red) dots show the ODF (LDFA) results from previous work of Luntz et al.,[20,22] and the reaction coordinate is defined in the same way as in that work. The barrier and thus the transition state for dissociation are located at the vertical gray line (i.e. 0 Å) in (A–C) and indicated by the empty circle in (D). Negative numbers up to the transition state denote the approach from the gas phase (i.e., decreasing heights Z above the surface).

In order to study the effect of the different friction models on actual experimental observables, we perform MDEF calculations according to the quasi-classical trajectory method.[28] In view of the short interaction time of the molecules with the Cu(111) surface during all simulated trajectories, we neglect the fluctuating forces in eq .

Figure A,B shows the results for the dissociative chemisorption probability S0 for both H2 and D2 molecular beams based on the SRP48-PES, respectively. Due to the construction of the latter,[32,33] already the adiabatic calculations yield good agreement with the experimental data.[40,41] Inclusion of electronic friction reduces S0, leading to even better agreement with the experimental data, in particular, at high incidence energies. The reduction is strongest for ODF and weaker for LDFA and ODF-iso, which are very similar to each other. It can be rationalized by the differences of the friction models for the friction η along the reaction coordinate close to the dissociation barrier (see Figure B). This effect of ηODF on S0 for H2 and D2 on Cu(111) has not been reported for two-dimensional ODF calculations.[22] Consequently, a high-dimensional treatment of ODF in MDEF, on an equal footing with LDFA,[21] is important. However, the overall small effect of electronic friction on S0 makes this not an optimal observable for experimental validation of the different friction models.

Figure 3

Calculated reaction probabilities S0 based on the SRP48-PES for dissociative chemisorption of (A) H2 and (B) D2 molecular beams as a function of average normal incidence energy ⟨Ei⟩ for the indicated nozzle temperatures Tnozzle in comparison to experimental data (brown unfilled squares) from refs (40) and (41), respectively. The calculations are adiabatic (filled black squares) or employ the LDFA (red plusses), ODF (blue circles), and ODF-iso (purple crosses) models for the electronic friction coefficients.

Instead, we have identified vibrational de-excitation probabilities Ptransition to yield a clearly distinguishable difference between LDFA and ODF. On the basis of 50000 MDEF trajectories, we calculate Ptransition as a function of incidence energy Ei from the scattered trajectories by a conventional binning procedure. The concomitant average gain in translational energy ⟨ΔEtrans⟩ is calculated from the final center-of-mass velocities. As detailed in the Supporting Information, the error bars reflect the error due to statistical sampling of the initial conditions. Only by employing our continuous representation to compute a large amount of trajectories were we able to reduce these errors so that the different electronic friction models can be distinguished. We focus the discussion on de-excitation from vibrational state ν = 2, J = 1 (2) to ν = 1, J = 1 (2) for H2 (D2), respectively, as shown in Figure A,C (B,D). Unlike for other vibrational transitions,[42] for this transition, we obtain results that are not only qualitatively but even almost quantitatively identical to corresponding results obtained with the PW91-PES (see the Supporting Information).

Figure 4

Vibrational de-excitation probabilities Ptransition (lower row) and concomitant average gain in translational energy ⟨ΔEtrans⟩ (upper row) as a function of normal incidence energy Ei for state-to-state scattering using the SRP48-PES. Panels (A,C) [(B,C)] are for the transition from the rovibrational state ν = 2, J = 1 [2] to ν = 1, J = 1 [2] for H2 [D2]. Shown are results from adiabatic calculations (filled black squares), as well as those including electronic friction according to the LDFA (red pluses), ODF (blue circles), or ODF-iso (purple crosses) models. The error bars indicate the error due to statistical sampling as described in detail in the Supporting Information.

At low incidence energies, with increasing Ei, more and more molecules come close enough to the surface so that the curvature of the PES and electronic friction lead to an increase of Ptransition. Both effects are additive and result in de-excitation probabilities that are up to 6 (2) times larger for H2 (Figure C) and up to 3 (2) times larger for D2 with ODF (LDFA), respectively (Figure D). At high incidence energies, the dissociation channel (see Figure A,B) becomes more effective, which is why Ptransition decreases again in all cases. For the adiabatic simulations on the static surface, ⟨ΔEtrans⟩ is equal to the rovibrational energy loss of one vibrational quantum and thus by about √2 larger for H2 than that for D2 (Figure A,B). Electronic friction reduces the energy gain. The reduction is almost twice as large for ODF compared to that for LDFA. The fact that it does not very strongly depend on Ei for the energy range considered here suggests that it is dominated by η and thus directly reflects the differences observed along the minimum-energy path depicted in Figure A–D. Consequently, when comparing MDEF with other nonadiabatic models[43−45] for vibrational de-excitation, our results suggest that it is crucial to also take into account whether the friction coefficients include any molecular anisotropy. Unfortunately, because molecular beam experiments for this system have hitherto focused on rovibrational excitation rather than de-excitation,[46−48] experimental verification of this effect is still pending.

Although ODF-iso inherits the spurious memory effects as well as going beyond the independent atom approximation from ODF, quite remarkably, for Ei > 15 (20) kJ/mol for H2 (D2), we obtain results with ODF-iso that are almost identical to those with LDFA. That means that (at least in this energy range) these do not affect the dynamics and the molecular anisotropy is the most important difference. For lower incidence energies, scattering over the top site has been found to dominate vibrational de-excitation from ν = 1 in adiabatic calculations.[42] Indeed, for top sites, ηODF-iso is rather different from ηLDFA so that the difference in electronic structure inherent to LDFA and ODF also becomes visible in the dynamics in this case. If the molecular anisotropy could be experimentally validated, it would greatly encourage future theoretical work to develop extensions to LDFA that might be able to (at least approximately) account for it.

In summary, we have obtained different observables for H2 and D2 on Cu(111) from extensive MDEF trajectory calculations for the first time using full-dimensional friction tensors based on both LDFA and ODF. The molecular anisotropy as described by ODF and absent from LDFA leads to strongly enhanced friction for motion along the molecular axis when the molecules are close to the surface. The dissociative sticking probability is almost negligibly reduced compared to adiabatic simulations. The effect is slightly stronger for ODF compared to LDFA and improves the agreement with experimental data in both cases. For the state-to-state scattering of vibrationally excited molecules (from ν = 2, J = 1 (2) to ν = 1, J = 1 (2) for H2 (D2)), we predict up to six (two) times larger vibrational de-excitation probabilities with ODF (LDFA) compared to adiabatic simulations. Remarkably, isotropicalization of ODF yields results almost identical to LDFA for incidence energies larger than 15 (20) kJ/mol for H2 (D2). The predicted differences between the vibrational de-excitation probabilities are a “fingerprint” of the molecular anisotropy as described by ODF. Recently suggested techniques to prepare H2 molecular beams with 1 ≤ ν ≤ 4[49] should allow testing for this fingerprint. This would provide unprecedented insight into the accuracy of state-of-the-art electronic friction models for molecules and allow analysis of the importance of concomitant approximations.

Computational Details

In eq , we calculate the electron–phonon matrix elements from the change of the Kohn–Sham potential with respect to nuclear coordinate R, which is obtained from density functional perturbation theory (DFPT)[50] employing the PW91[51] exchange–correlation functional as implemented in the Quantum Espresso package.[52] Surfaces are modeled by 2 × 2 Cu(111) slabs with four layers and 10 Å of vacuum. A plane-wave cutoff energy of 816 eV is used, together with ONCV pseudopotentials[53] from the SG15[54] library and an 18 × 18 k-point grid. These settings reproduce the PW91-PES from ref (34) up to a few meV. They also enable an accurate evaluation of the sum over electronic states in eq at the Fermi level using an equivalent Gaussian envelope technique to broaden the δ-function with a width of 0.6 eV as suggested in ref (26). We note that this implies the possible presence of spurious electronic memory effects, as argued in ref (27).

The neural network fits for the 21 independent elements of ηODF(R) are based on ∼30000 ODF coefficients obtained from DFT calculations on the same seven lateral sites that have been used to construct the SRP48-PES.[33] For LDFA, we extract the background electron density ρ(R) from a DFT calculation with the same computational setup as described above. Employing the functional form for ηH(ρ) suggested in ref (38), we then construct three-dimensional neural network interpolations for ηH(ρ(R)) based on the symmetry-adapted coordinates[55,56] in order to obtain a continuous representation of ηLDFA(R).