Synergistic effect of cation substitution and p-type anion doping to improve thermoelectric properties in Zintl phases

Synthesis

Sample preparation was carried out within an argon-filled glovebox, where oxygen and moisture concentrations were kept below 0.1 ppm, or under vacuum. All reactants were bought from Alfa Aesar, and these included Ca (shot, 99.5%), Yb (ingot, 99.9%), Zn (shot, 99.99%), Cu (shot, 99.5%), Sb (shot, 99.9999%), and Pb (granules, 99.99%). The surfaces of Ca and Yb, which had become tarnished, were scrubbed using a metal brush inside an argon-filled glovebox before use. All reactant sets were introduced into alumina crucibles, accompanied by an excess quantity of Pb metal used as the metal-flux in the following ratios: Ca:Yb:Zn:Cu:Sb:Pb = 9:0:4.5:0:9:45, 8.5:0.5:4.5:0:9:45, 8:1:4.5:0:9:45, 7.5:1.5:4.5:0:9, 8:1:4.45:0.05:9:45, 8:1:4.4:0.1:9:45, and 8:1:4.35:0.15:9:45. Individual alumina crucibles were inserted into separate fused-silica tubes which were subsequently sealed under vacuum conditions to prevent oxidation of reactants during processing at elevated temperatures. Then, the reactants were annealed in a furnace ramped up to 1323 K at 200 K/h, kept there for 24 h, and then cooled down to 873 K at a rate of 5 K/h in 90 h. Afterward, the reaction ampoule was quickly taken out of the furnace and centrifuged for 3 min to separate the crystallized products from the molten Pb-metal flux. A small amount of remaining Pb metal after the centrifugation was further dissolved by soaking crystals in H₂O₂ for 48 h, and the needle/bar-shaped nicely grown single crystals were eventually obtained as shown in Figure 1.

Figure 1. OM images for the needle/bar-shaped single-crystals of (A) Ca₈YbZn_4.45Cu_0.05Sb₉, (B) Ca₈YbZn_4.4Cu_0.1Sb₉, and (C) Ca₈YbZn_4.35Cu_0.15Sb₉. Scale bars are also provided.

Four of the title compounds were additionally prepared via ball-milling followed by a hot-pressing method to obtain larger quantities of samples suitable for TE property measurements. Initially, each reactant mixture corresponding to the target composition was placed into a stainless-steel container along with two 0.5-inch and two 0.25-inch stainless-steel balls inside an argon-filled glovebox. The mixtures were then ball-milled using a SPEX 8000M milling machine for 3 h. During this procedure, milling was paused every 30 min to scrape and thoroughly remix the powders, ensuring sample homogeneity. After completion of ball-milling, the resulting powdered products were consolidated into disk-shaped samples by hot-pressing under the following conditions: first at 623 K under 60 MPa for 2 h, and secondly at 1,073 K under the same pressure for an additional 2 h.

X-ray diffraction analysis

PXRD analysis was initially used to evaluate the phase purity of the five title compounds in the Ca_9-xYb_xZn_4.5-yCu_ySb₉ system (0.57(2) ≤ x ≤ 1.28(2), 0 ≤ y ≤ 0.15). The diffraction data were collected using a Bruker D8 diffractometer with monochromatic Cu K_α1 radiation (λ = 1.54059 Å). The data acquisition was conducted for 30 min per sample, utilizing a step size of 0.02°, covering the angular range of 20° ≤ 2θ ≤ 80°.

The phase purities were first examined by comparing the measured diffraction patterns with a calculated pattern based on the SXRD refinement results of Ca₉Zn_4.48(2)Sb₉; this comparison is shown in Figure 2. In addition, the Rietveld refinement was conducted for the collected PXRD patterns of three Zn/Cu mixture-containing quinary compounds using PDXL2 software to obtain the lattice parameters and volumes of the unit cells. The refinement results and lattice parameters are presented in Supplementary Table 1 and Supplementary Figure 1.

