Impact of the dopant-induced ensemble structure of hetero-double atom catalysts in electrochemical NH₃ production

Abstract

Using spin-polarized density functional theory (DFT) calculations, we examined electrochemical N₂ reduction (N₂RR) toward NH₃ production on hetero-RuM (M = 3d transition metals) double atom catalysts supported on defective graphene by means of analysis on the geometric ensemble structure, the N₂RR mechanism, the decoupling of strain, dopant and configurational effects and the d-orbital resolved density of states (ORDOS) (d_z², d_xz, d_yz, d_xy, and d_x²–y²) on the hetero-double atoms. In addition, we computationally screened novel catalysts by exploring 4d, 5d and p block metals as the hetero-M metals in the RuM system. First, we found the significantly enhanced N₂RR activity of inclined pentagon M (Fe, Mn, and Sc) double atom catalysts (RuFe has the highest activity) compared to the homo-Ru₂ double atom catalyst. Our DFT calculations on the interplay of strain, dopant and configurational effects in the inclined pentagon M (Fe, Mn, and Sc) double atom catalysts predicted that (1) the dopant effect was the promoter to improve the N₂RR activity of RuSc and RuMn, (2) the tensile strain (RuSc) tended to reduce the NH₃ productivity via the N₂RR, while the effect of compressive strain (RuFe and RuMn) was insignificant, and (3) the dopant-support interaction induced a unique inclined pentagon M double atom ensemble structure, which leads to the large reduction of the N₂RR activity of the hetero-RuSc double atom but the activity increases for the hetero-RuFe and RuMn cases. Finally, our DFT calculation on the analysis of the p–d (d_z², d_xz, d_yz, d_xy, and d_x²–y²) orbital overlap identified the key d orbitals in determining the descriptor (NH₂ adsorption energy) for representing the N₂RR. That is, the orbitals (d_z², d_xz, and d_yz) having an orientation toward the z direction in the horizontal Ru₂ double atom played an important role in determining the NH₂ adsorption process, while for the inclined pentagon M double atoms (RuFe, RuSc, and RuMn), the d_xz and d_xy orbitals were found to be essential for the modification of NH₂ adsorption energy. Finally, a descriptor based DFT search additionally discovered promising hetero-RuOs and RuIr double atom catalysts. This study highlights that the dopant engineering of hetero-double atom catalysts supported on defective graphene can significantly modify the electrochemical reactivity, particularly by the dopant type and geometric ensemble structure.

---

Introduction

To reduce global warming, many studies are being actively conducted around the world to reduce anthropogenic carbon dioxide (CO₂) emissions and find new clean energy sources to replace fossil fuels. Among the candidates for new energy media, hydrogen (H₂) is being spotlighted.¹ H₂ is the most bountiful chemical substance in the universe, and it is colourless, non-toxic and highly combustible. It can be produced without carbon emission through electrolysis of water using renewable electricity and used as a power source for fuel cells that can produce electricity for ships, airplanes, automobiles, and drones.^2–4 When H₂ is consumed in a fuel cell, it only emits water as a byproduct.⁵ However, its low volumetric energy density, wide range of flammability limits, and high transportation cost are obstacles to using it directly as an energy medium.

Recently, ammonia (NH₃) has been prominent as a promising candidate for storage and transportation of H₂ due to its several advantages. It can be readily liquefied under 8 bar at ambient temperature,⁶ and its relatively high autoignition temperature (651 °C, compared to 254 °C for diesel) enables us to use it safely from fire and explosion.^7–9 Zamfirescu et al. reported that NH₃ is competitive compared to other common fuels such as gasoline, liquefied petroleum gas (LPG) and methanol due to its high volumetric and gravimetric H₂ storage density and low energy costs.⁷ It is already widely used for industrial use and has a distributional infrastructure to transport it in amounts larger than 100 million tons yearly or more.¹⁰

In industrial NH₃ production so far, the Haber–Bosch process constitutes the dominant route.^11,12 However, this process consumes a huge amount of energy due to its high operating temperatures (400–500 °C) and pressure (150–300 bar) to break the intermolecular N≡N bond, which is equivalent to about 2% of the worldwide energy use.^13,14 Moreover, natural gas (mainly methane) is used as a source of H₂, releasing massive amounts of CO₂ as a by-product.

In contrast, electrochemical NH₃ production through the N₂ reduction reaction (N₂RR) can synthesize NH₃ directly from N₂ and H₂ without carbon emission, and it can be operated by supplying renewable electricity such as solar and wind power. In general, electrolytic cells based on a solid electrolyte membrane are used for this process as they allow easy separation of the H₂ feed from the NH₃ product.^15,16 The basic equations for electrochemical NH₃ production using electrolytic cells under acidic conditions are as follows:

Anode:

Cathode:

N₂ + 6H⁺ + 6e⁻ → 2NH₃

Overall:

In the anode, water (H₂O) is oxidized through the oxygen evolution reaction, giving protons and electrons and releasing oxygen. These protons and electrons are transferred to the cathode via the membrane and electric potential, respectively. In the cathode, NH₃ is produced by the N₂RR, in which gaseous N₂ is sequentially combined with protons and electrons.

Group VIII elements (Fe, Ru and Os) have been known to have excellent activity in electrochemical NH₃ production.^17–22 In addition, the catalyst research has been expanded to noble metals (e.g. Rh, Pd, Au, Ru, Pt and their oxides, nitrides, and sulphides), non-noble metals (e.g. Ti, Fe, Ni, Mo, V, W and their oxides, nitrides and sulphides), and carbon-based non-metal catalysts.^23–39 Among those species, Ru has been reported to be the most active catalyst. For example, according to Nørskov's study, there is a volcanic relationship between the number of valence electrons and the NH₃ productivity, indicating that Ru has optimal valence electrons for NH₃ production through N₂ reduction.²⁰ Ru has exhibited a high reactivity toward electrochemical NH₃ production under room temperature and pressure conditions.^{11,15,16,21,22} However, Ru still shows poor NH₃ productivity for the commercialization of electrochemical NH₃ production.

To deal with this issue, the concept of single atom catalysts (SACs) has been recently introduced to make a breakthrough in activity for electrochemical NH₃ production, since it can improve activity beyond the existing volcanic curve, maximize the number of active sites, and achieve economic feasibility by reducing catalyst loading.^40–47 Zhao et al. theoretically reported that single Mo exhibited the lowest N₂RR onset potential (−0.35 V) among SACs embedded in defective boron nitride sheets,⁴⁰ and Lü et al. claimed that the active site for the N₂RR on an Fe catalyst anchored to an N-doped carbon framework is Fe–N–C, in which a single Fe atom is bonded to four N atoms.⁴² For Ru-based SACs, Geng et al. reported that Ru–N–C SACs showed a lower overpotential [ΔG(Ru₁–N₃) = 0.73 eV and ΔG(Ru₁–N₄) = 0.77 eV] compared to the Ru (101) catalyst [ΔG = 0.91 eV] by using calculations and experiments.⁴⁷ This study showed that Ru SACs supported on N-doped graphene have better activity than metallic Ru catalysts for the N₂RR. However, they did not discuss the origin of the enhanced N₂RR activity in detail.

In the catalytic design, the introduction of hetero-atoms into a pure metal (so-called alloy or multi-metallic catalysts) has been known as a salient approach for the enhancement of catalyst activity and durability: for example, core@shell, surface alloys, and intermetallic alloys.^48–50 Nonetheless, a study on the heteroatom addition strategy has not yet been conducted to tune the activity of a single atom catalyst.

In this study, the electrochemical activity of Ru-based SACs, homo-Ru₂ double atom catalysts (DACs), and hetero-RuM DACs having Ru and M (M = 3d transition metal) was investigated using quantum mechanics computation based on spin-polarized density functional theory (DFT). We calculate the structural stability, reaction mechanism and activity of the N₂RR on homo- and hetero-Ru DACs. In addition, the decoupling of strain, dopant and configurational effects and the d-orbital resolved density of states (ORDOS) (d_z², d_xz, d_yz, d_xy, and d_x²–y²) on the hetero-double atoms are presented for a clear understanding of enhanced catalysis. Finally, using the descriptor we have identified, we computationally screen for a novel catalyst composition by exploring more search spaces (4d, 5d, and p block metals) for boosting NH₃ production via N₂ reduction.

Computational details

All of the calculations were performed on the basis of spin-polarized density functional theory (DFT) as implemented in the Vienna Ab initio Simulation Package (VASP).^51–53 The projector augmented wave (PAW) method using a plane wave basis set was implemented to describe the interaction between core and valence electrons.⁵⁴ The Perdew–Burke–Ernzerhof (PBE) exchange-correlation functional within the generalized gradient approximation (GGA) was employed.⁵⁵ For the plane wave expansion of the electronic eigenfunctions, an energy cutoff of 400 eV was used.

