Vibrational excitation in plasma catalysis: how important are dynamical effects?

Abstract

Plasma catalysis offers a promising alternative to current ammonia production processes, due to the combination of high selectivity of heterogenous catalysis and efficient activation of nitrogen in the plasma. However, the theoretical understanding of how various plasma processes contribute to efficiency improvements remains limited. The pioneering work of Metha et al. (Nat. Catal., 2018, 1, 269) extended the standard formulation of transition state theory by making it vibrational state-specific through the use of the Fridman–Macheret α model. The resulting microkinetic model accounted for vibrational contributions under the non-equilibrium conditions of a plasma reactor. In this work, we critically examine the prototypical chemical process of activated N₂ reactivity on ruthenium through explicit rate coefficient calculations using state-of-the-art molecular dynamics, based on a potential energy surface previously validated against molecular beam experiments. Our findings reveal that vibrational activation is significantly more effective in promoting surface reactivity than predicted by the Fridman–Macheret α model, which fails to capture the full complexity of state-specific contributions. Furthermore, our calculations indicate that vibrational activation is also the primary driver of highly activated thermal catalytic reactions. These results provide a valuable benchmark to guide the development of future state-specific microkinetic models for heterogeneous and plasma catalysis.

---

Broader context

Pre-activating reactants by a plasma holds great potential to boost heterogeneously catalysed processes, especially in scenarios where dissociative chemisorption (DC) reactions are rate-limiting. However, the fundamental atomic-scale understanding and quantification of this enhanced performance are still lacking. The pioneering work on plasma-enabled ammonia synthesis for DC of N₂ introduced the Fridman–Macheret (FM) α model, which is widely used to extend standard micro-kinetic modelling to account for vibrational excitations in plasma catalysis. While this model is computationally inexpensive, it relies on assumptions based on gas-phase chemistry. In this article, we scrutinise the FM-α model for N₂ dissociation on Ru(0001) using explicit molecular dynamics simulations on an accurate potential energy surface. We find that vibrational excitation plays a far more dominant role in N₂ dissociation than the FM-α model predicts, and the reasons behind this discrepancy are analysed. This is a crucial step towards developing accurate effective models to replace the FM-α model.

1 Introduction

The synthesis of ammonia is a pivotal chemical reaction for the world's food supply. To date, it has been dominantly driven by the Haber–Bosch process,¹ with the dissociative chemisorption of N₂ on a catalyst's metal surface being the rate controlling step.^2,3 This elementary reaction step has been the subject of many studies – both experimental^4–6 and theoretical^7–13—and ruthenium has been found to be the optimal single element catalyst for the Haber–Bosch process.^14–17 Despite optimisations, this process's power consumption remains high even under optimal reaction conditions. Current efforts to improve the energy efficiency of the process involve applying plasma-enhanced catalysis. This approach aims to overcome the kinetic limitations of traditional thermal (heterogeneous) catalysis, potentially significantly reducing the operational temperature and pressure.¹⁸ With reactants being excited in the plasma, numerous (synergistic) effects can increase reactivity, such as rovibrational and electronic excitation, ionisation and dissociation, modification of the catalyst surface through etching and charging, and the presence of electric and magnetic fields.^19–21 Moreover, the combination of a non-equilibrium plasma with a heterogeneous catalyst has been demonstrated to be able to surpass the sum of its parts.^22,23 However, plasma catalysis is also considerably more complex than heterogeneous catalysis. The aforementioned plasma-induced effects might all play a role in plasma catalysis, or some not at all, depending on the operating conditions, the catalyst, and the reactor design.²⁰ This makes it challenging to investigate and disentangle these effects through experiments alone. Therefore, simulations are necessary to elucidate the role of individual plasma-induced effects and thus enhance our understanding of the entire process.

Vibrational excitation of N₂ is known to increase the dissociative chemisorption (DC) rate.^24–27 The ground-breaking work of Mehta et al.²⁸ for plasma-enhanced ammonia production has made the DC of N₂ a prototypical reaction and an important showcase for the importance of vibrational excitations in plasma catalysis.²⁰ This notion has been reached based on simulations for the turnover frequency (TOF) for ammonia synthesis under realistic catalytic reaction conditions. The underlying microkinetic models (MKMs) rely on the Fridman–Macheret (FM) α model²⁹ to quantify the effect of vibrational excitation on the reaction rate coefficient of the DC of N₂.^28,30 The FM α model is rooted in gas-phase reactions, where it models the vibrational contribution in diatom–atom reactions by computing the relative ratio of the forward and backward barrier heights. This model has become the workhorse approach for modelling vibrational excitation in plasma catalysis.^31,32

Recently, questions have arisen regarding the relative importance of vibrational excitation compared to other plasma effects.^21,33 For instance, the concentration of vibrationally excited molecules is typically much lower compared to the concentration of radicals in dielectric barrier discharges, which are the plasmas most commonly used in plasma catalysis.²¹ It is important to note that the FM α parameter cannot be directly measured experimentally. Instead, it requires models fitted to experimental data or theoretical calculations to validate its value. While a simple gas-phase reaction may be reasonably well described by a model like TST and the FM α, molecule–metal surface reactions typically involve significantly more complex potential energy surfaces (PESs) and necessitate high-dimensional models. The applicability of the FM α model to these complex reactions has already been questioned by Kedalo et al.³⁴ for N₂ + Ru(0001) based on comparison with MD simulations. While their study should be considered with caution (see Section 3.3 and Section S2.2 in the ESI†), they found that the effect of vibrational excitation exceeds the FM α model's predictions but is not as significant as the effect of translational energy. Recently, some of us have suggested that the FM α approach performs poorly, both qualitatively and quantitatively, for DC in general. This is based on an extensive analysis of theoretical molecular dynamics (MD) studies available in the literature.³⁵ For instance, the effect of the vibrational excitation of polyatomic molecules on the reactivity depends on the specific vibrational mode(s) being excited.^35–40 For methane, excitation of the vibrational stretch modes is more effective than exciting the bend modes.^36–38 Similarly, exciting overtones generally leads to complex distributions of near-degenerate vibrational states that again yield considerably different reaction probabilities.^39,40 Furthermore, the FM α model seems to significantly underestimate the effect of vibrational excitation for a variety of DC reactions.³⁵ In late barrier systems, the bobsleigh effect can make molecules with high incidence energy deviate from the minimum energy path (MEP), and consequently they experience much higher barrier heights, lowering reactivity^40–42 (also often referred to as Polanyi's rules⁴³). Additionally, sterical hindrance in the MEP combined with dynamical effects can have a major influence on the effect of rotational and vibrational excitation of reactants.⁴⁴ Unfortunately, the FM α model is not able to capture these effects, since it only considers the forward and backward reaction barrier heights of the ground state PES and the absolute vibrational energy.

MD approaches are able to capture these complex effects of vibrational excitation, offering an intriguing alternative to the FM α model.⁴⁵ In this study, we explore how and why the FM α model and some simple extensions deviate considerably from full-scale MD simulations. For the latter, we build on our previous work that established an accurate MD model for N₂ dissociation on a Ru(0001) surface, meticulously validated against the gold standard in gas-surface dynamics provided by molecule-beam experiments.^11,12,46 This model employs quasi-classical trajectory (QCT) calculations using a machine-learned high-dimensional potential energy surface (PES) based on DFT calculations, including all relevant degrees of freedom of both the metal surface and the molecule. We utilise this model to compute the effect of vibrational excitation on the reaction probabilities and uncover a nontrivial relationship between vibrational excitation and reactivity. Crucially, we show that vibrational excitation plays a significantly more important role in plasma catalysis than previously anticipated, influencing the predicted TOF of ammonia in plasma catalysis. If the rate-controlling reaction step is highly activated, our findings suggest that the same can hold for conventional thermal catalysis. This could have serious implications for the modelling of heterogeneous catalysis; to the best of our knowledge, vibrational excitation is consistently neglected in MKMs for thermal catalysis. Finally, we analyse the reaction dynamics to elucidate the disparity in reactivity between conventional TST methods and our dynamics-based approach, concluding that dynamical effects play an important role in plasma catalysis.

2 Methods

2.1 Microkinetic model

Our study builds on the MKM originally developed by Mehta et al.²⁸ for modelling plasma-enabled catalysis. In this MKM, each vibrational level is treated as a distinct reactant, with its own specific reaction rate and a concentration directly proportional to the vibrational distribution. We consider three distributions for the vibrational states of the reactant N₂ molecules:

(1) A ‘ground-state-only’ distribution (ν = 0), as typically used to model thermal catalysis.

(2) A Boltzmann distribution to account for the population of vibrational states in thermal equilibrium at a given vibrational temperature T_vib that is the same as the gas temperature T_gas.

(3) A Treanor distribution⁴⁷

to describe the non-equilibrium population of vibrational states, which is commonly used in plasma catalysis.^28,48,49 The Treanor probability density f(ν; T_vib, T_gas) is a function of the vibrational state ν at given temperatures T_vib and T_gas, which usually are significantly different from each other. B is a normalization constant, ω is the vibrational frequency, and x_e is the anharmonicity coefficient.

Throughout this work, we use this MKM to compute relative TOFs of ammonia synthesis, which are independent of pressure. This highlights the fact that different reaction rate coefficients for the N₂ dissociation step cause the same relative differences for the TOFs for both the industrial high-pressure (100 atm for the Haber–Bosch process) and the typical low-pressure conditions in plasmas (1 atm). We assume a vibrational temperature of T_vib = 3000 K for the plasma, like in the work of Mehta et al.²⁸ As mentioned in Section 1, other plasma-induced effects (e.g., the presence of radicals) can also play a role, which we neglect here, in order to focus on the effect of vibrational excitation.