Figure 2. (A-F) PXRD patterns of five title compounds in the Ca_9-xYb_xZn_4.5-yCu_ySb₉ (0.57(2) ≤ x ≤ 1.28(2), 0 ≤ y ≤ 0.15) system. A PXRD pattern calculated from the SXRD data of CaZn_4.48(2)Sb₉ is also provided for a comparison purpose.

SXRD data was collected at room temperature for one ternary and three quaternary compounds: Ca₉Zn_4.48(2)Sb₉, Ca_8.42(2)Yb_0.57Zn_4.46(1)Sb₉, Ca_8.01(3)Yb_0.99Zn_4.46(1)Sb₉, and Ca_7.72(2)Yb_1.28Zn_4.45(2)Sb₉. The measurements were performed using a Bruker SMART APEX2 CCD-based diffractometer with Mo K_α1 radiation (λ = 0.71073 Å). Several single crystal specimens were initially screened for quality, and the most suitable specimen was selected for a thorough analysis using Bruker's APEX2 program^[20]. Data reduction, integration, determination of unit cell parameters, and atomic site refinement were performed using the SAINT program^[21]. Semi-empirical absorption corrections based on equivalents were applied using the SADABS program^[22]. The entire reflection sets from these four products nicely matched with the orthorhombic Pbam space group (no. 55), and the isotypic crystal structures eventually crystallized in the Ca₉Mn₄Bi₉-type phase. The crystal structures were determined using a direct method and refined to convergence using full-matrix least-squares methods on F². The refinement included several parameters: the scale factor, atomic positions with anisotropic displacement parameters (ADPs), extinction coefficients, and occupancy factors for the Ca/Yb mixed sites. The STRUCTURE TIDY program^[23] was employed for the standardization of atomic positions in the final refinement step. Comprehensive crystallographic data, atomic positions with ADPs, and selected interatomic distances are presented in Tables 1-3. All these crystallographic data can be obtained free of charge from The Cambridge Crystallographic Data Centre via www.ccdc.cam.ac.uk/data_request/cif. The depository numbers are as follows: CCDC- 2422809 for Ca₉Zn_4.48(2)Sb₉, CCDC-2422810 for Ca_8.43(2)Yb_0.57Zn_4.46(1)Sb₉, CCDC- 2422811 for Ca_8.01(3)Yb_0.99Zn_4.46(1)Sb₉, and CCDC- 2422812 for Ca_7.72(2)Yb_1.28Zn_4.45(2)Sb₉.

DFT calculation

To understand how atomic substitutions influence the electronic structure of the studied system, DFT calculations were carried out by employing the TB-LMTO-ASA method^[24-27] with two hypothetical models. For practical reasons, these two hypothetical models were designed to have the monoclinic subgroup Pc (no. 7) rather than the refined space group orthorhombic Pnma (no. 55) to apply the idealized compositions of Ca₉Zn_4.5Sb₉ and Ca₈YbZn₄Cu_0.5Sb₉. Lattice parameters and atomic sites information necessary to build these models were taken from the SXRD refinement results of Ca₉Zn_4.48(2)Sb₉ and Ca_8.01(3)Yb_0.99Zn_4.46(1)Sb₉, respectively. The obtained DOS, COHP curves, band structures, and ELF illustrations were carefully investigated. Detailed crystallographic information on these models is presented in Supplementary Table 2.

All relativistic effects were incorporated using a scalar relativistic approximation, except for spin-orbit coupling. The potential inside each Wigner-Seitz (WS) sphere was treated as spherically symmetric^[28], and a combined correction was applied to account for the overlapping regions. Empty spheres were automatically included, and the space was filled with overlapping WS atomic spheres^[28]. The WS sphere radii were automatically calculated to optimize the overlapping potential's approximation to the full potential. The WS radii used for the two models are as follows: Ca, 1.891-2.158 Å; Zn, 1.335-1.504 Å; Sb, 1.631-1.911 Å for Ca₉Zn_4.5Sb₉; Ca, 1.891-2.158 Å; Yb, 1.975-2.040 Å; Zn, 1.335-1.504 Å; Cu, 1.465 Å; Sb, 1.631-1.911 Å for Ca₈YbZn₄Cu_0.5Sb₉. The used basis sets were 4s, 4p and 3d orbitals for Ca, 6s, 6p and 5d orbitals for Yb, 4s, 4p and 3d orbitals for Zn, 4s, 4p and 3d orbitals for Cu, 5s, 5p, 5d and 4f orbitals for Sb. The Löwdin downfolding technique^[29] was employed for the Ca 4p, Yb 6p, and Sb 5d and 4f orbitals. Integration in k-space was performed using the tetrahedron method, and the self-consistent charge density was computed across the first Brillouin zone using 320 irreducible k-points for both models^[30].

