Analysis of Energy Dissipation Channels in a Benchmark System of Activated Dissociation: N2 on Ru(0001).

The excitation of electron-hole pairs in reactive scattering of molecules at metal surfaces often affects the physical and dynamical observables of interest, including the reaction probability. Here, we study the influence of electron-hole pair excitation on the dissociative chemisorption of N2 on Ru(0001) using the local density friction approximation method. The effect of surface atom motion has also been taken into account by a high-dimensional neural network potential. Our nonadiabatic molecular dynamics simulations with electronic friction show that the reaction of N2 is more strongly affected by the energy transfer to surface phonons than by the energy loss to electron-hole pairs. The discrepancy between the computed reaction probabilities and experimental results is within the experimental error both with and without friction; however, the incorporation of electron-hole pairs yields somewhat better agreement with experiments, especially at high collision energies. We also calculate the vibrational efficacy for the N2 + Ru(0001) reaction and demonstrate that the N2 reaction is more enhanced by exciting the molecular vibrations than by adding an equivalent amount of energy into translation.

Introduction

Accurately describing chemical reactions often requires the knowledge of phenomena controlling the energy exchange between different degrees of freedom. Among these phenomena, electronic excitations are of particular importance to reactions involving metal surfaces, where the lack of an energy band gap paves the way for creation of electron–hole pairs (ehps). In this respect, experimental evidence has been obtained, confirming the excitation of charge carriers during interactions of a diverse range of atoms and small molecules with metal surfaces.[1−4] Chemicurrent measurements, for instance, revealed the presence of electronic excitations during adsorption of a variety of molecular and atomic species on polycrystalline silver.[1,2] There has also been another set of experiments indicating that the electron transfer from metal surfaces to molecules can be promoted by collisions of highly vibrationally excited molecules with a metal surface.[3,4] This suggests that the ability to accurately describe chemical processes at metal surfaces would necessitate describing the effect of electronic excitation.[5−8] This is especially true for some systems: as shown for the inelastic scattering of H from Au(111),[8] only dynamical simulations taking ehps explicitly into account provided scattering probabilities in good agreement with experimental results.

Apart from ehp excitation giving rise to nonadiabatic effects, energy transfer to surface phonons (i.e., the thermal motion of surface atoms) can also affect the reaction dynamics strongly.[9−13] The role of surface phonons is highlighted in reactions where the collision of a heavy projectile with the surface allows a considerable amount of energy to be transferred to the surface. The reactivity of the projectile is therefore determined by the excess of kinetic energy remaining in the molecular motion to cross the reaction barrier.[14] Another notable aspect of the thermal motion of surface atoms is its direct influence on the reaction barrier height, modifying it dynamically during the reaction process.[15−17] This clearly shows that for a better understanding of surface reactions, it is often mandatory to go beyond static surface models, which have been the dominant theoretical approach for decades.[18−24]

Describing the energy exchange with the surface during N2 reaction on Ru(0001), which is a prototypical case of highly activated diatomic molecule dissociation, here we address the combined effect of surface phonons and ehp excitation on the reaction probability S0 of N2. The N2 + Ru(0001) system has attracted a lot of attention both experimentally[14,25−32] and theoretically.[13,23,24,31] Earlier theoretical approaches indicated that adiabatic simulations of N2 + Ru(0001) compared to experiment overestimate the value of S0 by 2 orders of magnitude, while adding the coupling to ehp excitation results in a semiquantitative agreement with experiments.[31] However, these calculations were strictly limited to low-dimensional models where only the vibrational and normal translational motion of N2 was included and the coupling to surface phonons was approximated by a single harmonic oscillator.[33] Subsequent studies[23,24] combining an evaluation of the static properties of the potential energy surface (PES) with quasi-classical dynamics simulations performed within the Born–Oppenheimer static surface (BOSS) model suggested that the discrepancy between the experiments and earlier adiabatic calculations[31] can largely be resolved by accounting for all the molecular degrees of freedom. Specifically, a study of the transition state of N2 + Ru(0001) suggested that the reaction barrier represents a very narrow bottleneck to the reaction. However, the role of surface temperature Ts and also the dependence of S0 on electronically nonadiabatic effects remained unclear, and it has not yet been demonstrated with dynamics calculations that the reaction indeed only proceeds for a narrow range of impact sites and incidence angles. Very recently, it was demonstrated by large-scale molecular dynamics (MD) simulations that the interaction of N2 with Ru(0001) and the motion of Ru atoms can be described accurately by a high-dimensional neural network potential (HD-NNP), which was constructed based on a set of density functional theory (DFT) data.[13] The computed S0 showed reasonable agreement with experimental results to within the experimental error, and it was suggested that the remaining discrepancies might be due to creation of ehps.