2.2 Reaction rate coefficients

2.2.1 Transition state theory. Usually, reaction rate coefficients in an MKM are obtained based on transition state theory (TST) in the form of an Arrhenius equation:

In this equation, k_B is the Boltzmann constant, E_a is the activation barrier, including the effects of zero-point energy, and A is a frequency factor related to the entropy barrier and pressure. All activation barriers E_a in the work of Mehta et al.²⁸ are determined by scaling relationships (SR), i.e., using fitted trends between adsorption energy and reaction barriers across many metal surface species.⁵⁰ Focussing on N₂ dissociation in this work, we opt to use a more accurate activation barrier of E^TS_a = 1.83 eV on the basis of a DFT-based transition state (TS) search^11,51 for this rate-limiting step. This is a substantially higher barrier than the E^SR_a = 1.48 eV barrier from the SR. For all other elementary reaction steps, we keep the SR-based values as determined by Mehta et al.²⁸ for terraces on the Ru(0001) surface unchanged.

2.2.2 Fridman–Macheret α models. To make the rate coefficient for N₂ dissociation vibrational state specific, Mehta et al.²⁸ have suggested to use the Fridman–Macheret (FM) α model, which has originally been developed for gas phase reactions. According to this model, reaction rate coefficients are obtained from a simple extension of the Arrhenius equation:

Here, H is the Heaviside step function, and most importantly, the FM α quantifies how effectively the vibrational energy E_ν (excluding zero-point energy) reduces the effective barrier height. The FM α model aims to enforce the Polanyi rules⁴³—which state that, in a late-barrier system, vibrational excitation has a greater impact on reactivity—by linking the ratio between the forward and backward reaction barriers to the “lateness” of the barrier and, consequently, to α as follows:

Here, the forward and backward reactions are the DC and desorption of N₂, respectively, and can be obtained with quantum chemistry calculations, such as DFT, either directly through a transition state (TS) search, or indirectly using the SR. Using scaling relationship activation energies results in an FM α^SR of 0.37 and using the DFT activation energy for the forward dissociation reaction (keeping the adsorption energy equal) a slightly larger α^DFT = 0.39 is obtained. We refer to this rate calculation method as the ‘TST + FM’ level of theory, appending ‘@SR’ or ‘@TS’ when using the activation energy of N₂ dissociation from scaling relationships or a DFT transition state search, respectively.

In the field of gas-surface dynamics, vibrational efficacies (VEs)

can be considered the equivalent of α in the FM α model, i.e., they quantify how much more efficiently vibrational energy E^ν_vib increases the reaction probability P relative to the same amount of translational incidence energy E^ν_inc (taken to be perpendicular to the surface) for a molecule in vibrational state ν. While VEs have been measured for the dissociative chemisorption of a number of molecules,^52,53 no accurate data are available for N₂, since a suitable experimental technique to prepare molecular beams in specific vibrational states is yet to be found. To assess the accuracy of the FM α, we have obtained a mean [small eta, Greek, macron] (detailed in Section 2.4 and Section S4 in the ESI†) using N₂ dissociation probability curves P_ν(E_inc) as a function of incidence energy E_inc per vibrational state ν calculated through MD simulations. This [small eta, Greek, macron] we use as a substitute to the FM α in eqn (3) in order to compute rate coefficients. We refer to this method as the ‘TST + [small eta, Greek, macron]@TS’ level of theory in the following.

2.2.3 Kinetic gas theory. The last method for obtaining vibrational-state-specific rate coefficients forgoes the Arrhenius equation entirely by making use of the expressions derived from kinetic gas theory.⁵⁴ Instead of approximating the reactivity with a single barrier height, this requires reaction probability curves P_ν(E_inc) from MD to define rate coefficients:

k_ν = A〈P_ν〉(T_gas), (6)

where A is a frequency constant as in eqn (2) and (3), and 〈P_ν〉(T_gas) is the ensemble-averaged reaction probability at temperature T_gas, which is given by:

In this integral, f_inc(E_inc; T_gas) is the probability density of the 1D Maxwell–Boltzmann distribution for the incidence energy E_inc according to the gas temperature T_gas given by:

Here, , where m is the molecular mass of an N₂ molecule and v_inc is its center of mass's velocity component perpendicular to the surface. We calculate the N₂ dissociation probability P_ν(E_inc) using MD simulations, as detailed in Section 2.3.3. Rate coefficients obtained with this method are referred to as the ‘MD’ level of theory in the following.

2.3 Molecular dynamics simulations

2.3.1 High-dimensional neural network potential. All MD simulations were performed using the high-dimensional neural network potential (HDNNP) constructed by Shakouri et al.¹¹ according to the Behler–Parinello approach.⁵⁵ Briefly, this potential is based on 25 000 single-point RPBE-DFT⁵⁶ calculations sampling the interaction of N₂ molecules with the Ru(0001) surface, which is modelled by a slab of 7 layers using a 3 × 3 supercell, where the bottom layer is kept fixed throughout all simulations. Of these 25 000 single points, 5000 were for a relaxed surface at 0 K surface temperature, and the remaining 20 000 were generated taking lattice expansion and surface atom displacements into account. In this work, the potential is evaluated in the LAMMPS code⁵⁷ using the ML-HDNNP package's interface to the n2p2 library.^58,59

2.3.2 Initial conditions. The initial height of N₂ above the surface is 5.5 Å. Because normal energy scaling has been observed for N₂ dissociation probabilities on Ru(0001),⁶⁰ we only consider perpendicular incidence with a center of mass velocity . The other initial conditions of the N₂ molecule (i.e., orientation, phase of the vibration, and the position relative to the unit cell) are sampled quasi-randomly according to a low-discrepancy sequence as defined in ref. 61. The same set of initial conditions (θ, ϕ, x, y, and vibrational phase within its period) is used for all incidence energies and vibrational states. The surface temperature for this simulation was 575 K and we expanded the lattice constant of the slab according to experimental thermal expansion factors.⁶² A Nosé–Hoover canonical ensemble (NVT) simulation using LAMMPS for a total of 150 ps with a 0.5 fs time step and a 50 fs damping time is used to extract 100 different snapshots of surface configurations with displaced atoms.

2.3.3 Quasi-classical trajectory method. The initial conditions are sampled separately by MD trajectories at different incidence energies according to the QCT approach.⁶³ In the QCT method, the N₂ molecule is initialised with vibrational energy equal to the eigenenergies of the Hamiltonian of the free diatom. The vibrational states and corresponding energies are computed from the diatom potential via the Fourier grid Hamiltonian method.⁶⁴ Further details of the QCT method for N₂ on Ru(0001) are provided in ref. 12. The dissociation probabilities were calculated for incidence energies ranging from 0.25 to 10.0 eV in steps of 0.25 eV and for vibrational quantum numbers ν ∈ {0…10}. We only simulate the rotational ground state J = 0, because previous work has shown that rotational excitation has limited to no effect on the reaction probability of N₂ on Ru(0001).⁹ The number of trajectories used for computing the reaction probabilities was 10⁴ (10⁵) for total energies (i.e. incidence and vibrational energy, including the zero-point energy) larger (smaller) than 2.25 eV. All MD trajectories are calculated with the LAMMPS code using the standard NVE integrator with a time step of 0.33 fs.

2.3.4 Reaction probabilities. The reaction outcome of each MD trajectory is determined by the following criteria: the N₂ molecule is considered to have reacted (dissociatively) if its bond length is larger than 2.65 Å. The molecule is considered to be scattered if it is further than 5 Å away from the surface with the centre-of-mass velocity vector pointing away from the surface. If neither occurs within the maximum simulation time of 36.3 ps, the molecule is considered to be trapped on the surface. Since most (99.9%) reactive trajectories occurred within the first 3.3 ps, we neglect the possibility that trapped molecules might react (dissociatively) after 36.3 ps and count them as non-reactive. We calculate N₂ dissociation probability P_ν(E_inc) using an estimator:

Here, the total number of trajectories is the sum of dissociated, scattered and trapped MD trajectories, respectively, at each E_inc and vibrational quantum number ν (only explicitly denoted on the left-hand side of eqn (9) for brevity).

We have compared different approaches to numerically evaluate the integral in eqn (7) using the results from our MD calculations (see Section S2 in the ESI† for details). The most robust method uses the carefully constructed fitting function

to obtain a continuous representation of p̂_ν(E_inc). Here, α_ν, β_ν, and γ_ν are the optimal fitting parameters per vibrational state ν, which we have obtained from a least-squares fit as tabulated in Section S2.1 of the ESI.†

The HDNNP was designed with an emphasis on configurations where the total initial energy of the N₂ molecule (both translational and vibrational components) and the Ru surface does not exceed 15 eV. Consequently, the chance increases that the HDNNP needs to extrapolate towards parts of the PES that are less-well covered by the original DFT data set for MD trajectories at higher initial energies. We account for these uncertainties in p̂_ν(E_inc) and propagate them in the form of 95% confidence bounds. Using the exact Clopper–Pearson method,⁶⁵ we compute the 95% confidence intervals for the estimators. For the lower bound, we include the number of non-conclusive trajectories terminated because of extrapolation N_extr in the total as if they are non-reactive (N⁺_total = N_diss + N_scat + N_trap + N_extr), such that the lower bound is from the 95% confidence interval of a binom(N_diss, N⁺_total) distribution. Likewise, for the upper bound, we also include the inconclusive extrapolative trajectories in the reactive trajectories N⁺_diss = N_diss + N_extr, such that the upper bound is from the 95% confidence interval of a binom(N⁺_diss, N⁺_total) distribution. This allows both the uncertainty in the estimator due to finite sampling and the outcome of the trajectories due to the extrapolation of the HDNNP to be accounted for (see Fig. S1 in the ESI†). These errors are dominated by the former (latter) at low (high) incidence energies. Uncertainties in the resulting reaction rate coefficients k_ν and TOFs from the MKM calculations are accounted for by distinct fits of eqn (10) to the lower bound and upper bound of the estimators, giving us effective 95% confidence intervals (p₋ and p₊) for a given fit parameter p. Additionally, we account for the uncertainty introduced by the least-squares procedure σ_± at the confidence bounds p_± by expanding the individual bounds as follows:

