Determining the Effect of Hot Electron Dissipation on Molecular Scattering Experiments at Metal Surfaces.

Nonadiabatic effects that arise from the concerted motion of electrons and atoms at comparable energy and time scales are omnipresent in thermal and light-driven chemistry at metal surfaces. Excited (hot) electrons can measurably affect molecule-metal reactions by contributing to state-dependent reaction probabilities. Vibrational state-to-state scattering of NO on Au(111) has been one of the most studied examples in this regard, providing a testing ground for developing various nonadiabatic theories. This system is often cited as the prime example for the failure of electronic friction theory, a very efficient model accounting for dissipative forces on metal-adsorbed molecules due to the creation of hot electrons in the metal. However, the exact failings compared to experiment and their origin from theory are not established for any system because dynamic properties are affected by many compounding simulation errors of which the quality of nonadiabatic treatment is just one. We use a high-dimensional machine learning representation of electronic structure theory to minimize errors that arise from quantum chemistry. This allows us to perform a comprehensive quantitative analysis of the performance of nonadiabatic molecular dynamics in describing vibrational state-to-state scattering of NO on Au(111) and compare directly to adiabatic results. We find that electronic friction theory accurately predicts elastic and single-quantum energy loss but underestimates multiquantum energy loss and overestimates molecular trapping at high vibrational excitation. Our analysis reveals that multiquantum energy loss can potentially be remedied within friction theory whereas the overestimation of trapping constitutes a genuine breakdown of electronic friction theory. Addressing this overestimation for dynamic processes in catalysis and surface chemistry will likely require more sophisticated theories.

Introduction

The Born–Oppenheimer approximation gives rise to the notion of a single potential energy surface (PES) that governs chemical dynamics. Despite its great success, the breakdown near electronic degeneracies is well-known and corresponding nonadiabatic effects have profound implications in various fields such as photochemistry and single molecule electronics.[1] This is particularly true in elementary chemical reactions at metal surfaces, which are of fundamental and practical importance in heterogeneous catalysis, as there is virtually no energy threshold for electronic excitation in metals.[2] As a result, the gaseous species in the vicinity of a metal surface can easily dissipate their energy not only by exciting lattice vibrations but also through electron–hole pair excitations (EHPs).[3] Indeed, there has been growing experimental evidence of such nonadiabatic effects in surface chemistry[4] from quantum-state-resolved molecular beam scattering experiments,[5] chemicurrent measurements,[6,7] and ultrafast spectroscopy,[8] providing valuable benchmark data for testing first-principles theories of nonadiabatic gas-surface interactions.[9] However, a predictive quantitation of how nonadiabatic effects contribute to measurable dynamic properties remains elusive.

The continuum of electronic states in metallic systems is a daunting challenge to the first-principles simulation of nonadiabatic gas-surface scattering dynamics.[1] While a full-dimensional quantum treatment is at present unfeasible, several pragmatic mixed quantum-classical dynamics (MQCD) methods have been developed,[10−14] two of which stand out for their practical feasibility when combined with ab initio electronic structure theory. The first is the independent electron surface hopping (IESH) method[10,15,16] which is based on the popular surface hopping trajectory method,[17] that characterizes nonadiabatic effects via probabilistic electronic transitions between electronic states.[18] The IESH method describes the hopping of independent electrons with a Newns–Anderson Hamiltonian parametrized with density functional theory (DFT) data. However, it is difficult to determine excited states and their nonadiabatic couplings from first-principles for metallic systems and several ad hoc approximations are required in the parametrization.[15] An alternative is the molecular dynamics (MD) with electronic friction (MDEF) method which assumes weak nonadiabaticity[19,20] (eq ). Herein, electronic degrees of freedom (DOFs) are described via a frictional damping force that represents the nonadiabatic linear response of electrons to the motion of adsorbate atoms.[21] This force acts on the atoms in addition to the force arising from the PES (eq ).In eq , M and R are the mass and position of a nucleus, respectively, i and j are nuclear coordinates, V is the PES, Λ is the electronic friction tensor (EFT), and is a force associated with random white noise from the bath of electrons. In practice, the MDEF method is always further approximated by imposing the Markov approximation of instantaneous response in the constant coupling limit,[20] where some pragmatic assumptions are made in how the friction tensor is calculated from DFT in that limit.[12] Examples of approximations include the local density friction approximation (LDFA), which allows for an efficient calculation of scalar isotropic friction from the electron density of the metal[22,24−30] and the more realistic orbital-dependent friction (ODF),[12,23] which is calculated from Kohn–Sham DFT via time-dependent perturbation theory to provide a better description of the mode-selective nature of nonadiabatic molecule-metal energy transfer.[24−28]

