Catalytic resonance of ammonia synthesis by simulated dynamic ruthenium crystal strain.

Ammonia affords dense storage for renewable energy as a fungible liquid fuel, provided it can be efficiently synthesized from hydrogen and nitrogen. In this work, the catalysis of ammonia synthesis was computationally explored beyond the Sabatier limit by dynamically straining a ruthenium crystal (±4%) at the resonant frequencies (102 to 105+ Hz) of N2 surface dissociation and hydrogenation. Density functional theory calculations at different strain conditions indicated that the energies of NHx surface intermediates and transition states scale linearly, allowing the description of ammonia synthesis at a continuum of strain conditions. A microkinetic model including multiple sites and surface diffusion between step and Ru(0001) terrace sites of varying ratios for nanoparticles of differing size revealed that dynamic strain yields catalytic ammonia synthesis conversion and turnover frequency comparable to industrial reactors (400°C, 200 atm) but at lower temperature (320°C) and an order of magnitude lower pressure (20 atm).

INTRODUCTION

The emergence of affordable electrical power from distributed wind turbines and solar panels requires a complementary chemical energy storage technology (). Electrochemical conversion of water to hydrogen affords a carbon-free energy storage medium, which can then be converted to liquid ammonia as an energy-dense, fungible liquid fuel (). Currently manufactured using hydrogen from methane via steam reforming, ammonia becomes more sustainable when associated carbon emissions are captured (blue ammonia) or when hydrogen is generated via carbon-free methods, such as water electrolysis (green ammonia) (). These approaches to energy storage are best suited to small-scale modular microprocesses that can be located close to distributed renewable energy sources, providing immediate energy storage with the capability for energy-dense liquid transportation (, ). However, scaling down conventional ammonia synthesis is not economically viable primarily due to the poor performance of ammonia synthesis catalysts ().

Improving ammonia synthesis catalysts for new small-scale distributed applications requires improving the rate of reaction and single-pass reactor conversion. As a sequence of surface reactions including N2 dissociation and sequential hydrogenation, ammonia synthesis exhibits a Sabatier peak (Fig. 1A) in catalytic turnover frequency (TOF) consistent with a transition in rate-limiting surface reactions with ruthenium and iron near the peak (, ). The resulting rate limitations on either side of the Sabatier peak then derive from the reduced rate of N2 dissociation and the hydrogenation of nitrogen-containing surface species (–). Catalytic TOF closer to the Sabatier peak has been obtained by tuning the catalyst surface energy using surface promoters or bimetallic catalysts such as cobalt-molybdenum (). However, the limitations of the conventional Haber-Bosch catalytic ammonia synthesis have led to alternative ammonia synthesis technologies including electrocatalysis () and plasma catalysis ().

Fig. 1.

Catalytic synthesis of ammonia via strain.

(A) Sabatier volcano catalytic TOF at 0.31% conversion of N2 on supported metal nanoparticles with 2 or 8% step sites at 50 atm and 400°C; full details in the Supplementary Materials. (B) Two-dimensional strain changes the atom-to-atom spacing (−4, 0, and +4% spacing from left to right) across surface planes and steps. (C) Two-dimensional strain shifts the d-band center of Ru by ~0.23 eV. (D) The binding energy (BE) of NH(T) on Ru terrace sites (left) varies slightly and has a low γ (<1) relative to N(T) on a clean terrace (low surface coverage). The binding energy NH(S) on Ru step sites (right) varies less than terrace sites due to the anisotropic behavior of step sites. (E) Schematic of a typical Ru nanoparticle on alumina, designed via Wulff construction, showing step sites (orange) and various NH adsorbates. These nanoparticles typically have 2 to 8% step sites (remainder terrace sites) for particle sizes of 8.0 to 2.5 nm. The Wulff construction provides only an estimate of step sites for microkinetic calculations and is not used otherwise.

One strategy to overcome the Sabatier limit is by crystal engineering whereby one purposely creates defects on the surface to create bifunctional catalysis (, ) or preferred sites (). Another strategy for catalysts to achieve turnover frequencies beyond the Sabatier limit aims to oscillate the electronic state of the catalyst at frequencies conducive to the kinetics of surface chemistry. As recently described via simulation (, ), oscillation of the catalyst electronic state between sides of the Sabatier peak can accelerate the overall reaction rate, control the extent of reaction, and even select specific reaction pathways in a reaction network (, ). These dynamic catalysts exhibit a band of resonant-applied frequencies, yielding the highest catalytic TOF unique to each catalytic reaction. Dynamic catalytic rate acceleration has already been demonstrated with dynamic formic acid electro-oxidation () and dynamic photocatalysis for methanol decomposition ().

