Introduction

Thermoelectric (TE) materials are capable of converting heat energy into electricity without additional pollution, promising to alleviate the energy shortage and environmental pollution caused by energy use1,2,3. Conventionally, the conversion efficiency of TE materials is measured by the figure of merit \(zT=\frac{{S}^{2}\sigma T}{{\kappa }_{{{{\rm{e}}}}}+{\kappa }_{{{{\rm{L}}}}}}\), in which S, σ, κe, κL, and T are the thermopower, electrical conductivity, electrical thermal conductivity, lattice thermal conductivity, and absolute temperature, respectively. In principle, high-efficiency TE materials need to possess both a high TE power factor (PF = S2σ) and low thermal conductivities (κ = κe + κL). The former involves band tuning4,5, applying strain6, and doping7,8. Meanwhile, to capture an inherently high PF, the anisotropic electronic bands and multi-valley band structure have been proposed9,10,11,12,13, in which the dispersive part leads to high σ and carrier mobility μ, while the flat (multi-valley) part induces large S. However, due to the mutual coupling of S and σ, it is difficult to improve the zT value by tuning a single parameter. To address this problem, common strategies are to suppress the κL of existing TE materials through introducing defects14,15, alloying16,17, nanostructures18, and substructures19,20, or to search for TE materials with inherently low κL.

Recently, high zT values of a series of full-Heusler compounds have been reported, such as Ba2BiAu21, Sr2SbAu22, Ca2HgPb23, and Ba2AgSb24, which is mainly attributed to the ultralow κL. Particularly, a very high zT value of 5 has been reported in full-Heusler compounds Ba2BiAu21. Hence, understanding the lattice dynamic behavior and phonon thermal transport in full-Heusler compounds is critical to gain insights into which can further suppress κL and enhance zT. However, existing theoretical studies only focus on full-Heusler compounds with lattice stability (no imaginary phonon frequencies) at 0 K. However, many materials with ferroelectric-like lattice instability also have good TE properties, such as Rb3AuO25, SnSe26, Cu2Se27, and GeTe16,28. However, the comprehensive understanding of full-Heusler compounds is hampered by the existence of imaginary phonon frequencies and the failure to consider higher-order anharmonic effects. Especially for some TE materials with imaginary phonon frequencies, the four-phonon (4ph) scattering processes even exceed three-phonon (3ph) scattering processes16,29,30.

In view of the above, we perform fist-principles calculations, self-consistent phonon (SCP) theory, and Boltzmann transport equations to investigate the thermal transport and TE properties in full-Heusler compound Na2TlSb. We consider the contribution of quartic anharmonic renormalization of phonon frequencies to phonon group velocity υph as well as the influence of 3ph and 4ph scattering to the phonon lifetimes τph. The phonon spectrum, thermal conductivity spectrum κL (ωph), 3ph and 4ph scattering processes, etc. are investigated to uncover the microscopic mechanism of thermal transport. Our key finding is that a rational κL and temperature dependence can be obtained by including full quartic anharmonicity. Meantime, the ultralow κL is captured in Na2TlSb due to the soft Tl-Sb bonding and resonant bonding in the pseudocage composed of the Na and Sb atoms interaction. Additionally, the multi-valley band structure increases the band degeneracy, resulting in a high PF in Na2TlSb. High PF combined with ultralow κL means good TE performance, with the highest zT values of 2.88 (0.94) at 300 K in p-type (n-type) full-Heusler compound Na2TlSb.

Results

Structural stability

Na2TlSb crystallizes in the face-centered cubic structure (space group Fm\(\overline{3}\)m [225]), where Na, Tl, and Sb take up the 8c, 4a, and 4b sites, respectively, as shown in Fig. 1d. The calculated lattice constants of the crystal cell is 7.485 Å, and the bond length of Na-Sb/Tl (Tl-Sb) is 3.241 (3.743) Å. The dielectric tensors ϵ and Born effective charges Z* are listed in Supplementary Table 1. To investigate the dynamics stability, we calculated the HA and anharmonic phonon dispersion curves of Na2TlSb, as shown in Fig. 1a. In HA approximations, the imaginary frequencies indicate that Na2TlSb is dynamically unstable at 0 K. However, after considering the renormalization of phonon frequencies by the quartic anharmonicity, the imaginary frequencies disappear, indicating that Na2TlSb is dynamics stable between 100 and 700 K. Furthermore, as shown in Supplementary Fig. 6, Na2TlSb is on the convex hull in the ternary phase diagram, which indicates it is the minimum free energy structure under this component. Hence, the full-Heusler compound Na2TlSb is most likely to be experimentally synthesized to form a stable structure. Additionally, a 20,000-step AIMD is simulated to estimate the stability of Na2TlSb at high temperatures, as shown in Supplementary Fig. 1. There is no significant change in free energy, indicating that Na2TlSb is stable at 700 K. Furthermore, the mechanical stability is also estimated, and the elastic constants are list in Supplementary Table 1. Na2TlSb meets the mechanical stability criteria of the cubic lattice structure, as written in