Using the same HD-NNP, here, we address the two remaining questions discussed above. Most importantly, here, we go beyond the previous adiabatic model and study the effect of ehp excitations by performing MD simulations with electronic friction (MDEF). As already noted, earlier, it had been speculated that modeling ehp excitation is needed to bring theory in agreement with the experiment for N2 + Ru(0001) (ref (31)). Here, we employ the local density friction approximation (LDFA) to obtain electronic friction coefficients for MDEF and rely on the independent atom approximation (IAA) as originally suggested.[34] Unlike recent extensions of the LDFA,[35,36] LDFA-IAA has the advantage that the electronic friction coefficients are completely defined based on the density of the bare surface alone. Furthermore, LDFA-IAA has been used successfully in studies of atoms scattering from metals[7,8,68] and even of reactive scattering of molecules from metal surfaces,[34,37−41,69] although a disadvantage of the LDFA-IAA method is that it does not take into account molecular electronic structure effects on the friction coefficients. Methods using orbital dependent friction (ODF)[6,41] do take these effects into account, but using these methods is computationally expensive, with a challenge being to construct a continuous representation of the coordinate-dependent friction coefficients, as would also be required for a PES when sticking coefficients are to be calculated.[6] Furthermore, the ODF evaluation of friction coefficients may lead to spurious electronic memory effects which may affect the MDEF.[42] Finally, a recent study suggests that MDEF based on LDFA and ODF friction coefficients may lead to very similar computed sticking coefficients for late barrier systems, as found for H2 + Cu(111).[6] Our LDFA-MDEF simulations suggest that nonadiabatic effects only play a small role in the dissociative chemisorption of N2. This is in agreement with recent studies on low barrier N2–metal systems like N2 + W(100) and N2 + W(110),[37] but as pointed out in ref (43), studies on such systems are not necessarily predictive of nonadiabatic effects in systems with high late barriers, such as N2 + Ru(0001). Contrary to the seemingly small effect of ehps, surface phonons have a large effect on the reaction probability, as the energy transferred from N2 to the surface is on average between ≈0.4 and 1.1 eV, for the collision energy range 1.5–3.25 eV, which is roughly an order of magnitude larger than the energy lost to ehps.

The second question we address here is whether the MDEF simulations indeed explicitly show that the reactivity is further reduced through effects not concerning the dissipative degrees of freedom, but the molecular degrees of freedom. We show through explicit dynamics calculations that the effect of the transition state representing a narrow bottleneck to reaction[23,24] indeed leads to reaction only occurring for impact sites near the transition state, and for orientations of N2 close to that occurring in the transition state. Finally, we evaluate the state-selected reaction probability of N2 and show that adding vibrational energy to N2 promotes its reaction more efficiently than increasing the collision energy.

This paper is organized as follows. In Section , we give details of our theoretical model used in the MD and MDEF simulations. The constructed HD-NNP is briefly discussed, and we also explain how the Langevin equation of motion is integrated. Results of our MD and MDEF calculations are detailed in Section for both static and mobile surfaces of Ru(0001). Comparisons are made between molecular beam experiments and the dynamical results obtained with and without electronic friction. Energy exchange with surface phonons and energy loss to ehp excitation are all described in this section. We also discuss based on a statistical analysis the key dissociation geometries of N2 for which the molecule reacts with the Ru surface. We summarize our main results in Section .

Model and Method

We study the dissociative chemisorption of N2 on Ru(0001) in the presence of both ehp excitation and surface atom motion. Our MDEF simulations have been carried out within the LDFA model using a preconstructed HD-NNP of N2 + Ru(0001),[13] which accounts simultaneously for both molecular and surface atom motion. The molecular coordinate system used in our calculations and the high-symmetry points of the Ru(0001) surface are illustrated in Figure . The training data set used as reference for the neural network fit was collected from a set of DFT calculations performed at the generalized gradient approximation. All DFT calculations were done with the RPBE density functional[44] using the Vienna Ab initio Simulation Package code.[45−47] To model the N2 + Ru(0001) system, a 3 × 3 surface unit cell comprising 7 Ru layers was employed as our convergence tests for the barrier height have shown that the latter is converged to within 30 meV for the chosen surface slab model.[13] The surface Brillouin zone of Ru(0001) was sampled by a Γ-centered 7 × 7 × 1 k-point grid in combination with a broadening according to Methfessel and Paxton[48] using a width of 0.3 eV. The bottom Ru layer has been kept frozen in the training of the neural network fit (ref (13)), and in the MD and MDEF simulations reported below.

Figure 1

