The ultra-low lattice thermal conductivity dominated by the quartic anharmonicity in Bi-based binary compounds A₃Bi (A = K, Rb)

Thermal conductivity plays a vital role in the development of many modern sciences and technologies such as thermoelectric (TE) devices and transistors, and it can represent the heat conduction performance of materials under limited temperature gradients [1]. In the TE field, we need to minimize the thermal conductivity to get the best energy conversion efficiency, and high thermal conductivity is necessary for the heat dissipation of the devices in the reactor [2,3]. Therefore, semiconductors that exhibit extremely high or extremely low thermal conductivity, such as TE materials with low lattice thermal conductivity $(\kappa_{L})$, which have been widely studied in recent years [4–7], are worthy of our attention.

The Bi-based binary compounds are very interesting compounds and have recently attracted much attention due to their unusual electronic properties [8–11]. Studies have shown that in the low-temperature hexagonal crystal structure, these compounds possess 3D Dirac states, which are protected by crystal symmetry [12], which makes them widely used in the study of various topological properties [13–15]. At present, most of the researches are on the Bi-based binary compounds in the hexagonal phase, while the research on the cubic-phase Bi-based binary compounds is very rare, this is because the traditional harmonic approximation (HA) based on density functional theory (DFT) fails in these severe anharmonic materials.

Ultra-low $\kappa_{L}$ and inelastic neutron scattering spectra have proved that lattice anharmonic effects are universal in thermoelectric materials [16]. In the past few decades, by including lattice anharmonic effects, people have developed a variety of methods to calculate finite temperature anharmonic phonons [17,18]. We chose the self-consistent phonon (SCP) theory from these methods [19,20]. Combined with compressed sensing (CS) lattice dynamics, the interatomic force constants (IFCs) can be quickly trained according to the required displacement and force data sets [21,22]. Phonon scattering caused by lattice anharmonicity is also a key factor that affects $\kappa_{L}$. The $\kappa_{L}$ in the compound is obtained by the Boltzmann transport equation (BTE) solution with three-phonon (3ph) scattering (arising from cubic IFCs), which can successfully reveal the microscopic mechanism of high and low $\kappa_{L}$ [23–25]. However, current research has found that sometimes only considering 3ph scattering can cause $\kappa_{L}$ to be higher than the experimental value [26,27]. In this case, we further consider that four-phonon (4ph) scattering (arising from quartic IFCs) can reduce the obtained $\kappa_{L}$ [28], which is equivalent to the experimental result, indicating the necessity of fourth-order anharmonicity.

In summary, we use the first-principles calculation method combining SCP theory, CS lattice dynamics and BTE to calculate $\kappa_{L}$ containing 3ph and 4ph scattering in the Bi-based binary compounds A₃Bi (A = K, Rb), and the results show that the Bi-based binary compounds A₃Bi (A = K, Rb) of the cubic phase have extremely low $\kappa_{L}$. According to previous studies, K₃Bi and Rb₃Bi are hexagonal phases at low temperatures, and transform to cubic phases at about 550 K and 500 K [29,30].

All calculations are implemented through the density functional theory (DFT) software package Vienna Ab initio Simulation Package (VASP) [31, 32]. The projector augmented wave (PAW) [33] potentials and the exchange-correlation functional of Perdew-Burke-Ernzerhof revised for solids (PBEsol) [34] are used to emulate the ion cores and valence electrons with a cutoff energy of 600 eV. We set the criterion of the force acting on each atom to be less than $10^{-6}\ \text{eV/}\unicode{8491}$, and set the energy convergence criterion to $10^{-8}$ to fully optimize the primitive cell, where the k-point mesh is $14 \times 14 \times 14$ to sample the Brillouin zone.

