Stability, Elastic and Electronic Properties of Ta2N by First-Principles Calculations

Owing to exploring the influence of the N atoms ordering in Ta2N compounds on their properties, the stability, elastic, and electronic properties of Ta2N compounds (Ta2N-I: P3ml and Ta2N-II: P31m) were investigated using first-principles calculations based on density functional theory. Ta2N-II is energetically favorable according to the enthalpy of formation. Elastic constants were employed to reveal the stronger resistance to deformation, but weaker anisotropy, in Ta2N-II. A ductile-brittle transition was found between Ta2N-I (ductile) and Ta2N-II (brittle). The partial density of states showed a stronger orbital hybridization of Ta-d and N-p in Ta2N-II, resulting in stronger covalent bonding. The charge density difference illustrated the interaction of the Ta-N bond and electron distribution of Ta2N.


Introduction
The transition metal nitrides (TMNs) have widespread applications and motivate intense research because of their outstanding properties, such as good mechanical, thermal stability, electronic, and optical properties [1][2][3][4][5][6]. Tantalum nitrides (TaN x ) possess a series of equilibrium and metastable structures with various excellent properties, showing a remarkable richness of applications, including hard coatings for cutting tools, thinfilm resistors, and diffusion barriers in integrated circuits [7][8][9][10][11][12]. The change of growth conditions generates different phase and chemical compositions of polycrystalline TaN x , whose structures and properties vary greatly. Experimentally [7,[13][14][15][16][17] (No. 193, P6 3 /mcm), and η-Ta 2 N 3 (No.62, Pbnm) etc., have been synthesized by controlling growth conditions, and have attracted further study on their properties. Limited by the imperfection of X-ray diffraction (XRD), the structural difference (slight shift of N) between ε-TaN and π-TaN had not been distinguished until the application of neutron diffraction [18]. After detailed investigation, a similar phenomenon was found in Ta 2 N. In reality, there are two structures in Ta 2 N compounds (which are called Ta 2 N-I and Ta 2 N-II for convenience in this work). Terao et al. [14] detected hexagonal close-packed γ-Ta 2 N (Ta 2 N-I), also called β-Ta 2 N, with lattice parameters a = 3.05 Å, c = 4.92 Å (Space group, No.164) through XRD. Moreover, they also discovered the other structure of Ta 2 N (Ta 2 N-II) with lattice parameters a = 5.28 Å, c = 4.92 Å (Space group, No.162), due to a different ordering of N atoms, characterized by electron diffraction, which can display many weak peaks invisible in XRD. As the characterization techniques developed, neutron diffraction was employed by Conroy and Christensen to further confirm this structure (Ta 2 N-II) [19]. It is because of the ordering of vacancies that Ta 2 N-II come into being, which was presented by Friedrich et al. [20].
The same structure of Ta 2 N-II was also found in Cr 2 N [21], V 2 N [22], and Nb 2 N [23] compounds. By comparing the two structures of Ta 2 N, it is found that partial N atoms migrate to the gap between two Ta atom layers and leave corresponding vacancies in the N atom layer. Thus, the atomic stacking model changes from N-Ta-Ta-N-Ta-Ta-··· of Ta 2 N-I to N-Ta-N-Ta-··· of Ta 2 N-II. To sum up, the shift of N atoms generates the vacancies of the initial N atom layers in Ta 2 N-I, making a contribution to the formation of Ta 2 N-II.
Compared to Ta 2 N-II, Ta 2 N-I was very exciting for researchers. Cu diffusion studies [24] applied using density functional theory (DFT) showed that Ta 2 N-I can effectively prevent undesired Cu diffusion as an effective diffusion barrier material. The electronic and elastic properties of Ta 2 N-I were explored using first-principles calculation by Yu et al. [25]. Hence, it is significant to explore the properties of Ta 2 N-II, and clarify differences between Ta 2 N-I and Ta 2 N-II. In this work, the stability, elastic, and electronic properties of Ta 2 N-I and Ta 2 N-II were illustrated thoroughly using first-principles calculations based on DFT. The enthalpies of formation were calculated for stability comparison. Elastic constants were computed using a high-efficiency stress-strain method for elastic studies. Partial density of states (PDOS) and charge density difference (CDD) were applied for electronic research. Relevant results can promote the scientific study and commercial application of Ta 2 N.

First-Principles Calculations
The crystal structure information of Ta 2 N-I and Ta 2 N-II with space groups P3m1 and P31m, respectively, are shown in Table 1.  Figure 1 illustrates clearly the difference of crystal structure from different viewpoints. First-principles calculations were used for the Ta 2 N-I and Ta 2 N-II using the Vienna Ab initio Simulation Package (VASP) 5.4.1 version [26,27]. The generalized gradient approximation (GGA) parameterized by Perdew-Burke-Ernzerhof (PBE) [28] and the projector augmented wave (PAW) [29] method, with an energy cut-off of 600 eV, were employed for the exchangecorrelation effects and ion-electron interactions, respectively. The k-point meshes were Γ-centered 17 × 17 × 9, 9 × 9 × 11, 15 × 15 × 15 for the bulk Ta 2 N-I, Ta 2 N-II, bcc Ta; and 1 × 1 × 1 for the molecule (N 2 ), respectively. Herein, total energy calculations were performed in a large supercell of 14 Å × 15 Å × 16 Å for the N 2 molecule. The tetrahedron method with Blöchl corrections [30] was adopted for accurate total energy and electronic structure in integration of reciprocal space. The selected valence electrons were 5p 6 5d 6 6s 2 for Ta, and 2s 2 2p 3 for N. All systems were calculated with the convergence criteria 0.02 eV/Å for forces on ions, and 0.01 meV/atom for electrons, respectively.
In order to obtain the equilibrium structure, the four-parameter Birch-Murnaghan equation of state (EOS) was applied to fit the energy vs. volume (E-V) data points [31,32] where a, b, c, and d are fitting parameters. This EOS can estimate the equilibrium total energy E 0 , volume V 0 , bulk modulus B 0 , and its pressure derivative B 0 '. To compare Ta 2 N-I and Ta 2 N-II in terms of their thermodynamic stability, enthalpy of formation at 0 K, f E (kJ/mol), was calculated as: where ρ(Ta 2 Nself-consistent) is the total charge densities obtained by self-consistent calculation. The superposition of the atomic charge densities, ρ(Ta 2 Natomic), can be acquired by non-self-consistent calculation. The k-point meshes were Γ-centered 17 × 17 × 9, 9 × 9 × 11, 15 × 15 × 15 for the bulk Ta 2 N-I, Ta 2 N-II, bcc Ta; and 1 × 1 × 1 for the molecule (N 2 ), respectively. Herein, total energy calculations were performed in a large supercell of 14 Å × 15 Å × 16 Å for the N 2 molecule. The tetrahedron method with Blöchl corrections [30] was adopted for accurate total energy and electronic structure in integration of reciprocal space. The selected valence electrons were 5p 6 5d 6 6s 2 for Ta, and 2s 2 2p 3 for N. All systems were calculated with the convergence criteria 0.02 eV/Å for forces on ions, and 0.01 meV/atom for electrons, respectively.
In order to obtain the equilibrium structure, the four-parameter Birch-Murnaghan equation of state (EOS) was applied to fit the energy vs. volume (E-V) data points [31,32]: where a, b, c, and d are fitting parameters. This EOS can estimate the equilibrium total energy E 0 , volume V 0 , bulk modulus B 0 , and its pressure derivative B 0 . To compare Ta 2 N-I and Ta 2 N-II in terms of their thermodynamic stability, enthalpy of formation at 0 K, E f (kJ/mol), was calculated as: where E Ta 2 N is the total energy of Ta 2 N, and E Ta , E N are the energy per atom of bcc Ta and one N 2 molecule, respectively. n and N are the atom numbers of Ta and all atoms in bulk Ta 2 N. F is Faraday's constant. The charge density difference (also called deformation charge density) can be calculated as [33]: where ρ(Ta 2 N self-consistent ) is the total charge densities obtained by self-consistent calculation. The superposition of the atomic charge densities, ρ(Ta 2 N atomic ), can be acquired by non-self-consistent calculation.

Elastic Calculations
The elastic constants were calculated by the strain-stress method, as demonstrated by Shang et al. [6,[34][35][36]. For this methodology, the strains (ε = (ε 1 , ε 2 , ε 3 , ε 4 , ε 5 , ε 6 ), where ε 1 , ε 2 , ε 3 are normal strains and ε 4 , ε 5 , ε 6 are shear strains) of a crystal are imposed through specifying the lattice vectors in Cartesian coordinates, as shown in Equation (4). The lattice vectors, R, after deformation can be described as: , and → c are the same. In this work, linearly independent and isotropic strains were used, and the strain was set to ±0.01. By first-principles calculations, the corresponding set of stresses σ = (σ 1 , σ 2 , σ 3 , σ 4 , σ 5 , σ 6 ) of the deformed crystal with set of strains ε = (ε 1 , ε 2 , ε 3 , ε 4 , ε 5 , ε 6 ) could be obtained. The elastic constants (C ij ) are then acquired through the general Hook's law This equation can be solved via the singular value decomposition method to get the least square solutions of the elastic constants. For a trigonal crystal with space groups P3m1 and P31m, there are six independent single crystal elastic constants (C 11 , C 12 , C 13 , C 14 , C 33 , and C 44 ). The mechanical stability can be determined by the Born-Huang [37] criterion as: Then, the Voigt and Reuss bounds of bulk modulus and shear modulus are [38]: Through the Viogt-Reuss-Hill approximation [39], the bulk modulus B, shear modulus G, Young's modulus E, and Poisson's ratio ν can be further obtained by: Crystals 2021, 11, 445 5 of 11 The anisotropy of elastic properties can induce microcracks and reduce the mechanical durability. Thus, to characterize the elastic anisotropies, the universal elastic anisotropic index (A U ) [40] is described as follow: Besides, for trigonal crystals, the three-dimensional (3D) anisotropy of Young's modulus E can be obtained by the reciprocals of E [41]: where S ij are the elastic compliance constants and l i (i = 1, 2, 3) are the direction cosines: here, θ is the angle between a certain direction and crystal orientation [001]. φ is crystal orientation [100]. The greater the deviation of the 3D anisotropy diagram from a sphere (isotropy) is, the higher the degree of anisotropy is.

Structure Stabilities
The differences of structure between Ta 2 N-I and Ta 2 N-II are illustrated in Figure 1 and Table 1. It was revealed that the shift of N atoms between layers changes the symmetry and atomic stacking model of Ta 2 N-I, promoting the formation of Ta 2 N-II. To compare the structure and stability of Ta 2 N-I and Ta 2 N-II, the lattice parameters and enthalpies of formation were calculated, as shown in Table 2. With regard to lattice parameters, a good agreement is shown between our calculated values and the reference data [14,19,24,25,42,43]. Especially, the calculated lattice parameters (a = 3.110 Å, c = 4.896 Å) of Ta 2 N-I were very close to the results (a = 3.11 Å, c = 4.88 Å) [24] calculated through the same pseudo-potential (GGA-PBE), which confirms the reliability of this work. In Figure 2, the total energies of Ta 2 N-I and Ta 2 N-II are plotted as a function of volume through Birch-Murnaghan equation of state fitting. It is shown that Ta 2 N-II has lower energies in the whole variation range of volume. As for the enthalpies of formation, the negative values, −78.61 kJ/mol for Ta 2 N-I and −89.43 kJ/mol for Ta 2 N-II, declare that both structures are stable, agreeing well with experimental [42] and theoretical [23] results.
In Figure 2, the total energies of Ta 2 N-I and Ta 2 N-II are plotted as a function of volume through Birch-Murnaghan equation of state fitting. It is shown that Ta 2 N-II has lower energies in the whole variation range of volume. As for the enthalpies of formation, the negative values, −78.61 kJ/mol for Ta 2 N-I and −89.43 kJ/mol for Ta 2 N-II, declare that both structures are stable, agreeing well with experimental [42] and theoretical [23] results. It is worth noting that a more negative enthalpy of formation of Ta 2 N-II means it is more energetically favorable than Ta 2 N-I. Therefore, the ordering change of N atoms in Ta 2 N-II brings a more stable atomic stacking model (N-Ta-N-Ta-···) for Ta 2 N compound, comparing to the close-stacking (N-Ta-Ta-N-Ta-Ta-···) in Ta 2 N-I. That is because the N atoms between the two Ta atom layers can generate covalent bond (Ta-N) layer by layer in Ta 2 N-II but not in Ta 2 N-I. In addition, Figure 3 demonstrates the variation of bond length between the two structures. It is worth noting that a more negative enthalpy of formation of Ta 2 N-II means it is more energetically favorable than Ta 2 N-I. Therefore, the ordering change of N atoms in Ta 2 N-II brings a more stable atomic stacking model (N-Ta-N-Ta-···) for Ta 2 N compound, comparing to the close-stacking (N-Ta-Ta-N-Ta-Ta-···) in Ta 2 N-I. That is because the N atoms between the two Ta atom layers can generate covalent bond (Ta-N) layer by layer in Ta 2 N-II but not in Ta 2 N-I. In addition, Figure 3 demonstrates the variation of bond length between the two structures.  Only one kind of Ta-N bond with the bond length 2.184 Å exists in Ta 2 N-I. However, two types of Ta-N bond, with the shorter bond length (2.155 Å and 2.174 Å), are found in Ta 2 N-II, owing to the different ordering of N atoms. The covalent bond (Ta-N) with a shorter bond length means a stronger bond strength. This also makes a contribution to the higher stability of Ta 2 N-II.

Elastic Properties
To further investigate the difference in properties between Ta 2 N-I and Ta 2 N-II, elastic constants were calculated to reveal their elastic properties. There are six independent elastic constants in a trigonal system: C11, C12, C13, C14, C33, and C44, as shown in Table 3. According to the Born-Huang [37] criterion, both structures are mechanically stable. The difference between our calculated data of Ta 2 N-I and the other theoretical values [25] was acceptable because of the different pseudo-potential (ultrasoft pseudo-potential, USPP) used for calculations. Comparing the C ij of Ta 2 N-I and Ta 2 N-II, all the absolute values of Cij increased except for C13 and C14. The reduction of the C14 absolute value reveals that the Ta 2 N-II became less anisotropic than Ta 2 N-I.  Only one kind of Ta-N bond with the bond length 2.184 Å exists in Ta 2 N-I. However, two types of Ta-N bond, with the shorter bond length (2.155 Å and 2.174 Å), are found in Ta 2 N-II, owing to the different ordering of N atoms. The covalent bond (Ta-N) with a shorter bond length means a stronger bond strength. This also makes a contribution to the higher stability of Ta 2 N-II.

Elastic Properties
To further investigate the difference in properties between Ta 2 N-I and Ta 2 N-II, elastic constants were calculated to reveal their elastic properties. There are six independent elastic constants in a trigonal system: C 11 , C 12 , C 13 , C 14 , C 33 , and C 44 , as shown in Table 3. According to the Born-Huang [37] criterion, both structures are mechanically stable. The difference between our calculated data of Ta 2 N-I and the other theoretical values [25] was acceptable because of the different pseudo-potential (ultrasoft pseudo-potential, USPP) used for calculations. Comparing the C ij of Ta 2 N-I and Ta 2 N-II, all the absolute values of C ij increased except for C 13 and C 14. The reduction of the C 14 absolute value reveals that the Ta 2 N-II became less anisotropic than Ta 2 N-I. The calculated bulk modulus B, shear modulus G, Youngs' modulus E, Poisson's ratio ν, B/G ratio, the universal elastic anisotropic index A U , and mechanical stability are summarized in Table 4. A crystal with a higher bulk modulus B signifies a lesser compressibility. Shear modulus G describes the resistance of a material to shear stress. Young's modulus is considered as a measure of the stiffness of a solid. Table 4. The calculated bulk modulus B (GPa), shear modules G (GPa), Youngs' modulus E (GPa), Poisson's ratio ν, B/G ratio, universal elastic anisotropic index A U , and mechanical stability of Ta 2 N-I and Ta 2 N-II. It was shown that Ta 2 N-II possesses a higher bulk modulus B, shear modulus G, and Youngs' modulus E. This means Ta 2 N-II is harder to deform than Ta 2 N-I. In addition, the intrinsic ductility and brittleness of Ta 2 N can be estimated by Poisson's ratio ν [44] and B/G ratio [45]. ν> (<)0.26, B/G ratio> (<)1.75 can be regarded as a performance index to distinguish brittle (ductile) materials. Notably, as shown in Table 4, the ν (0.29) and B/G ratio (2.09) decreased to 0.26 and 1.73 in the structural change from Ta 2 N-I to Ta 2 N-II, respectively. Namely, Ta 2 N-I is ductile and Ta 2 N-II is brittle. Hence, a different ordering of N atoms in Ta 2 N-II leads to a ductile-brittle transition in Ta 2 N compounds, which can markedly impact the mechanical properties. Furthermore, the universal elastic anisotropic index (A U ) and the 3D Young's modulus anisotropy were calculated to characterize the elastic anisotropy of Ta 2 N compounds. Ta 2 N-II with A U = 0.55 behaved with a weaker elastic anisotropy than Ta 2 N-I (1.42). Isotropic materials show the shape of a sphere in the 3D diagram of Young's modulus anisotropy. Figure 4 clearly demonstrates the difference of anisotropy between Ta 2 N-I and Ta 2 N-II. Obviously, the shape is more like a sphere in Figure 4b, corresponding to the weaker anisotropy of Ta 2 N-II. A crystal with weaker anisotropy means the possibility of stress concentration is lower and the microcracks are harder to extend. Therefore, the higher resistance to deformation of Ta 2 N-II was further confirmed. intrinsic ductility and brittleness of Ta 2 N can be estimated by Poisson's ratio ν [44] and B/G ratio [45]. ν > (<) 0.26, B/G ratio > (<) 1.75 can be regarded as a performance index to distinguish brittle (ductile) materials. Notably, as shown in Table 4, the ν (0.29) and B/G ratio (2.09) decreased to 0.26 and 1.73 in the structural change from Ta 2 N-I to Ta 2 N-II, respectively. Namely, Ta 2 N-I is ductile and Ta 2 N-II is brittle. Hence, a different ordering of N atoms in Ta2N-II leads to a ductile-brittle transition in Ta 2 N compounds, which can markedly impact the mechanical properties. Furthermore, the universal elastic anisotropic index (A U ) and the 3D Young's modulus anisotropy were calculated to characterize the elastic anisotropy of Ta 2 N compounds. Ta 2 N-II with A U = 0.55 behaved with a weaker elastic anisotropy than Ta 2 N-I (1.42). Isotropic materials show the shape of a sphere in the 3D diagram of Young's modulus anisotropy. Figure 4 clearly demonstrates the difference of anisotropy between Ta 2 N-I and Ta 2 N-II. Obviously, the shape is more like a sphere in Figure 4b, corresponding to the weaker anisotropy of Ta 2 N-II. A crystal with weaker anisotropy means the possibility of stress concentration is lower and the microcracks are harder to extend. Therefore, the higher resistance to deformation of Ta 2 N-II was further confirmed.

Electronic Properties
To further probe the influence of the ordering change of N atoms in Ta 2 N compounds, the total and partial density of states were calculated, and are presented in Figure   Figure 4. The three-dimensional anisotropy diagram of Young's modulus E of (a) Ta 2 N-I and (b) Ta 2 N-II.

Electronic Properties
To further probe the influence of the ordering change of N atoms in Ta 2 N compounds, the total and partial density of states were calculated, and are presented in Figure 5. Comparing Figure 5a,b, the similarities of the two structures are as follows.

Electronic Properties
To further probe the influence of the ordering change of N atoms in Ta 2 N compounds, the total and partial density of states were calculated, and are presented in Figure  5. Comparing Figure 5a,b, the similarities of the two structures are as follows.  proving the formation of the covalent bond, which was confirmed to be favorable to the hardness and shear resistance of the materials [46]. However, the energy region of Ta-N hybridization broadens (from −6~−8eV to −5~−8eV) in Ta 2 N-II. It is implied that the ordering change of N atoms enhances the hybridization of Ta-d and N-p states and the overlaps of valence bonds between N and Ta atoms. The shorter bond lengths in Figure 3 of Ta-N bond in Ta 2 N-II also contribute to the hybridization enhancement. This can explain the elevation of shear modulus in Ta 2 N-II. Furthermore, the max energy of density of states distribution decreases from 7 eV in Ta 2 N-I, to 3 eV in Ta 2 N-II, which is beneficial to the higher stability of Ta 2 N-II.
For analyzing the bonding interaction and electron transfer of Ta and N atoms, charge density difference was calculated, as shown in Figure 6. Herein, Figure 6a,b illustrate the sectional view of charge density difference of the (110) lattice plane in Ta 2 N-I and the (100) lattice plane in Ta 2 N-II, respectively. The lattice plane (110) of Ta 2 N-I and (100) of Ta 2 N-II are equivalent. The red areas denote the gain in electrons, and the blue areas are the loss of electrons.
Moreover, Figure 6c-f exhibit the side and top view of charge density difference of Ta 2 N-I and Ta 2 N-II, in which the blue and yellow regions indicate the gain and the loss of electrons, respectively. Electrons transfer from Ta to N during the formation of the Ta-N bond. The electrons shared between the Ta and N atoms show the interaction of Ta and N, indicating the orbital hybridization of Ta-d and N-p. Moreover, in Figure 6a, the electron distribution is also found between two Ta atom layers in Ta 2 N-I. That means there is an interaction between Ta and Ta atoms in Ta 2 N-I. However, such a phenomenon becomes unapparent in Ta 2 N-II, owing to the ordering change of N atoms, as shown in Figure 6b. It is clear that the electrons no longer accumulate around the Ta atoms like Ta 2 N-I. It can be seen in Figure 6d,f that these electrons transfer and assemble around the new N atom layer formed in the Ta 2 N-II. Thus, the ordering change of N atoms has an outstanding influence on the electron distribution of Ta 2 N. As seen in Figure 6c-f, the electron distribution of Ta 2 N-II is more homogeneous in the unit cell than that of Ta 2 N-I, which may explain why the anisotropy of Ta 2 N-I is stronger.
For analyzing the bonding interaction and electron transfer of Ta and N atoms, charge density difference was calculated, as shown in Figure 6. Herein, Figure 6a,b illustrate the sectional view of charge density difference of the (110) lattice plane in Ta 2 N-I and the (100) lattice plane in Ta 2 N-II, respectively. The lattice plane (110) of Ta 2 N-I and (100) of Ta 2 N-II are equivalent. The red areas denote the gain in electrons, and the blue areas are the loss of electrons.

Conclusions
The stability, elastic, and electronic properties of Ta 2 N-I and Ta 2 N-II were investigated using DFT-based first-principles calculations. The calculated lattice parameters, enthalpies of formation, and elastic constants were a good match with the references, which confirms the reliability of our calculations. According to the enthalpy of formation, Ta 2 N-II is energetically favorable compared to Ta 2 N-I, due to the formation of a more stable atomic stacking model (N-Ta-N-Ta···) and the stronger covalent bond (Ta-N) in Ta 2 N-II. As for the elastic properties, Ta 2 N-II possesses a higher bulk modulus B, shear modulus G, and Youngs' modulus E; but a weaker anisotropy. Especially, a ductile-brittle transition was found between Ta 2 N-I (ductile) and Ta 2 N-II (brittle). It can be assumed from the partial density of states that both Ta 2 N-I and Ta 2 N-II are metallic. In addition, Ta 2 N-II revealed the stronger orbital hybridization of Ta-d and N-p, resulting in stronger covalent bonding, which can explain the elevation of shear modulus. The charge density difference clearly illustrates the interaction of the Ta-N bond and electron distribution difference between Ta 2 N-I and Ta 2 N-II. These analyses of property differences, because of the ordering change of N atoms in Ta 2 N compounds can further promote their application and research.

Data Availability Statement:
The data presented in this study are available within the article.