Transition metal nitrides are widely studied due to their excellent mechanical properties, thermal stability, wear resistant properties, corrosion resistance, varying electrical and magnetic properties and their existing and emerging applications as hard wear-resistant coatings, diffusion barriers in microelectronics, protective decorative coatings, energy storage and spintronic materials. [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17] The best known transition metal nitrides are rock-salt structure group IV B and V B nitrides, with TiN being the most studied and having well established properties. 8,9,[18][19][20][21] Increasing the valence electron concentration by moving towards the right in the periodic table results in an increasing occupation of d-bands which increases the ductility and may therefore result in supertoughened alloys. [22][23][24] Correspondingly, rock-salt structure molybdenum nitride MoN with 3 electrons per formula unit occupying 4d orbitals has the potential to exhibit high ductility and toughness, and therefore represents a promising hard coating candidate. However, the high dorbital occupation is expected to also lead to thermodynamic or even mechanical instability of the rocksalt phase, similar to what has been reported for tungsten nitride, [25][26][27][28] which has the same number of valence electrons as MoN. Therefore, it is useful to study the phase stability in the Mo-N binary system, since the structure-composition-property relations of Mo1-xNx are not established, yet.
Molybdenum nitride has been reported to crystallize in a variety of phases depending on synthesis conditions including nitrogen partial pressure, temperature, substrate material, and energetics of depositing species. [29][30][31][32][33] The cubic γ-Mo2N is the most common experimentally observed phase of molybdenum nitride, 29,30,[34][35][36] while tetragonal β-Mo2N, 33,37,38 hexagonal δ-MoN 31,34,[39][40][41][42] and cubic γ-MoN [42][43][44] have also been reported, as well as a high pressure rhombohedral MoN2 phase at 3.5 GPa. 45 Theoretical studies based on density functional theory (DFT) calculations [46][47][48] predict that the most stable phases of the Mo-N system are hexagonal MoN2, 47 orthorhombic Mo4N3, 47 and hexagonal δ-MoN, 46 suggesting a proposed convex hull with three compound phases, namely δ-MoN2, ε-Mo4N3, and δ-MoN. 47 However, others predict lowenergy cubic structures for MoN0.44 and MoN0.69, 49 such that cubic (or tetragonal) Mo1-xNx phases are most stable for x < 0.42 while hexagonal structures are preferred at larger x, 50,51 suggesting that the convex hull for 0.31 < x < 0.5 is defined by β-MoN0. 44 and hexagonal δ-MoN. 48 These studies illustrate the diversity of possible structures in the Mo-N system, raise the question if other Mo1-xNx phases may be thermodynamically stable, and motivate our study on the convex hull of the Mo-N system.
In this paper, we report on the results from our search for the thermodynamically stable phases of the Mo-N system, using the evolutionary structure predictor USPEX [52][53][54] in combination with DFT calculations. We found eight stable Mo1-xNx compounds including newly discovered phases like monoclinic σ-MoN and σ-Mo2N3. Computation of the mechanical properties of the stable phases indicates that the elastic modulus increases with increasing nitrogen content. This is attributed to an increasingly covalent bonding character, leading to a ductile-to-brittle transition at x = 0.44 -0.5, as determined from Poisson's and Pugh's ratios. The previously reported experimentally synthesized Mo-N phases at various N/Mo ratios differ significantly from the thermodynamically stable phases predicted in this study. This may be due to a preference for phases with high symmetry, such as cubic or hexagonal, during a kinetically limited phase formation process.
Thermodynamically stable phases were identified by employing the evolutionary algorithm USPEX [52][53][54] coupled with first principles calculations done with the Vienna ab initio simulation package (VASP) 55,56 at zero temperature and pressure. As a first step, USPEX was employed for a variable composition search, using single molybdenum and nitrogen atoms as building blocks and varying both the composition and the number of formula units for a given composition. For this purpose, a maximum of 16 atoms within the primitive unit cell was used and the structure search was carried out until the predicted convex hull remained unchanged for twentyfive consecutive generations, resulting in a total of 953 distinctly different simulated structures with fully relaxed atomic positions as well as relaxed unit cell volumes and shapes. The limit of 16 atoms per unit cell is required due to the steep increase in computational cost with increasing unit cell size, but also limits the possible compositions that are explored with the USPEX variable composition search. Therefore, secondly, a large number of manually created structures were also explored by relaxing their unit cell and atomic positions with the same computational parameters as is done within the USPEX approach. In particular, (i) 20 structures that have a low energy within the W-N system 26 were explicitly investigated. (ii) Over 200 cubic structures with various compositions were formed by introducing cation and anion vacancies in stoichiometric NaCl structure γ-MoN, employing 8 different symmetric and asymmetric super cell sizes with 8 to 32 atomic sites, following the procedure described in Ref. 43. (iii) Over 100 structures with hexagonal symmetry containing 2-25 atoms were formed by stacking close-packed Mo and N planes, varying both the sequence of Mo and N planes as well as their stacking positions A, B, and C, while keeping the number of atoms per basal plane fixed at one. The purpose of the manually generated structures is twofold: They allow discovery of phases and compositions that contain more atoms in the primitive unit cell than the maximum number of atoms in the USPEX variable composition search. They also bridge a gap between the simulations and the experimentally most commonly observed cubic and hexagonal phases. More specifically, (ii) and (iii) allow to determine the cubic and hexagonal convex hull that can be directly compared to the overall convex hull, providing both a measure of the unfavorability of these phases with respect to the thermodynamically stable phases as well as insight into the competition between the hexagonal and cubic phases.
The DFT calculations were done using a plane wave basis set, periodic boundary conditions, the projector augmented wave method, 57 and the Perdew-Burke-Ernzerhof generalized gradient approximation exchange correlation functional. 58 The computational parameters are chosen such that all reported energies are computationally converged to <1 meV/atom. This includes a 500 eV energy cut-off for the plane-wave basis set expansion and a k-point mesh that varies with unit cell size from a total of 1728 to 8000 k-points, corresponding to a k-point resolution of 0.04 Å -1 in reciprocal space. Mo 4s, 4p, 4d, and 5s electrons are explicitly calculated, that is, they are not included in the core of the pseudo potential. This corresponds to a total of 14 valence electrons per Mo atom and 5 valence electrons per N atom. For all configurations, atomic positions, unit cell volume and shape are relaxed until an energy convergence of 10 -5 eV is reached. The convex hull is then constructed by identifying phases at different compositions 0 ≤ x ≤ 1 that are thermodynamically stable against dissociation into any other Mo1-xNx phase, where x = 0 and 1 correspond to bcc Mo and molecular N2. The latter was calculated using a 15×15×15 Å 3 unit cell containing a single relaxed N2 molecule, that is, the x = 1 phase in this study corresponds to vapor phase N2 rather than solid nitrogen, which is chosen because of the more direct relevance to actual experimental nitride synthesis. Phases with an energy less than 5 meV/atom above the convex hull are considered to lie on the convex hull. The crystal structures are visualized using VESTA. 59 The single crystal stiffness matrix is obtained from calculated stress-strain relationships as implemented in VASP, 60 by straining unit cell vectors by 0.015 Å in each direction required for a given crystal symmetry, and calculating the stress-strain matrix. The bulk modulus is obtained from the Murnaghan equation of state by fitting the calculated energy vs unit cell volume of the structures with a second order polynomial. The isotropic elastic modulus E and the isotropic shear modulus G are obtained from the stiffness tensor using Hill's criterion. 61 The isotropic Poisson's ratio ν and Pugh's ratio k are obtained from the elastic, shear and bulk moduli using the relations ν = E/2G -1 and k = G/B. 62 Vickers hardness H is calculated using Tian's model, 63 that is Hv = 0.92k 1.137 G 0.708 . We expect, based on previous studies that evaluate the accuracy of DFT material property predictions, [64][65][66][67] that our calculated atomic volume has an uncertainty of 3%, while the coefficients of the elastic tensor as well as E, G, and B have an uncertainty of ~5 % and the predicted H may have a ~10 % uncertainty based on the study by Tian. 63 Predicting the hardness using the hardness models from Refs. 68-70 indicates deviations by up to 15% in comparison to Tian's model, suggesting an approximately 15% uncertainty in our hardness predictions. The mechanical stability of the eight stable compound phases is verified by checking if their elastic constants satisfy the Born-Huang criteria for structural stability which are based on the individual elastic constants representative of the crystal symmetry. 71,72 For example, the Born-Huang criteria for hexagonal crystals is C11 > C12; 2C13 2 < C33(C11+C12); C44 > 0; C66 > 0. The dynamical stability of these crystals is analyzed by calculating the phonon frequencies with PHONOPY, 73 using conventional unit cells and checking for imaginary frequencies along high symmetry crystal directions. This approach is based on calculating a force constant matrix and therefore requires calculations with large surpercell sizes to achieve reliable phonon dispersion curves. However, the computational cost considerably limits the supercell size. For example, the conventional cubic unit cell of the γ-Mo14N11 phase already contains 50 atoms. Thus, simply doubling the cell size in three dimensions leads to a supercell with 400 atoms, which is already beyond what we can reasonably simulate for transition metal nitrides without compromising computational accuracy.
Figure 1 is a plot of calculated zero-temperature formation enthalpies Hf for fully relaxed Mo1-xNx compounds, plotted as a function of composition x = 0.20 -0.68. The Hf values are provided with standard states being bcc molybdenum and molecular nitrogen, such that Hf = 0 for both x = 0 and x = 1, corresponding to bulk Mo and molecular N2, respectively. The convex hull, denoted by the solid line, indicates the lowest energy for a given composition. Thus, any phase above the convex hull is thermodynamically unstable while phases on the convex hull, plotted as open black circles, are stable against thermodynamic dissociation. However, as mentioned in Section II, phases which lie within 5 meV/atom above the convex hull are considered to be thermodynamically stable, leading to a total of ten stable Mo1-xNx phases, including metallic Mo and molecular nitrogen, at compositions x = 0, 0.25, 0.4, 0.421, 0.429, 0.44, 0.5, 0.6, 0.667 and 1. They correspond to bcc-Mo, β-Mo3N, δ-Mo3N2, γ-Mo11N8, ε-Mo4N3, γ-Mo14N11, σ-MoN, σ-Mo2N3, δ-MoN2, and molecular nitrogen, where β, δ, γ, ε, and σ indicate tetragonal, hexagonal, cubic, orthorhombic and monoclinic phases, respectively. There are no stable phases for compositions 0 < x < 0.25 and 0.67 < x < 1. Thus, the figure shows the zoomed-in "interesting region" with x = 0.2 -0.68, while a plot of the full convex hull over the entire compositional range can be found in the supplementary document. Figure 1 also includes the convex hulls for cubic and hexagonal structures indicated by the red and green dashed lines. These curves are defined by the red squares and green triangles that denote the calculated Hf of cubic and hexagonal Mo1-xNx phases that lie on or near these convex hulls. The end points of the cubic and hexagonal convex hulls are set to the enthalpies of fcc and hexagonal Mo at x = 0, respectively, and to molecular nitrogen at x = 1. We note that these Hf values are determined by retaining bcc Mo and molecular N2 as zero-enthalpy reference states, even though the bcc Mo reference state is not included in the cubic and hexagonal hulls since they are limited for cubic structures to vacancy containing rocksalt MoN and for hexagonal structures to stacking of close-packed Mo and N planes, as described in more detail in Section II.
In the following, we briefly describe and discuss stable and competing Mo1-xNx phases, while a discussion of the overall Mo1-xNx convex hull, and the relationship to mechanical properties as a function of x will follow in the Discussion Section IV. Figure 2 illustrates the stable phases and Table 1 summarizes the key structural parameters. A detailed description of all phases, including complete unit cell vector and atomic position information is provided as Supplementary Material.
β-Mo3N (x = 0.25):
The Mo3N phase has a tetragonal structure with a conventional unit cell that contains 12 Mo and 4 N atoms and can be described as three conventional cubic rock-salt structure MoN unit cells that are stacked on top of each other and from which 8 nitrogen atoms have been removed, as illustrated in Fig. 2(a). Its formation enthalpy of -0.232 eV/atom lies 2 meV/atom above the calculated line of the convex hull connecting other phases, namely ε-Mo4N3 and bcc Mo, as shown in Fig. 1. Thus, there is no thermodynamic driving force for the formation of Mo3N from ε-Mo4N3 and bcc Mo. This may explain why there are no reports of the experimental observation of β-Mo3N, as discussed in more detail in Section IV. A competing cubic phase at this composition (x = 0.25) is 23 meV/atom higher than Hf of the stable β-Mo3N phase. This cubic γ-Mo3N has a similar atomic arrangement as the β-Mo3N shown in Fig. 2(a) and is obtained by relaxing both atomic positions and the volume of the unit cell, but forcing the c/a ratio to remain 3 so that the structure is cubic. This structure defines the cubic convex hull. That is, the line for the cubic convex hull plotted in Fig. 1 between x = 0.25 and 0.421 is defined by the calculated Hf of γ-Mo3N and γ-Mo11N8. The formation enthalpy of the competing hexagonal phase at x = 0.25 is 192 meV above Hf of the stable β-Mo3N phase and is actually higher than the hexagonal convex hull which is 126 meV above Hf of β-Mo3N and is defined by a line that connects the calculated formation enthalpies of δ-Mo2N at x = 0.333 and hexagonal Mo at x = 0.
The hexagonal δ-Mo3N2 phase consists of hexagonally close-packed Mo and N planes, which are stacked along the z-direction as illustrated in Fig. 2(b) and described in more detail in the supplementary document. Its formation enthalpy Hf = -0.371 eV per atom is 3 meV/atom above the convex hull. Therefore, δ-Mo3N2 is considered to be stable within our 5 meV tolerance criterion. The competing cubic γ-Mo3N2 phase has a formation enthalpy that is just 7 meV per atom larger than Hf of the hexagonal δ-Mo3N2. This cubic phase is formed by introducing 6% and 37% of cation and anion vacancies in rock-salt structure MoN, respectively, suggesting stabilization of the cubic phase by vacancies. 43 The small ΔHf = 7 meV between the cubic and hexagonal Mo3N2 structures may explain why the majority of experimental studies report a cubic γ-phase at this composition, as discussed in more detail in Section IV.
The cubic Mo11N8 phase has a primitive unit cell that contains 11 Mo and 8 N atoms and is derived from rock-salt structure stoichiometric MoN by removing one Mo and three N atoms from a supercell with 12 cation and 12 anion lattice sites, defined by the lattice vectors [ 1 1 1], [ 0 1 1 ] and [ 1 10 ] as illustrated in Fig. 2(c). Its formation enthalpy of Hf = -0.390 eV per atom is 4 meV/atom above the convex hull defined by ε-Mo4N3 and bcc Mo. The competing hexagonal phase at the same composition x = 0.421 has a formation enthalpy that is 6 meV per atom larger than that of the cubic phase, indicating that the cubic phase is thermodynamically more stable. This is in contrast to x = 0.4, as discussed above, where the hexagonal phase is 7 meV per atom more stable than the cubic phase. Thus, there is a predicted hexagonal-to-cubic transition at x = 0.41, as indicated by the cross-over of the red and green curves in Fig. 1.
The body-centered orthorhombic ε-Mo4N3 phase is illustrated in Fig. 2(d). It defines, together with bcc Mo, the line of the convex hull for x = 0 -0.429. Correspondingly, in the absence of any kinetic barriers, all MoxN1-x compounds with 0 < x < 0.429 should in principle phase segregate into bcc Mo and ε-Mo4N3, if following purely thermodynamic arguments at zero temperature and assuming zero uncertainty in our calculations. The convex hulls of competing cubic and hexagonal phases at x = 0.429 exhibit formation enthalpies that are 3 and 15 meV/atom higher than Hf of ε-Mo4N3. That is, the enthalpy of a cubic phase with x = 0.429 is very close to the overall convex hull and is therefore expected to have a considerable chance to be synthesized, while hexagonal phases are less likely to develop at x = 0.429. γ-Mo14N11 (x = 0.44): Fig. 2(e) illustrates the cubic γ-Mo14N11 phase which is formed by removing Mo and N atoms from stoichiometric rock-salt structure MoN, similar to the γ-Mo11N8 phase discussed above. This cubic phase is 16 meV/atom more stable than the competing hexagonal phase, continuing the trend of an increasing thermodynamic preference of cubic over hexagonal Mo1-xNx structures with increasing x in the compositional range 0.4 -0.45, as evident in Fig. 1. We note that the γ-Mo14N11 and the γ-Mo11N8 phases have similar lattice constants (0.5% difference) and are both derived from the same rock-salt MoN structure by the introduction of vacancies. Thus, we expect that cubic Mo1-xNx phases with any composition between these two phases (i.e. 0.421 < x < 0.44) can be formed without phase separation into γ-Mo14N11 and γ-Mo11N8. That is, the red-dashed line in Fig. 1 that connects the calculated Hf of γ-Mo14N11 and γ-Mo11N8 does not just represent the formation enthalpy of a phase separated (two cubic phases) compound, but is an estimate of Hf of a single γphase with a specific composition x. For example, at x = 0.429, this line indicates an Hf = -0.398 eV, which is just 3 meV/atom above Hf of the competing ε-Mo4N3 at the same composition, as discussed above. In fact, all calculated cubic structures for which Hf is plotted in Fig. 1 are derived with the same approach as the γ-Mo14N11 and γ-Mo11N8 phases. Thus, the entire cubic convex hull plotted as red dash-dotted line in Fig. 1 can be considered a continuous single phase field, as discussed in more detail in section IV.
The σ-MoN monoclinic phase illustrated in Fig. 2(f) is the most stable phase at a stoichiometric 1:1 composition ratio. The competing hexagonal and cubic phases for stoichiometric MoN have formation enthalpies that are 7 and 36 meV per atom higher than Hf for σ-MoN. The hexagonal phase has a four-atom unit cell where alternating close-packed Mo and N planes form a ABAC stacking, such that the Mo sublattice forms a simple hexagonal lattice and the N sublattice a close-packed hexagonal lattice which is shifted relative to the Mo sublattice to be also close-packed relative to the Mo lattice. The cubic MoN phase exhibits the NbO structure with 3 Mo and 3 N atoms per unit cell with a lattice constant of 4.11 Å. In contrast to our results, previous studies have predicted the hexagonal phase to be the most stable stoichiometric MoN phase. 47,50,51 We attribute this disagreement to the monoclinic σ-MoN phase exhibiting a relatively uncommon structure and the energy difference to the hexagonal phase to be relatively small, such that σ-MoN is not easily "discovered". Also, we are not aware of any experimental study that has reported σ-MoN. As discussed in Section IV, this may be related to the relatively small enthalpy differences to the hexagonal or cubic phases which have higher symmetry and correspondingly may be easier to synthesize. We note here that stoichiometric MoN in the rocksalt structure has a formation enthalpy which is 387 meV/atom higher than the monoclinic phase and is expected to be mechanically unstable, similar to WN in the rocksalt phase, 26,27 and is thus very unlikely to form.
σ-Mo2N3 shown in Fig. 2(g) is a base centered monoclinic phase. It is similar to what has been reported for W2N3. 26 The competing hexagonal δ-Mo2N3 phase has a formation enthalpy that is 48 meV/atom higher than Hf of the monoclinc σ-Mo2N3 phase, while the cubic hull at x = 0.6 is even higher, 102 meV/atom above σ-Mo2N3. That is, the monoclinic phase is considerably more stable than hexagonal or cubic structures. Nevertheless, we are not aware of any experimental study reporting σ-Mo2N3, but multiple studies suggesting that x > 0.5 leads to hexagonal structures, 31,45,47,50,51 as also discussed further in Section IV. δ-MoN2 (x = 0.667):
The hexagonal close packed δ-MoN2 phase is shown in Fig. 2(h). Its formation enthalpy is 4 meV/atom above the line of the convex hull, which is defined at x = 0.667 by the enthalpies of σ-Mo2N3 and molecular nitrogen. The δ-MoN2 forms N-Mo-N units similar to the S-Mo-S layers in MoS2. However, contrary to the weak Van der Waals bonding between neighboring S-Mo-S layers with a resulting large S-S distance of 3.1 Å, 74 the N-N bonds in δ-MoN2 are covalent, strong, and only 1.38 Å long, which is just 25% larger than the 1.10 Å bond length in molecular nitrogen, thereby yielding a predicted high mechanical strength as previously reported 47 and discussed in more detail below. A competing hexagonal phase with a simpler ABBABB stacking 47 has a 14 meV/atom higher formation enthalpy, while the competing cubic structure at x = 0.667 has a formation enthalpy that is 95 meV/atom higher than that of the hexagonal δ-MoN2 phase, consistent with the trend shown in Fig. 1, indicating that increasing x beyond 0.5 leads to an increasing thermodynamic stability of the hexagonal over the cubic Mo1-xNx phases. Another metastable structure with x = 0.667 is a rhombohedral structure that has been experimentally synthesized at a high pressure of 3.5 GPa. 45 Its formation enthalpy is 0.241 eV above that for hexagonal δ-MoN2, which is more stable than the rhombohedral phase at all pressures, as previously reported. 47 Therefore, we attribute the need for high pressure for the synthesis of rhombohedral MoN2 primarily to the reduction in the nitrogen chemical potential to favor nitride formation over molecular N2.
Figure 3 shows plots of the average atomic volume Vo and coordination numbers CMo and CN of Mo and N atoms in stable Mo1-xNx compounds as a function of composition x. The atomic volume is obtained by dividing the relaxed unit cell volume by the total number of atoms per unit cell, irrespective of their type (Mo or N), and is plotted in Fig. 3(a) as well as provided in Table 1. It decreases from Vo = 13.17 Å 3 for β-Mo3N to a plateau with Vo = 11.59±0.03 Å 3 for hexagonal δ-Mo3N2 and cubic γ-Mo11N8 and γ-Mo14N11, to a second plateau with Vo = 10.95±0.02 Å 3 for orthorhombic ε-Mo4N3 and monoclinic σ-MoN and σ-Mo2N3, and drops to Vo = 9.56 Å 3 for δ-MoN2. The data clearly shows an overall trend of a decreasing atomic volume with increasing x, which can be attributed to the considerably smaller covalent radius of 0.75 Å for N in comparison to 1.45 Å for Mo. However, interestingly, the Vo vs x data does not show a continuous trend but resembles a step function. The first plateau is defined by hexagonal and cubic structures with x = 0.4 -0.44 and has a 5.8% larger atomic volume than the second plateau defined by orthorhombic εand monoclinic σ-phases with x = 0.429 -0.6. Particularly interesting is that the composition ranges of these two plateaus overlap. More specifically, the composition of the orthorhombic ε-Mo4N3 is between the composition of the cubic γ-Mo11N8 and γ-Mo14N11 phases, suggesting that the lower Vo of ε-Mo4N3 is associated with the higher and more uniform number of nearest neighbors (CN = 6) that allows for a denser atomic packing in ε-Mo4N3, while the cubic phases in the vicinity of this orthorhombic phase are created by the introduction of vacancies, leading to a reduction in the number of nearest neighbors and correspondingly a higher Vo. Similarly, the fact that Vo remains constant between σ-MoN and σ-Mo2N3 despite the relatively large increase in x is attributed to a decrease in the coordination number of N atoms. This is because a decrease in CN (and conversely also an increase in CMo) results in a more covalent bonding character and an increase in bond-directionality, leading to more open structures. That is, the general trend of a decreasing Vo with increasing x is counteracted by a transition to a more open structure, leading to the observed plateau as also discussed in Section IV.
Figure 3(b) is a plot of the average coordination number CMo and CN of molybdenum and nitrogen atoms in the stable Mo1-xNx phases. For all MoxN1-x compounds with the exception of the most nitrogen rich δ-MoN2, the nearest neighbors are always of opposite kind, that is the nearest neighbors for Mo atoms are nitrogen atoms and vice versa, as evident from the interatomic distance data in Table 1. The most nitrogen rich δ-MoN2 phase represents a special case, because it exhibits N2 dimers with a 1.38 Å bond length. Thus, we define CN for this compound separately, counting both the N-neighbor as well as the three Mo atoms as nearest neighbors, yielding CN = 4.
The plotted number of nearest neighbors of nitrogen atoms is 6 for x ≤ 0.4 (β-Mo3N and δ-Mo3N2), as well as for ε-Mo4N3. Figure 3(b) is a plot of the average coordination number CMo and CN of molybdenum and nitrogen atoms in the stable Mo1-xNx phases. For all MoxN1-x compounds with the exception of the most nitrogen rich δ-MoN2, the nearest neighbors are always of opposite kind. The most nitrogen rich δ-MoN2 phase represents a special case, because it exhibits N2 dimers with a 1.38 Å bond length. Thus, we define CN for this compound separately, counting both the Nneighbor as well as the three Mo atoms as nearest neighbors, yielding CN = 4. The plotted number of nearest neighbors of nitrogen atoms is 6 for x ≤ 0.4, for β-Mo3N and δ-Mo3N2, as well as for ε-Mo4N3. CN decreases steeply with the transition to cubic structures above x = 0.4, to CN = 5.25 and 4.9 for γ-Mo11N8 and γ-Mo14N11, remains nearly constant between x = 0.44 -0.5 and decreases further from CN = 5 for σ-MoN to CN = 4 for σ-Mo2N3 and δ-MoN2. Conversely, the average coordination number of molybdenum atoms increases from 2 for β-Mo3N, to approximately 4 for δ-Mo3N2, γ-Mo11N8 and γ-Mo14N11, 4.5 for ε-Mo4N3, 5 for σ-MoN, and 6 for σ-Mo2N3 and δ-MoN2. That is, the CMo increases approximately linearly with increasing x, which is a direct consequence of the increasing N-to-Mo ratio. Correspondingly, CN decreases with increasing x. However, this decrease is not linear but follows an approximate step-function with CN decreasing from 6 to 5 around x = 0.43 and from CN = 5 to 4 for x between 0.5 -0.6. These variations in CMo and CN are responsible for the plateaus and variations in Vo, as presented above, and also cause a transition from ductile to brittle with an associated increase in the predicted hardness, as presented below and discussed in section IV. We note that the C vs x data do not show a perfectly monotonic trend, particularly, for x = 0.4 -0.44. This is due to the competition between the cubic γ-Mo11N8 and γ-Mo14N11 and the orthorhombic ε-Mo4N3 phases. These cubic phases exhibit vacancycontaining NaCl structures and therefore have considerably smaller coordination numbers than the orthorhombic phase. The latter has a perfect CN = 6 while the cubic phases have CN = 5.25 and 4.9 and, similarly, CMo = 4.5 for ε-Mo4N3 is nearly one neighbor larger than CMo = 3.81 and 3.85 for γ-Mo11N8 and γ-Mo14N11. The smaller number of nearest neighbors is a direct consequence of the vacancies which also result in the larger atomic volume presented above.
Figure 4 is a plot of the calculated mechanical properties of the eight stable Mo1-xNx compounds as a function of composition x, including bcc Mo (x = 0) but excluding molecular nitrogen (x = 1). This data is also listed in Table 2. The calculated bulk modulus B = 264 GPa for pure Mo (x = 0) is in agreement with the previously measured 262 -265 GPa for bcc Mo. [75][76][77] B increases approximately linearly with x to 293 GPa for β-Mo3N and 347 GPa for δ-MoN2, resulting in an overall increase of 18% in the composition range x = 0.25 -0.667. We attribute the increase in B primarily to the decreasing atomic volume, which decreases by 27% over the same composition range.
The isotropic shear modulus G is determined from the calculated elastic tensor, as described in Section II. It increases monotonically with x, from 126 GPa for bcc Mo to 212 GPa for MoN2. This corresponds to a 68% overall relative increase, which is considerably larger than the increase in B. We attribute this difference to the fact that the bulk modulus is controlled by the bond density (and strength) while the shear modulus is determined by the directionality of the bonds. 78 We also note that the largest increase in G occurs at x = 0.41, which corresponds to the transition between hexagonal δ-Mo3N2 and cubic γ-Mo11N8 and is associated with a steep drop in CN from 6 to 5. That is, the hexagonal to cubic transition leads to a reduction in the number of nearest neighbors which causes more directional covalent bonds and a corresponding increase in G of 22%. For the same transition, B increases only by 4% since Vo remains nearly constant (0.5% increase). We note that the G for the ε-Mo4N3 remains nearly as high (2% drop) as for γ-Mo11N8, even though CN = 6 for the orthorhombic ε-Mo4N3. We attribute this to the 6% lower atomic volume of the particularly dense orthorhombic structure, which also results in a 7% higher B, as evidenced by an outlier in the plot of B in Fig. 4(a) at x = 0.429. The isotropic elastic modulus E increases by an overall 62% from 327 GPa for Mo to 530 GPa for δ-MoN2. Similar to G, the elastic modulus exhibits the steepest increase at x = 0.41, which corresponds to the transition from hexagonal Mo3N2 to cubic Mo11N8, and is attributed to an increase in the covalent character of the Mo-N bonds.
Figure 4(a) also includes the hardness H, determined from the calculated B and G using Tian's model, 63 as described in Section II. It increases gently from 12 GPa for bcc Mo to 14 GPa for δ-Mo3N2, exhibits a steep increase to 20 GPa for γ-Mo11N8 and 24 GPa for σ-MoN and σ-Mo2N3, and finally reaches 26 GPa for δ-MoN2. The overall H vs x data resembles the composition dependence of E, showing a steep increase at the hexagonal-to-cubic transition at x = 0.41, followed by a second increase as the structure transitions from cubic to monoclinic between x = 0.44 and 0.5. The steep increase at x = 0.41, is primarily attributed to the decrease in CN and the corresponding increase in bond directionality, as presented above. The second steep increase at x = 0.44 -0.5 is attributed to an increase in CMo from 3.85 to 5, while CN remains constant. This results in an increasingly densely packed structure and a higher density of Mo-N bonds, leading to a higher shear modulus and hardness and transition from ductile to brittle, as discussed below. The calculated Poisson's ratio remains approximately constant at ν = 0.29±0.01 for x = 0-0.4, which includes metallic Mo as well as tetragonal β-Mo3N and hexagonal δ-Mo3N2. However, it then drops to ν = 0.27±0.01 for the cubic γ-Mo11N8 and γ-Mo14N11 and orthorhombic ε-Mo4N3 phases, and to ν = 0.23±0.01 for the monoclinic σ-MoN and σ-Mo2N3 phases, before increasing again to 0.25 for δ-MoN2. The latter value is in agreement with ν = 0.25 from a previous computational study. 47 We use the Poisson's ratio as an indicator of bond character and ductility. More specifically, a material can be classified as brittle if ν < 0.25, and conversely as ductile if ν > 0.25. 62 In addition, a decreasing Poisson's ratio suggests an increasing resistance against changes in bond angles, indicating an increasing covalent bonding character. Correspondingly, the steep decrease from ν = 0.29 for δ-Mo3N2 to ν = 0.24 for σ-MoN indicates an increasing tendency for brittle fracture associated with a transition to an increasingly covalent bonding character at x = 0.4 -0.5. Thus, Mo1-xNx is predicted to exhibit a ductile-to-brittle transition at x = 0.44 -0.5, as the lowest energy phase transitions from cubic γ-Mo14N11 to monoclinic σ-MoN. This transition is also evident from the steep increase in shear and elastic moduli at x = 0.4 -0.5, presented above. Another indicator for a ductile-to-brittle transition is Pugh's ratio k, where k < 0.6 is typically used as a condition for ductility while k > 0.6 leads to a brittle material. 62 Pugh's ratio is plotted in Fig. 4(b). It remains approximately constant at k = 0.48 -0.51 for x = 0 -0.4, increases steeply to k = 0.64 for σ-MoN and finally decreases slightly to 0.61 for δ-MoN2. That is, Pugh's ratio exhibits a composition dependence that is approximately inverse to that of the Poisson's ratio. The steep increase in Pugh's ratio at x = 0.4 -0.5 indicates an increase in the covalent bonding character with increasing x, which is also reflected by the steep increase in the elastic modulus from 387 GPa for δ-Mo3N2 to 487 GPa for σ-MoN, and in the shear modulus from 150 GPa for δ-Mo3N2 to 196 GPa for σ-MoN. Using the above k < 0.6 criterion for ductility suggests a ductile-to-brittle transition between x = 0.44 and 0.5, in perfect agreement with the Poisson's ratio argument. Thus, in summary, the stable phases of Mo1-xNx are predicted to be ductile for 0 ≤ x ≤ 0.44 and brittle for x ≥ 0.5, with a ductile-to-brittle transition between the cubic γ-Mo14N11 and the monoclinic σ-MoN phases at x = 0.44 -0.5.
The calculated stiffnes tensors are also used to confirm mechanical stability of the Mo1-xNx phases that are on or near (within 5 meV) the convex hull. All eight compound phases satisfy all Born-Huang mechanical stability criteria. 71,72 That is, as expected for phases that are thermodynamically stable, the presented Mo1-xNx compound phases exhibit mechanical stability. In addition, the dynamical stability has been explored by calculating the phonon dispersion curves and checking for imaginary frequencies. These calculations indicate that four phases, namely δ-Mo3N2, ε-Mo4N3, σ-MoN, and δ-MoN2 exhibit real frequencies throughout the entire Brillouin zone, while the remaining four structures β-Mo3N, γ-Mo11N8, γ-Mo14N11 and σ-Mo2N3 show imaginary frequencies in small fractions of their Brillouin zone. An example is given in the supplementary document. These imaginary frequencies may indicate dynamical instability which may be associated with slightly more stable structures that, however, could not be found with the performed phase-search due to the limited number of atoms. However, more likely, the few imaginary frequencies are a computational artifact associated with the limited size of the supercells used for these phonon calculations, that is, they are due to limited accuracies caused by the large computational cost of calculating the phonon dispersion curves of these relatively asymmetric structures, as also discussed in Section II.
The calculated Mo1-xNx convex hull presented in Fig. 1 indicates a diverse set of stable compound phases, and well defined structural transitions with increasing N content from bcc Mo for x = 0 to tetragonal β-Mo3N for x = 0.25, hexagonal δ-Mo3N2 for x = 0.4, cubic and orthorhombic phases for x = 0.421 -0.44, monoclinic σ-phases for x = 0.5 -0.6, and hexagonal MoN2 for x = 0.667. As a first discussion point, we note that the calculated convex hull for 0 ≤ x ≤ 0.429 is a straight line which, within uncertainty, even remains straight to γ-Mo14N11. That is, the convex hull is not particularly "convex," which means that the six listed phases for x ≤ 0.44 effectively can co-exist in thermodynamic equilibrium, assuming contributions from the entropy to the free energy are negligible. Consequently, there is a negligible driving force for, for example, a mixture of bcc Mo and ε-Mo4N3 to form a β-Mo3N or a δ-Mo3N2 compound, and similarly, these latter compounds also have a negligible driving force for dissociation into neighboring compounds with higher and lower N-concentrations. As a consequence, a Mo1-xNx compound with 0 ≤ x ≤ 0.44 can exhibit multiple phases and/or spatial composition variations without increasing its enthalpy. This suggests, that kinetic considerations during synthesis will determine the developing structure as well as the composition of Mo1-xNx. That prediction is consistent with the diverse structures that have been reported from studies of Mo1-xNx thin film growth as well as prior theoretical investigations, as summarized in Section I.
The region of the convex hull between x = 0 and 0.4 exhibits only one phase: tetragonal β-Mo3N at x = 0.25. This phase has an Hf that is just 2 meV above the line connecting bcc Mo at x = 0 and ε-Mo4N3 at x = 0.429 and is therefore considered to lie on the convex hull, however, without a thermodynamic driving force for its formation from bcc Mo and ε-Mo4N3. Previous experimental reports indicate that compositions with x = 0 -0.33 lead to phase separation into bcc Mo and β-Mo2N 33,37,38,79 or γ-Mo2N 29,30,[34][35][36] with the γ-Mo2N being more commonly observed, while computational reports indicate variants of tetragonal β-Mo2N with nitrogen vacancies 50 and a lowenthalpy orthorhombic ε-Mo2N. 47 Our calculations of these phases with x = 0.33 show, however, that both ε-Mo2N with Hf = -0.304 eV and β-Mo2N with Hf = -0.293 eV are 7 and 18 meV/atom above the convex hull, as also indicated in Fig. 1, suggesting that they are not stable against dissociation into ε-Mo4N3 and bcc Mo or β-Mo3N. Nevertheless, we are not aware of any experimental report on the synthesis of β-Mo3N, which may be due to (i) the relatively low symmetry of this phase which limits its nucleation, (ii) the competing phase transition directly to bcc Mo and ε-Mo4N3, or (iii) some experimental ambiguity in the phase determination, as different cubic, tetragonal, and orthorhombic phases or phase mixtures exhibit similar diffraction pattern signatures. The predicted β-Mo3N phase exhibits a primarily metallic bonding character, which results in a relatively high Poisson's ratio ν = 0.28 and a low Pugh's ratios k = 0.51. It has a volume per atom Vo = 13.17 Å 3 which is high in comparison to other Mo1-xNx phases and is attributed to the high Mo-to-N ratio that also results in the high nitrogen and low molybdenum coordination numbers of CN = 6 and CMo = 2, respectively. The hexagonal and cubic convex hulls, indicated by the green and red curves in Fig. 1 suggest formation enthalpies at x = 0.25 that are 23 and 126 meV/atom higher than the tetragonal β-Mo3N phase, respectively. These convex hulls are 18 and 60 meV/atom above the total convex hull for at x = 0.3, respectively. The relatively high formation enthalpies of cubic and hexagonal phases in this composition region along with the negligible thermodynamic driving force for the formation of β-Mo3N from ε-Mo4N3 and bcc Mo may explain the lack of experimental observations of molybdenum rich MoxN1-x with x = 0 -0.33. In contrast, experimental studies report either tetragonal β-Mo2N 33,37,38,79 or cubic γ-Mo2N. 29,30,34,[36][37][38] For 0.33 < x < 0.4, both the cubic and hexagonal hulls approach the overall stable hull with increasing x. More specifically, the cubic hull is 14 meV above the convex hull for x = 0.33, and this difference decreases to 9 and 6 meV for x = 0.375 and x = 0.4. Similarly, the Hf of calculated hexagonal phases decreases from 15 meV above the convex hull at x = 0.333 to 8 and 3 meV above the convex hull for x = 0.375 and x = 0.4. This indicates that both cubic and hexagonal phases are expected to have an increasing likelihood for formation with increasing x in this composition range, because of the decreasing enthalpy difference ΔHf for forming these phases instead of the stable phase defined by the overall convex hull. Comparing the cubic and hexagonal hulls shows that cubic phases are clearly preferred over the hexagonal phases for x < 0.3, while hexagonal and cubic phases have similar enthalpies for x = 0.3 -0.4. The difference in Hf is less than 7 meV for x = 0.33 -0.4, such that the competition between hexagonal and cubic phases for x = 0.33 -0.4 is likely determined by kinetic factors. The transition from β-Mo3N with x = 0.25 to δ-Mo3N2 with x = 0.4 leads to a denser packed structure with a 12% lower atomic volume and an increase in CMo from 2 for β-Mo3N to 4 for δ-Mo3N2 which results in a more covalent bonding character. Experimental studies have reported variants of γ-Mo2N 29,32,80 and β-Mo2N 33,35,37,81 for this composition range. However, both of these phases have calculated Hf values that are higher than that of orthorhombic ε-Mo2N which in itself lies above the convex hull, as indicated in Fig. 1 for x = 0.33. In spite of the hexagonal δ-Mo2N phase having a formation enthalpy that is just 1 meV higher than the cubic γ-Mo2N, we are not aware of any experimental reports for the synthesis of δ-Mo2N. This further illustrates the fact that kinetic constraints strongly affect Mo1-xNx synthesis. We note that particularly the orthorhombic ε-Mo2N and ε-Mo4N3 phases have lower symmetries than the experimentally observed tetragonal β-Mo2N, which may suppress nucleation of the εphases, leading to metastable β-Mo2N or even less stable (but more symmetric) cubic γ-Mo2N.
The elastic modulus and hardness moderately increase with x in the compositional range x = 0 -0.4, with E increasing from 327 GPa for bcc Mo to 385 and 387 GPa for β-Mo3N and δ-Mo3N2, while H increases from 12 GPa for bcc Mo to 15 and 14 GPa for β-Mo3N and δ-Mo3N2. These values are consistent with the experimentally reported mechanical properties in this composition range which are, however, not from the equilibrium β-Mo3N and δ-Mo3N2 phases but instead from β-Mo2N and γ-Mo2N, exhibiting a wide range of measured values with E = 130 -420 GPa, 36,50,[80][81][82] and H = 6 -33 GPa. 36,50,[80][81][82] These wide ranges can be attributed to a variety of reasons including residual stress in thin films and variations in composition, density and grain size.
In the compositional range x = 0.4 -0.5, cubic and hexagonal Mo1-xNx phases are in close competition, as evident from the nearly overlapping cubic and hexagonal convex hulls shown in Fig. 1. The cubic convex hull has its minimum of Hf = -0.416 eV at x = 0.48, while the minimum of the hexagonal hull of Hf = -0.422 eV is at x = 0.5. However, at this stoichiometric composition, the monoclinic σ-MoN phase is even slightly (7 meV) more stable than the hexagonal δ-MoN. In addition, the orthorhombic ε-Mo4N3 phase has an Hf that is slightly (3 meV) below the cubic hull defined by γ-Mo11N8 and γ-Mo14N11 phases. That is, while hexagonal and cubic phases dominate the interesting composition range x = 0.4 -0.5, lower symmetry ε-Mo4N3 and σ-MoN phases at x = 0.429 and 0.5 contribute to the possible phase diversity. Based on the plotted cubic and hexagonal convex hulls, compositions between x = 0.4 and 0.421 should in principle phase segregate into hexagonal δ-Mo3N2 and cubic γ-Mo11N8. However, the small enthalpy difference of the competing δ and γ phases suggests that a single phase, either δ or γ should be expected, while the complete convex hull would even suggest the unlikely phase segregation into bcc Mo and orthorhombic ε-Mo4N3. A similar argument can be made for the x = 0.44 -0.5 range, for which Mo1-xNx in principle should segregate into γ-Mo14N11 and σ-MoN. However, the enthalpies of the cubic phases are so close to the convex hull, for example, just 6 meV above the calculated hull for x = 0.48, such that we expect cubic phases to dominate at x = 0.44-0.49, while x = 0.5 may lead to the monoclinic or hexagonal phase since Hf of the cubic phase increases steeply with increasing x ≥ 0.48. More generally, Hf of both δ and γ phases are so close (0-16 and 0-7 meV, respectively) to the overall convex hull for the entire range x = 0.4 -0.49, that we expect kinetic barriers during compound formation to determine which phase will ultimately form. This is consistent with reports from experimental studies, indicating primarily δ-MoN 32,34,42,79 and γ-MoN in the vicinity of x = 0.5. 32,43,83,84 Our results refine previous theoretical predictions which report a region of γ-phases for x = 0.375 -0.44, 50 and a region of δ-phases for x = 0.44 -0.5. 46,50,85 We attribute the different results in our study primarily to the much more extensive phase search, both in compositional and symmetry space, as the cubic and hexagonal phases that define the convex hull have relatively large unit cell sizes while the ε-Mo4N3 and σ-MoN phases that also contribute to the convex hull have low symmetries and are correspondingly difficult to discover.
The composition range x = 0.4 -0.5 is also interesting with regards to changes in bond character and mechanical properties. The atomic volume remains roughly constant with Vo = 11.59 ± 0.03 Å 3 for δ-Mo3N2, γ-Mo11N8 and γ-Mo14N11 but decreases to 10.97 ± 0.02 Å 3 for ε-Mo4N3 and σ-MoN. This steep decrease is explained by the molybdenum coordination number which remains nearly constant at approximately 4 for δ-Mo3N2, γ-Mo11N8 and γ-Mo14N11, but increases to CMo = 4.5 and 5 for ε-Mo4N3 and σ-MoN, leading to more densely packed structures. A slight increase in x from 0.4 for δ-Mo3N2 to 0.421 for γ-Mo11N8 leads to steep increases in the elastic and shear moduli from 387 -462 GPa and 150 -183 GPa, respectively. These increases occur despite
Vo remaining constant and are attributed to the decrease in CN from 6 to 5.25 which promotes a more covalent character and more directional bonding. A second increase in the stiffness of MoxN1-x is predicted between x = 0.44 and 0.5, corresponding to the transition from the γ-Mo14N11 to the σ-MoN phase. This is facilitated by an increase in CMo from 3.81 to 5 which suggests an increasingly covalent bonding character, since an increase in (covalent) bonds to nearest neighbor nitrogen atoms replaces the strength of metallic 2 nd -nearest-neighbor Mo-Mo bonds, leading to a drop in Poisson's ratio from 0.27 to 0.24 and an increase in Pugh's ratio from 0.55 to 0.64. These changes suggest a ductile-to-brittle transition at x = 0.44 -0.5, which is accompanied by a 10% increase in the elastic modulus and 33% increase in calculated hardness. In total, E increases by 26% and H increases by 71% over the compositional range x = 0.4 -0.5.
For x = 0.5, experimental studies report the formation of δ-MoN 32,34,42,79 and the γ-MoN 32,43,83,84 phases. Our calculations of the δ-MoN and γ-MoN (in the NbO phase) phases indicate formation enthalpies which are 7 and 36 meV higher than the most stable σ-MoN phase. These enthalpy differences are small, suggesting that kinetic barriers likely play an important role in determining the stable phase at x = 0.5, which might lead to the synthesis of a more symmetric hexagonal structure over a relatively uncommon monoclinic structure. We are not aware of any experimental observation of the σ-MoN phase, which may be attributed to its low symmetry that makes synthesis challenging and also to kinetic barriers associated with the significant atomic rearrangements required for the transition from symmetric γ-MoxN1-x at x = 0.44 to σ-MoN at x = 0.5, such that a more symmetric δ-MoN phase is experimentally synthesized. 34,79,82,86 The reported mechanical properties of the experimentally synthesized δ-MoN phase and its variants exhibit a wide range, with E = 270 -420 GPa 80,82 and H = 15 -27 GPa, 80,82 in reasonable agreement with our predicted E = 487 GPa and H = 24 GPa.
The compositional range x = 0.5 -0.6 of the MoxN1-x convex hull is populated by the monoclinic σ-MoN and σ-Mo2N3 phases. The atomic volume remains constant between x = 0.5 -0.6, which is attributed to competing effects from CMo which increases from 5 to 6 and from CN which decreases from 5 to 4. Similarly, the mechanical properties including moduli and Pugh's and Poisson's ratios are nearly identical for σ-MoN and σ-Mo2N3. Previous studies report hexagonal phases for 0.5 < x < 0.6, including the experimental synthesis of δ-Mo5N6 with x = 0.54, 31 and theoretical predictions of a hexagonal symmetry. 47 This is in direct contradiction to our results which show that cubic and hexagonal convex hulls deviate significantly from the convex hull of the stable σ-phases in this region, indicating that the monoclinic phases are clearly more stable. In particular, the cubic convex hull is ΔHf = 36 meV above the overall convex hull at x = 0.5 and this difference increases to 109 meV for x = 0.6. Similarly, ΔHf for the hexagonal hull increases from 7 meV at x = 0.5 to 45 meV at x = 0.6. Therefore, increasing x from 0.5 to 0.6 results in an increasing preference of the monoclinic phases over hexagonal and (even more so) cubic phases, such that σ-Mo2N3 is thermodynamically much more stable than δ-Mo2N3, while a hypothetical cubic γ-Mo2N3 would be even less stable than the hexagonal phase. However, we are not aware of any experimental reports of the σ-Mo2N3 phase yet, which may again be due to kinetic barriers which favor the more symmetric hexagonal phase over the monoclinic phase.
The most nitrogen-rich compound of the Mo1-xNx convex hull is hexagonal δ-MoN2. Our search algorithm could not find any other thermodynamically stable phases with higher nitrogen content, i.e., with 0.67 < x < 1. The δ-MoN2 phase is rather different from the other phases predicted for x = 0 -0.6, as it includes a relatively short 1.38 Å nearest neighbor bond between two nitrogen atoms, while the nearest neighbor atoms for all other structures (x ≤ 0.6) are always of opposite type. The nitrogen atoms form dumbbells which also result in a dense structure as evident from the 13% decrease in Vo during the transition from σ-Mo2N3 with x = 0.6 to δ-MoN2 with x = 0.667. These two phases have the same coordination numbers CMo = 6 and CN = 4, where the Mo atoms in both phases form 6 bonds with nearest neighbor nitrogen atoms, while the bonding environment of the N atoms changes from four Mo neighbors for σ-Mo2N3 to three Mo and one N neighbors for δ-MoN2. The shear modulus G = 212 GPa and hardness H = 26 GPa of hexagonal δ-MoN2 is approximately 8% above the values for σ-Mo2N3, which we attribute to the N2 dumbbells as well as possibly the 13% smaller atomic volume. Similar transition metal di-nitrides such as WN2, 87 TiN2 88 and CrN2 16 have previously been predicted to exhibit high mechanical strength. In contrast, experimental studies for Mo1-xNx with x > 0.6 report phases with hexagonal symmetry such as rhombohedral MoN2, which is different from the predicted hexagonal symmetry of MoN2, 47 WN2, 87 TiN2 88 and CrN2, 16 at x = 0.667. 45 We calculate a formation enthalpy for this rhombohedral phase that is 241 meV higher than that predicted for δ-MoN2, in good agreement with previously reported calculations that predict a 221 meV difference. 47 This large energy difference raises the question why the rhombohedral instead of the δ-MoN2 phase forms. We speculate that this is likely associated with kinetic barriers for the inclusion of the N2 dumbbells within close-packed Mo planes, or conversely, the expected low kinetic barrier for an N2 dumbbell to detach and desorb into the vapor phase, providing a considerable entropy benefit at finite temperatures, consistent with the reported tendency for a decreasing N-to-metal ratio in transition metal nitrides with increasing synthesis temperatures, including TiN, 89 HfN, 90 TaN, 91 WN, 25 MoN, 43 CrN, 92 NbN 4 and Ti1-xAlxN. 93
In summary, the convex hull of MoxN1-x contains ten stable phases which are bcc-Mo, β-Mo3N, δ-Mo3N2, γ-Mo11N8, ε-Mo4N3, γ-Mo14N11, σ-MoN, σ-Mo2N3, δ-MoN2, and molecular nitrogen. The convex hull is a straight line for 0 ≤ x ≤ 0.44, such that the first six phases in the above list can effectively co-exist in thermodynamic equilibrium, which may explain why some of these phases have never been synthesized yet. The calculations also predict new unexplored monoclinic σ-MoN and σ-Mo2N3 phases for x = 0.5 -0.6, and confirm the previous theoretical prediction of a hexagonal phase with N2 dumbbells for x = 0.667. In addition, we report the convex hull of cubic and hexagonal structures of MoxN1-x, since these high symmetry structures are often preferred in experimental growth conditions. The calculated formation enthalpies indicate that the cubic structures are clearly preferred over hexagonal structures in the molybdenum rich region (x = 0 -0.3) and hexagonal structures are clearly preferred in the nitrogen rich region (x > 0.5), while the cubic and hexagonal structures compete with each other closely in the range x = 0.3 -0.5 where their difference in enthalpy is less than 16 meV/atom, which implies that kinetic barriers during synthesis play an important role in the formation of either the hexagonal or cubic phase.
The average coordination number of nitrogen atoms CN for the stable structures in the convex hull decreases from 6 to 4, while the average coordination number of molybdenum atoms CMo increases from 2 to 6 as x increases from 0.25 to 0.667, leading to an increasing covalent bonding character with increasing x. The elastic modulus of the stable phases of MoxN1-x increases from 327 to 530 GPa and the shear modulus increases from 126 to 212 GPa as x increases from 0 to 0.667. Particularly steep increases in the moduli are predicted for x = 0.4 -0.5, which correspond to the composition range where phase transitions from hexagonal to cubic to monoclinic structures. The predicted hardness increases over the same composition range from 14 to 24 GPa, while a drop in Poisson's ratio and an increase in Pugh's ratio indicate a ductile to brittle transition between x = 0.44 and 0.5.
See supplementary material for (1) a figure of the complete Mo1-xNx convex hull, (2) detailed descriptions of the structure of the eight stable compound phases including lattice vectors and atomic positions, and (3) example phonon dispersion curves. Figure 1: Calculated formation enthalpy Hf per atom of Mo1-xNx compounds vs composition x. The labeled circles denote stable phases and the solid black line shows the convex hull. Open red squares and green triangles indicate enthalpies of cubic γand hexagonal δ-phases, respectively, while the red and green dashed lines show the corresponding convex hulls. The plot also includes tetragonal β-Mo2N and orthorhombic ε-Mo2N. A plot of the complete hull for x = 0 -1 is provided as supplementary document. Phase x Hf (eV) a, b, c (Å) α, β, γ (deg.) dMo-N (Å) dMo -Mo (Å) dN -N (Å) CMo CN β-Mo3N 0.25 -0 Table 2: Predicted atomic volume Vo, bulk modulus B, isotropic elastic modulus E, isotropic shear modulus G, hardness H, Poisson's ratio ν, and Pugh's ratio k of stable Mo1-xNx phases.