Molecular coordinates of N2 in the center of mass reference frame and high-symmetry points of the Ru(0001) surface. The light and dark surface atoms correspond to the topmost and second topmost surface layers, respectively.

Equations of Motion in the Presence of ehps

Neglecting nonadiabatic effects for the Ru atoms in the following, their dynamics is simply given by Newton’s equations. The effect of ehps on surface reactions can be incorporated into classical dynamics by using the concept of electronic friction.[49−51] In this case, as shown in ref (50), the equations of motion of the adsorbate atom i with mass m is given by generalized Langevin equationHere, V denotes the PES as given by the HD-NNP described in the previous section depending on the coordinates r1, r2, ... of the adsorbate and surface atoms. In general, η is a tensor with each element having the same atomic coordinate dependence as V. As described in more detail in the following section, in this work, η is simplified to a scalar η for each adsorbate atom that only depends on the coordinates of this particular atom. According to the second term on the right hand side of eq , η is responsible for the dissipative effect of electronic friction, whereas Rel corresponds to the randomly fluctuating force that originates from nonadiabatic scattering due to thermal surface electrons and is typically approximated by a Gaussian white noise.[7,38] Considering the friction-dependent terms in eq , the energy exchange Ω with the electronic degrees of freedom of the surface after a time t can be described according to

Friction Coefficients According to the LDFA Method

There is yet some controversy which theoretical model is most suitable to describe the electronic friction tensor η.[6,41,52,53,67] One of the computationally most efficient models proposed so far is the LDFA method in combination with the IAA.[34] In this method, the friction tensor is diagonal and isotropic and thus, a scalar friction coefficient for each adsorbate atom is independently obtained according to the atoms-in-jellium model.[54] The density of the corresponding homogeneous electron gas (i.e., the jellium) is taken to be the density of the ideal frozen surface at the same point where the atom is placed. Ways to account for surface atom displacements in the LDFA method have been suggested,[35] but their effect has been found to be negligible for the reaction probability of N2 on Fe(110),[36] which is why the frozen surface description has been adapted here.

In the preceding work, the friction coefficient for N has been calculated as a function of the electronic density using the LDFA model and has been fitted to the following analytical function[55]where rs is the Wigner–Seitz radius and is inversely related to the embedding electron density ρe, that is, rs = (4πρe/3)−1/3. The fitting parameters of eq , in atomic units, are as follows:[55] (A1,A2) = (39.298,–34.62)ℏ/r, (B1,B2) = (−0.127,0.333), and (C1,C2) = (0.838,0.999) rb–1, where rb is the Bohr radius. The friction coefficient η of a single N atom on Ru(0001) is shown in atomic units in Figure . The surface electron density ρe of the ideal frozen Ru(0001) (i.e., without including N2 in the simulation cell) has been calculated using DFT with the same computational setup as described at the beginning of Section . In doing this single DFT calculation, a Zn core pseudo-potential (Ru pv) has been used for Ru, allowing 14 electrons of Ru in its 4p65s14d7 configuration to be modeled. We have verified that this pseudo-potential yields converged values of ρe and η at dynamically relevant distances from the surface. As is evident from Figure b, by moving N away from the surface along the Z direction, the friction coefficient decays roughly exponentially. We also find from Figure c that the friction coefficient is much larger on the top site (i.e., X = Y = 0) than, for instance, the bridge site (X ≈ 1.37 Å, Y = 0) where the friction coefficient arrives at its minimum. As will be discussed subsequently, this could largely diminish the dependence of S0 on the ehp excitation because the reaction of N2 takes place mostly at the vicinity of the bridge site.

Figure 2

(a) Contour plot of the surface density-dependent friction coefficient η for a single N atom on Ru(0001), which is obtained from the LDFA model for Y = 0. (b,c) Variations of the friction coefficient, in atomic units, along the orthogonal directions Z and X which are perpendicular and parallel to the surface, respectively. Note that Y = 0.

Sampling Initial Conditions

In sampling the initial position of the surface atoms at the experimental Ts, we carried out NVT classical MD simulations (i.e., the number of atoms, unit cell, and the surface temperature are kept fixed) using the large-scale atomic/molecular massively parallel simulator (LAMMPS).[56,57] To ensure that the thermally distorted surface atoms mimic the thermodynamic properties of a real canonical ensemble, the Nosé–Hoover thermostat[58,59] with a damping parameter of 0.05 ps has been employed. In our NVT simulations, a time step of 0.5 fs has been used for the velocity-Verlet[60,61] integration. The surface equilibration process continued for 3 × 105 time steps, that is, more than the number of trajectories required to obtain S0 with good statistics. With the aid of the long lasting NVT simulations performed, distinct random snapshots have been selected for the initial input conditions of the surface atoms in the MD and MDEF simulations.