TGA

The TGA, performed on a Netzsch TG 209 F1 Libra instrument, was utilized to evaluate the thermal stability of three synthesized compounds. Each pulverized sample (around 20 mg) was positioned in an alumina pan. Heating proceeded to 1,123 K under a constant nitrogen flow, with the temperature ramped at 10 K/min. Following the heating cycle, the samples underwent natural cooling back to ambient temperature. The resulting TGA curves [Supplementary Figure 2] indicated that the three compounds maintained thermal integrity up to roughly 900 K.

Electrical transport property measurement

First, disk-shaped specimens were prepared using a two-step pressing process: initial hot-pressing at 623 K for 2 h under 60 MPa, followed by pressing at 1,073 K for 2 h. From these disks, square bars (3 mm × 3 mm × 10 mm) were cut and polished for the electrical measurements. It was verified through geometric density checks that all samples achieved densities above 96%. Under a helium atmosphere, the temperature-dependent electrical conductivity σ and Seebeck coefficient S of the four compounds were measured simultaneously from 323 to 861 K with a ULVAC-RIKO ZEM-3 system.

Thermal conductivity measurement

Using the flash diffusivity method, thermal diffusivity was evaluated for the four disk-shaped title compounds between 323 and 861 K. The measurements were performed with a Netzsch LFA 467 HyperFlash system in an argon environment, where a laser pulse heats the sample's front surface, and an Infrared (IR) detector tracks the subsequent temperature rise on the back surface. The thermal conductivity κ was determined using κ = DC_pρ, where D = thermal diffusivity, C_p = heat capacity, and ρ = density^[31]. The Dulong-Petit value (3R/atom, R = gas constant) was used for C_p^[31]. The total thermal conductivity κ_tot was calculated by the sum of the lattice κ_latt and the electronic κ_elec thermal conductivities^[32]. The electronic thermal conductivity was evaluated using the Wiedemann-Franz law: κ_elec = LσT, where L is the temperature-dependent Lorenz number, which was derived from the single-parabolic band model using the temperature-dependent Seebeck coefficient S^[32]. Subsequently, the lattice thermal conductivity was calculated using: κ_latt = κ_tot - κ_elec.

EDS and EPMA analysis

The elemental composition and distribution within the compounds Ca₈YbZn_4.35Cu_0.05Sb₉, Ca₈YbZn_4.4Cu_0.1Sb₉, and Ca₈YbZn_4.35Cu_0.15Sb₉ were investigated using EDS. The measurements were carried out using a ULTRA Plus FE-SEM, with the operating condition set to an acceleration potential of 20 kV. Well-grown single crystals synthesized via the molten Pb-flux technique were affixed to an aluminum sample holder under an argon atmosphere. Results from the EDS analysis, presented in Supplementary Figure 3, indicated a homogeneous distribution of Cu within these single crystals. In particular, the quantitative elemental analyses of these three title compounds were also conducted by the JEOL JXA-IHP200F EPMA using an acceleration voltage of 20 KeV with a current of 50 nA in a WDS mode for randomly selected locations. EPMA analysis resulted in the chemical compositions of Ca_8.25Yb_0.86Zn_4.51Cu_0.07Sb_8.81, Ca_8.36Yb_1.00Zn_4.49Cu_0.11Sb_8.54, and Ca_8.11Yb_1.01Zn_4.54Cu_0.14Sb_8.69 for Ca₈YbZn_4.45Cu_0.05Sb₉, Ca₈YbZn_4.4Cu_0.1Sb₉, and Ca₈YbZn_4.35Cu_0.15Sb₉, respectively.