Here, z = 1.96 corresponds to 95% confidence intervals of the least-squares uncertainties and we assume that the fit uncertainty is independent of the input data uncertainty allowing us to sum the variances associated with the confidence intervals. This results in 3 sets of parameters for eqn (10) per vibrational state: a lower bound, an estimator and an upper bound (see Section S2.1 in the ESI†), which are used to compute the error bars for results at the ‘MD’ level of theory. Although the least-squares uncertainties are substantial for high vibrational quantum numbers, they ultimately hardly affect the uncertainties of the derived k_ν and TOFs.

2.4 Vibrational efficacy calculations

As described in Section 2.2.2, we have used a mean VE [small eta, Greek, macron] to substitute the FM α in eqn (3) for vibrational state-specific rate coefficients. However, η_ν(P) as defined in eqn (5) depends on the vibrational state-specific reaction probability and is not trivially reduced to a single scalar [small eta, Greek, macron]. In this section, we will discuss methods for computing [small eta, Greek, macron] as done in the literature and compare to our approach. One method is to pick a fixed reaction probability and take the mean over the vibrational quantum number such that

This method considers only a single point of the reaction probability curves, which needs to be chosen carefully. In previous work for N₂ dissociation on Ru(0001), P has been chosen between 0.1 and 0.5, resulting in [small eta, Greek, macron] ≈ 1.6 for ν ≤ 3.^10,12

Here, we extend the original approach by aligning and incorporating the reaction probability curves for different vibrational states over a broader range of incidence energies. To achieve this, we express the numerator in eqn (5) as a sum of shifts between the successive vibrational states, . The δE^ν_inc values are obtained as optimal shifts by minimising the total square difference integrals

separately for each vibrational state ν. The ν-specific vibrational efficacies are then given by:

Finally, the mean value [small eta, Greek, macron] is obtained as the slope of a linear fit through the points (ΔE^ν_vib and ΔE^ν_inc). For the details of this procedure and estimation of its uncertainty, see Section S4 in the ESI.†

3 Results: ammonia synthesis on Ru(0001)

3.1 Turnover frequency for ammonia production

In Fig. 1, we present the relative TOFs of ammonia synthesis with the DC rate of N₂ computed at the aforementioned separate levels of theory: TST + FM@SR, TST + FM@TS, TST + [small eta, Greek, macron]@TS, and MD defined in Section 2.2. Additionally, for each level of theory, TOFs are calculated for all three vibrational distributions (see Section 2.1), i.e., a ground state (T_vib = 0 K), thermal (T_vib = 673 K) and plasma distribution (T_vib = 3000 K). The TOFs presented in Fig. 1 are given relative to the reference ammonia TOF reported as ‘plasma-off’ results by Mehta et al.,²⁸ which, in our notation, corresponds to TST + FM@SR with a ground-state vibrational distribution. To ensure an adequate comparison, the vibrational temperatures and the gas temperature (T_gas = 673 K) have been chosen identical to this previous work.

Fig. 1 Relative turnover rates for ammonia production on Ru(0001) under industrial reaction conditions for different levels of theory using the ‘plasma-off’ result from the study by Mehta et al.²⁸ as the reference. For each level of theory, we report the TOF for the vibrational ground state (T_vib = 0 K), a thermal gas (T_vib = 673 K) where the vibrational states are Boltzmann distributed, and a plasma (T_vib = 3000 K) where the vibrational states are Treanor distributed. T_gas is 673 K in all three cases. The levels of theory are (i) SR activation energies with the FM α model as described in ref. 28, (ii) DFT-calculated activation energy^11,51 with the FM α model, (iii) DFT[-calculated] activation energies with the mean VE [small eta, Greek, macron] from QCT and (iv) explicit reaction probabilities computed with MD using the QCT method.

At this level of theory, Mehta et al.²⁸ demonstrated that the plasma vibrational distribution produces a significantly higher TOF—four orders of magnitude greater—compared to the ground-state-only result. Here, we show that the TOF increases by only 3% when a thermal distribution of N₂ vibrational states is considered, as opposed to only the ground state. In other words, TST + FM@SR predicts that vibrationally excited species contribute minimally to ammonia production under thermal conditions—consistent with common expectations. As mentioned above, the first refinement of the reaction rate (i.e., TST + FM@TS) involves computing the barrier height directly instead of relying on scaling relationships. In general, the TOFs for all vibrational distributions have decreased roughly by the same amount (almost 3 orders of magnitude), which is directly caused by the 0.35 eV increase in activation energy. Although this increase of E^forward_a also increases the FM α (|E^forward_a| ≫ E^backward_a in eqn (4)), the vibrational enhancement of the reactivity is negligible in the ground state and for thermal conditions, and only slightly larger under plasma conditions compared to TST + FM@SR.

At the next higher level of theory, TST + [small eta, Greek, macron]@TS, we observe a qualitative shift, with the vibrational ground state and thermal distributions now differing by more than 5 orders of magnitude. This suggests that, when the FM α = 0.39 is replaced by the much higher [small eta, Greek, macron] = 1.8, vibrationally excited species dominate the TOF in the MKM not only under plasma conditions but also under thermal conditions. Given the increased contributions of vibrationally excited N₂, one might expect a larger relative difference in TOFs between thermal and plasma distributions. However, the effect remains similar, with the TOF increasing by 5 orders of magnitude compared to 4 orders of magnitude previously. In Section 3.2, we will show that, at higher vibrational levels, the reaction becomes effectively ‘barrierless’, leading to a saturation of the reaction rate.

Finally, at our highest level of theory—rate coefficients from MD—the difference between the ground state and thermal distributions is larger than in the FM α-based results but remains considerably smaller than the [small eta, Greek, macron]-based results. Conversely, the gap between the plasma distribution and the ground state or the thermal distribution is immense, spanning 10 orders of magnitude—a difference unlikely to be observed in experiments. However, we argue that the MD simulations capture a realistic trend, as discussed in Section 3.2. Instead, the likely overestimation of the population of highly excited vibrational states by the Treanor distribution, as also reported in ref. 30, may account for this discrepancy. This underscores how the partitioning of energy can have a substantial impact on the predicted reaction rates in plasma catalysis. Finally, the ground state TOF at the MD level of theory is significantly lower than at the TST-based levels, indicating a higher effective activation energy for N₂ dissociation. In Section 3.4, we will show that this effective increase arises from dynamical effects.

3.2 Analysis of N₂ dissociation probability

To better understand the overall ammonia TOFs computed via the MKMs, we analyse the DC of N₂ on Ru(0001) in detail, as this elementary step is known to have a pronounced effect for the conditions considered in the previous section (T_gas = 683 K). Since the prefactor A is identical in eqn (3) and (6), relative differences between the reaction rate coefficients k_ν are best compared by the thermally averaged reaction probabilities. In Fig. 2, we present 〈P_ν〉 as a function of the vibrational quantum number ν for the four levels of theory previously introduced and discussed. This figure reasserts that the change in activation energy between the SR and the explicit TS only induces a shift in the N₂ reaction probability. On a logarithmic scale, the exponent in eqn (3) simplifies to a linear dependence on E_ν—consequently on ν when anharmonicity is negligible. Therefore, the nearly identical slopes for TST + FM@SR and TST + FM@TS explain the negligible differences in the effect of vibrational excitation on the ammonia TOFs between these two levels of theory in Fig. 1.

Fig. 2 Thermally averaged reaction probabilities per vibrational level at different levels of theory: (i) activation energies from the SR²⁸ with the FM α model²⁹ (red stars), (ii) activation energy for N₂ dissociation from a DFT TS search^11,51 with the FM α model (green triangles), (iii) the same activation energies with the mean VE from MD (orange squares), and (iv) reaction probability explicitly computed using MD.

When the FM α is substituted with [small eta, Greek, macron], the slope is notably increased, reaching 〈P_ν〉 = 1 already for vibrational states ν ≥ 4. In other words, according to the FM + [small eta, Greek, macron]@TS level of theory, N₂ DC on Ru(0001) is effectively barrierless and thus is guaranteed to occur starting from the 4^th vibrationally excited state. For both TST + FM@SR and TST + FM@TS, the first vibrationally excited state ν = 1 is only one order of magnitude more reactive than ν = 0, much less than the decrease in population between the two states for a thermal Boltzmann distribution. On the other hand, for TST + [small eta, Greek, macron], ν = 1 is four orders of magnitude more reactive, more than enough to allow this state dominate reactivity over the vibrational ground state despite its lower population. The individual contributions of each vibrational state to the total reaction probability is further discussed in Section S3 of the ESI.† All of this explains why the increase in TOF from a thermal to a plasma distribution in Fig. 1 for TST + [small eta, Greek, macron]@TS is smaller than expected a priori: the reaction rate saturates at high vibrational quantum numbers, which is precisely where the plasma distribution differs most significantly from the thermal distribution in terms of vibrational state populations. Nevertheless, the substantial impact of the plasma on the ammonia TOF is maintained, as the total number of vibrationally excited species remains higher in the plasma distribution compared to the thermal one.

