Macroscopic Elastic Properties of Textured ZrN--AlN Polycrystalline Aggregates: From Ab initio Calculations to Grain-Scale Interactions

Despite the fast development of computational materials modelling, theoretical description of macroscopic elastic properties of textured polycrystalline aggregates starting from basic principles remains a challenging task. In this communication we use a supercell-based approach to obtain the elastic properties of random solid solution cubic ZrAlN system as a function of the metallic sublattice composition and texture descriptors. The employed special quasi-random structures are optimised not only with respect to short range order parameters, but also to make the three cubic directions $[1\,0\,0]$, $[0\,1\,0]$, and $[0\,0\,1]$ as similar as possible. In this way, only a small spread of elastic constants tensor components is achieved and an optimum trade-off between modelling of chemical disorder and computational limits regarding the supercell size is achieved. The single crystal elastic constants are shown to vary smoothly with composition, yielding $x\approx0.4$-0.5 an alloy constitution with an almost isotropic response. Consequently, polycrystals with this composition are suggested to have Young's modulus independent on the actual microstructure. This is indeed confirmed by explicit calculations of polycrystal elastic properties, both within the isotropic aggregate limit, as well as with fibre textures with various orientations and sharpness. It turns out, that for low AlN mole fractions, the spread of the possible Young's moduli data caused by the texture variation can be larger than 100 GPa. Consequently, our discussion of Young's modulus data of cubic ZrAlN contains also the evaluation of the texture typical for thin films.


I. INTRODUCTION
Quantum mechanical calculations using Density Functional Theory (DFT) of structural as well as elastic properties of materials have become a standard tool in modern computational material science. Recently, also the alloying trends have been heavily investigated, which in the area of hard protective coatings addressed predominantly issues related to the phase stability (see e.g. Refs. 1-6). This has been possible due to the increased computational power, and the development of theories for treating random solid solutions. These include effective potential methods 7 (e.g., coherent potential approximation or virtual coherent approximation), cluster methods 8 (e.g., cluster expansion method), or supercell-based approaches, such a quasi-random structures (SQSs) 9 technique employed in this paper.
While the bulk modulus is relatively easy to obtain from the Birch-Murnaghan equation of state 10 as used during the structure optimisation, the full tensor of elastic constants, C ij , requires additional calculations. The two common methods to calculate C ij from the first principles are the total energy method and the stress-strain method. The latter relies on availability of stress tensor and uses Hooke's law to evaluate C ij directly. On the other hand, the total energy method assigns an energy difference between a ground and a deformed state to the strain energy. This is a function of applied strain and a specific combination of the elastic constants. The advantage of this method is that the total energy is always available from ab initio calculations, and it furthermore allows for estimation of higher order elastic constants 11 . The disadvantage is that it usually takes more CPU resources than the stress-strain method as more deformation modes need to be applied. It has been also recently proposed that the stress-strain method is a more robust technique 12 .
When it comes to the elastic constants of materials without any long-ranged periodicity, the supercell approach faces an apparent problem: on one hand, the distribution of atoms on the lattice sites is required to be as random as possible to mimic solid solutions, hence often leading to supercells with only primitive symmetry (space group P 1). On the other hand, the material is expected to exhibit certain symmetry based on its underlying lattice, for example the cubic symmetry of nitride coatings with the B1 (NaCl) structure. A combined ab initio and molecular dynamics study 13 has shown, that when the supercell is large enough, the differences between macroscopically equivalent directions or deformation modes (e.g., tension along the x, y, and z direction in the cubic systems) vanish. Although this is promising, the idea is not in line with the original purpose of SQS which was to simulate random alloys with as small supercells as possible. Moakher and Norris 14 provided a rigorous mathematical theory on how to project a tensor of elastic constants with an arbitrary symmetry onto a tensor with a desired crystallographic symmetry. This has been applied to the cubic Ti 1−x Al x N system 15,16 with a satisfactory agreement to available experimental data, however still requiring supercells with around 100 atoms and averaging over crystallographically equivalent directions.
In this work we investigate a possible trade-off between the randomness and the overall effective symmetry by introducing directionally-optimised SQSs (do-SQS, detailed description is given in Section II A) with the aim that the resulting tensor of elastic constants exhibits as small as possible deviations between the equivalent elastic constants. This in turn can lead to a significant reduction of computational resources by applying only a reduced set of deformations (similar to what is done to perfectly ordered and fully symmetric compounds, e.g. Ref. 11). The second part of this work is devoted to establishing the impact of a texture on the elastic constants of the polycrystalline aggregate. This is an important step towards a quantitative comparison of theoretical and experimental data, as well as theory-guided prediction of thin-film growth directions that provide extremal mechanical properties. Additional improvements towards modelling of real materials would be finite temperature effects and inclusion of grain boundaries, neither of which is addressed here.
To assess the performance of the here developed supercells, we have chosen cubic Zr 1−x Al x N system (NaCl prototype, F m3m space group). It is an isovalent system with well investigated and widely used Ti 1−x Al x N. Compared with TiN, ZrN has a lower coefficient of friction and has been suggested to have better oxidation resistance 3,17 . Additionally, calculated elastic constants of this system have not yet been published, and experimental values are only scarce.

