Within a 2-D scattering model, we investigate the vibrational relaxation of an idealized molecule colliding with a metal surface. Two perturbative nonadiabatic dynamics schemes are compared: (i) electronic friction (EF) and (ii) classical master equations. In addition, we also study a third approach, (iii) a broadened classical master equation that interpolates between approaches (i) and (ii). Two conclusions emerge. First, even though we do not have exact data to compare against, we find there is strong evidence suggesting that EF results may be spurious for scattering problems. Second, we find that there is an optimal molecule-metal coupling that maximizes vibrational relaxation rates by inducing large nonadiabatic interactions.

## I. INTRODUCTION

Non-adiabatic dynamics are known to play an essential role in photochemistry, and excited state dynamics play an essential role in the gas phase and in solution. For such experiments,^{1–6} it is obvious that photo-excitation is followed by energetic relaxation as electrons relax and nuclei heat up: after all, the density of states for nuclear motion is much larger than the density of states for electronic motion and thus, after a long time, all electronic energy must be converted into heat (or nuclear motion). Thus, in solution or gas phase photochemistry, nonadiabatic dynamics are paramount.

Now, using the same logic, consider the role of nonadiabatic effects at a metal surface. On the one hand, a bulk metal carries an enormous density of phonon states; these states will usually be the final acceptors of any excess energy and thus act as a driving force for nonadiabatic transitions. On the other hand, a metal also carries a large density of electronic states, and thus electronic transitions are possible (even without nuclear motion). Thus, there is no guarantee that one will observe nonadiabatic dynamics near metal surfaces, i.e., a nontrivial coupling of nuclear motion with electronic transitions. In a series of recent papers, however, the Wodtke group has given very convincing evidence that nonadiabatic effects *are ubiquitous* when studying scattering processes for molecules off of metal surfaces. In the most famous experiments, the Wodtke group has scattered NO molecules across Au(111) surfaces^{7–13} and found clear evidence that vibrational relaxation is mediated by nonadiabatic processes (in this case, transient electron transfer). Thus, there is currently a great deal of interest in the physics and chemistry communities regarding how to model these difficult experiments.

Now, obviously, with a metal substrate, a fully quantum description of the nuclear and electronic degrees of freedom would be prohibitively expensive.^{14} For realistic, multi-atom simulations of the Wodtke experiments, semiclassical treatments are the only possible way forward. And, in this context, there are two well-known perturbative limits.^{15} (i) On the one hand, for weak nonadiabatic effects (i.e., strong metal-molecule couplings and weak electron-phonon couplings), the usual semiclassical framework is to assume that the molecular motion on the metal surface feels the so-called “electronic friction (EF)” from the bath of metallic electrons. This concept of electronic friction has been used many times in the past,^{16–21} most famously by Head-Gordon and Tully to study the relaxation of CO on a Cu substrate;^{16} recent work by Juaristi and Reuter *et al.* has also successfully investigated the relaxation behavior for CO on Cu,^{20} as well as N on Ag,^{21} using a local density approximation for the electronic friction tensor. (ii) On the other hand, for strong nonadiabatic effects (i.e., weak metal-molecule couplings and strong electron-phonon couplings), another approach is a classical master equation (CME), whereby molecules move as if they are charged or uncharged with stochastic hops between different charge states. This master equation approach describes the processes known as dynamics induced by multiple electronic transitions (DIMET) by the surface-science community.^{22,23} In the context of standard electron transfer theory, one imagines that the former approach should describe the inner-sphere heterogeneous electron transfer and the latter approach should describe the outer-sphere heterogeneous electron transfer.^{24}

Unfortunately, for the case of NO scattering off of gold, neither EF nor CME may be valid. First, the NO–Au interaction is not small at short distances, and so the CME approach is likely inapplicable. Second, far from the metal, the NO–Au interaction is weak, and the Wodtke group has given three pieces of evidence suggesting that electronic friction is inapplicable for scattering problems:^{7–9}

For NO incoming in a highly excited state (e.g.,

*n*_{vib}= 15), vibrational relaxation shows no barrier (as one would expect with adiabatic dynamics). Furthermore, the most probable exit channel has*n*_{vib}= 7. However, if the metal Au is replaced by an insulator, LiF, a barrier does appear, and the most probable exit channel is the original vibrational state (and the second most probable exit channel is the original vibrational state minus one). These data obviously suggest strong coupling between electronic transitions in the metal and vibrational transitions in the molecule, and the resulting dynamics demonstrate sudden changes as opposed to gradual changes of state (not as we would expect with electronic friction).For NO incoming in a highly excited state (e.g.,

*n*_{vib}= 15), one can observe hot electron emission from the metal at very large kinetic energies. Such emission precludes simple electronic friction descriptions based on fast electronic equilibration.For NO incoming in the ground state, Wodtke

*et al.*have observed that multiple quanta can be excited directly, which would not agree with a frictional description, i.e., a golden-rule picture of the dynamics assuming small electron-phonon couplings.

