Introduction

Thermal conductivity and radiative properties are important material properties. Thermal conductivity (κ) is a key parameter for thermal management1 and thermoelectrics2. Thermal radiative properties, represented by the dielectric function (ε), are essential to applications including polaritonics3, thermal-photonic devices4, radiative energy converters5,6 and radiative cooling7,8,9. Both of these properties are related to phonon-phonon scattering, the process by which atom vibrations in the material interact with each other. First principles perturbation theory can be used to predict the phonon scattering rates (\({\tau }_{\lambda }^{-1}\)) and subsequently these two properties10,11,12,13,14. While the phonon-phonon interaction up to the lowest order, i.e., three-phonon (3ph) scattering, had long been considered to be adequate for the prediction of \({\tau }_{\lambda }^{-1}\)13,15,16,17, the importance of four-phonon (4ph) scattering, a higher order intrinsic scattering, was recently predicted by Feng et al.18,19 and confirmed by lattice thermal conductivity measurements of BAs20,21,22, as well as many studies on thermal conductivity23,24,25,26,27,28, Raman linewidth29,30,31, and infrared (IR) spectra7,32,33. For accurate predictions of κ and ε, especially when dealing with materials exhibiting ultra-high or ultra-low thermal conductivity, or when operating at high temperatures, or investigating optical phonons, both three-phonon and four-phonon (3ph+4ph) scattering mechanisms should be assessed.

However, the calculation of phonon scattering requires significant computational resources, especially for 4ph scattering. As a result, the first-principles predictions of κ and ε have only been done for a limited number of materials. For various important materials used in solar cells, thermal barrier coatings, and thermoelectric devices, the complex crystal structures can lead to complex phonon dispersion and numerous phonon branches. This complexity can result in billions of phonon scattering processes, making 4ph or even 3ph scattering computationally infeasible.

In this paper, we propose an approach based on sampling and maximum likelihood estimation (MLE) to reduce the computational cost of phonon scattering calculations and accelerate the predictions of κ and ε. Under the relaxation time approximation (RTA), the sampling method aims to estimate the scattering rates of 3ph and 4ph scattering for each phonon mode λ (denoted as \({\tau }_{\lambda ,3{{{\rm{ph}}}}}^{-1}\) and \({\tau }_{\lambda ,4{{{\rm{ph}}}}}^{-1}\), respectively) by sampling a subset from all phonon scattering processes. The concept is simple but works surprisingly well. After that, κ and ε are determined by using the scattering rate of all phonon modes and the IR-active phonon mode only, respectively. We demonstrate that the sampling method can significantly reduce the computational cost of the predictions of thermal conductivity and radiative properties. This allows us to revisit the thermal conductivity prediction of Si with an unprecedented q-mesh of 32 × 32 × 32, resulting in a converged thermal conductivity value that closely aligns with experimental data. The accuracy and efficiency of our approach make it ideal for high-throughput screenings of materials for thermal and optical applications.

Results

Overview of the sampling method

We start with an introduction to the rigorous calculation of \({\tau }_{\lambda }^{-1}\) and how we came up with the idea of using MLE to estimate it. The rigorous calculation of \({\tau }_{\lambda ,3{{{\rm{ph}}}}}^{-1}\) and \({\tau }_{\lambda ,4{{{\rm{ph}}}}}^{-1}\) requires an exhaustive calculation of the scattering rates of all possible 3ph and 4ph scattering processes for each phonon mode λ, denoted as \({{{\Gamma }}}_{\lambda {\lambda }^{{\prime} }{\lambda }^{{\prime}{\prime} }}^{3{{{\rm{ph}}}}}\) and \({{{\Gamma }}}_{\lambda {\lambda }^{{\prime} }{\lambda }^{{\prime}{\prime} }{\lambda }^{{\prime}{\prime} {\prime} }}^{4{{{\rm{ph}}}}}\), respectively. Then, under RTA, \({{{\Gamma }}}_{\lambda {\lambda }^{{\prime} }{\lambda }^{{\prime}{\prime} }}^{3{{{\rm{ph}}}}}\) and \({{{\Gamma }}}_{\lambda {\lambda }^{{\prime} }{\lambda }^{{\prime}{\prime} }{\lambda }^{{\prime}{\prime} {\prime} }}^{4{{{\rm{ph}}}}}\) are summed up separately to obtain \({\tau }_{\lambda ,3{{{\rm{ph}}}}}^{-1}\) and \({\tau }_{\lambda ,4{{{\rm{ph}}}}}^{-1}\) for phonon mode λ. To account for symmetry, multiple counting needs to be removed from the summation process (described in Method Section). At last, based on spectral Matthiessen’s rule11, we calculate the total phonon scattering rate: \({\tau }_{\lambda }^{-1}={\tau }_{\lambda ,3{{{\rm{ph}}}}}^{-1}+{\tau }_{\lambda ,4{{{\rm{ph}}}}}^{-1}\). Detailed formulations of 3ph and 4ph scattering rates can be found in12,18.

The huge computational cost of calculating \({\tau }_{\lambda }^{-1}\) originates from the substantial number of scattering processes involved. Calculating \({{{\Gamma }}}_{\lambda {\lambda }^{{\prime} }{\lambda }^{{\prime}{\prime} }}^{3{{{\rm{ph}}}}}\) and \({{{\Gamma }}}_{\lambda {\lambda }^{{\prime} }{\lambda }^{{\prime}{\prime} }{\lambda }^{{\prime}{\prime} {\prime} }}^{4{{{\rm{ph}}}}}\) of all of them can be expensive or even unaffordable, especially for 4ph scattering. To mitigate this computational cost, we start to consider the possibility of reducing the number of processes that we calculate. As \({\tau }_{\lambda ,3{{{\rm{ph}}}}}^{-1}\) and \({\tau }_{\lambda ,4{{{\rm{ph}}}}}^{-1}\) can be considered as the population sums of \({{{\Gamma }}}_{\lambda {\lambda }^{{\prime} }{\lambda }^{{\prime}{\prime} }}^{3{{{\rm{ph}}}}}\) and \({{{\Gamma }}}_{\lambda {\lambda }^{{\prime} }{\lambda }^{{\prime}{\prime} }{\lambda }^{{\prime}{\prime} {\prime} }}^{4{{{\rm{ph}}}}}\) under RTA, we can estimate the summations using MLE based on a randomly selected subset of scattering processes. \({\hat{\tau }}_{\lambda ,3{{{\rm{ph}}}}}^{-1}\) and \({\hat{\tau }}_{\lambda ,4{{{\rm{ph}}}}}^{-1}\) (‘\(\hat{}\)’ is used to denote estimated values) are calculated by multiplying the average values of \({{{\Gamma }}}_{\lambda {\lambda }^{{\prime} }{\lambda }^{{\prime}{\prime} }}^{3{{{\rm{ph}}}}}\) and \({{{\Gamma }}}_{\lambda {\lambda }^{{\prime} }{\lambda }^{{\prime}{\prime} }{\lambda }^{{\prime}{\prime} {\prime} }}^{4{{{\rm{ph}}}}}\) for scattering processes in the subsets and the total number of scattering processes for each phonon mode λ. In this way, only the \({{{\Gamma }}}_{\lambda {\lambda }^{{\prime} }{\lambda }^{{\prime}{\prime} }}^{3{{{\rm{ph}}}}}\) and \({{{\Gamma }}}_{\lambda {\lambda }^{{\prime} }{\lambda }^{{\prime}{\prime} }{\lambda }^{{\prime}{\prime} {\prime} }}^{4{{{\rm{ph}}}}}\) of the small subset are calculated, resulting in significant savings of computational cost. The uncertainty of our estimation can be evaluated by calculating the confidence interval using statistical information of our sample. After calculating \({\hat{\tau }}_{\lambda ,3{{{\rm{ph}}}}}^{-1}\) and \({\hat{\tau }}_{\lambda ,4{{{\rm{ph}}}}}^{-1}\), we can use them to calculate \(\hat{\kappa }\) and \(\hat{\varepsilon }\). The derivation of the maximum likelihood estimator of the population sum is shown in the Supplementary Material. A detailed explanation of the workflow can be found in the Method Section.

