Fast evaluating phonon life time and thermal conductivity determined by Grüneisen parameter and phase space size of three-phonon scattering

Abstract

Efficiently and fast seeking specific lattices with targeted phonon thermal conductivity {\kappa }_L plays an important role in the thermal design and thermal management of materials. How to efficiently and accurately evaluate the phonon lifetime determined by anharmonicity becomes a critical bottleneck when high-throughput measuring {\kappa }_L. Here, we propose a method of fast evaluating three-phonon scattering induced lifetime based on the many-body theory of phonon gas. In the high temperature limit, the phonon scattering rate is simply determined by the product of only two anharmonic parameters: the square of Grüneisen parameter and the phase space size of three-phonon scattering, both of which can be quickly derived from the harmonic phonon properties. We demonstrate the effectiveness of the method in high-throughput evaluating the {\kappa }_L in first-principles calculation, which exhibits a good consistence with our collected experimental data. This method shows promising application potential in exploring material screening of the targeted {\kappa }_L, which by improving the ability of characterizing phonon anharmonicity will further enhance the performance of {\kappa }_L prediction.

Graphical abstract

Keywords

many-body theory / thermal conductivity / anharmonic properties / three-phonon scattering / first-principles calculation / high-throughput calculations

1 Introduction

Phonons are quasi-particles, commonly known as the quantized unites of lattice vibrations [1, 2]. These vibrations result from the superposition of acoustic waves, carrying energy and thus conducting heat, with varying frequencies and each frequency corresponds to a specific phonon energy. As a result, phonon propagation can be considered as the primary source of {\kappa }_L. Notably, {\kappa }_L serves as a fundamental physical property that characterizes a material or system, and can be utilized as a critical physical quantity for measuring its heat transfer capacity. Moreover, {\kappa }_L plays a crucial role in the control and management of thermal energy [3–5]. Materials with high or low {\kappa }_L demonstrate great diverse potential for applications in optoelectronics, thermoelectrics, and energy conversion and storage devices [6–9]. For instance, high {\kappa }_L materials can significantly enhance the heat dissipation of materials, leading to improved cooling performance [10]. On the other hand, low {\kappa }_L materials are commonly utilized in thermoelectric materials and devices, which can directly convert energy into useful electrical energy under specific temperature differences [11]. In summary, materials exhibiting high or low {\kappa }_L are of paramount importance in enhancing energy utilization efficiency and thermal management capabilities.

The phonon lifetime (\tau) plays a crucial role in determining {\kappa }_L. In harmonic crystals, where phonon−phonon interactions are absent, \tau is infinite, irrespective of any defects or grain boundaries. In contrast, non-harmonic waves that enable the interaction between phonons with different modes lead to phonon scattering and the subsequent decay of phonon energy. As a result, the inclusion of non-harmonic waves decreases \tau. This significant anharmonic property of phonons effectively enhances phonon scattering and reduces {\kappa }_L, as supported by numerous theoretical calculations and experimental studies [11–17].

Fast and accurately evaluating the \tau thereby {\kappa }_L is crucial for exploring the thermal functional materials theoretically. Based on the phonon Boltzmann transport equation (BTE) [18, 19], the {\kappa }_L could be accurately calculated with the third or higher order interatomic force constants (IFCs). This approach have undergone extensive experimental verification and demonstrate significant relevance [20–22]. However, it necessitates substantial computational resources for calculating \tau, which is not suitable to high-throughput calculate the {\kappa }_L. With reasonable assumptions and simplifications, there are some semi-empirical models that enable to rapidly predict the {\kappa }_L [23–30]. Based on Debye linear phonon dispersion approximation, the anharmonicity of is generally captured by the Grüneisen parameter (\gamma). The Slack model [23–25] and Debye−Callaway model [26, 27] has concise form and clear physical meaning, and is able to reasonably predict the {\kappa }_L. The Cahill [28] model and the Clarke model [29] are primarily utilized to measure the lower limit of {\kappa }_L in solid materials. Moreover, with the sine type phonon dispersion approximation, the PET model [30] could accurately predict the lattice thermal conductivities of crystalline materials covering seven crystal systems. However, the phonon BTE theory shows that the anharmonicity of phonon is determined by both frequency-dependent \gamma and phase space size of phonon scattering P_3, and the P_3 plays an important role in anharmonicity. There are also more works to emphasize the importance of the P_3 to determine the {\kappa }_L [31–35]. Therefore, further considering both the \gamma and the P_3 with phonon dispersion could achieve fast and accurate prediction of the {\kappa }_L.

In this work, we introduce an effective approach for rapidly evaluating the \tau induced by three-phonon scattering of the {\kappa }_L. With the high-temperature limit, the phonon scattering rate is determined by anharmonic parameter Grüneisen parameter \gamma and the phase space size of three-phonon scattering P_3. The results of high-throughput calculation with this method are consistent with the experimental {\kappa }_L. Our approach shows the potential in exploring material screening of targeted {\kappa }_L, as it improves the ability to characterize phonon anharmonicity and enhances the performance of {\kappa }_L prediction. Meanwhile, by enabling rapid and accurate prediction of the {\kappa }_L, this method has important applications for the development of new materials with desired {\kappa }_L properties.

2 Methods

2.1 Three-phonon scattering process

From the perspective of the phonon gas model, the {\kappa }_L depends on the summation of heat capacity (C), the group velocity (\upsilon) and the relaxation time (\tau), of every individual phonon mode. Commonly, the above physical quantities are related to different phonon frequencies (\omega), branches and wave vectors (to simplify notation, we use \lambda to comprise both a phonon branch index and a wave vector). Therefore, one can calculate the {\kappa }_L of materials in the \alpha and \beta directions through the following relation [2, 36–41]:

(1){\kappa }_{L,\alpha \beta }={\displaystyle \frac{1}{V}}{\displaystyle \sum _{\lambda ,\omega }{C}_{\lambda }(\omega ){\upsilon }_{\lambda ,\alpha }(\omega ){\upsilon }_{\lambda ,\beta }(\omega ){\tau }_{\lambda }(\omega ),}