Thus, the Wodtke experiments present a clear challenge to theoretical chemistry and physics. Since the relevant dynamics have strong electron-phonon couplings, and because the metal-molecule couplings can be weak or strong (depending on the distance to the metal), and since transient electron transfer cannot be ignored, the Wodtke experiments simply do not sit in any simple perturbative regime.^{15} To model these difficult dynamics, a few years ago, Shenvi *et al.* suggested discretizing the metal continuum and they developed the so-called independent-electron surface hopping (IESH) approach.^{25,26} Quantum master equations have also been proposed for this purpose.^{27,28} More recently, by rederiving the origins of electronic friction, our research group showed how to extrapolate between the CME and EF regimes so that one could develop a universal, semiclassical nonadiabatic dynamics algorithm for strong or weak coupling near a metal surface. We labeled the resulting algorithm a broadened CME (BCME) approach.

With this background in mind, we have two goals for the present article. First, we would like to investigate the *consequences and signatures of nonadiabatic effects* for a diatomic molecule scattering off of a metal surface. Experimental signatures of nonadiabatic dynamics have been suggested by Wodtke *et al.*, and we would like to see how many of these signatures can be studied theoretically. To isolate these dynamical effects, we will work with a 2D model that will allow a thorough analysis. Second, to guide our understanding of the relevant process, we would also like *to compare and contrast three different nonadiabatic dynamics approaches:* (i) EF, (ii) CME, and (iii) BCME. Because we do not have an exact propagator, it is essential that we analyze multiple approaches. Naturally, since the EF and CME algorithms are based on perturbation theory, these algorithms must be accurate within their own, respective, parameter regimes. However, in a non-perturbative regime, it will be crucial to have different approaches so that we can make the best guess for the correct answer. In the course of our results, we will point out several surprising features that arise from these different methods. At present, our hypothesis is that, of the three methods above, BCME dynamics are the most reliable.

This paper is organized as follows: in Sec. II, we review all three dynamics schemes discussed above; in Sec. III, we define our 2D model Hamiltonian and provide details of the simulation; in Sec. IV, we show the simulation results for different sets of parameters; in Sec. V, we discuss the results and highlight why sometimes EF can yield vibrational relaxation rates that are too small, while at other times EF can yield vibrational relaxation rates that are too large; in Sec. VI, we conclude with a few suggestions for future work.

### A. Notation

The notation used in this paper is as follows: bold characters (e.g., ** r**) are vectors, bold characters with a left-right arrow (e.g., $\Lambda \u2194$) are tensors, plain characters (e.g., H) are either scalars or operators; for indices, we use Greek letters (

*α*,

*β*, etc.) for nuclei and Roman letters (

*i*,

*j*, etc) for electrons;

**and**

*r***always represent position and momentum vectors for molecules, respectively.**

*p*## II. THEORY

Henceforward, we will consider an idealized molecule (or impurity) on a metal surface using the Anderson-Holstein (AH) model in a nuclear *D*-dimensional space,

Here and below, the Fermi level of the metal is always chosen to be zero. *d* and *d*^{$\u2020$} are annihilation and creation operators for the impurity site, *c*_{k} and $ck\u2020$ are annihilation and creation operators for the *k*th orbital in the electronic bath. When $d\u2020d=1$, the impurity is occupied and we will speak of the molecule as being an anion; when $d\u2020d=0$, the impurity is unoccupied and we will speak of the molecule as being neutral. We define *U*_{i}(** r**) as the potential energy surface for the electronic state $i$ at position

**so that**

*r**U*

_{1}(

**) ≡**

*r**h*(

**) +**

*r**U*

_{0}(

**) is the potential energy surface for the anion. Vice versa,**

*r**h*(

**) ≡**

*r**U*

_{1}(

**) −**

*r**U*

_{0}(

**) is the energy gap between the anion state and the neutral state.**

*r**V*

_{k}(

**) denotes the coupling between the**

*r**k*th bath orbital and the impurity site.

*ϵ*

_{k}is the energy of the

*k*th bath orbital. Whenever possible, we apply the wide band approximation (WBA) and assume that the self-energy of the impurity has only an imaginary part which does not depend on

*k*or

*ϵ*,

where *ρ*(*ϵ*) is the density of states in the metallic bath.

### A. Electronic friction

In this paper, we will study the dynamics of the AH model using several approaches, especially electronic friction (EF).^{21,29–32} According to this model, all nonadiabatic effects are wrapped up into stochastic Langevin dynamics on a potential of mean force. Thus, the equations of motion are

where $M\u2194\u22121\u2261m1\u22121m2\u22121..$ is the inverse mass tensor, $\Lambda \u2194$ is the friction tensor, and *δ*** f** is the random force.

**is the mean force acting on the nuclear degrees of freedom,**

*F*and, for future reference, the potential of mean force is given by

