Trends in the elastic response of binary early transition metal nitrides

Motivated by an increasing demand for coherent data that can be used for selecting materials with properties tailored for specific application requirements, we studied elastic response of nine binary early transition metal nitrides (ScN, TiN, VN, YN, ZrN, NbN, LaN, HfN, and TaN) and AlN. In particular, single crystal elastic constants, Young's modulus in different crystallographic directions, polycrystalline values of shear and Young's moduli, and the elastic anisotropy factor were calculated. Additionally, we provide estimates of the third order elastic constants for the ten binary nitrides.


I. INTRODUCTION
Nitride compounds are a prominent class of materials with applications spanning from protective hard coatings (mostly transition metal nitrides, TMNs, of the IIIB-VIB group but also e.g. BN or SiN) 1,2 , to optoelectronic devices (mostly IIIA and VA groups, but also e.g. ScN or TiN) 3 , to potential hydrogen storage materials such as Li 3 N 4 .
When it comes to their superior mechanical properties such as high hardness and Young's modulus, TMNs, and in particular TiN are often the industrial material of choice for surface coating of e.g. cutting tools. The properties of simple binary compounds can be successfully enhanced by forming metastable alloys, e.g. Ti 1−x Al x N 1 , which at higher temperatures age-hardens before decomposing into its stable constituents 5 . Recent studies have shown beneficial effects (such as increased oxidation resistance or retardation of the final decomposition step to higher temperatures) of additional alloying elements in Ti 1−x Al x N 5-8 and other systems 9 . Another approach how to improve material properties is via multilayer design, where individual layers are typically simple binary or ternary systems [10][11][12] .
A modern way designing new and improve current materials is to combine experiment with modelling. For simple and/or small systems, quantum mechanical first principle approaches can be used. However, these are typically limited to several hundreds of atoms and, e.g. multilayers, crack propagation or nanoindentor tip-layer contact become difficult topics to handle. Here, the continuum mechanics employing finite element method (FEM) proves to be a successful tool [13][14][15] . A key prerequisite to perform FEM calculations is the knowledge of the elastic properties of the studied materials, which are not always experimentally available (e.g., because some phases are stable only in the multilayer arrangement, but not as a bulk material). In such cases the elastic constants can be provided by the first principle calculations.
The literature on first principle calculations of early TMNs (group IIIB-VB) is vast. The main focus of those papers is the electronic structure and related material chemistry problems, while the calculation of the elastic properties is often only a minor part of the results. Additionally, a lot of those reports focus only on one (or a few) systems. There are some exhaustive reports on the chemical trends in the early TMNs [16][17][18][19][20][21] , but apart from Refs. 20 and 21 they do not discuss the elastic properties. In addition, there are some discrepancies between reported values (e.g., C 11 for ZrN between 304 GPa 22 and 616 GPa 20 ) which are worth cross-checking.
The aim of the present paper is to give a comprehensive overview on elastic properties of early transition metal nitrides and AlN, since these materials are or have the potential to be used as protective coatings 23 . In particular, we investigate ten binary systems, AlN, ScN, TiN, VN, YN, ZrN, NbN, LaN, HfN, and TaN. We focus on the cubic variant (B1, F m3m, NaCl prototype), which is the stable configuration of all of them apart from AlN, NbN, and TaN being metastable in this configuration. The single crystal elastic constants are validated by several independent approaches as well as by a comparison with available theoretical and experimental data. Subsequently, we calculate directionally resolved Young's modulus, anisotropy factors and polycrystalline elastic properties of these compounds, and rationalise the trends in terms of their electronic structure and bonding.

A. Deformation modes
The linear elastic response of cubic materials is fully described by three independent components c xxxx , c xxyy , and c xyxy of the fourth rank tensor of the second order elastic constants (SOECs). It is convenient to represent this tensor with a 6 × 6 matrix where C 11 = c xxxx , C 12 = c xxyy , and C 44 = C 66 = c xyxy .
Here we make use of the Voigt notation xx ∼ 1, yy ∼ 2, zz ∼ 3, yz ∼ 4, xz ∼ 5, and xy ∼ 6. An additional relationship links SOECs with the bulk modulus, B, (2) B describes the elastic response of materials to volume changes, and it is obtained as a fitting parameter from the Birch-Murnaghan equation of state 24 . Consequently, two other deformation modes are needed to obtain all independent components of the cubic elastic tensor. The first pair consists of orthorhombic and monoclinic deformations. The orthorhombic mode results in a strain tensor and the corresponding strain energy density, U (δ), is Here, E tot (δ) and E eq are the total energies per unit volume, corresponding to ε orth (δ) and ε orth (0), respectively. A monoclinic deformation yielding a strain tensor is used to evaluate the C 44 elastic constant from the corresponding strain energy density One should note that these two modes keep the unit cell volume constant. The second pair of deformations, which is also often used, is a pair of tetragonal and triclinic distortion. The tetragonal deformation corresponds to a strain matrix producing a strain energy density The C 44 elastic constants is obtained from a trigonal distortion with a strain tensor and a strain energy density These two deformations are volume non-conserving.
Recently, Zhao et al. 25 proposed a set of six deformation matrices which allow for estimation of second and third order 26 elastic constants of cubic materials at the same time. These are: The strain energy density in these cases is where the coefficients A and B are specific combinations of C ij and C ijk as given in Table I.

