Simultaneous enhancement of power factor and suppression of thermal conductivity in bulk TlFe_1.6Se₂ via embedded atomically thin FeSe layers

Abstract

FeSe in the monolayer limit exhibits extremely large thermoelectric power factors (PFs). Extending the high-PF concept from two-dimensional FeSe to bulk materials, together with the suppression of lattice thermal conductivity, enables higher-performance thermoelectrics. Here, layered TlFe_1.6Se₂ is identified as a model system consisting of atomically thin two-dimensional FeSe layers separated by Tl atoms; i.e., FeSe monolayers are naturally confined within a bulk crystal. This compound uniquely exhibits a transition from Fe-vacancy (V_Fe)-ordered to-disordered states at around 200 °C. Although the V_Fe-disordered phase exhibits high electrical conductivity, carrier compensation suppresses the Seebeck coefficient and limits PF. In contrast, the V_Fe-ordered phase shows an enhanced Seebeck coefficient associated with Mott gap formation, resulting in improved PF which is much higher than that of bulk FeSe. The lattice thermal conductivity of the V_Fe-ordered phase is lower than that of representative thermoelectric chalcogenides and that of the V_Fe-disordered phase further decreases to ∼0.2 W (m⁻¹ K⁻¹) at 500 °C due to the V_Fe-induced bond heterogeneity. Consequently, the dimensionless figure of merit (ZT) of TlFe_1.6Se₂ reaches ∼0.2 at 50 °C in the V_Fe-ordered phase, which is two orders of magnitude higher than that of bulk FeSe. These results demonstrate that confining FeSe monolayers within a bulk crystal, alongside vacancy order–disorder control, is an effective design strategy for next-generation thermoelectrics.

---

1 Introduction

Thermoelectric energy conversion is considered a promising renewable technology for converting waste heat into electricity via the Seebeck effect.^1–3 The energy-conversion efficiency of thermoelectric materials is governed by the dimensionless figure of merit (ZT), defined as ZT = S²σT/κ, where S is the Seebeck coefficient, σ is the electrical conductivity, T is the absolute temperature, and κ is the thermal conductivity. The electrical output power is given by S²σ, known as the power factor (PF), and κ includes the electronic thermal conductivity (κ_ele) and lattice thermal conductivity (κ_lat). Achieving a high ZT requires simultaneously attaining a large S and high σ while minimizing κ. However, this requirement is fundamentally constrained by well-known trade-offs among these parameters. For example, σ and S exhibit an inverse dependence on the carrier concentration (n): increasing n enhances σ but reduces S, resulting in a maximum PF at an optimal n. Moreover, increasing σ typically raises κ_ele according to the Wiedemann–Franz law (κ_ele = LσT, where L is the Lorenz number). Owing to this intricate interdependence of transport properties, realizing thermoelectric materials with high ZT remains a major challenge.⁴ In contrast, κ_lat is largely independent of these electronic transport parameters, and materials with intrinsically low κ_lat are therefore well suited for achieving high ZT.

Recently, a giant PF has been reported in ultra-thin films of layered FeSe,⁵ one of the high-critical-temperature (T_c) Fe-based superconductors.⁶ Ultra-thin FeSe film grown on the SrTiO₃ substrate has attracted considerable research attention because it exhibits superconductivity with T_c up to 100 K,⁷ dramatically higher than the bulk T_c of ∼9 K. When the film thickness is reduced to ∼1 nm, σ increases due to interface-induced carrier doping from the SrTiO₃ substrate, while S also increases due to the enhanced two-dimensionality and electron correlation effects. As a result, the PF reaches 260 µW (cm⁻¹ K⁻²), a value much higher than that of practical Bi₂Te₃ at room temperature (RT).⁵ Furthermore, the PF of ultra-thin FeSe increases to 13 000 µW (cm⁻¹ K⁻²) at a reduced temperature of 50 K. Assuming κ ≈ 5 W (m⁻¹ K⁻¹), similar to bulk Fe-based superconductors,^8,9 the ZT of the ultra-thin FeSe is estimated to be ∼1.5 at 280 K. However, the excellent PF and ZT values appear only in the ultra-thin limit, which makes these films unsuitable for thermoelectric power generation. It would therefore be more interesting to achieve such a high PF throughout the bulk phase of FeSe. In addition, FeSe has been reported to exhibit a relatively high κ_lat ≈ 4.5 W (m⁻¹ K⁻¹) at RT in polycrystalline samples,¹⁰ which is nearly three times higher than that of Bi₂Te₃.¹¹ Therefore, reducing the κ_lat in FeSe is also crucial for realizing its potential as a high-performance thermoelectric material.

To overcome the limitations of ultra-thin FeSe, we herein focus on a layered A_1−xFe_2−ySe₂ (A = alkali metals and Tl), which can be regarded as a natural bulk analogue of FeSe monolayers. This compound crystallizes in a tetragonal ThCr₂Si₂-type structure composed of alternately stacked FeSe layers and A-ion layers along the c-axis (Fig. 1a).¹² This crystal structure intrinsically contains atomically thin FeSe monolayers within a bulk crystal. The FeSe layers in A_1−xFe_2−ySe₂ possess an identical crystal lattice to that of FeSe, in which Fe atoms form a square lattice and Se atoms occupy apical positions of the edge-shared FeSe₄ tetrahedra. The thickness of the FeSe layer in TlFe_1.6Se₂ is ∼3.0 Å, corresponding to a single FeSe slab, which is significantly thinner than the ∼1 nm minimum thickness reported for ultrathin FeSe films.⁵ Notably, the Fe sites intrinsically contain vacancies (V_Fe), forming a √5 × √5 × 1 supercell of the ThCr₂Si₂-type structure (Fig. 1b). Upon heating, the arrangement of V_Fe undergoes a transition between the ordered and disordered states at ∼500 K. Theoretical calculations have suggested that the V_Fe-ordered A_1−xFe_2−ySe₂ is a Mott insulator with a Mott gap of ∼100 meV,^13,14 and the gap has been experimentally confirmed to be ∼430 meV.¹⁵ The V_Fe-ordered Mott insulator phase also exhibits long-range antiferromagnetic (AFM) order with a Néel temperature as high as 470−560 K.¹² In contrast, by tuning deficiencies at the A site and the Fe site, A_1−xFe_2−ySe₂ can exhibit superconductivity with T_c = 30 K.^16,17 Note that the electronic transport in A_1−xFe_2−ySe₂ is mainly governed by FeSe layers, where the electronic states near the Fermi level are dominated by Fe 3d orbitals.¹⁴ As a result, the contribution of A-site ions to the electron transport is relatively limited. Based on these characteristics, A_1−xFe_2−ySe₂ is expected to exhibit a high PF owing to its two-dimensional FeSe layers naturally embedded within the bulk crystal. In particular, strong electron correlation may lead to an enhanced S, while the presence of the ordered and disordered V_Fe defects is expected to promote phonon scattering and thereby reduce κ_lat.