Another method to dynamically alter the electronic state of catalysts is physical strain. Tensile or compressive strain yields physical and electronic rearrangement, described by the d-band center (, ), leading to variable adsorbate binding energy (, ). Evaluated computationally, the binding energy of common adsorbates, such as CO or O2, can vary as much as half of an electron volt for large surface strain of 3 to 5% (, ), with unique variation in binding energy for combinations of adsorbate, surface composition, and site structure.

Achieving strains sufficiently large to affect binding energy variation requires both a perturbation method and a strategy to prevent catalyst surface rearrangement under strain. Static strain can be implemented by lattice mismatch of overlayers with catalyst supports (–), including the formation of subsurface alloys or core shell particles (, ). Dynamic strain can potentially be implemented by supporting a thin catalyst layer on temporally stressed supports. For example, deposition of a Pt overlayer on LiCoO2 exhibits a strain of ~3% as lithium ions intercalate into and out of the cobalt oxide layers with applied potential (). Another variation is deposition on a physically deformed support, such as tungsten carbide achieving >3% strain on a bendable poly(methyl methacrylate) diaphragm (). Strain is also dynamically applied via propagating surface acoustic waves, producing immense compressive and tensile strain regions when propagating through metal films (–).

The challenge with implementing plastic deformation of catalytic materials emerges at high static strain, where rearrangement via dislocation motion can relax materials obviating the catalytic benefit of strained surfaces. However, in an oscillating stress field at sufficiently high frequency (~10−6 s), the thermal fluctuations of atoms required for dislocation motion are statistically improbable, thereby allowing for substantially higher dynamic stress and strain (>5%) at short time intervals (e.g., millisecond to microsecond) (, ). Large strain at high frequencies (kilohertz to megahertz) is consistent with the requirements of resonant catalysis over this range of applied frequencies.

In this work, the ammonia synthesis reaction is computationally evaluated on ruthenium (Ru) undergoing substantial dynamic strain (±4%) to understand the impact on the catalytic TOF and extent of reaction relative to equilibrium. The energies of surface intermediates and transition states were calculated under different magnitudes of compressive and tensile surface strain on both step sites [required for N2 dissociation (, )] and Ru(0001) terrace sites using density functional theory (DFT). The microkinetic model (MKM) was then evaluated under dynamic strain with variable frequency and reactor temperature (320° to 400°C) and stoichiometric feed gas pressure (20 to 200 atm) while accounting for the effect of surface coverage and surface diffusion. Simulations indicate that dynamic ammonia synthesis at 20 kHz and 320°C can match industrial reactor performance at an order of magnitude lower pressure (20 atm).

RESULTS AND DISCUSSION

The straining of Ru (±4%) in two dimensions across the surface, depicted in Fig. 1B, alters the physical and electronic state of the catalyst. Interatomic spacing across the surface was adjusted from its relaxed state of 2.74 Å down to 2.63 Å under compressive strain and 2.84 Å under tension. Step site interatomic spacing exhibits anisotropic behavior, expanding from the relaxed condition of 2.74 to 2.63 and 2.84 Å under compression and tension, respectively, along the step edge while expanding from 2.64 to 2.63 and 2.66 Å under compression and tension perpendicular to the edge. Perpendicular to the step, edge atoms are not restrained by nearest-neighbor atoms and, therefore, relax closer to an unstrained position. The shift in relative atomic position expands and narrows the d-band due to changes in the strength of atomic interactions, altering the density of states (Fig. 1C) and shifting the d-band center by about ±0.4 eV for terrace sites and ±0.3 eV for step sites over the full range of strain. The resulting shift in the d-band energy affects adsorbate binding energy [fig. S30 and “Density of states (DOS) and d-band center” and “DFT calculations” sections in the Supplementary Materials].