where *ζ*(*r*_{0}) is some arbitrary reference potential.

In Eq. (4), the spectral function *A* and Fermi function *f* are

The reader may be well surprised that the bandwidth *W* appears in Eq. (4), given that we would like to take the wide-band limit. In fact, a finite *W* is required in this case to make sure that the integral in Eq. (4) does not diverge, given the form of ** K** in Eq. (7).

^{34}In practice, we choose

*W*$\u226b$ Γ or, to be specific,

*W*is at least 10 times larger than Γ

_{0}

^{50}[see Eq. (20)].

For a two-state model, with multiple nuclear degrees of freedom, the proper electronic friction tensor is^{33}

where ⊗ denotes an outer product. Equation (8) can be derived from many different approaches, including perturbation theory on top of a Meyer-Miller transformation,^{29} non-equilibrium Green’s functions,^{31} path integrals,^{35} or most generally, a projection approach applied to the quantum-classical Liouville equation.^{36} Although the original Head-Gordon Tully (HGT) friction model applies only at zero temperature,^{29,33} Eq. (8) is equivalent to Tully’s recent extrapolation for friction at finite temperature.^{37} Finally, *δ*** f**(

**,**

*r**t*) is the Markovian random force,

From the expressions above, one immediately finds a troubling attribute of electronic friction tensors. For some Hamiltonians, it is possible to encounter geometries where Γ(** r**) → 0 but $\u25bf$

*h*(

**) $\u2260$ 0. In such a case, the corresponding matrix elements in $\u2194\Lambda (r)$ [Eq. (8)] will diverge to infinity because $\u222bA2(\u03f5,r)d\u03f5\u221d1\Gamma $ as Γ → 0 [see Eq. (6)]. To avoid such a numerical instability, below we will choose a small artificial parameter Γ**

*r*_{cutoff}for our simulations such that for Γ(

**) < Γ**

*r*_{cutoff}(corresponding to very small molecule-metal coupling), we will ignore any effect from the electronic bath and set the friction and random force to 0,

We must always check whether or not our final results depend on Γ_{cutoff}. Unless stated otherwise, all data presented below are independent of Γ_{cutoff}.

### B. Classical master equations (CMEs)

Apart from electronic friction, classical master equations (CMEs) represent an entirely different approach for modeling nonadiabatic dynamics at metal surfaces. The CME approach^{38,39} treats electronic states explicitly and proposes stochastic trajectories. More specifically, nuclear trajectories are propagated either along *U*_{0} or *U*_{1} and, for each time step, the particle may hop from one surface to the other. The probability to hop is decided by the hybridization function Γ. This scheme is summed up by the following equations of motion for the probability densities:

Here *P*_{i}(** r**,

**,**

*p**t*) denotes the probability density to find a particle at phase point (

**,**

*r***) in the electronic state $i$ at time**

*p**t*. The CME in Eq. (11) can be derived by assuming (i) a high temperature such that classical nuclear motion suffices and (ii) a small hybridization function Γ <

*kT*.

#### 1. Relating classical master equations to electronic friction

For large enough Γ, with many hops back and forth between surfaces, Ref. 39 demonstrates that CME dynamics become equivalent to EF dynamics. Because this equivalence is key for deriving BCME dynamics (as discussed in Sec. II C), we will briefly review how such an equivalence can be demonstrated. First, we change variables to two new densities *P*_{A}(** r**,

**,**

*p**t*) and

*P*

_{B}(

**,**

*r***,**

*p**t*) as follows:

or, vice versa

Next, if we consider the equation of motion for the total density *P*_{A}(** r**,

**,**

*p**t*) and we assume that

*P*

_{B}(

**,**

*r***,**

*p**t*) is changing slowly, we find

### C. Broadened classical master equations

Finally, the last dynamics approach studied here will be an extrapolation of EF and CME dynamics, denoted a broadened CME (BCME)^{33,40} approach. The BCME approach is a natural extension of the CME if one wants to study large Γ. To include the broadening effects, one merely modifies the potential energy surfaces (PESs) *U*_{0} and *U*_{1} in Eq. (11) so that, in Eq. (14), the Fermi population *f*(*h*(** r**)) is replaced with the correctly broadened population

*n*(

*h*(

**)) [for the definition of**

*r**n*, see Eq. (17)]. Thus, the diabatic surfaces in Eq. (11) will now depend on both Γ(

**) and temperature. Although there is no single, unique means to modify Eq. (11)—because we specify only how the equation of motion for**

*r**P*

_{A}(

**,**

*r***,**

*p**t*) should change and not for

*P*

_{B}(

**,**

*r***,**

*p**t*)—the simplest BCMEs of motion are as follows:

^{33,40}

Here, the diabatic forces have been modified by the following correction:

The broadened population *n*(** r**) is defined as

For future reference, the broadened diabatic potentials of mean force are