A. Supercells
Warren-Cowley short-range order (SRO) parameters, α j , are commonly used to quantify randomness of an atom distribution on lattice sites. For binary alloys (or pseudo-binary, e.g. where the mixing happens only on one sublattice, as in the case of Zr 1−x Al x N), they are calculated as 18 where x A and x B (x A + x B = 1) are the mole fractions of atoms A and B, respectively, N is the number of sites in the supercell, M j is the site coordination in the jthneighbour distance, d j , and N j AB is the total number of {A, B} pairs of atoms separated by the d j (number of A-B bonds of length d j ). This definition implies that α j > 0 and α j < 0 correspond to tendency for clustering and ordering, respectively, while α j = 0 describes an ideal statistically random alloy. When constructing special quasi-random structures (SQSs), the aim is to minimise values |α j | for several first coordination shells (typically between 5 and 7). Tasnádi et al. 16 recently concluded that relatively large (around 100 atoms and more) supercells are needed to accurately describe the elastic response of a cubic Ti 0.5 Al 0.5 N. Nevertheless, somewhat smaller cells with 64 atoms and overall cubic shape do performed with an acceptable accuracy, too 16 . Moreover, we have applied a following additional constrain during the SQS generation: the number of bonds, N j AB , is divided into three subsets N j AB,x , N j AB,y , and N j AB,z , depending on which projection of the vector AB into x, y, and z direction is the longest (Fig. 1). Since the three directions x, y, and z are crystallographically equivalent in the cubic systems, the projected SRO parameters are calculated as (2) This way, the number of A-B bonds is optimised with respect to the three equivalent directions. We applied this requirement also to the A-A and B-B bonds. The resulting supercells, hereafter called directionally-optimised SQSs (do-SQSs), are summarised in Table II. They were generated using a script which randomly distributes atoms A and B (considering a required chemical composition) on the (sub)lattice, hence providing a large ensemble of various atomic arrangements. For every one of them, projected SROs α j,ξ up to j = 5 were evaluated, and a supercell with α j,ξ closest to 0, i.e. an ideal solid solution, was chosen. The projected SROs of the resulting supercells are listed in Table III. It is worth noting that the compositions x = 0.125 (x = 0.875) and x = 0.375 (= 0.625) are worse optimised than the other two x = 0.25 (x = 0.75) and x = 0.5, a behaviour consistent with the analysis of standard SQSs reported in Ref. 19. Finally, although the do-SQS method is developed here for high symmetry cubic systems, it can be applied also to other crystallographic classes by requesting the number of different bonds to be as similar as possible along equivalent directions.

