Doping effects on mechanical and thermodynamic properties of zirconium carbide systems: a first-principles study

Zirconium carbide (ZrCx) is an important high temperature structural material, whose wide engineering applications are limited by carbon vacancies. Doping various impurity elements (O, B, etc) into ZrCx may lead to a significant change in its mechanical properties and thermodynamic properties behaviors. In this paper, based on the density functional theory, the effects of carbon vacancy contents and dopant on mechanical properties and deformation behaviors of zirconium carbide were discussed. With the increase of the carbon vacancy contents, the Young’s modulus, bulk modulus, and shear modulus decrease gradually. When the tensile strain is greater than 0.4, ZrC0.75 has stronger plasticity than ZrC0.875, ZrC0.9375 and ZrC. Furthermore, the mechanical properties of ZrC, ZrC0.75O0.25, ZrC0.75B0.25 and ZrC0.75 were studied. Compared with ZrC0.75, the mechanical properties of ZrC0.75O0.25 and ZrC0.75B0.25 are improved, and the mechanical properties of the systems are improved the most by doping O atoms. Based on the quasi-harmonic approximation, the influence of doping atoms on thermodynamic properties of ZrC0.75O0.25, ZrC0.75B0.25 and ZrC0.75 was also investigated. Doping O and B atoms in ZrC0.75 can improve the thermal conductivity at high temperature, and ZrC0.75B0.25 has the highest thermal conductivity. The results also show that the thermal properties of ZrC0.75 can be improved by doping O and B atoms. With the increase of temperature, ZrC0.75O0.25 has the largest thermal expansion.


Introduction
As non-stoichiometric compounds, zirconium carbide (ZrC x ) systems have high temperature durability and chemical stability and can be used in extreme environments [1][2][3]. Ultra-high-temperature ceramic material comprising zirconium carbide is particularly outstanding and is the most promising candidate for extreme environments [4,5]. Because of its excellent thermal conductivity, high-temperature durability and low density, zirconium carbide can be widely used in cutting tools and aerospace [6][7][8][9].
The properties of ZrC x have been investigated via experimental methods [10][11][12]. For example, Guillermet et al [13] have studied non-stoichiometric NaCl-type-structure ZrC x carbides, calculated the heat capacity of ZrC x with different carbon contents, and found that its heat capacity deviates at high temperatures. However, it is difficult and time-consuming to study the mechanical and thermal properties via such methods, especially at high temperatures and pressures. Due to the superiority of theoretical and computational methods, it is possible to determine the properties of ZrC x in extreme environments by first-principles calculation [14][15][16][17][18]. For example, Jhi et al [18] have evaluated the shear modulus of ZrC x by ab-initio pseudo potential calculation, showing that its excellent mechanical properties can be understood via the electronic band structure. Moreover, Kim et al [19] have calculated the mechanical properties and thermal expansion of ZrC, ZrN, and Zr(C 0.5 N 0.5 ) under high temperatures and pressures.
Although many researchers have studied the mechanical and thermodynamic properties of ZrC x [20,21], few have considered the effect of vacancy contents and dopant on these properties [22,23]. The ZrC x system is stable across a wide composition range, with the x ranging from 0.98 to 0.63 [24]. In practical applications, carbon vacancies are easily occupied by other elements, which changes the properties of ZrC x [25][26][27][28][29]. Xie et al [28] have studied the influence of different carbon vacancies on the structural stability and mechanical properties of ZrC x and found that the mechanical properties of systems with carbon vacancies are improved by doping N and O atoms. As the vacancies are easily occupied by other elements (O, B, etc), this is favorable under high temperatures, and the performance of ZrC x systems can change greatly after the vacancies are occupied. Therefore, it is necessary to study the mechanical and thermodynamic properties of doping ZrC x systems under high temperatures.
In this paper, we study the mechanical properties of ZrC x systems with different carbon vacancy contents as well as the mechanical and thermodynamic properties of doping systems. The effects of O and B doping by in situ substitution on the properties of ZrC x systems were examined via first-principles calculations. Specifically, the mechanical properties in a small elastic range were predicted from the elastic constants, elastic moduli, ductility and elastic anisotropy. In addition, the fundamental thermodynamic quantities as a function of temperature were calculated systematically by full quasi-harmonic model.

