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) surfaces7–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., nvib = 15), vibrational relaxation shows no barrier (as one would expect with adiabatic dynamics). Furthermore, the most probable exit channel has nvib = 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., nvib = 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., ) 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; r and p always represent position and momentum vectors for molecules, respectively.
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 are annihilation and creation operators for the impurity site, ck and are annihilation and creation operators for the kth orbital in the electronic bath. When , the impurity is occupied and we will speak of the molecule as being an anion; when , the impurity is unoccupied and we will speak of the molecule as being neutral. We define Ui(r) as the potential energy surface for the electronic state at position r so that U1(r) ≡ h(r) + U0(r) is the potential energy surface for the anion. Vice versa, h(r) ≡ U1(r) − U0(r) is the energy gap between the anion state and the neutral state. Vk(r) denotes the coupling between the kth bath orbital and the impurity site. ϵk is the energy of the kth 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 is the inverse mass tensor, is the friction tensor, and δf is the random force. F is the mean force acting on the nuclear degrees of freedom,
and, for future reference, the potential of mean force is given by
where ζ(r0) 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 Γ or, to be specific, W is at least 10 times larger than Γ050 [see Eq. (20)].
For a two-state model, with multiple nuclear degrees of freedom, the proper electronic friction tensor is33
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 h(r) 0. In such a case, the corresponding matrix elements in [Eq. (8)] will diverge to infinity because as Γ → 0 [see Eq. (6)]. To avoid such a numerical instability, below we will choose a small artificial parameter Γ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 approach38,39 treats electronic states explicitly and proposes stochastic trajectories. More specifically, nuclear trajectories are propagated either along U0 or U1 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 Pi(r, p, t) denotes the probability density to find a particle at phase point (r, p) in the electronic state at time 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 PA(r, p, t) and PB(r, p, t) as follows:
or, vice versa
Next, if we consider the equation of motion for the total density PA(r, p, t) and we assume that PB(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) U0 and U1 in Eq. (11) so that, in Eq. (14), the Fermi population f(h(r)) is replaced with the correctly broadened population n(h(r)) [for the definition of n, see Eq. (17)]. Thus, the diabatic surfaces in Eq. (11) will now depend on both Γ(r) and temperature. Although there is no single, unique means to modify Eq. (11)—because we specify only how the equation of motion for PA(r, p, t) should change and not for PB(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 r. By contrast, n represents the correctly broadened equilibrium population of the impurity site at position r. Thus, 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 (Ui) and broadened () 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 mx 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 and . A0 represents well the depth of the Morse potential S0, and following Ref. 41, we fix A0 as 300 meV. B1 is the difference between the surface work function and NO electron binding strength, and thus we fix B1 to be 5.55 eV (the work function of gold is 5.1 eV and the electronic affinity of NO is 0.45 eV41). x0 and x1 are chosen to mimic the equilibrium bond length differences between NO (1.15 Å) and NO− (1.25 Å) so that x1 − x0 = 0.1 Å. Other parameters such as C0, z0, and C1 are chosen such that the potential surfaces look reasonable, given that the energy surfaces in the z direction (S0(z), S1(z)) resemble the electron mediated model for NO proposed by Newns.41 The second term in the expression for S1(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 state (i.e., x0); in the z direction, the coupling Γ(r) decreases exponentially as the distance between the particle and surface increases, and Γ(r) goes to 0 as z → −∞.
Almost all of the parameters listed above are defined in Table I (except for the hybridization function Γ0 and the displacement x1). 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 nvib = 15 (which makes the classical vibrational energy Evib ω) 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 x1 and Γ0 (as well as in the incoming momentum in the z direction, ).
Parameters used in the simulation.
Parameter . | Value (a.u.) . | Comment . |
---|---|---|
m | 55 000 | Mass of particle |
mx | 14 000 | Reduced mass |
kT | 0.001 | Temperature |
ω | 0.008 | Harmonic frequency |
x0 | 0 | Parameter in Eq. (19) |
A0 | 0.011 | Parameter in Eq. (19) |
C0 | 0.64 | Parameter in Eq. (19) |
z0 | −3.5 | Parameter in Eq. (19) |
B1 | 0.2 | Parameter in Eq. (19) |
C1 | 0.67 | Parameter in Eq. (19) |
Kg | 4 | Parameter in Eq. (20) |
cg | 0.64 | Parameter in Eq. (20) |
zg | 0 | Parameter in Eq. (20) |
W | 1.5 | Bandwidth in Eq. (4) |
Parameter . | Value (a.u.) . | Comment . |
---|---|---|
m | 55 000 | Mass of particle |
mx | 14 000 | Reduced mass |
kT | 0.001 | Temperature |
ω | 0.008 | Harmonic frequency |
x0 | 0 | Parameter in Eq. (19) |
A0 | 0.011 | Parameter in Eq. (19) |
C0 | 0.64 | Parameter in Eq. (19) |
z0 | −3.5 | Parameter in Eq. (19) |
B1 | 0.2 | Parameter in Eq. (19) |
C1 | 0.67 | Parameter in Eq. (19) |
Kg | 4 | Parameter in Eq. (20) |
cg | 0.64 | Parameter in Eq. (20) |
zg | 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).
A plot of the 2D model surfaces used in our simulation. The energy surfaces in the x direction (left) are both harmonic wells. The equilibrium positions for the neutral state and the anion state are x0 = 0 and x1 = 0.2 [see Eq. (19)], respectively. In the z direction (right), there is an energy barrier at z = −1.8 and an energy well at z = −0.9 (see arrows in the plot), which can potentially trap incoming particles.
A plot of the 2D model surfaces used in our simulation. The energy surfaces in the x direction (left) are both harmonic wells. The equilibrium positions for the neutral state and the anion state are x0 = 0 and x1 = 0.2 [see Eq. (19)], respectively. In the z direction (right), there is an energy barrier at z = −1.8 and an energy well at z = −0.9 (see arrows in the plot), which can potentially trap incoming particles.
The CME and BCME surfaces in the z direction. Here x = 0.0 and Γ0 = 0.03. Ui and are the CME and BCME surfaces for , respectively. UPMF is calculated according to Eq. (5). For a particle incoming from z = −∞, on the lower surfaces U0, there is an energetic barrier to reach the crossing point. When broadening is taken into account (e.g., through the BCME), this barrier is lowered. The arrows show the crossing point with or without broadening.
The CME and BCME surfaces in the z direction. Here x = 0.0 and Γ0 = 0.03. Ui and are the CME and BCME surfaces for , respectively. UPMF is calculated according to Eq. (5). For a particle incoming from z = −∞, on the lower surfaces U0, there is an energetic barrier to reach the crossing point. When broadening is taken into account (e.g., through the BCME), this barrier is lowered. The arrows show the crossing point with or without broadening.
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, px) satisfying equally. In the z direction, trajectories were initialized at position z = −15 and the momentum pz was chosen from a Gaussian distribution with and . We used a time step dt = 0.25 a.u. and propagated trajectories for 2 × 106 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).
Fraction of scattered particles as a function of t. Here x1 = 0.2, Γ0 = 0.03,, Γcutoff = 0.075ω. Note that the CME dynamics result in no trapping, the BCME dynamics result in a modest amount of trapping, and the EF dynamics result in the most trapping (relatively). We have checked that these reported fractions are unchanged for a long time after t = 49 × 104 such that the percentage of trapped trajectories are very meaningful plateau values.
Fraction of scattered particles as a function of t. Here x1 = 0.2, Γ0 = 0.03,, Γcutoff = 0.075ω. Note that the CME dynamics result in no trapping, the BCME dynamics result in a modest amount of trapping, and the EF dynamics result in the most trapping (relatively). We have checked that these reported fractions are unchanged for a long time after t = 49 × 104 such that the percentage of trapped trajectories are very meaningful plateau values.
Scattered particle statistics.
Parameters . | Percentage trapped . | Percentage reflected on state . | |||||
---|---|---|---|---|---|---|---|
x1 . | Γ0 . | . | 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 . | |||||
---|---|---|---|---|---|---|---|
x1 . | Γ0 . | . | 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 . For this reason, the vibrational state of each particle can be calculated as follows: (a) we compute the kinetic energy in the x direction, , (b) we compute the potential energy in the x direction, Epx = U0(x), and (c) we compute nvib = (Ekx + Epx)/ω − 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 U0(x) [rather than ] because at z = −20, Γ → 0, and the U0(x) and the surfaces have negligible differences.
1. Relaxation dependence on
In Fig. 4(a), we plot the vibrational distributions as a function of different Γ0 for fixed incoming z momentum . 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.
Vibrational distribution analysis for the scattered trajectories for different metal-molecule couplings. Here x1 = 0.2 (left), 0.4 (right) [see Eq. (19)]. The incoming z-momentum , Γcutoff = 0.075ω, Γ0 = 0.01, 0.03, 0.05, 0.08. The CME dynamics give no relaxation because without broadening, the trajectories never reach the diabatic crossing point. However, both EF and BCME give relaxation, and the methods agree for large Γ0. Note that both of these methods predict a turnover in relaxation: vibrational relaxation is maximized for Γ0 = 0.05. Note also that, for x1 = 0.4, the BCME gives significantly more relaxation than that of EF when Γ0 is small.
Vibrational distribution analysis for the scattered trajectories for different metal-molecule couplings. Here x1 = 0.2 (left), 0.4 (right) [see Eq. (19)]. The incoming z-momentum , Γcutoff = 0.075ω, Γ0 = 0.01, 0.03, 0.05, 0.08. The CME dynamics give no relaxation because without broadening, the trajectories never reach the diabatic crossing point. However, both EF and BCME give relaxation, and the methods agree for large Γ0. Note that both of these methods predict a turnover in relaxation: vibrational relaxation is maximized for Γ0 = 0.05. Note also that, for x1 = 0.4, the BCME gives significantly more relaxation than that of EF when Γ0 is small.
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 () 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 x1 = 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 x1 = 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 (x1 − x0) 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 with a fixed value of the hybridization (Γ0 = 0.03) and x1 = 0.2. Here, for , 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.
Vibrational distribution analysis for the scattered trajectories for different incoming momenta. Here x0 = 0.2 (left), 0.4 (right), Γ0 = 0.03, Γcutoff = 0.075ω, . The agreement between EF and BCME increases as increases. When x1 = 0.4, more relaxation is predicted for all three schemes.
Vibrational distribution analysis for the scattered trajectories for different incoming momenta. Here x0 = 0.2 (left), 0.4 (right), Γ0 = 0.03, Γcutoff = 0.075ω, . The agreement between EF and BCME increases as increases. When x1 = 0.4, more relaxation is predicted for all three schemes.
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, x1 = 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.
A histogram of hopping energies. Here x1 = 0.2, 0.4, Γ0 = 0.01, 0.05, , Γcutoff = 0.075ω. Positive ΔE means energy transfers from the metal to a trajectory, while negative ΔE means energy transfers from a trajectory to the metal. Note that the y axis has different scales for different subfigures. The plot suggests that, for all 4 cases, most hops have small energy changes, but a large energy transfer is not prohibited. Note the sensitivity of the BCME dynamics to x1.
A histogram of hopping energies. Here x1 = 0.2, 0.4, Γ0 = 0.01, 0.05, , Γcutoff = 0.075ω. Positive ΔE means energy transfers from the metal to a trajectory, while negative ΔE means energy transfers from a trajectory to the metal. Note that the y axis has different scales for different subfigures. The plot suggests that, for all 4 cases, most hops have small energy changes, but a large energy transfer is not prohibited. Note the sensitivity of the BCME dynamics to x1.
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 x1 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 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 x1 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 nvib 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).
BCME surfaces in the x direction for x1 = 0.2 (left), 0.4 (right), with a fixed z = −2.0. Here Γ0 = 0.01. The dotted line is the vibrational energy for particles at nvib = 15, and the shaded areas are the active regions for hopping events (darker) and (lighter), assuming only downward hops (suggested by Fig. 6). Because particles move more slowly in the z direction than that in the x direction when nvib = 15, this cartoon representation (with fixed z position) gives a reasonable explanation for why vibrational relaxation is faster with x1 = 0.4 (as opposed to x1 = 0.2). Obviously, hop1 can be triggered more easily when x1 = 0.4, and hop1 also releases more vibrational energy to the metal for larger x1.
BCME surfaces in the x direction for x1 = 0.2 (left), 0.4 (right), with a fixed z = −2.0. Here Γ0 = 0.01. The dotted line is the vibrational energy for particles at nvib = 15, and the shaded areas are the active regions for hopping events (darker) and (lighter), assuming only downward hops (suggested by Fig. 6). Because particles move more slowly in the z direction than that in the x direction when nvib = 15, this cartoon representation (with fixed z position) gives a reasonable explanation for why vibrational relaxation is faster with x1 = 0.4 (as opposed to x1 = 0.2). Obviously, hop1 can be triggered more easily when x1 = 0.4, and hop1 also releases more vibrational energy to the metal for larger x1.
Now, EF dynamics also predict stronger relaxation for large displacements—after all, the EF friction tensor is proportional to 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. x1 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.
Modified surfaces used in Γcutoff analysis (Fig. 9). x1 = 0.2. The only difference between this figure and Fig. 1 is the shape of T(z), which is relatively sharper here. Here Γ(r) can be very small (but still >Γcutoff) even in the surface crossing region.
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).
Vibrational distribution analysis for the scattered trajectories moving along the modified parameters in Fig. 8 (left) and the original surfaces in Fig. 1 (right). Here x1 = 0.2, Γ0 = 0.03, . For all the Γcutoff parameters used here, the cut-off region [where Γ(r) = Γcutoff] is to the left of the surface crossing point. In the original model, the vibrational distributions from EF are consistent for different Γcutoff parameters, but for the modified parameters, the distributions are quite sensitive to the artificial parameter Γcutoff.
Vibrational distribution analysis for the scattered trajectories moving along the modified parameters in Fig. 8 (left) and the original surfaces in Fig. 1 (right). Here x1 = 0.2, Γ0 = 0.03, . For all the Γcutoff parameters used here, the cut-off region [where Γ(r) = Γcutoff] is to the left of the surface crossing point. In the original model, the vibrational distributions from EF are consistent for different Γcutoff parameters, but for the modified parameters, the distributions are quite sensitive to the artificial parameter Γcutoff.
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 Γ(r) 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.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 IESH25 (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).