Fig. 1 Crystal structure of the layered AFe_1.6Se₂ (A = alkali metals and Tl) with an Fe-vacancy (V_Fe) ordered structure, viewed along the (a) [120] and (b) [001] directions. The spheres represent the Tl (gray), Fe (red), Se (blue) atoms, and V_Fe sites (green). The structure consists of FeSe layers, crystallographically equivalent to the high-power-factor (PF) FeSe monolayers, sandwiched by the A-ion layers. The Fe sites intrinsically host V_Fe, which undergo an order–disorder transition upon heating to ∼500 K, which is expected to enhance phonon scattering. The coexistence of FeSe layers with intrinsic V_Fe defects is expected to enable a high PF, together with low lattice thermal conductivity.

In this study, we investigate the electronic and phonon transport properties of bulk polycrystalline layered TlFe_1.6Se₂ as a model system to clarify the effectiveness of embedded two-dimensional FeSe layers in enhancing thermoelectric performance. This compound is much more stable in air than the other A_1−xFe_2−ySe₂ phases (A = K, Rb, and Cs),¹⁸ and a fully V_Fe-ordered phase with high chemical homogeneity can be obtained due to the lower vapor pressure of Tl compared with those of alkali metals.¹⁹ We demonstrate that the PF of TlFe_1.6Se₂ is significantly higher than bulk FeSe; furthermore, TlFe_1.6Se₂ exhibits an exceptionally low κ_lat, with both the electronic and phonon transport properties strongly modulated by the V_Fe ordering and disordering, resulting in an enhanced ZT.

2 Synthesis and characterization

2.1 Synthesis

FeSe was synthesized by a solid-state reaction using Fe powder (purity: 99.999%, Kojundo Chemical Lab.) and Se powder (purity: 99.999%, Kojundo Chemical Lab.), mixed at the molar ratio of Fe : Se = 1 : 1. The mixed powders were sealed in an evacuated silica-glass ampule and annealed at 400 °C for 5 h, followed by annealing at 800 °C for 10 h. TlFe_1.6Se₂ was prepared by a solid-state reaction of fine pieces of the Tl metal (purity: 99.9%, Aldrich) and powders of FeSe and Se. The precursors were mixed in a stoichiometric atomic ratio of Tl : FeSe : Se = 1 : 1.6 : 0.4 and sealed in an Ar-filled stainless-steel tube. The mixture was reacted at 400 °C for 5 h and subsequently at 650 °C for 10 h.²⁰ The resulting FeSe and TlFe_1.6Se₂ powders were thoroughly ground and sintered into pellets by a hot-press method at 500 °C for 10 min. All reagents and products were handled in a glove box filled with dry Ar gas (dew point < −80 °C and oxygen concentration < 1 ppm).

2.2 Characterization

Crystalline phases were examined by X-ray diffraction (XRD) in the Bragg–Brentano geometry using a Cu Kα radiation source (MiniFlex600, Rigaku Co.). The lattice parameters were refined by the Pawley method using the TOPAS ver. 4.2 program (Bruker AXS GmbH). The microstructure and chemical compositions (i.e., atomic ratios of Tl, Fe, and Se) of the samples were characterized using an electron probe micro-analyzer (EPMA).

σ and S were simultaneously measured using a ZEM-3 (ULVAC-RIKO, Inc.) under a He atmosphere. S was obtained from the thermoelectric voltage generated under an applied temperature difference (ΔT); the thermoelectromotive force (ΔV) and ΔT were measured simultaneously, and S was determined from the slope of the ΔV–ΔT plots. Pt electrodes deposited by sputtering were used to form ohmic contacts. κ was obtained from κ = D × C × d, where the thermal diffusivity (D) was measured in an Ar atmosphere by the laser flash method (LFA 457, NETZSCH), the specific heat capacity (C) was measured by differential scanning calorimetry (DSC), and the sample density (d) was determined from the mass and geometric dimensions of the pellets. κ_ele was calculated using κ_ele = LTσ. L was calculated based on a single-parabolic-band model as follows:²¹

where k_B is the Boltzmann constant. Here, the reduced Fermi energy (η) was obtained based on the free carrier model using the measured S as with the Fermi integral defined as where r = −1/2 is the scattering factor.²² κ_lat was obtained by subtracting κ_ele from the total κ, i.e., κ_lat = κ − κ_ele. σ and S were measured in the in-plane direction of the hot-pressed pellets, while κ was measured along the out-of-plane direction.

2.3 Theoretical analysis

First-principles calculations were conducted using the projector-augmented wave (PAW) method as implemented in the Vienna Ab initio Simulation Package (VASP).^23,24 Tl [5d 6s 6p], Fe [3d 4s], and Se [4s 4p] orbitals were treated as valence states. Variable-cell structural optimizations were performed using the Perdew–Burke–Ernzerhof functional adapted for solids (PBEsol) functional,²⁵ with a plane-wave cutoff energy of 550 eV, the convergence criteria of 10⁻⁶ eV for the energy and 0.01 eV Å⁻¹ for the force, and a Γ-centered k-mesh with a k-spacing of 0.1 Å⁻¹. Phonon calculations were performed using the ALAMODE codes.^26,27 For FeSe, a 4 × 4 × 2 supercell containing 192 atoms was employed to compute the interatomic force constants (IFCs). For TlFe_1.6Se₂, a 2 × 2 × 2 supercell with 368 atoms was used to obtain the harmonic IFCs, whereas a 1 × 1 × 1 unit-cell with 46 atoms was used to evaluate the anharmonic IFCs. The harmonic IFCs were fixed to the values determined by the finite-displacement approach,^28,29 and the anharmonic IFCs up to the sixth order were estimated by the compressive sensing lattice dynamics. We included all allowed interactions for the harmonic IFCs, all third-order IFCs within a 12-bohr cutoff radius, and all fourth-to sixth-order IFCs within an 8-bohr cutoff radius. Density functional theory (DFT) calculations to obtain energies and forces were performed using the PBEsol functional with a plane-wave energy cutoff of 400 eV and an energy convergence criterion of 10⁻⁸ eV. Finite-temperature phonon dispersions were computed by using the self-consistent phonon (SCPH) approach.³⁰ κ_lat was determined by considering the Peierls term, obtained by solving the Peierls–Boltzmann transport equation (PBTE) under the single-mode relaxation time approximation, and the coherent term using the unified formula.³¹ The phonon-isotope scatterings were considered.³² The convergence of κ_lat with respect to the q-point mesh was verified using a 7 × 7 × 7 grid for TlFe_1.6Se₂ and a 15 × 15 × 10 grid for FeSe in the PBTE calculations.

3 Results and discussion

3.1 Electrical transport properties

