Open Access Article

This Open Access Article is licensed under a

Creative Commons Attribution 3.0 Unported Licence

DOI: 10.1039/C8FD00140E
(Paper)
Faraday Discuss., 2019, Advance Article

Reinhard J. Maurer†
*^{a},
Yaolong Zhang†^{b},
Hua Guo^{c} and
Bin Jiang*^{b}
^{a}Department of Chemistry, Centre for Scientific Computing, University of Warwick, Gibbet Hill Road, Coventry, CV4 7AL, UK. E-mail: R.Maurer@warwick.ac.uk
^{b}Hefei 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
^{c}Department of Chemistry and Chemical Biology, University of New Mexico, Albuquerque, New Mexico 87131, USA

Received
1st October 2018
, Accepted 30th October 2018

First published on 10th December 2018

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

The proposed technological applications of hot electron effects in plasmonic sensors,^{12} plasmon-enhanced catalysts,^{13} and locally-enhanced cancer radiation treatment,^{14} are driven by exciting fundamental experiments that provide quantitative evidence for hot electron effects. These include molecular beam scattering experiments that explicitly measure the state-to-state scattering probability and the energy loss profile of a gas–surface reaction,^{15} but also the measurement of electric currents,^{11,16} so-called chemicurrents, that occur due to hot electrons that are generated upon adsorption, desorption, or chemical transformation of molecules.

The ability to computationally predict the molecular details that underlie hot electron-mediated processes and the intricate coupling of electrons and atomic motion is essential for the design of new devices and catalysts based on plasmonics. As such, it is imperative to develop first-principles computational simulation methods that correctly capture how molecular adsorbate motion generates hot electrons or how hot electrons affect measurable dynamical properties of chemical reactions. The intrinsic dimensionality of systems involving many hundreds of atoms and thousands of electrons make exact nonadiabatic quantum dynamics simulations close to unfeasible and therefore ask for efficient approximate mixed quantum-classical dynamics methods.^{17,18} Such methods are typically based on stochastic surface hopping models based on chemical master equations, where hot electron effects on the reaction dynamics are modelled as explicit electronic transitions.^{19,20} An alternative is the molecular dynamics with electronic friction (MDEF) method,^{21–23} where the classical motion of nuclei follows a single potential energy surface and hot electron effects are described as electronic friction forces that act on the adsorbate nuclei. These forces are governed by the so-called electronic friction tensor,^{24} which arises from nonadiabatic coupling and the resulting equations of motion follow a stochastic Langevin model. The electronic friction tensor is in general frequency-dependent and anisotropic but is often approximated to be frequency independent and spatially isotropic, i.e., nonadiabatic relaxation rates are equal in all spatial directions for rigid atomic and molecular motion. The most common is the Local Density Friction Approximation (LDFA),^{25,26} where nonadiabatic coupling is expressed as a simple function of the metallic electron density within which the adsorbate atom is embedded. This model is spatially isotropic for rigid molecular motion, but does account for some anisotropic effects along intramolecular degrees of freedom.^{26,27} Due to its simplicity and ease of computation, LDFA-based MDEF is a particularly appealing approach that, in the past, has been applied to describe hot electron-mediated dynamics in ambient^{28} and laser-driven conditions.^{29,30} However, especially for molecular scattering and light-driven plasmonics, the underlying neglect of frequency dependence and the role of molecular anisotropy of electronic friction forces and their effects on measurable observables have been the reason for exciting debates in recent literature.^{24,31–33}

In recent work, we have achieved a description of electronic friction that captures the directional and the explicit vibrational mode dependence and tensorial nature of electronic friction using an efficient implementation of first-order time-dependent perturbation theory (TDPT) based on Density Functional Theory (DFT).^{24,31} On the examples of vibrational energy loss of metal-adsorbed diatomics^{31} and the dissociative chemisorption of molecular hydrogen on Ag(111),^{34} we were able to show that the overall magnitude of friction is, on average, similar to the LDFA prediction, but the directional and vibrational mode dependences and their effect on the dynamic energy loss are different. Before us, Spiering and Meyer have reached a similar conclusion for reactive scattering of H_{2} on Cu(111),^{35} where mode-dependent electronic friction forces as described by TDPT provide a measurable signature of energy loss during state-to-state scattering of vibrationally excited molecules. We have recently developed a general method to construct an analytical representation of the electronic friction tensor as a function of the adsorbate coordinates that enabled us to simulate dynamic scattering probabilities and to discuss the subtle differences in electronic friction that give rise to the observation that hot electron effects are much more pronounced during reactive scattering of H_{2} on Cu(111) than on Ag(111).^{36}