In Eq. (16), *f* is the Fermi function and *f*(*h*(** r**)) represents the unbroadened, equilibrium population for the impurity site at position

**. By contrast,**

*r**n*represents the correctly broadened equilibrium population of the impurity site at position

**. Thus,**

*r**n*−

*f*indeed represents a broadening correction. We note that, for large enough Γ, the total probability density for BCME dynamics evolves on the same potential of mean force as EF dynamics [in Eq. (5)].

^{40}In Sec. IV A, we will plot and compare the unbroadened (

*U*

_{i}) and broadened ($Uib$) diabats.

## III. SIMULATION DETAILS

To study the methods above, we will simulate vibrational relaxation for a model two-dimensional (2D) system. Our 2D system has been roughly designed to mimic a scattering event, whereby a diatomic molecule impinges on a metal surface. The first dimension *x* corresponds to the vibrational degree of freedom (DoF) of the molecule, and the second dimension *z* is the molecular center-of-mass position. We assume that the metal surface is located at *z* = 0 and that an NO molecule approaches the surface from *z* = −∞. The energy surfaces we use are

Here *m*_{x} is the reduced mass for vibrational motion. The energy surfaces along the *x* direction are harmonic wells, where the eigenfrequency *ω* is chosen to be the same for $0$ and $1$. *A*_{0} represents well the depth of the Morse potential *S*_{0}, and following Ref. 41, we fix *A*_{0} as 300 meV. *B*_{1} is the difference between the surface work function and NO electron binding strength, and thus we fix *B*_{1} to be 5.55 eV (the work function of gold is 5.1 eV and the electronic affinity of NO is 0.45 eV^{41}). *x*_{0} and *x*_{1} are chosen to mimic the equilibrium bond length differences between NO (1.15 Å) and NO^{−} (1.25 Å) so that *x*_{1} − *x*_{0} = 0.1 Å. Other parameters such as *C*_{0}, *z*_{0}, and *C*_{1} are chosen such that the potential surfaces look reasonable, given that the energy surfaces in the *z* direction (*S*_{0}(*z*), *S*_{1}(*z*)) resemble the electron mediated model for NO proposed by Newns.^{41} The second term in the expression for *S*_{1}(*z*) does not appear in the Newns model but has been added to ensure that the impinging NO particles scatter back (rather than penetrate the metal). The metal surface is effectively located around *z* = 0.

For the hybridization function Γ(** r**), we choose

In the *x* direction, Γ(** r**) has a maximum near the equilibrium position of the $0$ state (i.e.,

*x*

_{0}); in the

*z*direction, the coupling Γ(

**) decreases exponentially as the distance between the particle and surface increases, and Γ(**

*r***) goes to 0 as**

*r**z*→ −∞.

Almost all of the parameters listed above are defined in Table I (except for the hybridization function Γ_{0} and the displacement *x*_{1}). Note that the temperature here is relatively low and should not satisfy the “high temperature” prerequisite for the CME dynamics. That been said, the experiments start in a hot vibrational state *n*_{vib} = 15 (which makes the classical vibrational energy *E*_{vib} $\u226b$ $\u210f$*ω*) such that classical dynamics may well still be valid. Furthermore, in this paper, we will also study the dynamics with the BCME to include broadening. In a future publication, we will consider these dynamics with a broadened version of a quantum master equation (QME) to include broadening plus nuclear quantum effects. For now, our major concern is how will the dynamics depend on different values of *x*_{1} and Γ_{0} (as well as in the incoming momentum in the *z* direction, $p0$).

Parameter . | Value (a.u.) . | Comment . |
---|---|---|

m | 55 000 | Mass of particle |

m_{x} | 14 000 | Reduced mass |

kT | 0.001 | Temperature |

ω | 0.008 | Harmonic frequency |

x_{0} | 0 | Parameter in Eq. (19) |

A_{0} | 0.011 | Parameter in Eq. (19) |

C_{0} | 0.64 | Parameter in Eq. (19) |

z_{0} | −3.5 | Parameter in Eq. (19) |

B_{1} | 0.2 | Parameter in Eq. (19) |

C_{1} | 0.67 | Parameter in Eq. (19) |

K_{g} | 4 | Parameter in Eq. (20) |

c_{g} | 0.64 | Parameter in Eq. (20) |

z_{g} | 0 | Parameter in Eq. (20) |

W | 1.5 | Bandwidth in Eq. (4) |

Parameter . | Value (a.u.) . | Comment . |
---|---|---|

m | 55 000 | Mass of particle |

m_{x} | 14 000 | Reduced mass |

kT | 0.001 | Temperature |

ω | 0.008 | Harmonic frequency |

x_{0} | 0 | Parameter in Eq. (19) |

A_{0} | 0.011 | Parameter in Eq. (19) |

C_{0} | 0.64 | Parameter in Eq. (19) |

z_{0} | −3.5 | Parameter in Eq. (19) |

B_{1} | 0.2 | Parameter in Eq. (19) |