B. Elastic properties
The single crystal elastic constants were obtained using the total energy method, discussed in detail in Ref. 11. The applied deformation matrices were (3) The volumetric density of the total energy increase after application of such deformation matrix is, in the first approximation, a quadratic function of δ The coefficient A depends on the applied deformation matrix A. For example, for A1 it is equal to 1 2 C P 1 11 . By changing the position of the non-zero component δ in the matrix A1 to position (2, 2) and (3, 3), the coefficient A changes to 1 2 C P 1 22 and 1 2 C P 1 33 , respectively. By fitting the U A (δ) curves for various deformation matrices (coefficients A are listed in Table I), a full tensors C P 1 ij of single crystal elastic constants was obtained for every composition. In doing so, the off-diagonal components for rows i = 4, 5, 6 and for columns j = 4, 5, 6 in C P 1 ij were fixed to 0.
Finally, application of the symmetry-based projection technique 14 as described in Ref. 16 yields the three cubic elastic constants C 11 , C 12 and C 44 : The anisotropicity of the material is quantified using Zener's anisotropy ratio, A: The stress-strain method 20 was used to confirm the total energy calculations of elastic constants. Additional
Orientation distribution function (ODF) is a convenient way to quantify the texture of polycrystals. ODF(α, β, γ) is a function of three Euler angles, α, β, γ, and for particular values it gives a fraction of grains with that orientation 21 . The Voigt (constant strain in all grains) and Reuss (constant stress in all grains) polycrystalline averages of elastic constants are defined as 22 C ijkl (α, β, γ) and S ijkl (α, β, γ) in the above expressions are stiffness and compliance tensors, respectively, in a coordinate frame rotated by angles α, β, and γ with respect to the reference coordinate frame, in which both ODF and the polycrystalline elastic constants are defined. In fact, Voigt and Reuss elastic averages represents elastic response of a polycrystalline material with grain boundaries oriented parallel and perpendicular with respect to the applied stress direction. A single-valued polycrystalline elastic properties are obtained from Hill's average A commercial package LaboTex 23 was used to generate the ODFs describing cubic fibre texture with preferred 1 0 0 or 1 1 1 orientations, different sharpness of the distribution (quantified by the full-width at half maximum (FWHM) of the distribution 22 ), and varying isotropic fraction.

C. First principle calculations
The quantum mechanical calculation within the framework of Density Functional Theory (DFT) were performed using Vienna Ab initio Simulation Package (VASP) 24,25 .
The exchange and correlation effects were treated using generalised gradient approximation (GGA) as parametrised by Perdew, Burke and Ernzerhof 26 and implemented in projector augmented wave pseudopotentials 27,28 . We used plane wave cut-off of 700 eV (500 eV) and 7 × 7 × 7 (6 × 6 × 6) Monkhorts-Pack k-point mesh for the 64 (96) atom supercells, yielding the total energy accuracy in the order of meV/at. The slightly different parameter sets are a consequence of combining results of two research groups; nevertheless, additional tests revealed that the changes in elastic constants induced by increasing the plane-wave cut-off energy from 500 eV to 700 eV and/or altering the k-mesh are not larger than a few GPa. All supercells were fully structurally optimised yielding energies and lattice parameters as discussed in Ref. 29.

