Computational Insight into Defective Boron Nitride Supported Double-Atom Catalysts for Electrochemical Nitrogen Reduction

: Designing highly selective and efﬁcient double-atom electrocatalysts (DACs) is essential for achieving a superior nitrogen-reduction reaction (NRR) performance. Herein, we explored the defective boron nitride–supported cage-like double-atom catalysts to rummage the qualiﬁed NRR catalysts. Based on a systematic evaluation of the stability, N 2 adsorption, NRR selectivity and activity of 10 DACs of TM 1 -TM 2 @V B -BN, we predicted Ru-Ti@V B -BN to be the NRR candidate with a limiting potential of − 0.40 V. Compared to the corresponding single-atom catalysts, the introduction of Ti/Mo modulates the d -band center of the active metal atom, which improves the NRR performance. Moreover, the magnetic Ru-Ti dimer can facilitate the transfer of charge to molecular N 2 , ensuring a signiﬁcant activation of the inert N ≡ N bond. This research not only opens up new avenues for designing boron nitride–supported DACs for NRR, but also deepens the understanding of DACs in N 2 activation.


Introduction
Ammonia (NH 3 ) plays an extremely important role in the modern agriculture and chemical industry [1,2].Currently, industrial NH 3 synthesis requires the reduction of nitrogen (N 2 ) to NH 3 with the help of iron-based catalysts at a high temperature (300-500 • C) and high pressure (150-300 atm) by means of the Haber-Bosch process [3,4].Inspired by the biological nitrogen fixation occurring under mild conditions, the electrochemical nitrogen reduction reaction (NRR) has attracted great attention as a promising NH 3 production method [5,6].However, efficient N 2 reduction has been proven to be extremely tough owing to the chemically inert N≡N triple bond with a high bond energy of 940.95 kJ mol −1 [7,8].Meanwhile, most catalysts exhibit poor selectivity toward NRR, mainly due to the existence of a competing hydrogen evolution reaction (HER) in aqueous solution and a low overpotential of HER [9][10][11][12].To overcome these challenges, searching for high-performance catalysts with low working overpotential, good stability and high selectivity remains a highly rewarding task.
As an emerging class of atomically dispersed catalysts, single-atom catalysts (SACs) have been widely used in NRR due to their high atom utilization and excellent catalytic performance [13].However, the low stability of the anchoring individual metal atom on the substrate material makes them prone to agglomerate during the reaction [14].In addition, NRR involves multiple reaction steps, and it is difficult for a single active site to simultaneously regulate the adsorbed intermediates, resulting in its low Faradaic efficiency and NH 3 yield [15].The introduction of bimetallic pairs can effectively solve this problem.For example, Chen et al. [16] systematically explored the performance of double transition metals (TMs) anchored on C 2 N (TM 2 -C 2 N) for NRR by first principles calculations.In contrast with the performance on corresponding SACs (TM 1 -C 2 N), the studies highlighted In this work, we adopted four commonly used transition metal atoms (TM = Ti, Fe, Mo and Ru) and designed cage-like diatomic TMs N2-fixing catalysts on defective BN (VB-BN).The justification for choosing these TMs lies in the fact that numerous studies have shown that Fe, Mo and Ru could act as catalytic centers for NRR [24][25][26][27].More recently, Xia et al. predicted that Ti doping could improve the stability of defective BN and enhance the NRR performance compared with the pristine BN [28]; hence, we adopted the earthabundant transition metal Ti to anchor on BN for NRR.Based on density functional theory In atomic catalysts, appropriate substrates also play major roles by stabilizing and dispersing DACs at a specific site.Boron nitride (BN) nanosheets are widely used in the field of energy storage due to their high chemical stability and unique physical and chemical properties [18][19][20].It has been proven that boron vacancies (V B ) are easier to form in h-BN sheets than nitrogen vacancies (V N ).Moreover, the h-BN sheets containing V B have been widely used as substrates to support metal atoms, which serve as a promising catalyst for the NRR process.Zhao et al. [21] showed that single Mo anchored on the defective BN exhibited high catalytic activity, along with extremely low overpotential (−0.19 V).However, the study of the DACs based on defective BN remains in its infant stage.Recently, researchers have designed DACs imitating the configuration of cage-like nitrogenase active centers, such as Fe-TM double atoms anchored on graphene [22,23].They have been reported to not only show superior NRR activity, but also effectively retard the agglomeration of metal atoms and facilitate electron transfer between N 2 and the graphene surface.It is evident that the cage-like double metal active centers have a unique application potential in NRR catalytic reactions.These pioneering works in fabricating cage-like DACs are of critical importance for the expansion of broad electrochemical reactions applying BN as supported catalysts.
In this work, we adopted four commonly used transition metal atoms (TM = Ti, Fe, Mo and Ru) and designed cage-like diatomic TMs N 2 -fixing catalysts on defective BN (V B -BN).The justification for choosing these TMs lies in the fact that numerous studies have shown that Fe, Mo and Ru could act as catalytic centers for NRR [24][25][26][27].More recently, Xia et al. predicted that Ti doping could improve the stability of defective BN and enhance the NRR performance compared with the pristine BN [28]; hence, we adopted the earth-abundant transition metal Ti to anchor on BN for NRR.Based on density functional theory (DFT), we systematically evaluated the stability of the substrate, N 2 adsorption, NRR selectivity and activity of 10 TM 1 -TM 2 @V B -BN candidates.Three promising DACs were proposed (Ru-Mo/Ru/Ti@V B -BN) as excellent NRR catalysts, and the best one was predicted to be Ru-Ti@V B -BN, with a low overpotential of −0.40 V, along the enzymatic pathway.We also performed an in-depth interaction mechanism and electronic structure analysis to reveal the inherent regularity of catalytic activity.Our work not only enriches existing BN-supported NRR electrocatalysts, but also offers a feasible clue for the experimental explorations ahead.