$${C}_{11}-{C}_{12} \,>\, 0,\qquad {C}_{11}+2{C}_{12} \,>\, 0,\qquad {C}_{44} \,>\, 0,$$
(1)

indicating Na2TlSb is mechanically stable31.

Fig. 1: Lattice vibration properties.
figure 1

a The temperature-dependent phonon spectrum from T = 0 to 700 K for Na2TlSb. The cyan line denotes the HA phonon dispersion (T = 0 K), and the orange, dark cyan, red, and blue lines represent the SCP dispersion at 100, 300, 500, and 700 K, respectively. b, c The element-decomposed and total phonon density of states (PDOS) at T = 0 and 300 K, respectively. d The crystal structure of Na2TlSb, as visualized in VESTA code62. The red, blue, and orange balls denote Na, Tl, and Sb atoms, respectively. e The temperature-dependent atomic mean square displacements (MSDs) from 100 to 700 K along the [100] direction. The MSD from 100 to 700 K are the result of the SCP approximations.

Phonon transport

For strongly anharmonic materials, quartic anharmonicity needs to be considered to obtain reasonable κL and temperature dependence32. First, we investigate the impact of quartic anharmonicity on phonon frequencies ωph. We find that the phonon modes below 50 cm−1 shift up significantly with increasing temperature, as shown in Fig. 1a. Figure 1b, c reveal that the vibrations of Tl atoms dominate these phonon modes with imaginary frequencies. While the optical phonon branches, mainly contributed by the vibrations of Sb and Na atoms, is shifted upward only slightly. The vibrational behavior associated with the Tl atom is similar to that of the rattling modes in clathrates and skutterudites, whose phonon frequencies exhibit small values and strong temperature dependence33. Actually, avoiding the overlap of acoustic and optical phonon modes is a clear sign of the rattling behavior of Tl atoms34. The large mean square displacements (MSDs) and weak bonding of Tl atoms also confirm the above conclusion, as shown in Figs. 1e and 5b. Additionally, the phonon frequency shifts caused by cubic anharmonicity is not included in our calculations. Generally, the frequency shifts due to cubic anharmonicity is negative (Δωq < 0) and much smaller than the frequency hardening induced by quartic anharmonicity33. If the frequency shifts caused by cubic anharmonicity is considered, the hardening of low-lying phonon mode will be slightly suppressed, resulting in a slight decrease in the κL.

The imaginary HA phonon frequency will induce the phase transition of Na2TlSb at low temperatures. To predict the phase transition temperature Tc of Na2TlSb, we also calculated the temperature dependence of the squared frequency of the lowest acoustic phonon mode at X point, which is the softest mode in the region away from the Brillouin zone center, as shown in Supplementary Fig. 4. Above approximately room temperature, the temperature dependence of squared frequency can be rationally fitted by equation \({\Omega }_{{{{\bf{q}}}}}(T)={a}{({T}-{{T}}_{{{{\rm{c}}}}})}^{2}\)35, where Tc is the phase transition temperature. Applying the above equation, we obtain the phase transition temperature Tc as 83 K.

