Realizing a thermoelectric conversion efficiency of 10% with long-term stability in a kesterite Cu₂ZnSn(S_1−xSe_x)₄ single-leg device

Abstract

The application of thermoelectric (TE) technology, contributing to carbon neutrality by the recovery of global waste heat, is hindered by low conversion efficiencies and high costs, demonstrating a strong demand for high-performance TE materials with environmentally friendly characteristics. Multinary Cu₂ZnSnS₄ (CZTS)-based materials have attracted considerable attention for TE power generation owing to their cost-effectiveness and earth-abundant features. However, their dimensionless figure of merit (ZT), which is the critical measure of TE performance, is relatively lower than that of other TE-material-based systems. Herein, we report the discovery of p-type Cu₂ZnSn(S_1−xSe_x)₄ (CZTSSe) single crystals with a record-high ZT of 1.9 at 800 K and high TE conversion efficiency of ∼10% at a temperature difference of 473 K with long-term stability over 1000 h. The exceptional performance originates from a high power factor and intrinsically low thermal conductivity owing to the pseudo-cubic structure, alkali element doping, and short-range lattice distortion. This study demonstrates the potential of CZTS-based TE materials for intermediate-temperature TE harvesting.

---

Introduction

Carbon neutrality, that is, the realization of net-zero CO₂ emissions, is attracting global attention owing to the depletion of non-renewable energy sources driven by increasing energy consumption and the resultant environmental deterioration. The generation of clean energy is critical to achieving carbon neutrality by the mid-21st century. Thermoelectric (TE) technology shows promise in enabling reversible conversion between thermal energy, such as industrial waste heat and automobile exhaust gas, and electricity. The total amount of global waste heat is roughly estimated at the order of ‘exajoules’,^1,2 which is comparable to the energy output of several hundreds of nuclear power plants. In particular, high-temperature waste heat in the range of >500 °C is emitted from the iron and steel industry. This energy- and emission-intensive sector accounts for ∼8% of annual global anthropogenic CO₂ emissions (>2800 Mt),³ to whose decarbonization the recovery of high-temperature waste heat contributes.

TE applications with high conversion efficiency have been developed for the recovery of waste heat over various temperature ranges, which are based on Te-related compounds. High-performance and cost-effective TE materials are required for industrial applications toward carbon neutrality. However, state-of-the-art TE materials typically comprise expensive, scarce, or toxic elements, such as Pb, Bi, Sb, Ge, and Te.⁴ Recently, sulfide and selenide TE materials have been receiving much attention because the earth abundances of S and Se (420 ppm and 0.05 ppm, respectively) are higher than that of Te (0.001 ppm).⁵ A comparison of the costs of chalcogenides (Fig. 1A) indicates that the cost of sulfide materials is an order of magnitude lower than that of Te element-based materials; however, their performance is lower than that of Te element-based state-of-the-art chalcogenides.⁶ Te-free technologies must be developed for large-scale fabrication of TE devices.

Fig. 1 High TE conversion efficiency of a kesterite-based single-leg device. (A) comparison of the cost of state-of-the-art chalcogenide compounds and kesterite Cu₂ZnSn(S_1−xSe_x)₄ (CZTSSe) in this study. (B) Schematic of a CZTSSe-based single-leg TE system. (C) Demonstration of high TE efficiency of ∼10% and long-term stability over 1000 h after an aging test at 800 K. TE performance of the single-leg device was characterized over 30 cycles at a temperature difference of 473 K after each aging time. The inset depicts the photograph of the CZTSSe single-leg after a 1000 h aging test.

The high TE conversion efficiency is quantified using the dimensionless figure of merit ZT = S²σT/κ, where S, σ, T, and κ are the Seebeck coefficient, electrical conductivity, temperature, and thermal conductivity, respectively. Over the past few decades, great efforts have been devoted to improving the ZT value for various TE materials such as PbTe, GeTe, Bi₂Te₃, filled skutterudites, and half-Heuslers.⁴ Current proven strategies are typified by the concepts of phonon-glass electron-crystals⁷ and phonon-liquid electron-crystals,⁸ electronic band engineering for improvement of the Seebeck coefficient,^9–12 and nanostructure modification to reduce thermal conductivity or enhance electrical properties.^13–18 The pseudo-cubic structure has great potential to maximize the ZT value, as it leads to low crystal-field splitting and, consequently, a high power factor (PF) owing to valley degeneracy at the valence band maximum.¹⁹ We demonstrated a high ZT value of 1.6 realized in a totally environmentally benign p-type kesterite Cu₂ZnSnS₄ (CZTS) single crystal based on a pseudo-cubic strategy.^20,21

Doping techniques have contributed to the development of various materials including high performance TEs and photovoltaics (PVs). In particular, it has been recognized that the incorporation of alkali elements (alkali effect) is crucial to obtaining high efficiency for improving the electrical properties of PV devices based on chalcogenide materials such as CZTS and Cu(In,Ga)Se₂.^22,23 In the field of TE materials, the higher Na concentration leads to band convergence of valence bands in PbTe suppressing bipolar conduction and increasing the Seebeck coefficient.²⁴ Furthermore, alkali element doping in CZTS leads to high TE performance owing to controlling point defects and enhancement in the electronic density of states (DOS), respectively.^21,25 Although the preliminary results of the alkali effect indicate a positive influence on the properties of PV and TE materials, the details such as mechanism and local structure are still unknown.

Considering the remarkable potential of thermal energy utilization, research on practical aspects of TE devices/modules is attracting greater attention than traditional studies on achieving higher ZT. Developing a TE device is highly challenging and thus, from an industrial viewpoint, progress in this area remains slow because it is difficult to optimize device structure. A practical TE module consists of two or more p- and n-type TE legs that are electrically connected by a metal electrode in series. Currently, practical TE device application requires overcoming some issues such as TE material coupling, contact resistance between the TE material and electrodes, mechanical properties including thermal expansion, and long-term stability at operating temperatures. Considering current TE device progress, the development of the single-leg device is useful to evaluate the potential output power and conversion efficiency leading to optimal device structure for the module. However, there are few practical TE single-leg devices even though TE materials with a high ZT value are applied. The p-type GeTe alloy-based TE single-leg with a high ZT of 2.5 indicates a conversion efficiency of 14% at a temperature gradient of 440 K.²⁶ The n-type Mg₃Sb_1.5Bi_0.5-based TE single-leg with a high ZT of 1.5 indicates the conversion efficiency of 13% at a temperature gradient of 495 K.²⁷ In addition, several TE single-legs demonstrate the conversion efficiency of 3–10% at a temperature gradient (ΔT) range of 164–600 K based on p-type CZTS with ZT of 1.6,²⁸ p-type SnS with ZT of 1.6,⁵ p-type SnSe with ZT of 2.4,²⁹ p-type AgSbTe₂ with ZT of 1.77,³⁰ p-type α-MgAgSb with ZT of 1.2,³¹ p-type half-Heuslers ZrCoBi with ZT of 1.42,³² p-type PbSe with ZT of 1.1,³³ p-type PbTe with ZT of 2.2,³⁴ p-type Bi_0.4Sb_1.6Te₃ with ZT of 1.35,³⁵ and n-type PbTe with ZT of 1.3.³⁴ These conversion efficiencies lag behind their theoretical equivalents because of the relatively large interfacial resistance of the TE material–electrode interface.^4,27 Thus, interface properties significantly affect the application of newly fabricated TE materials.

