Kinetic and mechanistic analysis of NH3 decomposition on Ru(0001), Ru(111) and Ir(111) surfaces

We investigated the catalytic NH3 decomposition on Ru and Ir metal surfaces using density functional theory. The reaction mechanisms were unraveled on both metals, considering that, on the nano-scale, Ru particles may also present an fcc structure, hence, leading to three energy profiles. We implemented thermodynamic and kinetic parameters obtained from DFT into microkinetic simulations. Batch reactor simulations suggest that hydrogen generation starts at 400 K, 425 K and 600 K on Ru(111), Ru(0001) and Ir(111) surfaces, respectively, in excellent agreement with experiments. During the reaction, the main surface species on Ru are NH, N and H, whereas on Ir(111), it is mainly NH. The rate-determining step for all surfaces is the formation of molecular nitrogen. We also performed temperature-programmed reaction simulations and inspected the desorption spectra of N2 and H2 as a function of temperature, which highlighted the importance of N coverage on the desorption rate.

. The flow diagram of thermodynamic calculation To improve the accuracy of the energy, the zero-point energies (ZPE) are used to correct the static DFT electronic energy. ZPE refers to the vibrational energies that exist even at 0 K and is calculated as Eq S1 where accounts for the vibrational modes of the species. partition function, , is used to depict the energy as a function of temperature for the intermediates on the surface or in the gas phase, and the basic thermodynamic characters such as entropy ( ), specific heat at constant pressure ( ) and enthalpy ( ) can be derived by as the following equations: Where k is Boltzmann constant, T is the temperature. The global partition function is calculated as Eq S6.
Translational, rotational, vibrational, electronic and nuclear contributions are considered. Normally, the electronic systems are in a single electronic state, and the nuclear partition functions are unity, i.e.
, equal constant 1. The vibrational partition function of a system is obtained via Eq S7, Where is a specific vibrational mode and is the number of vibrations. The vibrational partition function in the gas phase, , is also calculated using the equation above for 3 -6 and 3 -5 vibrational degrees of freedom of a non-linear and linear molecule in the gas phase, respectively, where is the number of atoms in the molecule. The 2D-translational partition function for a free molecule is derived by the Eq S8.
Where is the average area of one active site on a catalyst.
The 3D-translational partition function for a molecule is calculated by the Eq S9, where is derived by , is the pressure of the gas phase.
( , ) The Rotational partition function for a free molecule is calculated using Eq S10 and Eq S11, depending on its symmetry and linear type.

Eq S11
Where is the symmetry factor and I is the moment of inertia defined as Eq S12, Where the sum is over the atoms in the molecule, is the mass of atom and is its distance from the rotation axis. The thermodynamic properties in our system were calculated according to the above equations. 1 To evaluate the accuracy of our method, a comparison between calculated and standard thermodynamic properties of NH 3 , N 2 , H 2 was carried out. The reference properties are from the NIST database and Thermochemical Data of Pure Substances. 2,3 Table S1 contains the data we used in the thermodynamic calculations, including the gas pressure, the electronic energy, the moment of inertia and vibration frequency from VASP calculation and the mass and symmetry of molecules.  Figure S2, Figure S3 and Figure S4 show both the calculated and reference thermodynamic data between 300-1000K. According to the Gibbs energy calculated, the delta Gibbs energy of ammonia synthesis reaction was given. The delta Gibbs energy of calculation is corrected by experimental formation enthalpy of ammonia and pure gas, which is -45.94 kJ/mol and 0 kJ/mol. Since the absolute values of the delta Gibbs energy are small at 400 and 500K, the errors turned a little bit. Generally, our method based on VASP optimized results to calculate thermodynamic properties was proved accurately.     300 K, while it will increase to 0.75 eV when temperature up to 900K. Figure S7. Energy profile of ammonia decomposition on Ru(111) under 300K, 600K and 900K Figure S8 shows the energy profile of ammonia decomposition on Ir(111) at 300K, 600K and 900K. It suggests that nitrogen recombination, R7, is the rate-determining step as well, with barrier energy ( ) of 1.94 eV at studied temperature, which is close to previous research (1.75 eV for 4 Ir(110)). 29 The first (R1) and second (R3) dehydrogenation of NH 3 needs to overcome similar barriers of around 1.55 eV at 300 K, which are increasing with the temperature. The last dehydrogenation step, R5, yields N * and H * by surpassing a 0.91 eV barrier at 300-900K, which is a moderate elementary step. The recombination of H * , R9, is impaired by the rise of temperature since the barrier energy ( ) rise from 0.62 eV at 300 K to 0.66 eV at 900K.

Ⅲ. Equilibrium constant and Reaction rate constant calculations
A script reading the partition functions, thermodynamic properties was implemented to calculate the equilibrium constant and reaction rate constant of all elementary steps. Figure S9 shows the process of the script.

Figure S9. The flow diagram of reaction rate constant and equilibrium constant calculation
The barrier energy of adsorption and desorption processes are calculated by Eq S 13, Eq S 14 and Eq S 15: The adsorption energy is calculated by Eq S 16: The classical Hertz-Knudsen equation was employed to estimate the rate of adsorption, as following Eq S 17 to Eq S 18. Where is the pre-exponential factor. The sticking coefficient, Sticky, is a 0 measure of the fraction of incident molecules which adsorb upon the surface and is calculated via Eq S 20 and Eq S 21.
As for the surface reactions in the heterogeneous catalytic system, which is considered in our research, the constant rate (k) of each surface elementary step is commonly computed using the transition state theory (TST) approximation of Eyring and Evans and Polanyi, as follows: Where is the Plank constant, kB is the Boltzmann constant, is the temperature, ℎ is the pre-exponential factor, ΔG* is the reaction activation energy, and and 0 are the partition functions of reactants and transition states respectively. The translations and rotations of the adsorbed species are frustrated on the surface and therefore we considered only vibrational modes.
We have considered an active site as a hexagonal site where the reactants and products in every elementary step occupy only one site on the surface. Consequently, the coverage of free sites, , is defined by: Where the represents the coverage of the intermediates present in the reaction ( ) system.

IV. The reaction and rate equations
The adsorption and first dehydrogenation of NH 3 are highly exothermic and the desorbed hydrogen and nitrogen molecules are assumed to leave from the surface immediately, therefore, the dissociated adsorption model is applied in the adsorption and desorption process of molecules. All the elementary steps in the ammonia decomposition and their rate equations are listed below.
Where is the time-dependent ratio of the molecule and free sites. ( )

Ⅳ. Differential equations in the TPD simulation of hydrogen and nitrogen
Temperature programmed reaction model start from pre-adsorbed NH 3 , the temperature increase at a different rate from 200 to 1000 K while any gas was extracted to avoid the re-adsorption of gases. It is applied to examine the adsorption properties of N 2 and H 2 over different surfaces in our research. = 5 -6 -2 * 7 + 2 * 8 = 1 -2 + 3 -4 + 5 -6 -2 * 9 -2 * 10 *

=-( + )
The initial parameter is input into the program by the command:

V. Differential equations in the reactor simulation
A batch reactor model under a variety of conditions is employed to investigate the catalytic properties when the metallic surface is in contact with a given pressure of NH 3 .