A. Single crystal elastic constants
The single crystal elastic constants calculated using the total energy method as a function of the composition are shown in Fig. 2. C 11 , describing the uni-axial elastic response, decreases from 522 GPa for ZrN to 377 GPa for Zr 0.25 Al 0.75 N, and then increases again to 421 GPa for pure cubic AlN. On contrary, off-diagonal shear-related components C 12 and C 44 increase with the AlN mole fraction from 118 GPa (ZrN) to 165 GPa (AlN), and from 105 GPa (ZrN) to 306 GPa (AlN), respectively. As a results, also the Zener's anisotropy ratio, A, monotonically increases with the AlN mole fraction. This corresponds to a qualitative change of the directional Young's modulus distribution: the stiffest direction is 1 0 0 for ZrN and 1 1 1 for AlN 11 .
The error bars show standard deviation as obtained by averaging the three elastic constants equivalent for perfect cubic material 16 . The relative error is below 3% for C 11 , around 5% for C 12 and below 7% for C 44 (with the exception of Zr 0.5 Al 0.5 N where it is 11%). Consequently, one can in principle rely on values obtained only for deformation e.g. along the x direction to get an estimate for the C 11 with an accuracy better than 3% (a value usually regarded as an acceptable deviation between theory and experiment due to various exchange-correlation effects, temperature of measurement/calculation, material quality, etc. 30 ).
The Zener's anisotropy ratio reaches 1 for Zr 0.5 Al 0.5 N, implying that for this composition the alloy should have an isotropic elastic response independent of texture or crystallite orientation. This idea is further confirmed by evaluating the polycrystalline elastic constants in the next sections.

B. Elastic response of isotropic polycrystal aggregates
We begin with presenting the polycrystalline properties of isotropic aggregates, i.e. when all crystal orientations are equally probable. Figure 3 shows the compositional dependence of bulk modulus, B, Young's modulus, E, and shear modulus, G. The spread of the data as shown by the shaded area corresponds to the Voigt (upper) and Reuss (lower) limits. In the case of isotropic aggregates of cubic materials, i.e. when the ODF is a constant function, the Eqs. 9 and 10 simplify to 11 where α = V or R and S ij are elastic compliances corresponding to the elastic constants C ij 31 . Bulk modulus changes only a little with the composition of Zr 1−x Al x N, and the alloy is predicted to be some 10% softer than the binary ZrN (B = 244 GPa) and AlN (B = 250 GPa). E and G exhibit the same behaviour, being almost constant for AlN mole fractions up to x ≈ 0.6, and only then significantly rising. The range between Voigt and Reuss limits is largest for the binary nitrides, and becomes almost zero for x ≈ 0.4−0.5. Hence for these compositions, the microstructure (lamellar or columnar) does not play a role for the resulting elastic response, as has been already stated in the previous subsection based on the Zener's anisotropy ratio.