Results and Discussion
2.1.Structural and Electronic Properties of TM 1 -TM 2 @V B -BN Defective BN (V B -BN) with one B atom removed was adopted as the substrate to anchor the transition metal atoms (TM = Ti, Fe, Mo and Ru).Then 10 TM 1 -TM 2 @V B -BN structures were locked by embedding a cage-like TM dimer (dopant pairs) in the V B -BN substrate adjacent to the B vacancy.All the optimized structures of the TM 1 -TM 2 @V B -BN system are plotted in Figure 1a and Supplementary Figure S2.

N2 Adsorption and Selectivity on TM1-TM2@VB-BN
Firstly, the capture and adsorption of N2 are prerequisites, meaning that the Gibbs free energy change (ΔG) for stable N2 adsorption should be lower than 0 eV (ΔG(*N2)< 0)) [30].Considering the adsorption of N2, the heteronuclear and homonuclear DACs provide two adsorption sites and one adsorption site, respectively.Consequently, we totally considered 16 adsorption configurations and calculated the ΔG(*N2) values of N2 terminated on TM1-TM2@VB-BN in end-on/side-on manners (see Table 1).Moreover, the NRR catalyst should effectively inhibit competitive HER, which is also an important indicator of the catalyst's efficiency.The surface of catalysts must provide enough active sites for N2 activation, so as to facilitate the subsequent reactions.A stronger adsorption activity of the N2 molecule is essential for high NRR selectivity, which ensures that the catalyst selectively adsorbs N2 molecular prior to H [31,32].As a matter of convenience, we determined the catalytic selectivity of TM1-TM2@VB-BN through a comparison of the adsorption free energies of the *N2 and *H atom.To verify the experimental synthetic feasibility, we evaluated the structural and thermodynamic stability of all TM 1 -TM 2 @V B -BN candidates.As Supplementary Figure S2 shows, all the optimized DACs maintain their planar structures with the cage-like TM dopants protruding out of the BN monolayer, implying high structural stabilities.The formation energy of TM 1 -TM 2 @V B -BN was calculated based on the formula as follows: where E catalysts is the total energies of the TM-doped BN; µ B , µ N , µ TM1 and µ TM2 are the chemical potentials of single B, N and TM atoms in their bulk structure; n is the total number of atoms in TM 1 -TM 2 @V B -BN; and x and y are the number of B and N atoms in the substrate, respectively.The difference between the adsorption energies of the TM atoms and the cohesive energy of the corresponding base species can be used to assess whether the metal can be stably anchored on the base material or aggregates to form clusters.
where E catalysts and E slab are the total energies of the TM-doped BN and supports; and E TM1 , E TM2 , E bulk(TM1) and E bulk(TM2) are the energies of the isolated atom and the single TM atoms in their bulk structure.The calculated corresponding energies can be found in Figure 1b and Supplementary Table S1.As revealed in Figure 1b, it is interesting to note that all the E f and E ads − E coh values are well below zero, indicating the high stability of the dopant TM atoms.
With respect to DACs, the insertion of dopant TM pairs leads to the change of the local environment and electronic property [29].We hence examined the TMs-substrate interactions by performing the Bader charge analysis.From Figure 1c and Supplementary Table S2, we found significant electron transfer from the TM atoms to the V B -BN substrate, ranging from 1.15|e| to 1.97|e|, which supports the strong interactions between the TM atoms and the V B -BN substrate.However, the electronegativity also affects the electron transfer of the TM pairs embedded in V B -BN substrate.For instance, the Ti atom, which has the smallest electronegativity (1.54 e) among all the candidate metal atoms, exhibits a strong ability to bind to the substrate, accompanying large amounts of electron transfer.However, for Ru (2.2 e for electronegativity), the electron transfer is not promising among all the candidates.