where V is the volume of the unit cell, C_{\lambda }, {\upsilon }_{\lambda \alpha } and {\upsilon }_{\lambda \beta } indicate the phonon harmonic properties. And {\tau }_{\lambda } is the phonon lifetime in the relaxation time approximation (RTA), it is the reciprocal of phonon scattering rate with considering the phonon−phonon interactions, indicating anharmonic properties. However, it conventionally consumes too many computational resources on evaluating anharmonic properties, i.e., the phonon lifetime {\tau }_{\lambda }, based on the strictly iterative method [18, 19], which becomes the bottleneck of the prediction of {\kappa }_L in high-throughput calculations. Therefore, we consider simplifying the expression of {\tau }_{\lambda } in calculating the phonon scattering process based on the perturbation theory [36, 42].

When considering the phonon−phonon interaction, different numbers of phonons are involved in the phonon scattering processes. In general, the intrinsic phonon scattering process refers to three-phonon, four-phonon, and higher-order phonon scattering process. Here, in order to balance the computational efficiency and accuracy, we mainly consider the three-phonon scattering process, as shown in Fig.1. Considering the conditions of energy ({\omega }_{\lambda }\pm {\omega }_{{\lambda }^{\prime }}={\omega }_{{\lambda }^{″}}) and quasi-momentum conservation (\mathit{q}+{\mathit{q}}^{\mathbf{\prime }}+{\mathit{q}}^{\mathbf{″}}=\mathit{G}) in the simplest case, three phonons interact and one of which decays to create the other two, and vice versa. The cases of \mathit{G}=0 and \mathit{G}\ne 0 are known as normal and Umklapp processes that correspond to the contributions of harmonic and anharmonic interactions between phonons. Hence, we can then write the total Hamiltonian (\widehat{H_L}) of the lattice according to the perturbation theory [39, 42]:

(2)\widehat{H_L}=\widehat{H_0}+{\displaystyle \sum _{i=3,4,...}{\hat{H}}_i,}

where \widehat{H_0} and \widehat{H_i} are the harmonic and anharmonic term in the phonon scattering process, respectively. Here, we mainly consider the three-phonon scattering process and i is equal to 3. In the language of second quantization, the first term in Eq. (2) can be expressed as \widehat{H_0}=\sum _{\lambda }\hslash {\omega }_{\lambda } \left(a_{\lambda }^{†}{a}_{\lambda }+\frac{1}{2}\right), and the second term in Eq. (2) can be expressed as [39, 42]

(3)\widehat{H_3}={\displaystyle \sum _{\lambda {\lambda }^{\prime }{\lambda }^{″}}{H}_{\lambda {\lambda }^{\prime }{\lambda }^{″}}^{(3)}(a_{-\lambda }^{†}+a_{\lambda })(a_{-{\lambda }^{\prime }}^{†}+a_{{\lambda }^{\prime }})(a_{-{\lambda }^{″}}^{†}+a_{{\lambda }^{″}}),}

(4)H_{\lambda {\lambda }^{\prime }{\lambda }^{″}}^{(3)}={\displaystyle \frac{1}{3!}}{\displaystyle \sum _{\lambda {\lambda }^{\prime }{\lambda }^{″}}\sqrt{{\displaystyle \frac{{\hslash }^3}{2^{3}N}}}{\displaystyle \frac{V(\lambda ,{\lambda }^{\prime },{\lambda }^{″})}{\sqrt{{\omega }_{\lambda }{\omega }_{{\lambda }^{\prime }}{\omega }_{{\lambda }^{″}}}}}\mathrm{\Delta }({\mathit{q}}^{\mathbf{\prime }}+{\mathit{q}}^{\mathbf{″}}+\mathit{q},\mathit{G}),}

where a_{\lambda }^{†} and a_{\lambda } are the creation and annihilation operators of phonons, which satisfy a_{\lambda }^{†}|n_{\lambda }⟩=\sqrt{n_{\lambda }+1}|n_{\lambda }+1⟩ and a_{\lambda }|n_{\lambda }⟩=\sqrt{n_{\lambda }}|n_{\lambda } -1⟩, respectively. Expansion of ( a_{-\lambda }^{†}+a_{\lambda } )( a_{-{\lambda }^{\prime }}^{†}+a_{{\lambda }^{\prime }})(a_{{\lambda }^{″}}^{†}+a_{{\lambda }^{″}}) in Eq. (3) represents different processes of creation and annihilation of the three phonons, which are mainly shown in Fig.1. n_{\lambda } is the Bose−Einstein distribution n_{\lambda }=[{\mathrm{e}}^{\hslash {\omega }_{\lambda }/(k_{B}T)}- 1{]}^{-1}. \hslash, k_B, T and N are the reduced Planck constant, the Boltzmann constant, the temperature and the total number of \mathit{q} points. Besides, the Kronecker delta function \mathrm{\Delta } guarantees the conservation of momentum in the scattering process. V(\lambda ,{\lambda }^{\prime },{\lambda }^{″}) is the phonon scattering matrix element due to anharmonic interactions, whose specific expression is as follows [18, 39, 42]:

(5)\begin{array}{rl}V(\lambda ,{\lambda }^{\prime },{\lambda }^{″})=& \sum _{k{k}^{\prime }{k}^{″}}\sum _{l{l}^{\prime }{l}^{″}}\sum _{\alpha \beta \eta }{\mathrm{\Phi }}_{lk,l^{\prime }{k}^{\prime },l^{″}{k}^{″}}^{\alpha \beta \eta }(\lambda ,{\lambda }^{\prime },{\lambda }^{″})\\ & \cdot \frac{{\mathit{e}}_{\alpha }(k,\lambda ){\mathit{e}}_{\beta }(k^{\prime },{\lambda }^{\prime }){\mathit{e}}_{\eta }(k^{″},{\lambda }^{″})}{\sqrt{m_k{m}_{k^{\prime }}{m}_{k^{″}}}}\\ & \cdot {\mathrm{e}}^{\left[\mathrm{i}\left(\mathit{q}\cdot \mathit{R}(l)+{\mathit{q}}^{\mathbf{\prime }}\cdot \mathit{R}(l^{\prime })+{\mathit{q}}^{\mathbf{″}}\cdot \mathit{R}(l^{″})\right)\right]},\end{array}