C_{1} | 0.67 | Parameter in Eq. (19) |

K_{g} | 4 | Parameter in Eq. (20) |

c_{g} | 0.64 | Parameter in Eq. (20) |

z_{g} | 0 | Parameter in Eq. (20) |

W | 1.5 | Bandwidth in Eq. (4) |

In Fig. 1, we plot the individual components making up the diabatic potential energy surfaces in Eq. (19). In Fig. 2, we plot the total potential surfaces in the *z* direction (for one fixed *x*), and we show the effects of broadening. We also plot the hybridization function Γ(** r**).

For each calculation reported below, we have run 5000 trajectories. To roughly simulate the Wodtke experiments,^{8} each trajectory was initialized to the 15th vibrational state. In the *x* direction, trajectories were initialized with a microcanonical ensemble: we weighted all (*x*, *p*_{x}) satisfying $Ex\u2009=\u2009Evib0\u2009=\u200915.5\u210f\omega $ equally. In the *z* direction, trajectories were initialized at position *z* = −15 and the momentum *p*_{z} was chosen from a Gaussian distribution with $p0$ and $\sigma =mkT2$. We used a time step *dt* = 0.25 a.u. and propagated trajectories for 2 × 10^{6} steps. Dynamics were carried out with the velocity-Verlet propagator. Unless stated otherwise, the trapped particles were not considered, and we analyzed exclusively reflected particles (which were collected at *z* = −20).

## IV. RESULTS

We will now report our results, focusing mostly on the overall amount of predicted vibrational relaxation.

### A. Dynamics

In Fig. 3, we plot the number of particles collected at *z* = −20 as a function of time. In this case, for the CME dynamics, we find very few particles trapped near the surface. For BCME dynamics, ∼10% of the particles are trapped, and for EF dynamics, more than 20% of the particles trap in the well near *z* = −0.9. In general, we find that this trend holds for most calculations below with different parameters: the EF results usually result in far more trapping that CME or BCME dynamics (see Table II).

Parameters . | Percentage trapped . | Percentage reflected on state $1$ . | |||||
---|---|---|---|---|---|---|---|

x_{1}
. | Γ_{0}
. | $p0$ . | CME (%) . | BCME (%) . | EF (%) . | CME (%) . | BCME (%) . |

0.2 | 0.01 | 20 | 0.02 | 0.04 | 0.02 | 0.04 | 0.08 |

0.2 | 0.03 | 20 | 0.02 | 1.36 | 0.72 | 0.14 | 0.08 |

0.2 | 0.05 | 20 | 0.02 | 10.24 | 17.44 | 0.06 | 0.11 |

0.2 | 0.08 | 20 | 0.04 | 37.06 | 55.34 | 0.08 | 0.03 |

0.2 | 0.03 | 40 | 0.08 | 7.44 | 20.98 | 0.06 | 0.09 |

0.2 | 0.03 | 60 | 3.60 | 6.26 | 18.50 | 0.10 | 0.09 |

0.2 | 0.03 | 80 | 1.22 | 0.78 | 0.32 | 0.06 | 0.10 |

0.4 | 0.01 | 20 | 7.98 | 8.52 | 2.06 | 0.04 | 0.07 |

0.4 | 0.03 | 20 | 3.26 | 4.42 | 2.24 | 0.10 | 0.08 |

0.4 | 0.05 | 20 | 1.34 | 3.24 | 31.34 | 0.14 | 0.12 |

0.4 | 0.08 | 20 | 0.34 | 1.98 | 7.60 | 0.08 | 0.10 |

0.4 | 0.03 | 40 | 9.38 | 7.48 | 23.22 | 0.09 | 0.06 |

0.4 | 0.03 | 60 | 9.68 | 4.90 | 8.24 | 0.04 | 0.06 |

0.4 | 0.03 | 80 | 3.74 | 1.88 | 5.48 | 0.04 | 0.04 |

Parameters . | Percentage trapped . | Percentage reflected on state $1$ . | |||||
---|---|---|---|---|---|---|---|

x_{1}
. | Γ_{0}
. | $p0$ . | CME (%) . | BCME (%) . | EF (%) . | CME (%) . | BCME (%) . |

0.2 | 0.01 | 20 | 0.02 | 0.04 | 0.02 | 0.04 | 0.08 |

0.2 | 0.03 | 20 | 0.02 | 1.36 | 0.72 | 0.14 | 0.08 |

0.2 | 0.05 | 20 | 0.02 | 10.24 | 17.44 | 0.06 | 0.11 |

0.2 | 0.08 | 20 | 0.04 | 37.06 | 55.34 | 0.08 | 0.03 |

0.2 | 0.03 | 40 | 0.08 | 7.44 | 20.98 | 0.06 | 0.09 |

0.2 | 0.03 | 60 | 3.60 | 6.26 | 18.50 | 0.10 | 0.09 |