Calculation model and method
2.1. Computational details As shown in figure 1(a), zirconium carbide has a face-centered cubic structure in which Zr and C atoms occupy lattice points and octahedral interstitial sites, respectively [10]. As shown in figure 1(b), the carbon vacancy model of ZrC 0.75 was built by removing one C atom from the unit cell of ZrC, which resulted in ordered vacancies for ZrC 0.75 at the unit cell and supercell scale. There were eight vacancies in a 2×2×2 supercell of ZrC 0.75 . The ordered vacancy structures of ZrC 0.875 and ZrC 0.9375 were built by removing four and two C atoms from the 2×2×2 supercell of ZrC, respectively. The carbon vacancies were occupied by O and B atoms, figure 1(c) shows the unit cell doping with these atoms. The typical electronic configurations of Zr, C, O, and B are 4s 2 4p 6 4d 2 5s 2 , 2s 2 2p 2 , 2s 2 p 4 , and 2s 2 p 1 , respectively. Optimization and static calculations were performed using 11×11×11 K-point grids in the Monkhorst-Pack scheme. The density of the phonon states and mechanical and thermodynamic properties of the 2×2×2 supercell structure were calculated using 5×5×5 K-point grids in the Monkhorst-Pack scheme. First-principles calculations based on density functional theory were carried out using the Vienna Ab-Initio Simulation Package [30,31]. The projector augmented wave potential was employed to describe the interactions between ionic cores and valence electrons [32]. The generalized gradient approximation in Perdew-Burke-Ernzerhof scheme [33] was used for the exchange dependent function [34]. To obtain accurate results, a high energy of 600 eV was used in all calculations. The tolerance for geometry optimization was set as the difference in total energy withiń -1 10 eV atom, 6 / the maximum atomic force was less than 1×10 −2 eV Å −1 .

Computational theory of mechanical and thermodynamic properties
The elastic constants of ZrC x were C , 11 C , 12 and C . 44 The bulk modulus (B), Young's modulus (E), and shear modulus (G) were calculated by the Voigt-Reuss-Hill [35] method: where S ij is the compliance, which is the inverse matrix of elastic constant matrix C ij .
The thermodynamic properties of ZrC, ZrC 0.75 O 0.25 , ZrC 0.75 B 0.25 , and ZrC 0.75 were calculated by the full quasi-harmonic model of quasi-harmonic approximation [36], in which the Helmholtz free energy, F V T , ( )is: ) is the phonon free energy of the lattice vibration, and F V T , el ( )is the thermal electronic contribution to the free energy. In the quasi-harmonic approximation, where k B is the Boltzmann constant, and w l k V , q ( ) is the frequency of the phonon mode for wave vector k and volume V. The phonon densities of states at several volumes were calculated around the equilibrium volume, which was used to derive the thermodynamic properties of the system [37].
The calculation of the thermodynamic properties included heat capacity at constant volume C , V ( ) coefficient of volume thermal expansion (α), heat capacity at constant pressure C , P ( ) and adiabatic bulk modulus B S ( ) can be expressed as where U is the internal energy of the system, and T is the temperature, and where B T is the isothermal bulk modulus, and V is the volume affected by the temperature.