According to the density functional perturbation theory (DFPT) [35], we can calculate Born effective charge and dielectric tensor required by the non-analytic part of the dynamic matrix. By encoding in the ALAMODE software package, we obtained the harmonic IFCs in the $2 \times 2 \times 2$ supercell using the finite displacement method, with a displacement of 0.01 Å [20,36,37], and further through the CS lattice dynamics approach to obtain cubic and quartic IFCs. To be more specific, we first simulated a 4000-step ab initio molecular dynamics (AIMD) at 300 K with a time step of 2 fs in a $2 \times 2 \times 2$ supercell to sample 80 snapshots. In order to obtain quasi-random configurations, each sample will move all atoms by 0.1 Å in a random direction. Then, we use static DFT in a $7 \times 7 \times 7$ k-point mesh to calculate the Hellmann-Feynman force of each quasi-random configuration. Finally, the least absolute shrinkage and selection operator (LASSO) technique is used to estimate the anharmonic IFCs using a data set composed of displacements and forces of the 80 quasi-random configurations [38,39]. In the process of calculating harmonics and third-order IFCs, all interatomic interactions in the supercell are considered, and the quartic IFCs is considered to the fifth nearest neighbor interaction; since the fifth- and sixth-order IFCs will not change the results too much, the third and second nearest neighbor interactions are considered, respectively.

Based on the SCP theory, the anharmonic frequencies are obtained through the calculated harmonic and anharmonic IFCs. Since we used $2 \times 2 \times 2$ supercell when generating IFCs, a fixed $2 \times 2 \times 2$ q-mesh is utilized in the reciprocal-to-real-space Fourier interpolation, and a $8 \times 8 \times 8$ q₁-mesh was further applied to the inner loop of the SCP equation [37], which is used to achieve the convergence of anharmonic phonon eigenvalues. The FourPhonon package is a revised version of the ShengBTE code [28,40,41], which is used to solve the phonon BTE. In the calculation process, a $12 \times 12 \times 12$ q-mesh was used to simulate the phonon wave vector. We use the iterative solution of 3ph scattering to solve the BTE, while the iterative solution of 4ph scattering is processed by single-mode relaxation time approximation (SMRTA) due to the huge computational cost.

The Bi-based binary compounds A₃Bi (A = K, Rb) in the cubic phase have the same lattice structure as shown in fig. 1(a) and (b), and their face-centered cubic unit cell belongs to the space group of $Fm\bar{3}m$. The primitive cell of A₃Bi (A = K, Rb) is composed of three A atoms (gray sphere) and one Bi atom (blue sphere) as shown in fig. 1(b), it is worth noting that only two of the three A atoms are equivalent, so we use A1 and A2 to distinguish them. The calculated lattice constants are 8.72 and 9.14 Å for K₃Bi and Rb₃Bi, respectively, as listed in table S1 of the Supplementary Material Supplementarymaterial.pdf (SM) which are similar to the experimental values [29,42] and the result of Rusinov et al. [8]. In addition, we also calculated the dielectric tensor $(\epsilon)$ and Bonn effective charges $(Z^{*})$ of A₃Bi and listed them in table S1 of the SM.

Zoom In Zoom Out Reset image size

In order to better study the vibration characteristics of the Bi-based binary compounds A₃Bi (A = K, Rb), the phonon dispersion curves of finite temperature are calculated by solving the SCP equation [39],

$$\begin{equation} V_{q j j^{\prime}}=\omega_{q j}^{2} \delta_{j j^{\prime}}+\frac{1}{2} \sum_{q_{1}} \Phi\left(\bm{q} j;-\bm{q} j^{\prime}; q_{1};-q_{1}\right)\left\langle Q_{q_{1}}^{*} Q_{q_{1}}\right\rangle. \end{equation} \tag{ 1 }$$