Crystal structure analysis

Seven title compounds in the Ca_9-xYb_xZn_4.5-yCu_ySb₉ (0 ≤ x ≤ 1.50, 0 ≤ y ≤ 0.15) system were prepared using the molten Pb-flux method. Initial structural assessments, including the determination of phase purities and lattice parameters of all compounds, were conducted through Rietveld refinements for the collected PXRD data. As displayed in Figure 2 and Supplementary Figure 1, all seven compounds crystallize as single-phase products. The refined lattice parameters were consistent with those obtained from SXRD data refinement. Further detailed structural refinements for one ternary and three quaternary Zn-containing title compounds were conducted by SXRD analysis, and the resultant crystallographical details are listed in Tables 1-3. The studied compounds adopted the Ca₉Mn₄Bi₉-type phase having an orthorhombic Pbam space group (Pearson code oP44, Z = 2)^[13], and thirteen crystallographically independent atomic sites including the five Ca or Ca/Yb mixed-sites, three Zn or Zn/Cu mixed-sites, and five Sb sites were refined and presented in Table 2. However, the detailed chemical compositions of the three quinary compounds containing both Zn and Cu were not accurately refined using SXRD data since X-ray scattering factors of Zn and Cu were similar. Instead, the detailed chemical compositions of these compounds were evaluated by EPMA analysis, and the results are as follows: Ca_8.25Yb_0.86Zn_4.51Cu_0.07Sb_8.81, Ca_8.36Yb_1.00Zn_4.49Cu_0.11Sb_8.54, and Ca_8.11Yb_1.01Zn_4.54Cu_0.14Sb_8.69 for Ca₈YbZn_4.45Cu_0.05Sb₉, Ca₈YbZn_4.4Cu_0.1Sb₉, and Ca₈YbZn_4.35Cu_0.15Sb₉ (all nominal compositions).

Since all seven title compounds are isostructural, the detailed overall crystal structure can be explained using a representative quaternary compound Ca_8.43(2)Yb_0.57Zn_4.46(1)Sb₉ in Figure 3. The presented crystal arrangement can basically be elucidated as an assembly of (1) a complex 2-dimensional (2D) layered anionic structure; and (2) two different types of electron-donating cationic elements. Interestingly, the 2D anionic layered structure can be built up from its basic building block of the [ZnSb₄] tetrahedron, which initially forms the 1D $${^{1}_{\infty}}$$[ZnSb₃] chains propagating along the c-axis direction by sharing Sb2 with its neighboring tetrahedra [Figure 3B]. These 1D chains are linked together via Sb4 to build the $${^{1}_{\infty}}$$[Zn₂Sb₅] double chains. Then, two such double chains are further joined via Sb5 in the form of the inversion symmetry resulting in building the $${^{1}_{\infty}}$$[Zn₄Sb₉] cluster.

Figure 3. (A) A crystal structure of Ca_8.43(2)Yb_0.57Zn_4.46(1)Sb₉ illustrated by a combination of ball-and-stick and polyhedral representations viewed down from the c-axis direction. (B) A schematic presentation of the assembly of the 2D anionic layered structure. Color codes: A (Ca/Yb mixed-site): gray; Zn: orange; Sb: blue.