Maximum likelihood estimation of \({\tau }_{\lambda }^{-1}\)

We take Si as an example to show the accuracy of our approach. Figure 1a, c show the model performance at sample size n = 5 × 104 and n = 5 × 105 for 3ph and 4ph, respectively. Since the values of \({\tau }_{\lambda ,3{{{\rm{ph}}}}}^{-1}\) and \({\tau }_{\lambda ,4{{{\rm{ph}}}}}^{-1}\) vary orders of magnitude, the figures and the calculated R2 values are based on scattering rates in logarithmic scale, while the result in normal scale is shown in Supplementary Fig. 1. The high R2 values indicate that our estimation of \({\tau }_{\lambda ,3{{{\rm{ph}}}}}^{-1}\) and \({\tau }_{\lambda ,4{{{\rm{ph}}}}}^{-1}\) are quite accurate. Additionally, Fig. 1b, d shows the consistency of estimation and rigorous results across the frequency range, as well as the satisfaction of the physical scaling law \(\mathop{\lim }\limits_{\omega \to 0}{\tau }_{\lambda }^{-1}=0\) for both 3ph and 4ph, which further supports the accuracy of our estimation. Note that the average sample sizes we use for 3ph and 4ph scattering represent only 6.25% and 0.013% of the total 3ph and 4ph scattering processes for each mode, respectively. This implies that our model is expected to yield significant time savings, as we will discuss in the subsequent section.

Fig. 1: Maximum likelihood estimation of \({\tau }_{\lambda ,3{{{\rm{ph}}}}}^{-1}\) and \({\tau }_{\lambda ,4{{{\rm{ph}}}}}^{-1}\).
figure 1

a The scatter plot of \({\hat{\tau }}_{\lambda ,3{{{\rm{ph}}}}}^{-1}\) with respect to \({\tau }_{\lambda ,3{{{\rm{ph}}}}}^{-1}\). b \({\tau }_{\lambda ,3{{{\rm{ph}}}}}^{-1}\) and \({\hat{\tau }}_{\lambda ,3{{{\rm{ph}}}}}^{-1}\) as a function of phonon frequency. c The scatter plot of \({\hat{\tau }}_{\lambda ,4{{{\rm{ph}}}}}^{-1}\) with respect to \({\tau }_{\lambda ,4{{{\rm{ph}}}}}^{-1}\). d \({\tau }_{\lambda ,4{{{\rm{ph}}}}}^{-1}\) and \({\hat{\tau }}_{\lambda ,4{{{\rm{ph}}}}}^{-1}\) as a function of phonon frequency. The estimation is for all phonon modes of Si at 300 K, with sample size 5 × 104 and 5 × 105 for 3ph and 4ph, respectively. The reported R2 value is calculated from scattering rates in log scale.

Full size image

With the increase in sample size, the estimation accuracy is improved (Supplementary Fig. 2). This prompts two natural questions: What level of precision can we achieve with a given sample size? What is the appropriate sample size that yields an estimation with an acceptable level of uncertainty? While it’s possible to determine the uncertainty through multiple random samplings and estimations, it conflicts with our goal of saving computational time. Consequently, it is preferable for us to obtain the uncertainty of our estimation from a single run.

We take a phonon mode #315 as an example to illustrate how we evaluate the uncertainty, where the index of phonon is defined within ShengBTE34. The position of this phonon mode on the phonon dispersion curve is shown in Supplementary Fig. 3. According to the Central Limit Theorem, when we sample a sufficiently large number of phonon scattering processes, \({\hat{\tau }}_{\lambda ,3{{{\rm{ph}}}}}^{-1}\) and \({\hat{\tau }}_{\lambda ,4{{{\rm{ph}}}}}^{-1}\) are approximately following normal distributions with means close to \({\tau }_{\lambda ,3{{{\rm{ph}}}}}^{-1}\) and \({\tau }_{\lambda ,4{{{\rm{ph}}}}}^{-1}\), and variances decreasing as sample sizes increase (Fig. 2a, b). Based on this distribution, we can derive a confidence interval that specifies where the true values of \({\tau }_{\lambda ,3{{{\rm{ph}}}}}^{-1}\) and \({\tau }_{\lambda ,4{{{\rm{ph}}}}}^{-1}\) fall with a certain confidence level ((1 − α)100%), serving as a measure of the estimation’s uncertainty. Since each scattering process is independently and randomly sampled from the underlying distribution of \({{{\Gamma }}}_{\lambda {\lambda }^{{\prime} }{\lambda }^{{\prime}{\prime} }}^{3{{{\rm{ph}}}}}\) or \({{{\Gamma }}}_{\lambda {\lambda }^{{\prime} }{\lambda }^{{\prime}{\prime} }{\lambda }^{{\prime}{\prime} {\prime} }}^{4{{{\rm{ph}}}}}\) of phonon mode λ, our method satisfies the assumption of the Central Limit Theorem that each random variable should be independent and identically distributed (i.i.d.). We provide a rigorous derivation of the confidence interval for \({\tau }_{\lambda ,3{{{\rm{ph}}}}}^{-1}\) and \({\tau }_{\lambda ,4{{{\rm{ph}}}}}^{-1}\) in Method Section.

Fig. 2: Calculating confidence interval of \({\hat{\tau }}_{\lambda ,3{{{\rm{ph}}}}}^{-1}\) and \({\hat{\tau }}_{\lambda ,4{{{\rm{ph}}}}}^{-1}\).
figure 2