C. Influence of fibre texture
The elastic behaviour of real materials with always unique microstructure is, however, different from that of isotropic polycrystalline aggregates. Hard ceramic coatings typically exhibit a 1 0 0 or 1 1 1 fibre textures 32 , with a fibre axis oriented perpendicular to the substrate surface, which usually develops due to the minimization of the strain energy during nucleation or as a result of the surface/interface energy minimization, respectively. Hence we have investigated the influence of a particular fibre texture containing a certain fraction of isotropic background (isotropic aggregate of grains) on the elastic Young's modulus in a direction perpendicular to the film surface (z direction). Sharpness of the texture (width of the grain orientation distribution along the preferred orientation) is quantified by the full width at half maximum (FWHM) parameter 22 . Fig. 4 shows representative results of the Hill's average of Young's modulus in the z direction.
Increase of the FWHM leads to decrease (increase) of the Young's modulus in the z-direction for ZrN with the 1 0 0 ( 1 1 1 ) texture, as the contribution from softer (stiffer) oriented grains increases. Similar but opposite trend is obtained also for AlN. Increase of the isotropic content in the texture has qualitatively the same impact as increasing FWHM (decreasing sharpness of the texture). The Young's modulus values are independent of the FWHM parameter for 100% isotropic content, since in these cases no texture fibre is present. Hence also the values in the 1 0 0 and 1 1 1 plots are the same for 100% isotropic texture, and correspond to the values shown in Fig. 3.
The shape of the Young's modulus profile changes continuously with the composition. It is worth noting that for x = 0.5, almost flat dependence is obtained. This underlines our earlier estimates that the elastic response of a Zr 1−x Al x N solid solution around this composition is isotropic and independent of the texture (assuming the grain boundaries influence is negligible). To illustrate the compositional dependence further, we plotted in Fig. 5 the spread between the single crystal Young's moduli in the 1 0 0 and 1 1 1 directions. It follows that the elastic response is most sensitive to the texture for the binary ZrN and AlN. Both theory 3,29 and experiment [33][34][35] suggest that the maximum AlN mole fraction in a supersaturated cubic single phase is about x ≈ 0.4. Consequently it can be stated, that addition of Al into ZrN up to its solubility limit results in the solid solution becoming steadily more elastically isotropic, hence decreasing the impact of film microstructure on the elastic behaviour. with the previously discussed elastic constants in Fig. 6. It follows that there is up to ≈ 13% difference in the C 11 values, while the C 12 and C 44 are practically unchanged. As a consequence, the SQS-based calculation predicts Zener's anisotropy ratio A to be significantly higher than the value based on the do-SQS supercell. A test calculation using do-SQS and the stress-strain method yields the same result as the total energy approach applied to the do-SQS. Another test calculation using an ordinary 64 atom SQS cell yielded C 11 ≈ 2% smaller than the corresponding do-SQS value. It can be therefore concluded that the discrepancies shown in Fig. 6 originate from the different supercell sizes (and shapes). This is similar to the behaviour of an isovalent system Ti 1−x Al x N, where C 11 from 64 atom supercell is about 7.5% smaller than the corresponding value calculated using a 96 atom supercell.
Moreover, the shown 96 atom SQS data are C P 1 11 , C P 1 12 , and C P 1 44 instead of the projected elastic constants C 11 , C 12 , and C 44 (for which 9 elastic constants would be needed, see Eqs. 5-7). As shown in the Ref. 16, when only C P 1 11 , C P 1 12 , and C P 1 44 are used to calculate A (a value labelled A x ), A x is overestimated by ≈ 18% with respect to the value A based on the project elastic constants for x = 0.5. The other set of elastic constants, C P 1 22 , C P 1 23 , and C P 1 55 (A y ), and C P 1 33 , C P 1 13 , and C P 1 66 (A z ) underestimate A by 9% and 7%, respectively. On the contrary, do-SQS results in a significantly reduced the spread of the Zener's ratio. For example, for x = 0.5 our data yield A x /A = 1.026, A x /A = 0.997, and A x /A = 0.978. A similar trend was observed also for the standard 64 atom SQS in Ref. 16.
The supercell size together with using/not using the projection technique for C ij add up, and cause the Zener's anisotropy ratio, A, of the 96 atom supercell increase steeper (than the do-SQS value) for low AlN mole fractions, yielding x ≈ 0.35 as the composition with the nearly-isotropic response. As a consequence of the overestimated A for the 96 atom SQS, the isotropic concentration is expected to be shifted to higher Al concentrations in Fig. 6, hence get closer to the do-SQS predictions, when the corrected projected C ij values are used. It can be therefore concluded, that the here proposed do-SQS cells are more appropriate for a direct estimation of C ij (i.e., without projecting the 9 P 1 values on the 3 cubic C 11 , C 12 , and C 44 ) than the 96 atom SQS supercell.
Focusing on the single crystal elastic constants for compositions x < 0.4 (i.e. those experimentally accessible in a single phase), the relative error from averaging C 11 from deformations along the x, y, and z directions is maximum ≈ 1%, while it is below ≈ 6% for C 12 and C 44 36 . When such accuracy is acceptable, the here proposed do-SQSs can be used to perform only one of the three symmetry equivalent deformations (i.e., x, y, or z for the C 11 ) to obtain the respective elastic constant of the alloy.
The off-diagonal components in rows i = 4, 5, 6 and columns j = 4, 5, 6 set to 0 GPa during the fitting procedure of the total energy method were confirmed to be negligibly small (below 3 GPa) using the stress-strains method. Therefore, the used assumption does not influence our results.
As pointed out in the recent publications (see e.g. Refs. 12 and 37), the stress-strain method proofs to be more robust than the total energy approach. It is also faster as one deformation yields several linearly independent equations for obtaining the full 6×6 matrix of elastic constants. It follows, that the stress-strain method performed on a standard SQS (possibly with more atoms) should be the preferred method. When a reliable interatomic potential exists, a molecular dynamics based approach on large SQSs 13 may be an acceptable approach. However, in the case when neither stress tensor nor good non-DFT based method is available, we propose that our do-SQSs together with the total energy method are a suitable (CPU-time affordable) approach.