These electronic variations manifest in variations of the binding energy (∆G) of the surface species involved in ammonia synthesis. As depicted in Fig. 1D, the binding energy of each adsorbate NH on a terrace site at standard temperature and pressure (taken here as a negative quantity) increases with the degree of hydrogenation and with the N(T) binding energy on a clean (low coverage) Ru terrace site. In sharp contrast to N(T), the binding energy of NH surface species varies less over the full range of strain. While NH species binding energies on step sites are similar to terrace site, the anisotropic behavior of the step edge under strain results in a correspondingly lower strain impact on the binding energies. The relative change is defined as γ < 1, where γ is the ratio of binding energies of surface adsorbates (γ = ΔΔHNH/ΔΔHN). The change upon straining is larger for adsorbates that interact closer and stronger with the surface.

The binding energy also varies with the binding site (terrace versus step sites) and surface coverage. As depicted in Fig. 1E, nanoparticles of ruthenium exhibit step sites (orange), which promote N2(S) dissociation to N(S); the fraction of step sites varies from 2 to 8% for nanoparticles of diameters between 8.0 and 2.6 nm (fig. S26B). Binding energies for adsorbates are affected by nearest-neighbor adsorbates. Mean field models reflect the effect of adsorbate-adsorbate interactions by including an average interaction between adsorbates versus the fractional monolayer coverage of a species. When binding energies are high, as those of N and NH, the lateral interaction between adsorbates is also high (table S1). As surface coverages vary with changing strain levels (fig. S16), the impact on all adsorbates is substantial (fig. S17), especially when dynamic square wave oscillations result in surface vacancies changing as much as 40% (Fig. 2B and “Dynamic simulations” section in the Supplementary Materials for additional results).

Fig. 2.

Dynamic oscillation of Ru strain for catalytic ammonia synthesis compared to a static volcano.

(A) Sabatier static strain volcano catalytic TOF at 0.31% conversion of N2 on Ru nanoparticles at step site densities indicated. (B) Instantaneous catalytic TOF and terrace and step site vacancy at various frequencies. (C) TOF under applied dynamic strain frequency (0.31% conversion at steady state and 2% step sites). All calculations are done at 320°C and 50 atm. a.u., arbitrary units.

The kinetics of ammonia synthesis were evaluated within a multisite MKM using the calculated surface energies of ammonia-forming species under varying conditions of strain, pressure, and temperature. Reactions of N2 dissociation and N* hydrogenation occur on both step and terrace sites with molecules diffusing between sites. Reactions occur within a perfectly mixed continuous flow tank reactor, whose space velocity was varied to reach a steady-state conversion of N2 of 0.31% (differential conditions). Full details are provided in the Supplementary Materials and Methods (“Reaction mechanism” and “Thermodynamics, kinetics, and coverage effects” sections in the Supplementary Materials).

The kinetics of the simulated ammonia synthesis reaction on static Ru exhibits strain volcanoes as depicted in Fig. 2A at 320°C and 50 atm. Unique from the volcanoes for various metal surfaces in Fig. 1A, the strain volcano peaks shift with the density of steps to 140.7, 138.2, and 137.1 kcal/mol [on the N(T) terrace binding energy scale] for 2, 6, and 10% step sites, respectively. Larger nanoparticles (fewer step sites) shift the strain volcano peak toward conditions of higher tensile strain [a lower N(T) binding energy].

Dynamic simulations using a symmetric square wave oscillatory strain (±4%) between both sides of the strain volcano (Fig. 2B for a Ru surface composed of 2% step sites and “Terrace and step site densities” and “MATLAB kinetic model” sections in the Supplementary Materials) exhibit complex frequency response in TOF and terrace and step site vacancies. Upon transition from static to dynamic conditions, the catalytic rate increases before achieving a periodic solution (i.e., a stable limit cycle). Under these conditions, the surface coverages are changing substantially for each applied oscillation, with terrace and step site vacancy changing as much as 40 and 10%, respectively, with even higher surface coverage changes for individual adsorbates [N(T) coverage changes as much as 90%]. At 2 kHz, the catalytic rate increases one order of magnitude. Rate acceleration begins above about 1 Hz and achieves resonance at 100 to 105 Hz and beyond (Fig. 2C).