All the results of \({\hat{\tau }}_{\lambda ,3{{{\rm{ph}}}}}^{-1}\) and \({\hat{\tau }}_{\lambda ,4{{{\rm{ph}}}}}^{-1}\) in this figure are for phonon mode #315 at (0.5, 0, 0.5) of Brillouin zone and on the 3th phonon branch, with a phonon frequency of 9.62 THz. We perform random sampling and estimations 50 times for each sample size. a, b The probability density of \({\hat{\tau }}_{\lambda ,3{{{\rm{ph}}}}}^{-1}\) and \({\hat{\tau }}_{\lambda ,4{{{\rm{ph}}}}}^{-1}\) with different sample sizes. c, d The confidence interval for \({\hat{\tau }}_{\lambda ,3{{{\rm{ph}}}}}^{-1}\) and \({\hat{\tau }}_{\lambda ,4{{{\rm{ph}}}}}^{-1}\) with respect to sample sizes, respectively. The green and blue shaded areas are the upper and lower bound of confidence intervals and the green and blue lines are the means of the corresponding bounds, respectively. The orange dots are the estimation with different samples. The red lines are for the rigorously calculated scattering rate.

Full size image

We select a confidence level of 90% (α = 0.1) to calculate the confidence interval for each sample size. Then we perform 50 rounds of random samplings and estimations, and compare them with the confidence interval to verify the accuracy of our estimation (Fig. 2c, d). The results show that most of the estimated values fall within the confidence interval, indicating that our estimation is reliable. As the sample size increases, the upper and lower bounds of the confidence interval approach the true scattering rate, suggesting that the uncertainty decreases with more sampled phonon processes. With the confidence interval, we can determine whether the sample size is sufficiently large. The confidence interval goes below 10% of the scattering rate at n = 20,000 for \({\tau }_{\lambda ,3{{{\rm{ph}}}}}^{-1}\) and n = 800,000 for \({\tau }_{\lambda ,4{{{\rm{ph}}}}}^{-1}\), which indicates that these sample sizes are sufficient for phonon mode #315. We also verified our confidence interval on phonon mode #6 and #52, whose results are shown in Supplementary Fig. 4. Overall, the analytical derivation of the confidence interval serves to aid users in understanding the level of accuracy that can be achieved and in selecting an appropriate sample size.

Time-saving in phonon scattering calculations

Since we only rigorously compute a small subset of all the scattering processes, we anticipate significant computational cost savings. Figure 3 and Supplementary Fig. 5 illustrate the nearly linear relationship between the sample size and computational cost for Si and MgO, respectively. For phonon mode #315 that we discussed in the previous section, we only need to calculate 2.5% of all 3ph scattering processes to achieve an accurate \({\hat{\tau }}_{\lambda ,3{{{\rm{ph}}}}}^{-1}\). This would lead to a 90.6% reduction in computational cost compared to a full calculation. Comparing with \({\hat{\tau }}_{\lambda ,3{{{\rm{ph}}}}}^{-1}\), \({\hat{\tau }}_{\lambda ,4{{{\rm{ph}}}}}^{-1}\) requires an even smaller fraction of the sample (0.022% of all 4ph scattering processes), leading to a significant reduction in computational cost (99.97%). These substantial computational savings will significantly accelerate the calculation of \(\hat{\kappa }\) and \(\hat{\varepsilon }\), which we shall now discuss in the following two sections.

Fig. 3: Average computational cost of estimating scattering rate with respect to sample sizes for Si.
figure 3

a \({\hat{\tau }}_{\lambda ,3{{{\rm{ph}}}}}^{-1}\), b \({\hat{\tau }}_{\lambda ,4{{{\rm{ph}}}}}^{-1}\). For each sample size, sampling and estimation are repeated 50 times, and the CPU time is represented with error bars based on one standard deviation. The upper x-axis displays the ratio of the sample size to the total number of phonon scattering processes, while the right y-axis displays the ratio of CPU time for the sampling method compared to the rigorous calculation.

Full size image

The estimation of lattice thermal conductivity

To predict \({\hat{\kappa }}_{3{{{\rm{ph}}}}+4{{{\rm{ph}}}}}\), we first obtain samples of 3ph and 4ph scattering processes for each phonon mode and calculate \({\hat{\tau }}_{\lambda ,3{{{\rm{ph}}}}}^{-1}\) and \({\hat{\tau }}_{\lambda ,4{{{\rm{ph}}}}}^{-1}\). \({\tau }_{\lambda }^{-1}\) is then determined using the spectral Matthiessen’s rule11: \({\tau }_{\lambda }^{-1}={\tau }_{\lambda ,3{{{\rm{ph}}}}}^{-1}+{\tau }_{\lambda ,4{{{\rm{ph}}}}}^{-1}\). Finally, we calculate \({\hat{\kappa }}_{3{{{\rm{ph}}}}+4{{{\rm{ph}}}}}\) considering the spectral contribution of every phonon mode. When obtaining \({\hat{\kappa }}_{3{{{\rm{ph}}}}}\), only 3ph scattering is considered and \({\tau }_{\lambda }^{-1}\) contains only \({\tau }_{\lambda ,3{{{\rm{ph}}}}}^{-1}\) term.

Figure 4a, b, d, e shows the relation between the sample size for each mode and \(\hat{\kappa }\) for Si and MgO. It is worth noting that for \({\hat{\kappa }}_{3{{{\rm{ph}}}}+4{{{\rm{ph}}}}}\), the 3ph term \({\tau }_{\lambda ,3{{{\rm{ph}}}}}^{-1}\) is calculated rigorously. Only the 4ph contribution \({\tau }_{\lambda ,4{{{\rm{ph}}}}}^{-1}\) is estimated with our sampling method. This is to show the error and the saving of computational cost solely brought by the estimation of \({\tau }_{\lambda ,4{{{\rm{ph}}}}}^{-1}\). As the sample size increases, the uncertainty of estimation decreases, and the mean of \(\hat{\kappa }\) gradually approaches the true κ with a convergence rate proportional to \(\frac{1}{n}\) for all temperatures and materials. The converged sample sizes are 9 × 103 for \({\hat{\kappa }}_{3{{{\rm{ph}}}}}\) and 8 × 102 for \({\hat{\kappa }}_{3{{{\rm{ph}}}}+4{{{\rm{ph}}}}}\) for both materials, where the relative error for estimations with two consecutive sample sizes goes below 10%. Notice that \({\hat{\kappa }}_{3{{{\rm{ph}}}}}\) and \({\hat{\kappa }}_{3{{{\rm{ph}}}}+4{{{\rm{ph}}}}}\) tend to be higher than the corresponding true results, which is different from what we observe for \({\hat{\tau }}_{\lambda ,3{{{\rm{ph}}}}}^{-1}\) and \({\hat{\tau }}_{\lambda ,4{{{\rm{ph}}}}}^{-1}\) where the means of estimations remain close to the true result. This discrepancy arises because κ is inversely proportional to \({\tau }_{\lambda }^{-1}\), which amplifies the negative error and leads to an overestimation of \(\hat{\kappa }\). Besides, we observe that \(\hat{\kappa }\) converges much faster than \({\hat{\tau }}_{\lambda ,3{{{\rm{ph}}}}}^{-1}\) and \({\hat{\tau }}_{\lambda ,4{{{\rm{ph}}}}}^{-1}\) as the sample size increases. The required sample sizes for reaching convergence are much smaller than those for estimating scattering rates. This is due to the existence of error-canceling effects when accumulating the spectral contribution of κ of all phonon modes. As a result, even if there are errors in \({\hat{\tau }}_{\lambda ,3{{{\rm{ph}}}}}^{-1}\) and \({\hat{\tau }}_{\lambda ,4{{{\rm{ph}}}}}^{-1}\), the final results of \(\hat{\kappa }\) would not be significantly impacted. Moreover, We find that for both Si and MgO, the converged sample size for 4ph is smaller than for 3ph. This is because both materials are primarily dominated by 3ph scattering, and the error of \({\hat{\tau }}_{\lambda ,4{{{\rm{ph}}}}}^{-1}\) only contributes a small portion to the error of total thermal resistance, which suggests that a small sample size is enough to reach convergence for the estimation of κ3ph+4ph. To further demonstrate the performance on materials with strong 4ph effects, different symmetry and anisotropic thermal conductivity, we estimated κ3ph and κ3ph+4ph of LiCoO2, which is shown in Fig. 5. The converged sample sizes are 2 × 103 and 4 × 105 for 3ph and 3ph + 4ph, respectively, with relative errors of less than 10%.

