DFT Investigation of Ammonia Formation via a Langmuir-Hinshelwood Mechanism on Mo-Terminated δ-MoN(0001).

In this work, we employed density functional theory to elucidate the energetics associated with elementary steps along a Langmuir-Hinshelwood mechanism for the Haber-Bosch synthesis of ammonia from N2 and H2 on a hexagonal, Mo-terminated molybdenum nitride surface. Using nudged elastic band calculations, we determined the energy barriers involved in the reaction processes. An active site consisting of four nearest-neighbor Mo atoms, previously identified as an active site on similar surfaces, was chosen to investigate the reaction processes. Using this approach, we calculate a barrier of ∼0.5 eV for the dissociation of N2. The superior activity of the dissociation of the strong N2 bonds is rationalized based on the unique geometric and electronic configurations present at these active sites. Despite the favorable energetics for nitrogen dissociation, the energy cost for hydrogenation of NH x (0 ≤ x ≤ 2) species is shown to be energetically limiting for the formation of ammonia through the Langmuir-Hinshelwood mechanism at these sites, with elementary step activation barriers calculated to be as large as ∼2 eV. A comparison to Haber-Bosch results derived from a similar γ-Mo2N model system suggests the relative independence of surface chemistry and bulk stoichiometry for rhombic Mo4 active sites present on molybdenum nitrides.

Introduction

Ammonia is a very important compound used in the production of many useful goods, for example, fertilizers, other nitrogenous compounds, cleaning agents, and antiseptics, and for fermentation. Industrial synthesis of ammonia is carried out via the Haber–Bosch process, which consists of the direct reaction of nitrogen with hydrogen over iron-based catalysts. Efficiently achieving usable yields of ammonia requires ultra-ambient reaction conditions (P ≥ 100 atm. and T ≥ 400 °C), which makes this process extremely energy intensive. For perspective, it is estimated that about 2% of the world’s energy usage is consumed by the Haber–Bosch process,[1] with most of that supply being extracted from fossil fuels. Apart from energy utilization, this reaction also currently depends upon fossil fuels to provide hydrogen gas feedstocks. In an effort to produce relatively “green” ammonia, many researchers have devoted themselves to the search for better catalysts amenable to this process. Though Ru-based catalysts are found to be more efficient than iron ones,[2] their use is economically hindered by limited supplies of terrestrial ruthenium. Search for a more cost-effective catalyst is required to decrease the energy requirements for this process and bring down the cost of ammonia (in terms of energy consumption, materials costs, and environmental impact).

Thanks in large part to the Nobel Prize winning work of Ertl,[3] conventional Haber–Bosch reactions are known to follow a Langmuir–Hinshelwood mechanism by activating N≡N dissociation on fourfold Fe sites following electron transfer from the substrate. A similar mechanism is also predicted for Ru-based catalysts when modeling activity at five-fold-coordinated Ru steps.[4,5] For these active monometallic catalysts, adsorption and dissociation of triple-bonded N2 is the most energy-intensive step in the reaction process.[6] More generally, Nørskov and his coworkers have shown that N2 adsorption strength can be used as a proxy to predict the Haber–Bosch activity on different catalytic interfaces.[7] If nitrogen adsorbs too weakly on a material, then the reaction is limited by the inefficient dissociative adsorption of N. On the other hand, if N2 adsorbs too strongly, then the reaction becomes limited by inefficient nitrogen hydrogenation and subsequent NH3 product desorption.[8] Hence, optimal catalysts should present intermediate N2 adsorption energies to best balance the needs of these opposing processes in accordance with the Sabatier principle.

In the past 2 decades, early transition-metal nitrides (TMNs) have been shown to activate ammonia formation.[9−12] Hybridization between orbitals of the metals and non-metals in these materials has been shown to lower the Fermi energy of these systems (relative to the pure metals), such that some typically less-active metals, like molybdenum, can be used to create binary complexes with electronic properties more closely approximating those exhibited by more catalytically active late transition metals.[13] Introduction of non-metal atoms (e.g., N) may also potentially change interatomic spacings in manners structurally favorable for improved adsorption and activation of N2 molecules on such surfaces. Beyond Haber–Bosch, TMNs also show promising results as active sites for electrochemical conversion of dinitrogen to ammonia.[14,15] Within the TMN family, molybdenum nitrides have received considerable attention, with several recent studies suggesting enhanced catalytic properties important to the Haber–Bosch reaction process.[16−20] In addition to presenting promising catalytic properties, the molybdenum nitride family of materials also exhibit a host of other physical properties of broad scientific interest. For example, several crystalline phases are known to be low-temperature superconductors,[21,22] and hexagonal δ-MoN exhibits a hardness comparable to diamond and boron nitride.[23]