To model the RuM/C catalysts, we first built a periodic geometry of monolayer graphene consisting of 60 carbon atoms from a P6₃/mmc-structured graphite supercell by removing all layers except one and stretching the z vector to 20 Å.⁵⁶ After modelling of pristine graphene, two adjacent carbon atoms were removed to create unsaturated carbon to provide embeddable sites for Ru and 3d transition metal (M) atoms.

Here, we chose defective graphene as the supporting material for DACs instead of pristine graphene. Note that the interaction between pristine graphene and DACs is not strong enough to prevent the aggregation of DACs.⁵⁷ According to Jamie H. Warner et al., an Fe double atom in defective graphene has been successfully incorporated via the defect-assisted doping method by electron beam irradiation.⁵⁸ That is, the Fe precursor (FeCl₃) solutions are added on the surface of graphene via drop-casting and in turn vacancy sites are created through irradiation with a focused electron beam using an aberration-corrected transmission electron microscope (AC-TEM), leading to an Fe dimer (Fe double atom) embedded into defective graphene.^59,60 In addition, Pt dimers can be prepared on defect-rich graphene using atomic layer deposition (ALD), through the creation of nucleation sites, single Pt atom deposition and attachment of a second Pt atom selectively on the single Pt one.⁶¹ Following these considerations and computational validation, the RuM/C catalysts were modelled by using mono- and diatomic clusters anchored to defective graphene with double vacancies (see Fig. S1†).

To find the optimized geometry and the total energy of hetero-RuM DACs, all of the atoms were fully relaxed using the conjugate gradient method until residual forces on all the constituent atoms became smaller than 5 × 10⁻² eV Å⁻¹.⁶² For Brillouin zone integration, we chose a (4 × 4 × 1) Monkhorst–Pack mesh of k-points to determine the equilibrium geometries and total energies of RuM/C catalysts.⁶³ To calculate the electronic structure, average energy and occupancy of the d-band, and atomic charge density of the catalysts, we increased the mesh size to (10 × 10 × 1).

In the total N₂RR process (N₂ + 6(H⁺ + e⁻) → 2NH₃), N₂ is reduced to NH₃ through six net coupled proton and electron (CPE) transfer steps, and each step involves the transfer of a CPE from solution to an adsorbed intermediate.³⁹ To calculate the Gibbs free energy change (ΔG) of every step for the N₂RR, we introduced a computational hydrogen electrode (CHE) model as pioneered by Nørskov.⁶⁴ In this model, the free energy of a CPE is equivalent to half that of gaseous H₂ [G(H⁺ + e⁻) = 1/2G(H₂)] under standard reaction conditions (T = 298.15 K, P = 1 bar, and pH = 0) with no external potential. According to this method, the ΔG value can be determined using the following equation;

ΔG = ΔE − TΔS + ΔZPE − neU

where ΔE, ΔS, ΔZPE, n, and U are the total energy difference directly obtained from DFT calculations, the entropy change, the change in zero-point energies, the number of electrons transferred during the reaction, and the operating electrochemical potential in the SHE, respectively. The entropies and zero-point energies of the N₂RR intermediates were computed from the vibrational frequencies, in which the adsorption vibrational mode was calculated explicitly with all atoms fixed except for N₂RR intermediates, metal atoms, and carbon atoms bonded with metal atoms.⁶⁵ The entropies of gaseous molecules under reaction conditions were taken from the NIST Chemistry WebBook.⁶⁶

Results and discussion

Mechanism of electrochemical NH₃ production on Ru single and double atom catalysts

According to previous studies,^{13,45,46,67–69} the NH₃ formation paths from the electrochemical N₂RR on a SACs can be mainly classified into two reaction routes: the dissociative and associative mechanisms, which are shown in Fig. 1. The first step for both the dissociative and associative mechanisms is the adsorption of a gaseous N₂ molecule in the catalyst [(1) N₂(g) → N₂*] (the asterisk * denotes the end-on adsorption state). In the dissociative mechanism, the adsorbed N₂ molecule (N₂*) is dissociated into two N* atoms on the catalytic surface [(2) N₂* → 2N*], which can then undergo three successive protonation reactions, leading to NH₃ production [(3) N* + (H⁺ + e⁻) → NH*, (4) NH* + (H⁺ + e⁻) → NH₂*, and (5) NH₂* + (H⁺ + e⁻) → NH₃(g)]. On the other hand, in the associative mechanism, N₂* is protonated to form N₂H* [(6) N₂* + (H⁺ + e⁻) → N₂H*], which undergoes protonation or protonolysis (the cleavage of the N–N bond by the addition of H⁺) for NH₃ formation. Here, depending on the type of catalyst, the reduction of N₂H to NH₃ follows one of the following three reaction pathways: (A) alternating path [(7) N₂H* + (H⁺ + e⁻) → HN₂H*, (8) HN₂H* + (H⁺ + e⁻) → HN₂H₂*, (9) HN₂H₂* + (H⁺ + e⁻) → H₂N₂H₂*, and (10) H₂N₂H₂* + (H⁺ + e⁻) → NH₂* + NH₃(g)], (B) distal path [(11) N₂H* + (H⁺ + e⁻) → N₂H₂* and (12) N₂H₂* + (H⁺ + e⁻) → N* + NH₃(g)] and (C) enzymatic path (unlike alternating and distal paths, N₂ is adsorbed on the catalyst surface via a side-on configuration) [(13) *N₂H* + (H⁺ + e⁻) → *HN₂H*, (14) *HN₂H* + (H⁺ + e⁻) → *HN₂H₂*, (15) *HN₂H₂* + (H⁺ + e⁻) → *H₂N₂H₂*, and (16) *H₂N₂H₂* + (H⁺ + e⁻) → NH₂* + NH₃(g) (note that double asterisks ** denote the side-on states of adsorbates)]. In addition, the N₂RR can also follow a mixed mechanism of alternating, distal and enzymatic pathways [(17) *HN₂H₂* + (H⁺ + e⁻) → NH* + NH₃(g) and (18) N₂H₂* + (H⁺ + e⁻) → HN₂H₂*].

Fig. 1 Detailed N₂RR pathway considered in this study.

Here, we considered the NH₂* + (H⁺ + e⁻) → NH₃(g) reaction as the final reaction step in the N₂RR, where the NH₃ desorption energy is indirectly included in the NH₂* protonation reaction. Note that the NH₂* + (H⁺ + e⁻) → NH₃(g) reaction is considered to be the combination of NH₂* + (H⁺ + e⁻) → NH₃* and NH₃* → NH₃(g), which makes it possible to understand the electrode potential effect on the NH₃ desorption process.⁶⁹

In this study, we first attempt to determine the preference of the dissociative mechanism in the N₂RR by calculating the free energy change for the N₂(g) → N* + N* reaction (ΔG_diss) using the following equation.

ΔG_diss = 2G(N*) − [G(N₂) + 2G(bare)]

where G(N₂), G(bare) and G(N*) indicate the free energies of the gaseous N₂ molecule, the bare catalyst without N*, and the adsorbed N*, respectively.

If the ΔG_diss value (which can be used to understand what type of N₂RR mechanism takes place more favourably) on the catalytic surface is positive (endothermic process), the dissociative mechanism may not be possible for the N₂RR to NH₃. A negative ΔG_diss value means that N* is more stable than N₂*, indicating that the N₂RR may follow the dissociative mechanism. Here, we further confirm whether the catalyst having a negative ΔG_diss value goes through the dissociative mechanism by calculating the barrier (ΔG_barr) for the N₂ dissociation reaction into 2N*.

Fig. 2 displays the most favourable N₂RR pathways and optimized geometric structures of intermediates in each reaction step on the Ru₁ SACs and homo-Ru₂ DACs. For the Ru₁ SAC case, the N₂RR is found to follow the associative ‘Distal’ mechanism (note that ΔG_diss = 2.07 eV, indicating that the dissociative mechanism is unfavourable). First, a N₂ molecule is adsorbed on top of the Ru atom in an end-on configuration [(1) N₂(g) → N₂*] and then undergoes two successive protonation reactions to form *N₂H* with a side-on configuration [(6) N₂* + (H⁺ + e⁻) → *N₂H*] and N₂H₂* [(11) *N₂H* + (H⁺ + e⁻) → N₂H₂*] with an end-on configuration. Here, the N–N bond length of N₂(g), N₂*, *N₂H*, and N₂H₂* is 1.12 Å, 1.13 Å, 1.25 Å, and 1.31 Å, respectively. Next, the adsorbed N₂H₂* dissociates into N* and NH₃(g) by protonolysis [(12) N₂H₂* + (H⁺ + e⁻) → N* + NH₃(g)]. Finally, NH₃(g) is formed by three successive protonation reactions of adsorbed N*. This is consistent with prior reports on single atom catalysis.^{40–42,45–47}

Fig. 2 The optimized geometric structures of N₂RR intermediates along the most favourable reaction pathway on (a) the Ru SAC and (b) the homo-Ru₂ DAC.