We calculate the temperature-dependence lattice thermal conductivities \({\kappa }_{{{{\rm{3ph}}}}}^{{{{\rm{SCP}}}}}\) and \({\kappa }_{3,4{{{\rm{ph}}}}}^{{{{\rm{SCP}}}}}\) from 100 to 700 K for Na2TlSb, as plotted in Fig. 2a. Due to the existence of imaginary frequencies, the solution of κL under HA approximations is numerically invalid. Hence, only the κL results of SCP approximations are provided. Compared with \({\kappa }_{{{{\rm{3ph}}}}}^{{{{\rm{SCP}}}}}\), \({\kappa }_{3,4{{{\rm{ph}}}}}^{{{{\rm{SCP}}}}}\) is significantly reduced due to the effect of 4ph scatterings. The calculated \({\kappa }_{{{{\rm{3ph}}}}}^{{{{\rm{SCP}}}}}\) and \({\kappa }_{3,4{{{\rm{ph}}}}}^{{{{\rm{SCP}}}}}\) are 0.91 and 0.44 Wm−1 K−1 at 300 K for Na2TlSb, respectively. The values of \({\kappa }_{3,4{{{\rm{ph}}}}}^{{{{\rm{SCP}}}}}\) is ultralow, half that of quartz glass (κL ~0.9 WK−1 m−1). To further explore the microscopic mechanism of heat conduction, we use the κL ~ Tα to analysis the temperature dependence of κL. We evaluated the temperature dependence of κL from 200 to 700 K, since κL at extremely low temperature is mainly determined by the lattice-specific heat CV following the Debye T3 law29. The \({\kappa }_{{{{\rm{3ph}}}}}^{{{{\rm{SCP}}}}}\) exhibits anomalously temperature dependence, \({\kappa }_{{{{\rm{3ph}}}}}^{{{{\rm{SCP}}}}} \sim {T}\)−0.52, which deviates far from the universal law κL ~ T−1. By including 4ph scattering, the temperature dependence of the \({\kappa }_{3,4{{{\rm{ph}}}}}^{{{{\rm{SCP}}}}}\) is enhanced and the value of α becomes 0.84. The increased temperature dependence of \({\kappa }_{3,4{{{\rm{ph}}}}}^{{{{\rm{SCP}}}}}\) further confirms the importance of the quartic anharmonicity for the materials with imaginary frequencies. Additionally, thermal expansion and additional phonon-grain boundary scattering also play a crucial role in determining κL and corresponding temperature dependence32,36. If the above factors are taken into account when estimating the phonon thermal transport properties, κL will decrease and temperature dependence of κL will be closer to the experimental results. Furthermore, we verify the feasibility of excluding the thermal expansion and cubic anharmonicity, we calculated the \({\kappa }_{3,4{{{\rm{ph}}}}}^{{{{\rm{SCP}}}}}\) considering the thermal expansion and discussed the effect of the cubic anharmonicity on the phonon frequencies. As shown in Supplementary Figs. 7, 8, \({\kappa }_{3,4{{{\rm{ph}}}}}^{{{{\rm{SCP}}}}}\) considering thermal expansion at 700 K is 4.1% lower than that without considering thermal expansion. Moreover, we observe that the phonon frequency shift caused by cubic anharmonicity at X point is very small, only 4.21 cm−1 at 700 K, and this value will not significantly change the phonon spectrum. Hence, it is rational to exclude the above factors to calculate κL.

Fig. 2: Phonon transport properties.
figure 2

a The calculated temperature-dependent κL from 100 to 700 K by using SCP+3ph and SCP+3,4ph, respectively. b The calculated 3ph (blue empty circles) and 4ph (orange empty circles) scattering rates (SRs) in the SCP approximations. The purple line denotes 1/τph = ωph/2π, that is, the scattering rates equal to phonon frequencies. c The κL spectrum (The filled area below the lines) and cumulative κL at 300 K by using SCP+3ph and SCP+3,4ph, respectively. d Same as c but for the maximum mean free path (MFP) cumulative κL at 300 K.

