Hot-electron effects during reactive scattering of H2 from Ag(111): assessing the sensitivity to initial conditions, coupling magnitude, and electronic temperature

Low-lying electronic excitations in metals, so-called hot electrons, efficiently mediate molecule-metal energy transfer and contribute to energy loss during molecular reactions at surfaces. They furthermore play an important role in plasmon-driven chemistry. Electronic friction represents a simple and effective concept to model hot electron-induced energy loss under ambient conditions. Different methods exist that vary in their description of magnitude, coordinate and directional dependence of friction during reactive molecular scattering at metal surfaces. Using molecular dynamics simulations with electronic friction, we systematically study the effect of hot electrons on measurable state-to-state scattering probabilities of molecular hydrogen from a (111) surface of silver. We assess the ability of ab initio electronic friction methods to accurately describe hot electron-mediated energy loss as a function of initial reaction conditions and electronic temperature. We furthermore find that dynamic scattering results and the ensuing energy loss are highly sensitive to the magnitude of electronic friction. Therefore, existing approximate models of electronic friction, which exhibit inherent uncertainties with respect to the magnitude of electronic friction, may not be applicable for a quantitative prediction of plasmon driven hot electron effects in their current state. We outline a development direction to potentially overcome these limitations.


Introduction
Plasmonic science 1 and hot-electron chemistry [2][3][4] are based on the idea of controlling the flow of energy between electronic excitations in metallic surfaces or nanostructures, 5 so-called hot electrons, and molecular adsorbate degrees of freedom, such as molecular vibrations. A high level of control of this molecule-surface energy flow can enable to facilitate chemical bond breaking events with light by funneling energy into the reaction coordinate or to impair a reaction by efficiently dissipating energy away from a reaction centre. 6,7 The underlying mechanism for this energy exchange is based on the efficient coupling of molecular adsorbate motion and metallic electrons due to the vanishing band-gap in metallic systems. 8,9 The consequence of this coupling is a deviation from thermal chemistry as described by the Born-Oppenheimer approximation and a failure of existing methods to describe chemical reaction dynamics. 10,11 The proposed technological applications of hot-electron effects in plasmonic sensors, 12 plasmon-enhanced catalysts, 13 and locally-enhanced cancer radiation treatment, 14 are driven by exciting fundamental experiments that provide quantitative evidence for hot electron effects. These include molecular beam scattering experiments that explicitly measure the state-to-state scattering probability and the energy loss profile of a gas-surface reaction, 15 but also the measurement of electric currents, 11,16 socalled chemicurrents, that occur due to hot electrons that are generated upon adsorption, desorption, or chemical transformation of molecules.
The ability to computationally predict the molecular details that underlie hot-electron-mediated processes and the intricate coupling of electrons and atomic motion is essential for the design of new devices and catalysts based on plasmonics. As such, it is imperative to develop first-principles computational simulation methods that correctly capture how molecular adsorbate motion generates hot electrons or how hot electrons affect measurable dynamical properties of chemical reactions. The intrinsic dimensionality of systems involving many hundreds of atoms and thousands of electrons make exact nonadiabatic quantum dynamics simulations close to unfeasible and therefore ask for efficient approximate mixed quantum-classical dynamics methods. 17,18 Such methods are typically based on stochastic surface hopping models based on chemical master equations, where hot electron effects on the reaction dynamics are modelled as explicit electronic transitions. 19,20 An alternative is the molecular dynamics with electronic friction (MDEF) method, [21][22][23] where the classical motion of nuclei follows a single potential energy surface and hot-electron effects are described as electronic friction forces that act on the adsorbate nuclei. These forces are governed by the so-called electronic friction tensor, 24 which arises from nonadiabatic coupling and the resulting equations of motion follow a stochastic Langevin model. The electronic friction tensor is in general frequency-dependent and anisotropic but is often approximated to be frequency independent and spatially isotropic, i.e., nonadiabatic relaxation rates are equal in all spatial directions for rigid atomic and molecular motion.
The most common is the Local Density Friction Approximation (LDFA), 25,26 where nonadiabatic coupling is expressed as a simple function of the metallic electron density within which the adsorbate atom is embedded. This model is spatially isotropic for rigid molecular motion, but does account for some anisotropic effects along intramolecular degrees of freedom. 26,27 Due to its simplicity and ease of computation, LDFA-based MDEF is a particularly appealing approach that, in the past, has been applied to describe hot-electron-mediated dynamics in ambient 28 and laser-driven conditions. 29,30 However, especially for molecular scattering and light-driven plasmonics, the underlying neglect of frequency dependence and the role of molecular anisotropy of electronic friction forces and their effects on measurable observables have been the reason for exciting debates in recent literature. 24,[31][32][33] In recent work, we have achieved a description of electronic friction that captures the directional and the explicit vibrational mode dependence and tensorial nature of electronic friction using an efficient implementation of first-order time-dependent perturbation theory (TDPT) based on Density Functional Theory (DFT). 24,31 On the examples of vibrational energy loss of metal-adsorbed diatomics 31 and the dissociative chemisorption of molecular hydrogen on Ag(111), 34 we were able to show that the overall magnitude of friction is, on average, similar to the LDFA prediction, but the directional and vibrational mode dependences and their effect on the dynamic energy loss are different. Before us, Spiering and Meyer have reached a similar conclusion for reactive scattering of H2 on Cu(111), 35 where mode-dependent electronic friction forces as described by TDPT provide a measurable signature of energy loss during state-tostate scattering of vibrationally excited molecules. We have recently developed a general method to construct an analytical representation of the electronic friction tensor as a function of the adsorbate coordinates that enabled us to simulate dynamic scattering probabilities and to discuss the subtle differences in electronic friction that give rise to the observation that hot electron effects are much more pronounced during reactive scattering of H2 on Cu(111) than on Ag(111). 36 In this work, we use this analytical representation of electronic friction for H2 on Ag(111) further to systematically assess the validity and robustness of the current frequency-independent MDEF method based on TDPT electronic friction. We do this by studying the sensitivity of dynamic observables such as the vibrational excitation/deexcitation probability during molecular scattering and the integrated energy loss along a specific degree of freedom as a function of various initial conditions and the type of electronic friction description. The former includes the effects of the initial and final vibrational excited states and the electronic temperature of the substrate. The latter includes a comparison between LDFA and TDPT and a discussion of their limitations, and sensitivity to intrinsic numerical parameters. We find that the energy loss predicted by LDFA and TDPT shows different trends with respect to the parameters we studied and that dynamic observables are particularly sensitive to the absolute magnitude of tensorial friction. The intrinsic sensitivity of friction to numerical parameters such as the broadening width will represent a particular challenge to quantitative predictions of hot-electron mediated reaction dynamics at high electronic temperatures that correspond to light-driven conditions. We therefore stress the importance to overcome this limitation introduced by assuming frequency independence of the friction tensor.
This will enable to further study the applicability of electronic friction approaches for light-driven hot electron chemistry.