The initial states of N2 in the gas phase are randomly sampled according to molecular beam experiments performed in ref (14). The energy spread of the measured collision energy was ΔE/E̅ = 0.15, where ΔE denotes the full width at half maximum of the energy distribution and E̅ is its average collision energy. To sample the initial translational velocities of N2 from the aforementioned energy distribution, we employed the following flux-weighted velocity distributionwhere C is a constant, v is the translational velocity of N2, vs is the so-called stream velocity, and α determines the width of f(v). Considering that C is only a normalization factor, the average collision energy can be evaluated analytically from eq , which becomesprovided that α ≪ vs. The unknown parameters vs and α in eq have to be determined such that the experimental condition ΔE/E̅ = 0.15 is satisfied. After that, the translation velocity v can be sampled randomly using the distribution function f(v) given in eq .

To determine the initial rovibrational states of N2, a nozzle temperature of Tn = 1100 K has been used.[14] For each set of the ν and j rovibrational quanta, the Tn-dependent Boltzmann weight is then given by

In eq , Nf is the normalization coefficient, g = 2j + 1 the rotational degeneracy factor, and w(j) is the weight due to the ortho–para ratio of N2. Given that the nucleus spin of the 14N isotope is s = 1, the total spin of N2 with even integers takes place twice as often as the odd ones, meaning that w(j) is equal to 2/3 for even j and 1/3 for odd j. Also, similar to ref (13), the vibrational and rotational temperatures in eq are taken equal to Tvib = Tn and Trot = 0.1Tn, respectively. Equation is necessary to calculate the mono-energetic reaction probability of N2, which is directly expressed as the sum of Boltzmann-weighted contributions from all state-selected reaction probabilities Sν, that is

To compare with measured reaction probabilities in molecular beam experiments (ref (14)), the reaction probability S(E,Tn) also has to be averaged over the collision energy distribution

QCT Calculations of N2 on Ru(0001)

Dissociation probabilities for the reaction of N2 with both mobile and rigid surfaces of Ru(0001) are obtained by performing a large number of MD and MDEF simulations, where the latter go beyond the Born–Oppenheimer approximation as described before. Our dynamics calculations were done at the quasi-classical limit by taking into account the zero point energy and the additional initial rovibrational energy of N2 in the gas phase. The total energy of N2 + Ru(0001) and interatomic forces during the time evolution of the dynamics are evaluated on the fly by the HD-NNP. The use of an HD-NNP instead of ab initio MDs (AIMD) is computationally inevitable because the low reaction probability of N2 (S0 ≈ 10–5 to 10–3)[14] necessitates running many trajectories. The number of trajectories used in our MD and MDEF simulations ranges from 5 × 104 up to 2.5 × 105 which corresponds, respectively, to the highest (E̅ = 3.25 eV) and lowest (E̅ = 1.5 eV) collision energies considered. With such a large number of trajectories determined, it is possible to describe S0 with a reasonably good statistics.

For the numerical integration of the classical equations of motion, either with or without friction, a MD time step of Δt = 0.3 fs has been used. The energy error created because of this value of Δt is on average less than 1 meV for both the energy drift and energy oscillations in adiabatic simulations. To solve the Langevin equation of motion, that is, eq , the stochastic algorithm of Ermak and Buckholz[62] has been employed.[39,40] In this method, the dynamic variables of the mobile atom are propagated in time using the following recurrence relations[63]where F⃗(t) = −∇V is the adiabatic force acting on the adsorbate atom at position r⃗(t), and R⃗g is a vector of three standard Gaussian random numbers (normal distributions centered around 0 with a standard deviation of 1), the components of which are independent of one another. The Gaussian random numbers are generated in our simulations by Marsaglia’s algorithm as implemented in the LAMMPS code.[64] The additional parameters in eqs and 11 are defined as follows[63]

Note that in the zero friction limit (η → 0), one obtains a0 = a1 = 1 and a2 = 0.5; therefore, eqs and 11 directly yield the well-known velocity-Verlet algorithm[60,61] in adiabatic calculations. The position and velocity of the Ru atoms were updated according to the latter, whereas the N atoms were propagated using an in-house implementation of the algorithm presented.

All trajectories started with the N2 molecule approaching the metal surface at a distance of Z = 5.5 Å [at this distance, the interaction potential between N2 and Ru(0001) is less than 1 meV]. Also, each trajectory calculation ran until one of the following criteria was met: (i) the N–N interatomic distance is larger than 2.7 Å, in which case the N2 is assumed to have reacted with the surface and (ii) the molecule moves away from the surface at a distance of Z > 5 Å, in which case the molecule is assumed to be scattered by the surface.

