High melting point [1,2], high hardness [3], good thermal conductivity [4,5], excellent chemical stability, good electrical conductivity [5], high strength and stiffness [6,7], and resistance to oxidation even at high temperature [8] makes ZrC a suitable material for various high-temperature engineering applications [9][10][11][12][13] such as rocket nozzles, turbine blades, heat shields, cutting tools, refractory materials, armor materials, etc. Furthermore, ZrC is a common material used in the nuclear industry due to its ability to withstand high temperatures and radiation exposure, making it an ideal choice for fuel elements [14,15]. In addition to these applications, ZrC is also used to make electrodes in batteries and fuel cells [16,17]. To effectively design and optimize these applications, a thorough understanding of the underlying heat transfer mechanism on ZrC is necessary.
The study of heat transfer mechanisms in ZrC has been a topic of great interest. Many experimental thermal conductivity data of ZrC from room temperature up to its melting point (~3700 K) have been reported [4,5,[18][19][20][21][22][23]. However, the measured thermal conductivity data are scattered across the literature. For example, at room temperature, the values vary from as low as 20 W m -1 K -1 to as high as 40 W m -1 K -1 , and at 2000 K, they are scattered from 22 to 46 W m -1 K -1 . ZrC is a semimetallic material, with the total thermal conductivity (κ) arising from the combined effect of electronic (κel) and lattice thermal conductivity (κph). However, since direct measurement of lattice thermal conductivity is challenging, it is typically derived from the total thermal conductivity by subtracting the electronic contribution by the Wiedemann-Franz law [24,25]. This approach leaves some uncertainty because the Lorenz number (L) is not necessarily 𝐿 0 = 2.44×10 -8 W Ω K -2 . Moreover, κph is found to be more sensitive to external factors such as defects, grain boundaries, and impurities than κel, making it further challenging to unveil the intrinsic thermal conductivity directly from the existing experimental data.
Different theoretical studies [26][27][28][29] based on molecular dynamics (MD) and density functional theory (DFT) have been done to understand the thermal conductivity of ZrC. Crocombette [26] studied the phonon and electronic thermal conductivity of ZrC from 1000 to 3500 K using molecular dynamics (MD). They further found that the impacts of defects are strong. Although they account for the temperature effect and predict thermal conductivity correctly using MD, they did not discuss intrinsic scattering and fundamental thermal transport mechanisms. Zhou, Fahrenholtz, Graham, and Hilmas found that carbon vacancy [27] and Hf additive [28] greatly suppress the phonon thermal conductivity by using first-principles calculations. However, their results are based on the Debye-Callaway model (of which accuracy can be questionable) and ignored the higher order scattering, which is found important at elevated temperatures [30,31].
Moreover, their study is focused only on lattice thermal conductivity. Mellan, Aziz, Xia, Grau-Crespo and Duff [29] predicted the thermal conductivity of ZrC by incorporating both three phonon scatterings (3ph) and four phonon scatterings (4ph), as well as phonon renormalization for the 2 nd order force constant. However, their calculations were only based on ground state force constants (GSFC), which may not be accurate for high-temperature calculations as found in Ref [32]. In this work, we predict the thermal conductivity of ZrC by including high-order phonon scattering, lattice expansion, and temperature-dependent force constants (TDFC). We also calculate the off-diagonal (Wigner) thermal conductivity, which is found to be significant for various materials at high temperatures. Additionally, we have analyzed the impact of extrinsic factors such as grain boundary scattering, defect scattering, and impurity scattering on phonon thermal conductivity at different temperatures.
ZrC is a semi-metallic material with thermal conductivity (κ) being 𝜅 = 𝜅 𝑝ℎ + 𝜅 𝑒𝑙 .
(1)
Based on the Boltzmann transport equation (BTE), the κph can be calculated as:
where λ = (q,υ) denotes the phonon mode with wave vector q and polarization υ. 𝑐 𝜆 is the specific heat per mode, 𝑣 𝜆 𝛼 and 𝑣 𝜆 𝛽 are phonon group velocities along α and β directions. 𝜏 𝜆 𝑝ℎ is the phonon relaxation time and can be calculated using Refs. [33,34] Extrinsic scattering rates due to grain boundary (𝜏 𝑔𝑏,𝜆 -1 ) and vacancy (𝜏 𝑑,𝜆 -1 ) are calculated using the respective formulation shown in Equation (6) and Equation ( 7) respectively.
Here, 𝜆 represents phonon mode (𝐪, 𝑗) with 𝐪 and 𝑗 labeling the phonon wave vector and dispersion branch, respectively. 𝑣 𝑝ℎ represents the phonon group velocity at 𝜆, and 𝐷 𝑔𝑟𝑎𝑖𝑛 is the grain size.
Similarly, and 𝑓 𝑣 represents the concentration of the vacancy. Likewise, 𝜔 represents the angular velocities and pDOS(𝜔) represents the partial density of states of a basis atom. The coefficient 9
in Equation ( 7) accounts for the mass and bond loss associated with the defect vacancy. [35,36] More detail on the extrinsic scattering rates calculation can be found in Ref. [37] and Supplemental
Material [38].
Figure 1 shows the computational workflow of the study. The first principles calculations based on DFT are performed by using Vienna Ab initio Simulation Package (VASP) [39,40], using projected augmented wave (PAW) [41] method and Perdew-Burke-Ernzerhof (PBE) [42] functional. ZrC belongs to the Fm3m space group and exhibits a cubic FCC structure. The planeenergy cutoff is 500 eV and the energy and force convergence threshold are 2×10 -8 eV and 2×10 - 8 eV Å -1 , respectively. The relaxed lattice constant is 4.710 Å, which closely resembles the experimental value [2] of 4.694 Å. Harmonic (2 nd -order) force constants (HFC or 2FC) are extracted using Phonopy [43]. The anharmonic force constants (AFC), including 3 rd (3FC) and 4 th order (4FC), are calculated using Thirdorder and Fourthorder packages built inside ShengBTE [44], considering the 4 th and 2 nd nearest atoms, respectively. In all DFT calculations, supercells of 4×4×4 (128 atoms) are used with a k-point grid of 6×6×6. Temperature-dependent 2 nd , 3 rd , and 4 th order force constants are obtained by the TDEP method [45,46] with the input being the energy, forces, and stresses of atoms in randomly a displaced supercell lattice at provided temperature. The effect of the temperature is factored in as a thermal expansion as well as in the displacement of the generated supercells. 300 random configurations are found to be sufficient at lower temperatures (300 to 1000 K), while 500 configurations are needed at higher temperatures (1500 to 3500 K) to obtain converged force constants.
ShengBTE [44] is used to solve the BTE and calculate the 3ph and 4ph rates and κph. The convergence of ShengBTE is tested and provided in the Supplemental Material [38]. Specifically, the 3ph calculation converges at a q-mesh density of 36×36×36, and the 3ph+4ph calculation converges at 12×12×12, respectively. Different extrinsic scattering rates, such as grain boundary scattering, defect scattering, and isotope scattering are calculated using their respective formulations. Phonon-electron scattering is calculated using density functional perturbation theory (DFPT) [47][48][49] and maximally localized Wannier functions using EPW [50][51][52][53]. Finally, total scattering rates are obtained by adding all the phonon scattering rates and are used to calculate the intrinsic κph. Electronic contribution to κ is calculated using EPW [50][51][52][53] interfaced with Quantum Espresso [54,55]. The structure is relaxed with a k-mesh of 12×12×12, a kinetic energy cutoff of 200 Ry, and Gaussian smearing with a spreading parameter of 0.002 Ry. Self-consistency is achieved using Davidson iterative diagonalization with a convergence threshold of 10 -12 . Phononelectron scattering rates are calculated using a 12×12×12 k-mesh. Electron-phonon scattering rates are calculated using a coarse 6×6×6 q-mesh and k-mesh and a fine 60×60×60 q-mesh and k-mesh for convergence. The sp 3 and d entanglement is used for C and Zr, respectively. The electrons' velocities are calculated using the Wannier90 [56] under the EPW package. The Fermi window is selected so that the electron band dispersions obtained using the EPW package and DFT match with each other. Likewise, the outer window is selected so that it includes all the bands of interest.
Finally, σ, κel, and S are calculated. The thermal expansion coefficient (TEC) is calculated using quasi-harmonic approximation (QHA) with the formalism found in Ref. [59]. As shown in Fig. 2, the predicted TEC (blue curve) deviates from experimental data at high temperatures, which is commonly seen for QHA. To resolve the discrepancy, we replace the constant bulk modulus (229 GPa [4]) in the formalism with temperature-dependent (TD) bulk modulus, and the predicted TEC (black dashed curve) agrees better with experimental data. After we include TDFC, the TEC (black solid curve) at high temperature agrees even better with experimental data. It is important to note that TD bulk modulus tends to increase the slope of TEC with temperature, while TDFC tends to decrease it. Overall, these findings suggest that utilizing TD bulk modulus and TDFC is a promising approach for the accurate prediction of TEC at high temperatures.
The phonon dispersion relations of ZrC at 0, 300, 1000, and 3000 K are shown in Fig. 3, which match well with the experimental data [60][61][62]. The temperature softening effect in ZrC is not significant. The phonon dispersion calculated from TDEP at finite temperature deviates slightly from the ground state calculations. The partial density of states (PDOS) shows that the heavier metal Zr dominates the acoustic phonon of frequency (0-10 THz), while the lighter element C dominates the optical phonon of frequency (12)(13)(14)(15)(16)(17)(18)(19)(20). Since the acoustic phonon modes contribute the most to lattice heat transfer, we can expect Zr vibration dominates the heat transfer.
The 3ph, 4ph, and ph-el rates, at three different temperatures of 300, 1000, and 3000 K are presented in Fig. 4. 3ph dominates throughout all temperatures, and 4ph becomes important at high temperatures, while ph-el is insignificant at all temperatures (even though a very small portion of phonons show high ph-el rates). 4ph plays a more important role in optical phonons than acoustic phonons, similar to the findings in Si, BAs, and diamond [30]. The large 4ph rates are due to the large acoustic-optical phonon band gap, which restricts 3ph processes. Figures 4 (a-c) show the impacts of TDFC on 3ph rates at various temperatures. It is seen that TDFC decreases the 3ph rates at all temperatures. One reason is that the TD 2FC changes the phonon dispersion and reduces the 3ph phase space by making the 3ph processes' energy and momentum conservation rules harder to be satisfied, which has been widely discovered in many other materials [63,64]. We explicitly check the 3ph phase spaces of ZrC at 3000 K using GS and TD 2FC separately, as shown in Fig. 5 (a). It is seen that most modes are reduced even though some are increased. We have also compared the 3ph rates using GS and TD 2FC and found that TD 2FC gives smaller 3ph rates than GS 2FC, as shown in Fig. 5 (e). The other reason is that the 3FC is softened by temperature [32]. As seen in Fig. 5 (f), using TD 3FC gives slightly smaller 3ph rates than using GS 3FC. Figures 4 (d-f) show the impacts of TDFC on 4ph rates at various temperatures. Interestingly, acoustic 4ph rates decrease while the optical 4ph rates increase when using the TDFC instead of the GSFC. This trend becomes more prominent as the temperature increases. To understand the physical reason behind this phenomenon and tell whether this is induced by the temperature dependence of 2FC or 4FC, we have done the following two comparisons at 3000 K. First, we compare the "GS 2FC + TD 4FC" and "TD 2FC + TD 4FC" to isolate the impact of TD 2FC. As shown in Fig. 5 (g), TD 2FC results in higher 4ph rates for both acoustic and optical phonons than GS 2FC. This is because TD 2FC flattens both acoustic and optical phonon dispersions, which makes the 4ph energy conservation rule easier to be satisfied, especially for the 4ph redistribution processes (𝜆 1 + 𝜆 2 → 𝜆 3 + 𝜆 4 ) where the modes 𝜆 1 and 𝜆 3 sit on the same flat band while the modes 𝜆 2 and 𝜆 4 sit on the other same flat band. This effect is clearly shown in Figs. 5 (b-d), which shows the scattering phase space for splitting 4ph processes (𝜆 1 → 𝜆 2 + 𝜆 3 + 𝜆 4 ), labeled as "--", redistribution 4ph processes (𝜆 1 + 𝜆 2 → 𝜆 3 + 𝜆 4 ), labeled as "+-", and recombination 4ph processes (𝜆 1 + 𝜆 2 + 𝜆 3 → 𝜆 4 ), labeled as "++", respectively. It is seen that TD 2FC significantly increases the scattering phase space. This effect is stronger for optical phonons than acoustic phonons as optical branches are flattened more. Second, we compare the "TD 2FC + GS 4FC" and "TD 2FC + TD 4FC" calculations to isolate the impact of TD 4FC. As seen in Fig. 5 (h), using TD 4FC gives smaller 4ph rates for both acoustic and optical phonons. This agrees with the fact that TD 4FC reduces the scattering cross section, as found for UO2 [32]. In summary, TD 2FC and TD 4FC have competing impacts, with the former increasing 4ph rates while the latter decreases them.
The impact of TD 2FC dominates for optical phonons while the impact of TD 4FC dominates for acoustic phonons, which results in an increase of 4ph rates in optical phonon modes and a decrease in acoustic phonon modes.
Figure 6 shows the κph as a function of temperature, using various force constants. In the following, we track the changes of κph at a low (300 K) and a high temperature (3500 K) when we gradually increase the calculation comprehensivity. First, we calculate the basic 3ph thermal conductivity using GSFC, which gives 56.9 and 6.3 W m -1 K -1 at 300 and 3500 K, respectively. When we replace the GS 2FC with TD 2FC, κph increases by 9% and 33% to 62.1 and 8.4 W m -1 K -1 , respectively.
Then, we include 4ph, and κph decreases significantly by 2% and 75% to 61 and 2.1 W m -1 K -1 , respectively. After that, we replace the GS AFC with the TD AFC, and κph increases by 0% and 52% to 61 and 3.2 W m -1 K -1 , respectively. This increase in κph is attributed to the decrease in scattering cross-section with increasing temperature [32]. In the end, we add the ph-el scattering, and κph slightly decreases by 10% to 55.0 and 2.9 W m -1 K -1 at 300 and 3500 K, respectively. The inset illustrates the reduction of thermal conductivity by 4ph as a function of temperature. Further analysis on 4ph using different force constants is provided in the Supplemental Material [38]. The off-diagonal term, calculated from Wigner formalism, is not found to make a significant contribution, unlike in some other materials [65,66]. Although the Wigner contribution increases with temperature, it only reaches a maximum value of 0.2 W m -1 K -1 at 3500 K, which is much lower than the standard Peierl κph of 2.9 W m -1 K -1 at that temperature. Considering all these intrinsic effects, the κph is found to follow a temperature dependence of ̴ T -1.5 rather than ̴ T -1 . rates.
The predicted electrical conductivity (σ) as a function of temperature is shown in Fig. 7. The experimental data along with the DFT prediction from Ref. [29] are also shown for comparison.
The σ decreases with temperature monotonically as a result of the increase of electron-phonon scattering [34,[67][68][69]. The σ calculated in this study is consistent with the results of stoichiometric ZrC obtained in Ref. [26] using ab initio molecular dynamics simulations. However, it should be noted that our σ prediction is higher than most of the experimental data, particularly at low temperatures. This could be because the DFT calculations assume a perfect crystal and do not consider the effects of defects (especially the carbon vacancies) and porosity that are present in experimental samples. This idea is supported by the findings in Ref. [26], where an increase in electrical resistivity (equivalently decrease in σ) is observed with impurities and defects in the sample. Furthermore, the Lorenz number (L) using the Wiedemann-Franz law: 𝐿 = κ el /𝜎𝑇 is calculated, which is found to deviate significantly from the Sommerfeld value of 2.44×10 -8 W Ω K -2 and vary from 1.6 to 3.3×10 -8 W Ω K -2 at different temperatures. This deviation from the standard value highlights the limitations of using the Lorenz number to calculate κel from σ, as such an approach can lead to overprediction at lower temperatures and underprediction at higher temperatures. The κel is compared to κph as a function of temperature in Fig. 8. Phonons dominate thermal transport at lower temperatures whereas electrons dominate at higher temperatures, similar to results found in Refs. [26,29]. As temperature increases, e.g., from 300 K to 3500 K, κel increases from 30 W m -1 K -1 to 54.7 W m -1 K -1 , whereas κph decreases from 55.0 W m -1 K -1 to 2.9 W m -1 K -1 .
At room temperature, κph accounts for roughly 70% of κ, which declines to just about 35% and 10% at 1000 and 3500 K respectively. In Fig. 9, the predicted κ (= κph + κel) is compared to various experimental data. The blue-shaded and yellow-shaded region represents the contribution of κel and κph, respectively. We can find that the experimental data are scattered and considerably lower than our prediction. This should be due to the presence of various defects, porosity, and vacancies in the experimental samples, which significantly decrease κ. Carbon vacancy is known to inevitably present in ZrC due to the intrinsic thermodynamic instability [72,73]. This decrease is much more significant at lower temperatures compared to higher temperatures due to the suppression of κph. κel also gets suppressed due to impurities, defects, and vacancy scattering; however, the effect is smaller as found by Crocombette [26] through molecular dynamics. In that study, the predicted κ for stoichiometric ZrC was much higher than the reported experimental data. However, when some vacancy was introduced, both κel, and κph decreased, and the predicted κ matched the experimental data. It is worth noting that our κ, as well as κel, matches Crocombette's [26] results with good accuracy.
We also compare our results to those predicted in Ref. [29] by first principles. We find that they underpredict κel and overpredict κph. As a result, their total κ prediction agrees with ours at low and ultra-high temperatures, but is lower than ours in the intermediate temperature range, from 500 to 3000 K. In Ref. [29], at the temperature of 1500 K, κph is overestimated by ~27.5% (9.93→~12.67 W m -1 K -1 ), and κel is underestimated by ~38% (43.21 → ~26.78 W m -1 K -1 ), leading to κ underestimation of ~25.81% (53.15 → ~39.43 W m -1 K -1 ). At ultra-high temperatures of 3500 (in terms of %) at various grain sizes and temperatures. For instance, if the grain boundary size is 100 nm, the κph at 300, 1000, 2000, and 3000 K is roughly 55%, 76%, 87%, and 94% of the bulk values, respectively. The inset shows the grain size required to reduce κph to 80% of the bulk values. At 300 and 3500 K, such grain sizes are 445 and 25 nm, respectively. The difference in grain boundary scattering between TDFC and GSFC is not substantial. Since the grain size in the experimental samples [2,4] is in the order of µm ( ~2 to 20 µm), the grain boundary scattering might not be the primary reason behind κph suppression. The impact of carbon and zirconium vacancies on κph of ZrC is calculated and found to be strong, using perturbation theory, as shown in Fig. 11 (a). For example, at 300 K, 1% and 2% carbon vacancies can decrease κph by 39% and 56% respectively. The impact of Zr vacancy is even stronger, e.g., at 300 K, even 0.5% of Zr vacancy can lead to a decrease in κph by 80%. This can be explained by the partial density of states presented in Fig. 3. Zr has a much larger mass than C and dominates the acoustic frequency bands, which dominate the heat transfer. A small concentration of Zr vacancy can lead to a significant decrease in κph. and various experimental data collected by Jackson and Lee [4] are presented for comparison.
In this study, the thermal transport of ZrC is predicted using the first principles calculations.
Various factors affecting thermal transport, especially at high temperatures, are considered, including high-order phonon scattering, lattice expansion, TDFC, and inter-band phonon conduction. The following conclusions can be drawn. (1) TD bulk modulus and TDFC are important for the accurate prediction of TEC at high temperatures. (2) Phonons and electrons dominate the heat transfer at lower and higher temperatures respectively. κph contributes roughly 70% at room temperature and declines to about 35% and 10% when the temperature increases to 1000 and 3500 K, respectively. (3) The 4ph is important at high temperatures and reduces κph by 2%, 34%, 59%, and 76% at 300, 1000, 2000, and 3500 K, respectively. (4) Although the electronphonon scattering increases with temperature but its impact is much smaller than 3ph and 4ph scattering. (5) Interestingly, TD 2FC decreases 3ph rates but increases 4ph rates by decreasing and increasing the scattering phase spaces, respectively. For 4ph phase space, the TD 2FC flattens phonon bands, and allow more redistribution 4ph processes (𝜆 1 + 𝜆 2 → 𝜆 3 + 𝜆 4 ) to happen. (6) We confirm that the TD 3FC and 4FC decrease the phonon scattering cross-section at elevated temperatures and increase the κph significantly in metals (by 52% at 3500 K), in addition to that found in insulators [32]. (7) The combination effect of TD 2FC and TD 4FC reduces 4ph rates of acoustic modes but increases those of optical modes. (8) The maximum Wigner off-diagonal (diffuson) contribution to κph (0.2 W m -1 K -1 ) is much lower compared to standard Peierls phonon contribution κph (2.9 W m -1 K -1 at 3500 K). ( 9) The σ decreases monotonically while κel increases as temperature increases. The Lorenz number (L) deviates significantly from the Sommerfeld value of 𝐿 0 =2.44×10 -8 W Ω K -2 , highlighting the limitations of using the 𝐿 0 to calculate κel. (10) No experimental data have reached the predicted intrinsic thermal conductivity values due to the presence of inherent defects in the experimental samples. The matching of κ predicted in the literature with experimental data is due to the error-cancellation effect where the prediction underpredicts the κ and the experimental samples' defects reduce κ. Overall, our study makes a critical revisit to the thermal transport from room to ultra-high temperatures.