One further step that has to be taken to realize efficient utilization of waste heat sources in future carbon-neutral society is the fabrication of high efficiency TE devices with environmentally friendly characteristics, using various material strategies. In this study, a noteworthy TE performance comparable to that of state-of-the-art TE materials was achieved by revealing the alkali effect and anion substitution to simultaneously optimize the electrical and TE properties in CZTS. A record-high ZT of ∼1.9 at 800 K was achieved using p-type K-doped Cu₂ZnSn(S_1−xSe_x)₄ (CZTSSe) single crystals. The CZTSSe/Mo(SSe)₂-barrier-layer/Mo-contact single-leg TE device shown in Fig. 1B exhibited decent ohmic properties and chemical stability between the TE material and Cu electrode, enabling high energy-conversion efficiency and long-term stability. A maximum conversion efficiency (η_max) of 10.4% and output power (P_max) of 160 mW at a temperature difference of 473 K was achieved, demonstrating the remarkable potential of CZTSSe for low-cost TE power generation (Fig. 1C). The reduction in efficiency and output power within 10% was found after the thermal aging test of the CZTSSe single-leg over 1000 h, demonstrating excellent long-term stability (Fig. 1D).

Results and discussion

Fundamental properties of CZTSSe as a TE material

In the present study, 5% Cu-deficient and 0.1 mol% K-doped samples were fabricated with different anion compositions (x), including non-doped Cu_1.9ZnSn(S_1−xSe_x)₄ (x = 0.2, 0.5, and 0.8 for samples 1, 2, and 3, respectively) and K-doped Cu_1.9ZnSn(S_1−xSe_x)₄ (x = 0.5 and 1.0 for samples 4 and 5, respectively), using the Sn-solvent-based traveling heater method (THM).³⁶ The single crystals were cut in the transverse (⊥) and longitudinal (∥) directions (perpendicular and parallel to the c-axis, respectively) for TE measurements. Their structural and compositional properties are shown in Fig. S1–S3 and Tables S1, S2. High-quality p-type CZTSSe single crystals with a homogeneous composition were obtained. The pseudo-cubic approach associated with the Cu/Zn disordered structure in CZTSSe is described in Fig. S4.

The degrees of order in kesterite CZTSSe samples must be determined for a comprehensive study because the order–disorder transition likely affects TE properties, the details of which are discussed in the SI. Consequently, the order–disorder transition was found to not significantly affect the TE properties of CZTSSe, and all investigated samples had a disordered structure in this study.

TE properties

The temperature-dependent electrical and TE properties of the various compositional CZTSSe single crystals from 300 to 800 K were examined (Fig. 2). Here, σ systematically increased with increasing Se content over the entire temperature range (Fig. 2A). The larger size of the selenide anions than that of the sulfide anions led to a higher degree of covalency and, consequently, a lower band effective mass and higher carrier mobility, resulting in enhanced σ values with increasing Se content. Therefore, the Se alloying and K doping of sample 5 performed to enhance its electrical properties led to a metallic behavior. Results related to the Hall-effect measurements and electrical conduction parameters are listed in Table S3. The conductivity data of samples 1–4 were quantitatively described based on typical thermal activation: where σ_B, E_A, and k_bare are the pre-factor, band-conduction-associated activation energy, and Boltzmann constant, respectively.³⁷ As an intrinsic defect common to Cu-based kesterites, the acceptor level introduced by V_Cu is always shallow and dominant under the Cu-deficient conditions employed in this study. The E_A values of CZTSSe (23–75 meV; Table S3) are consistent with the theoretically predicted shallow V_Cu acceptor level present 20–30 meV above the valence band maximum.³⁸ The Cu/Zn disorder readily leads to a high Cu concentration on the Zn antisite (Cu_Zn) acceptor defect in CZTSSe. The acceptor defect level corresponding to CZTSe (110 meV) is relatively lower than that of CZTS (150 meV) because of the weaker Cu–Se p-d hybridization and its higher valence band,³⁸ which corroborates the enhancement in hole concentration and electrical conduction with increasing Se content at a low thermal activation energy.

Fig. 2 TE characterization of various compositional CZTSSe single crystals. Temperature-dependent (A) electrical conductivity σ [fits were obtained using σ(T) = σ_B exp(−E_A/k_bT)]; (B) Seebeck coefficient S; (C) power factor PF; and (D) PFs of other TE chalcogenide compounds. Measurements were performed in the transverse (⊥) and longitudinal directions (∥) (perpendicular and parallel to the c-axis, respectively).

The S value monotonically increased with increasing temperature because more charge carriers were thermally activated to high-energy states originating from the pseudo-cubic structure (Fig. 2B); this tendency is also observed in TE kesterites.^39–42 The anisotropy of S was caused by the anisotropy of the effective hole mass, indicating higher S values along the a-axis than those along the c-axis.²¹S decreased with increasing Se content and K doping because of enhancements in the hole concentration. The combination of high σ and S values led to a PF value of 2.2 × 10⁻³ W mK⁻² at 800 K for sample 4 along the a-axis (Fig. 2C). A comparison of the PF values of various TE chalcogenides at 300–800 K (Fig. 2D) indicated that the aforementioned PF exceeded that of state-of-the-art chalcogenides, such as SnSe (ZT = 2.6 (ref. 8)), Cu₂Se (ZT = 2.6 (ref. 7)), PbTe (ZT = 2.2 (ref. 15)), AgSbTe₂ (ZT = 1.77 (ref. 30)), PbS (ZT = 1.7 (ref. 18)), Cu_2−xS (ZT = 1.7 (ref. 16)), Cu₁₂Sb₄S₁₃ (ZT = 1.15 (ref. 43)) and polycrystalline CZTS (ZT = 1.1 (ref. 42)).