On the other hand, our DFT calculation predicts that the homo-Ru₂ DAC follows a mixed mechanism of ‘Alternating’ or ‘Distal’ and ‘Enzymatic’ pathways (note that ΔG_diss = 1.24 eV). First, a N₂ molecule is adsorbed on top of a single Ru atom in the Ru₂ DAC with an end-on configuration, which leads to side-on *N₂H* formation by protonation. Next, *N₂H* undergoes two successive protonation reactions to form *HN₂H* [(7) *N₂H* + (H⁺ + e⁻) → *HN₂H*] and *HN₂H₂* [(8) *HN₂H* + (H⁺ + e⁻) → *HN₂H₂*]. The N–N bond length in N₂*, N₂H*, HN₂H*, and HN₂H₂* is 1.14 Å, 1.28 Å, 1.37 Å, and 1.46 Å, respectively. Finally, *HN₂H* is decomposed into NH₃(g) + NH* by the protonolysis of the N–N bond [(17) *HN₂H₂* + (H⁺ + e⁻) → NH* + NH₃(g)] and then NH* is converted into NH₃(g) via two successive protonation reactions.

The calculated free energy diagrams for the N₂RR to NH₃ on the Ru SAC and homo-Ru₂ DAC are displayed in Fig. 3. The adsorption energies (ΔE_ads) are calculated as follows: [ΔE_ads(X) = EX* − (E(bare) + E(X))], where EX*, E(bare), and E(X) denote the total energy of the X-adsorbed surface, the bare surface, and the isolated X species, respectively. E(X) was calculated in a cubic vacuum unit cell of size a = 20 Å, and the calculated E(X) values are shown in Table S1.† For both the Ru₁ SAC and the homo-Ru₂ DAC, the N₂ adsorption under thermodynamic equilibrium potential conditions (0.06 V) is predicted to be an endothermic process, which comes from the large loss of entropy due to the formation of a bound state.⁶⁸ Note that ΔG = 0.54 eV and ΔE_ads(N₂) = −0.46 eV for the Ru₁ SAC, and ΔG = 0.21 eV and ΔE_ads(N₂) = −0.84 eV for the homo-Ru₂ DAC. For the N₂RR under thermodynamic equilibrium potential conditions for both the Ru₁ SAC and homo-Ru₂ DAC, the first step [(6) N₂* + (H⁺ + e⁻) → *N₂H*] is calculated to be the most endothermic [note that ΔG = 1.14 eV (Ru₁ SAC) and 0.91 eV (Ru₂ DAC)] among all the steps for the most favourable N₂RR pathway. This indicates that the first N₂ protonation [which is defined as the potential limiting step (PLS)] determines the NH₃ production potential. Consequently, for the Ru₁ SACs and homo-Ru₂ DACs, we predict an over-potential of 1.19 V and 0.96 V, respectively, suggesting the higher N₂RR activity of a homo-Ru₂ DAC than a Ru₁ SAC. This is higher than previous single Ru atom calculation results,⁴⁷ which shows that the overpotentials for Ru₁–N₃ and Ru₁–N₄ are predicted to be 0.79 V and 0.83 V, respectively. This difference is related to the presence of N atoms in the defective graphene support. In our DFT study, a single Ru atom is embedded into a defective graphene support (forming a Ru–C bond), while, in the previous calculation, single Ru atoms are associated with three N atoms (Ru₁–N₃) or four N atoms (Ru₁–N₄) in the defective graphene support (forming a Ru–N bond). Due to the higher electronegativity of the N atom (3.04) than the carbon case (2.55), the Ru in the Ru–N bond loses more electrons than that in the Ru–C case, which may lead to the increase of adsorption strength in reaction intermediates and in turn the enhancement of N₂RR activity.

Fig. 3 Calculated free energy profiles at equilibrium (○) and limiting (▲) potentials for the N₂RR on (a) the Ru SAC and (b) the homo-Ru₂ DAC via the most favourable pathway.

Next, we examine the effect of the addition of a hetero-M atom (M = Sc, Ti, V, Cr, Mn, Fe, Co, Ni) into a Ru₁ SAC in order to enhance the N₂RR activity toward NH₃ formation.

Geometric ensemble structure of hetero-RuM double atom catalysts

First, we investigate the structural stability of hetero-RuM DACs supported on defective graphene by calculating the formation energy (indicated as ΔE_f) for the addition of a hetero-M atom onto a Ru SAC [a negative (exothermic) ΔE_f value means the stable formation of hetero-RuM DACs from Ru SACs, and vice versa] as follows:

ΔE_f(RuM/C) = E(RuM/C) − [E(Ru/C) + E(M)]

where E(RuM/C), E(Ru/C), and E(M) represent the total energy of the RuM/C surface, the bare Ru₁ SAC surface, and an isolated M atom, respectively. Here, the choice of an isolated M atom for ΔE_f is based on the synthesis method of single or double atom catalysts, where atomic layer deposition or sputtering has been used for the utilization of atoms.^46,69–71

Here, to find the most stable geometric ensemble of hetero-RuM DACs, we compare the calculated ΔE_f values of nine possible ensemble geometries (referred to as c1, c2, c3, c4, c5, c6, c7, c8, and c9) (see Fig. 4 and Table 1). For the comparison, the homo-Ru₂ DAC case is also presented. We see that the most stable geometric ensemble of the RuM/C system is classified into three groups. That is, (1) group I (c4 geometry) (we call horizontal double atom); the Ru–M dimer [M = Ru (ΔE_f = −3.69 eV), Ti (ΔE_f = −4.49 eV), V (ΔE_f = −3.96 eV), and Cr (ΔE_f = −1.72 eV)] is nearly horizontally adsorbed at two carbon vacant sites, leading to chemical bond formation between a Ru (or M) atom and two C atoms and (2) group II (c8 geometry) (we call inclined pentagon M double atom); the Ru–M dimer [M = Sc (ΔE_f = −4.23 eV), Mn (ΔE_f = −1.82 eV), and Fe (ΔE_f = −2.41 eV)] are connected to one carbon vacant and hollow carbon sites. Here, Ru and M atoms are associated with rectangular- and pentagonal-shaped four carbons, respectively, leading to the formation of an inclined RuM structure. (3) Group III (c9 geometry) (we call inclined hexagon M double atom); the configuration of the Ru–M dimer [M = Co (ΔE_f = −2.49 eV) and Ni (ΔE_f = −2.87 eV)] is similar to that in the group II case, except that M atoms are bound to hexagon-shaped five carbons, respectively. Here, group III also shows an inclined RuM structure. The optimized ensemble geometries of all hetero-Ru DACs are shown in Fig. S2.† In the following section, we investigate the electrocatalytic activity of hetero-RuM DACs toward NH₃ formation in the most stable models.

Fig. 4 A schematic diagram and possible geometries for RuM/C catalysts.

Table 1 Calculated formation energy (ΔE_f/eV) of RuM according to various geometries. In this table, the bold text indicates the formation energy of the most stable geometry

	c1	c2	c3	c4	c5	c6	c7	c8	c9
Ru2/C	−2.98	−3.22	−3.39	−3.69	−2.02	−3.23	—	—	—
RuTi/C	1.09	—	—	−2.58	−1.30	−4.14	−2.74	−4.23	−4.13
RuV/C	−0.92	−2.83	—	−4.49	−2.52	−4.04	−2.90	−4.42	−4.23
RuCr/C	−1.78	−2.97	−2.55	−3.96	−1.99	−2.93	−2.16	−3.09	−3.19
RuSc/C	−0.39	−1.07	−0.79	−1.72	−0.03	−1.27	−1.04	−1.46	−1.37
RuMn/C	−0.28	−0.98	−0.94	−1.75	0.20	−1.45	−1.20	−1.82	−1.71
RuFe/C	−0.90	−1.74	−1.64	−2.26	−0.71	−1.81	−1.83	−2.41	−2.35
RuCo/C	−0.80	−1.97	−2.36	−2.31	−0.56	−2.00	−1.91	−2.40	−2.49
RuNi/C	—	−2.16	−2.00	−2.15	—	−2.48	−2.51	−2.73	−2.87

---

Reactivity of hetero-RuM double atom catalysts