The large 3ph and 4ph SRs are one of the main reasons for the ultralow κL of Na2TlSb, as shown in Fig. 2b. Similar to 3ph scattering, 4ph scattering is also important in phonon thermal transport. It is evident from the calculations that 4ph SRs even exceeds 3ph SRs below 50 cm−1. The strong 4ph SRs can be attributed to the strict constraints on the 3ph scattering phase space by the selection rule37. Similar to BAs, the 3ph scattering phase spaces are restricted by the selection rule due to the large phonon band gap, resulting in a weakened 3ph SRs in Na2TlSb. Meantime, the selection rule has little effect on the 4ph scattering phase spaces, resulting in large 4ph SRs. The SRs curve equal phonon frequencies is also plotted by the purple line in Fig. 2b, which means that the phonon lifetime τph is equal to the vibrational period of the phonon quasiparticle. The phonon quasiparticle image is invalid if the τph is less than one vibrational period, i.e., the phonon annihilates before completing one vibration29. As shown in Fig. 2b, all 3ph and most 4ph scattering distributions are below the curve, indicating that our BTE scheme is valid. The κL spectrum and frequency cumulative κL indicate that the major contributions to the κL are the phonon branches below 50 cm−1. Meantime, the phonon branches below 50 cm−1 are the major contributions to the suppressed κL calculated with SCP+3,4ph relative to SCP+3ph, which is consistent with strong 4ph SRs below 50 cm−1. Figure 2d show the cumulative κL as the function of maximum MFP. The result shows that heat-carrying phonons have considerable MFPs at 300 K. The maximum MFP of \({\kappa }_{{{{\rm{3ph}}}}}^{{{{\rm{SCP}}}}}\) is about 210 nm, while the maximum MFP of \({\kappa }_{3,4{{{\rm{ph}}}}}^{{{{\rm{SCP}}}}}\) is only about 80 nm. Through nanostructuring, the \({\kappa }_{3,4{{{\rm{ph}}}}}^{{{{\rm{SCP}}}}}\) can be lower, e.g., as low as 0.27 WK−1 m−1 when the maximum MFP is limited to 10 nm.

The phonon group velocity υph at 100 and 300 K is shown in Fig. 3a. The small υph further confirms the ultralow κL of Na2TlSb. Relative to 100 K, the υph of 300 K is slightly larger below 50 cm−1, which is consistent with the change of the phonon dispersion curves form 100 to 300 K. The observed small specific heats CV further indicates the existence of ultralow κL, as shown in Fig. 3b. Figure 3c show the decomposed 3ph scattering rate. The combination (splitting) processes dominate the low-frequency (high-frequency) 3ph SRs below 100 cm−1. Due to the large phonon band gap from 100 to 130 cm−1, it is difficult for low-frequency phonon modes to be excited to high-frequency phonon modes through the combination processes. To obey the energy conservation constraint, high-frequency phonons can only be generated by the three-phonon scattering process of high-frequency phonons. Hence, the 3ph SRs of high-frequency phonons exhibit the same regularity as those of low-frequency phonons, e.g., the 3ph SRs of low-frequency (high-frequency) modes above 130 cm−1 are mainly determined by the combination (splitting) processes. Figure 3d show the decomposed 4ph scattering rates. The redistribution processes dominate the 4ph scattering processes because it is easier to satisfy the selection rule37. With respect to the 4ph scattering processes, due to the strong redistribution processes, the splitting processes at low frequencies is suppressed and weaker than the combination processes. It is only at very high phonon frequencies that the splitting processes are larger than the combination processes.

Fig. 3: Phonon transport parameter.
figure 3

a The phonon group velocity υph at 100 and 300 K. b The calculated specific heats CV. c Decomposed 3ph scattering rates into the slitting (\(\lambda \to {\lambda }^{{\prime} }+{\lambda }^{{\prime\prime} }\)) and combination (\(\lambda +{\lambda }^{{\prime} }\to {\lambda }^{{\prime\prime} }\)) processes at 300 K. d Decomposed 4ph scattering rates into the splitting (\(\lambda \to {\lambda }^{{\prime} }+{\lambda }^{{\prime\prime} }+{\lambda }^{{\prime\prime\prime} }\)), redistribution (\(\lambda +{\lambda }^{{\prime} }\to {\lambda }^{{\prime\prime} }+{\lambda }^{{\prime\prime\prime} }\)), and combination (\(\lambda +{\lambda }^{{\prime} }+{\lambda }^{{\prime\prime} }\to {\lambda }^{{\prime\prime\prime} }\)) processes at 300 K. All of the above results were obtained using SCP method.