In this work, we use this analytical representation of electronic friction for H_{2} on Ag(111) further to systematically assess the validity and robustness of the current frequency-independent MDEF method based on TDPT electronic friction. We do this by studying the sensitivity of dynamic observables such as the vibrational excitation/de-excitation probability during molecular scattering and the integrated energy loss along a specific degree of freedom as a function of various initial conditions and the type of electronic friction description. The former includes the effects of the initial and final vibrational excited states and the electronic temperature of the substrate. The latter includes a comparison between LDFA and TDPT and a discussion of their limitations and sensitivity to intrinsic numerical parameters. We find that the energy loss predicted by LDFA and TDPT shows different trends with respect to the parameters we studied and that dynamic observables are particularly sensitive to the absolute magnitude of tensorial friction. The intrinsic sensitivity of friction to numerical parameters such as the broadening width will represent a particular challenge for quantitative predictions of hot electron mediated reaction dynamics at high electronic temperatures that correspond to light-driven conditions. We therefore stress the importance of overcoming this limitation introduced by assuming frequency independence of the friction tensor. This will enable further study of the applicability of electronic friction approaches for light-driven hot electron chemistry.

(1) |

(2) |

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

where T(θ, φ) is the transformation matrix that rotates configuration (θ, φ) to the corresponding reference geometry with θ_{0} = φ_{0} = 90°. Although this decomposition of the full 6 × 6 friction tensor certainly depends on the choice of the reference orientation, Spiering and Meyer have shown that the dependence is relatively weak near the minimum energy path (MEP).

where q and Λ are the original coordinate vector and friction tensor and U is the transformation matrix. For Ag(111), U_{1}, U_{2} and U_{3} are given below and illustrated in Fig. 1,

Spiering and Meyer have very recently reported an approximate NN-based representation of the 6 × 6 friction tensor for the H_{2} + Cu(111) system.^{35} They first described individual friction coefficients in terms of four symmetry independent internal coordinates, as illustrated in Fig. 1, namely X, Y, Z, and r with angular coordinates fixed (e.g., θ_{0} = φ_{0} = 90°). Since the tensorial friction (Λ) is a symmetric matrix, only the upper triangular elements Λ_{ij}(i ≤ j) are to be determined in practice. This treatment avoids the complex symmetry involved in the full six-dimensional representation. Once Λ(X, Y, Z, r) with θ_{0} = φ_{0} = 90° is obtained, the friction tensor at an arbitrary configuration depending on all six internal coordinates can be approximately generated by the following transformation,^{35}

Λ(X, Y, Z, r, θ, φ) ≈ T(θ, φ)Λ(X, Y, Z, r)T^{−1}(θ, φ),
| (3) |

In a previous publication^{36} and in this present work, we propose an alternative strategy. We choose to fit the elements of the full 6 × 6 matrix form of the electronic friction tensor in Cartesian coordinates inside an irreducible triangle in the (111) surface unit cell (marked in red in Fig. 1). In this triangle, each friction tensor element is symmetry unique with a one-to-one mapping to a molecular geometry, thus posing no problems for fitting. In practice, since we have collected many data scattered inside and out the unit cell, we first move those data distributed outside this triangle into it. This can be done by first translating a data point outside the unit cell into this region, which does not change the friction tensor, followed by a series of reflections when necessary.^{46} For each reflection, the new Cartesian coordinate vector q′ and friction tensor Λ′ can be obtained by,

q′ = Uq, | (4) |

Λ′ = U^{T}ΛU,
| (5) |

(6) |

(7) |

(8) |

In practice, the number of reflections operated on a given point may be less than three and, in the case that the molecule is already situated inside the triangle, no transformation is necessary. It should be noted that this treatment does not account for exchange symmetry of the two H atoms. We solve this issue by ordering the positions of the two hydrogen atoms when necessary (i.e. with z_{1} ≤ z_{2}) so that the related matrix elements in the friction tensor are exchanged accordingly. Similar mapping schemes have been successfully applied in constructing PESs for gas–surface systems using either NNs^{47–49} or the modified Shepard interpolation.^{46} One drawback of these PESs is that the gradients of PESs (i.e. forces) may not be continuous at the boundary, resulting in non-conservation of total energy when running classical trajectories. However, this is not a problem for the friction tensor as we only need values of friction coefficients but no gradients in eqn (1). Similar schemes can be devised for other facets of the metal surface.