Several groups have successfully demonstrated different approaches allowing for the controllable synthesis of various phase-pure molybdenum nitride structures.[24−26] Of particular interest to this paper is the recent work of Khaniya and Kaden,[24] which reports the epitaxial growth of MoN films on Ru(0001) via an ion-assisted and physical vapor deposition approach. The films produced in that work appear ripe for direct comparisons between surface-science experiments and density functional theory (DFT) predictions due to their single-crystalline and atomically planar terminations, which exhibit a long-range order and elemental stoichiometry compatible with δ-MoN(0001).

A computational work done by Zhao et al. investigated the catalytic properties of a γ-Mo2N(111) surface used to activate the Haber–Bosch reaction.[8] In that case, a rhombus containing four nearest-neighbor Mo atoms (Mo4) was identified as an active site for the adsorption of dinitrogen, with N≡N bond dissociation requiring only 0.58 eV when binding in highly coordinated arrangements at these sites. By contrast, subsequent NH hydrogenation at the same site was shown to be rate-limiting, with step-wise barriers reaching as high as 2.01 eV. Enhanced surface nitridation (likely to occur under Haber–Bosch conditions) was shown to generally lower the hydrogenation barriers, while impacts on nitrogen dissociation were shown to depend more strongly (and anisotropically) on N adatom locations relative to adjacent rhombic Mo4 ensembles.

The first atomic layer in Mo-terminated γ-Mo2N(111) and δ-MoN(0001) is near-identical, suggesting that similar Mo4 site behavior might be expected from the δ-phase interface (Mo-terminated surfaces of γ-Mo2N(111) and δ-MoN(0001) are shown in Figure for comparison). Based on (i) previous predictions indicating advantageous conditions for Haber Bosch catalysis on rhombic configurations of the nearest neighbor Mo sites present on Mo-terminated γ-Mo2N(111), (ii) the existence of near-identical surface site configurations on Mo-terminated δ-MoN(0001), and (iii) the existence of an established recipe for controllable growth of single-crystalline, hexagonal MoN thin films, we aim here to compare and contrast predictions for the energetics associated with each of the elementary steps along a Langmuir–Hinshelwood mechanism for ammonia production on a δ-phase model to those previously reported on analogous γ-phase sites. Results from this comparison may help further disambiguate the degree to which surface versus bulk modification (i.e., N enrichment/depletion) might be expected to affect the resultant Haber–Bosch surfacechemistry in subsequent experimental investigations on similar model systems.

Figure 1

(a) γ-Mo2N(111) and (b) δ-MoN(0001) showing various adsorption sites considered in this work. Purple and blue spheres denote Mo and N atomic positions, respectively. H atoms (present in subsequent figures) are represented using white spheres.

Results and Discussion

Adsorption Characteristics

To inform our selection of initial state (IS) and final state (FS) configurations used in the nudged elastic band (NEB) calculations, we first considered several possible adsorption geometries for each of the sequentially produced atomic/molecular species required to form ammonia via a Langmuir–Hinshelwood mechanism (i.e., N2, N, H2, H, NH, NH2, and NH3). A minimum of four distinct adsorption sites (bridge, fcc, hcp, and top) were considered in each case (see Figure ). For clarity, the four sites can be distinguished as follows: top sites refer to adsorbates placed directly above surface Mo atoms, bridge sites refer to adsorbates placed directly between two neighboring surface Mo atoms, hcp sites refer to adsorbates placed in three-fold hollows between three neighboring surface Mo atoms in positions directly above second-layer N atoms, and fcc sites refer to adsorbates placed in threefold hollows between three neighboring surface Mo atoms in positions directly above areas without second-layer N atoms. In order to clearly represent the adsorption configurations of the adsorbates, both “nearest bond lengths d” and “adsorbate heights [h]” are mentioned along with adsorption characteristics. Nearest bond lengths and adsorbate heights are calculated as the distance from the vertical position of the lowest atom within the adsorbed species to the nearest surface Mo atom and the average vertical position of top-layer Mo atoms, respectively.

Adsorption of N and N2

Of the four N/MoN trial configurations, only three proved stable. N atoms placed in bridge sites spontaneously migrate to hcp sites upon relaxation (see Table ). Side and top views of the three stable binding motifs are shown in Figure S1. Of the stable configurations, individual N atoms were found to adsorb most strongly in hcp sites (Eads = 3.20 eV), indicating a preference for terminal adsorption in sites having continued registry with the remainder of the bulk structure. Top and fcc sites also provide stable, albeit less favorable configurations (Eads = 0.75 and 2.55 eV, respectively). Whereas weaker adsorption at top sites results in considerably less substrate-to-N electron transfer (0.67 e–), stronger adsorption at both hcp and fcc sites results in a transfer of ∼1.05 e– onto the newly bound N adatom. Charge transfer from the substrate to electronegative N atoms is indicative of a new bond formation between the two.

Table 1