B. Calculation details
Quantum mechanical calculations employing density functional theory (DFT) 27,28 were carried out using Vienna Ab initio Simulation Package 29,30 .
Projector augmented-wave pseudopotentials 31 together with the generalised gradient approximation (GGA), as parametrised by Wang and Perdew 32 , for the exchange and correlation potential are used. The plane-wave cutoff energies and the k-vector samplings of the Brillouin zone were carefully checked to provide a total energy accuracy in the order of 1 meV/at or better. They are listed in Table II together with the used pseudopotentials; the suffices sv and pv refer to the exact valence configuration taking into account explicitely also the s and p closed-shell electrons, respectivelly.

A. Equilibrium properties
The optimised lattice constants, a, formation energies, E f , and mass densities, ρ, are summarised in Table III. Since there have been a vast number of publications on experimental as well as calculated equilibrium structure parameters of these early TMN compounds (see e.g., Refs. 16,18,19,21,33, and 34, and references therein), we limit the comparison of the here calculated lattice parameters to the experimental data from Ref. 35. The calculated lattice constants are, as expected for GGA, slightly larger than the experimental values. The error is smaller then 1% except for the case of TaN, where a deviation of about 1.5% is obtained. This is likely to be related to the fact that cubic TaN is metastable and prefers N deficient configurations resulting in a significant decrease of the lattice parameter with respect to a stoichiometric configuration 36 .
The trends in the energy of formation, i.e. less negative values as one moves from the IIIB to the VB group, as well as the absolute numbers agree well with those presented by Rovere et al. 34 .
Lastly, from the calculated lattice parameters (equilibrium volume) and the atomic weights we computed mass densities. Again, apart from the TaN case, where the under-stoichiometry of the experimental compound is likely to play a role, the agreement with experimental values is satisfactory.

B. Single crystal elastic constants
When calculating the single crystal elastic constants as described in Sec. II A, one should check how the fitted C ij s depend on the maximum deformation, δ max , i.e., on the range of deformations applied to the unit cell. When δ max is too small, the accuracy of C ij s is likely to be influenced by the numerical inaccuracies of the DFT calculations, while non-linear elastic (and perhaps also plastic) effects are no longer negligible for too large δ max 25 . To illustrate this behaviour, we plot in Fig. 1 the C 11 and C 44 elastic constants of ZrN as a function of δ max and the order of the fitting polynomial. It follows, that with increasing order of the fitting polynomial, the plateau region where the specific elastic constant is independent of δ max , enlarges. At the same time, the onset of the plateau shifts to higher values δ max . The reason is that a high fitting order leads to an over-fitting of the too few data points for a small δ max . In extreme cases, such overfitting may lead to incorrect plateaus, as shown e.g. for C 44 (monoclinic deformation) using a 13 th order fitting polynomial.
In general, the combination of tetragonal and trigonal deformations gives more robust results for the early TMN in the cubic structure then the combination of orthorhombic and monoclinic deformations. Nevertheless, the most robust results in terms of plateau values scatter and the dependence on the order of the fitting polynomial, were obtained when employing the deformation matrices A1 and A6 (see Fig. 1).
Recently, Udyansky et al. 39 showed that the elastic constants of α-Fe are also highly sensitive to the value of the smearing parameter, σ, in the Methfessel-Paxton scheme. We have therefore checked the convergence of C ij s also with respect to σ. It turns out that the cubic early TMNs are not hugely sensitive to σ, but in some cases, e.g., ZrN or NbN, the elastic constant values change by up to 5% when sigma is increased from 0.02 eV to 0.8 eV. Nevertheless, these variations typically take place only for small values of σ, and a converged behaviour is obtained for σ ≈ 0.6 eV.
The single crystal elastic constants are summarised in Table IV. When possible to evaluate, we give the C ij s based on all three methods described here (i.e., orthorhombic+monoclinic, tetragonal+trigonal, and A1+A6 deformation modes). Since the A1+A6 deformation modes were the only ones to provide well converged results for all ten binary nitrides, we show them in Fig. 2 and we will use them in the following analysis for consistency. A comparison with previous DFT-GGA literature data 21,22,40,41,43,44,49 yields, apart from a few exceptions, a good agreement with our results. The local density approximation (LDA) based C 11 elastic constants from the literature 20,42,47 are higher than our GGA-based data. This is a consequence of over-and under-binding of LDA and GGA, respectively, resulting in too small lattice constants and consequently too hard elastic constants in LDA. Finally, although many of the here calculated elastic constants agree well with the available experimental data, in a few cases the discrepancy is as large as 20% (e.g. C 11 of VN).
In all cases, the obtained elastic constants fulfil the stability criteria for cubic crystals The elastic constant C 11 is significantly stiffer than the other two elastic constant. Within each row the C 11 and C 12 elastic constants monotonically increase with increasing atomic number at the same time, C 11 decreases from Sc to Y to La (isovalent IIIB group) while it increases from V to Nb to Ta (isovalent VB group). It has been suggested in the literature 50,51 , that a negative   IV. Single crystal elastic constants. "o", "t", and "A" stand from values calculated by orthorhombic and monoclinic, tetragonal and trigonal, and using the A1 and A6 deformation modes, respectively. "c" and "e" stand for calculated and experimental data from the literature, respectively. The results in bold are used for the further analysis of polycrystalline elastic properties.
Cauchy pressure C 12 − C 44 < 0 corresponds to more directional bonding while positive values indicate predominant metallic bonding. Indeed, the calculated Cauchy pressure is most negative for AlN in which a significantly larger charge transfer from cation to anion takes place as e.g., in TiN 52 . The Cauchy pressure increases to positive values with increasing number of valence electrons within each periodic table row, as those contribute mainly to the metal-metal d-d interactions 19,53 . These trends may be used in the materials selection process to realize specific requirements.