Results and Discussion

Reaction Probability

We first examine the reaction probability S0 and compare our MD and MDEF results with one another and with molecular beam measurements (ref (14)). The theoretical and experimental data of S0 are all summarized in Figure . In comparison with adiabatic simulations, the dissociation probability obtained within the LDFA method (dashed lines) show slightly better agreement with the experimental S0, especially for 2.5 < E̅ eV. This holds true for both rigid (Figure a) and mobile (Figure b,c) surfaces, in the sense that the reaction probabilities computed with the LDFA method were somewhat smaller than the computed adiabatic reaction probabilities, whereas both were larger than the experimental values. However, the computed reaction probabilities for the experimental surface temperature, that is, Ts = 575 K, reproduce the experimental results within experimental error both with and without friction (see Figure b). To better understand why the effect of ehp excitation is mostly visible at high E̅ we evaluated the energy dissipated to ehps, that is, Ω, using eq (for sufficiently small MD time steps, the difference between the total energies at the last and first time steps can also accurately yield Ω). The resulting Ω for reacted trajectories is roughly two times larger for E̅ ≥ 3 eV than for low incidence energies E̅ < 2 eV. Another notable aspect that is directly seen by comparing Figure b,c is that increasing the surface temperature Ts from 0 to 575 K has a much larger effect on the reaction than the relatively minor contribution of ehps to the dissociation probability S0.

Figure 3

(a–c) Reaction probability S0 as a function of the average collision energy. Our QCT calculations have been performed with and without friction for both rigid (a) and mobile surfaces at Ts = 575 K (b) and Ts = 0 (c). The randomly fluctuating force in (a,b) is sampled according to the experimental surface temperature (i.e., Ts = 575 K), whereas in (c), the random force is zero [that is, Rel(Ts = 0) = 0]. The experimental S0 is related to ref (14).

Energy Exchange with Surface Phonons and ehps

Keeping track of the kinetic (EkRu) and potential (EpRu) energies of the surface at the beginning and the end of trajectories in which the molecule is scattered by the surface, the average energy exchanged between the molecule and surface phonons has been obtained from , where i and f represent the initial and final surface structures in each trajectory, respectively, and N is the number of trajectories (see Figure a). For a meaningful analysis in which the direction of energy flow is also taken into account, the energies transferred from the molecule to the surface and vice versa are distinguished by δE+ and δE–, respectively. On the basis of this notation, δE+ only takes the average over trajectories in which (EkRu + EpRu)f > (EkRu + EpRu)i and δE– does the same for (EkRu + EpRu)i > (EkRu + EpRu)f. As a result, the average energy exchanged with the surface atoms can be calculated directly from δE = δE+ – δE–. By comparison between δE+ and δE– in Figure a, we find that the transfer of energy from the surface to N2 takes place very rarely, particularly at high collision energies (for E̅ ≥ 2.5 eV, the mean energy δE– makes almost no contribution so that δE ≈ δE+). In addition, as seen in Figure a, a significant portion of the collision energy of N2 (roughly 1/3) is transferred to the surface atoms, indicating that the surface phonons play a key role in decreasing the reactivity of N2.

Figure 4

(a) Average energy exchanged between the surface atoms and molecules scattered from the surface. The red and blue bars correspond to unidirectional energy flows from N2 to the surface and vice versa, respectively. The MDEF simulations have been performed for Ts = 575 K. (b) Comparison of the computed energy transfer to surface phonons and ehps with experiments (ref (32)).

The energy losses δE and δE + Ω we computed in the calculations for Ts = 575 K using MDEF simulations is compared to the energy loss accompanying rotationally elastic scattering of j = 8 N2 as extracted from experiments at Ts = 610 K (ref (32)); see Figure b. In view of the approximations that needed to be made to extract energy loss to the surface from the experiments in ref (32), the agreement between theory and experiment is remarkably good.

As for the competition between the energy transfer to the surface phonons and ehps and to identify which one is the dominant mechanism affecting the N2 dynamics, the mean values of the energy exchange with ehps, ⟨|Ω|⟩, and δE are evaluated for a large number of trajectories for which the impinging N2 has been scattered by the surface; see Tables and 2, respectively. Comparison of the data obtained reveals that the average energy exchanged with surface phonons completely dominates the energy lost to ehps (on average, δE is an order of magnitude larger than Ω). This allows us to conclude that, within the LDFA and the IAA, the surface phonons influence the reactivity of N2 to a much larger extent than ehp excitation. This result is in agreement with results obtained earlier for N2 + Fe(110) using AIMD simulations.[11] Also, by direct comparison of the data in Table , we notice that the amount of energy transferred to the surface phonons is about ∼25–40% higher at 0 K than at Ts = 575 K. This can rationalize that the dissociation of N2 on Ru(0001) is promoted by increasing the surface temperature, assuming that reduced energy transfer to the surface better facilitates to overcome barriers to reaction.