All existing practical methods to study nonadiabatic dynamics introduce significant approximations which need to be scrutinized against experiment. This is of course mixed with the underlying errors of adiabatic PES itself in describing the energy landscape. Unfortunately, there is little quantitative data for realistic systems that describes under which conditions exactly which approximation breaks down and how this depends on the molecule-metal coupling strength. In other words, how do we know when the weak-coupling limit is satisfied and when the MDEF method is reliable for a particular system? While this question was partially addressed by Dou and Subotnik for simple model systems,[29] here we provide quantitative insights for state-to-state scattering of NO from Au(111), which has been considered a representative strong-coupling showcase for the breakdown of electronic friction theory.[30,31]

Over many years, Wodtke and co-workers have collected ample state-to-state experimental data that reveals unambiguous nonadiabatic characteristics of this benchmark gas-surface process for a wide range of scattering conditions,[5,9,31−36] stimulating many different theoretical studies.[10,15,37,38] While the aforementioned IESH model has partially accounted for the multiquantum vibrational relaxation/excitation of NO scattered from Au(111),[16,33] its predictions on the translational energy dependence of vibrational relaxation probabilities[36] and some other observables[31] were qualitatively inconsistent with experimental findings. These discrepancies have been largely attributed to the “too-soft” and “too-corrugated” adiabatic PES within the diabatic model Hamiltonian expressed by simple pairwise potentials.[36] Using the same adiabatic PES, earlier MDEF calculations have qualitatively failed to describe the nonadiabatic dynamics for this system, especially the vibrational excitation of NO(vi = 0) scattering from Au(111) in the vibrational ground state.[33] Interestingly, a reduced-dimensional quantum-mechanical version of electronic friction has been able to yield broad vibrational state distributions compatible with experiment.[38]

The recent emergence of high-dimensional machine-learning-based PESs has enabled the reduction of interpolation errors stemming from the fitting of the PES.[39] In this work, we use a recently developed embedded atom neural network (EANN) based adiabatic PES of this system with high-fidelity involving realistic surface DOFs.[40] With a more accurate description of repulsive NO–Au(111) interaction and potential energy topography shaped by the tight transition state, this PES has enabled much more adiabatic vibrational energy transfer than previously expected and provides a qualitatively correct translational energy dependence of vibrational inelasticity.[40] This allows us to largely reduce errors in the description of the adiabatic PES and to focus on scrutinizing the quality of nonadiabatic description at a level that was not possible before. By combining this highly accurate PES with a faithful multidimensional EANN representation of the full-rank ODF EFT derived by time-dependent perturbation theory,[41] in the present work, we study systematically the influence of EF on the state-to-state scattering dynamics of NO from Au(111) (see Figure for a schematic system definition). Impressively, incorporating both molecular and surface DOFs, the MDEF(ODF) model allows a quantitatively correct description of the single quantum vibrational relaxation dynamics of NO(vi = 3) and NO(vi = 2) and their dependence on translational energy. The MDEF(LDFA) model, by contrast, has little effect on the dynamics beyond the adiabatic description. We provide a detailed analysis to rationalize this. By comparing against experimental data with systematically increasing incidence vibrational energy, we further pinpoint the energetic regime in which MDEF and, specifically, Markovian MDEF break down.

Figure 1

Schematic plot of NO on Au(111) showing internal (X, Y, Z, r, θ, φ) coordinates of NO molecules and relaxation of surface atom.