Crystal structures and microstructure characterizations are summarized in Fig. 2, S1 and S2 in the SI. It was confirmed that TlFe_1.6Se₂ adopted the V_Fe-ordered crystal structure with a space group of I4/m (Fig. 2a) and FeSe crystallized in a PbO-type α-phase with a space group of P4/nmm at RT (Fig. S2). Upon heating the TlFe_1.6Se₂ bulk sample, a transition from the V_Fe-ordered structure to the V_Fe-disordered structure (space group: I4/mmm) was observed in the T range of 100 °C–200 °C (Fig. 2b). At higher T, the V_Fe-disordered phase became stable.³³ When cooled back to RT, the V_Fe-ordered phase reappeared, indicating a reversible structural transition. The relative density of the sintered bulk samples was ∼96%. A minor secondary phase of Fe₇Se₈ (∼2.6 mol%) was detected in the FeSe sample (Fig. S2). Although Fe₇Se₈ is a ferromagnetic metallic phase with a vacancy-ordered NiAs-type crystal structure³⁴ and may, in principle, contribute to electrical transport, its small fraction suggests that its influence on the overall transport properties is expected to be limited.

Fig. 2 (a) XRD pattern of the TlFe_1.6Se₂ bulk sample at RT. (b) Enlarged view of the XRD patterns measured at different T. In the T range from RT to 100 °C, superlattice diffraction peaks of 110 and 112, originating from V_Fe ordering, were clearly observed. On the other hand, these peaks disappeared between 200 °C and 300 °C, indicating that V_Fe disordering occurred in the T range. Upon cooling back to RT, the superlattice reflections reappeared, demonstrating that the V_Fe ordered structure reversibly recovered.

Fig. 3 compares the T dependences of electrical transport properties of TlFe_1.6Se₂ and FeSe bulk polycrystals. FeSe tends to transform to the NiAs-type β-phase at ∼300 °C,³⁵ and accordingly, all measurements were performed below this temperature to retain the α-phase FeSe. For FeSe, σ exhibited only a slight T dependence; it increased from 1.7 × 10⁵ S cm⁻¹ at 50 °C to 1.9 × 10⁵ S cm⁻¹ at 290 °C (upper panel of Fig. 3a). At 50 °C, FeSe showed a small positive S of 10 µV K⁻¹ (Fig. 3b), which is consistent with the previously reported p-type conduction inferred from Hall-effect measurements.³⁶ Upon heating, however, the absolute value of S (|S|) gradually decreased and the sign eventually became negative around 200 °C, reaching −4 µV K⁻¹ at 290 °C. This sign reversal suggests that the dominant carrier type changes with T, which is reasonable given that metallic FeSe is a multiple-band system.³⁶ In contrast, TlFe_1.6Se₂ showed a significant increase in σ in the T range of 150 °C–200 °C, rising from 402 S cm⁻¹ in the V_Fe-ordered phase at 50 °C to 1064 S cm⁻¹ at 595 °C upon transition to the V_Fe-disordered phase (lower panel of Fig. 3a). Upon cooling, σ reversibly returned to its original value, reflecting the reversibility of the V_Fe order–disorder transitions. The sign of S in TlFe_1.6Se₂ was negative over the entire T range, indicating dominant electron conduction (Fig. 3b). In TlFe_1.6Se₂, V_Fe ordering induces a reconstruction of the electronic structure, associated with antiferromagnetic ordering and the opening of a band gap. The n-type behavior should not be attributed solely to Tl intercalation because the nominal composition is charge-balanced. Rather, it is more likely related to the reconstructed electronic structure associated with V_Fe ordering, possibly together with slight off-stoichiometry such as reduced V_Fe concentration or Se deficiency. At 50 °C, TlFe_1.6Se₂ exhibited a large |S| value of 110 µV K⁻¹; however, once the V_Fe order was lost through phase transition, |S| decreased substantially around 150 °C and remained almost constant at ∼34 µV K⁻¹ above 200 °C. This behavior can be attributed to the collapse of the Mott gap upon V_Fe disordering, which leads to partial compensation between the electron and hole contributions and consequently reduces the S.

Fig. 3 T dependence of the (a) electrical conductivity (σ) and (b) Seebeck coefficient (S) for the FeSe and TlFe_1.6Se₂ polycrystals. (c) T dependence of the weighted mobility (µ_w) and weighted carrier concentration (n_w) for the TlFe_1.6Se₂ polycrystal. (d) T dependence of the power factor (PF) for the TlFe_1.6Se₂ and FeSe polycrystals. The green and yellow regions indicate the T ranges where the V_Fe-ordered and V_Fe-disordered phases are stabilized in TlFe_1.6Se₂, respectively.

As a reliable Hall voltage could not be obtained by Hall effect measurements, the carrier mobility was estimated as the weighted mobility (µ_w), calculated from σ and S using the following formula:

where m_e is the free electron mass, k_B is the Boltzmann constant, h is the Planck constant, and e is the elementary charge.³⁷ The µ_w is related to the drift mobility (µ) through , where is the effective mass of the density of state. Using µ_w, the weighted carrier concentration was calculated as . In the V_Fe-ordered phase, TlFe_1.6Se₂ exhibited a high µ_w of 60 cm² (V⁻¹ s⁻¹) at 50 °C, whereas µ_w markedly decreased through transition to the V_Fe-disordered phase, reaching 8 cm² (V⁻¹ s⁻¹) at 595 °C. This reduction is likely associated with enhanced carrier scattering induced by the V_Fe disorder. The n_w in the V_Fe-ordered phase was ∼4.3 × 10¹⁹ cm⁻³, but it increased sharply during the phase-transition and then continued to increase gradually in the V_Fe-disordered phase, reaching ∼8.0 × 10²⁰ cm⁻³. This behavior is consistent with the collapse of the Mott gap upon V_Fe disordering. These results suggest that both σ and S undergo pronounced changes in response to the transition between the V_Fe ordered and disordered states, accompanied by changes in the electronic structure. Fig. 3d summarizes the T dependences of PF (=S²σ) for TlFe_1.6Se₂ and FeSe bulk samples. Bulk FeSe exhibited a small S, resulting in low PF over the entire T range. In contrast, TlFe_1.6Se₂ showed a pronounced increase in PF in the V_Fe-ordered phase, reaching 5.2 µW (cm⁻¹ K⁻²) at 50 °C, which is 30 times higher than that of bulk FeSe. Above 200 °C, where the V_Fe-disordered phase was stabilized, PF was reduced to ∼1.0 µW (cm⁻¹ K⁻²). These results demonstrate that TlFe_1.6Se₂, which contains FeSe monolayers embedded within the bulk crystal, can realize a substantially larger PF than bulk FeSe, with its electrical transport properties strongly modulated by the transition between V_Fe ordered and disordered states.

3.2 Thermal and thermoelectric properties