Fig. 5 shows the calculated ΔG_diss for the N₂(g) → N* + N* reaction on the surface of the hetero-RuM DACs. We observe the endothermic process when M is Ti, V, Cr, and Fe in RuM/C [note that ΔG_diss = 1.29 eV (Ti), ΔG_diss = 0.85 eV (V), ΔG_diss = 0.58 eV (Cr), and ΔG_diss = 0.07 eV (Fe)], suggesting that the N₂RR via the dissociative mechanism is energetically unfavourable. On the other hand, the RuSc/C, RuMn/C, RuCo/C, and RuNi/C cases exhibit negative ΔG_diss values [−1.39 eV (Sc), −0.28 eV (Mn), −0.79 eV (Co), and −0.01 eV (Ni)], indicating that there is a possibility for the dissociative N₂RR on those catalysts. Even if the ΔG_diss value is less than zero, the dissociative NRR will not occur if the barrier (ΔG_barr) is insurmountably high under the reaction conditions. The ΔG_diss value on the stepped Ru (0001) surface, which is known to be capable of the dissociative N₂RR under standard temperature and pressure conditions, is 0.78 eV.^21,72 Zhao et al. reported that they did not consider the dissociative N₂RR because the dissociation barrier of N₂* in their catalyst (Mo SAC on BN) was as high as 4.25 eV.⁴⁰ To clearly determine which mechanism to follow on hetero-RuM/C, we calculate ΔG_barr for the RuSc/C, RuMn/C, and RuCo/C catalysts with substantial negative values of ΔG_diss, as well as for the RuNi/C and RuFe/C catalysts with near zero values of ΔG_diss (see Fig. S3†). We find that all the ΔG_barr values are above 2.2 eV [note that ΔG_barr(Sc), ΔG_barr(Mn), ΔG_barr(Fe), ΔG_barr(Co), and ΔG_barr(Ni) are 2.92 eV, 2.43 eV, 2.71 eV, 2.40 eV, and 2.22 eV, respectively], suggesting that the N₂RR via the dissociative mechanism cannot occur on the RuSc/C, RuMn/C, RuCo/C, RuNi/C and RuFe/C surfaces under the reaction conditions.

Fig. 5 Calculated N₂ dissociation energy (ΔG_diss) to determine the reaction mechanism.

Fig. 6 displays the Gibbs free energy diagram for the favourable pathway for the N₂RR on the hetero-RuM DACs. The PLS for the N₂RR on the Ru₁ SAC and homo-Ru₂ DAC is the protonation of N₂* to form N₂H* (see Fig. 3). However, in the case of hetero-RuM DACs, the PLS is changed to the protonation of NH₂* except for hetero-RuFe and RuTi DACs. We also find that the onset potential for NH₃ formation strongly depends on the type of M, and hetero-RuM DACs show a less negative onset potential than the Ru₁ SAC (−1.14 V) except for the hetero-RuV DAC. In particular, the hetero-RuSc, RuMn, and RuFe in group II (inclined pentagon M double atom) have a less negative onset potential [−0.89 V (RuSc/C), −0.81 V (RuMn/C), and −0.73 V (RuFe/C), where RuFe exhibits the best activity toward NH₃ production] than the homo-Ru₂ DAC [−0.91 V], indicating a reduced over-potential and in turn enhanced NH₃ productivity for the hetero-RuSc, RuMn, and RuFe DACs. In contrast, for the hetero-RuTi, RuV, and RuCr DACs in group I (horizontal double atom) [−1.12 V (RuTi/C), −1.27 V (RuV/C), and −1.06 V (RuCr/C)], and RuCo and RuNi DACs in group III (inclined hexagon M double atom) [−1.02 V (RuCo/C) and −1.06 V (RuNi/C)], the opposite is true.

Fig. 6 Calculated Gibbs free energy profiles at equilibrium (○) and limiting (▲) potentials for the N₂RR on hetero-RuM DACs of group I [((a) RuTi/C, (b) RuV/C, and (c) RuCr/V)], group II [(d) (RuSc/C, (e) RuMn/C, and (f) RuFe/C)], and group III [(g) RuCo/C and (h) RuNi/C] via the most favourable pathway.

Next, we display the change of the onset potentials of the Ru₁ SAC, homo-Ru DAC, and hetero-RuM DACs for NH₃ production as a function of the adsorption energy of N₂H and NH₂ (see Fig. 7). The calculated ΔE_ads(N₂H) and ΔE_ads(NH₂) are summarized in Table S2.†

Fig. 7 Relationship between the onset potential of the N₂RR and the adsorption energies of (a) N₂H and (b) NH₂. In (b), the PLS of the catalysts located on the left side of the peak position is (5) NH₂* + (H⁺ + e⁻) → NH₃(g), while for the right-side case, the PLS is (6) N₂* + (H⁺ + e⁻) → N₂H*.

We see a clear volcanic relationship between the onset potential of the N₂RR and ΔE_ads(NH₂), whereas there is no correlation between the onset potential and ΔE_ads(N₂H). The strong dependence of N₂RR activity in the RuM system on ΔE_ads(NH₂) [rather than ΔE_ads(N₂H)] is related to the PLS we have identified. For most RuM (M = V, Cr, Sc, Mn, Co, Ni) cases, the PLS is NH₂* protonation [(4) NH₂* + (H⁺ + e⁻) → NH₃(g)], where the level of ΔE_ads(NH₂) determines the reactivity toward NH₃ production. That is, the strong binding of NH₂ to the catalyst leads to a decrease in reaction potential and vice versa. Thus, ΔE_ads(NH₂) can be used for the representation of NH₃ productivity via the N₂RR, which is called as the descriptor.

Looking closely at the volcano plot of Fig. 7(b), the lowering of the descriptor values [ΔE_ads(NH₂)] compared to the Ru₁ SAC increases the activity (onset potential) of hetero-RuM DACs (M = Ru, Ti, and Fe) toward NH₃ production until the ΔE_ads(NH₂) value approaches about −3.68 eV (which is defined as the peak position of the volcanic activity-descriptor plot). After that, the N₂RR activity is reduced (onset potential is lowered) as descriptor values are further lowered. Here, the opposition correlation between the right and left sides of the volcano plot is related to the type of PLS. On the left leg of the volcano plot, the PLS is the protonation of NH₂* [(5) NH₂* + (H⁺ + e⁻) → NH₃(g)] regardless of its mechanism, where the level of NH₂ adsorption energy (descriptor value) determines the reactivity toward NH₃ production. That is, the strong binding of NH₂ to the catalyst leads to a decrease in the reaction potential. On the other hand, on the right leg of the volcano plot, the N₂RR occurs via the associative mechanism due to the weak ability to bind the reaction intermediates such as NH₂, where the increase of NH₂ affinity enhances the NH₃ productivity.

Origin of the enhanced N₂RR catalysis for hetero-RuM double atom catalysts

Decoupling strain, dopant and configurational effects. The reactivity of hetero-RuM DACs compared to the homo-Ru₂ DAC can be modified by three factors. The first one is related to the atomic size mismatch in the Ru–M system. In bimetallic catalysis, the lattice distance of the bulk metal determines the d-band structure and in turn the surface reactivity.^72–74 Thus, the atomic distance between Ru and M is one of the most important factors in determining the activity of the catalyst, which is called as the strain effect. The second one is the change of catalytic activity in both Ru and M atoms induced by the interatomic mixing between hetero-Ru and M atoms.^75,76 That is, compared to the respective pure Ru and M atoms, the d-band structures of both Ru and M atoms in hetero-RuM DACs are significantly changed by the electronic charge transfer between Ru and M atoms and in turn affects the binding energy of reaction intermediates on hetero-RuM DACs, which is referred to as the dopant effect. The last configurational effect is the modification of catalytic activity caused by the structural change of RuM DACs by the strong interaction between the hetero-RuM double atom and the graphene. That is, the configurational effect induces the geometric ensemble configuration change of the hetero-RuM DACs from the horizontal Ru₂ DAC to the inclined pentagon M double atom, which affects the binding energy of reaction intermediates in the N₂RR.

To clearly understand the relative role of strain, dopant and configurational effects in determining NH₃ productivity in the hetero-RuM DACs in group II catalysts (here, we choose group II catalysts which exhibit higher N₂RR activity compared to the group I and III cases), we decoupled these three effects as shown in Fig. 8 and Table S3.† We also displayed optimized N₂-, N₂H-, and NH₂-adsorbed geometries in Fig. S4.†

Fig. 8 Contribution of strain, dopant and configurational effects to the (a) onset potential, and adsorption energy of (b) N₂, (c) N₂H, and (d) NH₂ of the N₂RR on the hetero-RuM DACs. In (a), the bold solid lines indicate the change in the limiting step from [N₂* + (H⁺ + e⁻) → N₂H*] to [NH₂* + (H⁺ + e⁻) → NH₃(g)].

First, we attempt to only extract the strain effect from the total N₂RR catalysis by calculating the onset potential for strained homo-Ru₂ DACs having the same geometric configuration as the homo-Ru₂ DAC but having a different Ru–Ru bond distance (indicated by RuM_strained). To model RuM_strained, we use the optimized atomic position from hetero-RuM DACs with the geometry of horizontal double atom (group I). Compared to the homo-Ru₂ DAC case (2.276 Å), the Ru–Ru bond distance for the RuSc_strained, RuMn_strained, and RuFe_strained models is 2.300, 2.135 and 2.142 Å, respectively, which are under a tensile strain of +1.05% and compressive strain of −6.20% and −5.90%, respectively, (note that ‘+’ and ‘−’ signs represent the tensile and compressive strain, respectively). Our DFT calculation predicts a substantial decrease of onset potential by 0.15 V for the RuSc_strained catalyst, but a small variation for the RuFe_strained (onset potential decreases by 0.05 V) and RuMn_strained (onset potential decreases by ∼0 V) catalysts [see Fig. 8(a)], suggesting that the tensile strain tends to reduce the NH₃ productivity via the N₂RR, while the effect of compressive strain is insignificant. The large decrease of onset potential in the RuSc_strained catalyst is related to the decrease of N₂ adsorption energy (increase of N₂ affinity) [see Fig. 8(b)] from −0.84 V (homo-Ru₂ DAC) to −1.10 V and little change of N₂H adsorption energy [see Fig. 8(c)], which only raise the level of reactant energy for the potential limiting step [(6) N₂* + (H⁺ + e⁻) → N₂H*] and in turn increase the endothermicity of the reaction.