Results and Discussion

General Performance of Electronic Friction Methods

State-to-state quantum scattering of NO from Au(111) is a perfect system to scrutinize the performance of MDEF due to the availability of experimental data for various different incidence conditions. Figure shows final vibrational state (vf) distributions for NO scattering prepared in initial vibrational states (vi) of 2, 3, 11, and 16.[31,32,35,36] All experiments were performed at incidence energies ranging between 0.52 and 1.08 eV. Initial rotation states are chosen to closely match that employed in experiment. The scattering events are expected to be dominated by a single bounce due to the narrow angular distribution observed in experiment.[5,36] We first compare experiments with IESH and MDEF simulations performed on a previously published PES.[31] The IESH simulations, show strong overestimation of the elastic scattering contributions for vi = 11 and vi = 16. MDEF simulations on the same PES correctly predict the elastic scattering populations but deliver vibrational state distributions that only lose 2–4 quanta on average with almost no population at lower final vibrational states. The failure of both methods is evidence of an inaccurate PES.[44]

Figure 2

Experimental (Exp, golden histogram bars, respective references) vs BOMD vs LDFA vs ODF final vibrational state distributions for (a) vi = 2, ji = 2, Ei = 0.640 eV,[32] (b) vi = 3, ji = 0, Ei = 1.08 eV,[36] (c) vi = 11, ji = 0, Ei = 0.950 eV,[35] and (d) vi = 16, ji = 0, Ei = 0.520 eV.[35] Only single bounce trajectories are included in all models including the reference IESH and MDEF data;[31] additionally all data has been renormalized to the experimental limits (see Figure S7 for further clarification). The referenced MDEF friction coefficients have been calculated with a different approximate approach that is not comparable to ODF or LDFA. Lines are drawn between markers for visual clarity.

Figure further shows the results of adiabatic scattering simulations (labeled Born–Oppenheimer MD, BOMD), and MDEF simulations with LDFA and ODF using the new high-dimensional EANN PES and representation of EFT (the construction of which is described in the Computational Methods section and the Supporting Information (SI)). We find that the angular scattering distribution predicted by the EANN model is in good agreement with the experiment, with a low population of multibounce events (see Figure S6). Recent theoretical studies with other less accurate PESs[36,42] showed the necessity of excluding multibounce trajectories to acquire more realistic results. It is not necessary to exclude multibounce trajectories with our EANN PES, we do so anyway for the vibrational state distributions shown in the main text to ensure that vibrational energy loss only arises due to scattering and not due to equilibration on the surface. For completeness, Figures –5 are reproduced in the SI with multibounce events included (Figures S7–S10). As expected, the distributions do not differ significantly. We further exclude from our analysis any vibrational states whose populations were not measured in the corresponding experimental work[43] (see Figures S7–S10 for an analysis featuring all simulated final states).

Figure 5

(a) Mass-weighted internal stretch friction excitation spectrum for the adsorption structure defined in the SI. The gray dashed line is a 0.6 eV Gaussian curve used when evaluating the internal stretch element value (gray arrow), while the horizontal purple dotted line depicts the element multiplied by 4 as employed when scaling ODF. ODF with and without an anisotropic scaling of 4 vs LDFA with an isotropic scaling of 4 is shown for (b) vi = 11 and (c) vi = 16;[35] conditions are otherwise the same as those in Figure . Single bounce selected only. Lines are drawn between markers for visual clarity.

