Electronic, structural, and elastic properties of metal nitrides XN (X = Sc, Y): A first principle study

We utilized a simple, robust, first principle method, based on basis set optimization with the BZW-EF method, to study the electronic and related properties of transition metal mono-nitrides: ScN and YN. We solved the KS system of equations self-consistently within the linear combination of atomic orbitals (LCAO) formalism. It is shown that the band gap and low energy conduction bands, as well as elastic and structural properties, can be calculated with a reasonable accuracy when the LCAO formalism is used to obtain an optimal basis. Our calculated, indirect electronic band gap (E$^\mathrm{\Gamma-X}_g$) is 0.79 (LDA) and 0.88 eV (GGA) for ScN. In the case of YN, we predict an indirect band gap (E$^\mathrm{\Gamma-X}_g$) of 1.09 (LDA) and 1.15 eV (GGA). We also calculated the equilibrium lattice constants, the bulk moduli (B$_{o}$), effective masses, and elastic constants for both systems. Our calculated values are in excellent agreement with experimental ones where the latter are available.

Even with the vast number of experimental and computational studies already available for XN, a satisfactory description of their electronic, transport (effective mass), elastic, and structural properties is still an area of active research. Experimentally, ScN is known to be a semiconductor with a band gap in the range of 0.90 ± 0.1 eV to 1.32 ± 0.3 eV [4,7], while most DFT calculations utilizing local density approximation (LDA) and generalized gradient approximation (GGA) potentials found it to be a metal [3,18,30,[36][37][38][39]. Recent Green's function quasiparticle [32,40], exact exchange [36,41] and screened exchange [18,42] calculations have reported that ScN is a semiconductor with an indirect gap (E Γ−X g ) in the range 0.54 -1.70 eV. For YN, few theoretical calculations utilizing various forms of DFT potentials have reported that it is semi-metallic or semiconductor with an indirect band gap (E Γ−X g ) of 0.2 -0.3 eV [29,43,44], 0.54 eV [39], and 0.80 eV [45]. We are not aware of any experiments reporting the indirect band gap of YN.
Theoretical computations have had difficulties in predicting the correct band gap energy and other related electronic properties of XN from first principle. Indeed, the "band gap problem" is decades old. Several approaches to solve it have been proposed with significant successes. Density functional theory plus additional Couloumb interactions (DFT+U) formalism [46][47][48][49] has had good successes in obtaining correct energy bands and gaps of materials, but can only be applied to correlated and localized electrons, e.g., 3d or 4f in transition and rare-earth oxides. Hybrid functionals [50][51][52] have also been used in attempts to improve on the calculated energy bands and band gaps of materials. This approach involves a range separation of the exchange energy into some fraction of nonlocal Hartree-Fock exchange potential and a fraction of local spin density approximation (LSDA) or generalized gradient approximation (GGA) exchange potential. We should note that this range separation is not universal. There is always a range separation parameter ω which varies between 0 and ∞. While it is reasonably clear that there exists a value of ω that gives the correct gap for a given system, this ω is not universal as its value is always adjusted from one system to another [53,54]. For example, in HSE06 [50,55], ω = 0.11a −1 0 (a 0 is the Bohr radius) and in Perdew-Burke-Ernzerhof (PBEh) global hybrid [56], it is 25 % short-range exact exchange and 75 % short-range PBE exchange. Even though the HSE functional, in most cases, accurately reproduces the optical gap in semiconductors, it severely underestimates the gap in insulators [54,57] and the band width in metallic systems is generally too large [54,[57][58][59], The Engel and Vosko [60] (EV) GGA and the Tran and Blaha [61] modified Becke-Johnson (TB-mBJ) potentials have also provided some improvements to the calculated band gap of materials. For TB-mBJ, while the band gaps are considerably improved, the effective masses are severely underestimated. [59] In the case of arXiv:1206.4277v1 [cond-mat.mtrl-sci] 19 Jun 2012 the EV-GGA potential, the equilibrium lattice constants are far too large as compared to experiment and, as such, leads to an unsatisfactory total energy [62,63].
The theoretical underestimations of band gaps and other energy eigenvalues have been ascribed to the inadequacies of density functional potentials for the description of ground state electronic properties of XN [29]. Also, other methods [64,65] that entirely go beyond the density functional theory (DFT) do not obtain the correct band gap values of most semiconductors without adjustment or fitting parameters [29]. This unsatisfactory situation is a key motivation for our work. In light of the promising technological properties of these materials, parameter-free computations could aid significantly in the design and fabrication of devices. To this end, we have also investigated the elastic properties of ScN and YN within our parameter-free method. In YN, we predict several of these properties for which there are yet no experimental data.
In this paper, we utilized a simple and robust approach based on basis set optimization. The rest of this paper is organized as follows. After this introduction in Section I, the computational method and details used in our work are given in Section II. Section III shows our computed results. Calculated electronic structures are given in Subsection III A. The results for the chemical bonding and structural properties are presented and discussed in Subsections III B and III C. We will conclude with a summary in Section IV.

A. Method
We utilized the electronic structure package developed at Ames Laboratory of the U.S. Department of Energy (DOE), Ames, Iowa [66]. For the LDA computations, we used the Ceperley and Alder DFT exchange-correlation contribution [67] as parameterized by Vosko-Wilk-Nusair [68]. We refer to it as the CA-VWN potential. The GGA calculations were carried out using the Ceperley and Alder DFT exchange correlation contribution [67] as parameterized by Perdew and Wang [69]. We refer to it as the CA-PW potential.
We employed the linear combination of atomic orbitals (LCAO) approach where an unknown wave function, for the solid state, is expressed as a linear combination of atomic orbitals of the individual atoms in the system. The radial parts of these orbitals are generally exponential or Gaussian functions resulting from self consistent calculations of energy levels of the atomic or ionic species that are present in the solid under study. We use Gaussian functions and refer to that rendition of LCAO as the linear combination of Gaussian orbitals (LCGO).
In addition to the use of DFT potentials and of the LCAO formalism, our computational approach rests on the Bagayoko, Zhao, and Williams [70,71] method as recently enhanced by the work of Ekuma [72][73][74][75] and Franklin [76] (BZW-EF). The method heeds the statement by Kohn and Sham [77] that the equations defining LDA have to be solved self-consistently. Upon the selection of an exchange correlation potential, these equations reduce to (1) the one giving the ground state charge density in terms of the wave functions of the occupied states and (2) the Kohn-Sham equation. We describe below the aforementioned method.
In the BZW-EF method, we begin our calculations with the minimum basis set, the one that is just large enough to account for all the electrons in the system under study (in this case XN). Our self consistent calculation with this basis set is followed by another where the basis set is augmented with one additional orbital from one of the atoms or ions in the system. The comparison of the occupied energies from Calculations I and II generally shows that they are different, with those from Calculation II generally lower than their corresponding ones from I. A third, self consistent calculation is performed with a basis set that includes that for Calculation II plus another orbital from one of the atoms in the system. This process of augmenting the basis set and of carrying out self consistent calculations continues until a calculation, say N, is found to have the same occupied energies as calculation (N + 1) that follows it, within computational uncertainty of about 50 meV. This convergence of the occupied energies identifies the basis set of Calculation N as the optimal one. The optimal basis set is thus the smallest one with which the occupied energies verifiably reach their respective minima [73,78]. The resulting eigenvalues and the corresponding eigenfunctions provide the physics description of the system.
In calculations with basis set larger than the optimal one, the ground state charge density does not change, nor do the Hamiltonian and the eigenvalues of the occupied states. Consequently, these calculations do not lower any occupied energies (as compared to the results obtained with the optimal basis set), even though they generally lead to some lower, unoccupied energies [72,73] by virtue of the Rayleigh theorem [84,85]. The Rayleigh theorem can be stated as follows: if one solves an eigenvalue equation with two basis sets I and II, with set II larger than I and where I is entirely contained in II, then the eigenvalues obtained with set II are lower than or equal to their corresponding ones obtained with basis set I (i.e., E N +1 . This theorem explains the reasons that some unoccupied energies are lowered when the Kohn-Sham equation is solved with basis sets larger than the optimal one. Such a lowering of unoccupied energies with basis sets larger than the optimal one is fundamentally different from the one that occurs before the occupied en-ergies reach their minima. The latter lowering is ascribed to physical interactions given that both the charge density and the Hamiltonian change, from one calculation to the next, before one reaches the optimal basis set, while the former one is unphysical. Important points about this method include the fact that it provides the needed variational freedom for the restructuring of the electronic cloud in the material under study, as opposed to the case in isolated atoms or ions. Indeed, the methodical increase of the basis set remedies possible radial, angular symmetry, and size deficiencies by the time the optimal basis set is reached. By changing the basis set, i.e., the external input to the charge density equation, the method solves self-consistently both this equation and the Kohn Sham equation in a way that is quintessentially different from input changes limited to new expansion coefficients in the iterative process [86]. The latter input changes cannot remedy any significant deficiency of the basis set, i.e., in terms of radial orbital, angular symmetry, or the total number of functions (the dimension of the Hamiltonian). In light of the Rayleigh theorem, a deficiency could be over-completeness as well, in single trial basis set calculations. The critical importance of the BZW-EF method resides in part in the fact that, by its very variational derivation, the validity and the related physical content of the eigenvalues of the Kohn Sham equation directly hinge on the final, self-consistent charge density being as close as possible to that of the real material under study.
The thesis work of Williams [71] showed that the method may not be needed for the description of materials established, experimentally, to be metals. At least one band crosses the Fermi level in a metal and the lowestlaying conduction bands reach their minima at the same time as the valence bands do. Consequently, the basis set and variational effect [70,71] that explains the lowering of unoccupied energies, while the Hamiltonian does not change, does not affect these lowest-laying conduction bands in metals, even though it lowers some other conduction bands. This fact partly explains the early success of DFT for the description of metals as opposed to semiconductors and insulators. In the initial BZW method, we increased the basis set by adding orbitals in the order of increasing energies in the atomic or ionic species. Recent works by Ekuma [72][73][74][75] and Franklin [76] led us to the realization that, for the valence states, polarization (i.e., p, d, and f orbitals) have primacy over spherical symmetry (s orbital) in diatomic molecules and more complex systems with three or more atomic or ionic sites. Hence, while still useful, adding orbitals in the order of increasing energies in atomic or ionic species can be relaxed.

B. Computational Details
In the ground state, both ScN and YN have the rocksalt crystal structure. We utilized room temperature experimental lattice constant of 4.501Å [7,[87][88][89] and 4.837Å [88] for ScN and YN, respectively. Each Sc (or Y) atom is surrounded by six N atoms, thereby providing the octahedral environment at the Sc (or Y) site, which leads to the splitting of the degenerate d orbital into t 2g and e g states. Every Sc (or Y) atom has 12 nearest neighbors Sc (or Y) atoms and 6 next nearest neighbors Sc (or Y) atoms. Preliminary calculations indicated that scandium and yttrium are closer to Sc 3+ and Y 3+ than to the neutral Sc and Y, respectively. Similarly, nitrogen is N 3− as opposed to the neutral N. Therefore, we first carried out self consistent calculations of the electronic properties of X 3+ (X = Sc, Y) and N 3− . Atomic orbitals utilized in these calculations, for the valence states, are given between parentheses: Sc 3+ (3s3p3d4s4p4d) and N 3− (2s2p3s3p) for ScN and Y 3+ (4s4p4d5s5p) and N 3− (2s2p3s3p) for YN. Other atomic states with higher binding energies were treated as deep core states. In the optimal basis set for the valence states of ScN, (4p4d) and (3p) are unoccupied for Sc 3+ and N 3− , respectively. For YN, the unoccupied orbitals in the optimal basis set for the valence states are (5p) and (3p) for Y 3+ and N 3− , respectively. Nevertheless, these orbitals are included in the self-consistent LCAO calculations to allow for a reorganization of electronic cloud in the solid environment, including polarization.
In the self-consistency calculation, both the potential and charge density are also expanded in terms of even tempered Gaussian orbitals (Sc: 15, 15, and 13; Y: 15, 15, and 13; N: 17, 17, and 15 for s, p, and d orbitals, respectively). The exponents, α, of the Gaussian basis sets range from 0.19 to 10 5 . The charge fitting error using the Gaussian functions in the atomic calculation was about 10 −5 . Since the deep core states are fully occupied and are inactive chemically in the materials, the charge densities of the deep core states are kept the same as in the free atoms. However, the core states of low binding energy are still allowed to fully relax, along with the valence states, in the self consistent calculations. The computational error for the valence charge was 1.2 × 10 −5 and 5.9 × 10 −6 per valence electron for ScN and YN, respectively. The self consistent potential converged to a difference of 10 −5 after about 60 iterations.
The Brillouin zone (BZ) integration for the charge density in the self consistent iterations was based on 28 special k points in the irreducible BZ (IBZ). The energy eigenvalues and eigenfunctions were then obtained at 161 special k points in the IBZ for the band structure. A total of 152 weighted k points, chosen along the high symmetry lines in the IBZ of ScN and YN, respectively, were used to solve for the energy eigenvalues from which the electron densities of states (DOS) were calculated using the linear, analytical tetrahedron method [90]. The partial density of states (pDOS) and the effective charge of each atom were calculated using the Mulliken charge analysis procedure [91]. We also calculated the equilibrium lattice constant (a o ), the bulk modulus (B o ) and the electron effective masses for different directions with respect to the Γ point. For the calculation of the equilibrium lattice constant, we utilized the Murnaghan equation of state [92,93]. By applying an appropriate set of strains to the undeformed unit cell lattice, we calculated the elastic constants from the resulting change in total energy on the deformation using the strain-energy method. For a typical cubic crystal, three independent elastic moduli are of importance; they are usually denoted as C 11 , C 12 , and C 44 .

III. RESULTS AND DISCUSSION
The results of the electronic structure computations are given in Figs. 1 to 2. Figure 3 shows the contour plot of the distribution of the electron charge density while Figure 4 depicts the total energy curves of XN. We discuss the electronic structure and the effective masses in III A. Densities of States, chemical bonding, and electron distribution are described in Subsection III B. The structural and elastic properties are presented in Subsections III C.

A. The Electronic Structure, Band gap and Effective mass
Figures 1 to 2 exhibit the energy bands, and the related total (DOS) and partial (pDOS) densities of states of XN. All energies are referred relative to zero energy at the top of the valence band (VB). The electronic structures of the valence bands, the low energy conduction bands, and the band gap determine the most important properties of these materials in device applications.
Our ab-initio method shows that the fundamental gaps of both ScN and YN are indirect ones, with the maximum of the valence band (V B max ) occurring at Γ and the minimum of the conduction band (CB min ) at the X point (see Figs. 1(a) and 2(a) for ScN and YN, respectively). Our calculated, indirect gaps for ScN, using the equilibrium lattice constant, is 0.79 eV and 0.88 eV for CA-VNW (LDA) and CA-PW (GGA) potentials, respectively. The difference in the band gaps and energy eigenvalues between CA-VWN and CA-PW may be attributed to the enhancement factor (s = | n| 2kFn ) [94,95] in the GGA functional. The band gap of ScN has been reported to suffer from the Burnstein-Moss shift [96,97] due to large background carrier concentration. As a consequence of this effect, several experimental indirect band gap values ranging from 0.90 ± 0.1 eV [4] to 1.32 ± 0.3 eV [98] have been reported. Our computed band gap values are in good agreement with the lower, experimental values that result from relatively low carrier concentrations. For YN, our computations predict an indirect band gap value of 1.09 eV and 1.15 eV for LDA and GGA, respectively. To our knowledge, no measurement of the YN indirect gap has been reported. We expect YN to suffer from the same Burnstein-Moss effect as ScN, due to the close similarity of their structures.
Our calculated valence bands for ScN (see. Fig. 1(a)) and YN (see Fig. 2(a)) resemble those from other calculations available in literature. A major difference is the size of the band gap. For both ScN and YN, there is significant N 2p -X id (X = Sc, Y; and i = 3, 4 for Sc and Y, respectively) hybridization in the valence bands. The low energy conduction bands up to about 5.24 eV are mainly from X id (see Figs. 1(b) and 2(b) for ScN and YN, respectively). The following can further be confirmed from the partial density of states (see Figs. 1(c) and 2(c) for ScN and YN, respectively). The upper valence bands originated from the bonding between N 2p states and an admixture of 3d states (for Sc) and 4d states (for Y). The lowest conduction bands are mainly the antibonding t 2g states with 3d (for Sc) and 4d (for Y) character, respectively. The two nonbonding e g bands are higher in energy. These bonding characters are in perfect agreement with the bonding analysis of Harrison and Straub [99]. The authors showed that rock-salt XN have 3 p-like bonding, 3 d-like antibonding t 2g , and two d-like nonbonding e g bands formed by the hybridization of three valence p states of N with the five d states of X. The assignment of the bondings is consistent with the partial density of states analysis of the screened exchange LDA FLAPW calculations of Stampfl and co-workers [18,31] and of the exact exchange based quasiparticle calculations of Qteish et al. [36], respectively.
The density of states (DOS) of ScN is shown in Fig.  1(b). In the valence bands, we found sharp peaks at -2.46 ± 0.1 eV and -2.96 ± 0.1 eV, respectively. A shallow minimum is observed at -3.05 ± 0.1 eV, followed by a shoulder at -3.2 ± 0.1 eV. A broad peak is located at -4.05 ± 0.2 eV. These peaks are formed by a strong hybridization between N p states and Sc d states with little contribution from Sc p and s states, respectively. In the conduction bands, the first observed, significant feature is a shoulder at about 3.80 ± 0.1 eV. It is due to the hybridization between N p and Sc d states. A sharp peak is observed at 4.7 ± 0.1 eV; it is formed from all the states of the atomic species except N s states. The density of states of YN (see Fig. 2(b)) is similar to that of ScN, except that it has a pronounced shoulder at 0.92 eV below the Fermi energy. In the DOS of the valence bands, we predicted a sharp peak at -1.72 ± 0.1 eV, followed by two smaller ones at -2.1 ± 0.2 eV and -2.82 ± 0.1 eV, respectively. These peaks are mainly of N p states and Y d states. In the conduction bands, we observed a small peak at 3.6 ± 0.1 eV, which is mainly of N p states and Y d states. A sharp peak is observed at 5.94 ± 0.1 eV. This peak is of N s, N p, and Y d states character.
To determine some specific properties of XN relevant to transport, we calculated the effective mass in different directions with respect to the Γ point. The effective masses are calculated from curves that fit the bands in the immediate vicinity of minima (for electrons) or a maxima (for holes). The calculated electron effective masses (in units  Fig. 1(a). The Fermi energy (EF) has been set equal to zero. (c) The calculated partial density of states (pDOS) of c-ScN, as obtained from the bands shown in Fig. 1(a). The Fermi energy (EF) has been set equal to zero.  Fig. 2(a) The Fermi energy (EF) has been set equal to zero. (c) The calculated partial density of states (pDOS) of c-YN, as obtained from the bands (one with solid line) shown in Fig. 2(a). The Fermi energy (EF) has been set equal to zero.  Figure 3 shows the contour plots of the distribution of the electron charge densities of XN along the (100) plane, cutting through the atoms. Away from the atomic centers (i.e., in the interstitial region) the electron charge density distributions are not spherically symmetric. The chemical bonding in ScN and YN appears to be of inter- mediate character between ionic bonding in the nonconducting calcium nitride or strontium nitride and that of the metallic bonding prevailing in the conducting group IV transition nitrides. As a consequence of this complex bonding, there seems to be a cooperative/competing bonding mechanisms in XN. As can be seen from Fig. 3, the bonding between X and N is covalent, with a significant ionic character. The bond length of Sc -N is 2.25Å while that of Y -N is 2.42Å. The experimental Sc -N bond length is 2.24Å [7,88] while the experimental bond length, Y -N, is 2.44Å [100]. The covalent bonding in ScN can be seen to be stronger than that in YN. This difference is understandable in light of the larger Y -N bond.

C. Structural and Elastic Properties
Total energy versus lattice constant data are as shown in Fig. 4 for ScN and YN. The computation was done with both the CA-VWN and CA-PW potentials. However, the reported data are for CA-PW, as we noted that the total energy difference between results from the two potentials were small (± 0.12 eV  Using the elastic constant relations as in Ref. [101], we calculated the independent elastic constants C 11 , C 12 , and C 44 . The calculated independent elastic constants is used to determine the compliance constants S ij (S 11 , S 12 , and S 44 ) [101,102]. The Poisson ratios ν are measure of materials tendency to react to applied strain. For a cubic system (averaged over transverse directions in the [100], [110], and [111] directions), they are given by: [103,104] ν 100 = -S12 S11 , ν 110 = -S11+2S12−S44/2) S11+2S12+S44 , and ν 111 = -S11+3S12−S44/2 2S11+2S12+S44 . The isotropic shear modulus G [105] can be estimated from the average of the shear modulus in the Voigt approximation (G V = 1/5(C 11 -C 12 + 3C 44 )) [106] and Reuss approximation (G R = 5C44(C11−C12) 4C44+3(C11−C12) ). [107] It is given as G = (G V + G R )/2. Note that the Voigt and Reus approximations represent the upper and lower bounds of the isotropic shear modulus. In a typical cubic system, the isotropic bulk modulus K, can be calculated using the relation K = 1 3 (C 11 + 2C 12 ). Using these quantites, we obtain the isotropic Young's modulus, E = 9KG 3K+G and the isotropic Poisson ratio,ν = 3K−2G 2(3K+G) ≈ 1 3 (ν 100 + ν 110 +ν 111 ). The equilibrium elastic constants, compliance constants (S 11 , S 12 , and S 44 ), the elastic moduli corresponding to the determined stiffness constants and other elastic parameters are shown in Table I. The calculated, independent elastic constants for ScN are C 11 = 452.55 GPa, C 12 = 98.82 GPa, and C 44 = 184.60 GPa. Using the Zener's anistropy index (ratio) [108,109], η = 2C 44 /(C 11 -C 12 ), we quantify the  [110], and Gall et al. [89] on epitaxial films of ScN. These authors reported values of E in the range 270 ± 25 to 388 ± 20 GPa and v value of 0.15, 0.188 ± 0.002 and 0.20 ± 0.04. Our calculated and measuredν values of ScN are in the same range. Our calculated value of E is close to the upper bound, given above, of the measured values for the films. The slight differences between some of our computed elastic constants and experimental ones may be due to the fact that these experiments were done on films of ScN. We should note that growth conditions and film thickness [7] are known to affect significantly the elastic constants, as is evident in the works of Gall et al. [7], Gall et al. [89], and Moram et al. [110] which showed a wide range of values. For YN,Ē = 357.94 GPa, G = 149.53 GPa, andν = 0.20. We note that our computed Poisson ratios,ν are comparable to known Poisson ratios for other transition-metal nitrides that are in the range 0.19 -0.22 [101,111,112]. We are not aware of any experimental measurements of E, G, and v for YN; hence, our calculated results are predictions. In both ScN and YN, our computed isotropic bulk modulus, K is in basic agreement with the bulk modulus we obtained for the equilibrium structural optimization.
Our calculated elastic stiffness constants, obey the stability conditions for a cubic crystal system, i.e., B o = (C 11 + 2C 12 )/3 > 0, G = (C 11 -C 12 )/2 > 0, and C 44 > 0 [113]. The calculated η values are both close to 1, implying that the elastic behavior is almost isotropic. The Cauchy pressure (C 12 -C 44 ), in both ScN and YN, is significantly less than 0. It corresponds to a directional bonding [74,114,115], leading to large charge transfers from cations to anions as has been observed experimentally in TiN [116].

IV. SUMMARY AND CONCLUSION
We have performed first principle calculations of electronic and related properties of ScN and YN, based on basis set optimization with the BZW-EF method. Accurate, band gaps, structural parameters, lattice constants, bulk moduli, and elastic constants were calculated. We have provided detailed analyses of the electronic and related properties of both ScN and YN. Our computed properties are in good agreement with experiment. In particular, we predict various electronic (indirect band gap value), effective masses, elastic and structural properties of YN. It is our hope that the present results will motivate further experimental and theoretical studies of the properties of these materials, with emphasis on measurements of electronic and elastic properties of YN.