Yaolong
Zhang
a,
Reinhard J.
Maurer
*bc,
Hua
Guo
d and
Bin
Jiang
*a
aHefei National Laboratory for Physical Science at the Microscale, Department of Chemical Physics, University of Science and Technology of China, Hefei, Anhui 230026, China. E-mail: bjiangch@ustc.edu.cn
bDepartment of Chemistry and Centre for Scientific Computing, University of Warwick, Gibbet Hill Road, Coventry, CV4 7AL, UK. E-mail: R.Maurer@warwick.ac.uk
cDepartment of Chemistry, Yale University, New Haven, Connecticut 06520, USA
dDepartment of Chemistry and Chemical Biology, University of New Mexico, Albuquerque, New Mexico 87131, USA
First published on 8th November 2018
The breakdown of the Born–Oppenheimer approximation gives rise to nonadiabatic effects in gas-surface reactions at metal surfaces. However, for a given reaction, it remains unclear which factors quantitatively determine whether these effects measurably contribute to surface reactivity in catalysis and photo/electrochemistry. Here, we systematically investigate hot electron effects during H2 scattering from Ag(111) using electronic friction theory. We combine first-principles calculations of tensorial friction by time-dependent perturbation theory based on density functional theory and an analytical neural network representation, to overcome the limitations of existing approximations and explicitly simulate mode-specific nonadiabatic energy loss during molecular dynamics. Despite sizable hot-electron-induced energy loss, no measurable nonadiabatic effects can be found for H2 scattering on Ag(111). This is in stark contrast to previous reports for vibrationally excited H2 scattering on Cu(111). By detailed analysis of the two systems, we attribute this discrepancy to a subtle interplay between the magnitude of electronic friction along intramolecular vibration and the shape of the potential energy landscape that controls the molecular velocity at impact. On the basis of this characterization, we offer guidance for the search of highly nonadiabatic surface reactions.
Given the large number of DOFs in molecule–surface systems, exact quantum nonadiabatic dynamics calculations are currently infeasible to interrogate these systems, yet a variety of approximate models have been developed in recent years to computationally simulate nonadiabatic dissipation during gas-surface dynamics.18 In the weak coupling limit, molecular dynamics with electronic friction (MDEF) models have been widely applied to account for the excitation of low-energy EHPs by the motion of adsorbate molecules.19–22 The electronic friction approach replaces the explicit account of EHPs by coupling of the molecular dynamics with an electronic bath that exerts frictional forces on the adsorbate dynamics within a generalized Langevin dynamics framework.20 This approach is approximate, yet highly efficient as it reduces the description of the EHP spectrum to a single coordinate-dependent quantity, namely the electronic friction tensor. Combined with the local density friction approximation (LDFA),23 which assumes a single friction coefficient for each non-substrate atom depending on the local electron density of the substrate, MDEF simulations have successfully captured the non-adiabatic energy loss of atomic scattering9,24,25 and hot atom relaxations,26,27 at least qualitatively. Within the independent atom approximation (IAA), the LDFA-based MDEF model has also been extended to describe the reactive and non-reactive scattering of molecules,23,28–32 Eley–Rideal reactions,33–35 and laser-induced desorption processes36 at various metal surfaces. In such cases, the electrons in the metal surface are represented by a homogeneous free electron gas within which each atom of the adsorbate is assumed to be embedded. However, the validity of LDFA-IAA has been continuously questioned due to the lack of intrinsic anisotropy in the friction tensor.15,37–41 The latter issue has been partially addressed via the Hirshfeld partitioning of the overall system density to contributions of relevant atoms, which has recently been applied to the relaxation of hot atoms and molecules on metal surfaces.27,38,42
The approximations in the LDFA-IAA model can be overcome using first-order time-dependent perturbation theory (TDPT) that fully accounts for the electronic structure of the interacting molecule–surface system.19,20,43 Because of the high computational costs of calculating the nonadiabatic matrix elements, however, TDPT calculations of electronic friction have been mostly restricted to low-dimensional problems20,37 until very recently. Maurer and coworkers have recently reported an efficient implementation of TDPT based on the Kohn–Sham (KS) orbitals of Density Functional Theory (DFT).40,41 Applications of TDPT to nonadiabatic vibrational relaxation rates for a number of diatomic molecules on metal surfaces revealed the significance of friction-induced mode-coupling.40,41 Based on MD simulations with an on-the-fly calculated TDPT friction tensor, we recently found that tensorial friction induces more energy loss of rovibrational over translational energy in scattering of H2 on Ag(111), resulting in a significant impact on the fate as well as the final energy distributions of individual H2 scattering trajectories.39 Due to the limited number of trajectories, unfortunately, state-to-state scattering probabilities were not obtained. This prediction was later confirmed by Spiering and Meyer for a similar system,15 namely for H2 on Cu(111). These authors overcame the high cost of the on-the-fly friction tensor calculations by using an analytical representation of the friction tensor and obtained an up to six times higher vibrational deexcitation probability of H2(ν = 2, j = 1 → ν = 1, j = 1) upon inclusion of hot electron effects. Interestingly, mode-dependent friction as described by TDPT provides a deexcitation probability that is three times higher than an isotropic description based on LDFA. The latter finding, if experimentally verified, could serve as an important experimental fingerprint to evidence dynamical signatures of the mode-selectivity and anisotropy of electronic friction effects.35
In this work, we report state-to-state scattering and dissociation probabilities of vibrationally excited H2 and D2 from Ag(111) using MDEF with both the TDPT and LDFA models. To this end, we develop an alternative and more general analytical representation of the friction tensor as a function of Cartesian coordinates using neural networks (NNs). Interestingly, with a larger number of trajectories on an analytical potential energy surface (PES), the MDEF simulations with TDPT and with LDFA both yield scattering-induced vibrational deexcitation probabilities of H2(ν = 2, j = 0 → ν = 1, j = 0) that are not significantly increased compared to the adiabatic case, which is in stark contrast to findings on Cu(111). Deeper analysis reveals large differences in electronic friction between the two methods and between the two substrates (Ag vs. Cu). We find that, for nonadiabatic effects to significantly affect experimental measurements, the velocity profile as controlled by the PES is equally important as the magnitude of electronic friction with respect to the intramolecular stretching motion. The strong interplay between the two quantities is in fact another manifestation of mode-selective electronic friction and the analysis of the two can provide guidance on whether nonadiabatic effects can produce quantum state-resolved measurable differences in experimental attributes.
(1) |
Herein, the three force components on the right side of eqn (1) 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 vector of the atoms Ṙ, and a temperature- and friction-dependent random force term that restores detailed balance. To arrive at a GLE, one must assume that the electron-nuclear coupling acts instantaneously, which equates to the Markov approximation. The result is that, within 1st order TDPT, we evaluate the electronic friction tensor using Fermi's golden rule:40,43
(2) |
Herein, the factor 2 accounts for spin multiplicity in the case of non-spin-polarised calculations and ψkν and εkν correspond to KS eigenvectors and eigenstates. In eqn (2), we have not explicitly addressed 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 between k states to obtain a continuous electronic structure representation giving rise to friction tensor elements, which correspond to relaxation rates due to electron-nuclear coupling along the (mixed) Cartesian directions i and j. We are using local atomic orbital basis functions and evaluate nonadiabatic coupling matrix elements by numerical finite difference in terms of Hamiltonian and overlap matrices in real-space basis representation that are subsequently transformed into k-space. As detailed in our previous work39 and by Spiering and Meyer,15 we do this by replacing the delta function with a normalized Gaussian function of finite width 0.6 eV. It has recently been pointed out that such broadening can introduce contributions that are not rigorously contained in a zero-frequency 1st order approximation within many-body perturbation theory44,45 and leads to relaxation rates that depend on the choice of broadening. This pragmatic choice has been discussed in detail in our previous work39 and further methodological developments to overcome the existing restrictions of Markovian MDEF with TDPT are in progress. All elements of the friction tensor are calculated as a function of adsorbate atom position using the all-electron, local atomic orbital code FHI-aims46 and the PBE functional.47 Our computational settings regarding model set-up, basis set, and Brillouin zone sampling have been detailed in our previous publication.39 In practice, we carry out MDEF calculations with the quasi-classical trajectory (QCT) method as detailed elsewhere48 and briefly described in the ESI.† Adiabatic QCT calculations for reactive and nonreactive scattering of H2 on several metal surfaces have been validated before to reproduce almost quantitatively the quantum results.49,50 We neglect the random force term in eqn (1) considering the short collision time in this process which will simplify the discussion.15
Nevertheless, the TDPT-based friction tensor is directionally-dependent on all coordinates of the molecule. Consequently, a proper transformation matrix must exist to link the correlated friction tensors of two symmetry equivalent molecular geometries due to the intertwined translational, rotational, and permutational symmetries. In this work, we explicitly account for the complicated covariance properties of the friction tensor with respect to the surface symmetry by a simple mapping scheme. We transform the original data of coordinates and friction tensor elements into an irreducible triangle in the surface unit cell, as depicted in Fig. 1b. The molecular geometry and friction tensor thus form a symmetry unique one-to-one mapping in this triangle, which can be represented by NNs. Since NNs are analytical functions of the nuclear coordinates, the calculation of the interpolated friction tensor at any given geometry can be computed efficiently. Our approach and the interpolation procedure are briefly described in the SI and will be more thoroughly discussed in a forthcoming publication.56
Fig. 1 (a) Internal (X, Y, Z, r, θ, ϕ) and Cartesian (x1, y1, z1, x2, y2, z2) coordinates used in the H2 + Ag(111) system, and MEP following red arrows along which the molecule dissociates from a bridge site to two hollow sites. (b) Mirror reflections that can move a molecular geometry outside the irreducible triangle (red) into the triangle. U1/U1−1, U2/U2−1, or U3/U3−1 corresponds to the transformation/inverse transformation matrix for each reflection, respectively. (c) Diagonal terms of nonadiabatic relaxation rates (effective friction coefficients) obtained by NNs (solid symbols) and TDPT (open symbols) along the MEP in terms of internal coordinates. (d) Comparison of the total, internal and translational energies as a function of time for an H2(ν = 1) molecule scattering vertically from a hollow site with a translational incidence energy of 0.6 eV (ref. 39) with the on-the-fly MDEF simulation (OTF, solid lines) and those with the simulation based on the interpolated friction tensor (NN, dotted lines). |
Despite the lack of EHP effects on the sticking probability, Spiering and Meyer found significant EHP effects for the vibrational de-excitation probabilities. Following their work, we analyse the de-excitation probabilities of scattered molecules H2 (D2) from the (ν = 2, j = 0) state to (ν = 1, j = 0) state (P2→1) as a function of incident energy via the standard histogram binning scheme. As shown in Fig. 2b and c, P2→1 shows an initial increase followed by a gradual decrease with the increasing Ei, peaking at Ei = 0.4 (0.5) eV by ∼1.5%. This dependence and magnitude of P2→1 are similar to what have been observed in H2 (D2) scattering from Cu(111).15 The vibrational deexcitation probability measures how effectively adsorbate vibrational energy is lost during molecular scattering at a given molecular impact velocity. The peak at 0.4 eV indicates that vibrational energy loss is most effective at intermediate kinetic energies. One may notice that the translational energy thresholds for both the dissociation and vibrational deexcitation of H2 and D2 differ by roughly 0.2 eV, which may be attributed to the difference of the vibrationally adiabatic barriers for the two isotopologues. This suggests that the vibrational deexcitation possibly occurs when the molecule has sufficient energy to climb up from the entrance channel to the transition state region but fails to dissociate. H2(ν = 2) carries more vibrational energy than D2(ν = 2) and hence requires less translational energy to access the high density and repulsive region, resulting in a lower threshold for vibrational deexcitation. However, surprisingly, we do not find a noticeable increase of P2→1 when accounting for nonadiabatic energy loss with MDEF compared to the adiabatic results, as reported in the study of Spiering and Meyer.15 This is true for both tensorial TDPT friction and isotropic LDFA friction. In contrast, for H2 (D2) scattering from Cu(111), the vibrational deexcitation probabilities in the presence of tensorial friction were found to increase by a factor of up to 6 (3) for H2 and 3 (1.5) for D2 when compared to the adiabatic (LDFA) ones.15 Surprisingly, the seemingly similar H2 + Cu(111) and H2 + Ag(111) systems behave quite differently with respect to the EHP effects for molecular scattering. We note in passing that this small nonadiabatic effect in our case is independent on the final rotational states and scattering angles (see Fig. S2†).
How can the effects of electronic friction on the vibrational deexcitation probability lead to a significant change of a measurable observable in one case and leave no measurable trace in the other? In the adiabatic case, the possibility of vibrational-to-translational energy transfer largely determines the vibrational deexcitation probability, since the total energy is conserved. When accounting for the nonadiabatic energy loss induced by electronic friction, vibrational energy can be dissipated to surface EHPs in addition to being transferred to translation, which increases the vibrational deexcitation probability. In this regard, the larger the vibrational-to-electronic energy dissipation, the larger the vibrational deexcitation probability is expected to be.
Following this line of thought, we first compare the most important components of the electronic friction tensor computed by the LDFA and TDPT methods, in terms of internal coordinates along the MEPs for H2 dissociative chemisorption on Cu(111) and Ag(111). The LDFA friction coefficients Λrr and ΛZZ shown in Fig. 3a are more or less the same for both metal substrates (note that the LDFA model does contain mode–mode coupling between internal coordinates, and ΛrZ values are close to zero here because of the parallel orientation of H2 along the minimum energy path). On the other hand, Fig. 3b presents TDPT results for Λrr, ΛZZ, and ΛrZ. Remarkably, although these three coefficients exhibit similar trends and peak at transition states for both systems, the absolute values of Λrr, ΛZZ, and ΛrZ for H2 + Cu(111) are more than twice of those for H2 + Ag(111).
Fig. 3 Comparison of several LDFA (a) and TDPT (b) friction coefficients in internal coordinates (Λrr, ΛZZ, and ΛrZ) for the H2 + Cu(111)15 (dashed–dotted curves) and H2 + Ag(111) (solid curves) systems. Note that the data for H2 + Cu(111) are extracted from Fig. 2 in ref. 15, and the reaction coordinates (s) are defined by the arc length in the (Z, r) plane starting from the transition state (s = 0) for both systems. |
Upon careful comparison of our computational approach and numerical settings with the calculations in ref. 15, we conclude that this difference is rooted in the differences in the underlying electronic structure of the two metals. In comparison of the two models, for H2 + Ag(111), LDFA gives up to 6 times larger friction coefficients with respect to translation parallel to the surface normal (ΛZZ) than TDPT, while both models predict a similar magnitude of the friction coefficient with respect to H–H vibration (Λrr). Interestingly, for H2 + Cu(111), the overall magnitudes of TDPT friction elements are quite large, with Λrr and ΛZZ being much larger and comparable to the counterparts based on LDFA, respectively. The difference between the TDPT and LDFA based electronic friction models can be largely attributed to their different physical contents. Whereas LDFA only depends on the magnitude of the electron density of the clean substrate – a quantity that is similar along the MEPs for both substrates, TDPT captures the electronic structure differences between the two molecule–surface systems.
The above analysis is based on the comparison of the friction tensor elements of the two systems along the global MEP on which the molecular center is more or less over the bridge site. However, the molecules could be scattered upon the impact at other sites in the unit cell. In Fig. S3a,† we show the distribution of initial lateral positions of the scattered trajectories with H2(ν = 2, j = 0) to H2(ν = 1, j = 0). Given the weak steering effect observed previously,55 it is clear that the vibrational deexcitation process occurs with higher probabilities near the top site than other places, due presumably to the lower reactivity and stronger repulsion there. We therefore include the friction tensor analyses for the reaction paths at fixed top, bridge, and fcc sites for H2 + Ag(111) in Fig. S3b and c.† Interestingly, the results are qualitatively similar for the three paths, along which the absolute values of TDPT-based friction tensor elements gradually increase up to the barrier region then go down. The molecules that follow the top site path could not reach too closely to the surface atoms (high density region) because of the strongly repulsive potential, resulting in the smallest friction tensor values. These additional results further suggest a minor role of EHPs in the H2 + Ag(111) case.
These findings on electronic friction lead us to a better understanding of the mode-specific nonadiabatic energy loss, which naturally explains our observations. In Fig. 4a, the total nonadiabatic energy losses as predicted by TDPT and LDFA for the state-to-state H2(ν = 2, j = 0 → ν = 1, j = 0) scattering on Ag(111) are compared. The LDFA energy loss is on average greater than that with the TPDT model using the same initial conditions. Interestingly, LDFA induces increasingly stronger energy loss with increased translational energy than TDPT, due apparently to the much higher ΛZZ values with LDFA along the MEP seen in Fig. 3 and S3.† This difference can be more clearly seen in Fig. 4b by comparing the mode-dependent energy losses due to EHPs, which can be separately evaluated by the integration of the friction force,27
(3) |
Fig. 4 (a) Average non-adiabatic energy loss with the TDPT and LDFA models as a function of translational energy for H2 scattering from Ag(111) and Cu(111) (approximately converted from Fig. 4 in ref. 15 by energy conservation). (b) Mean energy losses along translational (Z) and vibrational coordinate (r) as a function of translational energy for H2. The variation of square of velocity (solid lines) in Z (|Ṙz|2, panel c) and r (|Ṙr|2, panel d) direction along with the corresponding friction coefficients (dotted lines), as a function of time for a representative trajectory with the same initial conditions at translational energy of 0.3 eV. |
On the contrary, the amount of nonadiabatic energy loss for H2 + Cu(111), which has been extracted from ref. 15 and displayed in Fig. 4a, is much larger and the pattern is quite different. We find in that case that the TDPT induced energy loss is on average about 0.30 eV, which is twice as large than that for LDFA due mainly to the larger value of the Λrr component with TDPT. This large vibrational energy loss due to EHPs certainly affects the vibrational deexcitation probabilities, and as a result, P2→1 of H2 with TDPT is about 2–3 times higher than that with LDFA, both of which prove a significant increase compared to the adiabatic results.
It is important to note that the energy loss in eqn (1) arises from the product of electronic friction and the velocity of the particle. Therefore, the underlying potential energy landscape and the velocity profile along the scattering trajectories also play an important role. One notices that the dissociation barrier locations are different for H2 + Ag(111) and H2 + Cu(111). The transition state of H2 dissociation on Ag(111) is more “product-like”,53 which features a longer H–H distance (r = 1.26 Å) and higher barrier (Ea = 1.16 eV), compared to r = 1.03 Å and Ea = 0.63 eV on Cu(111).50 The large barrier and the late transition state for H2 dissociation on Ag(111) could cause strongly reduced velocities of atoms in a large region around the turning point, as illustrated in Fig. 4c and d for a representative scattering trajectory. In Fig. 4c, the velocity along Z direction shows a sharp dip as the molecule approaches the point of reflection at the surface (∼110 fs) where the corresponding friction coefficients are the largest for both TDPT and LDFA models. A similar reduction of velocity and increase of friction coefficients along the intramolecular stretching DOF can be also seen in the same region in Fig. 4d, although the velocity is highly oscillating as the molecule vibrates. Interestingly, the TDPT based Λrr also oscillates out of phase with the velocity. This is indeed a signature of the mode-specific electronic friction varying with the molecular structure, which is absent in the LDFA-IAA friction. As a result, even when the electronic friction near the transition state region close to the surface is quite large, the H2 molecule has a low velocity, which effectively diminishes the electronic friction effects. The nonadiabatic energy loss, i.e. the integral over the friction force that can be calculated by eqn (3), is therefore very small relative to the total energy of the molecule. The potential energy landscape during scattering, therefore, plays a critical role in enabling or inhibiting large nonadiabatic energy loss by its effect on the average velocity profile during a scattering event. In case of H2 on Ag(111), the velocity profile and friction profile are such that the effect on the dynamic deexcitation probabilities is small. On Cu(111), in contrast, it is expected that the molecule can more easily reach the region of high electron density, due the lower barrier, giving rise to much larger energy loss as seen in Fig. 4a and a more remarkable increase of P2→1. The important role of this interplay between the velocity and friction coefficients has been remarked by Juaristi et al. when using LDFA to understand EHP effects for dissociative chemisorption.23 We further emphasize in the present work that this interplay depends on the mode-specific friction coefficients and corresponding velocity profile, which is the key to understand the role of electronic excitation in state-to-state scattering of molecules from metal surfaces.
Our results emphasize the importance of simultaneously accounting for the mode-dependent magnitude and the tensorial nature of electronic friction and the underlying molecule–surface interaction as described by a realistic ab initio calculated PES. The latter is much more system-dependent, which means that strong intrinsic electronic friction components during a gas-surface reaction do not a priori lead to measurable electronic dissipation effects. For such measurable effects, an interplay between high adsorbate velocities along DOFs and strong electronic friction is required and should be a guiding principle when selecting systems with high propensity for nonadiabatic effects. To arrive at a clearer picture of nonadiabatic dissipation during gas-surface dynamics, and, more importantly, to design systems where nonadiabatic effects can be utilized to selectively control chemical reactions as envisioned in hot-electron chemistry, this interplay needs to be studied systematically in the future.
Footnote |
† Electronic supplementary information (ESI) available. See DOI: 10.1039/c8sc03955k |
This journal is © The Royal Society of Chemistry 2019 |