N 2 Adsorption and Selectivity on TM 1 -TM 2 @V B -BN
Firstly, the capture and adsorption of N 2 are prerequisites, meaning that the Gibbs free energy change (∆G) for stable N 2 adsorption should be lower than 0 eV (∆G(*N 2 )< 0)) [30].Considering the adsorption of N 2 , the heteronuclear and homonuclear DACs provide two adsorption sites and one adsorption site, respectively.Consequently, we totally considered 16 adsorption configurations and calculated the ∆G(*N 2 ) values of N 2 terminated on TM 1 -TM 2 @V B -BN in end-on/side-on manners (see Table 1).Moreover, the NRR catalyst should effectively inhibit competitive HER, which is also an important indicator of the catalyst's efficiency.The surface of catalysts must provide enough active sites for N 2 activation, so as to facilitate the subsequent reactions.A stronger adsorption activity of the N 2 molecule is essential for high NRR selectivity, which ensures that the catalyst selectively adsorbs N 2 molecular prior to H [31,32].As a matter of convenience, we determined the catalytic selectivity of TM 1 -TM 2 @V B -BN through a comparison of the adsorption free energies of the *N 2 and *H atom.
As Figure 2a displays below, all structures have a negative ∆G (*N 2 ), indicating that the N 2 molecule can be stably anchored on the catalyst for subsequent hydrogenation reactions.For all the DACs, the capture of N 2 prefers side-on mode with lower ∆G (*N 2 ) than that of end-on configurations.Fe-Fe@V B -BN was not displayed, since it has the most negative ∆G(*N 2 ), at −2.12 eV, and would be detrimental to production and desorption of NH 3 .For catalyst selectivity, TM 1 -TM 2 @V B -BN in the lower-right region of Figure 2a is favorable for NRR.Most of the dopant pairs decorated V B -BN are within the area of ∆G(*N 2 ) < ∆G(*H), with the exception of Ti-Mo and Ti-Ru (see Table 1).Our results imply that these candidates are more inclined to adsorb *N 2 and that the adsorption of *N 2 is stronger than that of *H.This means that the coverage of N 2 molecules on the surface of these catalysts is relatively large; that is to say, these candidates show high NRR selectivity.However, the candidate Ti-Mo/Ti-Ru@V B -BN displayed lower NRR activity when N 2 molecules were adsorbed in end-on mode compared to HER, as they were deleted from the ensemble of catalysts candidates.