To further understand the microscope mechanism of the ultralow κL, we calculated the Gr\(\ddot{u}\)neisen parameter γ and 3ph scattering phase space W3ph, as shown in Fig. 4a, b. Compared with the optical phonon modes, the acoustic phonon modes have larger γ, indicating that the acoustic phonon modes have stronger cubic anharmonicity. This result is consistent with our previous analysis that the rattling behavior of Tl atoms have stronger anharmonicity. Figure 4b, c show the large 3ph and 4ph scattering phase spaces, indicating that Na2TlSb have large anharmonic scattering rates. Meanwhile, The decomposed anharmonic (3ph and 4ph) scattering phase spaces exhibit a similar law to the decomposed anharmonic scattering rates. Figure 4d show the decomposed 4ph SRs into normal and Umklapp processes. The Umklapp processes dominate the 4ph SRs of the entire phonon modes, indicating that the 4ph scattering processes mainly suppress heat conduction, which is consistent with the above analysis of κL. The Umklapp processes are much larger than the normal processes, which also confirms the above analysis, as shown in the inset of Fig. 4d.

Fig. 4: Phonon scattering.
figure 4

a The Gr\(\ddot{u}\)neisen parameter γ of acoustic (violet) and optical (orange) phonon branches. b The decomposed 3ph scattering phase space W3ph into the splitting and combination processes. c The decomposed 4ph scattering phase space W4ph into the splitting, redistribution, and combination processes. d Decomposed 4ph scattering rates into the normal and Umklapp processes. The inset represents the ratios between Umklapp and normal processes.

Electronic structure

Figure 5a show the electronic band structure and density of states (DOS) of Na2TlSb. Na2TlSb is an indirect band gap semiconductor with the conduction band minimum (CBM) at the L point and the valence band maximum at the K − Γ high symmetry line. The calculated band gap is 0.08 eV using the PBEsol functional with SOC, which is in the range commonly found in TE materials. In Na2TlSb, the CBM is mainly contributed by Sb and Tl atoms, while the VBM is almost entirely contributed by Sb atoms. The CBM has high electronic band dispersion, which will lead to high electrical conductivity σ. Meantime, the observed remarkable band asymmetry at CBM further manifests a good TE performance in n-type Na2TlSb38. A multi-valley band structure at VBM results in high degeneracy. Additionally, the second (Γ − X line) and third (W point) highest valleys of the valence band are very close to the VBM, by doping, we can achieve a higher degeneracy. A large degree of degeneracy means that p-type Na2TlSb has a high Seebeck coefficient S. Meanwhile, multiple valleys create additional conduction channels, resulting in high σ. Figure 5b show the projected crystal orbital Hamilton population (pCOHP) analysis. At the VBM, Na−Sb and Tl−Sb show bonding states, while Na−Tl presents antibonding states. Obviously, the Tl atoms are weakly bound to other atoms, which is consistent with the MSD results.

Fig. 5: Electron structure.
figure 5

a The projected electronic structure and partial electronic density of states for Na2TlSb using the PBEsol functional. b The calculated projected crystal orbital Hamilton population (pCOHP) analysis for the Na-Sb, Na-Tl, and Tl-Sb of Na2TlSb. The negative values of pCOHP present antibonding states. The Fermi energy level is put in 0 eV.

Since the band gap of materials is usually underestimated using PBE functional, we recalculated the band gap accurately with HSE06 functional with SOC, as shown in Supplementary Fig. 2. The band gap of HSE06 functional with SOC is 0.30 eV. Since the electronic structures of HSE06 and PBEsol functional are not obviously different except for the band gap for Na2TlSb, the band gap of the HSE06 functional is used to obtain rational electronic transport properties. The carrier relaxation time τ of Na2TlSb is calculated by considering the ADP, IMP, and POP scattering mechanisms, as shown in Supplementary Fig. 3.

Electron transport