0.2 | 0.03 | 80 | 1.22 | 0.78 | 0.32 | 0.06 | 0.10 |

0.4 | 0.01 | 20 | 7.98 | 8.52 | 2.06 | 0.04 | 0.07 |

0.4 | 0.03 | 20 | 3.26 | 4.42 | 2.24 | 0.10 | 0.08 |

0.4 | 0.05 | 20 | 1.34 | 3.24 | 31.34 | 0.14 | 0.12 |

0.4 | 0.08 | 20 | 0.34 | 1.98 | 7.60 | 0.08 | 0.10 |

0.4 | 0.03 | 40 | 9.38 | 7.48 | 23.22 | 0.09 | 0.06 |

0.4 | 0.03 | 60 | 9.68 | 4.90 | 8.24 | 0.04 | 0.06 |

0.4 | 0.03 | 80 | 3.74 | 1.88 | 5.48 | 0.04 | 0.04 |

### B. Vibrational distribution

Let us now discuss the vibrational relaxation of the outgoing particles that are scattered backwards (and ignore all the trapped particles). With CME or BCME dynamics, because of the large energy penalty to emerge as an anion asymptotically, almost all (>99%) the reflected particles are found to lie on the neutral state $0$. For this reason, the vibrational state of each particle can be calculated as follows: (a) we compute the kinetic energy in the *x* direction, $Ekx=px2/2mx$, (b) we compute the potential energy in the *x* direction, *E*_{px} = *U*_{0}(*x*), and (c) we compute *n*_{vib} = (*E*_{kx} + *E*_{px})/$\u210f$*ω* − 0.5 and round it to an integer. This procedure can be applied for all the methods above (CME, BCME, and EF). Note that, for BCME dynamics, we may safely use *U*_{0}(*x*) [rather than $U0b(x)$] because at *z* = −20, Γ → 0, and the *U*_{0}(*x*) and the $U0b(x)$ surfaces have negligible differences.

#### 1. Relaxation dependence on $\Gamma 0$

In Fig. 4(a), we plot the vibrational distributions as a function of different Γ_{0} for fixed incoming z momentum $p0=20$. We observe the vibrational relaxation for both EF and BCME dynamics, while CME dynamics do not yield any relaxation. From these plots, it is straightforward to see that the CME dynamics fail for an obvious reason: nearly all particles are blocked by the energy barrier in the *z* direction, and they do not have enough energy to reach the surface crossing seam, see Fig. 1.

Focusing now on the EF and BCME dynamics, as one would expect, we find that relaxation becomes stronger as the molecule-bath coupling increases from zero, reaching a maximum ($nvib=7\u22128$) at Γ_{0} = 0.05. Note, however, that there is a turnover. As Γ_{0} increases even more, corresponding to the extreme adiabatic limit, vibrational relaxation slows down. Such a surprising turnover feature is not found in condensed phase electron transfer dynamics where the rate of electron transfer strictly increases as the molecule-metal coupling parameter Γ grows. See, e.g., Fig. 2 in Ref. 42 and Fig. 10.8 in Ref. 43. And yet this turnover is clearly analogous to Kramers’ theory,^{44} whereby the unimolecular escape rate from a well is maximized for a friction that is not too large or too small. Even though we are modeling transient vibrational relaxation (rather than activated nuclear barrier crossings), the same physics applies.

Regarding reliability, we note that the EF and BCME relaxation rates agree, especially in the adiabatic limit as Γ_{0} increases. Thus, even though we cannot propagate exact dynamics, we do calculate similar observables with two different and orthogonal methods. Furthermore, recent benchmarking of the BCME algorithm has suggested that BCME dynamics should be quite accurate with only two electronic states.^{51} For both of these reasons, we have a great deal of confidence in the data from Fig. 4(a), at least qualitatively. In Fig. 4(b), however, we plot the same result for the displacement *x*_{1} = 0.4, and we show that the agreement between EF and BCME does not always hold. While both methods predict more relaxation than the case of *x*_{1} = 0.2, the BCME approach predicts far more relaxation for small Γ_{0} than does EF. In this case, because CME dynamics can be derived with perturbation theory assuming small Γ, it is easy to argue that CME dynamics (and not EF) must be accurate here for Γ_{0} = 0.01. Furthermore, from the fact that BCME dynamics exactly agree with CME dynamics for small Γ_{0} and qualitatively agree with EF dynamics for large Γ_{0}, we hypothesize that the BCME dynamics should be meaningful over a wide range of parameter space. Our intuition is that EF dynamics will fail for large displacements (*x*_{1} − *x*_{0}) and small or moderately sized hybridization functions Γ_{0}.

#### 2. Relaxation dependence on incoming momentum

