Boron-Doped MXenes as Electrocatalysts for Nitrogen Reduction Reaction: A Theoretical Study

Electrocatalytic nitrogen reduction reaction (NRR) is a promising and sustainable approach for ammonia production. Since boron as an active center possesses electronic structure similar to that of transition metals with d-orbitals (J. Am. Chem. Soc., 2019, 141 (7), 2884), it is supposed to be able to effectively activate the triple bond in N2. MXenes can be applied as substrates due to the large specific surface area, high conductivity, and tunable surface composition. In this work, the catalytic performance of a series of MXenes-supported single boron atom systems (labeled as B@MXenes) has been systematically studied by using density functional theory (DFT). B@Nb4C3O2, B@Ti4N3O2, and B@Ti3N2O2 were screened out owing to outstanding catalytic activity with limiting potentials of −0.26, −0.15, and −0.10 V, respectively. Further analysis shows that the unique property of boron that can intensely accept lone pair and back-donate the anti-bond of nitrogen contributes to the activation of inert triple bond. This work provides a new idea for the rational design of NRR catalyst and is of great significance for the future development of nitrogen reduction catalysts.


INTRODUCTION
NH 3 is one of the most important basic chemicals used in the downstream manufacturing of various products in the agricultural industry, synthetic fiber, and fine chemicals (Giddey et al., 2013;Zhan and Zhang, 2021). Moreover, it is supposed to be promising energy storage intermediate. Although the atmosphere consists of more than 70% of nitrogen and the reduction of nitrogen to produce NH 3 is an exothermic reaction under standard conditions, the nitrogen fixation from N 2 to NH 3 is still tough because of the inherent nature of the N≡N bond with extremely high bond energy (Gambarotta and Scott, 2004). Currently, industrial production of NH 3 relies on the high energy-intensive Haber-Bosch process, which results in a huge amount of energy consumption and carbon emissions around the world (Cui et al., 2018).
In contrast, inspired by the biological N 2 fixation in bacteria (Singh et al., 2021), the electrochemical nitrogen reduction reaction (NRR) provides a green and sustainable alternative to produce NH 3 at relatively mild conditions (Liu et al., 2021a). The development of highly active and stable catalysts for the electrochemical NRR is urgently demanded and challenging. In many cases (Kugler et al., 2015;Manjunatha and Schechter, 2018;Liu et al., 2020a;Wang et al., 2020a;Wang et al., 2020b;Cai et al., 2020;Yang et al., 2020;Zheng et al., 2020;Dai et al., 2021), electrocatalysts contain transition metals as the catalytic center. The empty d-orbitals accept the σ-electron pair and the filled d-orbitals back-donate the anti-bonding π orbitals, activating the inert bond (Rao and Rao, 1991;Légaré et al., 2018), while there is much room for improvement of Faraday efficiency and reaction selectivity. In recent years, nonmetallic boron-based catalysts, such as boron-doped graphene and dicoordinate borylene, have attracted great interest in N 2 reduction (Légaré et al., 2018;Yu et al., 2018). The boron sites can provide enhanced affinity to N 2 molecules, leading to high NH 3 productivity (Hering-Junghans, 2018).
Introducing boron dopant into an appropriate substrate is a strategy to construct boron-centered catalytic system (Ji et al., 2019;Liu et al., 2019;Zheng et al., 2019). The transition metal carbides/ nitrides/carbonitrides, known as MXenes, is a new class of twodimensional materials . MXenes materials are generally expressed in chemical formula as M n + 1 X n T x . M represents early transition metals (Hf, Zr, Nb, Ti, and Ta), X represents the nonmetal element C or N, and T refers to the surface terminal groups, like hydroxyl group, −F, or −O. Normally, n is 2 or 3. Because of the high conductivity, large specific surface area, and tunable surface composition, it has attracted extensive attention in the field of electrocatalysis and is regarded to be an ideal substrate to support the active center (Liu et al., 2020b;Wang et al., 2021). For example, Mo@Mo 2 CO 2 Mo 2 TiC 2 -based systems, Mo 2 C, and W 2 C have been reported as novel electrocatalysts for the NRR (Shao et al., 2018;Huang et al., 2019;Gao et al., 2020). Previous works mainly focused on three-and five-layered M 2 XO 2 as electrocatalyst for the NRR, such as Mo@Ti 2 NO 2 (Cheng et al., 2019), Mo@Mo 2 CO 2 (Huang et al., 2019), Ru@Mo 2 CO 2 (Peng et al., 2020), or B@W 2 CO 2 (Zheng et al., 2019). In this work, we further consider MXenes composed of more layers that possess high thermodynamic and dynamic stability. The nonmetal B atom is doped into M 3 X 2 O x or M 4 X 3 O x which were selected out from the C2DB database (Rasmussen and Thygesen, 2015;Gjerding et al., 2021). Fourteen boron-doped MXenes (B@MXenes) have been studied as potential electrocatalysts. Three candidates are screened out, labeled as B@Nb 4 C 3 O 2 , B@Ti 4 N 3 O 2 , and B@Ti 3 N 2 O 2 . The nitrogen reduction processes are explored to evaluate the catalytic performance. In order to understand the promoted activity of the boron-containing systems, charge density differences and projected density of states (PDOS) are analyzed. Moreover, AIMD simulation is employed to verify the stability. This work may broaden the application of MXene-based materials and has great significance for the future development of electrocatalyst for the NRR. We hoped our work can inspire more researches on M 3 X 2 O x or M 4 X 3 O x that have the potential to be implemented experimentally.