Fig. 4: Estimation of κ with the sampling method.
figure 4

Figures a and d show \({\hat{\kappa }}^{3ph}\) and figures b and e show \({\hat{\kappa }}^{3ph+4ph}\) as functions of sample size. Subplots a and b are for Si and subplots d and e are for MgO. For each sample size, sampling and estimation are repeated 50 times and the error bar is based on one standard deviation. Figures c and f show the convergence test of \({\hat{\kappa }}^{3ph+4ph}\) of Si and MgO at 300 K, respectively.

Full size image
Fig. 5: Estimation of the anisotropic κ of LiCoO2 with the sampling method.
figure 5

The figures show \({\hat{\kappa }}^{3ph}\) (a) and \({\hat{\kappa }}^{3ph+4ph}\) (b) with respect to sample size. The calculation is performed at 300 K. For each sample size, sampling and estimation are repeated 50 times and the error bar is based on one standard deviation.

Full size image

When calculating the thermal conductivity, the first Brillouin zone should be discretized into q-mesh which should be determined based on a rigorous convergence test. While a large and converged q-mesh can used in 3ph scattering calculation, many studies employ a small q-mesh in 4ph scattering calculation due to the large computational cost35,36,37, which may not be converged. To reduce the computational cost, a smaller broadening factor (scalebroad in ShengBTE34) was usually adopted in calculations38,39,40,41,42. However, this method is known to underpredict the scattering rates. To get the result under high q-mesh, some meaningful attempts have been made by calculating thermal conductivity using several small q-meshes and extrapolating the result to the large q-mesh34,43,44. However, it lacks physical evidence to support a specific mathematical relationship between thermal conductivity and the q-mesh size. Since our method dramatically reduces the computational cost, it is now possible to study the 4ph scattering with a dense q-mesh. Take Si as an example. In previous studies, while the 3ph calculation is based on a dense q-mesh (around 28 × 28 × 28), the rigorous 3ph + 4ph calculation can only be carried out on 16 × 16 × 16 q-mesh at most due to the limit of computational power18,19,43,45. We do a convergence test on 3ph + 4ph scattering at 300 K for Si using our sampling method, which is shown in Fig. 4c. The calculation of 3ph scattering rates is based on the iterative scheme and the calculation of 4ph scattering rates is based on RTA with a sample size 106 to ensure the accuracy. The uncertainty of our method is only 0.2% while the computational cost is dramatically reduced. We found that κ3ph+4ph converged at 32 × 32 × 32 q-mesh (133 W m−1 K−1), which is much higher than the result from 16 × 16 × 16 q-mesh (114 W m−1 K−1) and is closer to the experimental results (130-150 W m−1 K−1 46,47,48,49). For MgO, the 3ph + 4ph thermal conductivity converges at 15 × 15 × 15 q-mesh (Fig. 4f), which is the same as the mesh used in previous study50. In this paper, since we need to do a rigorous 3ph + 4ph calculation to get a ground truth value to verify our estimations, the q-mesh is set to 16 × 16 × 16 for Si when illustrating the model performance of estimating \({\tau }_{\lambda }^{-1}\), the confidence interval and the time-saving. The estimation of thermal conductivity of Si that we have shown in Fig. 4 is based on 32 × 32 × 32 q-mesh, while the result with 16 × 16 × 16 q-mesh is shown in Supplementary Fig. 6.

To evaluate the time-saving effect, we test our method on Si and MgO at 300 K for both 3ph and 4ph calculations. For a sample size that reaches a relative uncertainty of less than 10%, our method achieves three to four orders of magnitude acceleration for 4ph scattering calculations (Fig. 6). The reduction in computational time is significant.

Fig. 6: Time-saving of predicting lattice thermal conductivity.
figure 6

The average CPU times for calculating \({\hat{\kappa }}_{3{{{\rm{ph}}}}}\) and \({\hat{\kappa }}_{3{{{\rm{ph}}}}+4{{{\rm{ph}}}}}\) are based on 50 runs at their respective converged sample sizes. The speed-ups are 47 × for Si with 3ph scattering, 23,000 × for Si with 3ph+4ph scattering, 24 × for MgO with 3ph scattering and 44,000 × for MgO with 3ph+4ph scattering, 65 × for LiCoO2 with 3ph scattering and 1.700 × for LiCoO2 with 3ph+4ph scattering.

Full size image

The estimation of thermal radiative properties

For the prediction of ε, \({\tau }_{\lambda }^{-1}\) of IR-active phonon modes are required, which are used as the damping factor (γ) in the Lorentz oscillator model to determine ε. Since Si has no IR-active phonon modes, we take MgO, a polar material, to demonstrate the performance of our method.