Next, we investigated the thermal transport properties of TlFe_1.6Se₂ and FeSe polycrystals. The T dependences of thermal diffusivity and specific heat capacity are summarized in Fig. S3a and b, SI. Because the specific heat capacity in TlFe_1.6Se₂ exhibited a peak near the phase-transition temperature (125 °C–200 °C), data in this region were excluded from the analysis of κ. Fig. 4a shows the T dependence of the total κ. FeSe exhibited a κ value of 3.2 W (m⁻¹ K⁻¹) at 50 °C, which gradually increased with increasing T and reached 4.1 W (m⁻¹ K⁻¹) at 300 °C. In contrast, TlFe_1.6Se₂ showed a much lower κ of 0.8 W (m⁻¹ K⁻¹) at 50 °C, remaining nearly constant up to ∼125 °C. After the transition to the V_Fe-disordered phase, κ slightly increased to 1.0 W (m⁻¹ K⁻¹) and then increased gradually, reaching 1.5 W (m⁻¹ K⁻¹) at 500 °C. We next discuss the thermal transport properties by separating electronic and lattice contributions. As shown in Fig. S3c in the SI, κ_ele is very large in FeSe than TlFe_1.6Se₂ due to its higher σ, and therefore plays a dominant role in determining κ. Fig. 4b compares the T dependence of κ_lat (=κ − κ_ele), in comparison with the typical high-performance thermoelectric materials, including (Bi_0.5Sb_1.5)Te₃,³⁸ (Pb_0.98K_0.02)(Te_0.15Se_0.85),³⁹ (Ge_0.9In_0.015Cu_0.125)Te,⁴⁰ Cu₂Se,⁴¹ and SnSe⁴² polycrystals. FeSe exhibited a κ_lat value of 1.8 W (m⁻¹ K⁻¹) at 50 °C, which gradually decreased with increasing T and reached 1.4 W (m⁻¹ K⁻¹) at 300 °C. In contrast, TlFe_1.6Se₂ showed an extremely low κ_lat of 0.7 W (m⁻¹ K⁻¹) at 50 °C in the V_Fe-ordered phase, where κ_lat exhibited a weak T dependence and remained nearly constant. After the transition to the V_Fe-disordered phase, κ_lat decreased further to 0.4 W (m⁻¹ K⁻¹) at 200 °C and ∼0.2 W (m⁻¹ K⁻¹) at 500 °C. In addition, these values are lower than those of high-performance thermoelectric materials with intrinsically low κ_lat, and are comparable to the lowest κ_lat values reported for layered SnSe polycrystals.⁴² Therefore, TlFe_1.6Se₂, containing FeSe layers with V_Fe order–disorder structures, is demonstrated to achieve extremely low κ_lat.

Fig. 4 T dependence of (a) κ and (b) κ_lat for the TlFe_1.6Se₂ and FeSe polycrystals. For comparison, the κ_lat vs. T curves of typical high-performance thermoelectric materials, including the (Bi_0.5Sb_1.5)Te₃,³⁸ (Pb_0.98K_0.02)(Te_0.15Se_0.85),³⁹ (Ge_0.9In_0.015Cu_0.125)Te,⁴⁰ Cu₂Se,⁴¹ and SnSe⁴² polycrystals, are also shown in (b). (c) T dependence of the dimensionless figure of merit (ZT) for the TlFe_1.6Se₂ and FeSe polycrystals. (d) Decomposition of contributions in the ZT change across the V_Fe order–disorder transition, evaluated between the V_Fe-ordered phase at 50 °C and the V_Fe-disordered phase at 400 °C. Here, Δ log(X) = log(X_ordered/X_disordered); the positive (negative) values indicate parameters that enhance (reduce) ZT in the V_Fe-ordered phase.

Here, the Lorenz number (L) for κ_ele was estimated using a single-parabolic-band model with a scattering parameter of r = −1/2. Since TlFe_1.6Se₂ is a multiband Mott system with strong electron correlations, L depends on the Fermi level, effective mass and r, which are difficult to determine uniquely because of the complex electronic structure and the uncertainty in the dominant scattering mechanism. Then, we evaluated L over a broad physically reasonable range of r (−1/2 to 2) and Fermi level. For the V_Fe-ordered phase (n = 4.2 × 10¹⁹ cm⁻³), L varies within 1.5–2.8 × 10⁻⁸ WΩ K⁻², corresponding to an uncertainty of ±30–40%. This leads to a variation of κ_ele from 0.05 to 0.11 W (m⁻¹ K⁻¹) (nominal value: 0.084 W (m⁻¹ K⁻¹)), resulting in κ_lat ranging from 0.71 to 0.79 W (m⁻¹ K⁻¹) (nominal value: 0.74 W (m⁻¹ K⁻¹)), indicating a relatively limited uncertainty. In contrast, for the V_Fe-disordered phase with a higher carrier concentration (n = 6.0 × 10²⁰ cm⁻³), L becomes much less sensitive to model parameters and converges to a narrow range of 2.4–2.5 × 10⁻⁸ WΩ K⁻², close to the degenerate limit. In this case, κ_ele is ∼1.01 W (m⁻¹ K⁻¹), with a variation of less than ±0.05 W (m⁻¹ K⁻¹), leading to only minor alternations in κ_lat (∼0.27 W (m⁻¹ K⁻¹)). Although the exact value of L depends on the model assumptions, the resulting uncertainty in κ_lat remains limited. Therefore, even when this uncertainty is taken into account, TlFe_1.6Se₂ exhibits intrinsically low κ_lat.

Based on the electronic and thermal transport properties discussed above, we evaluated the overall thermoelectric performance of TlFe_1.6Se₂ and FeSe (Fig. 4b). Bulk FeSe exhibited both a small S and a high κ, resulting in low ZT over the entire T range, consistent with its metallic electronic structure. In contrast, TlFe_1.6Se₂ showed a pronounced increase of ZT in the V_Fe-ordered phase, reaching ∼0.2 at 50 °C, which is two orders of magnitude higher than that of bulk FeSe, whereas ZT was limited to ∼0.05 in the V_Fe-disordered phase at T > 200 °C. To clarify the origin of this strong contrast in ZT between the V_Fe-ordered and V_Fe-disordered phases, the relative change in ZT across the transition was analyzed by decomposing ZT into contributions from individual transport parameters. Taking the logarithm of ZT = S²σT/κ, we obtain Δ log(ZT) = Δ log(σ) + 2Δ log(S) + Δ log(T) − Δ log(κ), which are the approximations of relative differences, . Here, the contributions are evaluated using logarithmic differences, defined as Δ log(X) = log(X_ordered) − log(X_disordered), which represents the change in each parameter across the transition. Fig. 4c summarizes the changes in each transport parameter between the V_Fe-ordered phase at 50 °C and the V_Fe-disordered phase at 400 °C. These temperatures were chosen because they lie sufficiently far from the order–disorder transition, ensuring high phase purity in the fully V_Fe-ordered and fully V_Fe-disordered states, respectively. This allows us to compare the intrinsic transport properties of the pure phases. Positive values indicate quantities that are enhanced in the V_Fe-ordered phase and therefore contribute positively to the increase in ZT, whereas negative values indicate factors that reduce ZT. Note that κ appears in the dominator of ZT, an increase of κ (i.e., positive Δ log(κ)) contributes negatively to ZT. Although κ_lat is effectively reduced in the V_Fe-disordered phase, the overall contribution of the total κ to ZT is only slightly negative. This is because the transition to a metallic state simultaneously leads to a substantial increase in κ_ele, which compensates for the benefit of the reduced κ_lat. In contrast, S exhibits the largest positive contribution, indicating that the enhancement of ZT in the V_Fe-ordered phase is predominantly driven by the large S associated with the opening of the Mott gap. To assess possible effects arising from the large temperature difference in the above comparison, the ZT decomposition analysis was also performed with a closer temperature difference of 125 °C for the ordered phase and 225 °C for the disordered phase (Fig. S4). The results confirm that the overall conclusion remains unchanged; the dominant contribution to the ZT enhancement arises from the increase in S. These results demonstrate that V_Fe ordering in layered TlFe_1.6Se₂ enables a favorable balance between electronic and thermal transport, allowing relatively large PF and ZT values to be realized in the bulk form.

