Orbital-Dependent Electronic Friction Significantly Affects the Description of Reactive Scattering of N2 from Ru(0001).

Electron-hole pair (ehp) excitation is thought to substantially affect the dynamics of molecules on metal surfaces, but it is not clear whether this can be better addressed by orbital-dependent friction (ODF) or the local density friction approximation (LDFA). We investigate the effect of ehp excitation on the dissociative chemisorption of N2 on and its inelastic scattering from Ru(0001), which is the benchmark system of highly activated dissociation, with these two different models. ODF is in better agreement with the best experimental estimates for the reaction probabilities than LDFA, yields results for vibrational excitation in better agreement with experiment, but slightly overestimates the translational energy loss during scattering. N2 on Ru(0001) is thus the first system for which the ODF and LDFA approaches are shown to yield substantially different results for easily accessible experimental observables, including reaction probabilities.

In the dawning age of sustainability, chemical reactions on metal surfaces play a crucial role in heterogeneously catalyzed processes that feed and fuel our modern societies. The corresponding reaction rates are usually obtained based on the Born–Oppenheimer (BO) approximation and concomitant (adiabatic) potential energy surfaces (PESs).[1,2] It has been suggested that nonadiabatic effects in the form of electron–hole pair (ehp) excitations,[3−5] which are not captured within the BO approximation, may significantly affect the underlying dynamics of molecules on metal surfaces.[6−12] However, nonadiabatic effects cannot be quantified by experimental data alone. Instead, state-of-the-art first-principles based computer simulations are mandatory, in combination with measurements from well-defined molecular beam experiments under clean ultrahigh vacuum conditions.

The current workhorse model for including nonadiabatic effects in simulations of molecular beam experiments is molecular dynamics with electronic friction (MDEF),[13,14] with two different theoretical approaches for obtaining the electronic friction coefficients: The local density friction approximation (LDFA) determines the latter based on the surface electron density[15] according to the computationally inexpensive atoms-in-jellium model.[16] The LDFA enables the inclusion of all degrees of freedom of a molecule in dynamical simulations[15,17] but at the same time implies the independent atom approximation (IAA) in most practical applications, thus neglecting any potential molecular effects.[18,19] Obviously, this is no problem for atomic projectiles. Experiments have confirmed that the LDFA yields accurate results for hydrogen atoms scattering from metal surfaces.[20,21] However, orbital-dependent friction (ODF) invokes first order time-dependent perturbation theory (TDPT) for the Kohn–Sham orbitals resulting from density functional theory calculations of an atom or molecule interacting with the surface,[13,14,22−24] so that the effects of molecule’s electronic structure (i.e., no IAA) and its interaction with the surface are taken into account. ODF is thus expected to be important for reactive scattering of molecules from metal surfaces.[18,25−27] However, the pragmatic use of broadening techniques for the calculation of ODF coefficients, which is currently without any alternative for calculating these coefficients on a large scale in the context of dynamical simulations,[25,26,28,29] has recently been criticized to affect the values obtained for these coefficients in a way that is not well-defined within the quasi-static limit of TDPT.[30] Altogether, the LDFA and ODF as formulated and implemented at present both have advantages and disadvantages. It is thus not clear a priori which model better describes dissociative chemisorption. Therefore, empirical evidence (insights based on comparisons with experimental data) is needed. Due to the very high computational cost of ODF, so far it has been used for the simulation of reactive scattering in only two systems including all six molecular degrees of freedom, i.e., H2 and D2 from Ag(111)[25,28,29] and Cu(111).[26] For these two system, no significant differences were found between reaction probabilities computed with ODF and the LDFA.

Given this situation, other systems are required that offer the possibility to distinguish LDFA and ODF, ideally by benchmarking against data from molecular beam experiments. For this, Luntz and co-workers have suggested N2 on Ru(0001) on the basis of extensive experimental and pioneering low-dimensional computational studies.[31−33] This prototypical case of highly activated diatomic molecule dissociation has received much attention due to the relevance of N2 dissociation as rate-limiting step for ammonia production via the Haber–Bosch process.[34] Recent results from LDFA calculations indicate that electronic friction is not important for the dissociative chemisorption probability,[35] whereas experiments have demonstrated that N2 molecules associatively desorbing from Ru(0001) experience a large amount of vibrational quenching,[36,37] which cannot be explained using BO-based theory.[31,38]