Electronic friction based on time-dependent perturbation theory
In MDEF simulations, the nuclear dynamics are defined via a Langevin equation (LE) 21 : Herein, the three force components are the conservative force due to the potential energy V, the electronic friction force given as a product between the electronic friction tensor and the velocity of the atoms, and a temperature-and friction-dependent random force term that ensures detailed balance. The second and third terms describe the interaction of the adsorbate degrees of freedom with the bath degrees of freedom that represent the metal electrons. The two major assumptions to arrive at an LE are that the electron-nuclear coupling is weak and that it occurs instantaneously, i.e., coupling due to electron-hole pair (EHP) excitations has no memory of previous times and no dependence on the perturbing frequency (Markov approximation). The first assumption is problematic in cases where molecule-metal electron transfer significantly changes the underlying potential energy surfaces, such as in the case of nitrous oxide on Au(111). 37,38 The second approximation is equivalent to assuming that the electronnuclear coupling is constant regardless of nuclear velocities or vibrational frequencies.
The result is that, within 1 st order TDPT, we evaluate the electronic friction tensor using Fermi's Golden Rule 23,31 : Herein, the factor 2 accounts for spin multiplicity in the case of non-spin-polarised calculations and we have not included occupation factors that arise from finite temperature state populations. As our DFT calculations generate a finite number of states at discrete points k in momentum space, we need to interpolate to describe the friction tensor elements, which correspond to relaxation rates due to electron-nuclear coupling along the (mixed) Cartesian directions. This can be done in two ways: By choosing a broadening large enough to achieve convergence and small enough not to affect relaxation rates one approaches the zero-frequency limit, which is the formally correct limit within 1 st order TDPT. It was recently shown that this level of description struggles to capture vibrational energy loss for dense molecular overlayers 33 and it does not correspond to the correct limit for gas-surface scattering with incidence energies above several hundred meV. As detailed in our previous work 34 and by Spiering and Meyer, 35 we choose the second approach, namely to replace the delta function with a normalized Gaussian function of a finite width (0.6 eV). With such a broadening, we effectively model the inclusion of higher energy EHPs. It has recently been pointed out that such broadening introduces contributions that are not rigorously contained in a zero-frequency 1 st order approximation 32 and leads to relaxation rates that depend on the choice of broadening. This dependence remains a weak spot in this theory and has been discussed in detail in our previous work 31  Our computational settings regarding model set-up, basis set, and Brillouin zone sampling have been detailed in a previous publication. 34