Here, $\omega_{q j}$ is the harmonic phonon frequency with the momentum q and the branch index j; $\Phi\left(\bm{q} j;-\bm{q} j^{\prime}; q_{1};-q_{1}\right)$ is the reciprocal representation of the fourth-order IFCs; $\left\langle Q_{q_{1}}^{*} Q_{q_{1}}\right\rangle$ is the mean square displacement of the normal coordinate Q_q . The quartic anharmonic usually increases the frequency of low-lying rattling modes, which can be attributed to the dominant and positive contributions of the diagonal term of the quartic coefficients, for the low-lying optical modes j we have $\Phi(\mathbf{0} j; \mathbf{0} j; \mathbf{0} j; \mathbf{0} j)>0$. We have drawn the phonon dispersion curves of K₃Bi and Rb₃Bi as shown in fig. 1(c) and (d), where black curve represents the harmonic phonon dispersion, and the anharmonic phonon dispersion at 600 K, 800 K and 1000 K are represented by blue, purple, and gray curves, respectively. Since we mentioned earlier that the temperatures at which K₃Bi and Rb₃Bi transform into cubic phases are 550 K and 500 K, respectively, we have only given the anharmonic phonon dispersion curves above 600 K, and we will use the calculation data at 600 K in the subsequent descriptions. The dynamic stability of materials can usually be judged by observing whether there are imaginary frequencies $(\omega_{q}^{2} < 0)$ in the phonon spectrum, and through the harmonic phonon dispersion in fig. 1(c) and (d), it is obvious that there are multiple imaginary frequencies with the lowest value of $15.73i\ \text{cm}^{-1}$ at $X (\frac{1}{2}, 0, \frac{1}{2})$, $K (\frac{3}{8}, \frac{3}{8}, \frac{3}{4})$ and $U (\frac{5}{8}, \frac{1}{4}, \frac{5}{8})$ points. Comparing the results of HA phonon density of states (PDOS) in fig. 2(a) and (c), it can be found that the imaginary frequency is mainly caused by A1 (A = K, Rb). After considering quartic anharmonic renormalization based on the SCP theory, the imaginary frequency of the SCP phonon spectrum in A₃Bi disappears at the considered temperature because the low-lying phonon mode frequency in the materials are significantly hardened, as shown in fig. 1(c) and (d), this is mainly due to the positive and dominant contribution of the diagonal term of the quartic coefficient to the eigenvalue of the phonon, which is usually manifested as the quartic anharmonic effect increasing the frequency of the low-lying phonon mode [27,37]. In addition, the SCP phonon dispersion curve did not change significantly as the temperature increased to 1000 K. By plotting the anharmonic PDOS at 600 K in fig. 2(b) and (d), it can be found that the vibration frequency of A1 is significantly increased and that A1 and the heavier Bi atoms contribute to the low-frequency PDOS, and the PDOS at intermediate and high frequencies is mainly contributed by the lighter A2 atoms. Through the stable structure and corresponding phonon properties, the lattice thermal conductivity $\kappa_{L}$ can then be calculated.

Zoom In Zoom Out Reset image size

According to the solution of phonon BTE [43], the $\kappa_{L}$ along the α-axis can be expressed as

$$\begin{equation} \kappa_{\mathrm{L}}^{\alpha \alpha}=\frac{\hbar^{2}}{k_{B} T^{2} \Omega N_{q}} \sum_{\mathbf{q} \nu} n_{\mathbf{q} \nu}\left(n_{\mathbf{q} \nu}+1\right) \omega_{\mathbf{q} \nu}^{2} v_{\mathbf{q} \nu}^{\alpha} F_{\mathbf{q} \nu}^{\alpha}, \end{equation} \tag{ 2 }$$

where ℏ and k_B are the reduced Planck constant and Boltzmann constant; T, Ω and N_q are absolute temperature, volume of the unit cell and the number of wave vectors q, respectively; $n_{\mathbf{q} \nu}$, $\omega_{\mathbf{q} \nu}$ and $v_{\mathbf{q} \nu}$ are the population, frequency and group velocity of the phonon mode $|\mathbf{q} \nu\rangle;$ and $F_{\mathbf{q} \nu}^{\alpha}$ can be written as

$$\begin{equation} F_{\mathbf{q} \nu}^{\alpha}=\tau_{\mathbf{q} \nu}\left(v_{\mathbf{q} \nu}^{\alpha}+\Delta_{\mathbf{q} \nu}\right), \end{equation} \tag{ 3 }$$

where $\tau_{\mathbf{q} \nu}$ and $\Delta_{\mathbf{q} \nu}$ are the lifetime of phonon and correction term of the phonon mode $|\mathbf{q} \nu\rangle$.