2.2.2 Fitting and evaluation. After the mapping procedure, a single NN is utilized to describe the one-to-one mapping between the molecular geometry and each individual element of the friction tensor in the triangle. Six Cartesian coordinates of the molecule are set as the input layer of the NN, which pass geometric information to two hidden layers with 20 and 40 neurons via hyperbolic tangent activation functions. A hybrid algorithm combining the extreme learning machine and the Levenberg–Marquardt (ELM–LM) algorithm was used to train the NN.^{50,51} The data set consists of ∼1200 points sampled by previous AIMD trajectories^{34} and ∼2600 additional points used for constructing the PES^{44} that cover the dynamically relevant region. Since the friction coefficients are very small far above the surface where the electron density vanishes, we exclude data points in that region. Instead, a switching function was imposed to smoothly dampen the friction coefficients to zero when the molecule–surface interaction is negligible (Z → ∞, in Å),

(9) |

One should note that the NN fitting of each friction coefficient is essentially the same as that of the PES in previous works and very robust.^{47–49} However, for our purpose, to evaluate a friction tensor for any point outside this symmetry-unique region, we need to first find the symmetry identical geometry of that point in the irreducible triangle by the transformation described above, then calculate the NN interpolated friction tensor there, followed by inversely transforming the obtained full 6 × 6 friction tensor back to that corresponding to the original geometry.

Given an initial condition, eqn (1) is propagated using third-order Beeman’s algorithm with a time step of 0.1 fs. The MD and MDEF with LDFA simulations are easily implemented as friction coefficients only depend on the value of the electron density of the clean substrate.^{25,28} The tensorial nature of the TDPT based electronic friction makes the equations of all coordinates coupled. In practice, we need to transform the coordinates and velocities into the friction eigenspace and apply the friction eigenvalues as damping coefficients onto the transformed velocities, and then transform the products back to the Cartesian space. This requires a diagonalization of the friction tensor at each time step. In addition, the random force term on the right-hand side of the Langevin equation, given a finite surface temperature (T), is described as a Gaussian white noise of zero mean and a standard deviation of , where η_{i} is the ith friction eigenvalue (or simply the ith friction coefficient in LDFA) and k is the Boltzmann constant. However, we neglect surface temperature in most simulations because the contribution of the random force is typically small, except where we discuss temperature effects explicitly.

Furthermore, it is interesting to analyse whether the nonadiabatic effects stem from the anisotropy along different directions or the mode-coupling between these directions in the friction tensor. In Fig. 4c, we compare the P_{2→1} with TDPT electronic friction where the off-diagonal coupling between the intra-atomic coordinate and the molecule-adsorbate coordinate (TDPT) is included and where it is artificially set to zero (TDPT-Λ_{rZ=0}). To make the comparison clearer, a factor of two is also multiplied in both cases. Although the magnitude of Λ_{rZ} is not small as displayed in Fig. 2b, interestingly, the results with zero Λ_{rZ} are quite close to the original ones in the entire energy range, implying the negligible effect of the mode-coupling between the translational and vibrational DOFs in this case. In addition, to study the importance of directional anisotropy, we average the three diagonal components of the TDPT friction tensor (twice in magnitude) in the Cartesian space for each atom, then assign the same friction coefficient for each atomic Cartesian coordinate in the following way,

(10) |

Following this line of thought, we show the excitation probability (P_{2→3}) from H_{2} (v = 2, j = 0) to H_{2} (v = 3, j = 0) in Fig. 5b and deexcitation probability (P_{3→2}) from H_{2} (v = 3, j = 2) to H_{2} (v = 2, j = 2) in Fig. 5c. Both scattering processes involve a higher vibrational state which naturally increases the molecular velocity with respect to the intramolecular stretching mode. We see much stronger nonadiabatic effects in both cases, leading to the significant decrease of P_{2→3} and increase of P_{3→2} due to electronic friction at low translational energies. In such cases, the translational velocity is low enough for the molecule to remain close to the surface for an extended time, and the vibrational velocity is high enough to cause electronic excitations. While the results with TDPT and LDFA are very similar with respect to P_{2→3}, the former become qualitatively different from the latter in P_{3→2} below E_{i} = 0.2 eV. This interesting behavior is reflected by the nonadiabatic energy loss given in Fig. 5d, where the TDPT induced energy losses along the molecular vibration (r direction) are more significant than the LDFA ones in that range. Such a large vibrational energy loss with TDPT at low energy seems to be surprising given the similar magnitude of Λ_{rr} with TDPT and LDFA. We therefore monitor the velocity and friction coefficient profiles in the Z and r directions as a function of time for an exemplary trajectory at E_{i} = 0.1 eV in Fig. 6a and b. It appears that Λ_{rr} with TDPT is remarkably higher than those with LDFA in the strong interaction region, implying that the trajectory significantly departs from the MEP. Unlike our previous results (Fig. 4d in ref. 36) for H_{2} (v = 2 → v = 1), the H_{2} (v = 3) molecules do not necessarily follow the MEP. As a result, the multidimensional topography of the friction tensor determines the total nonadiabatic energy loss and thus the scattering probability. This reflects the importance of correctly capturing this tensorial topology.