Next, for the calculation of the contribution of the dopant effect to the NH₃ formation, we prepared RuSc_doped, RuMn_doped, and RuFe_doped models by replacing a Ru atom by a Sc, Mn, and Fe atom in RuSc_strained, RuMn_strained, and RuFe_strained models, respectively. Here, the change of onset potential from the RuM_strained model to the RuM_doped model indicates the contribution of the dopant effect to total catalysis. We find different aspects for all three catalysts due to the dopant effect. For the RuSc_doped case, the onset potential is significantly increased to −0.45 V compared to the RuSc_strained model, which is related to the dramatic increase of N₂H affinity [ΔE_ads(N₂H) = −3.15 eV] and almost no change of N₂ binding energy [ΔE_ads(N₂) = −1.12 eV] compared to the RuSc_strained model [ΔE_ads(N₂H) = −2.31 eV and ΔE_ads(N₂) = −1.10 eV] [see Fig. 8(b) and (c)]. Consequently, this induces the shift of the PLS from [(6) N₂* + (H⁺ + e⁻) → N₂H*] to [(5) NH₂* + (H⁺ + e⁻) → NH₃(g)]. In the case of RuMn_doped, the onset potential is decreased to −1.13 V by the dopant effect compared to the RuMn_strained model. Unlike the RuSc_doped case, this is mainly connected to the stronger N₂ adsorption and weaker N₂H adsorption [ΔE_ads(N₂) = −0.90 eV and ΔE_ads(N₂H) = −2.10 eV]. However, the N₂RR onset potential (−0.93 V) and the adsorption energies of N₂ [ΔE_ads(N₂) = −0.90 eV] and N₂H [ΔE_ads(N₂H) = −2.26 eV] on the RuFe_doped model are changed little by the dopant effect compared to the RuFe_strained models [−0.96 V (onset potential), ΔE_ads(N₂) = −0.83 eV] and [ΔE_ads(N₂H) = −2.22 eV].

For the contribution of the configurational effect to the total N₂RR catalysis, we calculate the difference of the onset potentials between the RuM_doped models and hetero-RuM DACs. Notice that the configurational effect induces the geometric ensemble configuration change of the hetero-RuSc, RuMn, and RuFe DACs from group I (horizontal double atom) to group II (inclined pentagon M double atom), which affects the N₂RR catalysis. We find that for the hetero-RuSc DAC, the N₂RR activity (onset potential) is reduced by 0.44 V since the NH₂ adsorption energy is decreased from −3.29 eV to −3.78 eV by the configurational effect. For the RuMn and RuFe DAC cases, the onset potential is increased from −1.13 V and −0.93 V to −0.81 V and −0.73 V, respectively, which is related to the increase of N₂H affinity by 0.67 eV and 0.51 eV compared to the respective RuMn_doped and RuFe_doped models. Especially for the RuMn case, the potential limiting step is shifted to (5) NH₂* + (H⁺ + e⁻) → NH₃(g) due to the strong NH₂ adsorption. These results demonstrate that the strain, dopant and configurational effects play an important role in determining NH₃ productivity via electrochemical N₂ reduction. In particular, the dopant and configurational effects are responsible for the enhancement of NH₃ production.

Electronic structure analysis. To better understand the fundamental reason why the hetero-RuM DACs boost the NH₃ production compared to the homo-Ru₂ DAC, we calculate the d-orbital resolved density of states (ORDOS) (d_z², d_xz, d_yz, d_xy, and d_x²–y²) of the Ru and M (M = Sc, Mn, Fe) atoms and the p density of states (DOS) of the N atom in the NH₂ for the NH₂-adsorbed homo-Ru₂ and hetero-RuM DACs [see Fig. 9]. Here, the introduction of a hetero-M atom into the homo-Ru₂ atom can affect the electronic structure of the catalyst through dopant and configurational effects and in turn modify the binding energy of the descriptor NH₂ for representing the NH₃ production catalysis. Thus, it is worthwhile to comprehend the change of the electronic structure induced by the-hereto M atom for the rational design of catalysts.

Fig. 9 Average d-orbital resolved density of states (ORDOS) (d_z², d_xz, d_yz, d_xy, and d_x²–y²) of the Ru/M atoms and the p-density of states (DOS) of the N atom in NH₂ (descriptor for representing the NH₃ production catalysis) for the NH₂-adsorbed (a) homo-Ru₂ DAC, and (b) hetero-RuSc, (c) hetero-RuMn, and (d) hetero-RuFe DACs. The DOS for the p orbital of N in NH₂ is shown as the grey area.

For the homo-Ru₂ DAC case, we see strong interatomic mixing between the p orbital of N in adsorbed NH₂ and the d_z², d_xz, d_yz, d_xy, and d_x²–y² orbitals of Ru₂ DAC in the energy ranges from −7.0 eV to −6.0 eV, from −4.5 eV to −3.5 eV and from −3.0 eV to −1.5 eV. To quantitatively understand the degree of interatomic mixing between each d orbital and p orbital, we calculate the fraction of p–d orbital overlap in the total p orbital of the N atom (referred to as χ_d) (see Fig. 10). Our DFT calculation shows that the order of overlap between p and d orbitals is d_z² (χ_dz² = 0.341) > d_xz (χ_dxz = 0.289), d_yz (χ_dyz = 0.290) > d_x²–y² (χ_dx²–y² = 0.250) > d_xy (χ_dxy = 0.208). This suggests that the orbitals (d_z², d_xz, and d_yz) having an orientation toward the z direction in the Ru₂ double atom [which is horizontally placed in the x–y plane; group I] have a higher impact in determining the adsorption of NH₂ with the catalyst than d_x²–y² and d_xy.

Fig. 10 Fraction of p–d orbital overlap in the total p orbital of the N atom (indicated by χ_d) according to the variation of the type of d orbital (d_z², d_xz, d_yz, d_xy, and d_x²–y²).

On the other hand, for the hetero-RuFe DAC case, we find a strong overlap between the p orbital and d orbitals in the range from −4.5 eV to −2.5 eV (see Fig. 9(d)). In particular, unlike the homo-Ru₂ DAC case (key orbitals in bonding are in the order of d_z² > d_xz > d_yz), the d_xz and d_xy orbitals exhibit larger interatomic mixing with p orbital cases, suggesting that the orbitals (d_xz and d_xy) having an orientation toward the x direction play an important role in stabilizing a NH₂ adsorbate on the hetero-RuFe DAC. Note that d_xz (χ_dxz = 0.336) > d_xy (χ_dxy = 0.315) > d_yz (χ_dyz = 0.296) > d_z² (χ_dz² = 0.261) > d_x²–y² (χ_dx²–y² = 0.182) (see Fig. 10). Here, the different nature of p–d orbital overlap between the homo-Ru₂ DAC and hetero-RuFe DAC is related to the geometric ensemble difference of double atoms induced by the strong hetero-dopant-support effect. That is, unlike the homo-Ru₂ DAC [whose two atoms are horizontally placed in the x direction; group I (horizontal double atom)], the hetero-RuFe DAC [group II (inclined pentagon M double atom)] has an inclined geometric ensemble structure in the y-direction where each Ru and M atom has a different z position in the normal direction of the x–y plane (see Fig. S2†). This makes the d_xz and d_xy orbitals more favourable for the interaction with the p orbital of NH₂ than the d_z², d_x²–y² and d_yz orbitals, leading to the tilted adsorption configuration of NH₂ on the RuFe DAC. For the RuSc and RuMn DACs, we see a similar p–d orbital overlap to the RuFe case, that is, the p orbital mixing mainly with d_xz and d_xy in the range from −5.0 eV to −2.7 eV (from −7.0 eV to −6.0 eV) for the RuSc DAC and from −4.7 eV to −2.6 eV for RuMn DAC. This demonstrates that the dopant and configurational effects in a hetero-RuM double atom induce an inclined RuM geometric ensemble configuration and in turn activate the d_xz and d_xy orbitals for NH₂ adsorption (rather than the d_z² and d_xy + d_x²–y² orbitals), leading to the modification of NH₂ binding energy.