Figure 7a, b show the relation between \(\hat{\gamma }\) for both LO (longitudinal optical) and TO (transverse optical) modes and sample size. Again, for \({\hat{\gamma }}_{{{{\rm{LO,3ph+4ph}}}}}\) and \({\hat{\gamma }}_{{{{\rm{TO,3ph+4ph}}}}}\), the 3ph scattering part is calculated rigorously and 4ph scattering part is estimated with our sampling method. Similar to \({\hat{\tau }}_{\lambda ,3{{{\rm{ph}}}}}^{-1}\) and \({\hat{\tau }}_{\lambda ,4{{{\rm{ph}}}}}^{-1}\), increasing the sample size leads to a decrease in the uncertainty of damping factor. We then employ the Lorentz oscillator model with these damping factors to estimate \({\hat{\varepsilon }}_{3{{{\rm{ph}}}}}\) and \({\hat{\varepsilon }}_{3{{{\rm{ph}}}}+4{{{\rm{ph}}}}}\), as shown in Fig. 7c, d. To emphasize the uncertainty brought by the sampling method, we only show the wavelength range near the resonance frequency. Our sampling method gives more and more accurate estimations as the sample size increases. The predictions of a wider wavelength range and under different temperatures are shown in Supplementary Figs. 7–9, together with other thermal radiative properties including refractive index (n) and normal reflectance from a material-air interface (R). We take 6 × 104 for 3ph scattering and 2 × 104 for 4ph scattering as sufficiently large sample sizes where the relative error of the damping factors for both LO and TO modes go below 5%. Regarding the time-saving effect, we can accelerate the calculation by 4 × and 34,000 × for 3ph and 3ph+4ph scattering calculations, respectively.

Fig. 7: Estimation of dielectric function for MgO at 300 K.
figure 7

a \({\hat{\gamma }}_{{{{\rm{LO}}}},3ph}\) and \({\hat{\gamma }}_{{{{\rm{TO}}}},3ph}\) with respect to different sample sizes. b \({\hat{\gamma }}_{{{{\rm{LO}}}},3ph+4ph}\) and \({\hat{\gamma }}_{{{{\rm{TO}}}},3ph+4ph}\) with respect to different sample sizes. The error bar is based on one standard deviation. c, d The dielectric function for 3ph (\({\hat{\varepsilon }}_{3ph}\)) and 3ph+4ph (\({\hat{\varepsilon }}_{3ph+4ph}\)), respectively. “Re" stands for the real part and “Im" stands for the imaginary part of the dielectric function, respectively. The rigorously calculated results are given by solid lines. The upper and lower boundaries of the shaded area are the maximum and minimum estimations of the 50 runs, so all estimations are within the shaded areas.

Full size image

Discussion

Our method is based on the RTA, which treats normal and umklapp scattering processes equally. This approximation is accurate for predicting thermal radiative properties because the optical responses only involve the excitation and decay of IR-active phonon modes, and Umklapp and normal processes are both effective decay channels for an excited (over-populated) zone-center optical phonon. Previous works have demonstrated the success of the RTA in predicting thermal radiative properties for such cases7,32. For predicting lattice thermal conductivity, the RTA is suitable for materials in which the Umklapp processes are the dominant contributors to thermal resistance.

We should highlight that the sample is not from the phonon phase space but from all possible combinations of 3ph/4ph no matter whether they obey the energy or momentum conservation or not. Since the phonon phase space calculation accounts for nearly a quarter of the total computational time (Supplementary Fig. 10), if we first determine the phonon phase space and then perform sampling within it, the time-saving effect will be at most 78%. Instead of taking this approach, we perform sampling based on all 3ph and 4ph scattering processes before the evaluation of phonon phase space. Scattering rates of scattering processes that violate the energy or momentum selection rules are set to zero. By adopting this method, we can bypass the calculation of the phonon phase space on the entire phonon dispersion but effectively only do this for the chosen sample, hence achieving a time reduction of over 99%. Furthermore, in a similar manner, we can also estimate the size of phonon phase space and the corresponding confidence interval using the sampling method, which is shown in the Supplementary Material.

In our previous study51, we introduced a machine learning-based approach that can reduce the computational cost of calculating phonon scattering by up to two orders of magnitudes. It’s important to note that this machine learning method and the sampling method proposed in this paper are fundamentally different. The machine learning method in our previous work predicts the scattering rate of each individual scattering process (\({{{\Gamma }}}_{\lambda {\lambda }^{{\prime} }{\lambda }^{{\prime}{\prime} }}^{3ph}\) and \({{{\Gamma }}}_{\lambda {\lambda }^{{\prime} }{\lambda }^{{\prime}{\prime} }{\lambda }^{{\prime}{\prime} {\prime} }}^{4ph}\)) with a machine learning model trained on a subset of scattering processes. The scattering rate of one phonon mode (\({\tau }_{\lambda ,3{{{\rm{ph}}}}}^{-1}\) and \({\tau }_{\lambda ,4{{{\rm{ph}}}}}^{-1}\)) is calculated subsequently. Theoretically, this method can work both under RTA and with the iterative scheme, since it retains all the details of phonon scattering processes. On the other hand, the sampling method does not predict \({{{\Gamma }}}_{\lambda {\lambda }^{{\prime} }{\lambda }^{{\prime}{\prime} }}^{3ph}\) and \({{{\Gamma }}}_{\lambda {\lambda }^{{\prime} }{\lambda }^{{\prime}{\prime} }{\lambda }^{{\prime}{\prime} {\prime} }}^{4ph}\), but directly estimate \({\tau }_{\lambda ,3{{{\rm{ph}}}}}^{-1}\) and \({\tau }_{\lambda ,4{{{\rm{ph}}}}}^{-1}\) based on the sample. This inherent simplicity enables the sampling method to achieve even greater computational efficiency, as it eliminates the need for a time-consuming training process. However, it does not retain detailed information about phonon scattering processes and is designed to operate exclusively under RTA. The comparison of these two methods is shown in Table 1, where we choose a sample size that leads to a similar relative error for better comparison.

Table 1 Comparison of the sampling method and machine learning method51 on Si and LiCoO2.
Full size table

For the prediction of thermal conductivity, there are mainly two computationally expensive steps: calculating force constants and calculating phonon scattering. Take Si as an example, calculating the 4th order force constants takes approximately 6000 CPU hours while computing the 4ph scattering takes around 7000 CPU hours with q-mesh 16 × 16 × 16. For larger q-mesh, the computational cost of phonon scattering would be higher. Our work aims to reduce the computational cost of the latter step, which by itself is a challenging task, while the former warrants further study in the future. There are some works that aim to accelerate the calculation of force constants such as using compressive sensing lattice dynamics52 and utilizing machine learning potentials to compute force constants53,54. Future work could be done on combining our approach with these methods to reduce the total computational cost of predicting thermal conductivity.

There is some room to further improve our method. In our study, we employ sampling with replacement during the sampling process as it is computationally faster and more memory-efficient, given the large number of phonon scattering processes involved in our analysis. With sampling with replacement, the chosen scattering process is returned to the population after being selected, allowing for the possibility of selecting the same process multiple times during the sampling. However, theoretically, sampling without replacement can provide even greater accuracy than sampling with replacement under the same sample size. By removing each selected sample from the population after selection, sampling without replacement reduces the potential for redundancy in the sampled data, leading to a more diverse and representative sample. The confidence interval can be corrected by incorporating a finite population factor55 to account for the finite size of the population. A detailed description is shown in the Supplementary Material. Future work can focus on implementing this sampling scheme. Besides, when estimating κ, we sample the same number of scattering processes from all phonon modes. However, since acoustic phonon modes have a larger contribution to κ compared to optical phonon modes, it would be more beneficial to focus on improving the accuracy of predictions for acoustic phonon modes. Thus, choosing different sample sizes for each mode based on their contributions to thermal conductivity can lead to a faster and more accurate estimation of κ. The same thing happens to the estimation of thermal radiative properties, where the convergence rate of LO and TO modes can be quite different. The model efficiency can be further improved by choosing different sample sizes for these two modes.