Symmetry mapping scheme
To solve Eq (1) efficiently, one has to evaluate the potential energy and friction tensor at each trajectory step analytically. Neural network (NN) based representations of the PES 44 and the embedding electron density 45 (directly associated with the LDFA friction coefficient) which fully comply with the surface symmetry, have been described previously. We thus focus here on the more complex representation of the anisotropic friction tensor computed by TDPT. Because the TDPT-based friction tensor is directionally-dependent on molecular coordinates, the symmetry equivalent molecular configurations that share the same potential energy and LDFA friction coefficient, would have totally different but symmetrically-correlated friction tensor elements. This feature immediately renders the extension of any standard methodology for PES construction difficult for representing the friction tensor as a function of molecular geometry. Specifically, we discuss the case of a diatomic molecule (H2) on a rigid metal surface, in which six coordinates are involved. These six coordinates can either be the Cartesian coordinates of (x1,y1,z1,x2,y2,z2) or internal coordinates (X,Y,Z,r,,). In the latter, the first three specify the center of mass of the molecule while the rest denote the internuclear distance and the polar and azimuthal angles to the surface normal.

Spiering and Meyer have very recently reported an approximate NN-based
representation of the 6×6 friction tensor for H2+Cu(111) system. 35 They first described individual friction coefficients in terms of four symmetry independent internal coordinates, as illustrated in Fig. 1, namely X, Y, Z, and r with angular coordinates fixed (e.g., θ0=φ0=90°). Since the tensorial friction ( Λ ) is a symmetric matrix, only the upper triangular elements () ij ij  are to be determined in practice. This treatment avoids the complex symmetry involved in the full six-dimensional representation. Once (X, Y, Z, r) with θ0=φ0=90° is obtained, the friction tensor at an arbitrary configuration depending on all six internal coordinates can be approximately generated by the following transformation, 35 where T(θ, φ) is the transformation matrix that rotates configuration (θ, φ) to the corresponding reference geometry with θ0=φ0=90°. Although this decomposition of the full 6×6 friction tensor certainly depends on the choice of the reference orientation, Spiering and Meyer have shown that the dependence is relatively weak near the minimum energy path (MEP).
In a previous publication 36 and in this present work, we propose an alternative strategy. We choose to fit the elements of the full 6×6 matrix form of the electronic friction tensor in Cartesian coordinates inside an irreducible triangle in the (111) surface unit cell (marked in red in Fig. 1). In this triangle, each friction tensor element is symmetry unique with a one-to-one mapping to a molecular geometry, thus posing no problems for fitting. In practice, since we have collected many data scattered inside and out the unit cell, we first move those data distributed outside this triangle into it. This can be done by first translating a data point outside the unit cell into this region, which does not change the friction tensor, followed by a series of reflections when necessary. 46 For each reflection, the new Cartesian coordinate vector  q and friction tensor  Λ can be obtained by,

Fitting and evaluation
After the mapping procedure, a single NN is utilized to describe the one-to-one mapping between the molecular geometry and each individual element of the friction tensor in the triangle. Six Cartesian coordinates of the molecule are set as the input layer of the NN, which pass geometric information to two hidden layers with 20 and 40 neurons via hyperbolic tangent activation functions. A hybrid algorithm combining the extreme learning machine and Levenberg-Marquardt (ELM-LM) algorithm was used to train the NN. 50,51 The data set consists of ~1200 points sampled by previous AIMD trajectories 34 and ~2600 additional points used for constructing the PES 44 that cover the dynamically relevant region. Since the friction coefficients are very small far above the surface where the electron density vanishes, we exclude data points in that region.
Instead, a switching function was imposed to smoothly dampen the friction coefficients to zero when the molecule-surface interaction is negligible (Z∞, in Å), One should note that the NN fitting of each friction coefficient is essentially the same as that of the PES in previous works and very robust. [47][48][49] However, for our purpose, to evaluate a friction tensor for any point outside this symmetry-unique region, we need to first find the symmetry identical geometry of that point in the irreducible triangle by the transformation described above, then calculate the NN interpolated friction tensor there, followed by inversely transforming the obtained full 6×6 friction tensor back to that corresponding to the original geometry.