Finally, MD-based reaction probabilities are the only ones that strongly deviate from a linear dependence on the vibrational level. While they initially exhibit a slope similar to [small eta, Greek, macron], this slope gradually decreases with increasing vibrational quantum number, appearing to asymptotically approach 〈P_ν〉 = 1. This change in slope explains the smaller difference between the TOFs for the ground state and thermal distribution, as only the first few vibrational quantum numbers contribute significantly to the N₂ dissociation rate (see Section S3 in the ESI†). The fact that the TOFs in the vibrational ground state at the TST + FM@TS and TST + [small eta, Greek, macron]@TS levels of theory are more than three orders of magnitude higher than those from MD is a direct consequence of the corresponding decrease in reaction probability at ν = 0 by a similar amount.

3.3 Vibrational efficacy

As outlined in Section 2.4, the mean vibrational efficacy, denoted by [small eta, Greek, macron], is determined by aligning reaction probability curves for subsequent vibrational levels, ν, and averaging over the resulting η_ν. Notably, ν-specific vibrational efficiencies, η_ν, exhibit minimal variations with respect to ν (as detailed in Section S4 in the ESI†). This behaviour is similar to the FM α parameter, which remains constant for each vibrational level. However, the mean vibrational efficacy calculated from our MD simulations, [small eta, Greek, macron] = 1.805, is significantly higher than the FM α = 0.387. [small eta, Greek, macron] > 1 implies that vibrational energy is more efficient than the same amount of translational energy in promoting the dissociation of N₂ on Ru(0001). In other words, vibrational excitation through plasma excitation is more efficient than thermal heating of the gas. This phenomenon cannot be captured by the FM α model, which, by definition (see eqn (4)), restricts α to the range of 0 to 1.

Going beyond the alignment of entire reaction curves, our molecular dynamics simulations enable us to delve deeper into the dependence of η_ν on the incidence energy according to eqn (5). This analysis allows us to examine the vibrational efficiencies in greater detail. Fig. 3a illustrates that the ν-specific vibrational efficacies are indeed significantly smaller at low incidence energies, exhibiting an increase with each vibrational level. Notably, η_ν(E_inc) rapidly increases until approximately E_inc ≈ 3 eV, where they attain a steady value of η_ν ≈ 2 for each vibrational quantum number. This observation suggests that the thermally averaged state-specific vibrational efficacies

strongly depend on the vibrational state. Fig. 3b reveals that the average efficacy 〈η_ν〉 varies from 〈η_ν=1〉 = 0.14 to 〈η_ν=10〉 = 1.0 for a thermal gas at 1000 K, exhibiting an almost linear increase with the vibrational quantum number. Notably, all these values are significantly lower than the mean efficacy [small eta, Greek, macron], because the thermal averages are primarily influenced by lower incidence energy contributions (see the inset of Fig. 3a).

Fig. 3 (a) VE as a function of incidence energy per vibrational state, with the inset zooming into 0 eV ≤ E_inc ≤ 0.5 eV. (b) Comparison of the vibrational state independent FM α^FM (eqn (4), orange line) and mean VE [small eta, Greek, macron] (eqn (14), blue line), as well as the state-specific VE 〈η_ν〉 (T_gas = 1000 K) (eqn (15), blue circles) and MD-derived α^MD_ν(T_gas = 1000 K) (eqn (16), orange circles). Results for α^MD_ν(T_gas = 1000 K) as obtained by Kedalo et al.³⁴ are also included (black crosses).

Previously, Kedalo et al.³⁴ proposed an approach to extract FM α values from the reaction probability curves generated through MD simulations. This involves solving for α in eqn (3), which yields

We have included the results from ref. 34 in Fig. 3(b) (black crosses) for α^MD_ν(T_gas = 1000 K) for N₂ dissociation on Ru(0001). In addition, we have also applied eqn (16) to the thermally averaged reaction probability obtained from our own MD simulations according to eqn (7) for the same temperature (orange circles in Fig. 3b). The results for α^MD_ν(T_gas = 1000 K) from the two different sets of calculations differ quite notably from each other. Most importantly, as ν increases, our results do not approach α^FM but rather converge to zero. We think that systematic shortcomings in the methodology of Kedalo et al.³⁴ are causing these differences: firstly, their barrier height for N₂ dissociation on Ru(0001) computed with the PBE density functional (1.9 eV) is larger than our value underlying our HDNNP computed with the RPBE density functional (1.84 eV¹¹). This is in large disagreement with the SBH17 database,⁵¹ where the PBE and RPBE barrier heights have been computed as 1.49 eV and 1.95 eV, respectively. The values from the SBH17 follow the common trend that RPBE yields larger barrier heights for molecule–metal surface reactions than PBE. The small deviation of the aforementioned RPBE values (0.1 eV) is caused by slight differences in the DFT setup used in both studies, whereas the much larger difference between the PBE values (0.4 eV) is clearly beyond this error margin. Secondly, the root mean square error (RMSE) for the forces of the neural network fit for the potential energy surface reported by Kedalo et al.³⁴ (0.2–0.4 eV Å⁻¹) is an order of magnitude larger than what is usual (0.02–0.03 eV Å⁻¹).⁴⁰ Thirdly, in their MD simulations, Kedalo et al.³⁴ shift the Maxwell–Boltzmann distribution for the velocities to considerably higher energies, effectively assuming that 80% of the total (translational) velocity of N₂ molecules is distributed into the component perpendicular to the surface (see Fig. S4 in the ESI†). Furthermore, their MD simulation yields higher N₂ dissociation probabilities for the flat Ru(0001) surface than that for the stepped Ru(113) surface—opposite to what has been obtained in other computational work.^7,66,67

In light of these concerns regarding the DFT calculations, machine-learned potential and MD simulations of Kedalo et al.,³⁴ we now focus our comparison of α^MD_ν with 〈η_ν〉 on the results obtained from our own MD simulations (orange circles and blue squares in Fig. 3(b)). Starting from α^MD₁ ≈ 1.4, there is a monotonous decrease with vibrational quantum number ν, which still leaves α^MD₂ significantly larger than the FM α given by eqn (4). The fundamentally opposite trends for 〈η_ν〉 and α^MD_ν with respect to ν arise from the different order of operations. For 〈η_ν〉, one first calculates the difference between the reaction probability curves and then integrates this difference over the velocity distribution. The definition of 〈η_ν〉 in eqn (15) consistently generalises the standard definition of vibrational efficacies used in gas-surface dynamics, completely analogous to 〈P_ν〉 in eqn (7). However, we do not consider 〈η_ν〉 as a suitable substitute to introduce ν-dependence to the FM α, as this would effectively double-count statistical effects in the TST and MD descriptions. On the other hand, for α^MD_ν, one first integrates P_ν(E_inc) over the Maxwell–Boltzmann distribution and then calculates the difference between the thermally averaged reaction probabilities 〈P_ν〉. Consequently, vibrational-state-related effects are not blended with a combination of (different) statistical treatments at the TST and MD levels. By construction, if the results for α^MD_ν were to be substituted for the constant [small eta, Greek, macron] in eqn (3), the resulting TOFs would closely align with the MD reference results shown in Fig. 1. Clearly, α^MD₁(1000 K) > 1 is the key enhancement that cannot be achieved with the FM α model.

Altogether, we present compelling evidence that capturing the impact of vibrational excitation on the reaction rate coefficient using the FM α approach, which currently stands as the state-of-the-art in plasma catalysis modelling, presents substantial challenges. According to the Polanyi rules,⁴³ this effect can be estimated based on the potential energy surface. Early (late) barriers typically exhibit low (high) vibrational enhancement. Therefore, we anticipate similar implications for the efficiency of plasma excitation for other highly activated late barrier reactions, such as CH₄,^40,68,69 HCl,^44,70 and CH₃OH,^71,72 where high VEs have also been observed.

3.4 Dynamical effects

We continue our analysis by focusing on the TOFs for the vibrational ground state as presented in Fig. 2. The results based on reaction probabilities obtained from our MD simulations are significantly lower than those based on TST-based methods. This difference is likely caused by dynamical effects that are not accounted in any of the latter approaches. The conventional formulation of TST, as employed here, primarily considers the minimum energy path (MEP) for the dissociation of N₂, with approximate vibrational corrections. This MEP is defined by a single reaction coordinate, which, in the simplest case for DC, encompasses only the bond length and the molecule–surface distance. In contrast, our MD simulations include a much higher number of degrees of freedom: 6 for the N₂ molecule and 162 resulting from the 54 moving atoms in the slab model for the Ru(0001) surface. Their dynamical interplay can facilitate dissociation along paths that differ substantially from the MEP.

To analyse these dynamical effects, we extracted statistics from MD trajectories of dissociated N₂ molecules in the vibrational ground state at a high incidence energy (E_inc = 3.25 eV) above the barrier along the MEP on our PES (E^TS_a = 1.83 eV). Even at such a high incidence energy the dissociation probability obtained from the MD simulations is less than 100% (see Section S1, ESI†). In other words, not all molecules are efficiently steered towards the MEP and can thus encounter barriers higher than E^TS_a. In our statistical analysis, we focus on the coordinates depicted in the insets of Fig. 4 that are traditionally not included in the reaction coordinate: the angular orientation of the molecular axis of N₂ by the polar angle θ and its lateral position (x and y) relative to the surface unit cell, all based on the N₂ centre of mass as the origin. As described in Section 2.3.2, our trajectory ensemble is initially uniformly distributed over both orientations and lateral positions. According to TST, all reactive trajectories will cross the same minimal energy barrier along the MEP, which implies that there should be no preferential initial conditions. In Fig. 4, we show the distributions for the initial orientation θ and xy position of the subset of dissociated N₂ molecules, which have been rescaled such that their maxima coincide with those of the distributions for all trajectories. A sharper distribution for reacted trajectories indicates a degree of inaccessibility of the MEP for initial conditions further from the MEP, leading to a lower reaction probability compared to TST. By computing the ratio between the integrals over these (rescaled) distributions and their original counterparts for all trajectories, we can estimate the reduction in reactivity in particular degrees of freedom due to dynamical effects.