Our results correctly predict that NO scattering is highly vibrationally inelastic featuring the loss of one or more vibrational quanta leading to very broad final state distributions of highly vibrationally excited NO molecules that are clearly dominated by multiquantum vibrational energy loss. Our tensorial orbital-based description of EF is a significant improvement over the adiabatic description and local-density description of EF for single vibrational quantum loss and elastic scattering for low initial vibrational states. Impressively, in all conditions, the new MDEF(ODF) results agree better with experiment than the reference IESH and MDEF data in terms of the broadness/shape and peak positions of final state distributions, despite some remaining discrepancies with experiment. This implies that previous studies on the failure of EF theory might have conflated PES artifacts with failings to describe nonadiabatic effects. Indeed, the current more accurate PES allows us to isolate the role of nonadiabatic effects by analyzing the remaining discrepancies of MDEF(ODF) simulations with experiment. We identify two major discrepancies, namely, (i) MDEF(ODF) fails to improve on the adiabatic description and continues to underestimate multiquantum vibrational energy loss for low and high incidence vibrational energies (e.g., the underestimation of vf = 1 population for vi = 3 shown in Figure. b) and (ii) MDEF(ODF) overestimates the trapping probability (see Figure ). The failure to reproduce the vf = 1 population for vi = 3 or make any significant improvement over the adiabatic description is particularly telling for the inability of MDEF to predict multiquantum vibrational energy loss, which we analyze in detail further below. This is further emphasized in Figure b,c where MDEF(ODF) only really modifies the populations of elastic and single quantum loss channels. We also note that the MDEF(ODF) description slightly overestimates the elastic population for vi = 11 but underestimates it for vi = 16. The adiabatic results (previously discussed by some of us[40]) capture a significant portion of multiquantum loss for vi = 11 and vi = 16, where the dominant vibrational scattering channels are 3 and 5 quanta loss respectively, though this is around half of what is predicted by experiment.

Figure 6

(a) Experimentally determined and BOMD, LDFA, and ODF predicted trapping probabilities for vi = 2 (ji = 2) over a range of incidence energies.[44] Also shown are BOMD and ODF with a rescaled potential surface, [RS]. The black dashed line represents an experimentally determined fit.[44] The absolute trapped populations for (b) vi = 11, ji = 0, Ei = 0.950 eV and (c) vi = 16, ji = 0, Ei = 0.520 eV are also shown. Lighter bars are recorded with the rescaled potential energy surface.

The behavior exhibited in Figure a,b for low initial vibrational states largely holds for a range of incidence translational energies as demonstrated in Figure . MDEF(ODF) performs well across the range of incidence energies for elastic and single-quantum inelastic scattering, capturing best the translational energy dependence of vibrational inelasticity among all theoretical models. However, it only slightly improves upon the adiabatic description of two-quantum inelastic scattering, with both only managing to capture the general qualitative trend that high incidence energies lead to more two-quantum loss. Interestingly, MDEF(ODF) performs slightly worse at lower incidence energies (Figure a,c). We can attribute this to an artifact of the underlying PES (Figure S8) which will be discussed below. MDEF(LDFA) fails to significantly improve upon the adiabatic description at any incidence energy, strongly suggesting an account of the molecular nature of the impinging NO and the directional dependence and intermode coupling of EF is required to describe this system. The adiabatic results capture the qualitative incidence energy dependence, though significantly overestimate the importance of the elastic channel. In the following, we will investigate the origin of the failures of our EF simulations.

Figure 3

Experimental (Exp, relevant references) vs BOMD vs LDFA vs ODF branching ratios (for population of final state, P(vf)) for vi = 2 (ji = 2)[32] and vi = 3 (ji = 0).[36] Each plot is labeled with an arrow from the initial state to the final state. Only single bounce trajectories are included in all models including the reference IESH and MDEF data.[36] BOMD predicts no vi = 2 to vf = 1 population at low incidence energies so is omitted. Lines are drawn between markers for visual clarity.

Failure to Predict Multiquantum Loss for Low Initial Vibrational States

We can further analyze the origin of the failure to capture multiquantum loss by breaking down the vi = 3, Ei = 0.950 eV experiment with respect to initial molecular orientation as the experiment shows a significant orientation dependence. In the experiments, the NO molecules are aligned with the nitrogen (N-down) or with the oxygen (O-down) pointing toward the surface. Figure shows that molecules that start with an N-down orientation experience significantly more inelastic energy loss than molecules that start with O-down. MDEF(ODF) simulations succeed in reproducing this effect qualitatively, whereas adiabatic and MDEF(LDFA) models do not.