Molecular dynamics simulations with tensorial friction
The molecular dynamics (MD) and MDEF trajectories have been computed with various initial conditions sampled in a quasi-classical way as described in our previous publication. 36 The initial molecular center is chosen to be far above the Ag (111)  and k is the Boltzmann constant. However, we neglect surface temperature in most simulations because the contribution of the random force is typically small, except where we discuss temperature effects explicitly.

Validation of the NN interpolation of the friction tensor
Based on our previously published analysis, 54 we have further examined the quality of our NN representation of the friction tensor. Fig. 2a shows the root mean square errors (RMSEs) averaged over all data for the 21 independent friction coefficients. The overall RMSEs uniformly range from ~0.01 ps -1 to ~0.03 ps -1 without any outliers, which are generally comparable or smaller than those reported for the H2+Cu (111) system. 35 In Fig. 2b Clearly, there are strong couplings between vibration (r) and translation (Z) directions, and between parallel motion (X/Y) and rotation (θ). Although the coupling terms can be both negative and positive, their absolute values increase as the molecule approaches the transition state. It has been previously shown that our NN diagonal tensor elements agree well with TDPT ones, 36 and we also find here the excellent agreement between NN and TDPT off-diagonal friction coefficients. As a further validation, Fig. 2c compares two-dimensional electronic friction contour plots, as a function of Z and r along the MEP, generated from NN interpolation and direct DFT calculations, respectively. Our NN representation faithfully captures the multi-dimensional topography of the friction tensor element, even though these data points were actually not included in the training set. These results indicate that the NN representation of the friction tensor is sufficiently accurate for subsequent MDEF simulations.

The physical differences between LDFA and TDPT
In our previous work, 36  This level of difference is position and direction dependent and can be directly seen in Fig. 3b. In the following we will study how these different descriptions translate into measurable dynamic observables and how they depend on system parameters.

The sensitivity of the vibrational deexcitation probability on the magnitude of friction
In our previous results on the vibrational deexcitation probabilities from H2(v=2, j=0) to H2(v=1, j=0) (i.e. P2→1) during H2 scattering on Ag(111), we found that both TDPT and LDFA yield comparable deexcitation probabilities that are slightly elevated compared to the adiabatic MD case (see Fig. 4a). This is in contrast to findings for H2 This can be related to a balanced effect between the longer interaction time and smaller molecular velocity at this incidence energy.
Furthermore, it is interesting to analyse whether the nonadiabatic effects stem from the anisotropy along different directions or the mode-coupling between these directions in the friction tensor. In Fig. 4c, we compare the P2→1 with TDPT electronic friction where the off-diagonal coupling between the intra-atomic coordinate and the moleculeadsorbate coordinate (TDPT) is included and where it is artificially set to zero (TDPT-Λ =0 ). To make the comparison clearer, a factor of two is also multiplied in both cases.
Although the magnitude of rZ  is not small as displayed in Fig. 2b, interestingly, the results with zero rZ  are quite close to the original ones in the entire energy range, implying the negligible effect of the mode-coupling between the translational and vibrational DOFs in this case. In addition, to study the importance of directional anisotropy, we average the three diagonal components of the TDPT friction tensor (twice in magnitude) in the Cartesian space for each atom, followed by assigning the same friction coefficient for each atomic Cartesian coordinate in the following way, and zeroing all off-diagonal terms (TDPT-ISO). This converts the TDPT friction tensor to a form comparable with LDFA. It is found in Fig. 4c that the calculated P2→1 in this condition become drastically decreased, but very similar to the unscaled LDFA ones, indicating that the anisotropy of the friction tensor plays an important role in this process. This is consistent with what Meyer and Spiering have found for the H2+Cu(111) system 35 and hints to the fact that larger translational energy loss as predicted by LDFA really stems from its neglect of the directional anisotropy of the electronic friction tensor. It should be noted that while statically isotropic, the LDFA (much as TDPT) description does provide a certain level of dynamic directional dependence during molecular motion within MDEF due to the directional dependence of the velocity profile.