In this work, we show that our high-dimensional ODF model[26] applied to N2 on Ru(0001), which includes frictional couplings and the motion in all six molecular degrees of freedom, reduces the dissociative chemisorption probability by about 50% compared to both adiabatic calculations and the LDFA, an effect that is almost as large as observed going from a frozen to a moving surface. Both ODF and LDFA yield good agreement for the normal translational energy loss during scattering with the corresponding experimental data, but the ODF description of concomitant vibrational excitations is better. N2 on Ru(0001) is thus the first system for which the ODF and LDFA approaches are shown to yield substantially different results for easily accessible experimental observables including reaction probabilities. The error bars on the experimental sticking probabilities still prevent an unequivocal verification of the quantitative performance of both methods. Nevertheless, our results pave the way for subsequent improved experimental studies, which can elucidate whether indeed ODF better describes the nonadiabatic reaction in this benchmark system.

We have performed MDEF[13,14] calculations according to the generalized Langevin equation (GLE)where V(rN, rRu) is the potential energy surface that describes the (electronically adiabatic) interaction between a N2 molecule and the Ru(0001) surface consisting of mobile surface atoms described by coordinates rN and rRu, respectively. We have used the high-dimensional neural network (HD-NNP) PES from Shakouri et al.,[39] which has been fitted to a DFT reference data set based on the RPBE functional[40] using the Behler–Parinello method.[41]

The friction tensor ηN and the random forces describe the nonadiabatic coupling of the N2 molecules with electron–hole pair excitations in the surface at the surface temperature T. We have calculated the ODF tensor from density functional perturbation theory in the same way as in our previous work.[26] This 6 × 6 tensor depends on the coordinates of the two nitrogen atoms, which are most conveniently described in the coordinate system shown in Figure A. Subsequently, we have constructed an accurate continuous representation using a neural network approach as detailed in the Supporting Information, which, together with the HD-NNP PES, allows calculating a large enough number of trajectories to obtain sticking probabilities that can be compared to experimental data.[39,42] We neglect the influence of surface atom displacements on the friction tensor. Previous work has shown that this hardly affects the results of ab initio MDEF calculations based on the LDFA for N2 on Fe(110).[43]

Figure 1

(A) Six-dimensional coordinate system for the description of N2 molecules on Ru(0001), consisting of the center of mass (COM) coordinates (X,Y,Z) and the N2 bond distance d as well as the polar angle θ and azimuthal angle ϕ. X,Y,Z = 0 corresponds to the position of a Ru atom in the surface plane (top site). (B) Top view of a N2 molecule with its molecular axis parallel to the surface over a bridge site in a bridge-to-hollow orientation (). a denotes the surface lattice constant. First- (second-) layer Ru atoms are shown in (transparent) green. Dashed black lines show the periodic boundary conditions of a 2 × 2 super cell.

In order to numerically integrate eq , we have adapted a recently suggested Liouville operator technique, denoted by OVRVO in ref (44), which simplifies to the conventional velocity-Verlet algorithm[45] in the absence of friction. This technique allows defining a conserved total energy,[46] which has enabled us to monitor and thus ensure the accuracy of the numerical integration of the trajectories. It has also greatly simplified the analysis of the energy exchange with the surface, in particular for the nonadiabatic energy, which is dissipated into the electron–hole pair excitations in the surface.

In the following we compare results obtained with our ODF model to LDFA and adiabatic simulations for a mobile surface (BOMS) without electronic friction (i.e., without the last two terms in eq ), focusing on a surface temperature Ts = 575 K, which is comparable to the experimental conditions for which data has been obtained[31,32,47] and relevant for catalytic conditions of the Haber–Bosch cycle. In the calculations labeled ODF and LDFA, the motion of the surface atoms is also modeled. In more approximate calculations we neglect electron–hole pair excitations and also freeze all the ruthenium atoms at their equilibrium positions, resulting in the so-called Born–Oppenheimer static surface (BOSS) model, which does not allow for any energy exchange with the surface.[42,48]

Figure shows the initial sticking probabilities for the dissociative chemisorption of N2 on Ru(0001). Except for the lowest incidence energy (1.50 eV), the BOSS model does not reproduce the experimental results.[47] Including surface atom motion (BOMS) reproduces the experiment within the large systematic error bars (see Supporting Information for a discussion), as has already previously been shown.[35,39] The LDFA does not yield any significant changes compared to the BOMS model. The ODF does yield significant changes, and brings the theory in better agreement with the best experimental estimates of S0. We are not aware of any other system where electronic friction has an effect on S0 with the same magnitude as surface motion. The effect of the ODF is quantified on a linear scale in the inset of Figure , which shows the decrease of the reaction probability of both electronic friction models relative to the BOMS results. Also on this scale, the results from the LDFA and BOMS model are hardly distinguishable. ODF, however, decreases the sticking probability, relative to BOMS, from lower to higher incidence energies by 61% to 41%.