Mechanical properties of different carbon vacancy contents
The elastic constants of ZrC 0.75 , ZrC 0.875 , ZrC 0.9375 , and ZrC are shown in table 1. All four systems satisfied the mechanical stability [38]: ( ) The systems' mechanical properties are shown in table 1 as well. The lattice constant of ZrC ( = a 4.7098Å) is in good agreement with the theoretical calculations [38] and experimental value [39]. The lattice constant of ZrC 0.75 reduced to = a 4.693Å, and is in good agreement with theoretical calculation [24]. The elastic constants, bulk modulus, shear modulus, and elastic modulus of polycrystalline ZrC x decrease as the carbon vacancy contents increase. The mechanical properties of ZrC and ZrC 0.75 are consistent with the results of previous studies [24,38,39]. It is found that the bulk modulus decrease from 226.3 GPa for ZrC to 214.8 GPa for ZrC 0.9375 , 206.7 GPa for ZrC 0.875 , and 193.7 GPa for ZrC 0.75 . Correspondingly, the shear modulus and elastic modulus of ZrC 0.75 , ZrC 0.875 , ZrC 0.9375 , and ZrC also decrease as the carbon vacancy contents increase.
To demonstrate the carbon content's effect on the systems' mechanical properties, the tensile strains of ZrC 0.75 , ZrC 0.875 , ZrC 0.9375 , and ZrC along the c direction were calculated by first-principles. As shown in figure 2, the stress-strain curves of the stretch ZrC 0.75 , ZrC 0.875 , ZrC 0.9375 , and ZrC overlapped at small strains. In the elastic stage (ε0.10), the systems' stress increases linearly as the strain increases. At the same tensile strain, the stress is in the order of ZrC>ZrC 0.9375 >ZrC 0.875 >ZrC 0.75 . When the strain is 0.2, the systems' stress reaches the ultimate strength. The stress peaks of ZrC 0.75 , ZrC 0.875 , ZrC 0.9375 , and ZrC are 26.45 GPa, 28.06 GPa, 29.16 GPa, and 30.01 GPa, respectively. Compared with the ideal tensile strength of ZrC, the tensile strength decreases as the carbon vacancy contents increase. As the strain increases (ε0.20), the stress of ZrC 0.75 , ZrC 0.875 , ZrC 0.9375 , and ZrC decreases gradually. When the strain is greater than 0.4, the stress of ZrC 0.75 is the largest of the four systems, indicating that its structural failure is slower.
The bonding properties of ZrC 0.75 , ZrC 0.875 , ZrC 0.9375 , and ZrC were analyzed by the electron localization function (ELF), which can be used to analyze information on the structure of atomic shells and show the location and size of bonding and lone electron pairs [40]. The ELF values ranges from 0 to 1. A value of 1 (shown in red) reflected perfect localization connects to the covalent bond. A value of 0.5 (shown in green) reflects a metallic bond, and 0ELF<0.5 (shown in blue) reflects ionic bonding.

Models
Lattice constant(Å)    Elastic properties are fundamental for structural materials and provide a lot of information about their mechanical properties [41]. The calculated C 11 of ZrC 0.75 O 0.25 is the largest of those systems here. As C 11 characterized the anti-compression ability of the main axis, this shows that ZrC 0.75 O 0.25 is more incompressible along the main axis than the other systems. Moreover, the π-π value of ZrC 0.75 O 0.25 is enhancive, indicating that its shear resistance is stronger than that of ZrC 0.75 . In addition, the bulk modulus (B VRH ) of ZrC 0.75 O 0.25 is significantly greater than that of ZrC 0.75 B 0.25 , indicating that the former has a higher resistance to deformation. The elastic constant and moduli of ZrC 0.75 are the smallest of all the systems. It is worth mentioning that the elastic and shear moduli of ZrC 0.75 O 0.25 are in good agreement with the experimental values [28].
The enthalpy of formation DH ( )is the energy between an alloy or compound and the most stable elemental solid in the standard state and is calculated according to the following equation: where E tot is the total energy, N i is the number of atoms in a cell, E i is the energy per atom of pure i elements in the elemental solid state. It can be found from In the experiments, the hardness (H V ) of the materials was measured by the nanoindentation method of Vickers hardness testing, predicted theoretically as follows: Poisson's ratio (μ) is an indicator of the bonding nature of materials. Generally, it varies between 0.2 and 0.3 for elastic deformation [42]. If the μ is about 0.1, this indicates the material is strongly covalent; if the μ is greater than or equal to 0.25, this indicates the material is ionic [43]. The μ of ZrC, ZrC 0.75 O 0.25 , ZrC 0.75 B 0.25 , and ZrC 0.75 is 0.201, 0.206, 0.217, and 0.24, respectively. This suggest that the bonds of these materials are a mixture of covalent and ionic bonds and that they have the characteristics of covalent and ionic compounds. In addition, the μ of ZrC 0.75 O 0.25 is smaller than that of ZrC 0.75 B 0.25, which indicates that the covalent property of the Zr-O bond is stronger than that of the Zr-B bond.
Anisotropy is a factor affecting the mechanical properties of structural materials. Generally, it directly affects lattice distortion, crack initiation, and crack propagation [25]. Therefore, information about the anisotropic properties of structural materials is essential. The elastic isotropy of a cubic lattice is evaluated by the Zener ratio:

= -
The shear anisotropy ratio is = -  the beginning of the strain. Then, all structures uniformly resist the uniaxial deformation, and leading to the monotonic increase in stress. When the strain of ZrC 0.75 O 0.25 , ZrC 0.75 B 0.25 , and ZrC 0.75 is 0.16, 0.22, and 0.20, respectively, the stress peak of all systems is reached, and the localization further weakened. Simultaneously, irrevocable deformation occurred. Furthermore, the electron delocalization of ZrC 0.75 O 0.25 is the largest, which leads to the greatest downward stress trend when the strain reaches the yield strength. Remarkably, the analysis of the atomistic structure indicates that stretching destroyed the strength of the covalent bonds. As the electron localization decrease, the properties of the covalent bonds decrease. These properties are obviously weakened at a strain of 0.46, and irreparable deformation remained in the systems.