Dependence on initial conditions and final states
Our previous analyses suggested that the velocity profile and potential energy landscape also play a significant role in determining the nonadiabatic effects in the molecular scattering process. 36 One can of course change the velocity profile by choosing various initial conditions. We first study the influence of the angular velocity by performing MDEF simulations with initially different rotational and orientational states at the same (v=2) vibrational state. Fig. 5a compares the vibrational deexcitation probabilities for H2 (v=2, j=5, mj=0/5) to H2(v=1, j=5) with and without electronic friction. These two initial states correspond to the cartwheel (mj=0) and helicopter (mj=j) orientations, respectively, which would impact the surface with similar vibrational and translational, but distinctly different angular velocities for a given incidence energy.

Electronic temperature
Within the electronic friction picture, plasmon driven or laser-heated hot-electron effects are typically simulated via an increase of the electronic temperature that defines the heat of the electronic bath and the strength of the random force contribution. This can be done in form of a constant elevated temperature or by simulating a timedependent temperature profile that is coupled with a two-temperature model of the coupled lattice-electron dynamics of the metallic substrate. 29,30 Electronic temperatures that arise from laser-heating a surface can reach transient temperatures up to several thousands of Kelvin. It should be noted that this approach will assume that all hot electrons are thermalised and the weak coupling assumption will neglect any possible resonant electron-nuclear coupling effects due to high-lying excitations. Furthermore, the use of a temperature ramp is not strictly valid within an equilibrium Langevin theory and corresponding implications on detailed balance need to be carefully considered. In a previous work, we studied the effect of including electronic temperature in eq. (2) and found that over a range of up to 6000 K, electronic relaxation rates within our currently employed level of theory depend very weakly on temperature. 31 We therefore describe the electronic friction tensor as temperature independent. However, within the dynamic formalism of eq. (1) temperature still enters the Langevin expression and it is interesting to study how dynamic state-to-state scattering is affected by varying the electronic bath temperature. Whereas previous sections studied scattering at low temperature conditions and neglected the random force, Fig. 7 shows the dependence of vibrational deexcitation as a function of different electronic temperatures (0K, 1000K, and 6000K).
We find that, while overall deexcitation probabilities do not significantly change for both friction models, higher temperatures lead to higher integrated energy loss and increased deexcitation probabilities for very low translational energies. For the LDFA method, this effect becomes relevant at lower temperatures than for the TDPT method.
This can be related to the higher average friction magnitude of the LDFA method. The increase of deexcitation probability at low translational energies can be directly understood by the increased residence time of molecules close to the surface and the increased exposure to random force perturbations. From our results, it appears that state-to-state scattering is not strongly affected by plasmonic heating. This will, however, be different for rates of reactions such as recombinative desorption and dissociation reactions with high barriers and physisorbed precursor states where surface residence times are much larger.

Conclusions
We have presented extensive MDEF simulations within the Markov approximation for the state-to-state scattering of molecular hydrogen from a Ag(111) surface. These simulations are enabled by the construction of neural-network based potential energy surfaces and electronic friction tensor with data calculated from Density Functional Theory and 1 st order TDPT. We find that nonadiabatic dynamic state-to-state scattering is strongly dependent on the absolute magnitude of friction with respect to the relevant internal coordinate, whereas coupling between internal coordinates does not appear to play a significant role. Different vibrational and rotational states are affected very differently by electronic friction; however, the two different friction models did not yield drastically different results. The rotational state dependence hereby is mostly due to the changed molecular orientation and particular trajectory route with respect to the surface during impingement that boost vibrational or translational energy loss. We have furthermore shown that, the inclusion of elevated electronic temperatures as they appear in plasmon-driven hot-electron chemistry, dominantly affect reaction probabilities at low incidence energy, where surface residence times are larger than at high incidence energies. Especially the latter finding suggests that, plasmon-driven MDEF studies with electronic friction directly computed from DFT orbitals are possible and could provide a basis for further systematic improvement beyond the current state-of-the-art represented by LDFA. 29,30 As guided by recent work of Novko et al., 33 a word of warning on the Markov and zero-frequency approximation is in order. Our methodology is numerically robust but suffers from the choice of a broadening parameter that, despite a weak dependence in the chosen regime, does introduce some variation of frictional magnitude. Our results have shown that this dependence can strongly affect the dynamical scattering results on a level that would, at best, allow a qualitative prediction of reaction rate trends for hot electron chemistry. We conclude from this that, to investigate the ability of MDEF to act as a predictive tool for plasmonic science, the Markov approximation needs to be overcome. This endeavor is currently under way.

Conflicts of interest
There are no conflicts to declare.