Table 1

Average Energy Loss, ⟨|Ω|⟩, Due to ehp Excitation, in eV, for Molecules Scattered from the Surface

collision energy (eV)	1.50	1.75	2.00	2.25	2.50	2.75	3.00	3.25
rigid surface	0.039	0.044	0.049	0.054	0.059	0.064	0.069	0.074
mobile surface	0.041	0.046	0.052	0.057	0.062	0.068	0.074	0.080

Table 2

Average Energy Transferred from N2 to the Surface (δE+), in eV, for Ts = 575 and 0 K

collision energy (eV)	1.50	1.75	2.00	2.25	2.50	2.75	3.00	3.25
δE at Ts = 575 K	0.394	0.500	0.606	0.712	0.818	0.921	1.022	1.113
δE at Ts = 0 K	0.542	0.653	0.773	0.899	1.030	1.162	1.298	1.434

Statistical Analysis of the Reactive N2 Geometries

Our dynamical results show that all parts of the Ru(0001) surface containing face-centered cubic (fcc) sites (see Figure ) are more reactive than those containing hexagonal close-packed (hcp) sites. To illustrate this point, the distribution of impact sites on the surface is depicted in Figure a for molecules reacting with an ideal frozen surface. The collision energy of N2 is taken as E̅ = 3.25 eV corresponding to the highest reaction probability observed, in which case the dynamics results can be described with the best statistics possible. Note, however, that the same results as shown in Figure hold true for lower collision energies. The question of why N2 reacts mostly in the fcc triangle rather than in the hcp triangle of the surface unit cell can be answered by considering the transition state of the system (see the inset in Figure b and also Table 1 of ref (13)). In the transition state, the center of mass of N2 is close to a bridge site, but shifted somewhat to the fcc site, away from the hcp site. Another notable feature that can be inferred directly from Figure a is that the reaction of N2 predominantly takes place at the vicinity of the bridge site. To further verify this, we plot in Figure b the distribution of the distance dCOM from impact sites to their nearest bridge site for trajectories in which the molecule has reacted on the surface. As seen, the curves show peak at small distances of dCOM, indicating that the dissociation probability of N2 increases for the center of mass of N2 located close to the bridge site, which is consistent with the transition state geometry of N2 + Ru(0001) depicted. Similarly, this result holds true for a mobile surface (see Figure S1a,b).

Figure 5

(a) Distribution of initial molecular center of mass impact sites for molecules reacting with an ideal frozen Ru(0001). The inset shows the distribution of impact sites in two parts of the surface containing the fcc and hcp sites. The average collision energy is taken as E̅ = 3.25 eV at which the best statistics for the reaction of N2 is achieved. (b) Distribution of the distance from impact sites to the nearest bridge site for the parts containing the fcc and hcp sites (dCOM is the lateral distance of the center of mass of N2 to the nearest bridge site). The inset displays the transition state geometry of N2 on Ru(0001). Note that the center of mass of N2 is slightly displaced from the bridge site toward the fcc site. (c) Distribution of the lateral displacement of the center of mass of N2 (δCOM) from its initial impact site to the place of dissociation.

Note that it is appropriate to discuss the reaction site in terms of the initial impact site: as the distribution of the lateral distance between the impact site of the molecule at t = 0 and the impact site at the time of reaction shows (Figure c), the molecule hardly moves along the surface while reacting. Here, the time of reaction is defined as the time at which the N–N distance becomes 2.7 Å. The noted behavior is not unusual for a late barrier reaction, as shown also for CHD3 + Ni(111) (see Figure S13 of ref (65)). This result suggests that a sudden approximation should be applicable to the molecular motion along the surfaces in quantum dynamical studies of N2 dissociation on Ru(0001).

Figure a,b shows the angular distribution of the molecules reacting with the surface as a function of the azimuthal (Φ) and polar (θ) angles, respectively, for E̅ = 3.25 eV. The N2 impact sites corresponding to the angular orientations depicted alongside the distribution’s peaks are illustrated in the inset of Figure a using the same color code. By visualizing the depicted molecular orientations on top of their impact sites, one can conclude that the distribution peaks in Figure a are most relevant to the bridge-to-hollow (B2H) dissociation geometry. In fact, the molecules mainly react at orientations for which they can cross lower barriers and for N2 + Ru(0001), the minimum reaction barrier (transition state) occurs for the B2H dissociation geometry (see the inset in Figure b). Moreover, the distribution shown in Figure b indicates that the reaction probability increases for the molecules aligned approximately parallel to the surface.