Finally, these anionic clusters are connected to one another via the bridging Zn1 along the a-axis directions to complete the complete 2D layered structure. In particular, these two different types of Zn sites were refined in this anionic framework with distinctive local coordination geometries and different occupations: (1) the Zn1 site (Wyckoff 4h) was located at the center of the trigonal planar site with a partial occupation of ca. 24%; and (2) the Zn2 and Zn3 (both Wyckoff 4g) were found to be at the center of the tetrahedral site with full occupations. Interestingly, the Zn1 site showed a rather larger atomic displacement parameter (ADP) value than two other Zn sites as given in Table 2. This kind of large ADP value was already reported in other isotypic Ca₉Mn₄Bi₉-type compounds, including Ca₉Zn_4.10(1)Bi₉, Sr₉Cd_4.49(1)Sb₉, Ca₉Zn_4.478(8)Sb₉, and Yb₉Zn_4.380(12)Sb₉^[33,34]. Thus, it can be considered as an intrinsic structural feature of the title phase due to the partial Zn1 occupation. Furthermore, based on the Zintl-Klemm concept, we could rewrite the chemical formula of the title compounds as: [(Ca/Yb)²⁺]₁₈[(4b-Zn^2-)₈(3b-Zn^1-)(1b-Sb^2-)₃(2b-Sb^1-)₁₃(3b-Sb⁰)₂]. Therefore, we believe that the 1/4 occupation at the Zn1 site as shown in the formula should occur to optimize the overall charge balance between cations and anions. In addition, the charge balance evaluation based on the coordination geometry is also nicely supported by the ELF analysis, which will be explained further in the Electronic Structure Analysis section.

Five different types of cationic sites existed in the title system (See Figure 4): (1) the two intra-layer A3 and A5 sites located inside the “channel-like” voids within the anionic 2D layers; and (2) the three inter-layer A1, A2, and A4 sites situated in-between two such neighboring 2D layers as displayed in Figure 3. All these cationic sites were primarily coordinated by six Sb, which formed a slightly distorted octahedral geometry, and depending on the given local coordination environment, some of the triangular faces were capped by Zn atoms resulting in the various site volumes. Interestingly, in the two quaternary compounds having the Ca/Yb mixed-sites, Yb showed a specific affinity for the A3 site among the five possible cationic sites.

Figure 4. Illustration of the five distinct cationic environments identified within the Ca_8.43(2)Yb_0.57Zn_4.46(1)Sb₉ structure. Detailed labeling of atoms and chosen interatomic distance are also shown.

The observed site-preference behavior is typically influenced by one of the following main factors: (1) the size factor, which reflects the compatibility between the central cation and the site volume; or (2) the electronic factor represented by the QVAL^[17], which is the relationship between the electronegativity of a cation and an electron density at the site. In our title system, the difference in size between two divalent cations Ca²⁺ and Yb²⁺ is quite subtle [r(Ca²⁺) = 1.00 Å and r(Yb²⁺) = 1.02 Å]^[35]. Consequently, the electronic factor should be considered the determining criterion for Yb's observed preference for the A3 site. A detailed analysis of the correlation between the electronic factor and cationic site-preference will be a part of the Electronic Structure section.

Electronic structure and chemical bonding analysis

To understand how the co-substitution of Yb for Ca and Cu for Zn affects the electronic structure of the Ca_9-xYb_xZn_4.5-yCu_ySb₉ system, DFT calculations using the TB-LMTO-ASA method were performed^[24-27]. Therefore, two hypothetical structural models containing only Yb or both Yb and Cu were systematically designed with the idealized compositions of Ca₈YbZn_4.5Sb₉ and Ca₈YbZn₄Cu_0.5Sb₉, and their electronic structures were calculated and compared with that of the parental compound Ca₉Zn_4.5Sb₉. To incorporate these specific compositions and atomic arrangements in the models, the symmetry was reduced from the experimentally refined space group Pbam (no. 55), which contains 13 crystallographically independent atomic sites, to its subgroup Pc (no. 7), which contains 48 atomic sites. In particular, three out of four Zn1 sites were treated as vacancies to mimic the partial occupation observed in the SXRD refinements. For the quinary Ca₈YbZn₄Cu_0.5Sb₉ model, Cu with a higher electronegativity than Zn was assigned to partially occupy the Zn2 site (Wyckoff 4g site) showing the largest QVAL amidst the three available Zn sites. Additional crystallographic details, including lattice parameters and atomic coordinates, were obtained from the SXRD results of Ca₉Zn_4.48(2)Sb₉ and Ca_8.01(3)Yb_0.99Zn_4.46(1)Sb₉, respectively. Supplementary Table 2 presents further details for the two studied hypothetical models.