We now study how the incoming momentum affects relaxation. In Fig. 5(a), we plot the vibrational distributions for different $p0$ with a fixed value of the hybridization (Γ_{0} = 0.03) and *x*_{1} = 0.2. Here, for $p0\u226540$, the CME approach finally gives relaxation [compared against Fig. 4(a)]: there is enough energy to reach the diabatic crossing point. However, the CME does not agree with BCME or EF dynamics for small momenta. Regarding EF and BCME dynamics, we find that the relaxation rates are also in disagreement (though not completely different) for small incoming momenta.

For large momenta, however, we note that all the dynamic protocols (EF, CME, and BCME) roughly agree: apparently, because of the large incoming momenta, there are enough classical crossings such that friction results become meaningful, but this kinetic energy is also large enough such that broadening effects on the surface are unimportant. This agreement between the CME and EF dynamics has been seen before in 1D problems.^{42}

Finally, we consider the same dynamics now for the case of a larger displacement, *x*_{1} = 0.4. Here, we find again that there is no agreement between any of the methods for small incoming momentum. Because of the ability to interpolate, however, we hypothesize that the BCME dynamics are the most accurate. That being said, at larger incoming velocities, the methods do become more similar. Interestingly, though, at very large incoming velocities, all methods become very different again. These features yet cannot be easily explained. In general, we find that the BCME dynamics consistently predict more relaxation than the electronic friction as well as slightly wider vibrational distributions.

### C. Electronic energy released and hopping energy histograms

Recent experiments have measured the distribution of electronic kinetic energies excited in a metal surface as the result of molecular scattering.^{45,46} With this experimental fact in mind, we plot the energy distribution for hopping according to the CME/BCME dynamics for the four different simulations in Fig. 6. Note that such energy distributions cannot easily be extracted from EF calculations.

The results in Fig. 6 show that, for most hops, energy transfers from the incoming particle to the metal (i.e., the particle loses energy). Most hops occur near the surface crossing region with small energy gaps (|Δ*E*| < 0.02). Even so, large energy transfer events are possible within a single hopping event, which does explain the “multi-quanta relaxation” observed in Wodtke’s experiments.^{9,12} These results are discussed in Sec. V.

Finally, let us return to the turnover in vibrational rates as predicted by both EF and BCME dynamics in Fig. 4. To understand this turnover, Fig. 6 is quite useful. According to BCME dynamics, energy loss occurs only through hops. Thus, on the one hand, if Γ_{0} is small, the absolute number of hops is small and vibrational relaxation must be slow. On the other hand, if Γ_{0} is too large, most hops occur just after particles pass the energy crossing seam, i.e., when downward hops become energetically preferable. These hops yield very small energy loss and again vibrational relaxation must be slow. Thus Fig. 6 offers us some intuition for understanding the predicted turnover phenomenon.

## V. DISCUSSION

Thus far, we have found that both electronic friction and broadened classical master equations are able to capture many features of vibrational relaxation, and sometimes these two methods even agree. At this point, however, there are two key features which must be discussed in more detail.

### A. The BCME is more sensitive to displacement than EF

From Figs. 4 and 5, we observed that, although EF and BCME both yield larger relaxation rates when the displacement *x*_{1} is increased from 0.2 to 0.4, the BCME approach is obviously more sensitive to this change in parameter—especially for small Γ_{0} and small $p0$ cases. This sensitivity is obviously important because, for many molecules, the anionic and neutral potential energy surfaces can be very different. Furthermore, EF should be reliable only when these differences (i.e., electron-phonon couplings) are relatively small.

To explain the sensitivity of BCME dynamics, a figure will be very useful (Fig. 7). Here, we observe that, as the displacement *x*_{1} becomes larger, the surface crossing point as a function of the *x* coordinate drops in energy. As a result, if trajectories move along diabats, trajectories with a given *n*_{vib} will spend more time in the regions of large hopping probability. Furthermore, in these very regions, there is the chance to lose a larger amount of energy in one hop (see also Fig. 6).

Now, EF dynamics also predict stronger relaxation for large displacements—after all, the EF friction tensor is proportional to $\u25bf$*h* [see Eqs. (7) and (8)]. However, EF dynamics are not as sensitive to the displacement as are BCME dynamics because EF dynamics move along the adiabatic surface (rather than the diabatic surface) and a dramatic, sudden energy loss is impossible. Indeed, Wodtke and Tully *et al.* have argued that EF dynamics cannot produce multi-quanta relaxation because, by damping the nuclear motion, nuclear velocities change continuously in time, and thus any quantum mechanical extension of electronic friction must predict step-by-step dissipation of vibrational quanta.^{12,47}

### B. EF can be sensitive to Γ_{cutoff}

The very last feature that must be discussed is the artificial parameter, Γ_{cutoff}, which we have included earlier [in Eq. (10)], so as to determine when to apply or not apply the frictional damping and random force. The parameter Γ_{cutoff} can sometimes be crucial because, as explained in Sec. II A, in certain cases, one can find infinite friction for extremely small Γ(** r**). We must emphasize, however, that

**this divergence of friction is not an artifact.**In fact, this divergence in friction actually forces EF dynamics to recover Marcus’s theory of electrochemical charge transfer in the