The temperature-dependent thermal conductivity analysis of the CZTSSe single crystals from 300 to 800 K (Fig. 3A) showed a weak anisotropy of thermal conductivity, where κ⊥ was ∼3% greater than κ∥. Phonon-Boltzmann-transport-based theoretical calculations have indicated that the lattice thermal conductivity along the a-axis in kesterite is ∼5% higher than that along the c-axis.⁴⁴ All investigated samples exhibited an overall T⁻¹ dependence on thermal conductivity, indicating that the lattice thermal conductivity was primarily driven by phonon–phonon Umklapp scattering. Therefore, the CZTSSe system reported herein is suitable for studying the effects of anion substitution on the decrease in lattice thermal conductivity. Here, κ decreased with increasing Se content (Fig. 3A), with the lowest values among the examined specimens (2.56 and 0.93 W mK⁻¹ at 300 and 800 K, respectively) exhibited by Cu₂ZnSnS₂Se₂ (x = 0.5). The Callaway model was adopted to analyze the lattice thermal conductivity of CZTSSe to elucidate the phonon scattering mechanisms, including the mass and strain differences between S and Se.^45,46 Assuming that the Umklapp and point defect scattering were the dominant scattering mechanisms, the ratio of the lattice thermal conductivity of alloyed samples (κ_L) to that of the pure sample (κ_L0) can be expressed using eqn (1) and (2).

where u, θ_D, Ω, h, ν_m, and Γ are the disorder scaling parameter, Debye temperature, average volume per atom, Planck's constant, average sound velocity, and disorder scattering parameter, respectively. Γ is defined as follows:where x, ΔM/M, ε, and Δδ/δ are the Se content, relative change in atomic mass due to Se alloying, lattice anharmonic parameter, and substitution-induced local change in lattice parameters, respectively. Fig. 3B shows a Callaway model fit for thermal conductivity data of CZTSSe specimens with various Se contents at 300 K, in which lattice thermal conductivity is dominant; the parameters used for the calculations are listed in Table S4. The average thermal conductivities of CZTS and CZTSe measured at 300 K and corresponding to both axes were adopted as κ_L0, and the weighted average of both calculations was employed as κ_L. It is worth noting that the parameter ε corresponding to the optimal fit for CZTSSe (ε = 125), which was obtained by implementing mass and strain contrast, was greater than that of PbTe_1−xSe_x (ε = 65–92)⁴⁷ and Cu₂ZnGeSe_4−xS_x (ε = 36).⁴⁸ε is a function of the Grüneisen parameter γ and is consequently a measure of the anharmonic contributions to the lattice. The γ values for CZTSSe at 300 K (2.36–2.37)⁴⁹ are greater than those for PbTe (1.96 in ref. 50 and 0.8 in ref. 51), which is consistent with the tendencies of ε. Experimental evidence of lattice strain was obtained by detecting dislocation loops by transmission electron microscopy (TEM) analysis of sample 2 (x = 0.5; Fig. 3C) assuming high defect concentration related to Se alloying; the dislocation loops were observed in all the samples. Although quantitative analysis of the effects of dislocations on thermal conductivity was difficult, because the dislocation density could not be accurately determined, dislocation loops with diameters of several hundreds of nanometers were observed. A larger lattice strain could possibly lead to a more effective decrease in thermal conductivity.

Fig. 3 Thermal conductivity and lattice distortion in CZTSSe single crystals. (A) Temperature dependence of thermal conductivity κ. The profiles roughly follow a T⁻¹ relationship. (B) Thermal conductivity of CZTSSe at 300 K with various Se contents, featuring a Callaway model fit (dashed line). (C) TEM image of dislocation loops for observation of lattice strain in sample 2 (x = 0.5).

Overall, the outstanding electrical properties and relatively low thermal conductivity of K-doped CZTSSe (sample 4) resulted in a TE material with a ZT of 1.9 at 800 K (Fig. 4A). This value is comparable to the highest ZT reported for other related quaternary polycrystalline materials, because the PF reported herein (2.22 × 10⁻³ W mK⁻²) is 2–18 times greater than those reported in the literature (Table S5). The combined uncertainty of all measurements involved in the calculation of ZT was estimated to be ∼10%. Additionally, the high ZT was proven to be reproducible by conducting repeated measurements for ten thermal cycles (Fig. 4B); the values for each cycle were within 10%. Compared with state-of-the-art chalcogenides,^{5,7,8,15,16,18,30,43} CZTSSe enables the investigated class of materials to be highly efficient for intermediate-temperature waste heat recovery. The approach based on the nature of single crystals in multinary compounds is comparable with existing strategies for discovering low-toxicity earth-abundant TE materials.

Fig. 4 TE performance of CZTSSe single crystals. (A) Temperature-dependent dimensionless figure of merit ZT of CZTSSe single crystals with different Se contents; inset shows sample wafers with (100) and (001) planes. (B) Repeated measurements and a comparison with other TE chalcogenide compounds.

Potassium doping effect on TE properties

Revealing the location of the K atom in CZTSSe (sample 4) is crucial for discussing its doping effects on TE properties. Atomically resolved high-angle annular dark-field scanning TEM (HAADF-STEM) in Fig. S7 and integrated differential phase contrast (iDPC) images in Fig. 5A represent the positions of each element. The locations of Sn, Cu/Zn, and Se/S are identified byZ-contrast, aligning with the expected lattice sites of the pseudo-cubic structure. While revealing the positions of dopant atoms in TEM is challenging, some regions in the iDPC image show distortions around several Sn atoms (circled in yellow), which may suggest that subtle variations in the local electromagnetic environment are caused by impurity atoms, in this case, K.

Fig. 5 The location and electronic effects of K dopants in kesterite materials. (A) iDPC image of the CZTSSe sample from the [100] direction showing the positions of each element. (B) Atom distribution maps of Cu, Zn, Sn and K at the {013} planes. Unfolded effective band structure of the supercell for (C) pristine CZTS and (D) CZTS with the K atom at the interstitial site. ΔCF represents the energy difference between the topmost and the lower valence bands. The color scale reflects the unfolded spectral weight. (E) DOS near the VBM for pristine CZTS and CZTS with the K atom at the interstitial site.