Since there is an imaginary frequency in the HA phonon dispersion curve, which indicates that it is dynamically unstable, we will not further study its related properties. The calculated SCP lattice thermal conductivity $\kappa_{3ph}$ and $\kappa_{3ph,4ph}$ of the Bi-based binary compounds A₃Bi (A = K, Rb) are plotted in fig. 3. As mentioned above, the transition temperatures of cubic phase K₃Bi and Rb₃Bi are 550 K and 500 K, respectively, so we only carry out calculation research on A₃Bi above 600 K. It can be found from fig. 3 that A₃Bi has extremely low $\kappa_{L}$ and the $\kappa_{L}$ decreases with increasing temperature. For K₃Bi, $\kappa_{3ph}$ and $\kappa_{3ph,4ph}$ are 0.501 (0.367) and 0.299 (0.214) W/mK at 600 K (1000 K), and for Rb₃Bi, $\kappa_{3ph}$ and $\kappa_{3ph,4ph}$ are 0.322 (0.289) and 0.174 (0.134) W/mK at 600 K (1000 K), which are much lower than classical TE materials PbTe [44] and PdTe [45], and can even be comparable to the recently discovered full-Heusler semiconductors [6] with extremely low $\kappa_{L}$. In addition, compared to the $\kappa_{L}$ calculated by SCP + 3ph, SCP + 3,4ph can increase the total scattering rates (SRs) without changing the phonon group velocity (v^ph ), thereby reducing the phonon lifetime, resulting in a lower $\kappa_{L}$ (a reduction of nearly 40% and 44% for K₃Bi and Rb₃Bi at 600 K). In general, a lower Debye temperature indicates a smaller $\kappa_{L}$ due to the enhanced phonon scattering at a given temperature. The Debye temperature of K₃Bi and Rb₃Bi is 125.1 K and 93.2 K, respectively, which can further confirm the extremely low $\kappa_{L}$ of K₃Bi and Rb₃Bi.

Zoom In Zoom Out Reset image size

Generally, in the HA + BTE calculations, because of the changes in the Bose distribution function, the increase of temperature results in the enhanced 3ph scattering and reduces phonon lifetime [41]. There is a universal law of $\kappa_{L} \sim T^{-1}$ in ordinary semiconductors. However, in the present SCP + BTE calculations [37], since the effect of the changes in the Bose distribution function is offset by the frequency enhancement caused by the strong quartic anharmonic renormalization, the three-phonon scattering and phonon lifetimes of the soft modes are less sensitive to the increase of temperature. Therefore, considering that the three-phonon scattering of the soft modes plays a significant role in the final $\kappa_{L}$, the temperature dependence of $\kappa_{L}$ calculated from SCP + BTE is weaker than the power law of $\kappa_{L} \sim T^{-1}$ given by conventional HP + BTE method. As listed in table 1, the temperature dependences of the $\kappa_{3ph} (\kappa_{3ph,4ph})$ for K₃Bi and Rb₃Bi are $T^{-0.51} (T^{-0.67})$ and $T^{-0.25} (T^{-0.42})$, it can be found that the temperature dependence of A₃Bi has increased to a certain extent after considering 4ph scattering, making the result closer to the trend of $T^{-1}$, which indicates that 4ph scattering occupies a larger proportion in the calculation of the $\kappa_{L}$ for A₃Bi.

Table 1:. The contribution ratio of the three acoustic phonon branches (TA1, TA2 and LA) and the optical phonon branch to $\kappa_{3ph}$ and $\kappa_{3ph,4ph}$ at 600 K. The temperature dependence of $\kappa_{3ph}$ and $\kappa_{3ph,4ph}$ for A₃Bi (A = K, Rb) is in the last column.

Material		TA1	TA2	LA	Optic	Temperature dependence
K3Bi	$\kappa_{3ph}$	30.0%	18.8%	28.7%	22.5%	$T^{-0.51}$
	$\kappa_{3ph,4ph}$	27.2%	20.8%	26.7%	25.3%	$T^{-0.67}$