Figure 4

(a) Depiction of N (blue) down and O (red) down orientations. Experimental vs BOMD vs LDFA vs ODF final vibrational state distributions for vi = 3, ji = 0, Ei = 0.950 eV[45] with (b) all, (c) only nitrogen down, and (d) only oxygen down orientations included. Only single bounce trajectories are included. Lines are drawn between predicted distributions for visual clarity. Experimental results are shown as histogram bars. Note that in the experimental work, the vf = 0 population is not explicitly measured but rather assumed to be one-half of the corresponding vf = 1 population in all cases.

A closer look reveals that the final state distribution of O-down scattering is particularly well reproduced by MDEF(ODF), albeit with some level of underestimation of multiquantum loss. Experiments reveal that N-down dynamics undergo more single and double vibrational quanta loss; MDEF(ODF) reproduces the former well but not the latter while overestimating the elastic contribution. It appears that the inability to describe sufficient multiquantum loss from N-down dynamics is the major source of discrepancy between MDEF(ODF) and experiment for vi = 3 scattering shown in Figures and 3.

We note that, in agreement with other theoretical results[16] for this system, there is a strong dynamical steering effect, such that in our results an initial orientation does not guarantee a similar orientation when colliding with the surface. N-Down collision geometries are energetically preferred, such that even O-down initially orientated trajectories are predominately steered into an N-down collision geometry (Figure S11a). On average, we can see initial N-down orientations correspond to closer approaches to the surface and higher elongation of the N–O bond (Figure S11b). Already an EF model as simple as LDFA tells us that nonadiabatic effects increase exponentially as molecules come closer to the surface and this leads to stronger nonadiabatic molecule–metal coupling. From our previous work on H2 on Ag(111), we know that bond elongation leads to drastic increases in nonadiabatic coupling along the intramolecular stretch mode.[24] We further see that trapping predominately occurs for trajectories with a close surface approach of <2 Å (Figure S12). We reasonably expect the trapping behavior for vi = 3 to be similar to that experimentally determined for vi = 2, which is expected to be negligible for Ei = 0.950 eV.[44]

At close surface approach, several effects could contribute to the underestimation of multiquantum loss. First, molecules could be trapped that should in fact scatter with substantial energy loss. We discuss this effect further below. A second effect could lie in the current calculation of the ODF EFTs which only considers excitations that are both (i) first-order (single-electron excitation) and (ii) interband (ie. transitions that conserve momentum). It has been shown that phonon-assisted intraband excitations are the dominant contributor to the short vibrational lifetime of CO adsorbed on a Cu(100) surface,[46] though at a dense coverage of adsorbate molecules which is not the case here. A possible neglect of intraband contributions will particularly affect lower vibrational states and lower translational energies. We can test this effect by increasing the size of the unit cell, which effectively increases the number of electron–hole-pair excitations in the Brillouin zone that are accessible by momentum-conserving excitations. The effect is explored in detail in Figure S4. When changing unit cell size, the EFT elements do not change drastically over a range of energies when no broadening is used nor does the broadened EFT significantly change. This suggests that intraband contributions are sufficiently accounted for in our description.

Lastly, practical MDEF simulations are always performed within the Markov approximation. The time-dependent motion of the adsorbate excites EHPs and the ensuing energy dissipation between these DOFs is dependent on the energy of the perturbing molecular motion and the density of states (DOS) of the substrate. Due to the Markov approximation, here we assume that it is independent of both. In the following, we will explore how this affects our results for high initial vibrational states.

Failure to Predict Multiquantum Loss for High Initial Vibrational States

In the case of highly vibrationally excited molecules (vi = 11 and vi = 16), the failure of MDEF(ODF) to predict multiquantum loss upon scattering is even more evident (see Figure c,d). A decreased sensitivity of the final vibrational state distribution to both incidence energy and molecular orientation was observed in experiment for vi = 11 and further for vi = 16.[35] This was suggested to be due to the driving force of vibrational relaxation becoming very large for high vibrational states.[35]