Fig. 4 Distributions for initial (a) orientation θ and (b) xy position of N₂ molecules above the Ru(0001) surface for E_inc = 3.25 eV (>E^TS_a = 1.83 eV) and ν = 0. The coordinates are depicted in the insets. (a) The initial distribution for all trajectories (both reacted and non-reacted, black dashed line, eqn (S6) in the ESI†) is compared to the distribution of the subset comprised by the reacted trajectories (blue bars). To guide the eye, the latter has been fitted to the initial distribution multiplied by a Gaussian centred at θ = 90° (orange line, eqn (S7) in the ESI†). (b) The distribution for the initial xy-position of reacted trajectories in the irreducible wedge of the unit cell, as depicted in the inset. a denotes the surface lattice constant. Individual reacted trajectories are in red, the histogram is in greyscale, and the fitted bimodal Gaussian distribution is the contour lines. For better comparison and to quantify dynamical steering effects, the reacted distributions have been rescaled such that their maxima coincide with those of the distributions for all trajectories.

As shown in Fig. 4a and b, respectively, both the θ and xy distributions for the reacted molecules deviate considerably from the uniform distribution of all molecules. The very narrow distributions in Fig. 4(b) illustrate that dissociation dominantly occurs for molecules that are initially placed over the bridge site, in agreement with previous findings.^10,11 This significant narrowing of the distribution for xy due to surface corrugation leads to a substantial reduction of reactivity by more than a factor of 100. Fig. 4a shows that the distribution for θ narrows similarly for the reacted molecules, but the reactivity reduction in this case is only a factor of 2–3×. In the ESI,† we show that the effect of the azimuthal orientation of the molecular axis on the reactivity is even smaller (see Fig. S7, ESI†). There we also demonstrate that xy-related reactivity reduction strongly depends on the incidence energy and ranges from a 50× at low incidence energy to a 10× reduction at high incidence energy.

Altogether, assuming that the aforementioned angular and translational degrees of freedom are independent, the total reduction in reactivity due to dynamical effects is a factor of 20–250× compared to TST, depending on the incidence energy (see Section S5 in the ESI†). This accounts for more than half of the 10³–10⁴× reduction in the reaction probability shown in Fig. 2. Recrossing events in our MD simulations, i.e., molecules which do not dissociate despite crossing the barrier, could be responsible for a further reduction of the reaction probability compared to TST. Our recrossing analysis reveals that this is not the case (see Section S6 in the ESI†). Consequently, the bobsleigh effect is very likely responsible for the remaining reduction in reactivity due to dynamical effects, since it is closely related to the movement along the remaining molecular degree of freedom, the z-coordinate.^40–42 In principle, MEP deviations and barrier recrossing could be accounted for in TST by including entropic contributions and transmission factors in the reaction rate coefficients for the vibrational ground. However, accounting for these effects requires additional computational efforts, which is why we have not further investigated this here.

4 Discussion: implications and recommendations for thermal catalysis and plasma catalysis

For thermal catalysis, vibrational excitation of reactants is generally not considered to significantly affect the results. However, in this work, we predict that highly activated thermal catalysis (e.g., ammonia synthesis) is indeed significantly affected by vibrationally excited states, even if their population is comparatively low. Dynamical effects, such as the bobsleigh effect, considerably limit the effectiveness of translational energy, whereas vibrational energy is not affected as much. Therefore, for a chemical reaction where the TS is characterised by a large barrier height and a very elongated bond, we suspect that simulating only the vibrational ground state is insufficient. Interestingly, dynamical effects can also influence other aspects of the reaction. For instance, moving surface atoms can alter reactivity locally.^73,74

For plasma catalysis, we found that the effects of vibrational excitation are much more pronounced than previously anticipated based on the basis of the FM α model. Some of us have recently observed that, in general, the vibrational efficacies for DC reactions scale with the absolute barrier height.³⁵ In comparison to other metal surfaces, Ru(0001) has a relatively low barrier height for the DC of N₂. Other metals thus likely exhibit even higher vibrational efficacies than described by the FM α model. Consequently, the predicted optimal catalyst and operating conditions might change considerably, depending on the level of theory employed to account for the vibrational enhancement of DC reactions.

Similarly, it is often unclear which kind of plasma-induced effects are dominant.²⁰ Using different plasma catalytic models that try to disentangle the various effects might lead to qualitatively different conclusions. For example, the Langmuir–Rideal reaction mechanism (also, incorrectly, known as Eley–Rideal⁷⁵), where a reaction takes place directly between a gas phase reactant (in plasma catalysis typically a radical) and a surface adsorbate, is often studied along vibrationally excited reactants. For the plasma-catalytic non-oxidative coupling of methane, it has been predicted that the importance of radicals and vibrational excitation of reactants depends on the binding strength of the catalyst surface.^30,76,77 Similarly, for CO₂ hydrogenation, radicals were predicted to be much more important.³¹ In a kinetic investigation of ammonia synthesis on Fe, adsorption of radicals, Langmuir–Rideal reactions, and vibrational excitation were all found to be important.⁷³ Bayer et al.,³³ on the other hand, concluded on the basis of a combination of simulations and experiments that the loss of vibrationally excited N₂ on Fe catalysts was due to vibrational relaxation (i.e., non-reactive energy loss), and not DC. Furthermore, for dry reforming of methane on Ni, it was predicted that the Langmuir–Rideal reaction mechanism dominates, compared to vibrational excitation.³² Generally, it has been concluded that Langmuir–Rideal reactions and adsorption of radicals are much more important than vibrational excitation of reactants. However, in all cases, the FM α approach was employed. Based on our work, it can be expected that vibrational excitation plays a significantly more important role than has previously been accounted for by catalytic models. Furthermore, it is likely that the competing Langmuir–Rideal reaction mechanism occurs less often than typically modelled, especially when dynamical effects are considered.⁷⁸ This could have positive implications for the practical application of plasma catalysis. As discussed above, it is often believed that radicals are the driving force behind the enhanced chemical reactivity. However, designing a catalyst and a reactor in such a way that radicals are not ‘quenched’, which reduces the availability of radicals for desired reaction steps, is not facile. But if vibrational excitation is much more important than thus far thought, practical plasma catalysis might be more easily achieved. In short, the way dynamical effects are modelled in plasma catalysis can alter predictions considerably.

In this work, we have also observed that scaling relationships can lead to a considerable underestimation of the barrier height, necessitating explicit TS search calculations instead. Similarly, it can be expected that, for the majority of catalytically important reactions, standard DFT calculations will underestimate the barrier height, even when the TS is directly computed instead of employing the SR.⁷⁹ Self-interaction errors are likely the main culprit and are most prominent when the charge transfer at the TS is large.^80,81 Future improvements in electronic structure calculations are crucial for obtaining more reliable PESs, including barrier heights for MKMs. Fortunately, for N₂ + Ru(0001), the employed PES yields very good agreement with currently available experimental data from well-defined molecular beam experiments.^11,46

We also acknowledge that the employed QCT method may have certain limitations, such as neglecting nuclear quantum effects. Although tunnelling is likely insignificant for a molecule as heavy as N₂, zero-point energy violation does occur in our simulations (see Fig. S5 in the ESI†), leading to a slight overestimation of the reaction rate of lower vibrational states. In the future, by substituting the QCT approach with ring polymer MD, nuclear quantum effects can be simulated approximately, while keeping the simulations tractable.⁸² For the DC of methane on a moving Pt(111) surface, this approach has been demonstrated to accurately reproduce results from experiments, even for low incidence energies and high vibrational temperatures.⁸³ However, to the best of our knowledge, ring polymer MD simulations cannot yet be performed for specific vibrational states, necessitating methodological advancements to simulate vibrational state populations found in plasma catalysis.

Despite all methodological challenges, dynamical effects remain a crucial important aspect in determining reaction rates of molecules on metal surfaces. For instance, in the case of heterogeneous catalytic cracking of ammonia on Fe surfaces, which is the reverse of the Haber–Bosch process investigated in this work, dynamical effects play a significant role not only in the diffusion coefficient of ammonia but also in reaction steps involving surface adsorbates.⁸⁴ Therefore, we conclude that dynamical effects cannot be ignored in determining the reaction rates of molecules on metal surfaces. In the context of plasma catalysis, it appears that dynamical effects are even more important than in thermal catalysis, which unfortunately cannot be fully captured by TST models (at present). Accurate MD calculations for gas-surface dynamics are challenging to perform, as evidenced by the discrepancies between the present work and the work by Kedalo et al.³⁴ (see Section 3.3 and Section S2.2 in the ESI†). In the meantime, if available in the literature, vibrational state-resolved reaction probabilities measured or calculated for molecular beams can be utilised to guide catalytic modelling. Unfortunately, accurate data are scarce, particularly for higher vibrational states and very low reaction probabilities.