The calculated total DOS (TDOS) and partial DOS (PDOS) curves of Ca₉Zn_4.5Sb₉ and Ca₈YbZn₄Cu_0.5Sb₉ resembled each other since they adopted the identical Ca₉Mn_4.5Bi₉-type crystal structure [Figure 5A]. Strong orbital mixing among all elements was observed across the entire energy window, and the pseudogap noted at the Fermi level (E_F) was indicative of the poor metallic properties of these title compounds. In particular, the slightly increased metallic properties should be expected for the p-type Cu-substituting compounds since, as shown in Figure 5B, the DOS level at E_F increased in the Ca₈YbCu_0.5Sb₉ model. In addition, the wide p-orbital states ranging between ca. -6.5 and 0 eV can be interpreted by dividing into two sub-sections: (1) the lower section between ca. -6.5 and -4.5 eV, where Zn2 and Zn3 showed relatively larger contributions; and (2) the higher section between ca. -4.5 and 0 eV, where Zn1 displayed a relatively stronger contribution. Two or three COHP curves were evaluated to represent the averaged interatomic interactions between Zn and Sb atoms either in a triangular environment of the [ZnSb₃] unit or in the tetrahedral coordination of the [ZnSb₄] unit [See Figure 5C and D]. In particular, the Cu2-Sb COHP curve in the quinary Ca₈YbZn₄Cu_0.5Sb₉ model was also displayed in Figure 5D. Overall, both of these types of Zn-Sb and Cu-Sb COHP curves were nicely optimized at E_F in both structural models, and the non-bonding characteristics continued even below E_F indicating the possibility of a certain amount of p-type doping on the anionic sites. However, above E_F, these two Zn-Sb interactions showed quite strong antibonding characteristics. In particular, the Zn1-Sb COHP curves displayed a three-fold stronger antibonding character than that of Zn2/3-Sb. This strongly indicated that the central site of the [ZnSb₃] triangle was supposed to be partially occupied as experimentally observed in the SXRD refinement results.

Figure 5. TDOS, PDOS and COHP curves of Ca₉Zn_4.5Sb₉ (A, C) and Ca₈YbZn₄Cu_0.5Sb₉ (B, D). The E_F corresponding to 72 and 71.5 valence electron count (VEC) in (A and B) is set as the reference energy (0 eV) and marked with a dashed vertical line. An additional line corresponding to 72 VEC in (B) is shown with a dotted line. In the COHP diagram, "+" regions indicate bonding interactions, while "-" regions represent antibonding interactions. Color codes in DOS curves are as follows: TDOS: bold black outline; Ca PDOS: gray region; Yb PDOS: magenta region; Zn TDOS: yellow region; Cu PDOS: green region; Sb PDOS: sky blue region.

The ELF calculations were conducted as well for these two structural models to analyze the different distributions of lone-pairs of electrons on the bridging Sb4 in two structural models. As shown in Figure 6, two lone-pairs of electrons on Sb4 merged together producing the “C-shaped” electron densities. In addition, due to the stronger electron withdrawing ability of Cu than Zn, these lone-pairs in the Cu-containing compound shown in Figure 6B were further spread out toward the Cu-substituents in Ca₈YbZn₄Cu_0.5Sb₉. In general, this kind of substitution can affect the bond polarity on the anion frameworks. Given the electronegativity difference among anionic elements (Zn = 1.65, Cu = 1.90, and Sb = 2.05 in the Pauling scale)^[36], as Cu was introduced, the bond polarity of the Cu-Sb became smaller than that of Zn-Sb. As the bond polarity decreases, the particular bond becomes relatively more flexible, which leads to the high doping efficiency and the corresponding high carrier concentration^[37]. This enhancement in carrier concentration contributes to improved electrical conductivity, which eventually improves the electrical conductivity based on σ = neμ, where σ = electrical conductivity, n = carrier concentration, e = the charge of an electron, μ = carrier mobility^[6]. Detailed discussion will be discussed in the following Thermoelectric Properties section.

Figure 6. ELF analysis illustrating the electron lone-pairs associated with the bridging Sb4 in (A) Ca₉Zn_4.5Sb₉ and (B) Ca₈YbZn₄Cu_0.5Sb₉. The 3D isosurfaces identify regions where the ELF values are higher than 0.59.