Adsorption Characteristics of N on MoN

	Eads (eV)	d [h] (Å)	charge transfer (e–)
N
fcc	2.55	2.02 [1.30]	1.07
hcp	3.20	2.01 [1.21]	1.05
top	0.75	1.72 [1.72]	0.67

To accommodate the greater degrees of freedom associated with the molecular adsorption of dinitrogen, the number of trial configurations was increased from four to eight when searching for stable N2/MoN binding arrangements. A first set of molecules were placed in vertical configurations (with the N2 bond aligned along the surface normal) at top, hcp, and fcc sites. A second set were placed in horizontal configurations with the N2 molecule centered above the same series of positions in addition to the bridge site (oriented with the N≡N bond aligned along the Mo–Mo bridge bond). Within this second set, two azimuthally rotated N2 configurations (N≡N surface orientations rotated by 30° with respect to one another) were separately considered at top sites to differentiate the molecules oriented parallel to the hexagonal arrangement of Mo rows from those bridging between threefold hollows on the opposite sides of the host atom. Optimization of these trial structures converged upon four stable adsorption arrangements depicted and characterized in Figure and Table , respectively.

Figure 2

Side and top views of stable adsorption configurations of dinitrogen on δ-MoN(0001). (a,e) bridge, (b,f) fcc-top, (c,g) top-h, and (d,h) top-v.

Table 2

Adsorption Characteristics of N2 on MoN

	Eads (eV)	d [h] (Å)	d(N–N) (Å)	charge on first N (e–)a	charge on second N (e–)	charge transfer (e–)
N2
bridge	2.54	2.10 [1.30]	1.33	5.67	5.61	1.31
fcc-top	1.95	2.12 [2.00]	1.30	5.67	5.35	1.04
top-h	1.07	2.35 [2.30]	1.16	5.12	5.25	0.39
top-v	1.92	2.07 [2.07]	1.13	5.43	4.86	0.29

“First N” refers to nitrogen atoms positioned closer to the surface plane in Figure .

Figure 3

Relative potential energy profile associated with N2 dissociation on δ-MoN(0001). A labeled space-filling model of the “bridge” N2/δ-MoN structure is provided in the top left corner. Remaining insets show top (upper) and side (lower) views of IS, TS, and FS configurations associated with the selected reaction process.

Most trial structures significantly distort to take on one of two new binding arrangements following relaxation (distinguished as “bridge” and “fcc-top”), suggesting a preference for non-vertical adsorption in more highly coordinated configurations. By contrast, the remaining arrangements (distinguished as “top-v” and “top-h”) indicate the stable adsorption of molecules remaining essentially unchanged relative to their trial configurations following relaxation. Comparing across these two groups, we find greater adsorption energy, MoN-to-N2 electron transfer, and N≡N bond distortion for the statistically more likely arrangements resulting from N2 coordination to multiple Mo atoms. Consistent with these tendencies, the most strongly bound configuration (bridge; Eads = 2.54 eV) also features the greatest degree of N2–Mo coordination (both N atoms bind over three-fold hollow sites), the closest average molecular proximity to the surface (1.3 Å), the greatest N2 bond elongation (0.20 Å growth vs gas phase geometry), and the largest amount of charge transferred from the surface to the molecule (1.31 e– evenly distributed across both atoms), such that the atoms in the molecule most closely resemble those in the two most favorable atomic adsorption configurations. Perhaps unsurprisingly, this arrangement bares a near-perfect agreement to that previously identified as the active configuration for N2 dissociation on γ-Mo2N(111),[8] including a mild distortion of the Mo4 site upon adsorption (0.4 Å contraction of the Mo atoms positioned at the vertices spanning the long axis of the rhombic ensemble). By comparison, top site configurations transfer far less charge onto the adsorbing molecule and leave the N2 bond lengths unperturbed, suggesting poor activation toward dissociation. Adsorption in the fcc-top configuration results in intermediate effects, with the atom coordinating to the threefold hollow (designated as the “first N” in Table ) exhibiting charge localization characteristics similar to atoms in the bridge configuration, and the atom coordinating to the top site taking on characteristics similar to atoms in the top-h configuration by Bader charge analysis. Given similar bond elongation characteristics (vs bridge site adsorption), this latter N2/MoN adsorption configuration appears to provide an alternative activation site also likely to be suitable for subsequent dissociation.