*nonadiabatic regime*for a one-dimensional quantum Brownian oscillator model.

^{42}Thus, the friction tensor cannot be improved simply by smoothing away the divergence.

^{48}And yet, at the same time, the existence of an infinite frictional tensor must give one doubt about the overall applicability of EF dynamics. To investigate the practical consequences of this divergence, we will now modify the original

*T*(

*z*) model with the new parameters in Table III; this substitution forces

*T*(

*z*) to be much sharper than before.

*x*

_{1}is kept at 0.2. The modified parameters are plotted in Fig. 8. With these new parameters and surfaces, we report the relaxation rates from scattering simulations as a function of Γ

_{cutoff}.

From the data in Fig. 9(a), we find that the EF results are not equivalent for different values of Γ_{cutoff}. Indeed, for these parameters and such a small value of Γ_{0}, we find that the friction tensor is extremely large (nearly divergent). Thus, if Γ_{cutoff} is very small, we find that the EF dynamics can actually (and spuriously) predict more vibrational relaxation than BCME or CME dynamics. Luckily, this issue should not be important when Γ(** r**) is not infinitesimal near a surface crossing region, as shown in Fig. 9(b).

Note that, except for Fig. 9, all results reported in this paper using electronic friction can be considered reliable and converged with respect to Γ_{cutoff} (see Figs. 4 and 5). Obviously, looking forward, the fact that the BCME requires no such artificial parameter is a huge relative advantage. Unlike the case of EF, both the BCME and CME dynamics propose simple smooth dynamics along diabats in regions with Γ(** r**) → 0.

## VI. CONCLUSIONS AND FUTURE DIRECTIONS

In summary, we have investigated vibrational relaxation within a 2D scattering simulation where we expect the transient electron transfer for a variety of different approaches. Our conclusions are as follows:

We find that the CME approach is unable to predict accurate vibrational relaxation probabilities. Whenever the metal-molecule coupling Γ is large, CME dynamics along the simple diabatic curves are usually not accurate. In particular, the trajectory often misses the crossing region entirely. These dynamics usually disagree with BCME and EF dynamics (even for large Γ) and are a strong reminder that the propagating dynamics on raw diabatic surfaces is dangerous when Γ is really large.

We find that EF dynamics give reasonable probabilities of vibrational relaxation (and thus agree with BCME dynamics) when Γ(

) is reasonably large in the surface crossing region. However, there are clearly spurious effects when Γ becomes too small. We should emphasize that these effects likely represent the most severe failures possible for the EF approach. After all, the EF approach was originally designed to treat small electron-phonon couplings, whereas here the electron-phonon couplings can be large. Furthermore, Γ vanishes far away from the surface so that the assumption of large molecule-metal couplings is obviously violated. To address the shortcoming of the EF approach, the only future path forward would be to include non-Markovian effects.*r*^{49}Overall, the BCME approach appears to give the most sensible data. By construction, this algorithm mostly agrees with the CME algorithm (in the limit of small Γ) and with the EF algorithm (in the limit of large Γ). The BCME approach tends to be more sensitive to electron-phonon couplings, and the BCME approach usually results in more relaxation and a slightly wider distribution of vibrational quanta than do EF dynamics.

Finally, perhaps the most surprising conclusion of this work is the prediction that there is a turnover in the rate of vibrational relaxation for scattering experiments as a function of Γ. According to Fig. 4, we predict that the probability for vibrational relaxation peaks when Γ is neither too small nor too large. This turnover feature is not found in the condensed phase dynamics, where the rate of molecule-metal electron transfer is strictly increasing with the coupling parameter Γ: see Fig. 2 in Ref. 42. Instead it would appear that we have uncovered a situation where the Kramers’ turnover phenomenon reappears. As a practical matter, it would be very interesting to identify a series of different metal substrates with varying degrees of metal-molecule coupling (Γ) from which this trend could be confirmed experimentally.

Finally, looking forward, we have two clear next steps. First, given the simplicity of the BCME approach (which ignores electronic coherences for a two-state problem), it will be very interesting to compare the BCME algorithm with IESH^{25} (which includes coherences within the framework of a discretized metal). Such a comparison will tell us a great deal about when and why the BCME works/fails. Second, in order to apply the present dynamics to a real (and not model) system, it will be essential to extract (rather than conjecture) the relevant parameters from *ab initio* electronic structure calculations. This work is ongoing.

## ACKNOWLEDGMENTS

This material is based upon the work supported by the (U.S.) Air Force Office of Scientific Research (USAFOSR) PECASE award under AFOSR Grant No. FA9950-13-1-0157. J.E.S. acknowledges a Camille and Henry Dreyfus Teacher Scholar Award and a David and Lucille Packard Fellowship. We also acknowledge helpful suggestions from our JCP referees for interpreting the predicted turnover phenomenon in Fig. 4(a).