To further elucidate the underlying mechanism responsible for the low κ_lat in TlFe_1.6Se₂, we conducted first-principles anharmonic lattice dynamics (ALD) calculations, as implemented in the ALAMODE code.^26,27 We first calculated the harmonic phonon dispersion of TlFe_1.6Se₂, shown as the dotted black lines in Fig. 5a. Imaginary modes appear at the X and R points, indicating that the structure is dynamically unstable at T = 0 K within the harmonic approximation. To account for the temperature effect, SCPH calculations were performed. At T = 300 K, as shown by the solid blue lines in Fig. 5a, all phonon frequencies become positive, meaning that the previously imaginary modes are stabilized after including the temperature-induced anharmonic correction. TlFe_1.6Se₂ exhibits a large number of phonon branches with strong mode entanglement (left panel of Fig. 5a). The acoustic phonons and the low-frequency optical phonons up to ∼2.5 THz exhibit relatively flat dispersions, indicating low phonon group velocities and the suppression of heat transport. The partial phonon density of states (DOSs) projected on each element for the V_Fe-ordered TlFe_1.6Se₂ at T = 300 K is shown in the right panel of Fig. 5a. The phonon DOS reveals that vibrations of the heavy Tl atoms predominantly contribute to the low-frequency phonon modes below ∼2 THz. These Tl-derived optical modes exhibit extremely flat dispersions and strongly overlap with the acoustic branches, indicating highly localized vibrations associated with weak coupling between the Tl atoms and the Fe–Se framework. In contrast, the higher-frequency phonon branches are mainly associated with vibrations of Fe and Se atoms. For comparison, the phonon dispersion and partial phonon DOS of FeSe at T = 300 K are shown in Fig. 5b. FeSe exhibits a much simpler phonon spectrum with fewer but strongly dispersive branches (left panel of Fig. 5b), consistent with its relatively simple crystal structure. While the phonon dispersions along the out-of-plane Γ–Z direction are relatively weak, much stronger dispersions are observed along the in-plane Γ–X and Γ–M direction, extending up to ∼5 THz. The broad frequency distribution of the Fe–Se vibrational modes in the partial phonon DOSs (right panel of Fig. 5b) reflects the strong Fe–Se bonding and is consistent with the highly dispersive phonon branches. Compared with bulk FeSe, the Fe–Se phonon branches in TlFe_1.6Se₂ exhibit reduced dispersion and stronger mode entanglement, reflecting the reduced symmetry and structural complexity of FeSe layers caused by V_Fe and interlayer Tl atoms.

Fig. 5 Anharmonic phonon dispersions (left panels) and partial phonon density of states (DOS) projected on each element (right panels) for (a) TlFe_1.6Se₂ and (b) FeSe at T = 300 K. The dotted black lines in (a) represent the phonon dispersion calculated within the harmonic approximation. The k-path is Γ (0.0, 0.0, 0.0), X (0.0, 0.5, 0.0), M (0.5, 0.5, 0.0), Z (0.0, 0.0, 0.5), R (0.0, 0.5, 0.5), and A (0.5, 0.5, 0.5). T dependence of the calculated κ_lat along the xx and yy (parallel to FeSe layer) and zz (perpendicular to FeSe layer) axes for (c) TlFe_1.6Se₂ and (d) FeSe. The averaged κ_lat values obtained using the series and parallel conduction models are also shown. The experimentally measured κ_lat for the V_Fe-ordered TlFe_1.6Se₂ and FeSe polycrystals are also plotted for comparison. Comparison of the κ_lat spectra for (e) TlFe_1.6Se₂ and (f) FeSe at T = 300 K. Frequency-dependent cumulative κ_lat for (g) TlFe_1.6Se₂ and (h) FeSe at T = 300 K.

We subsequently calculated κ_lat by solving the Peierls–Boltzmann transport equation (PBTE) within the relaxation time approximation. Fig. 5c and d present the anisotropic κ_lat components along the in-plane (xx and yy) and the out-of-plane (zz) directions for TlFe_1.6Se₂ and FeSe as a function of T, considering the contribution from the coherent term. The effect of the coherent contribution on the T dependence of κ_lat is shown in Fig. S5 in the SI. κ_lat for TlFe_1.6Se₂ was calculated based on the ordered structure and is therefore shown only at low temperatures. In layered FeSe, the pronounced structural anisotropy leads to a strong crystallographic-direction dependence of κ_lat (Fig. 5d). At 27 °C, the in-plane κ_{lat_xx} (κ_{lat_yy}) reaches an exceptionally high value of 33.74 W (m⁻¹ K⁻¹) within the FeSe layers, whereas κ_{lat_zz} is only 1.09 W (m⁻¹ K⁻¹), resulting in an extremely large anisotropy ratio of ∼33. In contrast, TlFe_1.6Se₂ exhibits intrinsically suppressed κ_lat even along the in-plane direction (Fig. 5c). At 27 °C, κ_{lat_xx} (κ_{lat_yy}) is only 1.43 W (m⁻¹ K⁻¹), more than an order of magnitude lower than that of FeSe, while κ_{lat_zz} is further reduced to 0.42 W (m⁻¹ K⁻¹).

To facilitate comparison with the experimental data for anisotropic polycrystalline materials, we evaluated orientation-averaged κ_lat values using both a harmonic-mean (series-type) average, , and an arithmetic-mean (parallel-type) average, , which are plotted together with the measured κ_lat of V_Fe-ordered TlFe_1.6Se₂ and FeSe polycrystals (Fig. 5c and d). The parallel-type average corresponds to a transport regime in which heat flows independently along different crystallographic directions, whereas the series-type average represents a transport regime largely controlled by the lowest κ_lat component. Both TlFe_1.6Se₂ and FeSe polycrystals exhibit κ_lat close to κ^series_lat, indicating that heat transport in the polycrystalline form is largely governed by the low out-of-plane component. In FeSe, the intrinsically high in-plane κ_lat is strongly suppressed in the polycrystal, whereas in TlFe_1.6Se₂, the already low in-plane κ_lat is further reduced by the out-of-plane contribution.