To allow for more direct comparison to previous literature, we reoptimized the pristine MoN, gas-phase N2, and N2/MoN structures to calculate a new bridge site Eads value using the Perdew–Burke–Ernzerhof (PBE) functional in a manner consistent with that used for the γ-Mo2N(111) model system.[8] Changing from PBE to optB88 causes Eads(bridge) to drop from 1.61 to 2.54 eV, indicating a 0.93 eV increase in this value presumed to predominantly derive from the inclusion of van der Waals interactions in the optB88 calculations.[27] The PBE-calculated Eads for an analogous (bridge) N2/γ-Mo2N(111) configuration was previously reported to be 2.27 eV,[8] indicating a 0.66 eV decrease in the strength of the attractive N2-MoN interactions for this arrangement when transitioning from γ-Mo2N(111) to δ-MoN(0001). Whereas transition from the γ- to δ-phase model requires doubling the sub-surface concentration of N present in sixfold coordination sites between each subsequent Mo layer within the bulk, the introduction of a single N adatom into a threefold hollow adjacent to the bridge-bound nitrogen molecule has been shown to decrease the same configuration’s (PBE-calculated) Eads(N2) by as much as 1.23 eV on the γ-phase model system.[8] Consistent with this stark contrast between the impact from introducing excess N as undercoordinated terminal adatoms versus placement within fully coordinated sub-surface sites, we find ∼0.6 additional electrons localized at analogous surface versus bulk N sites in our N/MoN models, such that nearby N adatoms would be expected to more detrimentally limit the capacity for terminal Mo → N2 electron donation requisite for stronger molecule–support interactions.

Adsorption of NH, NH2, and NH3 Molecules

Adsorption energies, heights, and charge-transfer characteristics are given in Table for NH/MoN, NH2/MoN, and NH3/MoN adsorption in different configurations. Trial configurations were, respectively, limited to linear, trigonal, and tetragonal molecular arrangements designed to maximize surface-H and H–H distances within the MoN-supported NH, NH2, and NH3 structures. Optimization of the NH/MoN trial structures produces four stable adsorption configurations depicted in Figure S2. Adsorption energies associated with these structures range from 6.5 to 8.2 eV, with the most stable configuration resulting from adsorption at the hcp site. General trends in both surface-to-molecule charge transfer and molecular adsorption energy are similar to those observed for the adsorption of N atoms, with decreased surface coordination leading to decreased charge transfer and lower adsorption energies. Similarly, the fcc configuration again exhibits lower adsorption energy but slightly greater charge transfer to the adsorbed molecule relative to the hcp structure. A site-by-site Bader analysis comparison of analogous NH/MoN and N/MoN structures implies net electron redistribution from the NH adsorbate to the support upon N–H bond formation at all applicable N/MoN configurations. Unlike N, NH adsorption is also stable at bridge sites.

Table 3

Adsorption Characteristics of NH, NH2, and NH3 on MoN

	Eads (eV)	d [h] (Å)	charge transfer (e–)
NH
bridge	7.48	2.03 [1.30]	0.80
fcc	7.88	2.10 [1.34]	0.83
hcp	8.19	2.10 [1.26]	0.80
top	6.51	1.76 [2.11]	0.51
NH2
fcc	4.72	2.19 [1.60]	0.48
hcp	4.53	2.19 [1.60]	0.49
NH3
fcc	0.89	2.55 [2.00]	–0.01
top	1.51	2.28 [2.35]	–0.09

All NH2/MoN trial configurations converged into one of the only two stable adsorption structures (see Figure S3). In this case, similar characteristics are exhibited at both hcp and fcc sites, with a slight preference for fcc configurations indicated by a slightly greater adsorption energy (4.72 vs 4.53 eV) and the tendency for both top and bridge site migration to fcc rather than hcp sites upon relaxation.

H addition required for H–NH bond formation again leads to further NH adsorbate-to-support electron redistribution when comparing the Bader analyses of analogous NH2/MoN and NH/MoN configurations. Unlike the first H, the addition of the second H significantly weakens the NH adsorbate–surface interactions when comparing NH2 and NH adsorption energies and adsorbate–surface separations at common sites.

Like NH2/MoN, we again find only two stable NH3/MoN adsorption structures after relaxing several trial configurations (see Figure S4). Unlike NH and NH2, which exhibit preferences for highly coordinated adsorption at threefold hollow sites, NH3 instead exhibits a converse preference for less coordinated adsorption at top sites (Eads = 1.51 eV), with all but fcc trial configurations migrating to top sites upon relaxation. The general molecule surface bond weakening and electron-transfer trends continue when adding the third H, with Bader analysis showing NH3 charge localization roughly equivalent to that in free molecules for both stable adsorption configurations. NH bond lengths (1.03 Å) and orientations (∼109° H–N–H bond angles) also closely match those for isolated NH3 in both cases, suggesting favorable conditions for product desorption from either site once formed. NH bonds, within the adsorbed ammonia molecules, show an affinity for azimuthal alignment toward the nearest neighbor Mo atoms, while also exhibiting a modest repulsion away from nearby sub-surface N atoms. These tendencies combine to result in a slight (∼7.2°) rotation of top site ammonia relative to the hexagonal pattern of Mo atoms, such that each NH bond orients with a slight preference toward adjacent fcc rather than hcp hollows.

Adsorption of H and H2