Figures 6,7 show the electronic transport parameters for p-type and n-type Na2TlSb, respectively. We observe that σ is proportional to carrier concentration n and inversely proportional to temperature T. The former is due to an increase in the number of carriers participating in the conduction process due to an increase in the concentration n. The latter is attributed to the increased scattering rates dominated by the electron-phonon interaction caused by the temperature increase, which is consistent with the results in Supplementary Fig. 3. The n-type Na2TlSb exhibits higher σ relative to the p-type Na2TlSb, which is consistent with our previous analysis of the larger electronic band dispersion at the CBM. Unlike σ, the S decreases with increasing n at the same T, while the S increases with increasing T at the same n. Apparently, p-type Na2TlSb exhibits a large S due to the large degeneracy at the VBM. In general, κe and σ vary in the same law with n due to the increase in heat carrier, which is consistent with the trend in Fig. 6c and d. Due to the coexistence of a large S and σ, we obtain a high TE power factor. Particularly, S have single peak at 500 and 700 K due to the bipolar effect39. The above phenomenon is also observed in PbTe and PbSe with small band gaps40,41. The strong bipolar effect severely suppresses the S, leading to a decrease in TE performance. As a result, the highest PF is obtained at 500 K because the bipolar effect increases with temperature. At the optimal doping concentration n, the power factor is 2.26 (2.92) m Wm−1 K−2 for n-type Na2TlSb at 300 (500) K. Due to the multi-valley structure at the VBM, p-type Na2TlSb has a large power factor, e.g., 9.84 and 12.76 mW m−1 K−2 at 300 and 500 K, respectively. Moreover, the strong quartic anharmonicity of Na2TlSb affects not only phonon thermal transport properties but also electronic transport properties. Concretely, strong quartic anharmonicity leads to a hardening of phonon frequency, thereby reducing the electron-phonon coupling strength. Furthermore, the reduced electron-phonon coupling strength results in higher carrier mobility and larger electrical conductivity42. Hence, if the influence of quartic anharmonicity on electronic transport is included, the TE performance of Na2TlSb will be further improved.

Fig. 6: Electronic transport parameters.
figure 6

The calculated electronic transport parameters for p-type Na2TlSb at 300, 500, and 700 K. a Electrical conductivity σ. b The Seebeck coefficient S. c The electronic thermal conductivity κe. d TE power factor. The SOC is also included in the calculations.

Fig. 7: Electronic transport parameters.
figure 7

The calculated electronic transport parameters for n-type Na2TlSb at 300, 500, and 700 K, respectively. a Electrical conductivity σ. b The Seebeck coefficient S. c The electronic thermal conductivity κe. d TE power factor. The SOC is also included in the calculations.

The combination of ultralow κL and high TE power factor in the Na2TlSb captures an anomalously high zT ~4.81 for hole doping at nh ~ 7 × 1019 cm−3 and 500 K, as shown in Fig. 8a. Meanwhile, we predict the zT ~ 2.88 at nh ~ 4 × 1019 cm−3 and 300 K, which is also ultrahigh value for almost TE materials at room temperature. Due to the existing of large electronic band dispersion and remarkable band asymmetry at CBM, the zT is also very high for n-type Na2TlSb, e.g., 0.94 (1.48) at ne ~ 3 × 1018 (7 × 1018) cm−3 and 300 (500) K, as shown in Fig. 8b. The optimum ne for n-type Na2TlSb is far small related to nh in the p-type case, which means that n-type case is easier to achieve good TE performance. Hence, we also recommend n-type Na2TlSb as a potential TE material. Additionally, the zT value calculated based on the result of \({\kappa }_{{{{\rm{3ph}}}}}^{{{{\rm{SCP}}}}}\) is also given as the lower limit. Nonetheless, Na2TlSb also exhibits good TE properties, e.g., 1.05 (3.72) for n-type (p-type) case at 500 K. Additionally, thermal radiation43, air-induced thermal convection44, and the effect of grain boundary scattering on carrier mobility45 can lead to degradation of TE performance in experiments. Hence, if the above factors are taken into account, the calculated zT value will be slightly lower. Moreover, we only use the BTE to compute the electronic transport properties of Na2TlSb at specific doping concentrations without considering specific dopability, which also means that further experimental and theoretical explorations are required. Additionally, to determine the thermodynamic limit of achievable dopant concentrations, we estimated defect solubility for the formation of p-type and n-type semiconductors. We considered several possible neutral defects, including Na vacancy, Tl vacancy, Au in place of Tl (AuTl), Hg in place of Tl (HgTl), and Pb in place of Tl (PbTl). The calculated defect formation energy (\(\Delta {{E}}_{{{{\rm{F}}}}}^{{{{\rm{def}}}}}\)) of Na (Tl) vacancy in Na2TlSb is 0.10 (0.27) eV/defect in Na (Tl) poor condition. The above values are comparable to Na-doped PbTe (0.27 eV/defect)46, where the hole doping concentration can achieve 1020 cm−3 at room temperature41,47. The calculations indicated that the optimum hole concentration can be obtained in Na2TlSb. Furthermore, the \(\Delta {{E}}_{{{{\rm{F}}}}}^{{{{\rm{def}}}}}\) of AuTl (HgTl) is 0.75 (0.19) eV/defect. These values suggested that HgTl is more likely to reach the optimum value of hole concentration, while AuTl is difficult. Actually, since Au (Hg) are as heavy as Tl, their introduction is expected to maintain the Tl-dominated soft phonons and strong anharmonicity, thereby preserving an ultralow κL in p-type Na2TlSb. For these reasons, the use of AuTl and HgTl would be a strategic doping mechanism to best trigger p-type performance despite the toxicity of Hg and the low doping solubility of Au. Similarly, the introduction of PbTl is also expected to maintain ultralow κL in n-type Na2TlSb. Hence, we also calculated the \(\Delta {{E}}_{{{{\rm{F}}}}}^{{{{\rm{def}}}}}\) of PbTl to be −0.18 eV/defect, which indicates that PbTl is easier to achieve the optimal doping concentration. The above results strongly indicate that the doping concentration required for optimal zT can be achieved through defects.