To understand how the failure of MDEF(ODF) occurs, we must recall how the EFT is calculated. The EFT element Λ associated with adsorbate motion in directions Ri and R in first order perturbation theory can be expressed as[12,19]

The EFT for adsorbate motion with a frequency associated with energy ϵ = ℏω is calculated by summing over the product of relevant nonadiabatic coupling matrix elements over all possible excitations between effective single particle Kohn–Sham states ε and ε with respective occupation factors f(ε) and f(ε). By assuming a constant DOS around the Fermi level, Head-Gordon and Tully[20] were able to invoke a constant coupling assumption, which leads to the Markovian expression of MDEF. In practice, the EFT is evaluated at the Fermi level (zero excitation energy), replacing the delta function with a smearing function of 0.6 eV finite width.[12] Lifting the Markov approximation would lead to the inclusion of memory effects, which corresponds to the inclusion of EF at higher perturbing energies due to the modulation of particle velocity along the scattering trajectory. The inclusion of memory effects in the electronic friction force leads to a response between EHPs and adsorbate DOFs that draws contributions from the full EF spectrum. We expect that the importance of EF at higher perturbing energies will increase for higher incidence vibrational energies.

The friction excitation spectrum (Figure a) shows that for small broadening values and perturbing energies other than zero, the coupling may reach values several times higher than the Markovian EFT value indicated by the arrow. High vibrational states of NO lead to strong velocity oscillations and the excitation of EHPs further away from the Fermi level. As can be seen in the spectrum, the constant coupling approximation is not a good one in the case of NO on Au(111). No full memory-dependent implementation of MDEF exists at the moment, but it is clear that the inclusion of memory effects will lead to an increase in the magnitude of electronic friction and the Markovian EFT corresponds to a lower bound. We can investigate the potential effects of memory by scaling the ODF EFT internal stretch element by 4 (see the SI for methodology), which approximately represents the difference between the broadened Markovian friction value and the highest friction values present in the spectrum at nonzero frequencies. In this manner, we are studying close to an upper bound of the effects of memory on the strength of EF forces. Further in the SI, we demonstrate that the internal stretch element governs the nonadiabatic vibrational distribution with very little difference between an isotropic scaling of the whole ODF EFT or just the internal stretch element (Figure S10).

Though the individual state populations described by the scaled MDEF(ODF) model for vi = 11 and vi = 16, presented in Figure b,c, show deviations in relative contribution from experiment of about 0.05–0.10, the overall vibrational distribution is well represented. A similar scaling for the vi = 3 case also shows an improvement of the final state distribution (see Figure S9). Notably, scaling of the LDFA EFT does not provide any improvement on the results presented in Figure , which again confirms that the anisotropic nature of EF must be accounted for. Scaling the EFT is of course a primitive approach to account for the nonadiabatic coupling that arises from the excitation of EHPs at various energies present within this system, we instead use it for qualitative analysis of the shortcomings of MDEF(ODF) and its comparison to MDEF(LDFA). A more advanced representation of dynamical energy loss by including the memory-dependence of EF would likely provide a more accurate final state distribution. If confirmed, this would mean that the energy loss of gas-surface scattering, even when it involves high vibrational excitation, can be represented without having to abandon the conceptual basis of electronic friction theory. This stands in contrast to previous experimental and theoretical works,[30] which assign the inability to correctly represent the final state distribution to a direct metal-molecule electron transfer due to the presence of a transient anionic state of NO. This would, however, correspond to a clear departure from the weak coupling limit described by MDEF toward multistate dynamics as described by trajectory surface hopping techniques such as IESH.

Failure to Capture the Trapping Probability