Figure 6

(a,b) Angular distribution of the N2 molecules reacting with the surface for E̅ = 3.25 eV. The red solid lines are to guide the eye only. For each local peak in (a), the corresponding angular orientation of N2 is also illustrated. In addition, the inset in (a) shows the impact sites of N2 on the surface using the same color codes as the N2 configurations depicted. Note that all results are obtained within the LDFA model.

In general, there are two main mechanisms active in decreasing the reaction probability of N2: (i) energy transfer to the surface atoms in the entrance channel (this effect is discussed in detail above) and (ii) the transition state barrier exhibits a very narrow bottleneck to reaction. As for the latter, 2D cuts through the PES are plotted in Figure a,b for two representative dissociation geometries of N2, both of which have a large contribution to the reaction probability S0. The B2H dissociation geometry considered in Figure a yields the maxima of the angular distribution in Figure a and has been selected corresponding to the molecular orientation at the transition state. As is obvious from both elbow plots for the B2H and B90 (see Figure b) dissociation geometries, the PES of N2 + Ru(0001) exhibits a late and also a very high reaction barrier (the height of the minimum reaction barrier is 1.84 eV).[13] At this late barrier, the potential increases rapidly if the molecule is rotated in θ or moved along X or Y, leading to the small dynamical bottleneck.[23] We note, however, that the extremely narrow bottleneck of the PES near the transition state has no significant consequences for N2 energy loss to ehps, as the LDFA friction coefficients are small in any case along the reaction path (η < 0.42ℏ/rb2); see Figure c.

Figure 7

(a,b) 2D cuts through the PES of N2 on an ideal frozen Ru(0001) for Φ = 30° and θ = 84° (a), and Φ = 90° and θ = 90° (b), which correspond to the transition state (B2H) and bridge-90° (B90) dissociation geometries, respectively. The difference between successive energy lines is 0.25 eV. The reaction path for each elbow plot is characterized by a dashed line. (c) Friction coefficient of a single N atom along the reaction path coordinate q. The values q < 0 and q > 0 correspond to the entrance and exit channels, respectively (q = 0 is the location of the reaction barrier). For the B90 dissociation geometry, the difference between the value of η for each N is less than 1.5 × 10–3ℏ/rb2.

That the LDFA friction coefficient is small for N2 + Ru(0001) can also be seen from Figure . The diagonal components of the friction tensor shown in this figure, that is, η and η, which act on the motion in z and r, were previously obtained with the so-called ODF model[6] using DFT calculations (refs (52) and (66)), although one should bear in mind that both these coefficients are restricted to a 2D model of N2 on Ru(0001). In the neighborhood of the reaction barrier at q = 0, the friction coefficient calculated by the LDFA model is smaller than η and η by a factor of ∼4 and ∼2, respectively (the LDFA η in Figure indicates the sum of the friction coefficients for both N atoms). This simple comparison clearly confirms that the ODF model suggests a stronger effect of ehp excitation on the N2 + Ru(0001) reaction than the LDFA method does, as noted already by Luntz et al.[52] Accordingly, use of ODF coefficients in MDEF simulations is expected to lead to a larger reduction of the reaction probability than obtained here with the LDFA, and therefore to yield results in even better agreement with the experiment.

Figure 8

Friction coefficients along the minimum reaction path for dissociation of N2 on Ru(0001). |q| is the distance from the reaction barrier. The coefficients η, η, and LDFA* are reported in refs (52) and (66). The LDFA η (red curve) is the sum of friction coefficients for both N atoms.

Vibrational Efficacy of N2

To assess the importance of preexciting the molecular vibration for promoting the N2 reaction on Ru(0001), we examine the vibrational efficacy of N2 using the following conventional definition

To calculate E(Sν), it is often convenient to utilize the graph of the state-selected reaction probability Sν in terms of the collision energy E̅; see Figure for the Sν obtained with MDEF simulations at Ts = 575 K. Making use of this figure, the vibrational efficacy ζν is evaluated from eq for Sν = 0.02 and is then summarized in Table for different pairs of ν1 and ν2. The value ζ0,1 = 1.64 presented for the excitation of the molecular vibration from ν = 0 to ν = 1 is consistent with the value of 1.6 obtained in ref (23) from adiabatic calculations, although in ref (23) the MD simulations were limited to the BOSS approximation. Furthermore, for every set of ν1 and ν2 in Table , we observe that ζν > 1, suggesting that the reaction of N2 is promoted more efficiently by exciting initial molecular vibration than by increasing its translational energy. A similar conclusion has also been drawn by the experiments of ref (29), emphasizing the significant role of initial vibrational excitation in enhancing the dissociation probability of N2.