Fig. 5e and f compare the anisotropic frequency-resolved κ_lat at T = 300 K for TlFe_1.6Se₂ and FeSe, while the frequency-dependent cumulative κ_lat are shown in Fig. 5g and h. In FeSe, the in-plane κ_{lat_xx} (κ_{lat_yy}) exhibits a pronounced peak around ∼1.5 THz, whereas the out-of-plane κ_{lat_zz} shows a small peak at ∼0.5 THz (Fig. 5f). As shown in Fig. 5h, phonons below ∼5 THz account for about 90% of the total κ_{lat_xx} (κ_{lat_yy}), indicating that heat transport is dominated by acoustic phonon modes. The contribution to κ_{lat_zz} is confined to frequencies below ∼1.5 THz and arises almost exclusively from acoustic phonons, reflecting the intrinsically inefficient thermal transport along the out-of-plane direction. This anisotropic behavior is consistent with the steep phonon branches along the Γ–X and Γ–M directions, arising from highly dispersive in-plane acoustic phonons of the Fe–Se framework, and the much less dispersive phonon branches along the Γ–Z direction, reflecting the layered structural characteristic (Fig. 5b).

For TlFe_1.6Se₂, the frequency-resolved κ_lat shows a clear directional contrast (Fig. 5e). The in-plane κ_{lat_xx} (κ_{lat_yy}) spans a much broader frequency range extending up to ∼10 THz, mainly associated with the vibration of the Fe–Se framework (left panel of Fig. 5a). In contrast, the out-of-plane κ_{lat_zz} is dominated by phonons below ∼2 THz, primarily arising from the vibrations of the heavy Tl atoms (right panel of Fig. 5a). The cumulative κ_lat (Fig. 5g) further clarifies this difference. Below 2.5 THz, only ∼49% of κ_{lat_xx} (κ_{lat_yy}) is accumulated, whereas more than ∼96% of κ_{lat_zz} is contributed within this range. This indicates that higher-frequency optical modes contribute to more than half of the in-plane heat transport in TlFe_1.6Se₂. Such behavior is markedly different from that of FeSe, where the in-plane heat transport is dominated almost entirely by low-frequency phonons. Consequently, the contribution of low-frequency phonons to κ_{lat_xx} (κ_{lat_yy}) is strongly suppressed in TlFe_1.6Se₂.

Fig. 6 compares the phonon group velocity (ν_ph), phonon lifetime (τ_ph), and phonon mean free path (l_ph) as functions of phonon frequency. TlFe_1.6Se₂ exhibits much lower ν_ph than FeSe over a wide frequency range (Fig. 6a), reflecting the flatter acoustic and optical phonon branches shown in Fig. 5a. This reduction in ν_ph originates from the weak coupling between the heavy Tl atoms and the Fe–Se framework, which effectively softens the lattice vibrations and limits phonon propagation. More importantly, τ_ph in TlFe_1.6Se₂ is reduced over a broad frequency range, with an exceptionally strong suppression observed for phonons below ∼4 THz (Fig. 6b). In the lowest-frequency region below ∼1.5 THz, this reduction is primarily associated with Tl-related vibrational modes, whereas in the higher-frequency region up to ∼4 THz, it mainly affects phonons derived from the Fe–Se framework. As a result, l_ph in TlFe_1.6Se₂ is strongly suppressed, particularly because the τ_ph of phonon modes below ∼4 THz become extremely short (Fig. 6c). These results indicate that two distinct scattering mechanisms operate in different frequency regimes: low-frequency phonon scattering induced by the weakly bound heavy Tl atoms, and the additional scattering of the Fe–Se framework phonons in TlFe_1.6Se₂.

Fig. 6 (a) Phonon group velocity (ν_ph), (b) phonon lifetime (τ_ph), and (c) phonon mean free path (l_ph) in terms of the phonon frequency for TlFe_1.6Se₂ (pink plots) and FeSe (blue plots) at T = 300 K.

To further clarify the microscopic origin of the enhanced phonon scattering in the Fe–Se framework, we evaluated the Fe–Se bonding characteristics in TlFe_1.6Se₂ and FeSe by estimating their bonding energies using crystal orbital Hamilton population (COHP)⁴³ analysis, as implemented in the LOBSTER code.⁴⁴ The –iCOHP values, obtained by integrating –COHP up to the Fermi level and corresponding to the bond strength, are 2.35 eV per bond for the Fe–Se bonds in FeSe and range from 2.17 to 2.51 eV per bond in TlFe_1.6Se₂, with an average value of 2.34 eV per bond. The Fe–Se bond length in FeSe is 2.278 Å, while that in TlFe_1.6Se₂ spans a range of 2.242–2.290 Å, yielding an average bond length of 2.272 Å (Fig. S6 in SI). These results indicate that the average Fe–Se bond strength and bond length are nearly identical in TlFe_1.6Se₂ and FeSe. In contrast, TlFe_1.6Se₂ exhibits a pronounced distribution of Fe–Se bond strengths, reflecting bond heterogeneity induced by V_Fe ordering. Such bond heterogeneity provides an additional phonon-scattering channel, particularly for intermediate- and high-frequency Fe–Se vibrational modes, and plays a key role in shortening τ_ph and suppressing κ_lat. This is further supported by the interatomic force constant (IFC) distributions (Fig. S7 in SI), showing the broadened force-constant distribution in TlFe_1.6Se₂ compared to FeSe. In addition, the three-phonon scattering phase space (SPS) is significantly enhanced, particularly in the 1.5–4 THz region (Fig. S8 in SI), indicating an increased number of scattering channels in this frequency range. This frequency range is dominated by the Fe–Se vibrational modes with strong acoustic–optical mixing (Fig. 5a), making them sensitive, especially to the bond heterogeneity. The dominance of force-constant disorder is also inferred from the frequency-selective reduction of phonon lifetimes between 1.5 and 4 THz (Fig. 6b). Therefore, the lifetime reduction in this frequency window is primarily driven by emission channels activated by the force-constant disorder, while mass disorder provides a secondary background contribution. Such bond-heterogeneity-induced phonon scattering has also been reported in a variety of complex crystalline materials, where local fluctuations in bonding strength or force constants give rise to enhanced anharmonic phonon scattering and reduced κ_lat.^45–49 In the V_Fe-disordered phase, this bond heterogeneity is further amplified by the random distribution of V_Fe, introducing additional mass and force-constant disorder. This disorder-driven enhancement of phonon scattering accounts for the further reduction of κ_lat observed experimentally at elevated temperatures. Direct first-principles calculations for the V_Fe disordered phases were not performed due to the large configurational space associated with random vacancy distributions and the limitation of our computational resources. Random V_Fe distributions are expected to introduce spatially uncorrelated force-constant disorder, which is likely more effective in scattering phonons than the ordered structure, leading to additional suppression of κ_lat.

4 Conclusions

We investigated the thermoelectric properties of layered TlFe_1.6Se₂, motivated by the expectation that its embedded two-dimensional FeSe transport layers, together with the V_Fe order–disorder transition, could provide a unique combination of enhanced PF and suppressed κ_lat. TlFe_1.6Se₂ exhibited an exceptionally low κ_lat, reaching ∼0.7 W (m⁻¹ K⁻¹) in the V_Fe-ordered phase at 50 °C and further decreasing to ∼0.2 W (m⁻¹ K⁻¹) in the V_Fe-disordered phase at 500 °C, which is comparable to or even lower than those of high-performance thermoelectric materials with intrinsically low κ_lat. Notably, unlike layered FeSe, which shows extremely high in-plane κ_lat and strong anisotropy, κ_lat in TlFe_1.6Se₂ is intrinsically suppressed even along the in-plane direction. This ultralow κ_lat originates from the combined effects of weakly bound heavy Tl atoms, which reduce phonon group velocities, and enhanced phonon scattering caused by the V_Fe-induced bond heterogeneity within the Fe–Se framework.