Adsorption energies, heights, and charge-transfer characteristics are given in Table for H/MoN and H2/MoN adsorption in different configurations. Atomic H stably binds in four distinct adsorption structures following relaxation of each of our trial configurations (see Figure S5). Like N, H atoms also bind most stably at threefold hollow sites. Unlike N, H atoms adsorb more favorably at the fcc site (Eads = 1.03 eV), which again results in greater electron transfer to the adatom relative to the hcp configuration (0.46 vs 0.40 e–, respectively). Atoms placed at bridge sites were again found to migrate toward an adjacent hollow upon structural relaxation, but, unlike N, become trapped in a “bridge-hcp” configuration rather than fully relocating to the hcp site. A preference for adsorption in threefold hollows and accompanying a withdrawal of ∼0.4 e– was also exhibited by H on the γ-Mo2N(111) model system.[8] Different Eads values reported for analogous adsorption configurations in our work and that for γ-Mo2N(111) result from differences in how adsorption energy is defined for Hads in the two studies. Changing our definition to that used in the previous work [Eads(H) = EMoN + EH – EH/MoN] leads to an adsorption energy more comparable to that reported on the γ-phase model (Eads = 4.29 eV vs 3.55 eV on Mo2N). As shown with the N2/MoN(0001) versus N2/Mo2N(111) comparison, increased Eads is anticipated when comparing our optB88 values to those calculated using the PBE functional due to the added inclusion of van der Waals corrections in our work.[27]

Table 4

Adsorption Characteristics of H and H2 on MoN

	Eads (eV)	d [h] (Å)	charge transfer (e–)
H
bridge-hcp	0.68	1.94 [1.27]	0.41
fcc	1.03	2.00 [0.59]	0.46
hcp	0.72	2.02 [0.55]	0.40
top	0.30	1.75 [1.03]	0.33
H2
fcc-v	0.06	3.60 [3.26]	0.02
hcp-v	0.06	3.62 [3.27]	0.02
top-v	0.05	2.79 [2.79]	0.02
top-h	0.82	1.86 [1.85]	0.13
top-h-R30°	0.78	1.85 [1.82]	0.11

Of the five stable H2/MoN configurations depicted in Figure S6, only two show evidence of significant H2–MoN interactions (top-h and top-h-R30). By contrast, all three vertical adsorption configurations show near-zero adsorption energies but nonetheless present trapping barriers sufficient to prevent relaxation into one of the horizontal arrangements upon optimization. In the more favorable top-h configuration, H2 binds 1.85 Å above the surface with an adsorption energy of 0.82 eV and results in a 0.14 Å elongation of the H–H bond after a transfer of 0.13 electrons from the surface to the molecule. Azimuthally rotating the top-h H2 molecule by 30° about the surface normal results in similar, albeit slightly weaker, adsorption characteristics (see top-h-R30 in Table /Figure S6). While not discussed below, a representative barrier for the dissociation of top-h-R30 H2 into two nearest neighbor fcc sites was found to be 0.05 eV, which we take to imply the facile dissociative adsorption of hydrogen on Mo-terminated δ-MoN(0001).

Ammonia Formation Mechanism

N2 Dissociation

NEB was used to track changes in total energy when forcing a transition from the most stable N2/MoN configuration (bridge) to a dissociated configuration consisting of two singly occupied threefold hollow sites immediately adjacent to that bridge site (hcp + fcc). One might intuitively expect this process follows a bond elongation deformation mechanism (see Figure ). In the IS, a nitrogen molecule adsorbed in a bridge configuration is about 20 percent activated for dissociation (inferred by elongation in N–N bond length compared to gas-phase value). The highest energy image generated along the path connecting the IS and FS endpoints has been labeled TS and is calculated to lie 0.52 eV above the energy associated with the IS configuration. This value compares quite well with N2 dissociation barriers calculated for both stepped Ru[4] and an analogous rhombic Mo4 ensemble modeled as an active site on γ-Mo2N(111).[8] Consistent with individual hcp versus fcc N/MoN adatom adsorption characteristics, the N2 molecule tilts such that N(1) repositions slightly closer to the surface than N(2) throughout the dissociation process. For clarity, Mo(1) bounds the hcp site and Mo(3) bounds the fcc site.

Unlike the Mo2N results, we note far smaller changes in the dihedral angle formed between the fcc and hcp threefold hollows contained within the Mo4 active site following dissociation (∼5° here vs ∼20° in the previous work). Decreased torsional deformation could reflect a difference in active site rigidity resulting from changes in the bulk crystal structure. Despite near-negligible changes in a torsion angle, we do note significant lateral rearrangements within the Mo4 ensemble during dissociation. Specifically, Mo(2) and Mo(4) contract ∼0.5 Å toward one another and Mo(1) and Mo(3) migrate ∼0.3 Å away from one another while the system traverses the IS–TS activation barrier. The reconstructed Mo4 ensemble, which alleviates much of the surface strain introduced by the IS molecular adsorption configuration, is largely preserved upon relaxation from the TS into the fully dissociated FS structure.