Catalytic Activity of TM1-TM2@VB-BN
After clarifying the N2 adsorption strength and selectivity, we further studied the catalytic activity of DACs.Generally, the inevitable first hydrogenation step (*N2 + H -+ e − →*N2H) requires a large energy consumption to break the inert N≡N triple bond and, hence, induces a positive ΔG.Moreover, the last hydrogenation reaction (*NH2 + H − + e − →*NH3) is usually achieved with an uphill in ΔG due to the strong combination between We also measured the N-N bond length (d N-N ) and the charge N 2 gained from the base material, and the relationships between the two are shown in Figure 2b,c.Interestingly, Figure 2b,c show that the greater the charge accumulation on N 2 , the longer the N-N bond length.Upon N 2 adsorption, d N-N in terminated N 2 is elongated and within the range of 1.14-1.15Å(end-on)/1.18-1.26Å(side-on) (vs1.12Å in free gas phase), implying that TM 1 -TM 2 @V B -BN candidates are effective in the activation of N 2 molecule.In addition, these data also indicate that the bond length of N 2 adsorbed in the side-on configuration is significantly longer than that in the end-on configuration.This can be attributed to the fact that when the TM adsorbs N 2 in a side-on configuration, both N atoms will bond with the TM and will mutual repulsion, resulting in a longer N-N bond. Figure 2b also shows that the catalysts transferred 0.35-0.52|e| to N 2 adsorbed through end-on mode.In contrast, our results indicate that N 2 adsorbed through side-on mode can gain more charge (0.46-0.86 |e|) from the catalysts, thus activating N≡N triple bond more effectively (see Supplementary Table S3).

Catalytic Activity of TM 1 -TM 2 @V B -BN
After clarifying the N 2 adsorption strength and selectivity, we further studied the catalytic activity of DACs.Generally, the inevitable first hydrogenation step (*N 2 + H -+ e − →*N 2 H) requires a large energy consumption to break the inert N≡N triple bond and, hence, induces a positive ∆G.Moreover, the last hydrogenation reaction (*NH 2 + H − + e − →*NH 3 ) is usually achieved with an uphill in ∆G due to the strong combination between the *NH 2 species and TM atoms [33].Therefore, we set the Gibbs free energy change of the first (∆G *N2→*N2H ) and the last hydrogenation step (∆G *NH2→*NH3 ) at less than 0.6 eV as a criterion to ensure a low limiting potential.This criterion is half the rate-determining step of pristine BN (1.20 eV), which is superior to the commonly used metal catalyst Ru (0001) (0.98 eV).
As portrayed in Figure 2d, the TM 1 -TM 2 @V B -BN catalysts in the upper-left region, corresponding to the structures in the upper-right panel of Figure 2b or Figure 2c, can meet our restricted value only for the first hydrogenation step.Obviously, these catalysts exhibit a significant activation of the N 2 molecule, which makes the first hydrogenation step easier to occur.However, too-strong N 2 adsorption would lead to the unfavorable formation and desorption of NH 3 ; hence, the ∆G values of the last hydrogenation reaction for these systems deviate from our criterion (0.6 eV).In contrast, those catalysts located in the bottomleft corner of Figure 2b or Figure 2c require a higher energy demand in the formation of *N 2 →*N 2 H.Those catalysts are finally found in the bottom-right area of Figure 2d and ruled out at this step.Consequently, only the candidates located in the bottom-left area with moderate N 2 adsorption strengths can reach the criteria.Ru-Mo@V B -BN, Ru-Ru@V B -BN and Ru-Ti@V B -BN show superior NRR activity, along with the least energy requirements for both the first and last protonation process.In conclusion, Ru-Mo@V B -BN, Ru-Ru@V B -BN and Ru-Ti@V B -BN are considered to be potential NRR catalysts, and they are discussed in the following sections.Other structures' values of ∆G for the two hydrogenation steps are listed in Supplementary Table S4, respectively.