To understand the affinity of Yb for the A3 site, the QVAL of the five cationic sites in the ternary Ca₉Zn_4.5Sb₉ model was thoroughly evaluated by integrating electron densities within each corresponding WS sphere. As briefly mentioned in the previous section, in the Yb-containing quaternary compounds, Yb showed a discernible site-preference for the A3 site. This site-preference of any component is generally governed by two factors: (1) the size factor criterion taking into account the compatibility between the site volume and the size of the occupying element; and (2) the electronic factor criterion determined by the relationship between the intrinsic QVAL of each atomic site and the electronegativity of an occupying element. We thoroughly investigated those two factors to study the site preference of Yb and concluded that since Ca²⁺ and Yb²⁺ possess quite similar effective ionic radii (1.00 and 1.02 Å, respectively)^[35], any size-based preference could not be considered as a deciding factor. On the other hand, as presented in Table 4, the site with the highest QVAL of 2.336 is the A3 site. Therefore, this particular site should be preferred by the more electronegative Yb (Ca = 1.04, Yb = 1.06 on the Allred-Rochow scale)^[36].

Thermoelectric properties measurements

To investigate the effects of cationic substitution with divalent Yb and anionic p-type substitution with monovalent Cu, the temperature-dependent TE properties were measured across a temperature range of 323 to 861 K.

First, Figure 7A shows the temperature-dependent electrical conductivities σ of four title compounds. All compounds exhibited decreasing electrical conductivity σ with increasing temperature which is a characteristic typical of a heavily doped degenerate semiconductor. The measured σ values at 323 K were 243.4, 434.1, 542.8, 567.9, and 430.4 S/cm for Ca₉Zn_4.5Sb₉^[19], Ca₈YbZn_4.5Sb₉, Ca₈YbZn_4.45Cu_0.05Sb₉, Ca₈YbZn_4.4Cu_0.1Sb₉, and Ca₈YbZn_4.35Cu_0.15Sb₉, respectively.

Figure 7. For four compounds in the Ca₈YbZn_4.5-yCu_ySb₉ (0 ≤ y ≤ 0.15) system, temperature-dependent (A) electrical conductivity σ and (B) Seebeck coefficient S are plotted over the range of 323 to 861 K. Corresponding data for the ternary Ca₉Zn_4.5Sb₉^[19] compound are presented for comparison.

In particular, the σ values were higher in one quaternary compound with the Yb substitution and three quinary compounds with the Yb and Cu co-substitution compared to in the ternary parental compound Ca₉Zn_4.5Sb₉. This type of improvement of σ can be ascribed to the increased hole concentration in the quaternary and quinary compounds having the Yb and Cu substitutions as shown in Table 5.

The temperature-dependent Seebeck coefficients S of the title compounds are displayed in Figure 7B. The S values in all studied cases increased with temperature up to ca. 700 K, and then the values were saturated beyond 700 K. The maximum S values were 149.9, 163.2, 125.5, 144.5, and 147.2 μV/K for Ca₉Zn_4.5Sb₉^[19], Ca₈YbZn_4.5Sb₉, Ca₈YbZn_4.45Cu_0.05Sb₉, Ca₈YbZn_4.4Cu_0.1Sb₉, and Ca₈YbZn_4.35Cu_0.15Sb₉, respectively. The quaternary compound Ca₈YbZn_4.5Sb₉ showed the highest S, which typically showed the opposite trend to those of σ. In general, the correlation between S, the carrier concentration, and the effective mass can be expressed by = $$\frac{8\pi^{2} k_{B}^{2}}{3 e h^{2}} m^{*} T\left(\frac{\pi}{3 n}\right)^{\frac{2}{3}}$$, where k_B = the Boltzmann constant, e = the charge of electron, h = the Planck constant, n = the carrier concentration, and m^* = the effective mass^[6]. Therefore, according to this equation, the effective masses of our title compounds can be evaluated as follows: 0.74, 0.80, 0.82, and 0.96 m_e for Ca₉Zn_4.5Sb₉, Ca₈YbZn_4.5Sb₉, Ca₈YbZn_4.4Cu_0.1Sb₉, and Ca₈YbZn_4.35Cu_0.15Sb₉, respectively. In addition, the effective mass and mobility are inversely proportional to each other, but those hardly affect S directly. Therefore, this kind of observation implied the possibility of multiple electronic bands or nonparabolic band effects, which could influence the TE properties of the title compounds. In addition, Supplementary Figure 4A displays the results for the power factor (PF), which was computed using PF = S²σ.