Fig. 11 shows the relationship between ΔE_ads(NH₂) and χ_d for the RuM (M = Sc, Mn, Fe) DACs. We find that for the d_xz and d_xy orbitals, the increase of p–d orbital overlap (χ_d) enhances the NH₂ affinity linearly, while for the d_z², d_x²–y² and d_yz orbitals, no clear relationship is observed. This indicates that the d_xz and d_xy orbitals have a higher impact in determining the NH₂ adsorption process, which is in line with the larger p–d overlap of the d_xz and d_xy orbitals than the d_z², d_x²–y² and d_yz orbitals as discussed above.

Fig. 11 Relationship between NH₂ adsorption energy [ΔE_ads(NH₂)] and orbital overlap [χ_d] for the RuM (M = Sc, Mn, Fe) DACs (a) χ_d, (b) χ_dz², (c) χ_dxz, (d) χ_dyz, (e) χ_dxy and (f) χ_dx²–y².

HER activity

In addition to the N₂RR studies on the RuM DACs, we calculate the onset potential of the hydrogen evolution reaction (HER), which is one of the reasons for lowering the faradaic efficiency of NH₃ production.^13,42,45 The calculated onset potentials of the N₂RR (ΔG_NRR) and HER (ΔG_HER) and the difference of onset potentials between the N₂RR and HER (indicated by ΔΔG = ΔG_NRR − ΔG_HER) on the Ru₁ SAC, homo-Ru₂ DAC, and hetero-DACs are shown in Table 2. Our DFT calculation shows the enhanced selectivity (ΔΔG = −0.31 eV) of the most active hetero-ReFe double atom catalyst toward the N₂RR over the HER compared to the Ru₁ SAC (ΔΔG = −0.91 eV) and homo-Ru DAC (ΔΔG = −0.54 eV) cases. However, the onset potentials of the N₂RR for all RuM catalysts are still lower than the HER onset potential.

Table 2 Difference between the onset potential of the N₂RR and HER on the Ru SAC, and homo- and hetero-RuM DACs. Units are given in V

	Ru1/C	Group I				Group II			Group III
		Ru2/C	RuTi/C	RuV/C	RuCr/C	RuSc/C	RuMn/C	RuFe/C	RuCo/C	RuNi/C
ΔGNRR	−1.14	−0.91	−1.12	−1.27	−1.06	−0.89	−0.81	−0.73	−1.02	−1.06
ΔGHER	−0.22	−0.37	−0.20	−0.77	−0.55	−0.36	−0.29	−0.42	−0.81	−0.80
ΔΔG	−0.91	−0.54	−0.92	−0.51	−0.51	−0.53	−0.52	−0.31	−0.21	−0.26

---

Screening for additional hetero-RuM double atom catalysts

In the above investigation of the hetero-RuM DACs (whose M is a 3d transition metal) for NH₃ production, we concluded that a hetero-RuFe DAC is the best catalyst. To find out if there is a catalyst with better activity than a hetero-RuFe DAC, we further screen the candidate materials by engineering the chemical composition of RuM (M = 4d, 5d, and p blocks) using the descriptor (which is NH₂ adsorption energy and is indicated by λ) to predict the onset potential. First, the most stable geometry of the expanded hetero-RuM DACs was investigated using the same procedure as that in the case of 3d transition metals, and these results are shown in Table S4.† Using this optimized hetero-RuM DAC geometry, we further calculated the adsorption energy of NH₂, the descriptor for N₂RR activity we found. Fig. 12(a) displays the relative descriptor (denoted as Δλ) compared to the descriptor value of the peak position in the volcano plot (λ = −3.68 eV) identified in the above section (see Fig. 7(b)). We find that the descriptor value depends on the chemical composition and tends to increase as M moves to the right of the periodic table in the hetero-RuM DACs. Here, based on −0.1 eV < Δλ < 0.1 eV, we select four candidates (hetero-RuTc, RuRe, RuOs, and RuIr DACs), which can be expected to enhance NH₃ production via electrochemical N₂ reduction. The calculated onset potentials for these candidates are shown in Fig. 12(b). We see a higher NH₃ production reactivity (in particular, hetero-RuOs and RuIr DACs, whose onset potentials are −0.72 V and −0.74 V, respectively) than that of the homo-Ru₂ DAC case (−0.91 V), which is close to the onset potential of the hetero-RuFe DAC (−0.73 V).

Fig. 12 (a) Screening results using NH₂ adsorption energy as a descriptor in RuM/C, in which M is extended to 4d, 5d, 5p, and 6p elements, and (b) revised volcano plot by adding four more candidates. In (b), the PLS of the catalysts located on the left leg of the peak position is [(5) NH₂* + (H⁺ + e⁻) → NH₃(g)], while for the right leg case, the PLS is [(6) N₂* + (H⁺ + e⁻) → N₂H*].

We also calculate the onset potentials of the HER (ΔG_HER) on the additional four candidates and compared with the onset potential of the N₂RR (ΔG_NRR) to evaluate N₂RR selectivity (see Table S5†). However, ΔG_HER on all those candidates is also lower than ΔG_NRR. Improving the selectivity by increasing the N₂RR onset potential over the HER remains one of the biggest challenges to be addressed in the subsequent study.

Conclusions

In this study, using spin-polarized DFT calculations, we investigated NH₃ production catalysis via the electrochemical N₂ reduction (N₂RR) on hetero-RuM (M = 3d transition metal) double atom catalysts supported on defective graphene and computationally screened novel catalysts by exploring 4d, 5d and p block metals as the hetero-M metals.

First, we investigated the geometric ensemble structure of the hetero-RuM double atom catalysts on carbon vacant graphene. Depending on the type of hetero-atom, we identified three possible ensemble configurations; (1) group I: the horizontal double atom (Ru₂, RuTi, RuV, and RuCr), (2) group II: the inclined pentagon M double atom (RuSc, RuMn, and RuFe), and (3) group III: the inclined hexagon M double atom (RuCo and RuNi).

Second, our DFT calculation predicted the significantly enhanced N₂RR activity (reduced over-potential) of the inclined pentagon M double atom catalysts (group II) (such as RuFe, RuMn, and RuSc) (RuFe has the highest activity) compared to the homo-Ru₂ double atom case. By looking at the potential limiting step, we found that the NH₂ binding energy (so-called descriptor) was suitable for representing the N₂RR onset potential on the hetero-RuM double atom catalyst. In addition, from the volcanic activity-descriptor plot, we found that the peak position for a descriptor was around about −3.68 eV. We further searched for the catalyst candidate using the descriptor to represent onset potential by varying the chemical composition of RuM (M = 4d, 5d, and p blocks), which leads to the discovery of promising hetero-RuOs and RuIr double atom catalysts, whose onset potentials were close to the onset potential of the hetero-RuFe catalyst.

Third, the study on the relative role of strain, dopant and configurational effects in the inclined pentagon M double atom catalysts (RuFe, RuMn, and RuSc) demonstrated that the strain, dopant and configurational effects play an important role in determining NH₃ productivity via electrochemical N₂ reduction. That is, (1) the tensile strain (RuSc) tended to reduce the NH₃ productivity via the N₂RR, while the effect of compressive strain (RuFe and RuMn) was insignificant, (2) the dopant effect was found to have a significant effect on N₂RR activity, that is, a beneficial effect for the RuSc and RuFe cases and a detrimental effect for the RuMn case, and (3) the dopant-support interaction induced the change of the ensemble structure from the horizontal double atom to the inclined pentagon M double atom, where the N₂RR activity on the RuSc double atoms was substantially reduced, while for the RuFe, and RuMn cases, the change in the ensemble structure resulted in the enhancement of N₂RR activity.

Finally, our DFT calculation for the analysis of the p–d (d_z², d_xz, d_yz, d_xy, and d_x²–y²) orbital overlap identified the key d orbitals in determining the affinity of the descriptor NH₂. That is, the orbitals (d_z², d_xz, and d_yz) having an orientation toward the z direction in the horizontal Ru₂ double atom played an important role in determining the NH₂ adsorption process, while for the inclined pentagon M double atom [RuFe, RuSc, and RuMn, where each Ru and M (Fe, Sc, and Mn) atom has a substantially different z position in the normal direction of the x–y plane], the d_xz and d_xy orbitals were found to be essential for the modification of NH₂ binding energy.

Our theoretical study provides the fundamental mechanism of NH₃ production catalysis on hetero-double atom catalysts and gives physical and chemical intuition (in particular, interplay of strain, dopant and configurational effects) for next generation multi-metallic atom catalysts for hydrogen storage applications.

Author contributions

The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was financially supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (No. NRF-2020R1A2C1099711 and NRF-2018M1A2A2061975). This work was also supported by the National Supercomputing Center with supercomputing resources including technical support (KSC-2021-CRE-0264).

Notes and references