Electronic Structure of N2 during Dissociation

Bader charge differences tracking changes in electron accumulation at each of the four interacting molybdenums and both nitrogen atoms are provided in Table for the IS, TS, and FS N2 dissociation configurations described above. Transfer of electrons from the active site to the molecule is evident throughout the dissociation process. Beginning with the IS configuration, each N within the adsorbed molecule.

Table 5

Bader Charge Differences for Dinitrogen and Four Neighboring Mo Atoms in IS, TS, and FS Dissociation Configurationsa

	IS	TS	FS
Mo(1)	–0.26	–0.22	–0.37
Mo(2)	–0.24	–0.38	–0.44
Mo(3)	–0.27	–0.24	–0.36
Mo(4)	–0.24	–0.39	–0.45
N(1)	0.69	0.75	0.95
N(2)	0.62	0.72	0.96

Positive (negative) values signify electron accumulation (depletion) relative to the charge density localized on the same atom in the decoupled surface-molecule system.

M carries more than 0.6 excess electrons (relative to isolated N2), with the bulk of this excess charge balanced by ∼0.25 e– depletions from each of the four directly interacting Mo4 atoms (relative to pristine MoN). As the adsorbed configuration distorts, surface-to-molecule charge transfer results in increased electron accumulation on both N atoms. While approaching the TS, increased N2 electron accumulation is accompanied by anisotropically enhanced electron withdrawal from the Mo(2) and Mo(4) atoms, which both contract toward the mildly elongated molecule during this change in the reaction coordinate. After surpassing the TS, each N continues to sequester more charge away from the surrounding Mo4 ensemble until both adatoms carry ∼1 excess electron in the fully dissociated FS configuration. The majority of the additional charge-transfer noted while the system relaxes into the FS configuration arises from further electron withdrawal from the Mo(1) and Mo(3) atoms, which the N adatoms migrate toward following dissociation.

Figure a,b provides spatially resolved representations allowing for the direct visualization of charge accumulations and depletions described in the Bader analyses of the IS and FS configurations reported above. The overall increase in charge accumulation (depletion) surrounding the N (Mo) atoms upon adsorption and dissociation is readily apparent upon visual inspection. The general geometries and symmetries associated with the areas of enhanced/diminished charge density are consistent with those reported for analogous N2/Mo4 configurations on γ-Mo2N(111).[8] As such, we presume qualitatively similar charge-transfer characteristics, whereby Mo atoms positioned along (across) the N–N bond axis interact with a single N via σ-bonds (simultaneously interact with both N2 π-bonds). Interestingly, both structures (IS and FS) show the depletion of the electronic charge presumed to have been delocalized between Mo atoms prior to N2 adsorption and dissociation. While also present, similar effects are less pronounced for analogous configurations present on the N2/γ-Mo2N(111) model system.[8]

Figure 4

Electronic structure results: (a,b) charge density differences at 0.005 e/Bohr3 for N2 on a pristine MoN surface in IS and FS states. (c,d): DOS differences of N2 and four Mo atoms interacting with N2 molecules, respectively. The y-scale in (c) is 10 times larger than that in (d).

Figure c,d provides direct comparisons of the projected density of states (DOS) at each of the six directly interacting Mo and N atoms before and after dissociation by plotting DOS(FS) – DOS(IS) subtraction profiles for each atom common to both structures. For obvious reasons, dissociation-induced DOS changes are most apparent on the N atoms (note that the vertical axes in the Mo plots have been rescaled by a factor of ten to make the smaller changes more discernible for those sites). After dissociating, deep-lying states associated with N–N bonding interactions are completely lost (see negative density differences at energies below −6 eV) and replaced with new Mo-affiliated states located closer to the Fermi level. Increased N–Mo bonding interactions lead to mutually increased DOS near −5 eV at each of the six atoms after dissociation. Moreover, small hcp versus fcc related differences in the N(1) and N(2) plots correlate well with differences in the near-Fermi change in DOS at the correspondingly interacting Mo(1) and Mo(3) sites. General trends are again similar to those noted for Mo4-activated dissociation on γ-Mo2N(111);[8] however, N adatom-to-adatom asymmetries appear far less pronounced on our δ-MoN(0001) model, possibly reflecting a more homogeneous Mo-coordination to sub-surface N atoms in the δ-phase system.