COMPUTATIONAL METHODS
All density functional theory (DFT) computations were performed by using the Vienna ab initio simulation package (VASP) with the projector augmented wave (PAW) (Blöchl, 1994;Kresse and Furthmüller, 1996;Kresse and Joubert, 1999;Grimme, 2006). The van der Waals interaction was considered by the empirical correction of the Grimme's scheme (DFT-D2). The generalized gradient approximation (GGA) in Perdew-Burke-Ernzerh (PBE) functional was adopted (Perdew and Wang, 1992;Perdew et al., 1997). Energy cutoff of 450 eV for the plane wave basis was applied. The Brillouin zone was sampled by using a 3 × 3 × 1 Monkhorst-Pack k-points mesh. Both energy cutoff and k-points mesh have been checked based on calculations of a typical system, as shown in Supplementary Figures S1, S2. To avoid interactions between periodic images, a vacuum distance of 15 Å was imposed between two periodic units. All structures were optimized within energy and force convergences of 10 −5 eV/atom and 0.02 eV/Å, respectively. AIMD was performed under NVT ensemble by using Nosé-Hoover thermostat to control the temperature at 300 K. The time step is 2 fs, and each simulation run for 5,000 steps, namely, 10 ps. The energy barrier for B atom escaping from the vacancy is simulated by conducting the climbing-image-nudged elastic band (CI-NEB) methods. For each path, three intermediate images were inserted. The convergence criteria for energy and force were set to be 10 −7 eV and 0.05 eV/Å, respectively.
The computational hydrogen electrode model (CHE) was employed to calculate the Gibbs free energy change (ΔG) of each elementary step (Norskov et al., 2004;Skúlason et al., 2012). According to the CHE model, the ΔG value can be determined by the following Eq. 1: where ΔE is the internal energy change directly calculated from DFT, ΔE ZPE is the change of zero-point energy, T is the room temperature (300 K), ΔS is the entropy change, and ΔG U is the free energy contribution related to electrode potential U. ΔG pH is the free energy correction of pH, which can be expressed as follows: where k B is the Boltzmann constant and the value of pH is assumed to be zero. The activity was evaluated through the limiting potential (U L ) vs. RHE that was defined as follows: where ΔG is the free energy of the potential-determining step.
To evaluate the interaction between catalyst and nitrogen molecule, the adsorption energy (ΔE ad ) is calculated as follows: where E total , E B@MXene , and E N2 are the energies of the B@MXene-N 2 complex, B@MXene catalyst, and N 2 molecule directly obtained from DFT calculation, respectively. The binding energy of B atom (E b ) is calculated to evaluate the feasibility of the designed structures:

RESULTS AND DISCUSSION
In experiment, the oxygen-covered MXenes are stable, while the coverage of surface groups is normally far from uniform, that is, there can be plenty of vacancy (Xie et al., 2014;Karlsson et al., 2015). It is possible to introduce boron dopant into oxygen vacancy, shown as a B-doped M n + 1 X n O x system, or labeled as B@M n + 1 X n T x in this work. Herein, fourteen thermodynamically stable MXenes were selected from structure database (Groom et al., 2016) to construct corresponding B@MXenes models. In Figure 1, B@Ti 3 C 2 O 2 is taken as an example to show the typical atomic model. Intense adsorption of gas phase N 2 on the surface of the catalyst is suggested to be the prerequisite for an efficient NRR process, which may activate the inert triple bond. Among the nonmetal elements, B atom possesses the combination of empty orbitals and occupied orbitals in specific direction; thus, the incorporated B site can serve as an active center for the capture of N 2 . Compared to previous results, the adsorption of N 2 on B atom is preferable to that on O vacancy, indicating that the B center plays a key role in the electrocatalytic NRR process (Zheng et al., 2019). There are mainly four possible pathways for N 2 electrochemical reduction to NH 3 : distal, alternative, enzymatic, and mixed paths. Two typical N 2 adsorption configurations on B@MXenes systems are taken into account, including end-on and side-on configurations, which determine the following reduction path (Zhao and Chen, 2017). The simulated adsorption energies are shown in Figure 2. It is noticed that on all the fourteen MXene layers, the end-on pattern is considerably more stable than the side-on one and initializes all the following reduction paths. Furthermore, we also compare the hydrogen adsorption energy since the hydrogen evolution reaction (HER) is a competitive process. Both oxygen and boron sites are considered as the adsorption sites for some typical systems. As listed in Supplementary Table S1, the hydrogen adsorption on the O site is also weaker than that on B dopant. Although the desorption energy of hydrogen close to zero is normally related to high activity of the HER, the more negative adsorption energy of nitrogen indicates that the active center can be occupied by nitrogen, hindering the competitive hydrogen adsorption. The weaker hydrogen adsorption than nitrogen adsorption in the end-on configuration indicates a high selectivity of the NRR over the HER on all the considered B@MXenes.
Normally, the potential-determining step (PDS) of the NRR occurs in the first hydrogenation step (*N 2 →*NNH) or the last one (*NH 2 →*NH 3 ), which depends on the binding strength of the nitrogen-containing intermediate to the catalyst surface (Skúlason et al., 2012). If the bonding strength is too weak,  the breaking of the nitrogen triple bond is probable to be the PDS of the whole process. On the contrary, the strong bonding strength leads to the difficulty of the *NH 3 formation in the last hydrogenation step. We screened fourteen B@MXene systems according to the free energy changes of the first and the last hydrogenation steps as the criteria. As shown in Figure 3, Ti 4 N 3 O 2 , Ti 3 N 2 O 2 , and Nb 4 C 3 O 2 in the lower left corner exhibit relatively low energy barriers in both the first and the last key steps. Subsequently, we further explore the whole NRR processes of these three systems. Due to the high energy barrier of N 2 dissociation, the dissociative step was generally not considered in the hydrogenation process. As illustrated in Figure 4A, there are four possible associative NRR mechanisms, that is, distal, alternative, enzymatic, and mixed on the catalyst surface (Shi et al., 2019). The enzymatic path is excluded in our study since the side-on configuration is less stable than the corresponding endon one that generally initializes the enzymatic mechanism. Starting from the end-on configuration, free energy diagrams including distal, alternative, and mixed paths for three candidates (B@Nb 4 C 3 O 2 , B@Ti 4 N 3 O 2 , and B@Ti 3 N 2 O 2 ) are presented in Figures 4B-D.
On B@Nb 4 C 3 O 2 , the adsorbed N 2 (*N 2 ) can interact with the proton-electron (H + /e − ) pair to form NNH*. The Gibbs free energy change of the first hydrogenation step (*N 2 →*NNH) is 0.26 eV. Then, the H + /e − pair can attack 1 N atom of the adsorbed NNH species to form *NHNH or *NNH 2 intermediate, leading to the alternative or distal pathway, respectively. The calculated free energy change of the *NNH→*NNH 2 step is −0.70 eV, rather lower than that of the *NNH→*NHNH step (0.25 eV). Therefore, B@Nb 4 C 3 O 2 prefers to catalyze the NRR through the distal path. Once *NNH 3 is formed, the first NH 3 molecule can be readily desorbed from the surface, leaving a single nitrogen atom on the B site. In the subsequent reaction steps, three H + /e − pairs can continuously react with the remaining *N, forming *NH,  *NH 2 , and *NH 3 with energy changes of −1.01, −0.33, and −0.08 eV, respectively. Among all the elementary steps, the formation of *N 2 H (*N 2 →*NNH) is the PDS with the maximum free energy demand of 0.26 eV. In other words, when the external potential reaches to −0.26 V, the reaction can take place thermodynamically. The whole process proceeds the distal mechanism through the intermediates as *N 2 → *NNH → *NNH 2 → *N → *NH → *NH 2 → *NH 3 . Even though NH 3 desorption from the activation center requires a large amount of energy that may inhibit the regeneration of catalyst, we suppose that the strongly acidic electrolyte results in the formation of NH 4 + and promotes desorption of the products (Ling et al., 2018;Zheng et al., 2019;Guo et al., 2020).
Similarly, the NRR pathways on B@Ti 4 N 3 O 2 and B@Ti 3 N 2 O 2 are examined from the end-on adsorption configuration. The Gibbs free energy changes of the first hydrogenation step (*N 2 →*NNH) on B@Ti 4 N 3 O 2 and B@Ti 3 N 2 O 2 are both 0.04 eV. Afterward, the *NNH is easily protonated to the energetically more favorable *NNH 2 , releasing energies of 0.76 and 0.60 eV, respectively. Then, the H + /e − pair can attack *NNH 2 to form *NHNH 2 or *NNH 3 . . On B@Ti 4 N 3 O 2 and B@Ti 3 N 2 O 2 , the free energy changes of the *NNH 2 →*NHNH 2 step are 0.15 and 0.10 eV, respectively, while the free energy changes of the *NNH 2 →*N step are 0.30 and 0.22 eV, respectively. Thus, mixed paths through *NNH 2 →*NHNH 2 step are favorable on these two Ti-containing systems. In the following steps, three H + /e − pairs FIGURE 5 | Charge density differences of N 2 adsorbed on three B@ MXenes surfaces. The charge depletion and accumulation are depicted in cyan and yellow, respectively. The isosurface value is set to 0.001 e/Å 3 . continuously attack *NHNH 2 , forming *NH, *NH 2 , and *NH 3 . Among all the elementary steps, the formation of NHNH 2 * (*NNH 2 →*NHNH 2 ) is the PDS with maximum free energy changes of 0.15 or 0.10 eV on B@Ti 4 N 3 O 2 or B@Ti 3 N 2 O 2 , respectively. Overall, the NRR processes proceed within the mixed mechanism through the following intermediates: *N 2 → *NNH → *NNH 2 → *NHNH 2 → *NH → *NH 2 → *NH 3 . All the three candidate systems exhibit quite low barriers of PDSs, requiring less negative limiting potentials compared with other reported investigations including metals (−0.45 V for Ru) and single-atom catalyst (−0.69 V for Ti on graphene) (Skúlason et al., 2012;Liu et al., 2021b).
To gain deep insight into the high activity of the three B-doped systems, we explored the electronic structures. Charge density differences of N 2 adsorption configurations are analyzed as shown in Figure 5, indicating significant charge transfer from the B atom as the active site to the adsorbed N 2 molecule. According to the so-called push-pull hypothesis (Geri et al., 2017), B@MXenes can "pull" the lone-pair electrons from the N 2 and simultaneously "push" electrons back into the anti-bonding orbitals of N 2 . In other words, the N 2 molecule donates its lone-pair electron to the empty 2p z orbital of boron. Meanwhile, the antibonding of N 2 molecule can accept electrons from B@ MXenes catalysts, which helps trigger the NRR. As shown in Figure 5, the increased electron density between nitrogen and the surface is related to the electrons that are pulled from nitrogen by the metal atom. On the other hand, electrons also accumulate on the two nitrogen atoms, representing the π* orbital that accepts "pushed" electrons from surface. In addition, we plotted the projected density of states (PDOS) of the adsorption configurations and free N 2 molecules ( Figure 6) (Guo et al., 2020). Compared to the molecular orbitals of free nitrogen, the interaction with the active B center apparently downshifts the band of nitrogen, stabilizing the adsorbed N 2 . The empty p-orbital of B atom receives electrons from nitrogen molecule, increasing the adsorption strength. Especially, in these two Ti-containing systems, the B dopants possess obvious occupation under and close to the Fermi level, which feedbacks electrons to the antibonding orbital of nitrogen and results in further activated nitrogen as well as lower barriers than that of the B@Nb 4 C 3 O 2 system. Besides, to evaluate the feasibility and stability of the designed systems, we calculate the dopant energies of the three potential   Figure 7A. The kinetic barriers are 1.05, 1.59, and 1.80 eV for Nb 4 C 3 O 2 , Ti 3 N 2 O 2 , and Ti 4 N 3 O 2 , respectively. These relatively high barriers indicate that the boron localized in the defect hardly move out. Herein, these materials should be kinetically stable as catalysts. Moreover, ab initio molecular dynamic simulations are performed for these three B@ MXenes. In 10 ps at 300 K, there is no obvious structural deformation, and the B atom stays in the vacancy, also confirming the stability of these systems ( Figures 7B-D).