The electronic transport properties are strongly modulated across the transition between the V_Fe ordered and disordered states. In the V_Fe-disordered phase, TlFe_1.6Se₂ exhibits a metallic state with enhanced σ and ultralow κ_lat; however, strong carrier compensation significantly reduces S, limiting the ZT value to ∼0.05. In contrast, in the V_Fe-ordered phase, the opening of a Mott gap leads to a large S of ∼110 µV K⁻¹ at 50 °C, resulting in an enhanced PF of 5.2 µW (cm⁻¹ K⁻²) and a ZT of ∼0.2. These results demonstrate that the confinement of FeSe layers within bulk TlFe_1.6Se₂ enables a substantial PF enhancement compared with bulk FeSe, while V_Fe ordering and disordering play complementary roles by tuning carrier transport and strongly suppressing κ_lat. Looking forward, further enhancement of ZT may be achieved by optimizing carrier concentration in the V_Fe-ordered phase and by tuning the V_Fe concentration.

Although the incorporation of two-dimensional FeSe layers suggested the possibility of achieving a high PF comparable to that of FeSe monolayers, the obtained PF is substantially lower. This difference is primarily attributed to the much lower carrier mobility in TlFe_1.6Se₂, which is about four orders of magnitude smaller than that of FeSe monolayers, likely due to the polycrystalline nature and enhanced carrier scattering arising from V_Fe ordering/disordering and interactions involving Tl atoms. In contrast, the exceptionally high mobility in FeSe monolayers is realized for in-plane transport in the two-dimensional electronic structure. To clarify the intrinsic transport properties and to potentially realize higher PF in TlFe_1.6Se₂, studies using single crystals or epitaxial thin films will be essential, as they enable selective probing of in-plane transport in the embedded FeSe layers.

Extending this concept to other alkali-based A_1−xFe_2−ySe₂ (A = K, Rb, and Cs) is an attractive direction, as they consist of FeSe layers with V_Fe, forming quasi-two-dimensional frameworks, and exhibit high electrical conductivity. Moreover, these systems show high-T_c superconductivity, reminiscent of FeSe monolayers, suggesting a close connection between their electronic structures. In addition, the V_Fe concentration in these systems can be tuned over a wide range,^16,17 providing a versatile platform for controlling carrier concentration and optimizing thermoelectric performance. However, their poor air stability, precise control of vacancy ordering, and suppression of phase separation remain significant challenges.¹⁵ The growth of high-quality single crystals or epitaxial thin films will also be essential to clarify the intrinsic transport properties. A possible strategy to address air stability is partial substitution or solid-solution formation with air-stable Tl-based compounds,⁵⁰ which may enable simultaneous control of environmental stability and thermoelectric performance.

Taken together, the present study highlights the confinement of two-dimensional FeSe layers, combined with vacancy-order control, as a promising design strategy for bulk thermoelectric materials that combine large PF with intrinsically low κ_lat.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this study are included within the article and its supplementary information (SI). Supplementary information: the SI includes the bulk synthesis, characterization, and theoretical calculations for TlFe_1.6Se₂ and FeSe. See DOI: https://doi.org/10.1039/d6ta02075e.

Acknowledgements

This work was supported by the MEXT Program: Data Creation and Utilization Type Material Research and Development Project (Grant No. JPMXP1122683430), the Design and Engineering by Joint Inverse Innovation for Materials Architecture, and also by a project of Kanagawa Institute of Industrial Science and Technology (KISTEC). X. He was supported by Japan Society for the Promotion of Science (JSPS) through the Grant-in-Aid for Research Activity Start-up (Grant No. 25K23539). T. Katase was supported by the Special Award for Science Tokyo Advanced Researchers (STAR) funded by the Institute of Science Tokyo and JSPS through the Grant-in-Aids for Scientific Research (B) (Grant No. 26K01206), Scientific Research (A) (Grant No. 24H00314), Scientific Research (S) (Grant No. 22H04964), and Challenging Research (Exploratory) (Grant No. 24K21671). H. Hi. was supported by JSPS through the Grants-in-Aid for Scientific Research (A) (Grant No. JP20H00302, JP21H04612, and JP24H00376). The numerical calculations were carried out on the TSUBAME4.0 supercomputer at the Institute of Science Tokyo supported by the MEXT Project of the Tokyo Tech Academy for Convergence of Materials and Informatics (TAC-MI), and the supercomputer at the Research Center for Computational Science, Okazaki, Japan. The crystal structures in Fig. 1 were drawn using the VESTA code.⁵¹

References