Complementing the TEM result, atom probe tomography (APT) provides further insights into the location of the dopant atoms. The APT data was calibrated using a {013} crystal lattice plane, with the lattice spacing in the spatial distribution map corresponding to the theoretical distance of 0.302 nm (Fig. 5B). Cluster analysis applied to the region of interest in the calibrated data revealed 165 K – Sn clusters, in contrast to 130 K – Cu and 59 K – Zn pairs. This suggests that K tends to cluster with Sn or locate at interstitial sites near the Sn atoms.

The effect of potassium doping on the electronic band structure and transport properties of CZTS was investigated using density functional theory (DFT) calculations. To reproduce the experimentally observed K atom positions near Sn atoms, a K atom was introduced into an interstitial site in a 2 × 2 × 1 supercell (Fig. S8). Fig. 5C and D shows the unfolded effective band structures of the supercells without and with interstitial K atom in CZTS, respectively. The presence of an interstitial K atom induces changes in both the band gap and the fine structure near the valence band maximum (VBM). In pristine CZTS, the VBM is split into two bands by crystal field effects, and the energy difference between the two bands (ΔCF) is 0.058 eV. Upon interstitial K-doping, the reduced symmetry lifts the degeneracy of the second band, resulting in further splitting; however, the energy difference from the topmost band is reduced to 0.045 eV. In addition, the band dispersion along the Γ− (0, 0, 0.5) direction becomes slightly flatter. Consistent with the changes in the band structure, the DOS near the VBM increases in the energy range of approximately 0.05–0.10 eV below the band edge, and the slope of the DOS also becomes steeper (Fig. 5E). These features qualitatively explain the experimental results that K-doping enhances the electrical conductivity while maintaining a high Seebeck coefficient, thereby achieving a remarkably high PF. Similar trends are observed in the band structure of CZTSSe with 50% Se substitution for S upon interstitial K doping (Fig. S8B).

Several studies have reported simultaneous enhancements in hole concentration and carrier mobility induced by alkali–metal incorporation in kesterite compounds; however, the microscopic mechanism responsible for these improvements remains insufficiently understood.^21,25,52,53 A previous study on Na-doped CZTS proposed that Na substitution on the Cu site (Na_Cu) suppresses the dominant donor Zn on the antisite Cu (Zn_Cu) defect, thereby increasing the effective hole concentration and mobility through a reduction in the degree of carrier compensation. In this study, we found that K atoms tend to form clusters with Sn or occupy interstitial sites adjacent to Sn atoms, indicating that K preferentially resides in the local chemical environment around Sn. A recent study on K-doped CZTSe thin films similarly attributed the defect modification to doping-induced homogenization of local composition, which suppresses Sn-related deep defects such as Sn on antisite Sn (Sn_Zn) and Sn vacancies (V_Sn).⁵² Suppression of these defects reduces carrier compensation and ionized impurity scattering, resulting in increased hole concentration and mobility. The presence of K–Sn clusters or K in interstitial sites adjacent to Sn may influence Sn diffusion by increasing local diffusion barriers and stabilizing Sn at its original lattice sites. Nevertheless, the possible role of K in modifying Sn diffusion kinetics remains speculative and represents an important direction for future research.

High TE conversion efficiency in the CZTSSe single-leg

TE devices operating at high temperatures must be thermodynamically stable, with the TE-material–contact interface being particularly important in this regard. Therefore, a Mo(SSe)₂ layer was introduced as a diffusion barrier between CZTSSe and the electrode on the hot side of the TE device because of the presence of mobile atoms such as Cu and K in CZTSSe. Moreover, the Mo(SSe)₂/Mo layer is traditionally used as a back-contact in CZTSSe-based PV devices because of its decent ohmic properties and bonding with CZTSSe.⁵⁴ The TE single-leg device structure of the CZTSSe/Mo(SSe)₂ barrier layer/Mo-contact was applied in this study. The hot-side interface was investigated to confirm the effectiveness of the Mo(SSe)₂ layer as a diffusion barrier. A cross-sectional image of the interfaces in the aforementioned structure was obtained by scanning TEM (STEM) and energy-dispersive X-ray spectroscopy (EDX) (Fig. 6A). Chemical-diffusion-free uniform elemental distributions were observed by EDX-based compositional profiling across the CZTSSe/Mo(SSe)₂/Mo interfaces (Fig. 6B). The K-atom concentration in CZTSSe/Mo(SSe)₂/Mo was <1 at% over the entire measured region. The ∼5-µm-thick Mo(SSe)₂ diffusion barrier had a smooth interface with ∼100-nm-sized roughness undulations. STEM images obtained before and after a 1100 h aging test conducted at 800 K (Fig. S9) revealed that the CZTSSe/Mo(SSe)₂/Mo interfaces hindered chemical diffusion. Although the roughness undulations of the interface locally intensified after the aging test, no significant change of contact resistance was observed in the TE device properties.

Fig. 6 Design of the CZTSSe (x = 0.5) single-leg TE device. (A) and (B) Cross-sectional image of interfaces in the CZTSSe/Mo(SSe)₂-barrier-layer/Mo-contact device obtained by STEM and EDX. (C) Contact resistance at 300 K determined by line scanning using the four-probe technique.

TE devices can be optimized by minimizing the electrical contact resistance to ensure a low internal resistance. The total electrical resistance, including that of the Mo(SSe)₂ diffusion barrier, Mo contact, and Ag adhesive was estimated using a four-probe technique at 300 K to be as low as ∼0.55 mΩ (Fig. 6C), which corresponded to an interfacial contact resistivity of ∼88 µΩ cm², which is higher than the benchmark established by the Bi₂Te₃–Ni interface (1–5 µΩ cm²).⁵⁵

Output voltage (V), output power (P), and conversion efficiency (η) of a CZTSSe single-crystal-based single-leg TE device were monitored as a function of current (I) at different ΔT with a fixed cold-side temperature of 300 K (Fig. 7A–C). Current-dependent output heat flow (Q_out) from the cold side (Fig. S10B) was used to calculate the conversion efficiency, and the details of the calculation are discussed in the SI. The device's theoretical properties were simulated through the COMSOL Multiphysics software with a heat-transfer module using the measured Seebeck coefficient, electrical conductivity, and thermal conductivity in Fig. S11.