Activation Mechanisms of N 2 on Ru-Mo/Ru/Ti@V B -BN
To reveal the active origin of NRR, we analyzed the interactions between metal active sites and reaction intermediates.As shown in Figure 3, we took the most promising candidates, Ru-Mo@V B -BN, Ru-Ru@V B -BN and Ru-Ti@V B -BN, as the typical systems for further analysis.Their charge-density-difference plots and key bond length are displayed in Figure 3. Remarkably, N 2 tends to interact with the top Ru atom; hence, we mainly analyzed the electronic transfer situation between the Ru atom and N 2 .Obviously, Figure 3 reveals a significant charge redistribution, which mainly occupies the antibonding orbitals of the adsorbed N 2 , leading to the elongation of the N-N bond.Moreover, d N-N (1.18/1.19Å) of N 2 absorbed in the side-on mode is longer than that in the end-on mode (1.14 Å), and this resonates with the amount of charge transfer between the TM atom and *N 2 .
ering a model of single-atom Ru anchored on VB-BN.Similarly, N2 molecules are adsorbed on the Ru active site in two manners, and the charge density differences and Bader charge analysis were also analyzed (see Supplementary Figure S3).For Ru@VB-BN, the single Ru atom donates 0.34|e| and 0.28|e| to *N2, which is less than the amount of charge transferred in the case of Ru-Mo/Ru/Ti@VB-BN.Consequently, our simulations confirm that N2 adsorbed on the DACs can gain more charge from the catalysts, thus activating the inert N≡N bond more effectively.In addition, we also compared the activation effect of the TM dimer on N 2 by considering a model of single-atom Ru anchored on V B -BN.Similarly, N 2 molecules are adsorbed on the Ru active site in two manners, and the charge density differences and Bader charge analysis were also analyzed (see Supplementary Figure S3).For Ru@V B -BN, the single Ru atom donates 0.34|e| and 0.28|e| to *N 2 , which is less than the amount of charge transferred in the case of Ru-Mo/Ru/Ti@V B -BN.Consequently, our simulations confirm that N 2 adsorbed on the DACs can gain more charge from the catalysts, thus activating the inert N≡N bond more effectively.
To uncover the underlying mechanism of N 2 activation, the partial density of states (PDOS) of bare and N 2 -terminated Ru-Mo/Ru/Ti catalysts was analyzed.Particularly, Figure 4a shows four characteristic orbitals, namely 2σ*, 1π, 3σ and 1π*, in free N 2 molecular within the range of −5 ~7 eV.As shown in Figure 4d-f, the electrons were obviously transferred from the occupied Ru 4d orbital to the empty π* orbital of N 2 after N 2 adsorption.As a result, the antibonding π* orbitals of adsorbed *N 2 on these three catalysts reveal a significant down-shifting in comparison with isolated N 2 molecular.Correspondingly, the unoccupied d orbitals of the Ru atom accepted the electrons from *N 2 , forming the Ru-N bonding state to strengthen the N 2 capture, as evidenced by the evident hybridization between Ru_4d and N_2p orbitals around the Fermi level.In addition, the d-band centers (ε d ) of the active Ru sites before the adsorption of N 2 molecules for the three catalysts were also calculated and marked in Figure 4.It can be seen from Figure 4c that the Ru-Ti structure is easier to adsorb and activate N 2 because its ε d is closer to the Fermi level.Compared with the homonuclear diatomic catalyst Ru-Ru@V B -BN, the introduction of coordinating metal atoms Ti and Mo affects the ε d of Ru atom to different degrees.Interestingly, the Ti atom causes Ru to obtain a higher ε d , which can be attributed to the higher electronegativity of Ti, hence making the activation of Ru-Ti@V B -BN easier for N 2 .For Ru@V B -BN, we also calculated its d-band center (see Supplementary Figure S4).In contrast, since the second metal atom on the surface also interacts with the adsorbate, the d orbitals of both metal atoms are involved to some extent.Their activities can be reasonably described by the central location of the PDOS corresponding to the d state.It can be seen from the numerical value that the diatomic catalyst has higher catalytic activity.Recently, Guo et al. suggested that the magnetic moment, µ B , of a metal surface plays an important role in NH 3 synthesis [34].Based on this, we investigated the µ B of the Ru active site in Ru-Mo/Ru/Ti@V B -BN systems (see Supplementary Table S5).The calculated µ B values for Ru atoms in Ru-Mo/Ru/Ti@V B -BN catalysts are 1.55, 1.08 and 1.60, respectively, while the Ru-Ti@V B -BN system evidently shows the largest magnetic moment.For this reason, the N≡N tripe bond of the N 2 molecule on the magnetic Ru-Ti@V B -BN monolayer is weakened, thus leading to stronger N 2 adsorption.
metal atoms Ti and Mo affects the εd of Ru atom to different degrees.Interestingly, the Ti atom causes Ru to obtain a higher εd, which can be attributed to the higher electronegativity of Ti, hence making the activation of Ru-Ti@VB-BN easier for N2.For Ru@VB-BN, we also calculated its d-band center (see Supplementary Figure S4).In contrast, since the second metal atom on the surface also interacts with the adsorbate, the d orbitals of both metal atoms are involved to some extent.Their activities can be reasonably described by the central location of the PDOS corresponding to the d state.It can be seen from the numerical value that the diatomic catalyst has higher catalytic activity.Recently, Guo et al. suggested that the magnetic moment, μB, of a metal surface plays an important role in NH3 synthesis [34].Based on this, we investigated the μB of the Ru active site in Ru-Mo/Ru/Ti@VB-BN systems (see Supplementary Table S5).The calculated μB values for Ru atoms in Ru-Mo/Ru/Ti@VB-BN catalysts are 1.55, 1.08 and 1.60, respectively, while the Ru-Ti@VB-BN system evidently shows the largest magnetic moment.For this reason, the N≡N tripe bond of the N2 molecule on the magnetic Ru-Ti@VB-BN monolayer is weakened, thus leading to stronger N2 adsorption.2.5.Full Conversion of N 2 to NH 3 on Ru-Mo/Ru/Ti @V B -BN Based on the above discussion, we focused on three catalysts (Ru-Mo@V B -BN, Ru-Ru@V B -BN and Ru-Ti@V B -BN) to search for the most promising catalyst.We canvassed the full reaction path search through three possible mechanisms: distal, alternative and enzymatic mechanisms [35][36][37].All of the optimized intermediates along the proposed reaction mechanisms are illustrated in Figure 5a, and their relevant ∆G values are summarized in Supplementary Tables S6-S11.As Figure 5b demonstrates, the optimal NRR mechanism for the Ru-Mo@V B -BN and Ru-Ru@V B -BN systems are both the distal mechanism, while the Ru-Ti@V B -BN catalyst prefers the enzymatic mechanism.Importantly, the rate-determining step (RDS) of all catalysts is the first hydrogenation process.Our calculated RDSs for Ru-Mo/Ru/Ti@V B -BN catalysts are 0.47, 0.54 and 0.40 eV, respectively.Apparently, Ru-Ti@V B -BN possesses the highest activity, with a limiting potential of −0.40 V through an enzymatic mechanism.It is worth noting that Ru-Ti@V B -BN exhibits comparable and even better NRR activity than most of the reported NRR SACs, such as Mo 2 @PTI (−0.53 V) [38], Ru@Ti 3 C 2 (−0.40 V) [39], Ru 2 @GDY(−0.55V) [40], Fe-MoS 2 (−1.02V) [41] and Ti@Ru(0001) (−0.47 V) [42].In addition, we also considered the effect of van der Waals correction.As shown in Supplementary Figure S5, the van der Waals correction has no effect on the rate-determining steps and overall reaction trend of NRR.Taking the optimal Ru-Ti@VB-BN as an example, we further examined the evolution behaviors of the Bader charge of each intermediate along the favorable enzymatic mechanism.During the whole NRR process, each intermediate was divided into three parts, namely the defective BN monolayer, the TM dimer and the adsorbed NxHy species.Figure 6a shows the variations of the Bader charge of three parts, and all the data are summarized in Supplementary Table S12.The adsorbed NxHy and VB-BN can gain electrons from the Ru-Ti pairs for the N2 adsorption and the former two hydrogenation steps.Then the electrons of the adsorbed NxHy backflow to the Ru-Ti dimer for the *NNH2 formation.In the following hydrogenation processes, NxHy species continue to acquire electrons from Ru-Ti, and the backflow of electrons occur until the subsequent desorption of NH3.Thus, the Taking the optimal Ru-Ti@V B -BN as an example, we further examined the evolution behaviors of the Bader charge of each intermediate along the favorable enzymatic mechanism.During the whole NRR process, each intermediate was divided into three parts, namely the defective BN monolayer, the TM dimer and the adsorbed N x H y species.Figure 6a shows the variations of the Bader charge of three parts, and all the data are summarized in Supplementary Table S12.The adsorbed N x H y and V B -BN can gain electrons from the Ru-Ti pairs for the N 2 adsorption and the former two hydrogenation steps.Then the electrons of the adsorbed N x H y backflow to the Ru-Ti dimer for the *NNH 2 formation.In the following hydrogenation processes, N x H y species continue to acquire electrons from Ru-Ti, and the backflow of electrons occur until the subsequent desorption of NH 3 .Thus, the defective BN monolayer acts as an electron acceptor during NRR, while Ru-Ti pairs always reserve as the electron donors in the process of trapping the intermediates.Furthermore, we also examined the bond-length variations of N-N, Ru-N and Ru-Ti bonds (see Figure 6b and Supplementary Table S13).Markedly, the N 2 molecule can be properly activated on Ru-Ti@V B -BN, as is proved by the constant increase (decrease) in Afterward, the thermal stability of Ru-Ti@VB-BN was further evaluated by AIMD simulations, under 400 K for 4 ps, with the time step of 2 fs, as shown in Figure 6d.Obviously, the total energy converges quickly, and no significant geometrical distortion occurs over the whole simulation time, implying the relatively high thermodynamic stability of Ru-Ti@VB-BN.Ru-Ru/Mo@VB-BN exhibits similar behaviors under the same AIMD simulation conditions (Supplementary Figure S6).On this basis, we can predict that the established Ru-Ti@VB-BN catalyst is stable at ambient condition and holds great promise to be  We further performed the projected crystal orbital Hamiltonian population (pCOHP) analysis to illustrate the bonding-strength evolutions of the Ru-N bond for each NRR species [43].Figure 6c shows an increased orbital interaction of the Ru-N bond as the hydrogenation step continues, and this can be confirmed by the increase of bonding contributions below the Fermi level.To gain a quantitative indicator, we calculated the integral value of COHP (ICOHP, a more negative value corresponds to a stronger TM-N interaction).The obtained ICOHP values of *N 2 , *NNH, *NHNH, *NH 2 NH and *NH 2 NH 2 are −1.14, −1.22, −1.29, −1.45 and −1.92, respectively, further proving the gradually strengthened Ru-N bonding.After the first *NH 3 desorbed, the bonding state of Ru-N reached the strongest interaction, with an ICOHP value of −2.25.For *NH 3 , the ICOHP value decreased to −1.05, indicating the weakened Ru-N bond to expedite the formation of the second NH 3 .
Afterward, the thermal stability of Ru-Ti@V B -BN was further evaluated by AIMD simulations, under 400 K for 4 ps, with the time step of 2 fs, as shown in Figure 6d.Obviously, the total energy converges quickly, and no significant geometrical distortion occurs over the whole simulation time, implying the relatively high thermodynamic stability of Ru-Ti@V B -BN. Ru-Ru/Mo@V B -BN exhibits similar behaviors under the same AIMD simulation conditions (Supplementary Figure S6).On this basis, we can predict that the established Ru-Ti@V B -BN catalyst is stable at ambient condition and holds great promise to be synthesized experimentally ahead.