Accurate determination of the vibrational enhancement of DC reactions, as exemplified in this work, requires the computation of numerous reaction probability curves, each requiring many trajectories. Therefore, the remaining discrepancy between the most advanced TST model and the MD simulations for N₂ dissociation on Ru(0001) prompts further research to develop computationally efficient alternatives to account for the impact of vibrational excitation and dynamical effects on reaction rate coefficients. Although the mean VE does improve the FM α model by reproducing at least the initial trend of the MD results, there is room for improvement as the mean VE does overestimate the effect of vibrational excitations. The sudden vector projection model⁸⁵ is an alternative to performing full dimensional MD simulations, coming at a much lower computational cost. Unfortunately, this model still struggles with the same assumption as the FM α model that vibrational energy cannot be more effective than translational energy in promoting reactivity, and that dynamical effects such as the bobsleigh effect are not explicitly accounted for. Nevertheless, some of us have recently shown that a combination of the absolute barrier height, equilibrium and transition state bond lengths, and the sudden vector projection model manages to reproduce VEs of various molecule–metal surface reactions reasonably well at a similar computational cost as the FM α model.³⁵ Furthermore, it was observed that the curvature features of the MEP could significantly enhance the prediction of VEs with simplified static models.

5 Conclusions

In this study, we explored the impact of vibrational excitation on reactants in heterogeneous and plasma catalysis, and how this effect influences predictions obtained with microkinetic models. Our work marks a significant milestone in scrutinising the influence of vibrational effects on reaction rate coefficients for strongly activated dissociative chemisorption, which can be strongly driven by dynamical effects. Specifically, we computed vibrational-state-specific reaction rate coefficients for the dissociative chemisorption of highly vibrationally excited N₂ on Ru(0001), which is the rate-limiting step in ammonia synthesis. These reaction rate coefficients were obtained using molecular dynamics simulations at a level that is the state of the art in the field of gas-surface dynamics. They were subsequently employed in a MKM for ammonia synthesis on Ru(0001). The MD simulations revealed orders of magnitude lower reaction rate coefficients for N₂ and concomitant TOFs for ammonia compared to the reaction rate coefficients from transition-state-theory, which are the current workhorse for MKMs. Going beyond scaling relationships in the practical implementation of TST brings the results closer to the MD-based reference results, which strongly encourages investment of the additional computational effort to calculate transition states and their corresponding energy barriers, as this is still much easier than performing MD simulations.

Even under typical thermal catalytic conditions, the small fraction of vibrationally excited N₂ molecules dominate the rate coefficient for N₂ dissociation due to the very high vibrational efficacy. Consequently, typical MKMs for thermal catalysis, which only consider the vibrational ground state, prove insufficient for a highly activated chemical process like ammonia synthesis. Most importantly, we predict that the reactivity gap between thermal catalysis and plasma catalysis due to vibrational excitation is substantially larger than previously anticipated. MKMs used in plasma catalysis do commonly account for vibrational excitation with the FM α approach, which unfortunately leads to a large underestimation of the reactivity of vibrationally excited N₂. Employing vibrational efficacies determined from MD simulations yields a much larger impact of vibrational excitation on the reaction rate and improves the agreement between the TST model and the MD approach. Complex features from the PES, such as corrugation and anisotropy, along with dynamical effects like the bobsleigh effect, are responsible for the remaining discrepancy between the two approaches. It remains to be seen if a simple and computationally efficient modification of the Arrhenius equation can account for these effects. Our results emphasise the significance of how energy is partitioned and modelled in a chemical system, as it can significantly alter predictions. Since the energy distributions in plasma catalysis can be highly out of equilibrium, it is probable that dynamical effects play a substantial role not only in the DC of vibrationally excited reactants but also in other plasma-induced effects (e.g., Langmuir–Rideal reactions). For the time being, we conclude that the gold standard for determining the reaction rate coefficients for DC reactions in heterogeneous catalysis and plasma catalysis involves performing dynamical simulations.

Author contributions

FvdB: data curation, formal analysis, investigation, methodology, software, validation, visualisation, and writing – original draft preparation. NG: project administration, supervision, and writing – review and editing. JM: conceptualization, funding acquisition, resources, supervision, and writing – review and editing.

Conflicts of interest

There are no conflicts to declare.

Data availability