Figure 8A shows the temperature-dependent thermal conductivity κ_tot of the title compounds with minimum κ_tot values of 0.57, 0.82, 0.76, 0.70, and 1.00 W/mK for Ca₉Zn_4.5Sb₉^[19], Ca₈YbZn_4.5Sb₉, Ca₈YbZn_4.45Cu_0.05Sb₉, Ca₈YbZn_4.4Cu_0.1Sb₉, and Ca₈YbZn_4.35Cu_0.15Sb₉, respectively. Typically, the κ_tot is delineated as the sum of the electronic contribution κ_elec and the lattice contribution κ_latt. Given that the κ_elec value is determined through the application of the Wiedemann-Franz law (κ_elec = LσT; $$L=1.5+\exp \left[-\frac{|S|}{116}\right]$$, where L is in 10^-8 WΩK and S in μ/K), the determination of the term κ_latt involves the subtraction of the κ_elec from the κ_tot value (κ_latt = κ_tot - κ_elec)^[32].

Figure 8. For four compounds in the Ca₈YbZn_4.5-yCu_ySb₉ (0 ≤ y ≤ 0.15) system, temperature-dependent (A) total thermal conductivity κ_tot and (B) figure-of-merit ZT are plotted over the range of 323 to 861 K. Corresponding data for the ternary Ca₉Zn_4.5Sb₉^[19] compound are presented for comparison.

The magnitude of κ_tot for the three co-substituted compounds exhibited decreasing patterns down to ca. 800 K. The expected decrease in κ_tot due to the increased phonon scattering caused by Yb substitution was contrary to expectations, as it increased compared to Ca₉Zn_4.5Sb₉. This should be attributed to the significantly increased κ_elec due to the enhancement in σ. Nevertheless, despite the increase in κ_elec, the κ_tot did not significantly increase in most of the Ca₈YbZn_4.5-yCu_ySb₉ compounds, as κ_latt was reduced by the enhanced phonon scattering due to the Yb and Cu co-substitution. Based on this correlation, we have computed the κ_elec of 0.16, 0.28, 0.36, 0.27, and 0.27 W/mK at 325 K for Ca₉Zn_4.5Sb₉^[19], Ca₈YbZn_4.5Sb₉, Ca₈YbZn_4.45Cu_0.05Sb₉, Ca₈YbZn_4.4Cu_0.1Sb₉, and Ca₈YbZn_4.35Cu_0.15Sb₉, respectively. The κ_latt, derived by subtracting κ_elec from κ_tot, is also presented in Supplementary Figure 4B.

Ultimately, the temperature-dependent ZT was calculated and displayed in Figure 8B. Overall, the ZT for all title compounds showed an increasing trend up to 861 K. Among these title compounds, Ca₈YbZn_4.4Cu_0.1Sb₉ attained a maximum ZT of 0.81 at 843 K, representing an improvement over the ternary parental compound Ca₉Zn_4.5Sb₉. In addition, the previously reported Cu-doped Ca₉Zn_4.35Cu_0.15Sb₉^[38] showed a ZT of 0.72 at 873 K. Therefore, the highest ZT of our quaternary Ca₈YbZn_4.4Cu_0.1Sb₉ should be attributed to the enhanced σ and S resulting from the co-substitution of Yb and Cu compared to Ca₉Zn_4.5Sb₉. While the total thermal conductivity of Ca₈YbZn_4.4Cu_0.1Sb₉ displayed an increase in comparison to Ca₉Zn_4.5Sb₉, the enhanced σ was found to possess a disproportionately large influence over the overall TE properties.