1. F. J. DiSalvo, Science, 1999, 285, 703–706.
2. L. E. Bell, Science, 2008, 321, 1457–1461.
3. Q. H. Zhang, X. Y. Huang, S. Q. Bai, X. Shi, C. Uher and L. D. Chen, Adv. Eng. Mater., 2016, 18, 194–213.
4. G. J. Snyder and E. S. Toberer, Nat. Mater., 2008, 7, 105–114.
5. S. Shimizu, J. Shiogai, N. Takemori, S. Sakai, H. Ikeda, R. Arita, T. Nojima, A. Tsukazaki and Y. Iwasa, Nat. Commun., 2019, 10, 825.
6. F.-C. Hsu, J.-Y. Luo, K.-W. Yeh, T.-K. Chen, T.-W. Huang, P. M. Wu, Y.-C. Lee, Y.-L. Huang, Y.-Y. Chu, D.-C. Yan and M.-K. Wu, Proc. Natl. Acad. Sci. U. S. A., 2008, 105, 14262–14264.
7. J.-F. Ge, Z.-L. Liu, C. Liu, C.-L. Gao, D. Qian, Q.-K. Xue, Y. Liu and J.-F. Jia, Nat. Mater., 2014, 14, 285–289.
8. Y. Machida, K. Tomokuni, T. Isono, Y. Nakajima and T. Tamegai, J. Phys. Soc. Jpn., 2009, 78, 073705.
9. J. G. Checkelsky, R. Thomale, L. Li, G. F. Chen, J. L. Luo, N. L. Wang and N. P. Ong, Phys. Rev. B:Condens. Matter Mater. Phys., 2012, 86, 180502(R).
10. P. Devi Lodhi, N. Kaurav, K. K. Choudhary and Y. K. Kuo, J. Low Temp. Phys., 2019, 196, 494–509.
11. C. B. Satterthwaite and R. W. Ure, Phys. Rev., 1957, 108, 1164–1170.
12. F. Ye, S. Chi, W. Bao, X. F. Wang, J. J. Ying, X. H. Chen, H. D. Wang, C. H. Dong and M. Fang, Phys. Rev. Lett., 2011, 107, 137003.
13. R. Yu, J.-X. Zhu and Q. Si, Phys. Rev. Lett., 2011, 106, 186401.
14. X.-W. Yan, M. Gao, Z.-Y. Lu and T. Xiang, Phys. Rev. Lett., 2011, 106, 087005.
15. W. Li, H. Ding, P. Deng, K. Chang, C. Song, K. He, L. Wang, X. Ma, J.-P. Hu, X. Chen and Q.-K. Xue, Nat. Phys., 2012, 8, 126–130.
16. J. Guo, S. Jin, G. Wang, S. Wang, K. Zhu, T. Zhou, M. He and X. Chen, Phys. Rev. B:Condens. Matter Mater. Phys., 2010, 82, 180520.
17. Y. J. Yan, M. Zhang, A. F. Wang, J. J. Ying, Z. Y. Li, W. Qin, X. G. Luo, J. Q. Li, J. Hu and X. H. Chen, Sci. Rep., 2012, 2, 212.
18. B. C. Sales, M. A. McGuire, A. F. May, H. Cao, B. C. Chakoumakos and A. S. Sefat, Phys. Rev. B:Condens. Matter Mater. Phys., 2011, 83, 224510.
19. A. F. May, M. A. McGuire, H. Cao, I. Sergueev, C. Cantoni, B. C. Chakoumakos, D. S. Parker and B. C. Sales, Phys. Rev. Lett., 2012, 109, 077003.
20. T. Katase, H. Hiramatsu, T. Kamiya and H. Hosono, Proc. Natl. Acad. Sci., 2014, 111, 3979–3983.
21. H. J. Goldsmid, Introduction to Thermoelectricity, Springer, Heidelberg, 2009.
22. X. He, S. Kimura, T. Katase, T. Tadano, S. Matsuishi, M. Minohara, H. Hiramatsu, H. Kumigashira, H. Hosono and T. Kamiya, Adv. Sci., 2023, 11, 2307058.
23. G. Kresse and D. Joubert, Phys. Rev. B:Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
24. G. Kresse and J. Furthmüller, Phys. Rev. B:Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
25. J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou and K. Burke, Phys. Rev. Lett., 2008, 100, 136406.
26. T. Tadano and S. Tsuneyuki, Phys. Rev. B:Condens. Matter Mater. Phys., 2015, 92, 054301.
27. T. Tadano, Y. Gohda and S. Tsuneyuki, J. Phys.: Condens. Matter, 2014, 26, 225402.
28. G. Kresse, J. Furthmüller and J. Hafner, Europhys. Lett., 1995, 32, 729–734.
29. K. Parlinski, Z. Li and Y. Kawazoe, Phys. Rev. Lett., 1997, 78, 4063–4066.
30. T. Tadano and S. Tsuneyuki, Phys. Rev. B:Condens. Matter Mater. Phys., 2015, 92, 054301.
31. M. Simoncelli, N. Marzari and F. Mauri, Nat. Phys., 2019, 15, 809–813.
32. S. Tamura, Phys. Rev. B:Condens. Matter Mater. Phys., 1983, 27, 858–866.
33. H. Cao, C. Cantoni, A. F. May, M. A. McGuire, B. C. Chakoumakos, S. J. Pennycook, R. Custelcean, A. S. Sefat and B. C. Sales, Phys. Rev. B:Condens. Matter Mater. Phys., 2012, 85, 054515.
34. G. Li, B. Zhang, T. Baluyan, J. Rao, J. Wu, A. A. Novakova, P. Rudolf, G. R. Blake, R. A. de Groot and T. T. M. Palstra, Inorg. Chem., 2016, 55, 12912–12922.
35. T. M. McQueen, Q. Huang, V. Ksenofontov, C. Felser, Q. Xu, H. Zandbergen, Y. S. Hor, J. Allred, A. J. Williams, D. Qu, J. Checkelsky, N. P. Ong and R. J. Cava, Phys. Rev. B:Condens. Matter Mater. Phys., 2009, 79, 014522.
36. M. Kawai, F. Nabeshima and A. Maeda, J. Phys.: Conf. Ser., 2018, 1054, 012023.
37. G. J. Snyder, A. H. Snyder, M. Wood, R. Gurunathan, B. H. Snyder and C. Niu, Adv. Mater., 2020, 32, 2001537.
38. T. Zhu, Z. Xu, J. He, J. Shen, S. Zhu, L. Hu, T. M. Tritt and X. Zhao, J. Mater. Chem. A, 2013, 1, 11589–11594.
39. Q. Zhang, F. Cao, W. Liu, K. Lukas, B. Yu, S. Chen, C. Opeil, D. Broido, G. Chen and Z. Ren, J. Am. Chem. Soc., 2012, 134, 10031–10038.
40. Q. Zhang, Z. Ti, Y. Zhu, Y. Zhang, Y. Cao, S. Li, M. Wang, D. Li, B. Zou, Y. Hou, P. Wang and G. Tang, ACS Nano, 2021, 15, 19345–19356.
41. M. Li, S. M. Kazi Nazrul Islam, M. Yahyaoglu, D. Pan, X. Shi, L. Chen, U. Aydemir and X. Wang, INFOMAT, 2019, 1, 108.
42. C. Zhou, Y. K. Lee, Y. Yu, S. Byun, Z.-Z. Luo, H. Lee, B. Ge, Y.-L. Lee, X. Chen, J. Y. Lee, O. Cojocaru-Mirédin, H. Chang, J. Im, S.-P. Cho, M. Wuttig, V. P. Dravid, M. G. Kanatzidis and I. Chung, Nat. Mater., 2021, 20, 1378–1384.
43. R. Dronskowski and P. E. Bloechl, J. Phys. Chem., 1993, 97, 8617.
44. S. Maintz, V. L. Deringer, A. L. Tchougréeff and R. Dronskowski, J. Comput. Chem., 2016, 37, 1030.
45. T. Katase, N. Sato and T. Mori, Appl. Phys. Rev., 2025, 12, 041319.
46. D. T. Morelli, V. Jovovic and J. P. Heremans, Phys. Rev. Lett., 2008, 101, 035901.
47. M. Dutta, K. Pal, U. V. Waghmare and K. Biswas, Chem. Sci., 2019, 10, 4905–4913.
48. N. Sato, N. Kuroda, S. Nakamura, Y. Katsura, I. Kanazawa, K. Kimura and T. Mori, J. Mater. Chem. A, 2021, 9, 22660–22669.
49. X. He, S. Nomoto, T. Komatsu, T. Katase, T. Tadano, S. Kitani, H. Yoshida, T. Yamamoto, H. Mizoguchi, K. Ide, H. Hiramatsu, H. Kawaji, H. Hosono and T. Kamiya, Adv. Funct. Mater., 2023, 33, 2213144.
50. M.-H. Fang, H.-D. Wang, C.-H. Dong, Z.-J. Li, C.-M. Feng, J. Chen and H. Q. Yuan, Europhys. Lett., 2011, 94, 27009.
51. K. Momma and F. Izumi, J. Appl. Crystallogr., 2011, 44, 1272–1276.

---

Footnote

† X. H. and K. O. contributed equally to this work.

---