MD input files and trajectory outcomes including uncertainty analysis and scripts needed to reproduce all results from this article and the ESI† are available at zenodo.org (https://doi.org/10.5281/zenodo.15006139). Additional data can be obtained from the authors upon reasonable request.

Acknowledgements

We thank Annemie Bogaerts for fruitful discussions and encouragement in contributing our work to this special issue. This research is co-financed by Holland High Tech through a public–private partnership in research and development within the Dutch top sector of High-Tech Systems and Materials (HTSM).

References

1. F. Haber and R. Le Rossignol, Über Die Technische Darstellung von Ammoniak Aus Den Elementen, Z. Elektrochem. Angew. Phys. Chem., 1913, 19, 53–72,  DOI:10.1002/bbpc.19130190201.
2. P. H. Emmett and S. Brunauer, The Adsorption of Nitrogen by Iron Synthetic Ammonia Catalysts, J. Am. Chem. Soc., 1933, 55, 1738–1739,  DOI:10.1021/ja01331a507.
3. G. Ertl, M. Huber, S. B. Lee, Z. Paál and M. Weiss, Interactions of Nitrogen and Hydrogen on Iron Surfaces, Appl. Surf. Sci., 1981, 8, 373–386,  DOI:10.1016/0378-5963(81)90092-1.
4. C. T. Rettner and H. Stein, Effect of Translational Energy on the Chemisorption of N₂ on Fe(111): Activated Dissociation via a Precursor State, Phys. Rev. Lett., 1987, 59, 2768–2771,  DOI:10.1103/PhysRevLett.59.2768.
5. L. Diekhöner, H. Mortensen, A. Baurichter, E. Jensen, V. V. Petrunin and A. C. Luntz, N₂ Dissociative Adsorption on Ru(0001): The Role of Energy Loss, J. Chem. Phys., 2001, 115, 9028–9035,  DOI:10.1063/1.1413746.
6. Y. Wang and H. Guo, Post-Dissociation Dynamics of N₂ on Ru(0001): How Far Can the “Hot” N Atoms Travel?, J. Phys. Chem. C, 2023, 127, 4079–4086,  DOI:10.1021/acs.jpcc.2c08671.
7. R. van Harrevelt, K. Honkala, J. K. Nørskov and U. Manthe, The Reaction Rate for Dissociative Adsorption of N₂ on Stepped Ru(0001): Six-dimensional Quantum Calculations, J. Chem. Phys., 2005, 122, 234702,  DOI:10.1063/1.1927513.
8. C. Díaz, J. K. Vincent, G. P. Krishnamohan, R. A. Olsen, G. J. Kroes, K. Honkala and J. K. Nørskov, Multidimensional Effects on Dissociation of N₂ on Ru(0001), Phys. Rev. Lett., 2006, 96, 096102,  DOI:10.1103/PhysRevLett.96.096102.
9. C. Díaz, J. K. Vincent, G. P. Krishnamohan, R. A. Olsen, G. J. Kroes, K. Honkala and J. K. Nørskov, Reactive and Nonreactive Scattering of N₂ from Ru(0001): A Six-Dimensional Adiabatic Study, J. Chem. Phys., 2006, 125, 114706,  DOI:10.1063/1.2229197.
10. C. Díaz and R. A. Olsen, A Note on the Vibrational Efficacy in Molecule-Surface Reactions, J. Chem. Phys., 2009, 130, 094706,  DOI:10.1063/1.3080613.
11. K. Shakouri, J. Behler, J. Meyer and G.-J. Kroes, Accurate Neural Network Description of Surface Phonons in Reactive Gas–Surface Dynamics: N₂ + Ru(0001), J. Phys. Chem. Lett., 2017, 8, 2131–2136,  DOI:10.1021/acs.jpclett.7b00784.
12. K. Shakouri, J. Behler, J. Meyer and G.-J. Kroes, Analysis of Energy Dissipation Channels in a Benchmark System of Activated Dissociation: N₂ on Ru(0001), J. Phys. Chem. C, 2018, 122, 23470–23480,  DOI:10.1021/acs.jpcc.8b06729.
13. H. Shi, T. Liu, Y. Fu, X. Lu, B. Fu and D. H. Zhang, Quantum Effects in the Dissociative Chemisorption of N₂ on Fe(111): Full-Dimensional Quantum Dynamics and Quasi-Classical Trajectory Study, J. Phys. Chem. C, 2021, 125, 23105–23114,  DOI:10.1021/acs.jpcc.1c05334.
14. A. Ozaki, K.-I. Aika and H. Hori, A New Catalyst System for Ammonia Synthesis, Bull. Chem. Soc. Jpn., 1971, 44, 3216,  DOI:10.1246/bcsj.44.3216.
15. J. K. Nørskov, Electronic Factors in Catalysis, Prog. Surf. Sci., 1991, 38, 103–144,  DOI:10.1016/0079-6816(91)90007-Q.
16. C. J. H. Jacobsen, S. Dahl, B. S. Clausen, S. Bahn, A. Logadottir and J. K. Nørskov, Catalyst Design by Interpolation in the Periodic Table: Bimetallic Ammonia Synthesis Catalysts, J. Am. Chem. Soc., 2001, 123, 8404–8405,  DOI:10.1021/ja010963d.
17. N. Saadatjou, A. Jafari and S. Sahebdelfar, Ruthenium Nanocatalysts for Ammonia Synthesis: A Review, Chem. Eng. Commun., 2015, 202, 420–448,  DOI:10.1080/00986445.2014.923995.
18. Y. Wang, M. Craven, X. Yu, J. Ding, P. Bryant, J. Huang and X. Tu, Plasma-Enhanced Catalytic Synthesis of Ammonia over a Ni/Al2O3 Catalyst at Near-Room Temperature: Insights into the Importance of the Catalyst Surface on the Reaction Mechanism, ACS Catal., 2019, 9, 10780–10793,  DOI:10.1021/acscatal.9b02538.
19. E. C. Neyts, Plasma-Surface Interactions in Plasma Catalysis, Plasma Chem. Plasma Process., 2016, 36, 185–212,  DOI:10.1007/s11090-015-9662-5.
20. A. Bogaerts, X. Tu, J. C. Whitehead, G. Centi, L. Lefferts, O. Guaitella, F. Azzolina-Jury, H.-H. Kim, A. B. Murphy, W. F. Schneider, T. Nozaki, J. C. Hicks, A. Rousseau, F. Thevenet, A. Khacef and M. Carreon, The 2020 Plasma Catalysis Roadmap, J. Phys. D: Appl. Phys., 2020, 53, 443001,  DOI:10.1088/1361-6463/ab9048.
21. B. Loenders, R. Michiels and A. Bogaerts, Is a Catalyst Always Beneficial in Plasma Catalysis? Insights from the Many Physical and Chemical Interactions, J. Energy Chem., 2023, 85, 501–533,  DOI:10.1016/j.jechem.2023.06.016.
22. S. Futamura, H. Einaga, H. Kabashima and L. Y. Hwan, Synergistic Effect of Silent Discharge Plasma and Catalysts on Benzene Decomposition, Catal. Today, 2004, 89, 89–95,  DOI:10.1016/j.cattod.2003.11.014.
23. L. Wang, Y. Yi, Y. Zhao, R. Zhang, J. Zhang and H. Guo, NH3 Decomposition for H2 Generation: Effects of Cheap Metals and Supports on Plasma–Catalyst Synergy, ACS Catal., 2015, 5, 4167–4174,  DOI:10.1021/acscatal.5b00728.
24. C. T. Rettner and H. Stein, Effect of Vibrational Energy on the Dissociative Chemisorption of N₂ on Fe(111), J. Chem. Phys., 1987, 87, 770–771,  DOI:10.1063/1.453575.
25. L. Romm, G. Katz, R. Kosloff and M. Asscher, Dissociative Chemisorption of N₂ on Ru(001) Enhanced by Vibrational and Kinetic Energy: Molecular Beam Experiments and Quantum Mechanical Calculations, J. Phys. Chem. B, 1997, 101, 2213–2217,  DOI:10.1021/jp962599o.
26. M. J. Murphy, J. F. Skelly, A. Hodgson and B. Hammer, Inverted Vibrational Distributions from N₂ Recombination at Ru(001): Evidence for a Metastable Molecular Chemisorption Well, J. Chem. Phys., 1999, 110, 6954–6962,  DOI:10.1063/1.478601.
27. R. R. Smith, D. R. Killelea, D. F. DelSesto and A. L. Utz, Preference for Vibrational over Translational Energy in a Gas-Surface Reaction, Science, 2004, 304, 992–995,  DOI:10.1126/science.1096309.
28. P. Mehta, P. Barboun, F. A. Herrera, J. Kim, P. Rumbach, D. B. Go, J. C. Hicks and W. F. Schneider, Overcoming Ammonia Synthesis Scaling Relations with Plasma-Enabled Catalysis, Nat. Catal., 2018, 1, 269–275,  DOI:10.1038/s41929-018-0045-1.
29. A. Fridman, Plasma Chemistry, Cambridge University Press, Cambridge, 2008 DOI:10.1017/CBO9780511546075.
30. Y. Engelmann, K. van’t Veer, Y. Gorbanev, E. C. Neyts, W. F. Schneider and A. Bogaerts, Plasma Catalysis for Ammonia Synthesis: A Microkinetic Modeling Study on the Contributions of Eley–Rideal Reactions, ACS Sustainable Chem. Eng., 2021, 9, 13151–13163,  DOI:10.1021/acssuschemeng.1c02713.
31. R. Michiels, Y. Engelmann and A. Bogaerts, Plasma Catalysis for CO2 Hydrogenation: Unlocking New Pathways toward CH3OH, J. Phys. Chem. C, 2020, 124, 25859–25872,  DOI:10.1021/acs.jpcc.0c07632.
32. J. Sun, Q. Chen, W. Qin, H. Wu, B. Liu, S. Li and A. Bogaerts, Plasma-Catalytic Dry Reforming of CH4: Effects of Plasma-Generated Species on the Surface Chemistry, Chem. Eng. J., 2024, 498, 155847,  DOI:10.1016/j.cej.2024.155847.
33. B. N. Bayer, S. Raskar, I. V. Adamovich, P. J. Bruggeman and A. Bhan, Availability and Reactivity of N₂(v) for NH3 Synthesis by Plasma Catalysis, Plasma Sources Sci. Technol., 2023, 32, 125005,  DOI:10.1088/1361-6595/ad10f0.
34. Y. M. Kedalo, A. A. Knizhnik and B. V. Potapkin, Applicability of the Fridman–Macheret α-Model to Heterogeneous Processes in the Case of Dissociative Adsorption of N₂ on the Ru Surface, J. Phys. Chem. C, 2023, 127, 11536–11541,  DOI:10.1021/acs.jpcc.3c01967.
35. N. Gerrits and A. Bogaerts, Vibrationally Excited Molecule-Metal Surface Reactions in Heterogeneous and Plasma Catalysis: Going beyond the Fridman–Macheret α Model, EES Catal., 2025, 3, 733–742,  10.1039/D5EY00062A.
36. L. B. F. Juurlink, R. R. Smith, D. R. Killelea and A. L. Utz, Comparative Study of C-H Stretch and Bend Vibrations in Methane Activation on Ni(100) and Ni(111), Phys. Rev. Lett., 2005, 94, 208303,  DOI:10.1103/PhysRevLett.94.208303.
37. B. Jackson and S. Nave, The Dissociative Chemisorption of Methane on Ni(100): Reaction Path Description of Mode-Selective Chemistry, J. Chem. Phys., 2011, 135, 114701.
38. B. Jackson, F. Nattino and G.-J. Kroes, Dissociative Chemisorption of Methane on Metal Surfaces: Tests of Dynamical Assumptions Using Quantum Models and Ab Initio Molecular Dynamics, J. Chem. Phys., 2014, 141, 054102,  DOI:10.1063/1.4891327.
39. P. M. Hundt, M. E. van Reijzen, R. D. Beck, H. Guo and B. Jackson, Quantum-State-Resolved Reactivity of Overtone Excited CH4 on Ni(111): Comparing Experiment and Theory, J. Chem. Phys., 2017, 146, 054701,  DOI:10.1063/1.4975025.
40. N. Gerrits, K. Shakouri, J. Behler and G.-J. Kroes, Accurate Probabilities for Highly Activated Reaction of Polyatomic Molecules on Surfaces Using a High-Dimensional Neural Network Potential: CHD3 + Cu(111), J. Phys. Chem. Lett., 2019, 10, 1763–1768,  DOI:10.1021/acs.jpclett.9b00560.
41. E. A. McCullough, Jr. and R. E. Wyatt, Quantum Dynamics of the Collinear (H, H2) Reaction, J. Chem. Phys., 1969, 51, 1253–1254,  DOI:10.1063/1.1672133.
42. R. A. Marcus, On the Analytical Mechanics of Chemical Reactions. Quantum Mechanics of Linear Collisions, J. Chem. Phys., 1966, 45, 4493–4499,  DOI:10.1063/1.1727528.
43. J. C. Polanyi, Concepts in Reaction Dynamics, Acc. Chem. Res., 1972, 5, 161–168,  DOI:10.1021/ar50053a001.
44. N. Gerrits, J. Geweke, D. J. Auerbach, R. D. Beck and G.-J. Kroes, Highly Efficient Activation of HCl Dissociation on Au(111) via Rotational Preexcitation, J. Phys. Chem. Lett., 2021, 12, 7252–7260,  DOI:10.1021/acs.jpclett.1c02093.
45. G.-J. Kroes, Computational Approaches to Dissociative Chemisorption on Metals: Towards Chemical Accuracy, Phys. Chem. Chem. Phys., 2021, 23, 8962–9048,  10.1039/D1CP00044F.
46. P. Spiering, K. Shakouri, J. Behler, G.-J. Kroes and J. Meyer, Orbital-Dependent Electronic Friction Significantly Affects the Description of Reactive Scattering of N₂ from Ru(0001), J. Phys. Chem. Lett., 2019, 10, 2957–2962,  DOI:10.1021/acs.jpclett.9b00523.
47. C. E. Treanor, J. W. Rich and R. G. Rehm, Vibrational Relaxation of Anharmonic Oscillators with Exchange-Dominated Collisions, J. Chem. Phys., 2003, 48, 1798–1807,  DOI:10.1063/1.1668914.
48. T. Kozák and A. Bogaerts, Splitting of CO2 by Vibrational Excitation in Non-Equilibrium Plasmas: A Reaction Kinetics Model, Plasma Sources Sci. Technol., 2014, 23, 045004,  DOI:10.1088/0963-0252/23/4/045004.
49. A. Berthelot and A. Bogaerts, Modeling of CO2 Splitting in a Microwave Plasma: How to Improve the Conversion and Energy Efficiency, J. Phys. Chem. C, 2017, 121, 8236–8251,  DOI:10.1021/acs.jpcc.6b12840.
50. F. Abild-Pedersen, J. Greeley, F. Studt, J. Rossmeisl, T. R. Munter, P. G. Moses, E. Skúlason, T. Bligaard and J. K. Nørskov, Scaling Properties of Adsorption Energies for Hydrogen-Containing Molecules on Transition-Metal Surfaces, Phys. Rev. Lett., 2007, 99, 016105,  DOI:10.1103/PhysRevLett.99.016105.
51. T. Tchakoua, N. Gerrits, E. W. F. Smeets and G.-J. Kroes, SBH17: Benchmark Database of Barrier Heights for Dissociative Chemisorption on Transition Metal Surfaces, J. Chem. Theory Comput., 2023, 19, 245–270,  DOI:10.1021/acs.jctc.2c00824.
52. L. B. F. Juurlink, P. R. McCabe, R. R. Smith, C. L. DiCologero and A. L. Utz, Eigenstate-Resolved Studies of Gas-Surface Reactivity: CH4(v3) Dissociation on Ni(100), Phys. Rev. Lett., 1999, 83, 868–871,  DOI:10.1103/PhysRevLett.83.868.
53. M. P. Schmid, P. Maroni, R. D. Beck and T. R. Rizzo, Molecular-Beam/Surface-Science Apparatus for State-Resolved Chemisorption Studies Using Pulsed-Laser Preparation, Rev. Sci. Instrum., 2003, 74, 4110–4120,  DOI:10.1063/1.1599064.
54. K. Reuter and M. Scheffler, First-Principles Kinetic Monte Carlo Simulations for Heterogeneous Catalysis: Application to the CO Oxidation at RuO2(110), Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 045433,  DOI:10.1103/PhysRevB.73.045433.
55. J. Behler and M. Parrinello, Generalized Neural-Network Representation of High-Dimensional Potential-Energy Surfaces, Phys. Rev. Lett., 2007, 98, 146401,  DOI:10.1103/Phys-RevLett.98.146401.
56. B. Hammer, L. B. Hansen and J. K. Nørskov, Improved Adsorption Energetics within Density-Functional Theory Using Revised Perdew-Burke-Ernzerhof Functionals, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 7413–7421,  DOI:10.1103/PhysRevB.59.7413.
57. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton, LAMMPS – a Flexible Simulation Tool for Particle-Based Materials Modeling at the Atomic, Meso, and Continuum Scales, Comput. Phys. Commun., 2022, 271, 108171,  DOI:10.1016/j.cpc.2021.108171.
58. A. Singraber, J. Behler and C. Dellago, Library-Based LAMMPS Implementation of High-Dimensional Neural Network Potentials, J. Chem. Theory Comput., 2019, 15, 1827–1840,  DOI:10.1021/acs.jctc.8b00770.
59. A. Singraber, T. Morawietz, J. Behler and C. Dellago, Parallel Multistream Training of High-Dimensional Neural Network Potentials, J. Chem. Theory Comput., 2019, 15, 3075–3092,  DOI:10.1021/acs.jctc.8b01092.
60. D. Seets, M. Wheeler and C. Mullins, Kinetics and Dynamics of Nitrogen Adsorption on Ru(001): Evidence for Direct Molecular Chemisorption, Chem. Phys. Lett., 1996, 257, 280–284,  DOI:10.1016/0009-2614(96)00542-8.
61. P. G. Anderson, in Applications of Fibonacci Numbers: Proceedings of ‘The Fifth International Conference on Fibonacci Numbers and Their Applications’, The University of St. Andrews, Scotland, July 20–July 24, 1992, ed. G. E. Bergum, A. N. Philippou and A. F. Horadam, Springer Netherlands, Dordrecht, 1993, vol. 5, pp. 1–9 DOI:10.1007/978-94-011-2058-6_1.
62. J. W. Arblaster, Crystallographic Properties of Ruthenium: Assessment of Properties from Absolute Zero to 2606 K, Platinum Met. Rev., 2013, 57, 127–136,  DOI:10.1595/147106713X665030.
63. M. Karplus, R. N. Porter and R. D. Sharma, Exchange Reactions with Activation Energy. I. Simple Barrier Potential for (H, H2), J. Chem. Phys., 1965, 43, 3259–3287,  DOI:10.1063/1.1697301.
64. C. C. Marston and G. G. Balint-Kurti, The Fourier Grid Hamiltonian Method for Bound State Eigenvalues and Eigenfunctions, J. Chem. Phys., 1989, 91, 3571–3576,  DOI:10.1063/1.456888.
65. C. J. Clopper and E. S. Pearson, The Use of Confidence or Fiducial Limits Illustrated in the Case of the Binomial, Biometrika, 1934, 26, 404–413,  DOI:10.2307/2331986.
66. Z. Cao, H. Wan and Q. Zhang, Density Functional Characterization of N₂ Dissociation on the Step of Ruthenium Clusters, J. Chem. Phys., 2003, 119, 9178–9182,  DOI:10.1063/1.1615761.
67. Y. K. Kim, G. A. Morgan and J. T. Yates, Site-Specific Dissociation of N₂ on the Stepped Ru(109) Surface, Surf. Sci., 2005, 598, 14–21,  DOI:10.1016/j.susc.2005.07.043.
68. H. Yang and J. L. Whitten, Dissociative Chemisorption of CH4 on Ni(111), J. Chem. Phys., 1992, 96, 5529–5537,  DOI:10.1063/1.462690.
69. N. Gerrits, H. Chadwick and G.-J. Kroes, Dynamical Study of the Dissociative Chemisorption of CHD3 on Pd(111), J. Phys. Chem. C, 2019, 123, 24013–24023,  DOI:10.1021/acs.jpcc.9b05757.
70. T. Liu, B. Fu and D. H. Zhang, Six-Dimensional Potential Energy Surface of the Dissociative Chemisorption of HCl on Au(111) Using Neural Networks, Sci. China: Chem., 2014, 57, 147–155,  DOI:10.1007/s11426-013-5005-7.
71. R. García-Muelas, Q. Li and N. López, Density Functional Theory Comparison of Methanol Decomposition and Reverse Reactions on Metal Surfaces, ACS Catal., 2015, 5, 1027–1036,  DOI:10.1021/cs501698w.
72. J. Chen, X. Zhou, Y. Zhang and B. Jiang, Vibrational Control of Selective Bond Cleavage in Dissociative Chemisorption of Methanol on Cu(111), Nat. Commun., 2018, 9, 4039,  DOI:10.1038/s41467-018-06478-6.
73. J. Sun, Q. Chen, X. Zhao, H. Lin and W. Qin, Kinetic Investigation of Plasma Catalytic Synthesis of Ammonia: Insights into the Role of Excited States and Plasma-Enhanced Surface Chemistry, Plasma Sources Sci. Technol., 2022, 31, 094009,  DOI:10.1088/1361-6595/ac8e2c.
74. R. Moiraghi, A. Lozano, E. Peterson, A. Utz, W. Dong and H. F. Busnengo, Nonthermalized Precursor-Mediated Dissociative Chemisorption at High Catalysis Temperatures, J. Phys. Chem. Lett., 2020, 11, 2211–2218,  DOI:10.1021/acs.jpclett.0c00260.
75. R. Prins, Eley–Rideal, the Other Mechanism, Top. Catal., 2018, 61, 714–721,  DOI:10.1007/s11244-018-0948-8.
76. Y. Engelmann, P. Mehta, E. C. Neyts, W. F. Schneider and A. Bogaerts, Predicted Influence of Plasma Activation on Nonoxidative Coupling of Methane on Transition Metal Catalysts, ACS Sustainable Chem. Eng., 2020, 8, 6043–6054,  DOI:10.1021/acssuschemeng.0c00906.
77. P.-A. Maitre, M. S. Bieniek and P. N. Kechagiopoulos, Plasma-Catalysis of Nonoxidative Methane Coupling: A Dynamic Investigation of Plasma and Surface Microkinetics over Ni(111), J. Phys. Chem. C, 2022, 126, 19987–20003,  DOI:10.1021/acs.jpcc.2c03503.
78. R. Michiels, N. Gerrits, E. Neyts and A. Bogaerts, Plasma Catalysis Modeling: How Ideal Is Atomic Hydrogen for Eley–Rideal?, J. Phys. Chem. C, 2024, 128, 11196–11209,  DOI:10.1021/acs.jpcc.4c02193.
79. G.-J. Kroes and J. Meyer, Best-of-Both-Worlds Computational Approaches to Difficult-to-Model Dissociation Reactions on Metal Surfaces, Chem. Sci., 2025, 16, 480–506,  10.1039/D4SC06004K.
80. N. Gerrits, E. W. F. Smeets, S. Vuckovic, A. D. Powell, K. Doblhoff-Dier and G.-J. Kroes, Density Functional Theory for Molecule–Metal Surface Reactions: When Does the Generalized Gradient Approximation Get It Right, and What to Do If It Does Not, J. Phys. Chem. Lett., 2020, 11, 10552–10560,  DOI:10.1021/acs.jpclett.0c02452.
81. A. J. Cohen, P. Mori-Sánchez and W. Yang, Insights into Current Limitations of Density Functional Theory, Science, 2008, 321, 792–794,  DOI:10.1126/science.1158722.
82. I. R. Craig and D. E. Manolopoulos, Quantum Statistics and Classical Mechanics: Real Time Correlation Functions from Ring Polymer Molecular Dynamics, J. Chem. Phys., 2004, 121, 3368–3373,  DOI:10.1063/1.1777575.
83. N. Gerrits, B. Jackson and A. Bogaerts, Accurate Reaction Probabilities for Translational Energies on Both Sides of the Barrier of Dissociative Chemisorption on Metal Surfaces, J. Phys. Chem. Lett., 2024, 2566–2572,  DOI:10.1021/acs.jpclett.3c03408.
84. S. Perego, L. Bonati, S. Tripathi and M. Parrinello, How Dynamics Changes Ammonia Cracking on Iron Surfaces, ACS Catal., 2024, 14, 14652–14664,  DOI:10.1021/acscatal.4c01920.
85. B. Jiang and H. Guo, Relative Efficacy of Vibrational vs. Translational Excitation in Promoting Atom-Diatom Reactivity: Rigorous Examination of Polanyi's Rules and Proposition of Sudden Vector Projection (SVP) Model, J. Chem. Phys., 2013, 138, 234104,  DOI:10.1063/1.4810007.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5ey00132c

---