Fig. 8: TE performance.
figure 8

The calculated TE figure of merit with the \({\kappa }_{{{{\rm{3ph}}}}}^{{{{\rm{SCP}}}}}\) and \({\kappa }_{3,4{{{\rm{ph}}}}}^{{{{\rm{SCP}}}}}\) for a p-type and b n-type Na2TlSb at 300, 500, and 700 K. The SOC is also included in the calculations.

Discussion

In summary, we investigate the thermal and TE transport properties in Na2TlSb using the first-principles calculations combined with SCP theory and Boltzmann transport equations, which explicitly include the phonon frequency shift and 4ph scattering caused by quartic anharmonicity. Our results indicate that the ultralow κL of Na2TlSb can be explained by the small υph and strong 3ph and 4ph scattering. The Tl atom with rattling behavior has strong temperature dependence and cubic and quartic anharmonicity, which plays an important role in the phonon without imaginary frequencies. Additionally, the low-frequency four-phonon scattering rates in Na2TlSb can match or even exceed the three-phonon scattering rates. Meanwhile, the multi-valley band structure at VBM increases the band degeneracy, resulting in a high TE power factor in p-type Na2TlSb. Additionally, due to large electronic band dispersion and remarkable band asymmetry at CBM, the n-type Na2TlSb exhibit a high σ. By considering ADP, POP, and IMP scattering, we capture a rational electronic relaxation time and transport properties. As a consequence, the n-type and p-type Na2TlSb show a high TE figure of merit, whose values are 1.48 and 4.81 at the optimal carrier concentration and 500 K. Our study reveals the important role of quartic anharmonicity in phonon thermal transport, which contributes to our comprehensive understandings of ultralow κL microscopic mechanism in full-Heusler compounds. At the same time, we also provide ideas for the rational design of high-performance TE materials.

Methods

First-principles calculation and CSLD

We perform first-principle calculations employing the VASP code48,49, using the plane-wave basis set and projector augmented-wave method to simulate the potentials of ions cores and valence electrons50. The exchange-correlational interactions is dealt by the Perdew-Burke-Ernzerhof revised for solids (PBEsol) functional51 of the generalized gradient approximation (GGA)52. We use 520 eV as the kinetic energy cut-off, and Γ-centered 12 × 12 × 12 k-point meshes to sample the whole Brillouin zone. The structure of Na2TlSb is fully relaxed until the energy and Hellmann-Feynman force convergence criterion are less than 10−8 eV and 1 × 10−4 eV Å−1, respectively. Throughout the thermal transport calculations, we consider the nonanalytic part of the dynamics matrix, which is derived using the dielectric tensor ϵ and Born effective charges Z* calculated by the density functional perturbation theory (DFPT)53. All required harmonic (HA) and anharmonic interatomic force constants (IFCs) are calculated in the 2 × 2 × 2 supercell with 6 × 6 × 6 k-point meshes, implemented in the ALAMODE package35,54. Additionally, we also calculated the HA phonon dispersion curves within the 2 × 2 × 2 and 3 × 3 × 3 supercells to check the convergence of IFCs, as shown in Supplementary Fig. 4. It can be observed that the HA phonon dispersion curves agree well with each other, indicating the harmonic IFCs exacted form the 2 × 2 × 2 supercell are enough to capture convergence results. Furthermore, it can be deduced that anharmonic IFCs exacted from 2 × 2 × 2 supercell are also convergent, since anharmonic IFCs generally converge faster than harmonic IFCs. Specifically, the HA IFCs were extracted in three configurations produced by the finite-displacement method55. To obtain the displacement and force datasets required for anharmonic (cubic and quartic) IFCs, we use 4000-step ab initio molecular dynamics (AIMD) simulation with 2 fs time step and 300 K temperature to capture 80 snapshots first. On this basis, we obtain 80 quasi-random configurations by applying the random displacement of 0.1 Å to each atom in the 80 snapshots. Finally, we extract displacement and force information obtained from 80 quasi-random configurations to derive anharmonic IFCs via the compressive sensing lattice dynamics (CSLD) method56. In the CSLD calculations, we consider all (third-) nearest neighbor interactions for third- (fourth- to sixth-) order IFCs. Our calculations use 80 quasi-random configurations, more than used in previous work on Tl3VSe457 and cubic SrTiO335, are sufficient to capture accurate aharmonic IFCs that produce convergent results.