Finally, we wish to discuss the failure of MDEF to describe the trapping probability of NO scattering from Au(111). Trapping probabilities are overestimated from those that were experimentally determined for vi = 2; assuming the vf = 3 (excitation) and vf = 0 (double quantum loss) channels are negligible, the former has been experimentally recorded to be very small.[32] We employ the same methodology to calculate the model predicted trapping probabilities in Figure a. Indeed, the predicted vf = 0 and vf = 3 populations are very small (see Figure S7a) so that the trapping probability is very close to the absolute trapped population. Figure a shows the systematic overestimation of trapping over a range of incidence energies for adiabatic dynamics, with the application of either friction model not changing the picture significantly. The two possible origins of this are a potential overestimation of the adsorption well in the EANN PES rooted in the intrinsic errors of the semilocal PW91 functional or the presence of strong nonadiabatic effects such as transient ion formation that leads to a dynamical change in the energy landscape.[16,47] We expect the former to affect low incidence energy scattering more strongly and the latter to affect high incidence energy scattering more strongly.

Indeed, we find that the EANN PES for NO on Au(111), while substantially more accurate than previous models, does still overestimate the molecule–surface attraction. In order to identify if the overestimation of trapping is due to the energy landscape or due to the description of nonadiabatic effects, we have added a repulsive contribution to the PES to reduce the adsorption energy to the experimentally observed value of 0.24 eV[48] (see the SI for details). This results in a trapping probability at low incidence energies that is very close to the experimentally observed value (see ODF[RS] in Figure. ). In the SI, we show that the adjusted PES only has minor effects on final vibrational state distributions, with the exception of improving the agreement of MDEF(ODF) with experiment for low incidence energies (see Figure S8). Figures S9 and S10 show that the final state distributions for vi = 3, 11, and 16 at moderate incidence energies originally shown in Figures and 4 are not significantly affected, leaving our previous conclusions on multiquantum energy loss unaffected. This also suggests that our main conclusion in this work would not be significantly altered using different density functionals that may yield different adsorption well depts or barrier heights. Nevertheless, as can be seen in Figure a, the trapping probability as predicted by MDEF(ODF[RS]) remains too high at high incidence energies.

Figure b,c shows the absolute trapped population for both vi = 11 and vi = 16. The trapping probability has not been measured experimentally for high incidence vibrational states, though it is expected to be relatively insignificant, even at very low incidence energies (Ei = 0.05 eV).[5,49] On this basis, considering the high vibrational state and the moderate to high incidence energies employed, we should expect negligible trapped population. This is not the case for the adiabatic results and the application of EF which leads to an even larger trapping probability. Application of the rescaled PES significantly reduces the trapping (Figure b,c) but does not nullify it.

Contrary to our results, the low trapping probability at high vi has been correctly predicted by Shenvi et al.,[16] where IESH predicts a trapping probability far lower than their BOMD results and far lower than what our present BOMD and MDEF results suggest. The lowering of the trapping probability compared to adiabatic results has been related to transient nonadiabatic metal-to-molecule charge transfer which leads to an enhancement of vibration-to-translation energy transfer.[16] This is opposite to the effect that electronic friction has, which dominantly describes vibrational dissipation into EHPs, enhancing molecular trapping rather than reducing it.

Conclusion

We have presented a systematic analysis of the performance of state-of-the-art nonadiabatic simulation methods in describing hot-electron effects in vibrational state-to-state-scattering of NO on Au(111). To understand how nonadiabatic effects contribute to measurable dynamic reaction probabilities, we need to be able to isolate the role of nonadiabatic effects from other contributing factors. This is made possible with a newly created high-dimensional machine-learning-based potential energy landscape that resolves artifacts of previous PES models. While the model still overestimates the probability of trapping, readjustment of the PES to match experimental trapping at low incidence energies shows that the quantities of interest, namely final vibrational state distributions upon scattering, are not strongly affected by this. Using a rotationally covariant machine-learning model, we construct a high-dimensional model of ODF electronic friction calculated from DFT. We find that MDEF(ODF) provides excellent agreement with experiment for elastic scattering of various initial vibrational states and for single quantum vibrational energy loss of low initial vibrational states (vi = 2 and vi = 3), as well as the orientation dependence of vibrational state distribution, but otherwise underestimates multiquanta vibrational energy loss. Particularly in the case of high initial vibrational states such as vi = 11 and vi = 16, the width and the shape of the final state distribution is well described, but the average number of lost vibrational quanta is much smaller than in experiment. As we apply the Markov approximation, the high-lying EHPs of 1.5 eV and beyond that are excited by such high vibrational states are not included in our EF description. By analysis of the friction spectrum and rescaling of the friction tensor to account for this shortcoming, we find that we can reproduce the overall population distribution of inelastic scattering, albeit at a small remaining underestimation of the proportion of scattering outcomes with low vibrational energy. This finding is surprising as it suggests that a full account of memory effects within EF theory could potentially extend the remit of MDEF to describe the final vibrational state distributions of highly excited molecules without having to resort to the explicit inclusion of strong nonadiabatic effects on the energy landscape.