Materials and Methods
All of the spin polarized density functional theory (DFT) calculations were performed by using the projector augmented wave (PAW) formalism and implemented in the Vienna Ab initio Simulation Package (VASP) [44,45].Exchange-correlation potential was treated by the generalized gradient approximation (GGA) in the form of the Perdew-Bruke-Ernzerhof (PBE) functional [46].A plane-wave basis set with an energy cutoff of 520 eV was chosen for all computations, based on our convergence test (see Supplementary Figure S1).A (4 × 4 × 1) supercell was used to model the defective BN monolayer by removing one B atom.To avoid interactions between periodic images, a vacuum of 20 Å was set in the z direction.A Monkhorst-Pack grid of 3 × 3 × 1 was used to sample the first Brillouin zones.The convergence limits for energy and force were set to 10 −5 eV and 0.02 eV Å −1 , respectively.The long-range van der Waals interactions were treated by Grimme's scheme (DFT-D3) [47].A Bader charge analysis was carried out to account for the quantitative description of charge distribution and charge transfer.In addition, the thermodynamic stability of the catalysts was investigated by using ab initio molecular dynamics (AIMD) simulations under the NVT ensemble [48].The AIMD calculations lasted for 4 ps, with a time step of 2.0 fs at 400 K.
The change of Gibbs free energy (∆G) of each elementary step for NRR reaction was calculated by the following equation: where ∆E denotes the reaction energy, which could be directly obtained from DFT calculations; and ∆E ZPE and ∆S are zero-point energy and entropy change at 298.15 K, which are computed from the frequency analysis [49].Moreover, ∆G pH is the correction of the H + free energy in the aqueous solution, which can be calculated by using the formula: where the k B is the Boltzmann constant and the pH is set to be zero in this work; and e and U refer to the number of electrons transferred and the applied electrode potential, respectively.The applied electrode potential of each step is defined as follows:

Scheme 1 .
Scheme 1.Representative structures of single-and double-atom catalysts reported in the past five years and their corresponding advantages and disadvantages.

Catalysts 2022 , 14 Figure 1 .
Figure 1.(a) Schematic configurations of TM1-TM2@VB-BN.(b) Calculated formation energy and the difference between the adsorption energy and cohesive energy of the TM1-TM2@VB-BN systems.(c) Charge transfer from the TM pair to the VB-BN substrate.

Figure 1 .
Figure 1.(a) Schematic configurations of TM 1 -TM 2 @V B -BN.(b) Calculated formation energy and the difference between the adsorption energy and cohesive energy of the TM 1 -TM 2 @V B -BN systems.(c) Charge transfer from the TM pair to the V B -BN substrate.

Figure 2 .
Figure 2. (a) Comparison of the HER limiting potentials and NRR limiting potentials via side-on and end-on mode.The relationship between the N-N bond length (dN-N) and the charge *N2 gained from the TM1-TM2@VB-BN substrate through (b) end-on and (c) side-on pattern.(d) Free energy change of the first hydrogenation step (∆G*N2→*N2H) and the last hydrogenation step (∆G*NH2→*NH3) for different DACs.

Figure 2 .
Figure 2. (a) Comparison of the HER limiting potentials and NRR limiting potentials via side-on and end-on mode.The relationship between the N-N bond length (d N-N ) and the charge *N 2 gained from the TM 1 -TM 2 @V B -BN substrate through (b) end-on and (c) side-on pattern.(d) Free energy change of the first hydrogenation step (∆G *N2→*N2H ) and the last hydrogenation step (∆G *NH2→*NH3) for different DACs.