where l, k, \mathit{R}, and m refer to the ordinal number of the primitive cell, the specific atomic ordinal number in the primitive cell, the lattice vector, and the corresponding mass of the atom, respectively. \alpha, \beta and \eta are the atomic coordinates, {\mathit{e}}_{\alpha }(k, \mathit{q}\lambda ) is phonon eigenvector. And {\mathrm{\Phi }}_{lk,l^{\prime }{k}^{\prime },l^{″}{k}^{″}}^{\alpha \beta \eta }(\lambda ,{\lambda }^{\prime },{\lambda }^{″}) is the third-order interatomic force constant. In summary, by solving the third-order force constant matrix {\mathrm{\Phi }}_{lk,l^{\prime }{k}^{\prime },l^{″}{k}^{″}}^{\alpha \beta \eta }(\lambda ,{\lambda }^{\prime },{\lambda }^{″}) in the three-phonon scattering process and combining with Eqs. (2)−(5), we can obtain the V(\lambda ,{\lambda }^{\prime },{\lambda }^{″}) and the anharmonic phonon−phonon interaction \widehat{H_3}.

In addition, according to the \widehat{H_3}, the total phonon scattering rate \mathrm{\Gamma } in the three-phonon scattering process can be calculated from the Fermi’s golden rule [42, 43]:

(6)\mathrm{\Gamma }={\displaystyle \frac{2\pi }{\hslash }}{\left|⟨j|\widehat{H_3}|i⟩\right|}^2\delta ({\omega }_j-{\omega }_i),

where |i⟩ and |j⟩ are the initial and final states of the three phonon scattering states. {\omega }_i and {\omega }_j are the phonon energy corresponding to the initial and final states.

2.2 Deriving phonon lifetime {\tau }_{\lambda } and thermal conductivity {\kappa }_L

Here, we primarily consider the intrinsic phonon scattering rate in Fig.1, according to the single mode RTA, which is the reciprocal of {\tau }_{\lambda } in the phonon scattering process. Thus, it is given by [18]

(7){\displaystyle \frac{1}{{\tau }_{\lambda }}}={\displaystyle \frac{1}{N}}{\displaystyle \sum _{{\lambda }^{\prime }{\lambda }^{″}}\left({\mathrm{\Gamma }}_{\lambda {\lambda }^{\prime }{\lambda }^{″}}^{+}+{\displaystyle \frac{1}{2}}{\mathrm{\Gamma }}_{\lambda {\lambda }^{\prime }{\lambda }^{″}}^{-}\right).}

Based on the conservation of phonon energy and quasi-momentum, the first two terms on the right of Eq. (7) represent the phonon decay process (two phonons emit one phonon) and creation process (one phonon decomposes into two phonons). Here, the coefficient of 1/2 indicates that the decomposed two phonons are identical and cannot be distinguished, so as to avoid calculation duplication. Then, on the basis of the Fermi’s golden rule from Eq. (6) and Fig.1, the scattering rates of {\mathrm{\Gamma }}^{+} and {\mathrm{\Gamma }}^{-} can be further written as [18]

(8){\mathrm{\Gamma }}_{\lambda {\lambda }^{\prime }{\lambda }^{″}}^{+}={\displaystyle \frac{\hslash \pi }{4}}{\displaystyle \frac{n_{{\lambda }^{\prime }}-n_{{\lambda }^{″}}}{{\omega }_{\lambda }{\omega }_{{\lambda }^{\prime }}{\omega }_{{\lambda }^{\prime }}}}{\left|V_{\lambda {\lambda }^{\prime }{\lambda }^{″}}^{+}\right|}^2\delta ({\omega }_{\lambda } +{\omega }_{{\lambda }^{\prime }}- {\omega }_{{\lambda }^{″}}) ,

(9){\mathrm{\Gamma }}_{\lambda {\lambda }^{\prime }{\lambda }^{″}}^{-}={\displaystyle \frac{\hslash \pi }{4}}{\displaystyle \frac{n_{{\lambda }^{\prime }}+n_{{\lambda }^{″}}+1}{{\omega }_{\lambda }{\omega }_{{\lambda }^{\prime }}{\omega }_{{\lambda }^{″}}}}{\left|V_{\lambda {\lambda }^{\prime }{\lambda }^{″}}^{-}\right|}^2\delta ({\omega }_{\lambda }-{\omega }_{{\lambda }^{\prime }}-{\omega }_{{\lambda }^{″}}).

|V^{\pm }| represent the scattering matrix elements corresponding to the phonon creation and decay processes. Then, we substitute Eqs. (8) and (9) into Eq. (7), and further rewrite 1/{\tau }_{\lambda } as

(10)\begin{array}{rl}\frac{1}{{\tau }_{\lambda }}=& \pi \sum _{{\lambda }^{\prime }{\lambda }^{″}}{\left|V(\lambda ,{\lambda }^{\prime },{\lambda }^{″})\right|}^2[2(n_{{\lambda }^{\prime }}-n_{{\lambda }^{″}})\delta ({\omega }_{\lambda }+{\omega }_{{\lambda }^{\prime }}-{\omega }_{{\lambda }^{″}})\\ & +(n_{{\lambda }^{\prime }}+n_{{\lambda }^{″}}+1)\delta ({\omega }_{\lambda }-{\omega }_{{\lambda }^{\prime }}-{\omega }_{{\lambda }^{″}})].\end{array}

According to Eq. (5), we can further substitute it into Eq. (10) and obtain the following form:

(11)\begin{array}{rl}\frac{1}{{\tau }_{\lambda }}=& \frac{\pi \hslash }{8N}\sum _{{\lambda }^{\prime }{\lambda }^{″}}\frac{{\left|\mathrm{\Phi }(\lambda ,{\lambda }^{\prime },{\lambda }^{″})\right|}^2}{{\omega }_{\lambda }{\omega }_{{\lambda }^{\prime }}{\omega }_{{\lambda }^{″}}}\mathrm{\Delta }({\mathit{q}}^{\mathbf{\prime }}+{\mathit{q}}^{\mathbf{″}}+\mathit{q})\\ & \times [2(n_{{\lambda }^{\prime }}-n_{{\lambda }^{″}})\delta ({\omega }_{\lambda }+{\omega }_{{\lambda }^{\prime }}-{\omega }_{{\lambda }^{″}})\\ & +(n_{{\lambda }^{\prime }}+n_{{\lambda }^{″}}+1)\delta ({\omega }_{\lambda }-{\omega }_{{\lambda }^{\prime }}-{\omega }_{{\lambda }^{″}})].\end{array}

In order to simplify the calculation of 1/{\tau }_{\lambda } in Eq. (11), we firstly consider the approximate expression of Bose−Einstein distribution at the limitation of high temperature (\hslash {\omega }_{\lambda }\ll k_{B}T), based on the conservation of phonon energy. Therefore, we can finally simplify the last item on the right of Eq. (11), and the specific processes are

(12)\begin{array}{rl}{n}_{{\lambda }^{\prime }}-n_{{\lambda }^{″}}=& \frac{1}{{\mathrm{e}}^{\hslash {\omega }_{{\lambda }^{\prime }}/(k_{B}T)}-1}-\frac{1}{{\mathrm{e}}^{\hslash {\omega }_{{\lambda }^{″}}/(k_{B}T)}-1}\\ \approx & \frac{k_{B}T}{\hslash }\frac{{\omega }_{{\lambda }^{″}}-{\omega }_{{\lambda }^{\prime }}}{{\omega }_{{\lambda }^{\prime }}{\omega }_{{\lambda }^{″}}}=\frac{k_{B}T}{\hslash }\frac{{\omega }_{\lambda }}{{\omega }_{{\lambda }^{\prime }}{\omega }_{{\lambda }^{″}}},\end{array}

(13)\begin{array}{rl}{n}_{{\lambda }^{\prime }}+n_{{\lambda }^{″}}+1=& \frac{1}{{\mathrm{e}}^{\hslash {\omega }_{{\lambda }^{\prime }}/(k_{B}T)}-1}+\frac{1}{{\mathrm{e}}^{\hslash {\omega }_{{\lambda }^{″}}/(k_{B}T)}-1}+1\\ \approx & \frac{k_{B}T}{\hslash }\frac{{\omega }_{{\lambda }^{\prime }}+{\omega }_{{\lambda }^{″}}}{{\omega }_{{\lambda }^{\prime }}{\omega }_{{\lambda }^{″}}}=\frac{k_{B}T}{\hslash }\frac{{\omega }_{\lambda }}{{\omega }_{{\lambda }^{\prime }}{\omega }_{{\lambda }^{″}}}.\end{array}

Meanwhile, in the view of a simplified Grüneisen model based on the phonon coupling constant proposed by Klemens, the third-order interatomic force constant \mathrm{\Phi }( \lambda ,{\lambda }^{\prime } ,{\lambda }^{″}) in Eq. (11) can also be expressed generally as [44]

(14)\mathrm{\Phi }( \lambda ,{\lambda }^{\prime } ,{\lambda }^{″}) ={\displaystyle \frac{2\mathrm{i}\gamma (\omega ,\lambda )}{(3m{)}^{1/2}\upsilon }}{\omega }_{\lambda }{\omega }_{{\lambda }^{\prime }}{\omega }_{{\lambda }^{″}},

where \gamma (\omega ,\lambda ), \upsilon and m are the Grüneisen parameter, the sound velocity and the atomic mass, respectively. Importantly, \gamma (\omega ,\lambda ) usually characterizes the intensity of anharmonic properties in the process of three phonon scattering, and further directly determines the 1/ {\tau }_{\lambda }. Moreover, the \gamma is considered the dependence on \omega and \lambda, and it is given by [45]

(15)\gamma (\omega ,\lambda )=-{\displaystyle \frac{V}{{\omega }_{\lambda }}}{\displaystyle \frac{\mathrm{\partial }{\omega }_{\lambda }}{\mathrm{\partial }V}}.

The equation is based on the partial differentiation of the phonon frequency relative to the volume of the material. There are mostly spending too many computational resources and time on obtaining the \mathrm{\Phi } (\lambda , {\lambda }^{\prime },{\lambda }^{″}) for iteratively solving phonon BTE. Hence, Eq. (14) establishes the relationship between the \mathrm{\Phi }( \lambda ,{\lambda }^{\prime } ,{\lambda }^{″}) and the \gamma (\omega ,\lambda ), which can be applied to the estimation of the {\kappa }_L of materials using high-throughput calculation.

Then, we substitute Eq. (12), Eq. (13) and Eq. (14) into Eq. (11) and further obtain the following expression:

(16)\begin{array}{rl}\frac{1}{{\tau }_{\lambda }}=& \frac{\pi k_{B}T{\omega }_{\lambda }^2{\gamma }^2(\omega ,\lambda )}{6m{\upsilon }^2}\cdot \sum _{{\lambda }^{\prime }{\lambda }^{″}}\mathrm{\Delta }({\mathit{q}}^{\mathbf{\prime }}+{\mathit{q}}^{\mathbf{″}}+\mathit{q})\\ & \cdot \left[2\delta ({\omega }_{\lambda }+{\omega }_{{\lambda }^{\prime }}-{\omega }_{{\lambda }^{″}})+\delta ({\omega }_{\lambda }-{\omega }_{{\lambda }^{\prime }}-{\omega }_{{\lambda }^{″}})\right].\end{array}

Ulteriorly, in the process of three-phonon scattering, in addition to the fact that \gamma (\omega ) affects the anharmonicity of phonons and thus determines the phonon scattering rate 1/{\tau }_{\lambda }, the number of phonons allowed to pass through the three-phonon scattering channel is another important physical quantity that determines the scattering rate. Normally, we define the number of phonon involving into the three-phonon scattering channels as the size of the total phonon scattering phase space, P_3. The larger the P_3 is, the much phonon involved in three phonon scattering, the greater the scattering probability, the shorter the corresponding phonon relaxation time {\tau }_{\lambda }, and the lower the final {\kappa }_L. Therefore, according to the definition of the P_3 [46, 47], the second term on the right of equation in Eq. (11) can be replaced by P_3(\omega ,\lambda ), which is also related to \omega and \lambda in the phonon scattering space. Ultimately, we obtain the expression of phonon scattering rate 1/{\tau }_{\lambda }:

(17){\displaystyle \frac{1}{{\tau }_{\lambda }(\omega )}}={\displaystyle \frac{\pi k_{B}T{\omega }_{\lambda }^2}{6m{\upsilon }^2}}\cdot {\gamma }^2(\omega ,\lambda )\cdot P_3(\omega ,\lambda ).

At the same time, we derive the general expression for {\tau }_{\lambda }(\omega ) as follows:

(18){\tau }_{\lambda }(\omega )={\displaystyle \frac{6m{\upsilon }^2}{\pi k_{B}T{\omega }_{\lambda }^2}}\cdot {\gamma }^{-2}(\omega ,\lambda ) \cdot P_3^{-1}(\omega ,\lambda ) .

Combining Eq. (1) with Eq. (18), we then obtain the final expression of {\kappa }_L for the sum of all phonon modes and frequencies:

(19){\kappa }_L={\displaystyle \frac{1}{V}}{\displaystyle \frac{6m{\upsilon }^2}{\pi k_{B}T}}\cdot {\displaystyle \sum _{\lambda ,\omega }{C}_{\lambda }(\omega ){\upsilon }_{\lambda }^2(\omega ){\gamma }_{\lambda }^{-2}(\omega )P_3^{-1}(\omega ,\lambda ){\omega }_{\lambda }^{-2}.}

This formula directly relates the phonon harmonic (C_{\lambda }(\omega ), {\upsilon }_{\lambda }(\omega )) and anharmonic ({\gamma }_{\lambda }(\omega ), P_3(\omega ,\lambda )) physical quantities that affect the {\kappa }_L of materials during phonon thermal transport. By theoretically predicting the contribution of all phonon modes of these physical quantities at different phonon frequency \omega based on the first-principles calculations, {\tau }_{\lambda }(\omega ) determined by phonon anharmonicity with the \gamma and P_3 together and then the {\kappa }_L can be quickly obtained.

Meanwhile, not only does this method through Eq. (19) consider the three phonon scattering process from the phonon transport theory and have basic physical support, but it also avoids the waste of computing resources by approximately solving the third-order force constant, and then uses the two physical quantities {\gamma }_{\lambda }(\omega ) and P_3(\omega ,\lambda ) to measure the anharmonicity of phonons, which reasonably characterizes the {\kappa }_L of the materials. Therefore, we can further use this theoretical model to estimate the {\kappa }_L of a series of materials with high-throughput calculations, and finally explore more potential novel materials with high or low target {\kappa }_L.

3 First-principles results and discussion

Based on Eq. (19), we calculate all of the physical quantities through first-principles calculations. Importantly, on basis of the Vienna Ab-initio Simulation Package (VASP) [48, 49], the pseudo-potential density functional theory (DFT) of plane waves is calculated. Our geometric optimization and self-consistent energy calculations are implemented in VASP using the projected augmented wave (PAW) method [50, 51]. For the exchange correlation energy, Perdew−Burke−Ernzerhof (PBE) is used in the calculation of the generalized gradient approximation (GGA) [52]. To guarantee the complete convergence of the optimized structure, we use the 520 eV with a kinetic energy cutoff. A 9 \times 9 \times 9 Monkhorst−Pack K-point grid was used for Brillouin zone sampling. Meanwhile, we choose {10}^{-5} eV and 0.01 eV/Å as the energy convergence criterion and the maximum Herman−Feynman coefficient force convergence criterion, respectively. Furthermore, we calculate the phonon dispersions, the second-order IFCs, C_{\lambda }(\omega ), {\upsilon }_{\lambda }(\omega ) and \gamma (\omega ) of materials based on the the density functional perturbation theory (DFPT) [53], as implemented in the Phonopy package [54].

In the previous exploration, Yan et al. [30] directly estimated the {\kappa }_L through the phonon-elasticity-thermal (PET) model, by establishing the relationship between the elastic properties of the materials and the {\kappa }_L, and taking into account the contributions of the acoustic and optical phonons. As a result, the value of the {\kappa }_L is determined by the harmonic property dominated by the phonon group velocity and the anharmonic property dominated by the Grüneisen parameter. Among 226 materials covering all crystal systems, the predicted values of the model are in good agreement with the experimental values reported by the materials except for the poor prediction of anharmonicity in the trigonal crystal system [30]. Hence, we can further explore the novel way to improve the ability to predict the phonon anharmonicity of materials in this work, which is also crucial for us to excavate the target {\kappa }_L of functional materials based on the high-throughput calculations.

Commonly, in the process of phonon thermal transport of materials, the interactions between phonons directly determine the value of {\kappa }_L [22, 55]. And we can accurately describe the phonon scattering rate by calculating the harmonic and anharmonic properties. According to Eq. (19) derived by this work, \gamma and P_3 are usually used to characterize the intensity of phonon anharmonic properties, and the larger \gamma and P_3 indicate that the material has a stronger anharmonicity, which will eventually lead to the lower {\kappa }_L in theory [40, 41, 56]. For this purpose, in order to enhance the ability to evaluate the anharmonicity of material, we add the physical quantity \gamma (\omega ) and P_3(\omega ) to jointly measure the anharmonicity of phonons based on the phonon transport theory and the RTA. At the same time, these physical quantities take into account the contributions of all phonon modes to the {\kappa }_L at different frequencies, which can further improve the ability to characterize the anharmonicity. We can quickly evaluate the third-order force constants \mathrm{\Phi } (\lambda , {\lambda }^{{}^{\prime }},{\lambda }^{{}^{″}}) by using the approximate relationship between the \gamma (\omega ) in Eq. (14) to avoid spending massive computational resources. Ultimately, the accuracy of Eq. (19) will be verified through high-throughput first-principles calculation in view of the comparison with some experimental data.