To summarize, this study demonstrates a significant acceleration in the prediction of thermal conductivity and radiative properties, both of which are closely linked to phonon anharmonicity. By using a maximum likelihood estimation, we have successfully accelerated the calculation by orders of magnitude. This work removes the computational barrier associated with phonon scattering calculations, allowing for calculation with a converged q-mesh as well as high-throughput screening of materials’ thermal and optical properties.

Method

Estimate \({\tau }_{\lambda }^{-1}\)

We first describe how we use MLE to estimate \({\tau }_{\lambda ,3{{{\rm{ph}}}}}^{-1}\) and \({\tau }_{\lambda ,4{{{\rm{ph}}}}}^{-1}\). The analytical equation for calculating \({\tau }_{\lambda ,3{{{\rm{ph}}}}}^{-1}\) and \({\tau }_{\lambda ,4{{{\rm{ph}}}}}^{-1}\) are:

$${\tau }_{\lambda ,3{{{\rm{ph}}}}}^{-1}={({\tau }_{\lambda ,3{{{\rm{ph}}}}}^{(+)})}^{-1}+{({\tau }_{\lambda ,3{{{\rm{ph}}}}}^{(-)})}^{-1},$$
(1)
$${\tau }_{\lambda ,4{{{\rm{ph}}}}}^{-1}={({\tau }_{\lambda ,4{{{\rm{ph}}}}}^{(++)})}^{-1}+{({\tau }_{\lambda ,4{{{\rm{ph}}}}}^{(+-)})}^{-1}+{({\tau }_{\lambda ,4{{{\rm{ph}}}}}^{(--)})}^{-1},$$
(2)

where (+) and (−) terms stands for 3ph combination and splitting processes and the (++), (+−) and (−−) terms stand for 4ph combination, redistribution and splitting processes. Under RTA, the scattering rate of each type is given by12,13,18,19:

$${({\tau }_{\lambda ,3{{{\rm{ph}}}}}^{(+)})}^{-1}=\frac{1}{{N}_{{{{\rm{grid}}}}}}\left(\mathop{\sum }\limits_{{\lambda }^{{\prime} }{\lambda }^{{\prime}{\prime} }}^{(+)}{{{\Gamma }}}_{\lambda {\lambda }^{{\prime} }{\lambda }^{{\prime}{\prime} }}^{(+)}\right)$$
(3)
$${({\tau }_{\lambda ,3{{{\rm{ph}}}}}^{(-)})}^{-1}=\frac{1}{{N}_{{{{\rm{grid}}}}}}\left(\mathop{\sum }\limits_{{\lambda }^{{\prime} }{\lambda }^{{\prime}{\prime} }}^{(-)}\frac{1}{2}{{{\Gamma }}}_{\lambda {\lambda }^{{\prime} }{\lambda }^{{\prime}{\prime} }}^{(-)}\right)$$
(4)
$${({\tau }_{\lambda ,4{{{\rm{ph}}}}}^{(++)})}^{-1}=\frac{1}{{N}_{{{{\rm{grid}}}}}}\left(\mathop{\sum }\limits_{{\lambda }^{{\prime} }{\lambda }^{{\prime}{\prime} }{\lambda }^{{\prime}{\prime} {\prime} }}^{(++)}\frac{1}{2}{{{\Gamma }}}_{\lambda {\lambda }^{{\prime} }{\lambda }^{{\prime}{\prime} }{\lambda }^{{\prime}{\prime} {\prime} }}^{(++)}\right),$$
(5)
$${({\tau }_{\lambda ,4{{{\rm{ph}}}}}^{(+-)})}^{-1}=\frac{1}{{N}_{{{{\rm{grid}}}}}}\left(\mathop{\sum }\limits_{{\lambda }^{{\prime} }{\lambda }^{{\prime}{\prime} }{\lambda }^{{\prime}{\prime} {\prime} }}^{(+-)}\frac{1}{2}{{{\Gamma }}}_{\lambda {\lambda }^{{\prime} }{\lambda }^{{\prime}{\prime} }{\lambda }^{{\prime}{\prime} {\prime} }}^{(+-)}\right),$$
(6)
$${({\tau }_{\lambda ,4{{{\rm{ph}}}}}^{(--)})}^{-1}=\frac{1}{{N}_{{{{\rm{grid}}}}}}\left(\mathop{\sum }\limits_{{\lambda }^{{\prime} }{\lambda }^{{\prime}{\prime} }{\lambda }^{{\prime}{\prime} {\prime} }}^{(--)}\frac{1}{6}{{{\Gamma }}}_{\lambda {\lambda }^{{\prime} }{\lambda }^{{\prime}{\prime} }{\lambda }^{{\prime}{\prime} {\prime} }}^{(--)}\right),$$
(7)

where Ngrid stands for the total grid of q-mesh and Γ terms stand for the scattering rates of phonon scattering processes, which can be computed using Fermi’s golden rule56 and has been described in previous literature. The fractions in the summation account for the multiple counting of phonon scattering processes originated from symmetry. When illustrating the model performance of estimating \({\tau }_{\lambda }^{-1}\), the confidence interval and the time saving, the q-mesh for Si is set to 32 × 32 × 32 for 3ph and 16 × 16 × 16 for 4ph scattering. For MgO, both 3ph and 4ph scattering are using a meth of 15 × 15 × 15. When estimating κ and/or ε, the q-mesh is set to 32 × 32 × 32 for Si, 15 × 15 × 15 for MgO and 10 × 10 × 10 for LiCoO2, respectively.

We sample nλ scattering processes from all phonon scattering processes of a phonon mode λ. The total number of phonon scattering processes Nλ is given by:

$${N}_{\lambda }=\left\{\begin{array}{rlr}{N}_{{{{\rm{bands}}}}}^{2}{N}_{{{{\rm{grid}}}}}&&3{{{\rm{ph}}}}\,{{{\rm{scattering}}}}\\ {N}_{{{{\rm{bands}}}}}^{3}{N}_{{{{\rm{grid}}}}}^{2}&&4{{{\rm{ph}}}}\,{{{\rm{scattering}}}}\end{array}\right.$$
(8)

where Nbands is the number of phonon branches.

After sampling, we use the MLE to estimate \({\tau }_{\lambda }^{-1}\), which is given by:

$${({\hat{\tau }}_{\lambda ,3{{{\rm{ph}}}}}^{(+)})}^{-1}=\frac{1}{{N}_{{{{\rm{grid}}}}}}\left(\frac{{N}_{\lambda }^{(+)}}{{n}_{\lambda }^{(+)}}\mathop{\sum }\limits_{{\lambda }_{n}^{{\prime} }{\lambda }_{n}^{{\prime}{\prime} }}^{(+)}{{{\Gamma }}}_{\lambda {\lambda }_{n}^{{\prime} }{\lambda }_{n}^{{\prime}{\prime} }}^{(+)}\right)$$
(9)
$${({\hat{\tau }}_{\lambda ,3{{{\rm{ph}}}}}^{(-)})}^{-1}=\frac{1}{{N}_{{{{\rm{grid}}}}}}\left(\frac{{N}_{\lambda }^{(-)}}{{n}_{\lambda }^{(-)}}\mathop{\sum }\limits_{{\lambda }_{n}^{{\prime} }{\lambda }_{n}^{{\prime}{\prime} }}^{(-)}\frac{1}{2}{{{\Gamma }}}_{\lambda {\lambda }_{n}^{{\prime} }{\lambda }_{n}^{{\prime}{\prime} }}^{(-)}\right)$$
(10)
$${({\hat{\tau }}_{\lambda ,4{{{\rm{ph}}}}}^{(++)})}^{-1}=\frac{1}{{N}_{{{{\rm{grid}}}}}}\left(\frac{{N}_{\lambda }^{(++)}}{{n}_{\lambda }^{(++)}}\mathop{\sum }\limits_{{\lambda }_{n}^{{\prime} }{\lambda }_{n}^{{\prime}{\prime} }{\lambda }_{n}^{{\prime}{\prime} {\prime} }}^{(++)}\frac{1}{2}{{{\Gamma }}}_{\lambda {\lambda }_{n}^{{\prime} }{\lambda }_{n}^{{\prime}{\prime} }{\lambda }_{n}^{{\prime}{\prime} {\prime} }}^{(++)}\right),$$
(11)
$${({\hat{\tau }}_{\lambda ,4{{{\rm{ph}}}}}^{(+-)})}^{-1}=\frac{1}{{N}_{{{{\rm{grid}}}}}}\left(\frac{{N}_{\lambda }^{(+-)}}{{n}_{\lambda }^{(+-)}}\mathop{\sum }\limits_{{\lambda }_{n}^{{\prime} }{\lambda }_{n}^{{\prime}{\prime} }{\lambda }_{n}^{{\prime}{\prime} {\prime} }}^{(+-)}\frac{1}{2}{{{\Gamma }}}_{\lambda {\lambda }_{n}^{{\prime} }{\lambda }_{n}^{{\prime}{\prime} }{\lambda }_{n}^{{\prime}{\prime} {\prime} }}^{(+-)}\right),$$
(12)
$${({\hat{\tau }}_{\lambda ,4{{{\rm{ph}}}}}^{(--)})}^{-1}=\frac{1}{{N}_{{{{\rm{grid}}}}}}\left(\frac{{N}_{\lambda }^{(--)}}{{n}_{\lambda }^{(--)}}\mathop{\sum }\limits_{{\lambda }_{n}^{{\prime} }{\lambda }_{n}^{{\prime}{\prime} }{\lambda }_{n}^{{\prime}{\prime} {\prime} }}^{(--)}\frac{1}{6}{{{\Gamma }}}_{\lambda {\lambda }_{n}^{{\prime} }{\lambda }_{n}^{{\prime}{\prime} }{\lambda }_{n}^{{\prime}{\prime} {\prime} }}^{(--)}\right),$$
(13)

where \({\lambda }_{n}^{{\prime} }{\lambda }_{n}^{{\prime}{\prime} }\) and \({\lambda }_{n}^{{\prime} }{\lambda }_{n}^{{\prime}{\prime} }{\lambda }_{n}^{{\prime}{\prime} {\prime} }\) stand for the sampled 3ph and 4ph scattering process for phonon mode λ, nλ stands for sample size.

Confidence interval of \({\tau }_{\lambda }^{-1}\)

We take \({({\tau }_{\lambda ,3{{{\rm{ph}}}}}^{(+)})}^{-1}\) as an example to show the derivation of the confidence interval. The mean \({\bar{{{\Gamma }}}}_{\lambda }^{(+)}\) and variance \({({s}_{\lambda }^{(+)})}^{2}\) of the sampled \({{{\Gamma }}}_{\lambda {\lambda }_{n}^{{\prime} }{\lambda }_{n}^{{\prime}{\prime} }}^{(+)}\) is given by:

$${\bar{{{\Gamma }}}}_{\lambda }^{(+)}=\frac{1}{{n}_{\lambda }^{(+)}}\mathop{\sum }\limits_{{\lambda }_{n}^{{\prime} }{\lambda }_{n}^{{\prime}{\prime} }}^{(+)}{{{\Gamma }}}_{\lambda {\lambda }_{n}^{{\prime} }{\lambda }_{n}^{{\prime}{\prime} }}^{(+)}$$
(14)
$${({s}_{\lambda }^{(+)})}^{2}=\mathop{\sum }\limits_{{\lambda }_{n}^{{\prime} }{\lambda }_{n}^{{\prime}{\prime} }}^{(+)}\frac{{\left({{{\Gamma }}}_{\lambda {\lambda }_{n}^{{\prime} }{\lambda }_{n}^{{\prime}{\prime} }}^{(+)}-{\bar{{{\Gamma }}}}_{\lambda }^{(+)}\right)}^{2}}{{n}_{\lambda }^{(+)}-1}$$
(15)

According to Central Limit Theorem, the distribution of the sample mean \({\bar{{{\Gamma }}}}_{\lambda }^{(+)}\) is approximately a normal distribution with a variance given by:

$$Var({\bar{{{\Gamma }}}}_{\lambda }^{(+)})=\frac{{({s}_{\lambda }^{(+)})}^{2}}{{n}_{\lambda }^{(+)}}$$
(16)

The variance of \({({\hat{\tau }}_{\lambda ,3{{{\rm{ph}}}}}^{(+)})}^{-1}\) is then given by:

$$Var({({\hat{\tau }}_{\lambda ,3{{{\rm{ph}}}}}^{(+)})}^{-1})=Var({N}_{\lambda }^{(+)}{\bar{{{\Gamma }}}}_{\lambda }^{(+)})={N}_{\lambda }^{(+)}\frac{{({s}_{\lambda }^{(+)})}^{2}}{{n}_{\lambda }^{(+)}}$$
(17)

As the sample variance is approximated, we use t-distribution to calculate the confidence interval of \({({\hat{\tau }}_{\lambda ,3{{{\rm{ph}}}}}^{(+)})}^{-1}\). The (1 − α)100% confidence interval for \({({\hat{\tau }}_{\lambda ,3{{{\rm{ph}}}}}^{(+)})}^{-1}\) is given by:

$${({\hat{\tau }}_{\lambda ,3{{{\rm{ph}}}}}^{(+)})}^{-1}\pm {t}_{\alpha /2,{n}_{\lambda }^{(+)}}{N}_{\lambda }^{(+)}\sqrt{{{s}_{\lambda }^{(+)}}^{2}/{n}_{\lambda }^{(+)}}$$
(18)

where \({t}_{\alpha /2,{n}_{\lambda }^{(+)}}\) is the t-value of a t-distribution with \(({n}_{\lambda }^{(+)}-1)\) degrees of freedom and significant level α/2.

Similarly, we can derive the confidence interval of \({({\hat{\tau }}_{\lambda ,3{{{\rm{ph}}}}}^{(-)})}^{-1}\), which is given by:

$${({\hat{\tau }}_{\lambda ,3{{{\rm{ph}}}}}^{(-)})}^{-1}\pm {t}_{\alpha /2,{n}_{\lambda }^{(-)}}{N}_{\lambda }^{(-)}\frac{1}{2}\sqrt{{{s}_{\lambda }^{(-)}}^{2}/{n}_{\lambda }^{(-)}}$$
(19)

We denote the half length of the confidence interval with the symbol CI and we have:

$$C{I}_{\lambda ,3{{{\rm{ph}}}}}^{(+)}={t}_{\alpha /2,{n}_{\lambda }^{(+)}}{N}_{\lambda }^{(+)}\sqrt{{{s}_{\lambda }^{(+)}}^{2}/{n}_{\lambda }^{(+)}}$$
(20)
$$C{I}_{\lambda ,3{{{\rm{ph}}}}}^{(-)}={t}_{\alpha /2,{n}_{\lambda }^{(-)}}{N}_{\lambda }^{(-)}\frac{1}{2}\sqrt{{{s}_{\lambda }^{(-)}}^{2}/{n}_{\lambda }^{(-)}}$$
(21)

When summing up scattering rates in different categories, based on the rule of error propagation, the confidence interval of \({\hat{\tau }}_{\lambda ,3{{{\rm{ph}}}}}^{-1}\) is given by:

$${\hat{\tau }}_{\lambda ,3{{{\rm{ph}}}}}^{-1}\pm {\left({\left(C{I}_{\lambda ,3{{{\rm{ph}}}}}^{(+)}\right)}^{2}+{\left(C{I}_{\lambda ,3{{{\rm{ph}}}}}^{(-)}\right)}^{2}\right)}^{1/2}$$
(22)

The confidence interval of \({\hat{\tau }}_{\lambda ,4{{{\rm{ph}}}}}^{-1}\) can be derived in a similar manner:

$${\hat{\tau }}_{\lambda ,4{{{\rm{ph}}}}}^{-1}\pm {\left({\left(C{I}_{\lambda ,4{{{\rm{ph}}}}}^{(++)}\right)}^{2}+{\left(C{I}_{\lambda ,4{{{\rm{ph}}}}}^{(+-)}\right)}^{2}+{\left(C{I}_{\lambda ,4{{{\rm{ph}}}}}^{(--)}\right)}^{2}\right)}^{1/2}$$
(23)

where \(C{I}_{\lambda ,4{{{\rm{ph}}}}}^{(++)}\), \(C{I}_{\lambda ,4{{{\rm{ph}}}}}^{(+-)}\) and \(C{I}_{\lambda ,4{{{\rm{ph}}}}}^{(--)}\) are given by:

$$C{I}_{\lambda ,4{{{\rm{ph}}}}}^{(++)}={t}_{\alpha /2,{n}_{\lambda }^{(++)}}{N}_{\lambda }^{(++)}\frac{1}{2}\sqrt{{{s}_{\lambda }^{(++)}}^{2}/{n}_{\lambda }^{(++)}}$$
(24)
$$C{I}_{\lambda ,4{{{\rm{ph}}}}}^{(+-)}={t}_{\alpha /2,{n}_{\lambda }^{(+-)}}{N}_{\lambda }^{(+-)}\frac{1}{2}\sqrt{{{s}_{\lambda }^{(+-)}}^{2}/{n}_{\lambda }^{(+-)}}$$
(25)
$$C{I}_{\lambda ,4{{{\rm{ph}}}}}^{(--)}={t}_{\alpha /2,{n}_{\lambda }^{(--)}}{N}_{\lambda }^{(--)}\frac{1}{6}\sqrt{{{s}_{\lambda }^{(--)}}^{2}/{n}_{\lambda }^{(--)}}$$
(26)

Predicting lattice thermal conductivity

The total scattering rate of a phonon mode τλ is calculated based on spectral Matthiessen’s rule11:

$${\tau }_{\lambda }^{-1}=\left\{\begin{array}{rlr}{\tau }_{\lambda ,3{{{\rm{ph}}}}}^{-1}&&3{{{\rm{ph}}}}\,{{{\rm{scattering}}}}\\ {\tau }_{\lambda ,3{{{\rm{ph}}}}}^{-1}+{\tau }_{\lambda ,4{{{\rm{ph}}}}}^{-1}&&\,{{\mbox{3ph+4ph scattering}}}\,\end{array}\right.$$
(27)

The lattice thermal conductivity is calculated by:

$$\kappa =\frac{1}{V}\mathop{\sum}\limits_{\lambda }{c}_{\lambda }{v}_{\lambda }^{2}{\tau }_{\lambda },$$
(28)

where V is the unit cell volume, vλ and cλ are the group velocity and the specific heat of phonon mode λ, respectively.

Predicting thermal radiative properties

The complex dielectric function of polar dielectrics in the mid-IR range can be described by the Lorentz oscillator model57,58,59:

$$\varepsilon (\omega )={\varepsilon }_{\infty }\mathop{\prod}\limits_{m}\left(\frac{{\omega }_{{{{\rm{LO}}}},m}^{2}-{\omega }^{2}-i{\gamma }_{{{{\rm{LO}}}},m}\omega }{{\omega }_{{{{\rm{TO}}}},m}^{2}-{\omega }^{2}-i{\gamma }_{{{{\rm{TO}}}},m}\omega }\right)$$
(29)

where ε∞ is the dielectric constant at the high-frequency limit, ω is the photon frequency, ωLO,m and ωTO,m are frequencies of the zone-center IR active LO an TO phonon modes, respectively. γLO,m and γTO,m are the damping factors, which can be derived from the phonon scattering rate of zone-center IR active LO and TO phonon modes, respectively.

From the dielectric function, we can further derive many useful thermal radiative properties including the complex refractive index (m) and the normal reflectance of the air-materials interface (R).

$$m=\sqrt{\varepsilon }$$
(30)
$$R={\left| \frac{\sqrt{\varepsilon }-1}{\sqrt{\varepsilon }+1}\right| }^{2}$$
(31)