Fig. 7 Properties of the CZTSSe (x = 0.5)-based single-leg TE device. (A) Output voltage (V), (B) output power (P), and (C) conversion efficiency (η) as a function of current (I) at various temperature differences (ΔT). (D) Comparison of maximum conversion efficiencies (η_max) of the CZTSSe-based single-leg device and state-of-the-art p-type single-leg equivalents as a function of temperature difference.

The y-intercept and slope of linear V–I plots enable estimation of open-circuit voltage (V_OC) and internal resistance (R_in), respectively. All device parameters that are related to the maximum conversion efficiency are listed in Table 1, wherein all the experimental parameters show good agreement with the simulated results. The increase in V_OC and decrease in R_in with increasing temperature can be explained by the increase in S and σ of CZTSSe. The output power P increased with increasing current, reached a maximum when the external load resistance matched R_in, and subsequently decreased. The maximum output power (P_max) at the hot-side temperature of 773 K (P_max = 163 mW) corresponded to a high output-power-density of 6.52 kW m⁻². An η_max value of 10.4% was achieved by the CZTSSe single-leg device at a temperature gradient of 473 K. It is worth noting that this experimental value of η_max is plotted together with that of other reported single-leg TE devices^{5,26,29–35,43,56,57} as a function of ΔT in Fig. 7D. Although additional device-related developments are required, such as achieving low contact resistivity, these TE systems have remarkable potential for use at intermediate temperatures. The experimentally obtained value of η_max (10.4% at ΔT = 473 K) in this study is slightly lower than that of p-type GeTe (14% at ΔT = 438 K). However, our record value of 10.4% is higher than that of p-type half-Heusler ZrCoBi and other p-type chalcogenide-based single-leg devices such as PbTe, SnSe, AgSbTe₂ and Cu₁₂Sb₄S₁₃. In addition to its high conversion efficiency, the single-leg TE device operated at ΔT = 473 K showed decent thermal stability after various aging tests conducted at 800 K (Fig. S12). The maximum increase in R_in was ∼7% after the 1100 h aging test, and no distinct cracks were observed in the device by optical microscopy (Fig. S13). Such exceptional long-term service reliability at high temperature up to 800 K has rarely been reported in related-chalcogenide-based single-leg devices, marking a substantial advancement in TE device reliability.

Table 1 Device parameters of the CZTSSe single-leg from experimental and theoretical results. Open-circuit voltage (V_OC), the internal resistance (R_in), maximum output power (P_max), and maximum conversion efficiency (η_max) of the CZTSSe (sample 4) single-leg under different temperature gradients (ΔT) with a fixed cold-side temperature of 300 K

	ΔT (K)	V OC (mV)	R in (mΩ)	P max (mW)	η max (%)
Measured values	73	16.6	40.1	1.7	0.4
	173	42.3	38.9	11.6	1.4
	273	71.7	35.3	37.2	3.7
	373	93.2	32.1	66.7	5.2
	473	145.1	31.8	163.3	10.4
Simulated values	73	20.2	31.7	3.2	0.8
	173	48.4	28.7	20.3	2.6
	273	77.0	26.8	55.0	5.5
	373	108.9	25.3	116.3	9.0
	473	150.6	23.9	233.3	14.3

---

Conclusions

This study was aimed at achieving competitive ZT values for kesterite CZTSSe via a synergistic strategy involving a pseudo-cubic approach with a disordered structure and alloying. A record-high ZT of 1.9 was achieved at 800 K using a K-doped Cu_1.9ZnSn(S_1−xSe_x)₄ single crystal (x = 0.5), whose PF (2.22 × 10⁻³ W mK⁻²) can compete with that of state-of-the-art TE materials. Importantly, a high conversion efficiency of 10.4% was achieved using a CZTSSe-based single-leg TE device with long-term stability. The following conclusions were drawn from this study: (i) interstitial K-doping (alkali effect) in CZTSSe enhances the electrical conductivity while maintaining a high Seebeck coefficient due to the steeper slope of DOS; (ii) the CZTSSe/Mo(SSe)₂-barrier-layer/Mo-contact device structure led to decent ohmic properties and impeded chemical diffusion; (iii) thermal stability with small fluctuations of device parameters was ensured at aging times of up to 1100 h at 800 K. The results can pave the way for intermediate-temperature TE harvesting and accelerate the development of TE technology.

Materials and methods

Single-crystal growth and fabrication of a single-leg TE device

CZTSSe single crystals were grown using the Sn-solvent-based traveling heater method (THM at an axial temperature gradient of 40 °C cm⁻¹, a growth temperature of 850 °C, and a growth rate of 4 mm day⁻¹). Details are provided in the SI. A diamond blade was used to dice the CZTSSe (x = 0.5) single crystals to a leg size of approximately 5 mm × 5 mm × 10 mm. The grown crystals and legs were mechanically polished with 0.01 µm Al₂O₃ powder and subsequently etched with HCl solution for 5 min to remove the saw and polishing damage. First, an ∼100-µm-thick Mo-contact layer was deposited on both the top and bottom sides of the leg by RF sputtering at 300 K in an Ar atmosphere. The Mo(SSe)₂ barrier layer was formed by high-temperature annealing at 700 °C for 10 h under high vacuum (∼10⁻⁴ Pa), which enabled the Mo contact to react with the S and Se in CZTSSe. Subsequently, both sides and the Cu electrode were attached using Ag paste by heat treatment at 250 °C for 1 h in a N₂ atmosphere.

Characterization

The structural properties were analyzed by powder X-ray diffraction (XRD; Panalytical X'Pert PRO) and Raman spectroscopy (HORIBA T64000). The XRD measurements were performed at 40 kV and 40 mA using a Cu Kα radiation source. For the Raman analysis, a 532 nm laser excitation source was focused on the sample via an objective lens with a numerical aperture of 0.55, resulting in a laser power of 100 mW on the sample. The spectra were calibrated based on the 520 cm⁻¹ Si peak. The orientation and crystallinity of the CZTSSe single crystals were investigated by X-ray rocking curve (XRC) analysis (Panalytical X'Pert PRO). The XRC measurements were performed at 45 kV and 40 mA with a step width of 0.002° using a Cu Kα radiation source. The crystallographic orientation of all samples can be observed using a scanning electron microscope (SEM) equipped with an electron backscatter diffraction (EBSD; Quanta 600) system shown in Fig. S3. The samples were imaged at 30 kV using an SEM probe and that used an off-axis backscatter detector.