1. H.-J. Neef, International overview of hydrogen and fuel cell research, Energy, 2009, 34, 327–333.
2. L. van Biert, M. Godjevac, K. Visser and P. V. Aravid, A review of fuel cell systems for maritime applications, J. Power Sources, 2016, 327, 345–364.
3. J. W. Pratt, L. E. Klebanoff, K. Munoz-Ramos, A. A. Akhil, D. B. Curgus and B. L. Schenkman, Proton exchange membrane fuel cells for electrical power generation on-board commercial airplanes, Appl. Energy, 2013, 101, 776–796.
4. F. Ustolin and R. Taccani, Fuel cells for airborne usage: energy storage comparison, Int. J. Hydrogen Energy, 2018, 43, 11853–11861.
5. S. S. Kumar and V. Himabindu, Hydrogen production by PEM water electrolysis – a review, Mater. Sci. Energy Technol., 2019, 2, 442–454.
6. F. Jiao and B. Xu, Electrochemical ammonia synthesis and ammonia fuel cells, Adv. Mater., 2019, 31, 1805173.
7. C. Zamfirescu and I. Dincer, Ammonia as a green fuel and hydrogen source for vehicular applications, Fuel Process. Technol., 2009, 90, 729–737.
8. R. Metkemeijer and P. Achard, Ammonia as a feedstock for a hydrogen fuel cell; reformer and fuel cell behaviour, J. Power Sources, 1994, 49, 271–282.
9. W. Wang, J. M. Herreros, A. Tsolakis and A. P. E. York, Ammonia as hydrogen carrier for transportation; investigation of the ammonia exhaust gas fuel reforming, Int. J. Hydrogen Energy, 2013, 38, 9907–9917.
10. C. H. Christensen, T. Johannessen, R. Z. Sørensen and J. K. Nørskov, Towards an ammonia-mediated hydrogen economy?, Catal. Today, 2006, 111, 140–144.
11. V. Kyriakou, I. Garagounis, A. Vourros, E. Vasileiou and M. Stoukides, An electrochemical Haber-Bosch process, Joule, 2020, 4, 142–158.
12. J. Hymphreys, R. Lan and S. Tao, Development and recent progress on ammonia synthesis catalysts for Habor-Bosch process, Advanced Energy and Sustainability Research, 2021, 2, 2000043.
13. J. Kong, A. Lim, C. Yoon, J. H. Jang, H. C. Ham, J. Han, S. Nam, D. Kim, Y.-E. Sung, J. Choi and H. S. Park, Electrochemical synthesis of NH₃ at low temperature and atmospheric pressure using a γ-Fe2O3 catalyst, ACS Sustainable Chem. Eng., 2017, 5, 10986–10995.
14. N. Saadatjou, A. Jafari and S. Sahebdelfar, Ruthenium nanocatalysts for ammonia synthesis: a review, Chem. Eng. Commun., 2018, 202, 420–448.
15. N. Cao and G. Zheng, Aqueous electrocatalytic N₂ reduction under ambient conditions, Nano Res., 2008, 11, 2992–3008.
16. S. Giddey, S. P. S. Badwal and A. Kulkarni, Review of electrochemical ammonia production technologies and materials, Int. J. Hydrogen Energy, 2019, 38, 14576–14594.
17. J. K. Nørskov, Electronic factors in catalysis, Prog. Surf. Sci., 1991, 38, 103–144.
18. H. Iriawan, S. Z. Andersen, X. Zhang, B. M. Comer, J. Barrio, P. Chen, A. J. Medford, I. E. L. Stephens, I. Chorkendorff and Y. Shao-Horn, Methods for nitrogen activation by reduction and oxidation, Nature Reviews Methods Primers, 2021, 1, 56.
19. L. Forni, D. Molinari, I. Rossetti and N. Pernicone, Carbon-supported promoted Ru catalyst for ammonia synthesis, Appl. Catal., A, 1999, 185, 269–275.
20. G. Ertl, Primary steps in catalytic synthesis of ammonia, J. Vac. Sci. Technol., A, 1983, 1, 1247–1253.
21. V. Kordali, G. Kyriacou and C. Lambrou, Electrochemical synthesis of ammonia at atmospheric pressure and low temperature in a solid polymer electrolyte cell, Chem. Commun., 2000, 17, 1673–1674.
22. K. Imamura and J. Kubota, Electrochemical membrane cell for NH₃ synthesis from N₂ and H₂O by electrolysis at 200 to 250 °C using a Ru catalyst, hydrogen-permeable Pd membrane and phosphate-based electrolyte, Sustainable Energy Fuels, 2018, 2, 1278–1286.
23. K. Kugler, M. Luhn, J. A. Schramm, K. Rahimi and M. Wessling, Galvanic deposition of Rh and Ru on randomly structured Ti felts for the electrochemical NH₃ synthesis, Phys. Chem. Chem. Phys., 2015, 17, 3768–3782.
24. G. Marnellos, S. Zisekas and M. Stoukides, Synthesis of ammonia at atmospheric pressure with the use of solid state proton conductors, J. Catal., 2000, 193, 80–87.
25. J. Wang, L. Yu, L. Hu, G. Chen, H. Xin and X. Feng, Ambient ammonia synthesis via palladium-catalyzed electrohydrogenation of dinitrogen at low over-potential, Nat. Commun., 2018, 9, 1795.
26. S. Li, D. Bao, M.-M. Shi, B.-R. Wulan, J.-M. Yan and Q. Jiang, Amorphizing of Au nanoparticles by CeOx–RGO hybrid support towards highly efficient electrocatalyst for N₂ reduction under ambient conditions, Adv. Mater., 2017, 29, 1700001.
27. R. Manjunatha and A. Schechter, Electrochemical synthesis of ammonia using ruthenium–platinum alloy at ambient pressure and low temperature, Electrochem. Commun., 2018, 90, 96–100.
28. Y. Abghoui and E. Skúlason, Computational predictions of catalytic activity of zincblende (110) surfaces of metal nitrides for electrochemical ammonia synthesis, J. Phys. Chem. C, 2017, 121, 6141–6151.
29. M.-M. Shi, D. Bao, B.-R. Wulan, Y.-H. Le, Y.-F. Zhang, J.-M. Yan and Q. Jiang, Au Sub-nanoclusters on TiO2 toward highly efficient and selective electrocatalyst for N₂ conversion to NH₃ at ambient conditions, Adv. Mater., 2017, 29, 1606550.
30. S. Chen, S. Perathoner, C. Ampelli, C. Mebrahtu, D. Su and G. Centi, Electrocatalytic Synthesis of Ammonia at Room Temperature and Atmospheric Pressure from Water and Nitrogen on a Carbon-Nanotube-Based Electrocatalyst, Angew. Chem., Int. Ed., 2017, 56, 2699–2703.
31. X. Zhao, F. Yin, N. Liu, G. Li, T. Fan and B. Chen, Highly efficient metal–organic-framework catalysts for electrochemical synthesis of ammonia from N₂ (air) and water at low temperature and ambient pressure, J. Mater. Sci., 2017, 52, 10175–10185.
32. K. Kim, C.-Y. Yoo, J.-N. Kim, H. C. Yoon and J.-I. Han, Electrochemical synthesis of ammonia from water and nitrogen in ethylenediamine under ambient temperature and pressure, J. Electrochem. Soc., 2016, 163, F1523–F1526.
33. D. Yang, T. Chen and Z. Wang, Electrochemical reduction of aqueous nitrogen (N₂) at a low over-potential on (110)-oriented Mo nanofilm, J. Mater. Chem. A, 2017, 5, 18967–18971.
34. L. Zhang, X. Ji, X. Ren, Y. Ma, X. Shi, Z. Tian, A. M. Asiri, L. Chen, B. Tang and X. Sun, Electrochemical ammonia synthesis via nitrogen reduction reaction on a MoS2 satalyst: theoretical and experimental studies, Adv. Mater., 2018, 30, 1800191.
35. X. Zhang, R.-M. Kong, H. Du, L. Xia and F. Qu, Highly efficient electrochemical ammonia synthesis via nitrogen reduction reactions on a VN nanowire array under ambient conditions, Chem. Commun., 2018, 54, 5323–5325.
36. X. Ma, J. Hu, M. Zheng, D. Li, H. Lv, H. He and C. Huang, N₂ reduction using single transition-metal atom supported on defective WS2 monolayer as promising catalysts: A DFT study, Appl. Surf. Sci., 2019, 489, 684–692.
37. Y. Liu, Y. Su, X. Quan, X. Fan, S. Chen, H. Yu, H. Zhao, Y. Zhang and J. Zhao, Facile ammonia synthesis from electrocatalytic N₂ reduction under ambient conditions on N-doped porous carbon, ACS Catal., 2018, 8, 1186–1191.
38. C. Lv, Y. Qian, C. Yan, Y. Ding, Y. Liu, G. Chen and G. Yu, Defect Engineering Metal-Free Polymeric Carbon Nitride Electrocatalyst for Effective Nitrogen Fixation under Ambient Conditions, Angew. Chem., Int. Ed., 2018, 57, 10246–10250.
39. S. Mukherjee, D. A. Cullen, S. Karakalos, K. Liu, H. Zhang, S. Zhao, H. Xu, K. L. More, G. Wang and G. Wu, Metal-organic framework-derived nitrogen-doped highly disordered carbon for electrochemical ammonia synthesis using N₂ and H₂O in alkaline electrolytes, Nano Energy, 2018, 48, 217–226.
40. J. Zhao and Z. Chen, Single Mo atom supported on defective boron nitride monolayer as an efficient electrocatalyst for nitrogen fixation: a computational study, J. Am. Chem. Soc., 2017, 139, 12480–12487.
41. B. Huang, N. Li, W.-J. Ong and N. Zhou, Single atom-supported MXene: how single-atomic-site catalysts tune the high activity and selectivity of electrochemical nitrogen fixation, J. Mater. Chem. A, 2019, 7, 27620–27631.
42. F. Lü, S. Zhao, R. Guo, J. He, X. Peng, H. Bao, J. Fu, L. Han, G. Qi, J. Luo, X. Tang and X. Liu, Nitrogen-coordinated single Fe sites for efficient electrocatalytic N₂ fixation in neutral media, Nano Energy, 2019, 61, 420–427.
43. B. Qiao, A. Wang, X. Yang, L. F. Allard, Z. Jiang, Y. Cui, J. Liu, J. Li and T. Zhang, Single-atom catalysis of CO oxidation using Pt1/FeOx, Nat. Chem., 2011, 3, 634–641.
44. K. Liu, Y. Lei and G. Wang, Correlation between oxygen adsorption energy and electronic structure of transition metal macrocyclic complexes, J. Chem. Phys., 2013, 139, 204306.
45. C. Liu, Q. Li, J. Zhang, Y. Jin, D. R. MacFarlane and C. Sun, Conversion of dinitrogen to ammonia on Ru atoms supported on boron sheets: a DFT study, J. Mater. Chem. A, 2019, 7, 4771–4776.
46. C. Choi, S. Back, N. Kim, J. Lim, Y. Kim and Y. Jung, Suppression of hydrogen evolution reaction in electrochemical N₂ reduction using single-atom catalysts: a computational guideline, ACS Catal., 2018, 8, 7517–7525.
47. Z. Geng, Y. Liu, X. Kong, P. Li, K. Li, Z. Liu, J. Du, M. Shu, R. Si and J. Zeng, Achieving a Record-High Yield Rate of 120.9 μgNH₃ mgcat.-1 h-1 for N₂ Electrochemical Reduction over Ru Single-Atom Catalysts, Adv. Mater., 2018, 30, 1803498.
48. R. Burch, Importance of ligand effects in metal alloy catalysts, Acc. Chem. Res., 1982, 15, 24–31.
49. B. Li, J. Wang, X. Gao, C. Qin, D. Yang, H. Lv, Q. Xiao and C. Zhang, High performance octahedral PtNi/C catalysts investigated from rotating disk electrode to membrane electrode assembly, Nano Res., 2019, 12, 281–287.
50. S. Kim, J. Jung, E. Yang, K.-Y. Lee and D. J. Moon, Hydrogen production by steam reforming of biomass-derived glycerol over Ni-based catalysts, Catal. Today, 2014, 228, 145–151.
51. P. Hohenberg and W. Kohn, Inhomogeneous electron gas, Phys. Rev., 1964, 136, B864–B871.
52. W. Kohn and J. Sham, Self-consistent equations including exchange and correlation effects, Phys. Rev., 1965, 140, A1133–A1138.
53. G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 124, 11169–11186.
54. P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
55. J. P. Perdew, K. Burke and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett., 1996, 77, 3865–3868.
56. F. P. Bundy and J. S. Kasper, Hexagonal diamond—a new form of carbon, J. Chem. Phys., 1967, 46, 3437–3446.
57. H. Sevinçli, M. Topsakal, E. Durgun and S. Ciraci, Electronic and Magnetic Properties of 3d Transition-Metal Atom Adsorbed Graphene and Graphene Nanoribbons, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 77, 195434.
58. Z. He, K. He, A. W. Robertson, A. I. Kirkland, D. Kim, J. Ihm, E. Yoon, G.-D. Lee and J. H. Warner, Atomic Structure and Dynamics of Metal Dopant Pairs in Graphene, Nano Lett., 2014, 14, 3766–3772.
59. Y. Li, H. Su, S. H. Chan and Q. Sun, CO2 Electroreduction Performance of Transition Metal Dimers supported on Graphene: A Theoretical Study, ACS Catal., 2015, 5, 6658–6664.
60. A. W. Robertson, C. S. Allen, Y. A. Wu, K. He, J. Olivier, J. Neethling, A. I. Kirkland and J. H. Warner, Spatial Control of Defect Creation in Graphene at the Nanoscale, Nat. Commun., 2012, 3, 1144.
61. H. Yan, Y. Lin, H. Wu, W. Zhang, Z. Sun, H. Cheng, W. Liu, C. Wang, J. Li, X. Huang, T. Yao, J. Yang, S. Wei and J. Lu, Bottom-up precise synthesis of stable platinum dimers on graphene, Nat. Commun., 2017, 8, 1070.
62. H. Jiang and W. Yang, Conjugate-gradient optimization method for orbital-free density functional calculations, J. Chem. Phys., 2004, 121, 2030–2036.
63. P. E. Blöchl, O. Jepsen and O. K. Anderson, Improved tetrahedron method for Brillouin-zone integrations, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 49, 16223–16233.
64. J. K. Nørskov, J. Rossmeisl, A. Logadottir, L. Lindqvist, J. R. Kitchin, T. Bligaard and H. Jónsson, Origin of the over-potential for oxygen reduction at a fuel-cell cathode, J. Phys. Chem. B, 2004, 108, 17886–17892.
65. E. Skúlason, V. Tripkovic, M. E. Björketun, S. Gudmundsdóttir, G. Karlberg, J. Rossmeisl, T. Bligaard, H. Jónsson and J. K. Nørskov, Modeling the electrochemical hydrogen oxidation and evolution reactions on the basis of density functional theory calculations, J. Phys. Chem. C, 2010, 114, 18182.
66. NIST Standard Reference Database Number 69, https://webbook.nist.gov/chemistry/.
67. Y. Abghoui, A. L. Garden, J. G. Howalt, T. Vegge and E. Skúlason, Electroreduction of N₂ to ammonia at ambient conditions on mononitrides of Zr, Nb, Cr, and V: a DFT guide for experiments, ACS Catal., 2016, 6, 635–646.
68. X. Li, Q. Li, J. Cheng, L. Liu, Q. Yan, Y. Wu, X. Zhang, Z. Wang, Q. Qiu and Y. Luo, Conversion of dinitrogen to ammonia by FeN3-embedded graphene, J. Am. Chem. Soc., 2016, 138, 8706–8709.
69. H. Niu, X. Wang, C. Shao, Z. Zhang and Y. Guo, Computational Screening Single-Atom Catalysts Supported on g-CN for N2 Reduction: High Activity and Selectivity, ACS Sustainable Chem. Eng., 2020, 8, 13749–13758.
70. X. Zhai, L. Li, X. Liu, Y. Li, J. Yang, D. Yang, J. Zhang, H. Yan and G. Ge, A DFT screening of single transition atoms supported on MoS2 as highly efficient electrocatalysts for the nitrogen reduction reaction, Nanoscale, 2020, 12, 10035–10043.
71. Z. Xu, R. Song, M. Wang, X. Zhang, G. Liu and G. Qiao, Single atom-doped arsenene as electrocatalyst for reducing nitrogen to ammonia: a DFT study, Phys. Chem. Chem. Phys., 2020, 22, 26223–26230.
72. M. Mavrikakis, B. Hammer and J. K. Nørskov, Effect of Strain on the Reactivity of Metal Surfaces, Phys. Rev. Lett., 1998, 81, 2819–2822.
73. P. Strasser, S. Koh, T. Anniyev, J. Greeley, K. More, C. Yu, Z. Liu, S. Kaya, D. Nordlund, H. Ogasawara, M. F. Toney and A. Nilsson, Lattice-strain control of the activity in dealloyed core–shell fuel cell catalysts, Nat. Chem., 2010, 2, 454–460.
74. A. Khorshidi, J. Violet, J. Hashemi and A. A. Peterson, How strain can break the scaling relations of catalysis, Nat. Catal., 2018, 1, 263–268.
75. T. Bligaard and J. K. Nørskov, Ligand effects in heterogeneous catalysis and electrochemistry, Electrochim. Acta, 2007, 52, 5512–5516.
76. A. A. Jeffery, S.-Y. Lee, J. Min, Y. Kim, S. Lee, J. H. Lee, N. Jung and S. J. Yoo, Surface engineering of Pd-based nanoparticles by gas treatment for oxygen reduction reaction, Korean Journal of Chemical Engineering, 2020, 37(8), 1360–1364.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d1ta08358a

‡ Current address: Korea Institute of Energy Technology (KENTECH), Naju-si, Jeollanam-do, 58330, Republic of Korea. E-mail: E-mail: jhan@kentech.ac.kr

---