Figure 2

Reaction probability S0 as a function of the average incidence energy ⟨Ei⟩ calculated with the ODF model from this work for a surface temperature Ts = 575 K (purple diamonds). Corresponding results from Shakouri et al.[35] based on the LDFA (Ts = 575 K, red triangles), the BOMS (blue triangles), and the BOSS model (green squares). Experimental data from Diekhöner et al.[47] are shown for comparison (gray circles). The inset shows the ratio of reaction probabilities calculated with both electronic friction models relative to the corresponding adiabatic BOMS results.

We have also analyzed energy exchange for N2 scattering from Ru(0001) for two available experiments at different incidence energies (⟨Ei⟩). First, we compare with the translational energy loss of N2 along the surface normal (−⟨ΔE⊥⟩) in Figure A. At high ⟨Ei⟩, the BOMS model slightly overestimates this energy loss, while it underestimates at low ⟨Ei⟩. The effect of both ODF and LDFA is small for this observable and the same result is obtained as for BOMS for low ⟨Ei⟩, while ODF slightly overestimates at high incidence energies. In the second experiment, an upper bound of 0.05 eV has been obtained for the amount of vibrational excitation during N2 scattering from Ru(0001) at ⟨Ei⟩ = 2.8 eV. Earlier calculations within the BOSS model[48] using a different RPBE-based PES[42] have significantly overestimated this energy transfer ⟨ΔEvib⟩. As shown in Figure C, we reproduce this finding for BOSS-model-based simulations with our HD-NNP PES. Including surface mobility (BOMS) reduces the average vibrational excitation by up to 50% at the highest incidence energies, but the results are not yet compatible with the upper bound estimated from the experiments. LDFA does not yield any further improvement. Quite in contrast, ODF leads to a further reduction of 50–60% for all incidence energies, such that only this electronic friction model is compatible with the experimental upper bound.

Figure 3

(A) Average translational energy loss along the surface normal −⟨ΔE⊥⟩ and (B) average change of the vibrational energy ⟨ΔEvib⟩ as a function of the average incidence energy ⟨Ei⟩ for molecules scattered from the surface. Results from adiabatic calculations according to the BOSS model (green squares), the BOMS model (blue triangles), LDFA model (red triangles), and ODF model (purple diamonds) are plotted for a surface temperature Ts = 575 K. Experimental data from Mortensen et al.[32] (gray circles) are shown for comparison in (A). In (B), the maximum vibrational energy change of 0.05 eV at ⟨Ei⟩ = 2.8 eV estimated in the same study[32] is indicated (gray bar). Error bars in (A) and (B) are smaller than the line size and are consequently omitted.

We attribute the big effect of ODF on S0 and vibrational excitation to the extremely large electronic friction acting along the N2 bond axis. Figures A–C show the corresponding friction elements η, η, and η, respectively, along the minimum energy path q, at the bridge site as shown in Figure B, obtained by Shakouri et al.,[35] and also depicted in Figure D, for the HD-NNP PES used in this work. The ODF tensor elements for N2 on Ru(0001) are more than five times larger than for H2 on Cu(111).[26] Furthermore, as has been observed before,[18,25,26,33] ODF predicts increased friction along the N2 bond compared to LDFA, i.e., in contrast to ηODF ≈ 4 ηLDFA at the transition state. Moreover, the coupling between the d and Z modes is significantly smaller for LDFA compared to ODF along the reaction path that is dynamically most relevant for dissociation as discussed by Shakouri et al.[39] and in our Supporting Information. Hence, it is not surprising that dynamical processes that involve N2 bond activation, like dissociation on and vibrationally inelastic scattering from the Ru(0001) surface, are accompanied by a significantly larger concomitant nonadiabatic energy loss with ODF than with LDFA. Luntz and Persson have already pointed out large differences between ODF and LDFA,[18] but they did not explicitly compute sticking probabilities.[33] Juaristi et al. have emphasized that modeling the motion of a N2 molecule on a metal surface in all six degrees of freedom is a crucial advantage of the LDFA[15,17] over the dynamical ODF-based model in Luntz’s and Persson’s pioneering work, which only included Z and d.[33] Our results for η, η, and η slightly differ from theirs but still maintain the same essential features that distinguish LDFA from ODF. We show in the Supporting Information that these differences are related to a slightly different minimum energy path and the use of a different exchange-correlation functional.

Figure 4