The detailed composition of the crystals was determined by inductively coupled plasma-atomic emission spectroscopy (ICP-AES; SII NanoTechnology SPS3520UV). Samples were dissolved in a mixture of HNO₃ and HCl to facilitate this analysis.

Microstructure and composition of the CZTSSe single crystal and single-leg were determined by TEM and STEM (JEOL 2800) at 200 kV using two simultaneous solid-state EDX detectors. Spectral image acquisition was performed over a series of consecutive subsecond frames with drift correction between frames and a total acquisition time of less than 20 min. The EDX data were processed using Thermo Scientific Image Analysis software. The analytical certainty associated with EDX profiling was within 0.5 at%.

Atomic structure investigations were carried out by using TEM and APT. HAADF-STEM and iDPC images were acquired using a Themis-Z S/TEM (Thermo Fisher Scientific), operated at 300 kV, equipped with double-aberration correctors. The beam convergence angle in both HAADF and iDPC modes was 17.9 mrad, and the collection angles in HAADF and iDPC modes were 48–200 and 9–35 mrad, respectively. The iDPC images were acquired using a quadrant segmented dark field detector. APT was conducted using a CAMECA Local Electrode Atom Probe (LEAP) 4000X Si equipped with a picosecond-pulse ultraviolet laser. A laser pulse energy of 50 pJ and laser frequency of 200 kHz were used. APT was performed at ∼50 K and the data reconstruction and analysis were performed via CAMECA's AP Suite 6.3 software. The reconstruction was calibrated using the method described by Gault et al.^58,59 Site-specific TEM and APT specimens were prepared by standard lift-out techniques using a Helios G4 plasma focused ion beam scanning electron microscope (Thermo Fisher Scientific).

The electrical conductivity σ and Seebeck coefficient S of samples with dimensions of approximately 2 mm × 2 mm × 10 mm were simultaneously measured in a helium atmosphere at 300–800 K using a commercial system (ZEM3; ADVANCE RIKO). The thermal conductivity κ was calculated using the equation κ = λC_pD, where λ is the thermal diffusivity, C_p is the specific heat capacity, and D is the density. The thermal diffusivity of rectangular samples (8 mm × 8 mm × 1 mm) was measured using the laser flash method (Netzsch; LFA467); the samples were coated with a thin layer of graphite to minimize errors due to the emissivity of the material. The specific heat capacity of samples measuring 3 mm × 3 mm × 3 mm was determined by differential scanning calorimetry (DSC; Rigaku Thermo plus EVO2 DSCvesta). The mass of all samples was measured by using thermogravimetric and differential thermal analysis (TG-DTA; Rigaku Thermo plus EVO2 TG-DTA8122). Densities in the 4.4–4.5 g cm⁻³ range were measured using the Archimedes method at room temperature in ultrapure water.

The hole concentration p and mobility µ were determined based on the temperature dependence of the Hall effect (TOYO; ResiTest8400) in a 0.45 T magnetic field with the van der Pauw geometry in a He atmosphere at 300 K. For the Hall measurements, 200–300-nm-thick Au contacts with a diameter of 1 mm were evaporated onto the corners of each CZTSSe bulk single crystal with dimensions of 5 mm × 5 mm × 0.8 mm.

A four-probe technique was used to measure the electrical resistance of the single-leg CZTSSe device. A movable probe measured the voltage along the sample length when an alternating current was applied across the leg. The conversion efficiency of the CZTSSe-based single-leg TE device was evaluated by measuring electrical output power and output heat flow under vacuum (∼10⁻² Pa) using an M612-B tester (MOTTAINAI ENERGY).

Computational details

The DFT calculations were performed using the Quantum ESPRESSO (QE) package^60,61 with projector augmented wave pseudopotentials.⁶² The generalized gradient approximation functional with Perdew–Burke–Ernzerhof parametrization (GGA-PBE)⁶³ was chosen for the exchange and correlation potentials. We applied the Hubbard U correction to the strongly correlated Cu 3d and Zn 3d orbitals. We chose U_eff (= U − V) values as 5.2 eV and 6.5 eV for Cu and Zn, respectively. Here, U, V, and U_eff are the on-site Coulomb interaction strength, the on-site exchange parameter, and the resulting effective repulsion applied to the localized d states, respectively. Spin–orbit coupling was not considered in the present study. The conventional unit cell with a tetragonal lattice containing 16 atoms was fully relaxed until the residual forces became less than 10⁻⁶ eV Å⁻¹. A cutoff energy of 130 Ry for the plane wave basis and an 8 × 8 × 4 k-grid were used for structure relaxation. The 2 × 2 × 1 supercell was constructed based on the relaxed unit cell. Then, one K atom was introduced at an interstitial site identified using the pymatgen-analysis-defect module.⁶⁴ Although there are two candidate interstitial sites suggested by Voronoi tessellation, the one that reproduces the experimental atomic position obtained from STEM observation was adopted. To model the alloying in the (S,Se)-site in CZTSSe, the special quasirandom structure (SQS) method implemented in the Alloy Theoretic Automated Toolkit (ATAT) was used.⁶⁵ The experimentally determined lattice parameters were used and atomic positions in the supercell were fully relaxed until the residual forces became less than 10⁻³ eV Å⁻¹. A cutoff energy of 130 Ry for the plane wave basis and a 4 × 4 × 4 k-grid were used for supercell relaxation. The effective band structures of the supercells in the primitive Brillouin zone were obtained by using the band unfolding method implemented in BandUP code.⁶⁶

Author contributions

All authors discussed the results and provided feedback on the manuscript. Conceptualization: A. Nagaoka; methodology: A. Nagaoka, S. Miura, K. Nomoto, and N. Sato; investigation: A. Nagaoka, S. Miura, K. Nomoto, K. Chen, and N. Sato; visualization: A. Nagaoka, S. Miura, K. Nomoto, K. Chen, and N. Sato; supervision: K. Yoshino and K. Nishioka; writing—original draft: A. Nagaoka; writing—review & editing: A. Nagaoka, K. Nomoto, N. Sato, and K. Nishioka.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information: Single crystal growth by Traveling Heater Method (THM). Characterizations for structural properties, pseudo-cubic structure, degree of order in CZTSSe, compositional properties, and thermoelectric conversion efficiency. Fig. S1 to S13; Tables S1 to S5; references. See DOI: https://doi.org/10.1039/d5ta08640j.

Acknowledgements

