Thermal conductivity of complex crystals at high temperatures is critical for many applications, such as thermal barrier coatings [1,2], refractory materials [3,4], crucibles, and high-temperature thermal insulation [5,6] However, current state-of-the-art theoretical prediction based on phonongas model (PGM) [7][8][9] could not explain their intriguing thermal conductivity (κ). At intermediate temperatures, κ decays with temperature (T), being consistent with typical crystal behavior. However, at high temperatures, 𝜅 either increases or remains independent of temperature, displaying a peculiar and anomalous glass-like behavior. While the former behavior can be explained using PGM based on the Boltzmann transport equation (BTE), which primarily involves three-phonon (3ph) or four-phonon (4ph) scattering processes, the latter phenomenon remains a puzzle, presenting an unanswered question.
In this study, we take Al2O3 as an example to investigate the flatting or increasing trend of thermal conductivity at high temperatures for the following reasons. (1) Al2O3 has excellent mechanical strength, high-temperature thermal and chemical stability [10][11][12], large band gap [13], high dielectric constant [14], high melting point [15], and is widely used in many high-temperature engineering applications. (2) It does not have an electronic contribution to heat transfer, leaving the theoretical lattice thermal conductivity readily comparable to experimental data. (3) It can be grown to large-size single crystal with high purity, so grain boundary and defect scattering can be neglected. (4) Extensive experimental data [16][17][18][19][20][21][22][23][24][25] are available to validate our study. This study looks into several possible improvements to the current state-of-the-art κ prediction to explain the flatting or increasing trend of thermal conductivity in complex crystals at ultra-high temperatures. One improvement is the temperature correction on lattice structure and phononphonon scattering. The current phonon theory relies on the ground state lattice structure and force constant, where 3ph ∝ phonon population (n) ∝ T and 4ph ∝ n 2 ∝ T 2 [26,27], resulting to phonon thermal conductivity (κph) decay with power law T -1 and T -2 respectively for 3ph and 4ph scattering. However, phonon renormalization at high temperatures due to lattice expansion and temperature dependence of harmonic force constant changes 3ph and 4ph scattering phase space [28], which may lead to an increase in the κph. Recently, it has been reported that temperature dependence of anharmonic force constant reduces the scattering probability (or scattering crosssection) [28,29], which increases κph. These make the thermal conductivity flatter at higher temperatures compared to ground state calculations.
The second improvement could be the incorporation of the diffuson thermal conductivity (κdif or κc) along with κph. In PGM and BTE, the primary heat carriers are propagating phonons (particlelike phonons), which account only for the diagonal terms of the velocity operator. This is a good approximation for simple crystals [30], where the phonon branches are well separate, i.e., interband spacings are much larger than the phonon linewidths. However, this approximation fails in the complex crystals and disordered regime [31], where phonon bands are not well separated and many phonon modes might overlap with each other. In this scenario, phonon shows wavelike transport properties i.e., they can tunnel from one mode into another and conduct heat. Recently, Simoncelli et. al. [32,33] introduced the Wigner transport equation which encompasses particlelike and wavelike conduction mechanisms, providing a unified approach to heat-transport phenomena in solids, including simple crystals (where particle-like propagation dominates), glasses (where wavelike tunneling or diffusons dominates), and all intermediate cases (complex crystals where both particle-like and wavelike conduction mechanisms coexist). The significant κdif is reported for many complex crystals [34][35][36][37] and used to explain the flattening of κ at higher temperatures.
The temperature effect on both κph and κdif could also be studied using molecular dynamics (MD) simulations. In principle, MD should capture all orders of phonon-phonon interaction as well as diffuson, due to its ability to simulate the intricate behavior of individual atoms, thereby capturing all the underlying microscopic mechanisms governing lattice heat transfer. MD, while powerful, is limited by the accuracy of the interatomic potential functions they rely on. Classical potentials, based on empirical models, often involve approximations and force field parameters, leading to inherent inaccuracies in capturing the complex quantum mechanical behaviors of atoms and molecules. Machine learning interatomic potentials (MLIPs), on the other hand, have very high accuracy, comparable to that of density functional theory (DFT) calculations. In this study, we train moment tensor potential [38][39][40] (MTP) from ab initio molecular dynamics (AIMD) snapshots and used it to run Green-Kubo molecular dynamic (GKMD) [41,42] to calculate thermal conductivity (κGKMD).
The third improvement could be the incorporation of radiation thermal conduction (κrad) to overall thermal conductivity. Similar to phonon, photon also propagates inside and through the crystal and transports the heat, which contributes to apparent thermal conduction. Radiation contribution has been hypothesized as a primary reason for the increase in thermal conductivity at ultra-high temperatures in early literature [18,21,24,[43][44][45][46]. However, the radiation thermal conductivity has not been calculated rigorously from first principles. Here, we calculate the radiative properties using the Lorentz oscillator model [47,48] and estimate the radiation thermal conductivity based on the Rosseland model. This study presents the thermal transport of Al2O3 from room temperature to melting point (2200K) by incorporating the temperature-dependent force constants in κph calculation and considering the κdif and κrad. κph and κdif contributions are calculated using Wigner formalism and Green-Kubo molecular dynamics (GKMD) separately. The remaining section of the paper is structured as follows. In section II, we present the computations details and methodology used in the study. Section III and IV presents the main results and discussions respectively. Finally, section V presents the conclusions.
Figure 1 shows the computational workflow of the study. First of all, the structure is relaxed by using the Vienna Ab initio Simulation Package (VASP) [49,50], using generalized gradient approximational (GGA) [51] method and Perdew-Burke-Ernzerhof revised for solids (PBEsol) functional [52]. The plane-energy cutoff used in the calculations is 500 eV, and the energy and force convergence thresholds are 1×10 -8 eV and 1×10 -7 eV•Å -1 , respectively. Al2O3 belongs to the trigonal lattice structure system (space group 𝑅3𝑐 ̅̅̅ ), with a rhombohedral primitive unit cell For the ground state calculations, harmonic or 2 nd -order force constants (2FC) are extracted using Phonopy [55] considering the 4 th nearest neighbors. The anharmonic force constants (AFC), including 3 rd (3FC) and 4 th order (4FC), are calculated using Thirdorder and Fourthorder packages built inside ShengBTE [56], considering the 4 th and 2 nd nearest atoms, respectively. For the finite temperature calculations, the structure is expanded using thermal expansion coefficient (TEC).
The ab initio molecular dynamics (AIMD) is computed and snapshots are recorded, from which moment tensor potential (MTP) [39,40,57] based machine learning interatomic potential (MLIP) is trained. GKMD is performed by using the MLIP-LAMMPS package [39,58] to calculate κGKMD.
The temperature-dependent harmonic and anharmonic force constants are extracted using the temperature-dependent effective potential (TDEP) method [59,60] and MPT. Using the force constants, the κph and κdif are calculated by solving Wigner's formalism [32,33] under ShengBTE.
Radiative properties are calculated using the Lorentz oscillator model, from which κrad is calculated using the Rooseland model. Finally, the total thermal conductivity (κtot or κ) is calculated by summing up κph, κdif, and κrad as well as κGKMD and κrad separately. The details of each step are explained below.
The linear thermal expansion coefficient (TEC) is calculated using quasi-harmonic approximation (QHA) with the following formalism [61]:
Here, (𝐪, 𝑗) stands for a phonon mode with wavevector 𝐪 and dispersion branch 𝑗. 𝑁 𝐪 is the number of phonon q points. 𝑉 𝑐 is the volume of a primitive cell of Al2O3. 𝑘 𝐵 is the Boltzmann constant.
is the mode-dependent Grüneisen parameter. In this study, an 18×18×18 𝐪-mesh is used in the TEC calculations.
Ref. [17]. The DFT predicted data by Tohei et al. [62] and different experimental data [17,63,64] are included for comparison.
The bulk modulus at the ground state is calculated using VASP by finite difference method using
, which gives 𝐵 = 249 GPa. It matches well with the experimental data of 248.7
GPa [65] and 257 GPa [17]. The Grüneisen parameters are obtained by two methods. One is using finite difference method, i.e., 𝛾 𝐪,𝑗 = -𝑉 𝜔 𝐪,𝑗 𝜕𝜔 𝐪,𝑗
, implemented in Phonopy [55]. The other method is to use the third-order anharmonic force constant to predict 𝛾 𝜆 as implemented in ShengBTE [56].
As shown in Fig. 2, the predicted TEC using the ground state bulk modulus (blue-dash curve with Grüneisen parameter from Phonopy and black-dash curves with Grüneisen parameter from ShengBTE) match experimental data at low to medium temperatures but deviates at ultra-high temperatures. This deviation has been reported to be corrected by using the temperature-dependent bulk modulus [28]. With the temperature-dependent bulk modulus [17], the predicted TEC agrees with experimental data even at ultra-high temperatures [17,63,64]. QHA generally tends to overestimate thermal expansivity at high temperatures, which has been attributed to anharmonicity at high temperatures [66,67]. This deviation could be corrected by the phonon quasiparticle (PHQ) approach [68,69] or accurate calculation of mode Grüneisen parameters using temperaturedependent force constants [28]. TEC calculated using temperature-dependent force constants is shown in the dash-dot green line in Fig. 2, which shows that this correction is small in the case of Al2O3.
We employ a moment tensor potential (MTP) [39,40,57] based machine learning interatomic potential (MLIP) to characterize the temperature-dependent potential surface of Al2O3. The accuracy and effectiveness of this method have been demonstrated in previous studies [29,70]. A potential is trained for each temperature of study, i.e., at 500, 800, 1000, 1200, 1500, 1600, 1700, 1800, 1900, 2000, 2100, and 2200 K. The temperature interval is kept small at higher temperatures as the focus of the study is to study the flattening trend of κ at higher temperatures. The training database for each potential is prepared by AIMD with NVT ensemble for 500 steps with time step of 5 fs. The lattice structure for particular temperatures is expanded using TEC. A supercell of 3×3×3 primitive cell containing 270 atoms is used in the simulation domain. Four independent AIMDs with randomly displaced initial atomic positions are performed at each temperature to better sample the potential energy surface. Energies, forces, and stresses are recorded together with corresponding atomic configurations to construct the training and testing database. The database is separated randomly maintaining 75% (1500 snapshots) for training and 25% (500 snapshots) for testing. The initial MTP of level 22 is selected to train potential based on the accuracy and computational demand. The selected initial potential is trained for 1000 iteration steps with the minimum and maximum atomic interaction cutoff of 1.2 Å and 5.5 Å, respectively.
Once the MTP with a small error is developed, GKMD is performed by using the MLIP-LAMMPS package [39,58]. GKMD calculates the lattice thermal conductivity by integrating the heat current autocorrelation function based on the Green-Kubo formula [41,42]:
where 𝑘 𝐵 is the Boltzmann constant, 𝑉 is volume of total simulation domain, 𝑇 is temperature, 𝐽 ⃑ (𝑡) is the heat current, and the angular bracket represents an autocorrelation. In LAMMPS, the heat current vector is calculated by the energy and forces of the system, which is obtained from the MTP.
In this study, we employ a 7×4×3 supercell of the conventional cell, containing 5042 atoms. The size effect is studied. Periodic boundary conditions are implemented in all three dimensions. The time step of GKMD is set to 1 fs. First, an NVT ensemble is run for 200,000 steps (0.2 ns) to fully stabilize the temperature of the system. Then, an NVE ensemble is run for 200,000 steps (0.2 ns) to fully stabilize the system. Finally, another NVE ensemble is run for 800,000 steps with a correlation time of 200 ps, during which the heat current correlation is recorded. To mitigate the noise and intrinsic statistical error of GKMD, we conduct 16 independent runs with different initial velocities for each temperature. The ratio between the total running time and the correlation time is maintained larger than 300, which has been reported to be sufficient [71]. Since the temperatures in this study are relatively high, the difference between a classical statistic and a quantum statistic is neglected.
The temperature-dependent harmonic and anharmonic force constants are extracted using the temperature-dependent effective potential (TDEP) method [59,60] at 500, 1000, 1500, and 2000K.
The TDEP method extracts effective force constants at a certain temperature by fitting the potential energy of a series of atomic trajectory images at that temperature to the 2 nd , 3 rd , and 4 th orders.
1000 randomly displaced configurations are generated at each temperature to sample the potential surface. The forces, stress, and energies of these configurations are obtained from MTP at each temperature. The effect of the temperature is factored in by a thermal expansion, a temperaturedependent MTP trained from AIMD simulations at different temperatures, and a temperaturedependent displacement of atoms in the generated supercells.
Using the force constants, the thermal conductivity is calculated by solving the Wigner's formalism [32,33]:
𝜅 𝑝ℎ 𝛼𝛽 = ℏ 2 𝑘 𝐵 𝑇 2
𝜏 𝐪,𝑗 -1 = 𝜏 3𝑝ℎ -1 + 𝜏 4𝑝ℎ -1 + 𝜏 𝑝ℎ-𝑖𝑠𝑜 -1
.
(
Here 𝛼 and 𝛽 are cartesian directions. 𝜅 𝑝ℎ 𝛼𝛽 is the phonon particle thermal conductivity (or Peierls thermal conductivity), obtained from ShengBTE. 𝜅 𝑑𝑖𝑓 𝛼𝛽 is the diffuson thermal conductivity, calculated using an in-house code. 𝑗 and 𝑗 ′ are phonon branches, 𝜔 is the angular frequency, 𝑣 is the group velocity, and 𝜏 -1 is the phonon scattering rates, including three-phonon, four-phonon, and phonon-isotope scatterings. The convergence of ShengBTE at various q-mesh densities is tested. Specifically, the 3ph calculation converges at a q-mesh density of 18×18×18, and the 3ph + 4ph calculation converges at 6×6×6. The iterative solution to three-phonon and phonon-isotope scattering is included, as implemented in ShengBTE. Four-phonon scattering is taken at the relaxation time approximation (RTA) level. The formalism of 3ph and 4ph scattering rates can be found in Ref. [27].
Similar to phonon, photon can also transport in materials. Analogous to phonon creation and annihilation (i.e., phonon scattering events), which limit phonon mean free path, photon absorption and re-emission also limit photon mean free path inside a material. When a material's size is much larger than phonon mean free path, phonon transports diffusively, and the phonon thermal conductivity is approximately proportional to phonon mean free path. Similarly, when a material's size is much larger than its photon mean free path, the material is optically thick, and its radiation contribution to thermal conductivity is proportional to photon mean free path. Based on the Rosseland model, the radiation contribution to thermal conductivity (κrad) of an optically thick material [18,72,73] is:
where 𝑛(𝑇) and 𝛽(𝑇) are the temperature-dependent refractive index and extinction coefficient, respectively. 𝛽(𝑇) gives the attenuation of the electromagnetic waves inside the material and is the inverse of photon mean free path. 𝑛(𝑇) and 𝛽(𝑇) are given by
𝜎 𝑆𝐵 is the Stefan-Boltzmann constant. The spectral refractive index 𝑛(𝜆, 𝑇) and spectral absorption coefficient 𝛼(𝜆, 𝑇) can be calculated from the dielectric function 𝜖(𝜆, 𝑇) as:
𝜖 𝑟𝑒 is the real part of 𝜖. 𝑘(𝜆, 𝑇) is the spectral extinction coefficient given by
The dielectric function 𝜖(𝜔, 𝑇) can be predicted using the four-parameter Lorentz oscillator model [47,48] as:
where Γ is the same as phonon scattering rate, i.e., Γ = 𝜏 -1 . 𝜔 is the phonon or photon angular frequency, i.e., 𝜔 = 2𝜋𝑓. 𝑗 runs through all the infrared active transverse optical (TO) and longitudinal optical (LO) branches. 𝑖 is the imaginary unit number. 𝜖(∞) is the dielectric function at the high-frequency limit which is calculated from the density-functional-perturbation-theory (DFPT) [74]. The calculated values are 3.21 and 3.23 for ordinary and extraordinary rays, which match well with the experimental values of 3.2 and 3.1 [75]. The details of the calculation of 𝜖(𝜔, 𝑇) is similar to that of Ref. [76]. The MTP potential is trained by using the stress, forces, and energies obtained from AIMD calculations. The training and testing errors are less than 5% for each potential trained at various temperatures. To further verify the accuracy of the trained potentials, we extracted the forces and energies of the test dataset using MTP potential. Note that, this test dataset is not used to train the potential and is selected randomly. The extracted forces and energy from MTP are plotted against the forces and energies obtained from DFT calculations. As seen in Fig. 3, all the points lie along the diagonal with low root mean square error (RMSE) of 0.0426 eV/Å and 0.0664 eV for forces and energies respectively. This shows that MTP could accurately reproduce the forces and energies of the configurations and represent the actual potential surface with accuracy comparable to DFT calculation.
Using the MTP, GKMD is conducted to calculate the lattice thermal conductivity (κGKMD) as shown in Fig. 4. The average thermal conductivity is obtained by integrating the auto-correlated heat flux with correlation time. As seen in the inset of Fig. 4, the average thermal conductivity converges, suggesting that the parameters used in our calculations are appropriate. The κGKMD shows a flatting and slight increasing trend at ultra-high temperatures, which agrees reasonably well with experimental data. This further demonstrates the high accuracy of the trained MTP.
However, the experimental data shows much increasing thermal conductivity at ultra-high temperature, which could not be captured by the GKMD.
In the following, we track the changes of κ at low temperature (300 K), high temperature (1200 K), and ultra-high temperature (2200 K) when we gradually increase the calculation comprehensivity and include the contributions of phonon (κph), diffuson (κdif), and radiation (κrad).
First, we calculate the basic 3ph thermal conductivity using ground-state force constants (GSFC),
shown in the blue-dash line in Fig. 5(a), which follows κ~T -1 law. The κph is 28.59, 6.26, and 3.4 One interesting phenomenon seen in Fig. 5(a) is that the κph using GSFC and 3ph rates (dash-blue line) and κph using TDFC and 3ph+4ph rates (solid-blue line) agree with each other with reasonable accuracy at high temperatures. This is due to the competing effect of 4ph scattering and TDFC in thermal conductivity, considering 4ph rates decrease κph while TDFC increases it. In the case of Al2O3, these two opposite effects are found to cancel each other. This also introduces the possibility of an error-cancellation effect on other complex materials, where the effect of TDFC and 4ph rates cancel each other, and κph computed using GSFC and 3ph scattering matches with experimental data.
The diffuson thermal conductivity (κdif) increases with temperature at low and intermediate temperatures and becomes flat at ultrahigh temperatures as shown in Fig. 5(a). Using GSFC, the κdif is around 0.44 W•m -1 •K -1 at 300K which increases to 0.95 and 1.28 W•m -1 •K -1 at 1200 and 2200 K respectively. At high temperatures, the phonon linewidth increases, which increases the phonon tunneling probability and increases κdif. Using TDFC, the κdif decreases by 13%, 7%, and 16% to 0.38, 0.89, and 1.08 at 300,1200 and 2200 K respectively. Adding κdif to κph makes κ flatter at high temperatures as shown in dash-orange line in Fig. 7. κph + κdif matches experimental data at intermediate to high temperatures but could not explain the increase in thermal conductivity at ultra-high temperatures.
Also, κph + κdif obtained from the Wigner formalism matches with the GKMD results throughout the temperatures. This alignment underscores the capability of both GKMD and the Wigner formalism to effectively capture κph and κdif. Further, the good agreement between these two different approaches supports the Wigner formalism as well as the accuracy and consistency of our calculations. Experimental thermal conductivity data from Refs. [16][17][18][19][20][21][22][23][24][25] are also shown for comparison.
The radiation thermal conductivity (κrad) is negligibly small at room temperature, but it increases with the third power of temperature and reaches 0.22 and 1.05 Wm -1 K -1 at 1200 and 2200 K, respectively. The partial contributions of κph, κdif, and κrad are shown in Fig. 5
300, 1200, and 2200K respectively. This is shown in the solid-magenta line in Fig. 7, which matches well with the experimental data. This shows that the total thermal transport comes from the contribution of κph, κdif, and κrad. The red points in the graph show the κ obtained by summing up κGKMD and κrad, which also matches with experimental data.
The scaling laws of κph, κdif, and κrad with respect to temperature are shown in Fig. 5. For Al2O3, we find that κph decays approximately as ~T-1.14 after considering four-phonon scattering as well as finite temperature corrections to lattice constant, harmonic, and anharmonic force constants.
This is slightly different from κph~T -1.19 obtained from using ground state force constants. κdif increases roughly as ~T0.43 . κrad increases as ~T2.51 , being slightly smaller than ~T3 due to the increase of phonon linewidth with temperature, which increases photon extinction coefficient and reduces photon mean free path. These scaling laws are of interest when interpreting or predicting the thermal conductivity trends of other materials as well. Details of power law fittings can be found in the Supplemental Material [77].
Based on the diffuson theory, the phonons are characterized as either normal phonons or diffusonlike phonons based on Ioffe-Regel limit criteria [78,79]. The first criterion compares the lifetime (τ) of phonons with their period (P), stating that phonon modes with their τ smaller than P cannot be treated as particles anymore and should be treated as diffusons. In this condition, phonon modes exhibit wavelike nature, which allows them to tunnel between close eigenstates and transport heat diffusively. Figure 6(a) shows the 3ph and 4ph rates, along with P -1 (shown in the black-solid line).
As seen, both 3ph and 4ph rates, even at the higher temperature of 2000K are smaller than P -1
(equivalently τ > P) showing that all the phonon modes are normal phonons. Similarly, the second criterion states that phonon modes with mean free path (MFP) smaller than the minimum interatomic distance (𝐿 𝑎 ) become diffusons or diffuson-like phonons. As seen in Fig. 8
To look into the reason behind the radiation heat transfer, we calculate the spectral radiative properties of Al2O3 using the Lorentz oscillator model, which uses the frequency and damping of IR active phonon modes at the Γ point. Based on symmetry analysis [80][81][82], the phonon modes on Al2O3 are 2A1g + 2A1u + 3A2g + 2A2u + 4Eu+ 5Eg. Among these modes, the A2u (extraordinary ray with the electric field vector parallel to the z-axis) and Eu (ordinary ray with the electric field vector perpendicular to the z-axis) species are IR-active modes, A1g and Eg are Raman-active modes, and A2g and A1u are spectroscopically inactive. The details on determining TO and LO branch indexes of IR active modes are discussed in Ref. [76,83]. The experimental data from Querry [84] are also shown for comparison, which shows a close agreement. The reflectance is calculated from the dielectric function or refractive index as:
The spectral reflectance is shown in Fig. 9 (e), which matches well with the experimental data [48,75].
Extinction coefficient measures the attenuation of radiative waves inside the medium and is inversely correlated to photon MFP. It depends on the imaginary part of the dielectric function and is sensitive to the damping factor or phonon linewidth. Thus, temperature-dependent phonon linewidth is used to accurately calculate the extinction coefficient at higher temperatures. Note that the value for the extinction coefficient varies with the size of the material and was reported from 0.02 for bulk material [84] to 0.00008 for the thin film of 500 nm [85] at 1 μm wavelength. This results in significantly different κrad on bulk materials and thin film. The present calculation is based on the bulk material and is compared with the experimental data reported for bulk material by Querry [84]. From the spectral extinction coefficient, we can see that the Al2O3 is nearly transparent in the near-IR range, with a penetration depth of around 110 μm for ordinary rays and 60 μm for extra-ordinary rays in this range (Fig. 9(f)). In the case of thin film, the photon might pass through the material instead of interacting with it.
The present calculation is based on the Rosseland model, which assumes the materials to be optically thick, where the photon gets absorbed in the medium and re-emitted and re-absorbed. In this scenario, the medium behaves as a participating medium and leads to the radiation thermal transport, as discussed in the previous papers [18,21,24,[43][44][45][46]. The dominant radiation in the high-temperature region lies in the near-infrared range, as per Wein's law (𝜆 𝑑𝑜𝑚 𝑇 = 2898 μm ⋅ K), for which the penetration depth is in the order of ~100 μm. Since the experimental sample are in mm or higher [18,19,75,84], the optically thick medium approximation used in our calculations is justified. NAC consideration decreases κph from 30.58, 5.95, and 2.90 Wm -1 K -1 to 26.30, 5.09, and 2.48
Wm -1 K -1 at the temperature of 300, 1200, and 2200 K respectively. Whereas κdif changes from 0.48, 0.90, and 1.24 to 0.44, 0.95, and 1.28 Wm -1 K -1 respectively. Overall, the κ is overestimated if the NAC is not considered. This overestimation is due to the higher group velocities without NAC (Fig. 10(c)). The phonon scattering rates are not affected significantly. The change of group velocities is understandable since NAC causes TO-LO branch splitting, which flattens some branches in phonon dispersion as shown in Fig. 10 (c) and (d), and then reduces the group velocities.
In conclusion, this work presents the accurate first-principles prediction of the thermal conductivity of Al2O3 from room temperature to near melting point (2200 K). The lattice thermal conductivity is found to be composed of contributions of phonon, diffuson, and radiation. The nm at 300 K and 2200 K, respectively, indicating that it does not suffer from size effect for most experimental samples. (9) The photon penetration depth is about 100 nm, indicating that the ballistic effect of photon transport needs to be considered in the measurement of thermal conductivity of Al2O3 thin films when the film thickness is in the order of 100 nm at high temperatures. We hope this study has deepened the understanding of lattice thermal conductivity at ultra-high temperatures for complex crystals and will lead to more materials exploration for ultra-high temperature applications.