B. Comparison with experimental data
There are only a few experimental reports on the mechanical properties of Zr 1−x Al x N monolithic films in the literature. Lamni et al. 35 reported on magnetron sputtered thin films which keep cubic structure up to AlN mole fraction x = 0.43. The films have preferred 1 1 1 texture. The Young's modulus, as measured by nanoindentation, increases from 250 GPa for pure ZrN to 300 GPa for Zr 0.57 Al 0.43 N. The predicted Young's modulus for the preferred 1 1 1 texture increases from 307 GPa (335 GPa) for ZrN to 389 GPa (375 GPa) for Zr 0.625 Al 0.375 N for FWHM = 0.5 • (45 • ) (compare with Fig. 5). Therefore, the trend as well as the magnitude of the increase is correctly predicted by our calculations. The somewhat softer Young's modulus experimentally is likely to be related to the presence of grain boundaries in the real microstructure as well as finite temperature during experimental measurement. The same authors also published a value of nanoindentation Young's modulus 250 GPa for the as deposited Zr 0.57 Al 0.43 N 38 , which increased to 265 GPa after annealing the sample to 850 • C. Rogström et al. 34 also observed increase in Young's modulus after annealing their arc-deposited Zr 0.52 Al 0.48 N to 1400 • C, and argued that this is a consequence of increased grain size, decreased porosity and improved crystallinity of their sample. Although the grain size is not reflected in our calculations (ODF takes into account only the grain orientations, not their size nor shape), its effect, i.e. decreasing the volume fraction of the soft grain boundaries, can be intuitively foreseen. It therefore further underlines the importance of the grain boundaries (as well as the amorphous matrix present in the case of Zr 1−x Al x N) on the mechanical properties of nanocrystalline materials. This topic cannot be easily handled by the means of Density Functional Theory itself, and would require use of multi-method scale-bridging techniques (see e.g. 39,40 ).
Nevertheless, the influence of the grain boundary fraction has been experimentally proven by our results for reactively prepared Zr 0.65 Al 0.35 N and non-reactively prepared Zr 0.68 Al 0.32 N coatings possessing a single phase cubic structure and a mixed 1 1 1 -2 0 0 orientation. The reactively prepared coatings have grain sizes of 8 nm and Young's moduli of 347 GPa, whereas the non-reactively prepared coatings have grain sizes of 36 nm and also larger Young's moduli of 398 GPa 41 . The latter is in almost perfect agreement to our calculations for Zr 0.625 Al 0.375 N with FWHM = 0.5 • yielding 389 GPa.

V. CONCLUSIONS
The here proposed directionally-optimised SQSs seem to be a reasonable alternative to large standard SQSs for the estimation of elastic properties of alloys, in particular for systems with cubic symmetry. When well optimised, they can provide accurate single crystal elastic constants while significantly reducing the number of calculations needed by omitting the need for symmetry-based projection of C ij .
Based on the calculated single crystal elastic constants we propose that Zr 1−x Al x N with AlN mole fraction x ≈ 0.4-0.5 exhibits isotropic elastic behaviour. This in particular means, that any polycrystal with this composition will have texture-independent Young's modulus. This hypothesis has been supported by explicitly evaluating the compositional dependence of Young's modulus on fibre texture orientation, its sharpness and the amount of isotropic background. The comparison with experimental data showed decent agreement with our theoretically predicted values. The small discrepancy is ascribed mainly to the influence of grain boundaries. This phenomenon, however, goes beyond the capabilities of DFT and requires a multi-scale/multi-method approach.