In summary, N2 dissociation on Mo-terminated δ-MoN(0001) can be viewed as a two-step process. The process is first fully activated following Mo4 structural rearrangements associated with the donation of the electronic charge to the molecule via Mo(2,4) interactions across its π-bonds and concomitant N2 bond distortion. Once fully activated, the elongated molecule readily dissociates to relax into its FS structure, within which the resulting adatoms accumulate more electronic charge through σ-bond interactions with Mo(1,3) atoms aligned along the parent molecule’s primary axis. Despite the subtle differences in reaction-induced active-site reconstruction and qualitative electron redistribution characteristics, these results are otherwise consistent with those reported for N2 dissociation on Mo-terminated γ-Mo2N(111). This includes reaction energies, activation barriers, and net adsorbate-support charge transfers all falling within ∼10% of those in the former work (with all values skewing smaller in our work).[8] In addition to the differences in the bulk structure and N-content between previous studies and this, discrepancies in reaction energies and energy barriers can arise due to a lack of inclusion of van der Waals interactions (PBE vs optB88), which are shown to be important for the selection of suitable catalysts.[28,35]

NH Hydrogenation

To most efficiently draw comparisons to past works on γ-Mo2N(111), we have only considered hydrogenation pathways resulting from the sequential H adatom attack of isolated N(H) adsorbates. Nonetheless, we expect similar anisotropic impacts to overall Haber–Bosch catalysis in the presence of co-adsorbed N for our δ-phase model due to its similar N2/Mo4 active site geometry and interaction asymmetries. Using NEB, individual reaction barriers, energies, and transition states (TSs) were obtained for each of the three sequential hydrogenation iterations. IS configurations consist of NH and H adsorbates positioned in their most proximate combination of individually most stable adsorption geometries. On this basis, all IS configurations consist of NH species placed in threefold hollows (hcp for N and NH and fcc for NH2) immediately adjacent to an fcc-bound H adatom. Using the same logic, the most stable NH configurations were chosen as FS products for each elementary hydrogenation step. This overall reaction pathway is consistent with the more direct of the two ammonia formation mechanisms explored in the previous γ-Mo2N(111) study.[8]

Initial hydrogenation proceeds via simple H diffusion from an adjacent fcc- to top-site configuration (see Figure ). The energy cost associated with this process (Ea = 1.42 eV) is slightly larger than the difference in adsorption energy associated with top versus fcc H/MoN configurations on the pristine interface and again falls within ∼10% of the barrier found on γ-Mo2N for an analogous N hydrogenation mechanism.[8] Progressing from IS → TS → FS, N–H separations decrease from 3.32 → 1.82 → 1.03 Å, with the complete process resulting in an endothermic NH formation reaction energy of 0.58 eV.

Figure 5

Energy profile of N hydrogenation on the MoN surface. The N adatom is shown in dark blue, while the H atom is shown in red. Insets show top (upper) and side (lower) views of reaction.

As depicted in Figure , NH2 formation proceeds in a manner analogous to the first N hydrogenation step, with the barrier to reaction again stemming from the energy required for H diffusion from a hollow onto an adjacent Mo (Ea = 2.00 eV). In this case, the Mo–H TS bond tilts away from the adjacent NH, such that the H–NH distance only reduces from 3.30 → 2.70 Å when traversing the activation barrier. As the system relaxes out of the TS, H–H repulsion leads to the depicted NH2 FS configuration (θH–N–H = 101.6°), with the full hydrogenation step exhibiting an endothermic reaction energy of 1.68 eV. Unlike every prior step, the reaction barrier noted here is near-identical to (within 0.5%), while the NH–H TS geometry is significantly different from, that observed along an analogous reaction pathway in the γ-Mo2N(111) study.[8]

Figure 6

Energy profile of NH hydrogenation. The N adatom is shown in dark blue, while H atoms are shown in red. Insets show top (upper) and side (lower) views of the reaction.

Unlike N and NH hydrogenation, overcoming the final barrier to NH3 formation (Ea = 1.94 eV) includes the coordinated diffusion of both NH2 and H adsorbates from the nearest neighbor fcc sites onto a mutually common adjacent top site (see Figure ). During this process, the H and NH2 moieties contract about 1 Å toward one another while the latter rotates to maximize H–H–H separations in the TS geometry. Once formed, the TS readily transforms into the FS NH3/MoN configuration, with the full hydrogenation step again proceeding endothermically, with a reaction energy of 1.03 eV. In this case, a similar TS configuration was again noted on γ-Mo2N(111), but the reaction pathway leading to this structure is different due to a preference for NH2 adsorption at bridge sites within the previous model system.[8] The significantly larger barrier noted for this hydrogenation step on the γ-phase model (Ea = 2.56 eV) most likely at least partially results from this IS discrepancy between the two models. Like the previous study, progressing from IS → TS results in significant adsorbate-to-Mo4 electron transfer (not shown), with the diffusing H adatom contributing the majority of this charge while traversing the barrier to each sequential hydrogenation step.

Figure 7

Energy profile of NH2 hydrogenation. The N adatom is shown in dark blue, while H atoms are shown in red. Insets show top (upper) and side (lower) views of the reaction.