However, memory-dependent MDEF would likely not resolve the second failure of our MDEF results, namely, the significant overestimation of trapping at high incidence translational energies. Whereas experimental trapping probabilities at low incidence energies can be correctly predicted once the energy landscape matches the experimental binding energy, the same is not the case at high incidence energies. In agreement with previous literature, we conclude that this failure is likely due to the neglect of strong nonadiabatic effects that arise from transient charge- and/or spin-transfer between metal and molecule yielding a change in effective energy landscape, which goes beyond nonadiabatic energy dissipation that is described via electronic friction theory. To resolve this failure of electronic friction theory, stochastic surface hopping methods such as IESH will likely be required. However, such methods need to be integrated with more realistic first-principles determined charge transfer states,[50] which remains very challenging for periodic metallic systems. To fully understand the case of NO on Au(111), additional experiments that provide insight into the trapping probability of highly vibrationally excited molecules at high incidence energies will be useful in the future.

We believe that the here presented approach and simulation results, together with the extensive experimental data by Wodtke and co-workers, provide a firm baseline for the future development of more reliable and efficient nonadiabatic dynamics methods that will open the door to study nonadiabatic effects in thermal and photoelectrochemical reactions at catalyst surfaces in the future.

Computational Methods

In ref (40), we have performed spin-polarized DFT calculations for the NO + Au(111) system using the Vienna Ab initio Simulation Package[51,52] with the PW91 functional.[53] The Au(111) surface was represented by a four-layer slab model in a 3 × 3 unit cell with the top two layers movable. The Brillouin zone was sampled by a 4 × 4 × 1 Gamma-centered -point grid. A total of 2722 points with both energies and forces were collected mainly from direct dynamics trajectories to represent the adiabatic PES. More details can be found in ref (40). In this work, additionally, the 6 × 6 ODF electronic friction tensor (EFT) was evaluated for 1647 (+ 1052) training (+ test) points using our implementation within the all-electron numerical atomic orbital code FHI-Aims.[12,54] In the SI, we provide further details on numerical settings and evidence of the robustness of the friction tensor evaluation with respect to these settings (see Figure S3). Quasi-classical trajectory calculations were performed using a modified VENUS code.[55]

We employ the recently developed embedded atom neural network (EANN) approach to represent the scalar potential energy[56] and EFT[41] surfaces for NO on Au(111). In the EANN model, the potential energy is expressed as the sum of the embedded atomic energy, each of which is a complex function of the embedded density of the corresponding central atom. Different from potential energy, the EFT is covariant with respect to rotation (or reflection) of the molecule and permutation of identical atoms in the molecule, which is much more difficult to learn by neural networks. For the NO + Au(111) system, we start with a 6 × 3 first-order derivative matrix (corresponds to three neurons in the NN output layer) and a 6 × 6 second-order derivative matrix in terms of the partial derivatives of neural network outputs with respect to atomic Cartesian coordinates of the NO molecule. Multiplying the first- and second-order derivative matrices with their own transpose, respectively, yields two 6 × 6 matrices that naturally guarantee the rotational covariance and positive semidefiniteness of the EFT. The summation of the two 6 × 6 matrices is employed to approximate the EFT to account for additional symmetry of the EFT with respect to a symmetric mirror. More technical details can be found in our recent work[41,56−58] and in the SI.