Firstly, we have collected 63 structures with experimental values of the {\kappa }_L, which covering trigonal, hexagonal and cubic crystal systems based on the work reported by Yan et al. [30]. The structural numbers of these three systems are 16, 21 and 26, respectively. Then, we show the {\kappa }_L predicted by the theoretical model in this work are compared with the experimental results at room temperature in Fig.2. Among them, Pearson correlation coefficient between the predicted and the experimental values is 0.95, and the above crystal systems are corresponded to 0.87, 0.94 and 0.93. The results indicate the evaluated accuracy of this model by the good agreement between the assessment of our model with the experimental values, preliminarily proving the effectiveness of our model’s prediction.

Further, in order to demonstrate the predicted ability of the {\kappa }_L for these crystal systems, we also compare these values with the Slack model, as shown in Fig.3. The Slack model has been widely applied for the fast evaluation of the {\kappa }_L with minimal time and computational resources, showing the potential capability of high-throughput screening of {\kappa }_L. The expression of the {\kappa }_{\mathrm{S}\mathrm{l}\mathrm{a}\mathrm{c}\mathrm{k}} is as follows [22, 41, 57]:

(20){\kappa }_{\mathrm{S}\mathrm{l}\mathrm{a}\mathrm{c}\mathrm{k}}=A{\displaystyle \frac{\overline{M}{\mathrm{\Theta }}_D^3\delta }{{\gamma }^2{n}^{2/3}T}},

where \overline{M}, {\delta }^3, n, {\mathrm{\Theta }}_D, \gamma and T are respectively the average atomic mass, the atomic average volume in the primitive unit cell, the number of atoms in the primitive unit cell, the Debye temperature, the Grüneisen constant and the temperature. Moreover, A is defined as A=2.43\times {10}^{-8}\times (1-0.514 /\gamma +0.228/{\gamma }^2 {)}^{-1} [58]. {\mathrm{\Theta }}_D can be obtained by the average sound velocity {\upsilon }_m [59, 60]: {\mathrm{\Theta }}_D=\frac{h}{k_B}{\left[\frac{3n}{4\pi }(\frac{N_A\rho }{M})\right]}^{\frac{1}{3}}{\upsilon }_m, where h, k_B, M, n, N_A and \rho are the Plank’s constant, the Boltzmann’s constant, the weight of molecule, the number of atoms in the primitive cell, the Avogadro’s constant and the density, respectively. Meanwhile, the average sound velocity {\upsilon }_m stems from the transverse acoustic velocity {\upsilon }_s and the longitudinal acoustic velocity {\upsilon }_l: {\upsilon }_m={\left[\frac{1}{3}(\frac{2}{{\upsilon }_s^3}+\frac{1}{{\upsilon }_l^3})\right]}^{-\frac{1}{3}}, and they are based on the mechanical properties, namely bulk modulus B and shear modulus G. {\upsilon }_s and {\upsilon }_l are respectively defined as {\upsilon }_s=\sqrt{\frac{G}{\rho }} and {\upsilon }_l= \sqrt{\frac{B+\frac{4G}{3}}{\rho }} [41, 59]. Lastly, \gamma is calculated through the Poisson’s ratio \nu based on the Slack model: \gamma = \frac{3(1+\nu )}{2(2-3\nu )} [41, 58, 59].

The predicted values of {\kappa }_L by the Slack model and this work are compared with the corresponding experimental values for trigonal [Fig.3(a)], hexagonal [Fig.3(b)] and cubic [Fig.3(c)] crystal systems at 300 K, respectively. The predicted values of the Slack model are relatively high, mainly due to only considering the contributions of acoustic phonons. The results show that our model further improves the prediction accuracy of the {\kappa }_L compared with the Slack model. And the intuitive data show that the Pearson correlation coefficients of our model and experimental values are improved over those of the Slack model, especially in the prediction of trigonal system. The Pearson correlation coefficient between our model and the experimental value is increased to 0.87, which is significantly improved compared with 0.75 of the Slack model.

In order to analyze the causes for the gap in the {\kappa }_L prediction, we further show the difference of \gamma with phonon anharmonicity in Fig.4. Here, we simultaneously compare the predicted mean values of \gamma in our model and the relationship between the predicted values of the Slack model and the reference values provided. A total of 33 values for reference are collected in the three crystal systems, and the calculated results show that the Pearson correlation coefficient between this work and the experimental value is 0.72, which is much higher than 0.49 of the Slack model. It indicates that the work has high prediction accuracy for \gamma, and also improves the prediction ability of \tau and {\kappa }_L dominated by anharmonic properties. At the same time, the introduction of P_3 in our model can also characterize the evaluation of anharmonicity. In a word, enhancing the capability to evaluate the phonon anharmonicity of materials would facilitate precise estimation of their {\kappa }_L, thereby offering robust support for the exploration of materials with the target {\kappa }_L that hold significant practical value.

4 Conclusions

In summary, we have developed a theoretical model to fast and accurately evaluate the phonon lifetime \tau induced by the three-phonon scattering based on the many-body theory of phonon gas. In the process of three-phonon scattering, the phonon scattering rate is determined by the product of the square of Grüneisen parameter {\gamma }^2 and the size of scattering phase space P_3. Moreover, we perform the high-throughput calculations to evaluate the {\kappa }_L covering trigonal, hexagonal and cubic crystal systems. The results show that further including the phase space size of phonon in our model could achieve the more accurate prediction of lattice thermal conductivity than the Slack model. Furthermore, we have demonstrated that the \gamma in our model is more accurate than that in the Slack model, and the existence of P_3 further characterizes the evaluation of phonon anharmonicity. This study provides a rapid and accurate method for evaluating {\kappa }_L and can aid in the exploration of targeted {\kappa }_L materials with broad applications in the future.

References

[1] M. Born,K. Huang, Dynamical theory of crystal lattices, International Series of Monographs on Physics, Oxford University Press, Oxford, 1954

[2] G. P. Srivastava, The Physics of Phonons, Taylor & Francis, New York, 1990

[3] N. Li, J. Ren, L. Wang, G. Zhang, P. Hänggi, and B. Li, Phononics: Manipulating heat flow with electronic analogs and beyond, Rev. Mod. Phys. 84(3), 1045 (2012)