Figure 9

State-resolved reaction probability of N2 on a mobile surface of Ru(0001) for j = 0 and the lowest four vibrational quantum numbers. The results shown are all obtained with MDEF simulations at Ts = 575 K.

Table 3

Vibrational Efficacy of N2 for Different Pairs of ν1 and ν2a

ν1,ν2	0,1	0,2	0,3	1,2	1,3	2,3
ζν1,ν2	1.64 ± 0.15	1.50 ± 0.1	1.56 ± 0.1	1.56 ± 0.2	1.66 ± 0.1	1.60 ± 0.2

The state-resolved reaction probability is taken as Sν = 0.02.

Contributions of the vibrationally ground and excited states of N2 to the reaction probability S0 are given in Table for two collision energies, E̅ = 3.25 and 2.5 eV. The data provided in this table are all derived from our MDEF simulations for Tn = 1100 K and Ts = 575 K. As seen, for the experimental nozzle temperature Tn = 1100 K, the reaction probability is predominantly determined by the ground vibrational state ν = 0, and the additional contributions, especially from highly excited vibrational states ν ≥ 2, are nearly negligible. While ν = 1 N2 makes up 9 and 18% of the reacting molecules at E̅ = 3.25 and 2.5 eV, respectively, the population of the ν = 1 state at the nozzle temperature used (1100 K) is only 4.4%. This implies that to promote the dissociative chemisorption of N2 on Ru(0001) in molecular beam experiments, the nozzle temperature has to be increased to values much higher than Tn = 1100 K.

Table 4

Contributions of Different Initial Vibrational States of N2 in the Reaction Probability S0 at E̅ = 3.25 and 2.5 eV for Ts = 575 K

vibrational state ν	0	1	2	≥3
S0(ν) at E̅ = 3.25 eV	0.0197 ± 0.001	0.002 ± 4 × 10–4	(2 ± 1) × 10–4	<10–4
S0(ν) at E̅ = 2.50 eV	0.0034 ± 5 × 10–4	(8 ± 2) × 10–4	<10–4	<10–4

Unfortunately, we cannot yet comment on the controversial issue[27,31] of whether there should be a vibrational population inversion in associative desorption of N2 from Ru(0001). This would require running 1–2 orders of magnitude more trajectories for the present lowest energies and might require even more trajectories for lower incidence energies, as low energies are more heavily weighted in associative desorption experiments from moderately hot (∼1000 K) surfaces. Furthermore, it has been suggested[31] that modeling the experiment of ref (27) should require modeling reaction at steps on imperfect Ru(0001) surfaces. Finally, modeling the experiment of ref (31) would require modeling the reaction of N2 on a Ru(0001) surface with a precoverage of 0.6 monolayer of N-atoms.

Summary

In summary, we have carried out quasi-classical dynamics simulations for the dissociative chemisorption of N2 on Ru(0001). Nonadiabatic interactions of the N atoms with ehps created in the metal surface have been described by the LDFA model within the IAA. Our dynamical calculations have been performed using an HD-NNP of N2 + Ru(0001) which accounts for both the molecular and surface atom motions. Comparisons made between the energy exchanged with surface phonons and the energy lost to ehps showed that the reactivity of N2 is much more strongly influenced by energy transfer to the surface phonons than by ehps. In fact, our MD simulations with electronic friction (from the LDFA) turned on suggest much smaller nonadiabatic effects due to ehps than previously expected for N2 + Ru(0001).[31] Both reaction probabilities obtained with MD and MDEF simulations show good agreement with the experimental S0 (to within experimental error), although the MDEF results are in somewhat better agreement with the experiment at high collision energies. A small effect of LDFA electronic friction on dissociation of N2 on Ru(0001) is in line with earlier findings for N2 dissociation on transition metals (ref (53)). Future calculations will have to show whether a larger effect will be seen if friction coefficients computed at the ODF level of theory are used for N2 + Ru(0001), as Figure and ref (52) might be taken to suggest. Nevertheless, it is shown that the present MDEF simulations yield good agreement with the experiment for energy transfer to the surface. A statistical analysis has been provided for N2 dissociation geometries in trajectories in which the molecule reacts with the surface. It has been demonstrated that the reaction proceeds mostly with the N2 bond oriented and N2 positioned as in the minimum barrier B2H dissociation geometry. We have also evaluated the vibrational efficacy of N2 and showed that adding vibrational energy to initial gas-phase states promotes the N2 reaction more efficiently than increasing translational energy.