(A–C) Friction tensor elements related to the center of mass distance to the surface (η), the bond length (η), and the friction-induced coupling between these two (η), respectively, along the minimum energy path q for dissociative chemisorption over the bridge site in the bridge-to-hollow orientation with the molecular axis slightly tilted off parallel from the surface (θ = 84°, see Figure for the molecular coordinate system). This path is depicted in (D) together with the corresponding two-dimensional PES cut. The purple (red) lines indicate the electronic friction obtained for ODF (LDFA). Purple (red) dots show the ODF (LDFA) results from previous work of Luntz and Persson.[18,33] The transition state for dissociation is located at the vertical gray line in (A–C) (q = 0 Å) and indicated by the empty circle in (D). Negative numbers up to transition state denote the approach from the gas-phase (i.e., decreasing Z above the surface). We note that in ref (33)q is defined for the strictly parallel approach of the N2 molecule toward the surface (θ = 90°), but this does not correspond to the minimum energy path in our HD-NNP PES.[35]

In conclusion, for N2 on Ru(0001), our ODF approach, which includes nonadiabatic coupling of the motion in all six N2 molecular degrees of freedom due to ehp excitations, yields a reduction of the dissociative chemisorption probability by about 50%. Such a large effect on a reaction probability, which is of the same order of magnitude as the inclusion of surface atom motion, has never been observed for MDEF calculations before, most of which have been based on the LDFA model. Both ODF and LDFA agree with the currently available experimental data within its sizable systematic error bars, but ODF comes closer to the best experimental estimates for S0. ODF slightly overestimates the translational energy loss during scattering, but yields results for vibrational excitation in better agreement with the aforementioned experiments than LDFA. More accurate measurements for reactive scattering of N2 from Ru(0001) would in principle allow a verdict on whether ODF or LDFA better describes the effect of ehp excitation on reaction, or whether they are equally good. We believe it is possible to measure S0 with errors no greater than a factor 1.1, using a feasible strategy that would lead to a cancellation of systematic errors (or calibration factors), as discussed in the Supporting Information.[49] Likewise, absolute vibrational excitation probabilities could be measured more accurately. This might also allow to further develop theoretical modeling of nonadiabatic dynamics at metal surfaces, for example, by including higher order perturbation terms (electron-mediated phonon–phonon coupling), which Novko et al. demonstrated to play a crucial role for the nonadiabatic vibrational damping of CO on Cu(100).[30,50,51] Although not suggested by our present results, improvements of the exchange-correlation functional defining the PES might be required in order to achieve quantitative agreement[52,53] with the more accurate experimental data to be measured. Given the importance of N2 on Ru(0001) as a prototypical case of highly activated dissociative chemisorption, such future steps would be tremendously important to quantify how much nonadiabatic effects need to be accounted for in heterogeneous catalysis.

Computational Details

We have obtained orbital-dependent electronic friction tensors from density functional perturbation theory (DFPT)[54] results based on a computational setup similar to the one we used before for H2 on Cu(111).[26] Briefly, we have performed DFPT calculations as implemented in the Quantum Espresso package[55] for a 2 × 2 Ru(0001) slab with five layers employing the RPBE functional[40] as implemented in LibXC.[56] Using ONCV[57] pseudopotentials from the SG15[58] library together with a plane-wave cutoff of 816 eV, a 18 × 18 × 1 k-point grid, and a Gaussian envelope technique with a width of 0.6 eV for the sum over electronic states[23,30] yields converged results for the friction tensor elements. Tests with smaller widths (down to 0.15 eV) are provided in the Supporting Information. A continuous representation of the 6 × 6 frozen-surface friction tensor was obtained using neural networks constructed with the help of the TensorFlow package.[59] As detailed in the Supporting Information, we base this representation on a molecular orientation that corresponds to the most significant reaction path for dissociative chemisorption. Improving our previous approach,[26] special care has been taken in order to ensure positive definiteness of the friction tensor and to keep the amount of neural network weight parameters as small as possible (three hidden layers with 20 nodes each) by fitting all 21 independent friction tensor elements simultaneously.

Quasi-classical trajectory calculations with a time step of 0.3 fs were performed using the LAMMPS package,[60,61] into which we have implemented our adaptation of the OVRVO algorithm.[44] At every time step, we apply OVRVO by rotating to the six-dimensional coordinate system in which the ODF tensor is diagonal. The LDFA friction tensor, which is diagonal in the Cartesian representation, is the same as in ref (35), and our OVRVO implementation perfectly reproduces the results from that work. Likewise, we have taken over the equilibration procedure of the surface slab for generating initial conditions at Ts > 0 K.