[4] Y. Li, W. Li, T. Han, X. Zheng, J. Li, B. Li, S. Fan, and C. W. Qiu, Transforming heat transfer with thermal metamaterials and devices, Nat. Rev. Mater. 6(6), 488 (2021)

[5] J. Chen, J. He, D. Pan, X. Wang, N. Yang, J. Zhu, S. A. Yang, and G. Zhang, Emerging theory and phenomena in thermal conduction: A selective review, Sci. China Phys. Mech. 65, 117002 (2022)

[6] A. Polman, M. Knight, E. C. Garnett, B. Ehrler, and W. C. Sinke, Photovoltaic materials: Present efficiencies and future challenges, Science 352(6283), aad4424 (2016)

[7] R. Venkatasubramanian, E. Siivola, T. Colpitts, and B. O’quinn, Thin-film thermoelectric devices with high room-temperature figures of merit, Nature 413(6856), 597 (2001)

[8] Y. Ammar, S. Joyce, R. Norman, Y. Wang, and A. P. Roskilly, Low grade thermal energy sources and uses from the process industry in the UK, Appl. Energy 89(1), 3 (2012)

[9] C. Liu, C. Wu, T. Song, Y. Zhao, J. Yang, P. Lu, G. Zhang, and Y. Chen, An efficient strategy for searching high lattice thermal conductivity materials, ACS Appl. Energy Mater. 5(12), 15356 (2022)

[10] M. M. Waldrop, The chips are down for Moore’s law, Nature 530(7589), 144 (2016)

[11] L. D. Zhao, S. H. Lo, Y. Zhang, H. Sun, G. Tan, C. Uher, C. Wolverton, V. P. Dravid, and M. G. Kanatzidis, Ultralow thermal conductivity and high thermoelectric figure of merit in SnSe crystals, Nature 508(7496), 373 (2014)

[12] C. W. Li, J. Hong, A. F. May, D. Bansal, S. Chi, T. Hong, G. Ehlers, and O. Delaire, Orbitally driven giant phonon anharmonicity in SnSe, Nat. Phys. 11(12), 1063 (2015)

[13] J. He, M. Amsler, Y. Xia, S. S. Naghavi, V. I. Hegde, S. Hao, S. Goedecker, V. Ozoliņš, and C. Wolverton, Ultralow thermal conductivity in full Heusler semiconductors, Phys. Rev. Lett. 117(4), 046602 (2016)

[14] C. Chang and L. D. Zhao, Anharmoncity and low thermal conductivity in thermoelectrics, Mater. Today Phys. 4, 50 (2018)

[15] H. Xie,J. Yan,X. Gu,H. Bao, A scattering rate model for accelerated evaluation of lattice thermal conductivity bypassing anharmonic force constants, J. Appl. Phys. 125(20), 205104 (2019)

[16] J. Ding, T. Lanigan-Atkins, M. Calderón-Cueva, A. Banerjee, D. L. Abernathy, A. Said, A. Zevalkink, and O. Delaire, Soft anharmonic phonons and ultralow thermal conductivity in Mg₃(Sb,Bi)₂ thermoelectrics, Sci. Adv. 7(21), eabg1449 (2021)

[17] J. Zhang, H. Zhang, and G. Zhang, Nanophononic metamaterials induced proximity effect in heat flux regulation, Front. Phys. 19(2), 23204 (2024)

[18] W. Li, J. Carrete, N. A. Katcho, and N. Mingo, ShengBTE: A solver of the Boltzmann transport equation for phonons, Comput. Phys. Commun. 185(6), 1747 (2014)

[19] Z. Han, X. Yang, W. Li, T. Feng, and X. Ruan, Fourphonon: An extension module to ShengBTE for computing four-phonon scattering rates and thermal conductivity, Comput. Phys. Commun. 270, 108179 (2022)

[20] J. S. Kang, M. Li, H. Wu, H. Nguyen, and Y. Hu, Experimental observation of high thermal conductivity in boron arsenide, Science 361(6402), 575 (2018)

[21] S. Mukhopadhyay, D. S. Parker, B. C. Sales, A. A. Puretzky, M. A. McGuire, and L. Lindsay, Two-channel model for ultralow thermal conductivity of crystalline Tl₃VSe₄, Science 360(6396), 1455 (2018)

[22] X. Qian, J. Zhou, and G. Chen, Phonon-engineered extreme thermal conductivity materials, Nat. Mater. 20(9), 1188 (2021)

[23] G. A. Slack, Nonmetallic crystals with high thermal conductivity, J. Phys. Chem. Solids 34(2), 321 (1973)

[24] C. Toher, J. J. Plata, O. Levy, M. De Jong, M. Asta, M. B. Nardelli, S. Curtarolo, and H igh-throughput computational screening of thermal conductivity, Debye temperature, and Grüneisen parameter using a quasiharmonic Debye model, Phys. Rev. B 90(17), 174107 (2014)

[25] P. Nath, J. J. Plata, D. Usanmaz, C. Toher, M. Fornari, M. Buongiorno Nardelli, and S. Curtarolo, High throughput combinatorial method for fast and robust prediction of lattice thermal conductivity, Scr. Mater. 129, 88 (2017)

[26] J. Callaway, Model for lattice thermal conductivity at low temperatures, Phys. Rev. 113(4), 1046 (1959)

[27] J. Yan, P. Gorai, B. Ortiz, S. Miller, S. A. Barnett, T. Mason, V. Stevanović, and E. S. Toberer, Material descriptors for predicting thermoelectric performance, Energy Environ. Sci. 8(3), 983 (2015)

[28] D. G. Cahill, S. K. Watson, and R. O. Pohl, Lower limit to the thermal conductivity of disordered crystals, Phys. Rev. B 46(10), 6131 (1992)

[29] D. R. Clarke, Materials selection guidelines for low thermal conductivity thermal barrier coatings, Surf. Coat. Tech. 163–164, 67 (2003)

[30] S. Yan, Y. Wang, F. Tao, and J. Ren, High-throughput estimation of phonon thermal conductivity from first-principles calculations of elasticity, J. Phys. Chem. A 126(46), 8771 (2022)