Rate enhancement via dynamic strain is interpreted via sensitivity analysis of elementary reaction steps. The impact of seven elementary reactions on Ru step (S) and terrace (T) sites, shown in Fig. 3A, on the overall rate of ammonia synthesis was evaluated by modifying the preexponential factors (by ±2%). Under static conditions (zero strain), the rate-determining step is the scission of adsorbed N2(S) on the Ru steps (step 4). Under slow oscillation (20 Hz), the rate control shifts to NH2(T) and N(T) hydrogenation (steps 1 and 3). Higher frequencies at resonance conditions (>100 Hz) shift the rate control to N2(S) dissociation (step 4) and NH2(T) and NH2(S) hydrogenation (steps 1 and 5). This indicates that rate enhancement is associated with a shift in the rate-determining step. The shift in NH2 hydrogenation from steps to terraces also demonstrates how dynamics couples to individual reaction steps. The hydrogenation reaction, under static strain, is three orders of magnitude faster on steps than on terraces. At 20 Hz, the oscillations couple with the slower terrace reaction, increasing the ammonia TOF with little impact on the step site reaction, while the higher frequency oscillations (20 kHz) couple with the faster step site reaction increasing the TOF with negligible influence on the terrace site reaction.

Fig. 3.

Rate limitations of ammonia synthesis on dynamically strained Ru.

(A) Sensitivity analysis for the rate-relevant catalytic steps in ammonia synthesis for varying applied frequencies. S and SL indicate the upper and lower step sites, and * indicates a vacancy. (B) Dynamic flux, (C) terrace (T) surface coverage, and (D) upper step (S) surface coverage of ammonia synthesis at 350 Hz. All panels: 320°C, 50 atm, 2% step sites. a.u., arbitrary units, m.l., monolayer.

The transition in rate limitation also appears in the surface coverages at 350 Hz (Fig. 3, B to D). At this frequency, the catalytic system proceeds at each extreme strain condition (±4%) for 1.43 ms, allowing for the steps and terraces to almost transition to the new strain condition. However, the surface species never achieve steady state before reverting to the next strain condition, consistent with its overall catalytic TOF within the resonance range (102 to 105+ Hz). At +4% strain, N2(S) readily dissociates on the steps accumulating nitrogen on the step and adjacent terrace sites, N(S) and N(T). Simultaneously, hydrogen dissociates on both Ru sites, as H(T) and H(S), proportional to the hydrogenation fluxes. When reverting to negative (−4%) strain, the nitrogen and hydrogen atoms are depleted, and NH3(T) and NH2(S) accumulate on the surface. Oscillating between these two conditions produces ammonia from both step and terrace surface sites, unlike the common static finding of ammonia chemistry proceeding only on steps.

Further interpretation of dynamic rate enhancement requires analysis under static strain (Fig. 4). Sensitivity analysis identifies N2 dissociation on steps as rate-limiting at −4 and 0% strain (Fig. 4A). At +4% strain, the overall TOF is higher (Fig. 2A), and reactions 2, 3, 5, and 6 are all important; the N2 dissociation is no longer rate limiting. It follows that the higher stability of N2(S) associated with the higher adsorbate binding energy at +4% strain promotes the N2 dissociation, thereby shifting the rate-determining step further down the catalytic cycle. Individual elementary reaction rates (Fig. 4B) provide insight into the reaction fluxes contributing to the overall static strain (composite volcano curve). The adsorption, desorption, and surface diffusions are fast relative to the overall catalytic reaction (figs. S31 and S32). The N2 dissociation on steps fully defines the volcano between a strain of 3.25 and −4% (Fig. 3B), consistent with the sensitivity analysis (Fig. 4A). For +4% strain, the combination of reactions 2, 3, 5, and 6 defines the left side of the strain volcano. Thus, the strain volcano is associated with a shift in rate control between N2(S) dissociation and NH hydrogenation. The oscillatory operation discussed above overcomes the N2(S) dissociation limitation, seen under most conditions of the static strain, to enhance the reaction rate.

Fig. 4.

Rate limitations of ammonia synthesis on statically strained Ru.

(A) Sensitivity of ammonia synthesis at three strains indicated. (B) Static strain-affected TOF for potentially rate-determining reactions comprising the composite volcano curve on upper step (dashed) and terrace (line) sites. (C) Terrace and upper step site coverage (fraction of monolayer) as impacted by static strain. (D) Static strain impact on reaction barriers on upper step (dashed) and terrace (line) sites. Conditions of 320°C and 50 atm with 2% Ru step sites. S and SL indicate the upper and lower step sites, and * indicates a vacancy.