This study was supported by JSPS KAKENHI Grant Number JP23K26603 (A. Nagaoka), JSPS KAKENHI Grant Number JP22K14505 (N. Sato), TEPCO Memorial Foundation No. 23-006 (A. Nagaoka), and the Australian Research Council Discovery Early Career Researcher Award (DE220100527) (K. Nomoto).

References

1. C. Forman, I. K. Muritala, R. Pardemann and B. Meyer, Renewable Sustainable Energy Rev., 2016, 57, 1568–1579.
2. A. Firth, B. Zhang and A. Yang, Appl. Energy, 2019, 235, 1314–1334.
3. Y. Sun, S. Tian, P. Ciais, Z. Zeng, J. Meng and Z. Zhang, Nat. Commun., 2022, 13, 297.
4. X. L. Shi, J. Zou and Z. G. Chen, Chem. Rev., 2020, 120, 7399–7515.
5. W. He, D. Wang, H. Wu, Y. Xiao, Y. Zhang, D. He, Y. Feng, Y. J. Hao, J. F. Dong, R. Chetty, L. Hao, D. Chen, J. Qin, Q. Yang, X. Li, J. M. Song, Y. Zhu, W. Xu, C. Niu, X. Li, G. Wang, C. Liu, M. Ohta, S. J. Pennycook, J. He, J. F. Li and L. D. Zhao, Science, 2019, 365, 1418–1424.
6. Mineral Commodity Summaries 2024, United States Geological Survey.
7. A. A. Olvera, N. A. Moroz, P. Sahoo, P. Ren, T. P. Bailey, A. A. Page, C. Uher and P. F. P. Poudeu, Energy Environ. Sci., 2017, 10, 1668–1676.
8. L. D. Zhao, S. H. Lo, Y. Zhang, H. Sun, G. Tan, C. Uher, C. Wolverton, V. P. Dravid and M. G. Kanatzidis, Nature, 2014, 508, 373–377.
9. J. Sui, J. Li, J. He, Y.-L. Pei, D. Berardan, H. Wu, N. Dragoe, W. Cai and L.-D. Zhao, Energy Environ. Sci., 2013, 6, 2916–2920.
10. C. Fu, S. Bai, Y. Liu, Y. Tang, L. Chen, X. Zhao and T. Zhu, Nat. Commun., 2015, 6, 8144.
11. Y. Luo, J. Yang, Q. Jiang, W. Li, D. Zhang, Z. Zhou, Y. Cheng, Y. Ren and X. He, Adv. Energy Mater., 2016, 6, 1600007.
12. J. Zhang, L. Song, S. H. Pedersen, H. Yin, L. T. Hung and B. B. Iversen, Nat. Commun., 2017, 8, 13901.
13. K. F. Hsu, S. Loo, F. Guo, W. Chen, J. S. Dyck, C. Uher, T. Hogan, E. K. Polychroniadis and M. G. Kanatzidis, Science, 2004, 303, 818–821.
14. X. Shi, J. Yang, J. R. Salvador, M. Chi, J. Y. Cho, H. Wang, S. Bai, J. Yang, W. Zhang and L. Chen, J. Am. Chem. Soc., 2011, 133, 7837–7846.
15. K. Biswas, J. He, I. D. Blum, C. I. Wu, T. P. Hogan, D. N. Seidman, V. P. Dravid and M. G. Kanatzidis, Nature, 2012, 489, 414–418.
16. Y. He, T. Day, T. Zhang, H. Liu, X. Shi, L. Chen and G. J. Snyder, Adv. Mater., 2014, 26, 3974–3978.
17. S. I. Kim, K. H. Lee, H. A. Mun, H. S. Kim, S. W. Hwang, J. W. Roh, D. J. Yang, W. H. Shin, X. S. Li, Y. H. Lee, G. J. Snyder and S. W. Kim, Science, 2015, 348, 109–114.
18. M. Ibanez, Z. Luo, A. Genc, L. Piveteau, S. Ortega, D. Cadavid, O. Dobrozhan, Y. Liu, M. Nachtegaal, M. Zebarjadi, J. Arbiol, M. V. Kovalenko and A. Cabot, Nat. Commun., 2016, 7, 10766.
19. J. Zhang, R. Liu, N. Cheng, Y. Zhang, J. Yang, C. Uher, X. Shi, L. Chen and W. Zhang, Adv. Mater., 2014, 26, 3848–3853.
20. A. Nagaoka, T. Masuda, S. Yasui, T. Taniyama and Y. Nose, Appl. Phys. Express, 2018, 11, 051203.
21. A. Nagaoka, K. Yoshino, T. Masuda, T. D. Sparks, M. A. Scarpulla and K. Nishioka, J. Mater. Chem. A, 2021, 9, 15595–15604.
22. P. M. P. Salomé, H. Rodriguez-Alvarez and S. Sadewasser, Sol. Energy Mater. Sol. Cells, 2015, 143, 9–20.
23. D. Sharma, N. Nicoara, P. Jackson, W. Witte, D. Hariskos and S. Sadewasser, Nat. Energy, 2024, 9, 163–171.
24. P. Jood, J. P. Male, S. Anand, Y. Matsushita, Y. Takagiwa, M. G. Kanatzidis, G. J. Snyder and M. Ohta, J. Am. Chem. Soc., 2020, 142, 15464–15475.
25. A. Nagaoka, H. Miyake, T. Taniyama, K. Kakimoto, Y. Nose, M. A. Scarpulla and K. Yoshino, Appl. Phys. Lett., 2014, 104, 152101.
26. Z. Bu, X. Zhang, B. Shan, J. Tang, H. Liu, Z. Chen, S. Lin, W. Li and Y. Pei, Sci. Adv., 2021, 7, eabf2738.
27. X. Wu, Y. Lin, Z. Han, H. Li, C. Liu, Y. Wang, P. Zhang, K. Zhu, F. Jiang, J. Huang, H. Fan, F. Cheng, B. Ge and W. Liu, Adv. Energy Mater., 2022, 12, 2203039.
28. A. Nagaoka, K. Nagatomo, K. Nakashima, Y. Hirai, Y. Ota, K. Yoshino and K. Nishioka, Mater. Trans., 2023, 64, 2535–2541.
29. X. L. Shi, W. D. Liu, M. Li, Q. Sun, S. D. Xu, D. Du, J. Zou and Z. G. Chen, Adv. Energy Mater., 2022, 12, 2200670.
30. B. C. Chen, K. K. Wang and H. J. Wu, Small, 2024, 2401723.
31. Z. Liu, J. Mao, J. Sui and Z. Ren, Energy Environ. Sci., 2018, 11, 23–44.
32. H. Zhu, R. He, J. Mao, Q. Zhu, C. Li, J. Sun, W. Ren, Y. Wang, Z. Liu, Z. Tang, A. Sotnikov, Z. Wang, D. Broido, D. J. Singh, G. Chen, K. Nielsch and Z. Ren, Nat. Commun., 2018, 9, 2497.
33. S. Liu, Y. Wen, S. Bai, H. Shi, Y. Qin, B. Qin, D. Liu, Q. Cao, X. Gao, L. Su, C. Chang, X. Zhang and L. D. Zhao, Adv. Mater., 2024, e2401828.
34. P. Sauerschnig, P. Jood and M. Ohta, Adv. Mater. Technol., 2023, 8, 2201295.
35. X. Y. Tan, J. Dong, J. Liu, D. Zhang, S. F. D. Solco, K. Saglik, N. Jia, I. You, S. W. Chien, X. Wang, L. Hu, Y. Luo, Y. Zheng, D. X. Y. Soo, R. Ji, K. C. H. Goh, Y. Jiang, J. F. Li, A. Suwardi, Q. Zhu, J. Xu and Q. Yan, Adv. Sci., 2024, e2400870.
36. A. Nagaoka, R. Katsube, S. Nakatsuka, K. Yoshino, T. Taniyama, H. Miyake, K. Kakimoto, M. A. Scarpulla and Y. Nose, J. Cryst. Growth, 2015, 423, 9–15.
37. A. Nagaoka, H. Miyake, T. Taniyama, K. Kakimoto and K. Yoshino, Appl. Phys. Lett., 2013, 103, 112107.
38. S. Chen, A. Walsh, X. G. Gong and S. H. Wei, Adv. Mater., 2013, 25, 1522–1539.
39. M.-L. Liu, F.-Q. Huang, L.-D. Chen and I. W. Chen, Appl. Phys. Lett., 2009, 94, 202103.
40. Y. Dong, H. Wang and G. S. Nolas, Inorg. Chem., 2013, 52, 14364–14367.
41. Y. Dong, H. Wang and G. S. Nolas, Phys. Status Solidi, 2014, 8, 61–64.
42. Q. Jiang, H. Yan, Y. Lin, Y. Shen, J. Yang and M. J. Reece, J. Mater. Chem. A, 2020, 8, 10909–10916.
43. H. Hu, H. L. Zhuang, Y. Jiang, J. Shi, J. W. Li, B. Cai, Z. Han, J. Pei, B. Su, Z. H. Ge, B. P. Zhang and J. F. Li, Adv. Mater., 2021, 33, e2103633.
44. T. Gürel, C. Sevik and T. Çağın, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 205201.
45. J. Callaway and H. C. von Baeyer, Phys. Rev., 1960, 120, 1149–1154.
46. B. Abeles, Phys. Rev., 1963, 131, 1906–1911.
47. H. Wang, A. D. LaLonde, Y. Pei and G. J. Snyder, Adv. Funct. Mater., 2013, 23, 1586–1596.
48. C. P. Heinrich, T. W. Day, W. G. Zeier, G. J. Snyder and W. Tremel, J. Am. Chem. Soc., 2014, 136, 442–448.
49. X. He, J. Pi, Y. Dai and X. Li, Acta Metall. Sin., 2013, 26, 285–292.
50. S.-H. Lo, J. He, K. Biswas, M. G. Kanatzidis and V. P. Dravid, Adv. Funct. Mater., 2012, 22, 5175–5184.
51. W. G. Zeier, A. LaLonde, Z. M. Gibbs, C. P. Heinrich, M. Panthofer, G. J. Snyder and W. Tremel, J. Am. Chem. Soc., 2012, 134, 7147–7154.
52. S. Tao, L. Dong, J. Han, Y. Wang, Q. Gong, J. Wei, M. Zhao and D. Zhuang, J. Mater. Chem. A, 2023, 11, 9085–9096.
53. L. Wang, J. Ban, L. Han, Z. Zhou, W. Zhou, D. Kou, Y. Meng, Y. Qi, S. Yuan and S. Wu, J. Mater. Chem. A, 2024, 12, 25643–25677.
54. W. Wang, M. T. Winkler, O. Gunawan, T. Gokmen, T. K. Todorov, Y. Zhu and D. B. Mitzi, Adv. Energy Mater., 2014, 4, 1301465.
55. S. P. Feng, Y. H. Chang, J. Yang, B. Poudel, B. Yu, Z. Ren and G. Chen, Phys. Chem. Chem. Phys., 2013, 15, 6757–6762.
56. R. Chetty, Y. Kikuchi, Y. Bouyrie, P. Jood, A. Yamamoto, K. Suekuni and M. Ohta, J. Mater. Chem. C, 2019, 7, 5184–5192.
57. J. Jiang, H. Zhu, Y. Niu, Q. Zhu, S. Song, T. Zhou, C. Wang and Z. Ren, J. Mater. Chem. A, 2020, 8, 4790–4799.
58. B. Gault, F. de Geuser, L. T. Stephenson, M. P. Moody, B. C. Muddle and S. P. Ringer, Microsc. Microanal., 2008, 14, 296–305.
59. B. Gault, M. P. Moody, F. de Geuser, G. Tsafnat, A. La Fontaine, L. T. Stephenson, D. Haley and S. P. Ringer, J. Appl. Phys., 2009, 105, 034913.
60. P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari and R. M. Wentzcovitch, J. Phys.:Condens. Matter, 2009, 21, 395502.
61. P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. Buongiorno Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. Dal Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio Jr, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H. Y. Ko, A. Kokalj, E. Kucukbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H. V. Nguyen, A. Otero-de-la-Roza, L. Paulatto, S. Ponce, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu and S. Baroni, J. Phys.:Condens. Matter, 2017, 29, 465901.
62. P. E. Blochl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
63. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
64. S. P. Ong, W. D. Richards, A. Jain, G. Hautier, M. Kocher, S. Cholia, D. Gunter, V. L. Chevrier, K. A. Persson and G. Ceder, Comput. Mater. Sci., 2013, 68, 314–319.
65. A. van de Walle, Calphad, 2009, 33, 266–278.
66. P. V. C. Medeiros, S. Stafström and J. Björk, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 89, 041407.

---