- M. L. Brongersma, N. J. Halas and P. Nordlander, Nat. Nanotechnol., 2015, 10, 25–34 CrossRef CAS PubMed.
- W. Xie and S. Schlücker, Nat. Commun., 2015, 6, 7570 CrossRef PubMed.
- J. Y. Park, S. W. Lee, C. Lee and H. Lee, Catal. Lett., 2017, 147, 1851–1860 CrossRef CAS.
- J. Y. Park, S. M. Kim, H. Lee and I. I. Nedrygailov, Acc. Chem. Res., 2015, 48, 2475–2483 CrossRef CAS PubMed.
- S. D. Forno, L. Ranno and J. Lischner, J. Phys. Chem. C, 2018, 122, 8517–8527 CrossRef.
- E. Cortés, W. Xie, J. Cambiasso, A. S. Jermyn, R. Sundararaman, P. Narang, S. Schlücker and S. A. Maier, Nat. Commun., 2017, 8, 14880 CrossRef PubMed.
- S. Mukherjee, F. Libisch, N. Large, O. Neumann, L. V. Brown, J. Cheng, J. B. Lassiter, E. A. Carter, P. Nordlander and N. J. Halas, Nano Lett., 2013, 13, 240–247 CrossRef CAS PubMed.
- A. M. Wodtke, D. Matsiev and D. J. Auerbach, Prog. Surf. Sci., 2008, 83, 167–214 CrossRef CAS.
- S. P. Rittmeyer, V. J. Bukas and K. Reuter, Adv. Phys.: X, 2018, 3, 1381574 Search PubMed.
- A. M. Wodtke, J. C. Tully and D. J. Auerbach, Int. Rev. Phys. Chem., 2004, 23, 513–539 Search PubMed.
- A. Kandratsenka, H. Jiang, Y. Dorenkamp, S. M. Janke, M. Kammler, A. M. Wodtke and O. Bünermann, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, 680–684 CrossRef CAS PubMed.
- D. Sil, K. D. Gilroy, A. Niaux, A. Boulesbaa, S. Neretina and E. Borguet, ACS Nano, 2014, 8, 7755–7762 CrossRef CAS PubMed.
- P. Christopher, H. Xin and S. Linic, Nat. Chem., 2011, 3, 467–472 CrossRef CAS PubMed.
- X. Huang, P. K. Jain, I. H. El-Sayed and M. A. El-Sayed, Nanomedicine, 2007, 2, 681–693 CrossRef CAS PubMed.
- O. Bunermann, H. Jiang, Y. Dorenkamp, A. Kandratsenka, S. M. Janke, D. J. Auerbach and A. M. Wodtke, Science, 2015, 350, 1346–1349 CrossRef PubMed.
- H. Nienhaus, H. S. Bergh, B. Gergen, A. Majumdar, W. H. Weinberg and E. W. McFarland, Phys. Rev. Lett., 1999, 82, 446–449 CrossRef CAS.
- J. C. Tully, J. Chem. Phys., 2012, 137, 22A301 CrossRef PubMed.
- Y. Zhang, T. R. Nelson, S. Tretiak, H. Guo and G. C. Schatz, ACS Nano, 2018, 12, 8415–8422 CrossRef CAS PubMed.
- N. Shenvi, S. Roy and J. C. Tully, J. Chem. Phys., 2009, 130, 174107 CrossRef PubMed.
- H. Guo, P. Saalfrank and T. Seideman, Prog. Surf. Sci., 1999, 62, 239–303 CrossRef CAS.
- M. Head Gordon and J. C. Tully, J. Chem. Phys., 1995, 103, 10137–10145 CrossRef CAS.
- W. Dou and J. E. Subotnik, J. Chem. Phys., 2018, 148, 230901 CrossRef PubMed.
- B. Hellsing and M. Persson, Phys. Scr., 1984, 29, 360–371 CrossRef CAS.
- M. Askerka, R. J. Maurer, V. S. Batista and J. C. Tully, Phys. Rev. Lett., 2016, 116, 217601 CrossRef PubMed.
- J. Juaristi, M. Alducin, R. Muiño, H. Busnengo and A. Salin, Phys. Rev. Lett., 2008, 100, 116102 CrossRef CAS PubMed.
- D. Novko, M. Blanco-Rey, M. Alducin and J. I. Juaristi, Phys. Rev. B, 2016, 93, 245435 CrossRef.
- S. P. Rittmeyer, J. Meyer and K. Reuter, Phys. Rev. Lett., 2017, 119, 176808 CrossRef PubMed.
- M. Blanco-Rey, J. I. Juaristi, R. Díez Muiño, H. F. Busnengo, G. J. Kroes and M. Alducin, Phys. Rev. Lett., 2014, 112, 103203 CrossRef CAS PubMed.
- J. I. Juaristi, M. Alducin and P. Saalfrank, Phys. Rev. B, 2017, 95, 125439 CrossRef.
- I. Lončarić, M. Alducin, P. Saalfrank and J. I. Juaristi, Phys. Rev. B, 2016, 93, 014301 CrossRef.
- R. J. Maurer, M. Askerka, V. S. Batista and J. C. Tully, Phys. Rev. B, 2016, 94, 115432 CrossRef.
- D. Novko, M. Alducin, M. Blanco-Rey and J. I. Juaristi, Phys. Rev. B, 2016, 94, 224306 CrossRef.
- D. Novko, M. Alducin and J. I. Juaristi, Phys. Rev. Lett., 2018, 120, 156804 CrossRef CAS PubMed.
- R. J. Maurer, B. Jiang, H. Guo and J. C. Tully, Phys. Rev. Lett., 2017, 118, 256001 CrossRef PubMed.
- P. Spiering and J. Meyer, J. Phys. Chem. Lett., 2018, 9, 1803–1808 CrossRef CAS PubMed.
- Y. Zhang, R. J. Maurer, H. Guo and B. Jiang, Chem. Sci., 2019, 10, 1089–1097 RSC.
- S. Roy, N. A. Shenvi and J. C. Tully, J. Chem. Phys., 2009, 130, 174716 CrossRef PubMed.
- N. Shenvi, S. Roy and J. C. Tully, Science, 2009, 326, 829–832 CrossRef CAS PubMed.
- T. Olsen and J. Schiøtz, J. Chem. Phys., 2010, 133, 134109 CrossRef PubMed.
- M. Forsblom and M. Persson, J. Chem. Phys., 2007, 127, 154303 CrossRef PubMed.
- A. C. Luntz and M. Persson, J. Chem. Phys., 2005, 123, 074704 CrossRef CAS PubMed.
- V. Blum, R. Gehrke, F. Hanke, P. Havu, V. Havu, X. Ren, K. Reuter and M. Scheffler, Comput. Phys. Commun., 2009, 180, 2175–2196 CrossRef CAS.
- J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868 CrossRef CAS PubMed.
- B. Jiang and H. Guo, Phys. Chem. Chem. Phys., 2014, 16, 24704–24715 RSC.
- B. Jiang, M. Alducin and H. Guo, J. Phys. Chem. Lett., 2016, 7, 327–331 CrossRef CAS PubMed.
- C. Crespos, M. A. Collins, E. Pijper and G. J. Kroes, J. Chem. Phys., 2004, 120, 2392–2404 CrossRef CAS PubMed.
- T. Liu, B. Fu and D. H. Zhang, Sci. China: Chem., 2014, 57, 147–155 CrossRef CAS.
- T. Liu, Z. Zhang, B. Fu, X. Yang and D. H. Zhang, Chem. Sci., 2016, 7, 1840–1845 RSC.
- X. Shen, J. Chen, Z. Zhang, K. Shao and D. H. Zhang, J. Chem. Phys., 2015, 143, 144701 CrossRef PubMed.
- Y. Zhang, X. Zhou and B. Jiang, Chin. J. Chem. Phys., 2017, 30, 727–734 CrossRef CAS.
- M. T. Hagan and M. B. Menhaj, IEEE Trans. Neural Network., 1994, 5, 989–993 CrossRef CAS PubMed.
- X. Hu, W. L. Hase and T. Pirraglia, J. Comput. Chem., 1991, 12, 1014–1024 CrossRef CAS.
- M. C. Gutzwiller, Chaos in classical and quantum mechanics, Springer, New York, 1990 Search PubMed.
- P. Shi, Y. Yang, B. Ao, P. Zhang and X. Wang, J. Phys. Chem. C, 2014, 118, 26634–26640 CrossRef CAS.

## Footnote |

† These two authors contributed equally to this work. |

This journal is © The Royal Society of Chemistry 2019 |