C. Directional Young's modulus
The Young's modulus, E, in a certain direction, ξ, is defined as the ratio of longitudinal stress to longitudinal strain in this direction. The elastic compliances, S ij , are in the case of cubic crystals the solution of the following set of equations 54 : For a cubic crystals E ξ then reads 54 where l 1 , l 2 , and l 3 are the directional cosines of ξ. For 100 , 110 , and 111 directions this becomes The results are plotted in Fig. 3. The Young's modulus in 100 follows mostly the same trend as C 11 , since the C 11 elastic constant has the strongest contribution to E 100 . There is a considerable difference between the semiconducting compounds AlN, ScN, YN, and LaN, in which the 100 direction becomes the softest, and the metallic TiN, VN, ZrN, NbN, HfN, and TaN, where the 100 direction is clearly the strongest. In addition, the Young's modulus of AlN in the 111 direction is more than 1.5-times larger than in any other of the here investigated TMNs. This is mainly caused by the high value of C 44 of AlN, suggesting that AlN is much stronger in shear deformation than the other TMN.
To quantify the anisotropy, we employ the Zener's anisotropy ratio 55 , A, defined as The results, together with the ratio E 111 /E 100 which provides similar information, are shown in Fig. 4. The results suggest that AlN is clearly stiffer in the 111 than in the 100 direction. The opposite result is obtained for the group IVB and VB TMN where the 100 direction is the stiffest. The group IIIB semiconducting TMN exhibit values of both, A and E 111 /E 100 , very close to 1. This implies that their elastic behaviour is almost isotropic. The most isotropic response is predicted for YN with A = 1.05. The (an)isotropy of the Young's modulus is visualised in Fig. 5. Fig. 5b demonstrates the isotropic Some insight into these trends can be gained from considering the differences in bonding. The bonds in AlN are strongly ionic 52 while the TMNs contain a significant part of the covalent bonding 19 . Since the covalent bond is stronger than the ionic, this can rationalise why AlN has the smallest value of E 100 . When going from group IIIB to VB elements within each row, the extra electrons fill the bonding metal-metal orbitals, while the hybridised sp 3 d 2 states move to lower energies 19,53 . This can be interpreted as strengthening of the hybridised sp 3 d 2 bonds which are oriented along the 100 directions. As for the high C 44 elastic constant of AlN, one may argue that since there are no d electrons available to form the metalmetal bonds in the 110 directions (as it is the case for the IVB and VB group elements), upon a shear deformation a significantly increased repulsion between Al-Al and N-N ions as they get closer occurs, which causes the high value of E 111 .