The shift between N2 dissociation as the rate-determining step to the hydrogenation is accompanied by an N(S) coverage increase and H(S) decrease (Fig. 4C) when transitioning from negative to positive strain. At large positive strain (most negative binding energy), N(S) and N(T) are the dominant surface species, with H(S) no longer blocking the active sites. These shifts in surface coverage, combined with the moderate alterations in activation energy (Fig. 4D), change the rate-determining steps (summarized in Fig. 5) and enable the dynamic catalytic rate improvement.

Fig. 5.

Reaction paths of catalytic ammonia synthesis.

(A) +4%, (B) 0%, and (C) −4% strain. Numbers represent the percent of molar flux along that arrow relative to N2 gas (black text) or H2 gas (purple text) adsorbing on the surface. The normalized sensitivity coefficient of each path is indicated in teal text. Gray species and lines are reactions and species in the mechanism with de minimis reaction flux for a given static strain.

The dynamic strain results in a periodic conversion differing from equilibrium, consistent with other dynamic ratcheted systems (, ). In Fig. 6A, the ammonia synthesis reaction was simulated in a batch reactor at 320°C and 20 atm, starting with a stoichiometric mixture of H2 and N2 gases. Initial static catalysis proceeds to 33% conversion at equilibrium, after which dynamic strain (±4%) at varying frequency (5, 10, 15, and 20 kHz) was imposed. Dynamic operation increases the conversion to 45.7, 49.8, 51.2, and 51.8%, after which the catalyst returns to static operation and the chemistry proceeds to equilibrium. At higher pressure of 50 atm in Fig. 6B, variation of strain (±4%) at lower frequencies (350 Hz and 1, 3, and 5 kHz) profoundly increases the conversion above equilibrium.

Fig. 6.

Catalyst dynamic operation away from equilibrium.

(A) Surface dynamics (±4% square wave strain) initiated at time 0 at various frequencies promotes ammonia synthesis above the equilibrium conversion of 33% at 320°C and 20 atm to a periodic solution as high as 52% at 20 kHz. (B) Same as (A) but at 50 atm; the conversion is as high as 73% at 5 kHz. (C) The static reaction enthalpy diagram of ammonia synthesis on Ru at ±4% strain with steady-state coverage effects included. (D) Comparison of static (●) and dynamic (○, 10 and 20 kHz; ◊, 3.0 kHz; □, 1.0 kHz; Δ, 5.0 kHz; +, 15 kHz; ×, 10 kHz; blue points for 20 atm and red points for 50 atm) ammonia synthesis with equilibrium (lines) versus temperature at indicated pressures.

Interpreting the shift in conversion of ammonia synthesis away from equilibrium under dynamic strain is complicated by the number of sequential elementary steps in the overall mechanism. As depicted in Fig. 6B, strain oscillation can promote ammonia synthesis or decomposition (forward or reverse chemistry), depending on the applied frequency. Prior simulations of model catalytic systems (A to B) identified simple criteria for assessing the directionality of dynamic mechanisms (i.e., the Δ parameter of a single elementary reaction) (). However, ammonia synthesis exhibits a broad range of dynamic behaviors; as depicted in Fig. 6C, the reaction enthalpy diagrams of both positive and negative ±4% strain depict intermediate surface species energies across common transition states that vary uniquely (i.e., positive and negative γ). This includes reactions like N* hydrogenation, for which negative γ dynamic behavior is observed; the surface energies of N(T) and NH(T) change in opposite directions under variable strain. No overarching resonance theory yet exists to explain the overall directionality of complex reaction mechanisms.

Our proposed dynamic strain model demonstrates that catalytic resonance theory can be applied to reaction mechanisms more complex than the original simplified mechanism and still achieve rate and conversion enhancement. While the initial simulations of catalytic resonance only accounted for overall reactions determined by single elementary reactions (, ), the ammonia synthesis mechanism accounts for multiple pathways on two sites exchanging molecules via surface diffusion and lateral interactions of adsorbates, exhibiting linear scaling relationships in binding energy and transition state energy. Despite the inherent complexity of ammonia surface chemistry, the general predictions of catalytic resonance are retained.