Figure 3 .Figure 3 .
Figure 3. Optimized structures and change density differences of *N2 on Ru-Mo@VB-BN, Ru-Ru@VB-BN and Ru-Ti@VB-BN catalysts for (a) end-on and (b) side-on configurations.The key bond lengths are also given.

Figure 4 .
Figure 4. Partial density of states (PDOS) of (a) free N 2 molecule and Ru 4d orbitals in (a) Ru-Mo@V B -BN; (b) Ru-Ru@V B -BN and (c) Ru-Ti@V B -BN before N 2 adsorption.PDOS of N2p and Ru 4d orbitals in (d) Ru-Mo@V B -BN; (e) Ru-Ru@V B -BN and (f) Ru-Ti@V B -BN after N 2 adsorption in end-on (upper panel) and side-on manner (bottom panel).The purple dashed line represents the corresponding energy level of the d-band center (ε d ) of the Ru active site, and the Fermi level is set to 0 eV.

Figure 5 .
Figure 5. (a) Schematic diagram of the possible reaction mechanisms and geometry structures of the corresponding intermediate products during the NRR process.(b) Gibbs free energy diagram for NRR on Ru-Mo@VB-BN, Ru-Ru@VB-BN and Ru-Ti@VB-BN catalysts.

Figure 5 .
Figure 5. (a) Schematic diagram of the possible reaction mechanisms and geometry structures of the corresponding intermediate products during the NRR process.(b) Gibbs free energy diagram for NRR on Ru-Mo@V B -BN, Ru-Ru@V B -BN and Ru-Ti@V B -BN catalysts.