Thermodynamic properties
The phonon dispersion curves of the ZrC x systems at 0K were calculated first via a finite displacement method using Phonopy codes [45]. The harmonic approximation can describe well the phonon behavior at low temperature, but the influence of phonon anharmonic effects increase with the increase of temperature. The the stability of the ZrC, ZrC 0.75 O 0.25 , ZrC 0.75 B 0.25 and ZrC 0.75 at high temperatures was also studied because the zirconium carbide system is usually used in high temperature environments. The anharmonic phonon dispersion curves at high temperature were calculated by DynaPhoPy code [46]. As shown in figure 6, at 0K, it was found that phonon dispersion curves of ZrC, ZrC 0.75 O 0.25 , ZrC 0.75 B 0.25 and ZrC 0.75 have no imaginary frequency, and all phases are considered to be stable. The significant phonon broadening and shifts in the phonon dispersion curves at high temperatures come from many decay channels available to the phonons in the branches [47]. At 2000K, the longitudinal mode near X and M points decreases with the increase of temperature. It is observed that there are no imaginary frequencies along the Г-R high symmetry points, ZrC, ZrC 0.75 , ZrC 0.75 O 0.25 and ZrC 0.75 B 0.25 are dynamically stable. It demonstrates that the thermodynamic properties of the four systems are stable, and the thermodynamic calculations in this work are reliable. The phonon DOS curves soften with temperature. The acoustic branch, which represents the thermal properties, decreases with the increase of temperature. In the acoustic branch, the frequency of ZrC is the highest, which makes its thermal conductivity the best among the four systems. It can be found that the acoustic phonon frequency and thermal conductivity decreases after vacancies appear, and doping O and B atoms can alleviate this phenomenon.  Debye-Grüneisen model, the full quasi-harmonic model has certain advantages in the calculation of zirconium carbide systems [35].
The heat capacity is closely related to the thermodynamics and dynamics of materials and is essential for understanding vibration characteristics. The heat capacity at constant volume C V ( )and the heat capacity at constant pressure C P ( ) of ZrC, ZrC 0.75 O 0.25 , ZrC 0.75 B 0.25 and ZrC 0.75 are shown in figure 7. As the temperature increased, the C V was proportional to T 3 at low temperatures, so the C V values increased sharply with the increase in temperature. The values of C V rise slowly at higher temperatures, very close to a constant, and approximately close to Dulong-Petit limit  figure 7[a]). C P is very sensitive to temperature at low temperature, as shown in figure 7(b), is not restricted by the Dulong-Petit limit at high temperature. Due to the doping of O atoms, the C V and C P of ZrC 0.75 O 0.25 are slightly larger than those of the other systems.
As a hard coating, ZrC x is widely used in cutting tools, and its thermal expansion performance is of great value to the adhesion of coatings in the application environment. Volume thermal expansion refers to the change of the lattice volume with temperature and reflects the change of the lattice vibration frequency. The thermal expansions calculated for our systems are displayed in figure 8(a). In the low temperature environment (T<500K ), the thermal expansion of all systems is very sensitive to temperature, and increases sharply with the increase of temperature. Compared with other systems, the thermal expansion of ZrC 0.75 O 0.25 is most affected by temperature. In the high temperature stage (T>500K ), as the temperature increases, the rising trend of thermal expansion becomes slower, But the thermal expansion of ZrC 0.75 O 0.25 is always greater than the other systems. In figure 8(

Conclusion
In conclusion, the mechanical and thermodynamic properties of different zirconium carbide systems were studied by first principles method. The mechanical properties of ZrC 0.75 , ZrC 0.875 , ZrC 0.9375 and ZrC were studied. The results show that the elastic modulus, bulk modulus, and shear modulus increase with the decrease of carbon vacancy contents. When the strain is greater than 0.4, the tensile stress of ZrC 0.75 is greater than that of ZrC 0.875 , ZrC 0.9375 and ZrC. ZrC 0.75 has more metallic properties than ZrC 0.875 , ZrC 0.9375 and ZrC.