Implementation of dynamic strain will require optimization of the applied dynamic strain waveform beyond the single frequency square wave oscillations used here. The time scales in the ammonia synthesis steps (Fig. 6C) indicate that more complex waveforms (amplitudes and frequencies) could improve the reaction. Moreover, comparing the model with experiment will require including the transverse and longitudinal surface waves arising from different methods of surface perturbation () that introduce different types of local surface strain (one- or two-dimensional, varying local curvature). These improvements in applied dynamic strain waveform can further enhance ammonia synthesis beyond the order of magnitude rate increase and pressure decrease predicted here.

METHODS

Electronic structure calculations

DFT calculations were performed using the Vienna Ab initio Simulation Package (VASP) (, ) code with the projector augmented wave method. The generalized gradient approximation with the Perdew-Burke-Ernzerhof () functional was applied to describe exchange correlation effects. The cutoff energy was 520 eV, and the convergence of energy and forces for geometry optimizations were set to 1 × 10−5 eV and 0.01 eV/Å, respectively. For transition states searches, we used the climbing image nudged elastic band (NEB) method in the VASP transition state theory tools (). The convergence criteria of energy and forces for NEB were set to 1 × 10−5 eV and 0.05 eV/Å, respectively. The transition states were confirmed using vibrational frequencies analysis.

A four-layer Ru(0001) slab with a (3 × 3) unit cell and a four-layer Ru(0001) slab with a (3 × 6) unit cell were used, on which three of the six rows of top Ru atoms were removed to simulate the step site. A 3 × 3 × 1 Monkhorst-Pack () was used for k-point sampling for the (3 × 3) unit cells, and the k-point was set to 2 × 4 × 1 for (3 × 6) unit cells. All unit cells include a vacuum of 20 Å between slabs.

The surface strain was studied by changing the lattice constant of Ru slabs. Starting with the unstrained lattice constant of Ru bulk, compressive and expansive strain up to 4% were applied in the two principal directions parallel to the surface. For each model, the bottom two layers of Ru atoms were fixed, while the top two layers of Ru atoms with adsorbates were allowed to relax fully.

Electronic energies and vibrational data were computed for all surface species on terrace and step sites, and all transition states were located using NEB. Surface diffusions between step and terrace sites required multiple computations on adjacent binding sites away from the steps until terrace-like sites were reached. The diffusion step with the highest reaction barrier was used to model it as a single step.

Reaction mechanism and MKM

MKMs were constructed for 0% (unstrained), +4%, and −4% strain, including both Ru terrace and step sites. We chose 2:98 step–terrace site ratio based on an 8-nm Ru nanoparticle. The impact of strain on nanoparticle shape results in modest change in the step–terrace site ratio that we have omitted from this work. (See “Terrace and step site densities” section in the Supplementary Materials and figs. S21 to S23.) Scaling relationships were developed to model a continuum of strain levels. The MKM includes six adsorption/desorptions, eight elementary surface reactions on both terrace and step sites, and five surface diffusion steps between terrace and step sites. The DFT data were converted into thermodynamic and kinetic data using the Python Multiscale Thermodynamic Toolbox (). Gas species data were obtained from National Institute of Standards and Technology (), while reaction enthalpies and associated reaction barriers at each strain were used to determine Brønsted–Evans–Polanyi relationships for each surface reaction.

Haber-Bosch operating conditions (400°C, 200 atm, 4 s) with a stoichiometric feed of H2:N2 of 3:1 were compared to milder conditions at both static strain (steady-state operation at a single strain value) and dynamic strain (strain oscillated symmetrically with a square wave between +4% and −4%) at multiple frequencies in both a flow (continuous stirred tank reactor) and batch reactors (static and dynamic equilibrium studies). The simulations were carried out in a custom-made MATLAB code. This model accounts for multiple catalyst active sites, surface coverage effects, surface reactions on both step and terrace sites, and surface diffusion between catalyst sites. The results were analyzed using reaction path analysis [the Reaction Network Viewer visualization tool ()] to compute the reaction flux and the partial equilibrium index (), i.e., the distance from equilibrium, for all reactions and via sensitivity analysis to identify the rate-determining steps.