D. Polycrystalline properties
Several models exist which assess the isotropic polycrystalline elastic properties using the anisotropic single crystal elastic constants of a given material. Voigt's approach 56 of constant strains in all grains yields the upper limit, G V and E V , to the polycrystalline shear and Young's moduli, respectively. On the other hand, Reuss 57 proposed to apply constant stresses in all grains, which yields lower limits G R and E R . Taking B V = B R = B, where the bulk modulus, B, is obtained from the Birch-Murnaghan equation of state 24 , one gets Finally, Hershey 58 derived an equation for selfconsistently calculating the shear modulus, G H . In this approach, G H is the real positive root of the following fourth order polynomial This equation can be simplified by dividing it with (3C 11 + 6C 12 + 8G) to a third order polynomial 59 with the same positive real root Subsequently, Eq. 20 is used to estimate the Young's modulus within Hershey's approach.
The thus calculated polycrystalline elastic constants are shown in Fig. 6 together with the bulk modulus. They fit well with the few accessible experimental data-points included (black full symbols). The trends in E and G are akin: the maximum value in each row of the periodic table is obtained for the group IVB TMNs. The spread between G R and G V (shaded in Fig. 6), as well as between E R and E V is very small for the group IIIB and IVB TMNs, suggesting that the elastic properties of polycrystals of these materials will not be hugely influenced by the misorientations of individual grains. A different situation is obtained for AlN and group VB TMNs (in particular, for NbN and TaN), where the Reuss-Voigt range is quite large. As shown in Ref. 64, the ratio between the Voigt and Reuss bounds depends non-linearly on the anisotropy factor A. The ratio becomes particularly large when the anisotropy A approaches 0, as in the case of NbN and TaN. As a consequence, these materials are expected to be strongly affected by the actual microstructures (i.e., not only by the grain orientations, but also by the shape of the grain).
Based on an evaluation of a large experimental data set, Pugh 65 proposed that the higher (lower) the B/G ratio is, the more ductile (brittle) the material is. This ratio is plotted in Fig. 4. In general, the ductility increases from IIIB to VB group (e.g., with increasing number of valence electrons and thus increasing amount of metallic bonding), and within each group from lighter to heavier elements.

E. Third order elastic constants
The methodology employing the deformation matrices A1-A6 allows also to easily estimate third order elastic constants (TOECs), by following Eq. 12 and relations in Table I. TOECs, C ijk , appear in the Taylor series expansion of the strain energy where ε αβ = (1 + δ αβ )η i /2 is the relationship between components ε αβ of the Lagrangian strain tensor and six components η i of a corresponding vector in Voigt notation 59 . According to the above equation, TOECs give corrections when applying such large strains that linear elasticity no longer applies. TOECs are thus useful to describe the pressure dependence of second order elastic constants, C ij , or thermal properties of solids 66 . This can be of a particular interest for thin films where residual stresses in the range of several GPa can be realized.
The computed TOECs for the binary systems investigated in this work are summarized in Table V. For isotropic aggregates of cubic crystals, Lubarda 59 derived equations for Voigt-and Reuss-type averages of TOECs. These equations are equivalents to the elastic constants expressed by Eqs. 18-20. The corresponding formulae are briefly summarised in Appendix A. The three polycrystalline Voigt-and Reuss-type TOECs, C 123 , C 144 , and C 456 are presented in Fig. 7 as the upper and lower boundaries of the shaded areas. These boundaries provide an estimate for the expected spread of the data and depends on the actual microstructure. The lines in Fig. 7 represent the Hill's average of the Voigt-type, C ijk,V , and Reuss-type, C ijk,R TOECs. The six TOECs describing a crystal with the cubic symmetry are in the isotropic case related by The results suggest that in each row of the periodic table, C 123 and C 456 decrease to more negative values with increasing number of valence electrons from 3 to 5. For C 144 no clear trend is observed. Since the TOECs are mostly negative, second order elastic constants get stiffer with compressive stresses while they soften under tension (compare with Eq. 23).

IV. CONCLUSIONS
Calculating elastic properties using density functional theory is a powerful technique, in particular when material phases single crystals are not experimentally accessible. In this paper we provided a coherent description of the elastic behaviour of nine binary early transition metal nitrides (ScN, TiN, VN, YN, ZrN, NbN, LaN application-tailored properties. Single crystal elastic constants, C ij , and directionally resolved Young's moduli, E, in 100 , 110 , and 111 directions are provided. The results clearly indicate the special position of AlN. This material has the largest Young's modulus E along 111 , while all group IVB and VB nitrides exhibit the largest E value along 100 . These trends could be rationalised by analysing the bonding characteristics of these compounds. Computing the elastic anisotropy we find that YN followed by ScN and LaN are the materials closest to the elastically isotropic behaviour. Finally, the polycrystalline elastic properties (Young's and shear modulus) were calculated. Good agreement with the rather scarce available experimental data was obtained.
where A is the anisotropy ratio given by Eq. 17. It can be seen that in case of an isotropic materials (A = 1), the two approaches give the same results (compare e.g., Eqs. A3 and A4.