Rb3Bi	$\kappa_{3ph}$	27.2%	15.5%	27.3%	30.0%	$T^{-0.25}$
	$\kappa_{3ph,4ph}$	23.6%	15.5%	24.0%	36.9%	$T^{-0.42}$

We analyzed the contribution ratio of different phonon branches to the $\kappa_{L}$ in table 1, and plotted the normalized lattice thermal conductivity $\kappa_{C}/\kappa_{L}$ of K₃Bi and Rb₃Bi in fig. 4(a) and (b). The contribution of acoustic phonons to the $\kappa_{3ph} (\kappa_{3ph,4ph})$ are 77.5% (74.7%) for K₃Bi and 70% (63.1%) for Rb₃Bi, which indicates that the $\kappa_{L}$ in A₃Bi is mainly contributed by the low-frequency acoustic phonons. From the $\kappa_{C}/\kappa_{L}$ in fig. 4(a) and (b), we can also find that the slope of $\kappa_{C}/\kappa_{L}$ first increases and then decreases as the frequency increases, indicating that low-frequency phonons occupy the main component in the contribution of the $\kappa_{L}$. It is worth noting that the acoustic branches of K₃Bi are distinguished in fig. S5 in the SM. The results show that TA is not always the branch with the lowest frequency. However, the contribution of lattice thermal conductivity of each branch can be directly output by ShengBTE code, and the hybridization of LA and TA branches does not affect the final thermal conductivity results. In addition, in order to evaluate the influence of the size effect on the $\kappa_{L}$, we discussed the $\kappa_{L}$ of phonon maximum mean-free path (MFP) allowed in the SCP + 3ph and SCP + 3,4ph models in fig. 4(c) to show the effect of the 3ph and 4ph processes on the size-dependent $\kappa_{L}$. We can see that 1) in the SCP + 3ph (SCP + 3,4ph) approximation, the total $\kappa_{L}$ increases continuously with the increase of MFP, until $\kappa_{L}$ reaches the maximum value at MPF equal to 60 (38) nm for K₃Bi, and 29 (22) nm for Rb₃Bi; 2) in K₃Bi (Rb₃Bi), 4ph scattering significantly inhibits the heat conduction of carriers when the MFP is larger than 1.3 (1.6) nm; 3) the $\kappa_{L}$ below 10 nm occupies 70% (80%) and 76% (84%) of the total $\kappa_{3ph}$ and $\kappa_{3ph,4ph}$ for K₃Bi (Rb₃Bi). Therefore, in order to significantly reduce the $\kappa_{L}$ in A₃Bi, we need to reduce the sample size below 10 nm.

Zoom In Zoom Out Reset image size

The v^ph and phonon SRs are closely related to the $\kappa_{L}$. In order to explore the difference in $\kappa_{L}$ between K₃Bi and Rb₃Bi, we will conduct a detailed analysis for the v^ph and phonon SRs of A₃Bi (A = K, Rb). Figure 4(d) shows the v^ph of A₃Bi as a function of frequency. It can be found that low-frequency phonons have a larger v^ph , which is consistent with our previous analysis of low-frequency phonons contributing more $\kappa_{L}$. In addition, the v^ph of K₃Bi is significantly higher than that of Rb₃Bi, indicating that K₃Bi may have a larger $\kappa_{L}$, which is consistent with the trend of $\kappa_{L}$ in fig. 3. Figures 4(e) and (f) show the 3ph and 4ph scattering process of A₃Bi. Recent studies have shown that the acoustic mode is less affected by the 4ph scattering process compared to the optical mode [28,46], so it is surprising that the 4ph SRs of A₃Bi is higher than the 3ph SRs at low frequencies. We have previously analyzed that the $\kappa_{L}$ of A₃Bi is mainly contributed by low-frequency phonons, while the 4ph SRs plays a major role in the low-frequency phonons, it is not surprising that $\kappa_{L}$ is severely suppressed when 4ph scattering process is considered. The gray solid line indicates that the SRs equal phonon frequency $(1/\tau = \omega/2\pi)$ as shown in fig. 4(e) and (f), which means that the phonon lifetime is equivalent to the vibration period of the phonon quasiparticle. The 3ph and 4ph SRs are mainly distributed below the phonon quasiparticle picture curve, which supports the effectiveness of the BTE solution in this research. The only known limitation of the BTE approach is that the oscillation caused by the wave vector should have a long enough lifetime that its frequency is well defined $(\tau > 2\pi/\omega)$. However, the high temperature increases the scattering rate of phonons, so that phonons may not have a long enough phonon lifetime to complete an oscillation. In order to verify the effectiveness of BTE solution at high temperature, the SRs of A₃Bi at 1000 K and phonon quasiparticle picture curve are presented in fig. S3 of the SM. Most of the scattering rates in fig. S3 are lower than the phonon quasiparticle picture curve, ensuring that A₃Bi at 1000 K has sufficient phonon lifetime when the BTE method is used for calculation, which further verifies the accuracy of our calculation results. In addition, we show the Normal and Umklapp processes of the 4ph scatterings for A₃Bi in fig. S1 of the SM. Both the Umklapp processes of K₃Bi and Rb₃Bi are higher than the Normal process, which shows the rationality of the SMRTA scheme used in the BTE solution to treat the 4ph scattering process.