Calculation of transport properties

We calculate the temperature-dependence anharmonic phonon energy eigenvalues, including the off-diagonal elements of phonon self-energies. On top of SCP calculations, the thermal transport parameters is solved based on the phonon BTE, as employed in the FourPhonon package58,59,60. Generally, the computational cost of 4ph scattering is expensive in most materials. Hence, to ensure sufficient convergence, we use the available 4ph phase spaces as the criterion for calculating the 4ph SRs. We use 12 × 12 × 12 q-mesh to capture the κL, and the numbers of available 3ph and 4ph scattering processes have reached approximately 1.8 × 106 and 1.6 × 1010, respectively. The numbers of available 4ph scattering processes are much larger than LiCoO260, indicating that present calculations are sufficient to obtain an accurate κL. Finally, we obtained the κL, which is defined as

$${\kappa }_{{{{\rm{L}}}}}=\frac{{\hslash }^{2}}{{k}_{{{{\rm{B}}}}}{T}^{2}V{N}_{q}}\mathop{\sum}\limits_{{{{\bf{q}}}}\nu }{n}_{{{{\bf{q}}}}\nu }\left({n}_{{{{\bf{q}}}}\nu }+1\right){\omega }_{{{{\bf{q}}}}\nu }^{2}{\upsilon }_{{{{\bf{q}}}}\nu }{F}_{{{{\bf{q}}}}\nu },$$
(2)

where kB is the Boltzmann constant, is the reduced Planck’s constant, V is the volume of the unit cell, Nq is the number of wave vectors, q and ν are signs of the phonon modes, ωqν is frequency, υqν is the group velocity. The Fqν is defined as

$${F}_{{{{\bf{q}}}}\nu }={\tau }_{{{{\bf{q}}}}\nu }\left({\upsilon }_{{{{\bf{q}}}}\nu }+{\Delta }_{{{{\bf{q}}}}\nu }\right),$$
(3)

where τqν and Δqν are the phonon lifetime of single-mode relaxation time approximation and quantity displaying the population deviation of the iterative solution. To examine the stability of Na2TlSb at 700 K, we performed 20,000-step AIMD simulations with a time step of 2 fs.

The electronic band structure, crystal orbital Hamilton population (COHP), high-frequency dielectric constants ϵ, and deformation potentials are calculated using PBEsol functional with 12 × 12 × 12 k-point meshes. The elastic constants Cij, static dielectric constants ϵs, and effective polar phonon frequency ωpo are calculated using DFPT. The above materials’ parameters are listed in Supplementary Table 1. The accurate band gap is obtained using HSE06 functional and 12 × 12 × 12 k-point meshes. The electronic band structure is recalculated in uniform 135 × 135 × 135 k-point grids to obtain the electron relaxation time τe and electronic transport parameters, as performed in AMSET code61. Since Na2TlSb contains a heavy Tl element, the spin-orbit coupling (SOC) is also included in the calculations of electronic transport properties. The τe is calculated by including the fully anisotropic acoustic deformation potential (ADP) scattering, polar optical phonon (POP) scattering, and ionized impurity (IMP) scattering. The electron-acoustic phonon and electron-optical phonon interaction is treated by ADP and POP scattering. For details of scattering rates, please refer to the Supplementary Methods.