CONCLUSION
In summary, we conducted DFT calculations to evaluate the potential of B@MXenes systems in the electrochemical NRR. Based on energy changes of the first and last hydrogenation steps, B@Nb 4 C 3 O 2 , B@Ti 4 N 3 O 2 , and B@Ti 3 N 2 O 2 are screened out as the most possible candidates with high catalytic activity. Subsequently, the reaction mechanisms are revealed to predict the nitrogen reduction performance of these three catalysts. On B@ Nb 4 C 3 O 2 catalyst, transfers of electron-proton pairs proceed following the distal mechanism, with the limiting potential of −0.26 V, whereas the NRR processes on the two Ti-containing systems go through the mixed pathway, with quite low limiting potentials of −0.15 and −0.10 V, respectively. In order to understand the good catalytic activity of the three structures on the NRR, their electronic structures are further analyzed. We found that N 2 was captured by these three B@MXenes to trigger the NRR through the so-called push-pull hypothesis. The PDOSs show that the p-band of B atom close to the Fermi level in these Ti-containing systems contributes to the enhanced interaction with inert N 2 and the low reaction barriers. Furthermore, the stabilities of the B-doped systems are verified by the kinetic barriers of B-dopant diffusion and AIMD simulation. Our findings not only prove that the three catalysts, namely, B@Nb 4 C 3 O 2 , B@Ti 4 N 3 O 2 , and B@ Ti 3 N 2 O 2 , may exhibit outstanding performance on NRR catalysis but also provide further insight into the rational design of NRR catalysts. It is expected to stimulate deep research on exploration of novel catalysts in the field of nitrogen fixation and also expand the application range of MXenes-based materials.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
YW, XQ, and GZ performed DFT calculation. ZT, YW, and QZ wrote the manuscript. All authors discussed the results on the manuscript.

ACKNOWLEDGMENTS
The author thanks the high-performance computing center of Nanjing University and the Supercomputer Center of Ningbo Institute of Materials Technology and Engineering, CAS, for providing computing resources.