In order to further understand the 4ph scattering process, we show the 3ph and 4ph scattering phase spaces (P3 and P4) of A₃Bi and the three independent processes that make up the P4 in fig. S2 the SM. Figures S2(a) and (b) of the SM show the variation of P3 and P4 with frequency for K₃Bi and Rb₃Bi, respectively. Both K₃Bi and Rb₃Bi have quite high P4 that can be comparable to P3, which can explain the high 4ph SRs in fig. 4(e) and (f). As shown in fig. S2(c) and (d) of the SM, the P4 can also be decomposed into three separate processes, namely, splitting $(\lambda \rightarrow \lambda^{\prime} + \lambda^{\prime \prime} + \lambda^{\prime \prime \prime})$, redistribution $(\lambda + \lambda^{\prime} \rightarrow \lambda^{\prime \prime} + \lambda^{\prime \prime \prime})$, and recombination $(\lambda + \lambda^{\prime} + \lambda^{\prime \prime} \rightarrow \lambda^{\prime \prime \prime})$ processes. Taking into account the constraints of phonon energy conservation in the scattering processes, we estimate that the splitting processes is mainly related to the high-frequency phonon modes, while the recombination processes mainly exist in the low-frequency phonon modes. The redistribution processes, which are more flexible in satisfying the energy conservation, are found to contribute many more 4ph scattering events, thus dominating the scattering phase space [47].

In summary, we calculated the A₃Bi (A = K, Rb) phonon dispersion curve and the PDOS using the SCP theory considering the quartic anharmonicity, and found that the strong quartic anharmonicity, which is mainly caused by the alkali metal atoms A1 in A₃Bi (A = K, Rb), hardens the vibration frequency of the low-lying phonon branch and ensures the lattice dynamic stability. Through the RTA of the BTE, the calculated $\kappa_{L}$ considering the 3ph and 4ph scattering process shows an ultra-low value, which is 0.299 (0.214) and 0.174 (0.134) W/mK for K₃Bi and Rb₃Bi at 600 K (1000 K). In addition, we also conducted a detailed analysis of the normalized lattice thermal conductivity, $\nu^{ph}$, phonon SRs, and scattering phase space, revealing that the extremely low $\kappa_{L}$ in A₃Bi is mainly due to its relatively strong 4ph scattering process. In order to further study the method of reducing the $\kappa_{L}$ of A₃Bi, we also studied the size-dependent cumulative lattice thermal conductivity of A₃Bi. The results show that controlling the sample size below 10 nm can effectively reduce the $\kappa_{L}$ of A₃Bi. Our research prove the existence of strong quartic anharmonicity in A₃Bi, and we have studied in detail the influence of 3ph and 4ph scattering processes on the thermal transport properties of A₃Bi, which has never been reported before.

This research was supported by the National Natural Science Foundation of China under Grant No. 11974302, and No. 12174327.

Data availability statement: All data that support the findings of this study are included within the article (and any supplementary files).