A full summary of the energy profile associated with the most energetically favorable pathway for nitrogen dissociation and subsequent stepwise hydrogenation is provided in Figure S7. While the focus of this work is specific to the Haber–Bosch reaction, results may also provide some insights into related nitrogen reduction reaction (NRR) mechanisms, which have become increasingly important within the electrochemical catalysis community. Impacts resulting from simultaneous interactions between multiple N and H atoms are likely to be important to elucidating the improved understandings of both reaction processes (Haber–Bosch and NRR), which we aim to address (at least in part) in a forthcoming publication.

Conclusions

Comparing our results for δ-MoN(0001)-activated N2 dissociation and subsequent hydrogenation to those from γ-Mo2N(111) models exposing equivalent rhombic Mo4 active sites suggests predominantly bulk-invariant surface chemistry at such sites.[8] Despite drastic changes to bulk stoichiometry and local (sub-surface) Mo–N coordination, Mo4 sites present at both interfaces are predicted to exhibit extremely similar N2 adsorption/dissociation and NH hydrogenation characteristics, including adsorption configurations, electron redistributions, and overall reaction pathways. For example, when treated at the same level of accuracy, N2 adsorption at the Mo4 site only weakens by 30% when changing the structure from γ-Mo2N(111) to δ-MoN(0001); a change that increases the number of N atoms directly coordinating to the Mo4 active site from five to eight. Notably, the placement of as few as only one N adatom above one of these three contrastingly populated hollows (omitting sites that become occupied following N2 dissociation) is predicted to reduce Eads(N2) by as much as ∼40% on the γ-phase model, with increased co-adsorption expected to only further weaken the interaction; a discrepancy highlighting the relatively inconsequential impact of bulk (vs surface) changes in the Mo/N stoichiometry/structure on Haber–Bosch catalysis. Similar to the γ-phase results, pristine Mo4 sites are highly active toward N2 adsorption and dissociation (Eads = 2.54 eV and Ea = 0.52 eV), but rate limited by large barriers to subsequent hydrogenation [Ea(NH + H → NH2) = 2.00 eV]. Likewise, hydrogenation barrier height reductions will be expected when increasing N coverage around the δ-MoN-contained Mo4 sites, with effects expected to vary in a manner showing anisotropic dependence on the binding site locations of the co-adsorbates for reasons analogous to some of those previously described for the γ-Mo2N interface.[8] As such, improved Haber–Bosch activity on the δ-MoN-contained Mo4 sites should be achievable by controlling the reaction conditions to beneficially tailor the [N] surface concentration. Moreover, results from this work further suggest that such adjustments can include radical changes to bulk MoN stoichiometry without introducing significant changes to Mo4 chemistry, so long as a hexagonal δ-MoN(0001)-like terminal structure is preserved.

Computational Details

Calculations were performed in the framework of the DFT[29] using plane wave pseudo-potentials implemented in VASP (Vienna Ab-initio Simulation Package 5.4.4). The generalized gradient approximation[30−32] was used with the optB88 package[33] to account for van der Waals interactions. This approach to modeling van der Waals effects has previously shown improved agreement between modeled and experimental parameters for establishing the activity of catalysts.[34,35] TSs were isolated using the NEB method to interrogate configurations and energies associated with three or five images found along pathways linking fixed IS and FS reaction configurations. Relaxed structures were optimized until Hellman–Feynman forces on all atoms were less than 0.01 eV/Å. An energy cutoff of 400 eV was used along with a Monkhorst–Pack grid[36] of 2 × 2 × 1. Lattice constants of bulk hexagonal MoN were determined using a Monkhorst–Pack grid of 18 × 18 × 18 with a cutoff energy of 800 eV. The calculated lattice constants (a = b = 2.88 Å and c = 2.86 Å) were found to be close to the experimentally determined values (a = b = 2.87 Å and c = 2.81 Å).[37]

A molybdenum nitride surface was modeled using a 5 × 5 super cell with five layers. The bottom two layers of the slab were fixed to their bulk positions and a vacuum of more than 10 Å was maintained between periodic images of the slab in the z-direction (perpendicular to the surface). Adsorption energies are defined aswith Esub representing the energy of the pristine substrate, Emol representing the energy of the molecule in the gas phase, and Emol–sub representing the energy of the surface after adsorbing the molecule. Using this convention, positive adsorption energies indicate thermodynamically favorable configurations (and vice-versa). For single atom adsorbates (N and H), half of the gas-phase energy associated with the corresponding diatomic species is used to calculate Emol. VASPKIT[38] and VESTA[39] were used for post-processing and visualization, respectively. A few key structures were selectively reoptimized using the PBE functional to provide direct comparisons to previous results reported at that level of accuracy for analogous γ-Mo2N(111)-supported configurations.[8]