[31] T. Pandey, C. A. Polanco, L. Lindsay, and D. S. Parker, Lattice thermal transport in La₃Cu₃X₄ compounds (X = P, As, Sb, Bi): Interplay of anharmonicity and scattering phase space, Phys. Rev. B 95(22), 224306 (2017)

[32] M. Raya-Moreno, R. Rurali, and X. Cartoixà, Thermal conductivity for III−V and II−VI semiconductor wurtzite and zincblende polytypes: The role of anharmonicity and phase space, Phys. Rev. Mater. 3(8), 084607 (2019)

[33] L. Zhu, W. Li, and F. Ding, Giant thermal conductivity in diamane and the influence of horizontal reflection symmetry on phonon scattering, Nanoscale 11(10), 4248 (2019)

[34] M. K. Gupta, S. Kumar, R. Mittal, S. K. Mishra, S. Rols, O. Delaire, A. Thamizhavel, P. Sastry, and S. L. Chaplot, Distinct anharmonic characteristics of phonon-driven lattice thermal conductivity and thermal expansion in bulk MoSe₂ and WSe₂, J. Mater. Chem. A 11(40), 21864 (2023)

[35] Y. M. Zhao, C. Zhang, S. Shin, and L. Shen, Thermal conductivity of sliding bilayer h-BN and its manipulation with strain and layer confinement, J. Mater. Chem. C 11(32), 11082 (2023)

[36] A. Ward,D. Broido,D. A. Stewart,G. Deinzer, Ab initio theory of the lattice thermal conductivity in diamond, Phys. Rev. B 80(12), 125203 (2009)

[37] L. Lindsay, D. Broido, and T. Reinecke, First-principles determination of ultrahigh thermal conductivity of boron arsenide: A competitor for diamond, Phys. Rev. Lett. 111(2), 025901 (2013)

[38] D. A. Broido, M. Malorny, G. Birner, N. Mingo, and D. Stewart, Intrinsic lattice thermal conductivity of semiconductors from first principles, Appl. Phys. Lett. 91(23), 231922 (2007)

[39] T. Feng and X. Ruan, Quantum mechanical prediction of four-phonon scattering rates and reduced thermal conductivity of solids, Phys. Rev. B 93(4), 045202 (2016)

[40] Z. Gao, F. Tao, and J. Ren, Unusually low thermal conductivity of atomically thin 2D tellurium, Nanoscale 10(27), 12997 (2018)

[41] Y. Wang and J. Ren, Strain-driven switchable thermal conductivity in ferroelastic PdSe₂, ACS Appl. Mater. Interfaces 13(29), 34724 (2021)

[42] A. Maradudin and A. Fein, Scattering of neutrons by an anharmonic crystal, Phys. Rev. 128(6), 2589 (1962)

[43] S. Hepplestone and G. Srivastava, Phonon-phonon interactions in single-wall carbon nanotubes, Phys. Rev. B 74(16), 165420 (2006)

[44] P. Klemens, Thermal conductivity and lattice vibrational modes, in: Solid State Physics, Vol. 7, Elsevier, 1958, pp 1–98

[45] J. Xie, S. de Gironcoli, S. Baroni, and M. Scheffler, First-principles calculation of the thermal properties of silver, Phys. Rev. B 59(2), 965 (1999)

[46] L. Lindsay and D. Broido, Three-phonon phase space and lattice thermal conductivity in semiconductors, J. Phys.: Condens. Matter 20(16), 165209 (2008)

[47] S. Ju, R. Yoshida, C. Liu, S. Wu, K. Hongo, T. Tadano, and J. Shiomi, Exploring diamond-like lattice thermal conductivity crystals via feature-based transfer learning, Phys. Rev. Mater. 5(5), 053801 (2021)

[48] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54(16), 11169 (1996)

[49] G. Kresse and J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a planewave basis set, Comput. Mater. Sci. 6(1), 15 (1996)

[50] P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B 50(24), 17953 (1994)

[51] G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B 59(3), 1758 (1999)

[52] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77(18), 3865 (1996)

[53] S. Baroni, S. De Gironcoli, A. Dal Corso, and P. Giannozzi, Phonons and related crystal properties from density-functional perturbation theory, Rev. Mod. Phys. 73(2), 515 (2001)

[54] A. Togo, F. Oba, and I. Tanaka, First-principles calculations of the ferroelastic transition between Rutile-type and CaCl₂-type SiO₂ at high pressures, Phys. Rev. B 78(13), 134106 (2008)

[55] B. Wei, Q. Sun, C. Li, and J. Hong, Phonon anharmonicity: A pertinent review of recent progress and perspective, Sci. China Phys. Mech. Astron. 64(11), 117001 (2021)

[56] Z. Gao, G. Liu, and J. Ren, High thermoelectric performance in two-dimensional tellurium: An ab initio study, ACS Appl. Mater. Interfaces 10(47), 40702 (2018)

[57] X. Zhou, Y. Yan, X. Lu, H. Zhu, X. Han, G. Chen, and Z. Ren, Routes for high-performance thermoelectric materials, Mater. Today 21(9), 974 (2018)

[58] Y. Liu, D. Jia, Y. Zhou, Y. Zhou, J. Zhao, Q. Li, and B. Liu, Discovery of ABO₄ Scheelites with the extra low thermal conductivity through high-throughput calculations, J. Materiomics 6(4), 702 (2020)

[59] Y. Zhuo, A. Mansouri Tehrani, A. O. Oliynyk, A. C. Duke, and J. Brgoch, Identifying an efficient, thermally robust inorganic phosphor host via machine learning, Nat. Commun. 9(1), 4377 (2018)

[60] H. Zhu, C. Zhao, P. Nan, X. Jiang, J. Zhao, B. Ge, C. Xiao, and Y. Xie, Intrinsically low lattice thermal conductivity in natural superlattice (Bi₂)_m(Bi₂Te₃)_n thermoelectric materials, Chem. Mater. 33(4), 1140 (2021)

RIGHTS & PERMISSIONS

Higher Education Press