From synthesis to machine learning validation for the thermoelectric Zintl phase Ba_1−xEu_xZn_2−yAl_ySb₂ system

Abstract

Five new Zintl phase compounds in the Ba_1−xEu_xZn_2−yAl_ySb₂ system were synthesized using the molten Pb-flux method. All compounds crystallized in the orthorhombic BaCu₂S₂-type phase with a phase selectivity determined by the r⁺/r⁻ radius ratio criterion. Detailed structural analysis revealed a complex three-dimensional anionic framework of [(Zn/Al)Sb₄] tetrahedra combined with the Ba/Eu mixed-cations occupying the central voids within the frameworks. Density functional theory calculations confirmed structural stability and revealed changes in the electronic environment, resulting in reduced carrier mobility and ultra-low thermal conductivity. Thermoelectric property measurements showed that the enhanced Seebeck coefficients and the minimized thermal conductivities yielded a maximum ZT of 0.50 at 653 K. However, a central finding of this work is the elucidation of the complex chemical interplay governing the system's electronic properties. Despite intentional n-type substitution with Al, all compounds remained robustly p-type. We demonstrate this is a direct consequence of a chemical balance between intrinsic incomplete electron transfer and the competing electronic effects of Eu and Al substituents. Three machine learning models developed by XGBRegressor algorithms and trained using a customized dataset of 5798 thermoelectric real-world experimental outcomes demonstrated high predictive accuracy. These models not only matched the experimental trends but, more importantly, correctly predicted the persistent p-type behavior. This success demonstrates that the model learned the complex, real-world chemical physics, allowing it to make predictions that go beyond simple electron-counting rules. Herein, we established a powerful workflow for accelerating thermoelectric materials discovery by integrating experimental synthesis with density functional theory analysis and machine learning validation.

---

1. Introduction

Thermoelectric (TE) materials offer a sustainable solution to the global energy crisis by directly converting waste heat into useful electricity.¹ The conversion efficiency is governed by the dimensionless figure of merit, ZT = σS²T/κ, where σ is the electrical conductivity, S is the Seebeck coefficient, T is the absolute temperature, and κ is the thermal conductivity.² Achieving a high ZT value requires the difficult task of optimizing these inherently coupled transport properties—namely, maximizing the power factor (σS²) while simultaneously minimizing κ.³

Zintl phases have emerged as a highly promising class of TE materials, largely due to their intrinsic “phonon-glass electron-crystal” (PGEC) character. This has led to remarkable performance, with ZT values exceeding 2.0 in n-type Mg₃Sb₂-based systems.⁴ However, despite these breakthroughs, significant “bottlenecks” hinder widespread practical application. Many high-performing p-type Zintl compounds that achieve ZT > 1.0, such as those based on Yb or Cd, rely on elements that are either scarce or toxic, posing concerns for scalability and sustainability.^5–7 Furthermore, the field faces persistent challenges in material stability,⁸ device compatibility,⁹ and achieving effective doping.¹⁰ Crucially, while their complex structures are advantageous for achieving low thermal conductivity, they also create a fundamental challenge: the intricate coupling between thermal and electronic transport properties makes it difficult to suppress κ without also negatively impacting carrier mobility m*.¹¹ This difficulty is compounded because, owing to the complex crystal and electronic band structures, it is also not easy to understand the relationship between structures and properties in Zintl phases, which makes the optimization of TE properties even more challenging compared to those conventional semiconductors.¹¹ This landscape highlights a clear and pressing need for new p-type Zintl materials that are not only high-performing but are also composed of earth-abundant, benign elements, and for which novel optimization strategies can be rationally explored.

Furthermore, from a methodological standpoint, another gap has emerged. While artificial intelligence (AI) and machine learning (ML) methods are increasingly used to analyze complex data, the synergy between practice and validation remains challenging.^12,13 Many studies focus on using ML to predict novel compounds, but a critical gap still exists in using these powerful tools to interpret, validate, and extract physical insights from complex experimental results.¹⁴ To move the TE research field forward, it is essential to complete the synergistic loop not only by using predictions to guide experiments, but also by using experimental data to build robust ML models that explain and validate laboratory findings.

To address these concurrent challenges, we focus on the AM₂M′₂ Zintl family (where A = alkaline-earth/rare-earth, M = transition metal, M′ = pnictogen). These systems are known to draw special attention as they can crystallize in diverse and complex structures, such as the 2-dimensional (2D) layered CaAl₂Si₂-type or the 3-dimensional (3D) cage-like BaCu₂S₂-type, with the phase selection often being governed by factors like the ionic radii of the constituent elements.^15–20 Within this family, the BaZn₂Sb₂ system is a particularly noteworthy platform.²¹ Composed of earth-abundant and non-toxic elements, it is a compelling candidate to overcome the material sustainability bottleneck. The aforementioned radius ratio criterion was verified in it when we investigated the Ba_1−xSr_xZn_2−yCd_ySb_2, Ba_1−xSr_xZn_2−yCu_ySb₂, BaZn_2−xCd_xSb₂, and Ba_1−xEu_xZn₂Sb₂ systems over the past several years.^15,18–20 However, the performance of pristine BaZn₂Sb₂ is modest, limited by a non-optimal carrier concentration n and relatively high κ. Improving its performance requires complex multi-element doping, which generates detailed datasets that are difficult to interpret with conventional analysis alone. Therefore, this system provides a perfect opportunity not only to develop a new sustainable TE material, but also to demonstrate the power of using ML-based approaches to validate experimental findings and unravel complex structure–property relationships in a co-substituted Zintl phase.

Herein, we propose and investigate a synergistic co-substitution strategy in the BaZn₂Sb₂ system. Building on our previous work where Eu-substitution effectively scatters phonons, we introduce Al-substitution on the Zn-site.²⁰ As a donor impurity, Al is intended to optimize the p-type carrier concentration through compensation. However, BaZn₂Sb₂ possesses a robust intrinsic p-type character, which we attribute to the fundamental chemistry of Zintl phases, specifically, incomplete electron transfer from the cation to the anionic framework. This creates a complex scenario where the final electronic properties will be determined by the net outcome of these competing chemical influences. Our central goal was therefore to investigate this interplay, specifically to determine how the system balances its baseline p-type character against the electron-donating effect of Al and the electronegativity-driven influence of Eu.

In this work, we demonstrate the outcome of this investigation through a comprehensive study. Our results reveal that the final n is a predictable consequence of this chemical balance. We show that while Al-doping successfully reduces the hole concentration as intended, it is not sufficient to overcome the system's intrinsic p-type nature, which is further reinforced by the presence of the more electronegative Eu. This clear chemical model fully explains the persistent p-type behavior observed in all our compounds and accounts for the final ZT of ca. 0.5 at 653 K for the Ba_0.97Eu_0.03Zn_1.83Al_0.17Sb₂ composition. Crucially, the discrepancy between idealized Density functional theory (DFT) calculations (which predict n-type behavior) and our experimental results highlights the power of our ML models. By learning from real experimental data, the models correctly predicted the persistent p-type outcome, thereby validating our interpretation that the system is governed by complex chemical physics, not simple electron-counting rules.

2. Experimental section

2.1. Synthesis

The synthesis strategy was designed as a two-pronged approach to deliberately probe and modify the TE properties of the BaZn₂Sb₂ system. First, the molten Pb-flux method was chosen, as it is a well-established technique for growing the high-quality, low-defect single crystals necessary for accurate structural analysis. Second, the selection of Eu and Al as co-substituents was a targeted strategy. The isoelectronic substitution of Eu on the Ba site was based on our prior work, which established it as an effective method for reducing lattice thermal conductivity (κ_latt) via phonon scattering.²⁰ The aliovalent substitution of Al on the Zn site, explored here for the first time, was introduced as a classic n-type donor impurity with the scientific goal of investigating how it would interact with the system's intrinsic p-type character and influence the final carrier concentration.

Five title compounds in the Ba_1−xEu_xZn_2−yAl_ySb₂ (0 ≤ x ≤ 0.04(2); 0.07(2) ≤ y ≤ 0.39(2)) system were successfully synthesized by the molten Pb-flux method. All sample preparation processes were performed inside an Ar-filled glove box, where both H₂O and O₂ levels were below 0.1 ppm, or under a vacuum to avoid oxidation. The reactant elements were purchased from Alfa Aesar, and these included Ba (rod, 99+%), Eu (ingot, 99.9%), Zn (shot, 99.99%), Al (shot, 99.9%), Sb (shot, 99.9999%), and Pb (granules, 99.99%). The reactants were loaded in an alumina crucible with the excess amount of Pb-metals, and the ratios of these reactants which were systematically varied to achieve different final compositions were as follows: Ba : Zn : Al : Sb : Pb = 1 : 1.8 : 0.2 : 2 : 10, 1 : 1.7 : 0.3 : 2 : 10, 1 : 1.5 : 0.5 : 2 : 10 for BaZn_1.84(2)Al_0.16Sb₂, BaZn_1.75(2)Al_0.25Sb₂, BaZn_1.61(2)Al_0.39Sb₂, and Ba : Eu : Zn : Al : Sb : Pb = 0.8 : 0.2 : 1.9 : 0.1 : 2 : 10, 0.9 : 0.1 : 1.7 : 0.3 : 2 : 10 for Ba_0.96(2)Eu_0.04Zn_1.93(2)Al_0.07Sb₂, Ba_0.97(1)Eu_0.03Zn_1.83(1)Al_0.17Sb₂, respectively. The large excess of Pb, used as a metallic flux, is a standard practice to ensure complete dissolution of reactants and facilitate mass transport for high-quality crystal growth.

Each alumina crucible was enclosed in a fused-silica tube, with Al₂O₃ wool positioned at both the top and bottom and sealed under vacuum. The sealed ampoules were subjected to a carefully optimized multi-step thermal profile. They were first heated to 1273 K (120 K h⁻¹) and held for 24 h to ensure the formation of a homogeneous liquid solution. This was followed by slow cooling at 4 K h⁻¹ to 823 K, the critical step for nucleating and growing large single crystals. The ampoules were then held at 823 K for an extended period of 92 h to act as a crucial post-growth annealing step, allowing the crystals to reach thermodynamic equilibrium and minimizing point defects. Subsequently, the reaction containers were taken out of the furnace at 823 K, inverted, and instantaneously centrifuged for 3 min at 2400 rpm to separate the molten Pb metal from the needle-shaped crystal products (See SI Fig. S1). This entire synthesis protocol was held constant for all five compositions, with only the initial reactant loading ratios being changed. This represents a controlled experiment designed to isolate the effect of composition on the final structure.

These five title compounds were also prepared by ball-milling followed by hot-pressing to be subsequently exploited for the TE property measurement. First, a mixture of starting materials was loaded into a stainless-steel container with two 0.5-inch and two 0.25-inch stainless-steel balls inside an Ar-filled glove box. Then, the reactant mixture was ball-milled for 2 h 30 min using the SPEX 8000 M machine. During the process, this equipment was stopped every 30 min to scrape and remix the pulverized reactant mixture to ensure homogeneity. The resulting powder was then densified using a hot-pressing procedure. The powder was loaded into a 12.5 mm inner-diameter graphite die that had been pre-sprayed with a boron nitride (BN) coating to facilitate sample release. The loaded die was placed in the hot-press furnace, which was then evacuated and backfilled with high-purity argon three times. A uniaxial pressure of 60 MPa was applied, and the temperature was subsequently ramped to 823 K over a 2-hour period. The sample was held at 823 K and 60 MPa for a dwell time of 2 hours. After the dwell time, the furnace heater was turned off, the pressure was immediately released, and the sample was allowed to cool to room temperature. The final pellet was polished to remove the surface BN layer before its density and transport properties were measured. These final hot-pressing parameters were determined through a systematic optimization study. It was found that while a pressure of 40 MPa resulted in insufficient densification (ca. 80% relative density), increasing the pressure to 50 MPa and above consistently yielded dense pellets (>90% relative density). A pressure of 60 MPa was therefore selected to ensure maximum densification. The hot-pressing temperature was also optimized. It was observed that increasing the temperature from 823 K to 873 K did not lead to any significant improvement in the phase purity or crystallinity of the samples, as confirmed by powder X-ray diffraction (results not shown here). Given that thermal analysis (discussed later) indicates the material begins to degrade at temperatures near 900 K, the more conservative temperature of 823 K was chosen as it provides no disadvantage in sample quality while ensuring a wider, safer processing window. Consequently, the conditions of 823 K and 60 MPa were established as the optimal parameters to reliably produce dense, phase-pure samples for all compounds reported in this study.

2.2. Material characterization

The phase purities and detailed crystal structures of the isotypic five title compounds in the Ba_1−xEu_xZn_2−yAl_ySb₂ system were evaluated using powder and single-crystal X-ray diffraction (PXRD and SXRD) data analyses. First, the PXRD patterns of the samples were collected by the Bruker D8 diffractometer with monochromatic Cu Kα₁ radiation (λ = 1.54059 Å) in the range of 20° ≤ 2θ ≤ 80° with a step size of 0.02°, and the results are plotted in Fig. 1. In addition, the Rietveld refinement was conducted for the collected PXRD patterns of five title compounds using PDXL2 software to obtain the lattice parameters as well as volumes of the unit cells.²² The refinement results and lattice parameters are provided in SI Table S1 and Fig. S2. Second, the SXRD data of the three out of five title compounds: BaZn_1.75(2)Al_0.25Sb₂, Ba_0.96(2)Eu_0.04Zn_1.93(2)Al_0.07Sb₂, and Ba_0.97(1)Eu_0.03Zn_1.83(1)Al_0.17Sb₂, were collected at room temperature by using the Bruker D8 QUEST diffractometer equipped with Mo Kα₁ radiation (λ = 0.71073 Å). A full data collection was conducted for the best crystal using the APEX3 program.²³ The SAINT program was used for unit cell parameter refinement,²⁴ and the SADABS program was used for semiempirical absorption corrections based on equivalents.²⁵ The SXRD data of the two out of five title compounds: BaZn_1.84(2)Al_0.16Sb₂ and BaZn_1.61(2)Al_0.39Sb₂, were collected using synchrotron radiation (λ = 0.65000 Å) on an MX225HS detector (Rayonix, Evanston, IL, USA) at BL2D SMC in the Pohang Accelerator Laboratory. Data collection was conducted using the PAL BL2D-SMDC program,²⁶ and cell refinement, reduction, and absorption correction were performed using the HKL3000sm (version 715)²⁷ program. The detailed crystal structures including the atomic positions accompanied by anisotropic displacement parameters (ADPs), the extinction coefficients, and the mixed ratios of Ba and Eu, and Zn and Al, were refined to convergence using full-matrix least-squares methods on F². More detailed crystallographic data including atomic positions and some selected bond distances are presented in Tables 1–3. The depository numbers are as follows: CCDC-2470329 for BaZn_1.84(2)Al_0.16Sb_2, CCDC-2471602 for BaZn_1.75(2)Al_0.25Sb₂, CCDC-2470330 for BaZn_1.61(2)Al_0.39Sb₂, CCDC-2470328 for Ba_0.96(2)Eu_0.04Zn_1.93(2)Al_0.07Sb₂ and CCDC-2470327 for Ba_0.97(1)Eu_0.03Zn_1.83(1)Al_0.17Sb₂.

Fig. 1 (a)–(e) Collected PXRD patterns of five title compounds in the Ba_1−xEu_xZn_2−yAl_ySb₂ (0 ≤ x ≤ 0.04(2); 0.07(2) ≤ y ≤ 0.39(2)) system. For comparison purposes, (f) a simulated PXRD pattern of Ba_0.97(1)Eu_0.03Zn_1.83(1)Al_0.17Sb₂ is provided as a reference. An asterisk indicates a peak from the trace of remaining Pb metals used as inactive metal flux.

Table 1 SXRD data and structure refinement results for the Ba_1−xEu_xZn_2−yAl_ySb₂ (0 ≤ x ≤ 0.04(2); 0.07(2) ≤ y ≤ 0.39(2)) system

Empirical formula	BaZn1.84(2)Al0.16Sb2	BaZn1.75(2)Al0.25Sb2	BaZn1.61(2)Al0.39Sb2	Ba0.96(2)Eu0.04Zn1.93(2)Al0.07Sb2	Ba0.97(1)Eu0.03Zn1.83(1)Al0.17Sb2
a R1 = Σ||Fo| − |Fc||/Σ|Fo|; wR2 = {Σ[w(Fo^2 − Fc^2)/Σ[w(Fo^2)^2]]^1/2, where w = 1/[σ^2Fo^2 + (A − P)^2 + B − P], in which P = (Fo^2 + 2Fc^2)/3 and A and B are weight coefficients.
Crystal system	Orthorhombic
Space group	Pnma (no. 62)
Unit cell dimensions (Å)	a = 10.539(2)	a = 10.536(2)	a = 10.535(2)	a = 10.592(6)	a = 10.549(2)
	b = 4.504(1)	b = 4.503(1)	b = 4.504(1)	b = 4.505(2)	b = 4.509(2)
	c = 11.633(2)	c = 11.616(2)	c = 11.615(2)	c = 11.668(6)	c = 11.660(3)
Volume (Å^3)	552.2(2)	551.1(2)	551.1(3)	556.8(5)	554.6(2)
dcalcd (g cm^−3)	6.081	6.049	5.986	6.081	6.053
Data/restraints/parameters	686/0/33	1299/0/34	866/0/33	871/0/35	784/0/35
R indices^a (I > 2σ(I))	R1 = 0.0290	R1 = 0.0222	R1 = 0.0280	R1 = 0.0162	R1 = 0.0161
	wR2 = 0.0758	wR2 = 0.0428	wR2 = 0.0692	wR2 0.0329	wR2 = 0.0342
R indices^a (all data)	R1 = 0.0300	R1 = 0.0282	R1 = 0.0300	R1 = 0.0183	R1 = 0.0201
	wR2 = 0.0765	wR2 = 0.0437	wR2 = 0.0701	wR2 = 0.0334	wR2 = 0.0355
Goodness of fit on F^2	1.027	1.307	1.004	1.198	1.074
Largest diff. Peak/hole (e Å^−3)	2.303/−1.844	1.232/−2.019	2.139/−2.504	0.698/−1.136	0.840/−0.855

---

Table 2 Atomic coordinates and equivalent isotropic displacement parameters (U_eq^a) from the SXRD refinements for the Ba_1−xEu_xZn_2−yAl_ySb₂ (0 ≤ x ≤ 0.04(2); 0.07(2) ≤ y ≤ 0.39(2)) system

Atom	Wyckoff site	Occupation	x	y	z	Ueq^a (Å^2)
a Ueq is defined as one-third of the trace of the orthogonalized Uij tensor.b A = Ba for the three quaternary compounds; Ba/Eu mixed-site for the two quinary compounds.c M1 = Zn1/Al1 mixed-site.d M2 = Zn2/Al2 mixed-site.
BaZn1.84(2)Al0.16Sb2
A^b	4c	1	0.2461(1)	1/4	0.3209(1)	0.0070(2)
M1^c	4c	0.93(1)/0.07	0.0538(2)	1/4	0.6174(1)	0.0070(3)
M2^d	4c	0.91(2)/0.09	0.0934(1)	1/4	0.0484(2)	0.0072(3)
Sb1	4c	1	0.3462(1)	1/4	0.0354(1)	0.0026(2)
Sb2	4c	1	0.4759(1)	1/4	0.6641(1)	0.0036(3)
BaZn1.75(2)Al0.25Sb2
A^b	4c	1	0.24634(3)	1/4	0.3211(1)	0.0113(2)
M1^c	4c	0.90(1)/0.10	0.05335(6)	1/4	0.6176(1)	0.0118(3)
M2^d	4c	0.84(1)/0.16	0.09302(6)	1/4	0.0486(1)	0.0114(3)
Sb1	4c	1	0.34601(3)	1/4	0.0357(1)	0.0064(1)
Sb2	4c	1	0.47574(3)	1/4	0.6639(1)	0.0077(2)
BaZn1.61(2)Al0.39Sb2
A^b	4c	1	0.2469(1)	1/4	0.3209(1)	0.0120(2)
M1^c	4c	0.85(1)/0.15	0.0530(2)	1/4	0.6176(1)	0.0122(3)
M2^d	4c	0.76(1)/0.24	0.0921(2)	1/4	0.0486(2)	0.0118(3)
Sb1	4c	1	0.3454(1)	1/4	0.0354(1)	0.0073(2)
Sb2	4c	1	0.4754(1)	1/4	0.6641(1)	0.0085(2)
Ba0.96(2)Eu0.04Zn1.93(2)Al0.07Sb2
A^b	4c	0.96(2)/0.04	0.2461(1)	1/4	0.3202(1)	0.0124(2)
M1^c	4c	0.97(1)/0.03	0.0546(1)	1/4	0.6169(1)	0.0128(2)
M2^d	4c	0.97(1)/0.03	0.0931(1)	1/4	0.0480(1)	0.0127(2)
Sb1	4c	1	0.3459(1)	1/4	0.0348(1)	0.0082(2)
Sb2	4c	1	0.4760(1)	1/4	0.6648(1)	0.0089(1)
Ba0.97(1)Eu0.03Zn1.83(1)Al0.17Sb2
A^b	4c	0.97(1)/0.03	0.2462(1)	1/4	0.3208(1)	0.0155(2)
M1^c	4c	0.93(1)/0.07	0.0537(1)	1/4	0.6173(1)	0.0155(2)
M2^d	4c	0.90(1)/0.10	0.0932(1)	1/4	0.0483(1)	0.0155(2)
Sb1	4c	1	0.3460(1)	1/4	0.0354(1)	0.0107(2)
Sb2	4c	1	0.4759(1)	1/4	0.6642(1)	0.0117(2)

---

Table 3 Selected bond distances for the Ba_1−xEu_xZn_2−yAl_ySb₂ (0 ≤ x ≤ 0.04(2); 0.07(2) ≤ y ≤ 0.39(2)) system

Atomic pair	Bond distance (Å)
	BaZn1.84(2)Al0.16Sb2	BaZn1.75(2)Al0.25Sb2	BaZn1.61(2)Al0.39Sb2	Ba0.96(2)Eu0.04Zn1.93(2)Al0.07Sb2	Ba0.97(1)Eu0.03Zn1.83(1)Al0.17Sb2
a A = Ba for the three quaternary compounds; Ba/Eu mixed-site for the two quinary compounds.b M1 = Zn1/Al1 mixed site.c M2 = Zn2/Al2 mixed site.
A^a−M1^b (×2)	3.892(2)	3.888(1)	3.894(2)	3.894(2)	3.895(1)
A^a−M2^c (shorter)	3.558(2)	3.554(2)	3.566(3)	3.565(2)	3.563(2)
A^a−M2^c (longer, ×2)	3.863(2)	3.862(1)	3.872(2)	3.878(2)	3.872(1)
A^a−Sb1 (shorter)	3.486(2)	3.478(1)	3.483(2)	3.494(2)	3.491(2)
A^a−Sb1 (longer, ×2)	3.499(1)	3.497(1)	3.506(2)	3.506(2)	3.505(2)
A^a−Sb2 (shorter, ×2)	3.700(2)	3.698(1)	3.703(2)	3.711(2)	3.703(2)
A^a−Sb2 (longer, ×2)	3.725(1)	3.726(1)	3.733(2)	3.728(2)	3.729(1)
M1^b−Sb1 (shorter, ×2)	2.663(1)	2.664(1)	2.674(2)	2.665(1)	2.667(1)
M1^b−Sb1 (longer)	2.817(2)	2.818(1)	2.824(2)	2.832(2)	2.823(2)
M1^b−Sb2	2.672(2)	2.666(2)	2.671(2)	2.680(2)	2.676(2)
M2^c−Sb1	2.669(2)	2.670(2)	2.675(2)	2.681(2)	2.671(1)
M2^c−Sb2 (shorter, ×2)	2.724(2)	2.719(1)	2.722(2)	2.733(2)	2.727(1)
M2^c−Sb2 (longer)	2.763(2)	2.760(2)	2.766(2)	2.776(2)	2.771(1)

---

Elemental distributions were analyzed by energy-dispersive X-ray spectroscopy (EDS) using an ULTRA Plus field-emission (Carl Zeiss, Oberkochen, Germany) scanning electron microscope (SEM) system with an acceleration voltage of 30 kV. Detailed EDS analysis and mapping results are provided in SI Fig. S3–S7. Thermogravimetric analysis (TGA) was performed using a STA200 analyzer (HITACHI, NEXTA STA200, Japan) with a heating rate of 10 K min⁻¹ under continuous N₂ flow conditions (see SI Fig. S8).

2.3. Electronic structure calculations

To understand the overall electronic structure of the title compounds, a series of DFT calculations was conducted by the tight-binding linear muffin-tin orbital (TB-LMTO) method with the atomic sphere approximation (ASA)^28–30 using a total of three structural models (two hypothetical models and one reference model) with the idealized compositions of BaZn_1.5Al_0.5Sb₂, Ba_0.75Eu_0.25Zn_1.5Al_0.5Sb₂, and BaZn₂Sb₂. The resultant Density of states (DOS), crystal orbital Hamilton population (COHP) curves, band structures, and electron localization function (ELF) diagrams of these models were thoroughly studied. The symmetry of the two hypothetical models was lowered from the experimentally obtained space group Pnma (no. 62) to its subgroup P1m1 (no. 6) to accommodate the given chemical composition. The lattice parameters and the atomic positions of the two hypothetical models were extracted from the SXRD refinement results of BaZn_1.61(2)Al_0.39Sb₂ and Ba_0.97(1)Eu_0.03Zn_1.83(1)Al_0.17Sb₂, respectively. The structural details of these models are provided in SI Table S2. All the space in a unit cell is filled with overlapping Wigner–Seitz (WS) atomic spheres, and all relativistic effects (except spin–orbit coupling) were considered using a scalar relativistic approximation.³¹ Each WS sphere is regarded as containing a symmetrically spherical potential, and the combined correction was applied for the overlapping regions. The radii of each WS sphere were calculated automatically to make the overlapping potential the best approximation to the full potential.³¹ The WS radii used for the three models are as follows: Ba, 2.468 Å; Zn, 1.494–1.540 Å; Sb, 1.738–1.782 Å for BaZn₂Sb₂; Ba, 2.455–2.477 Å; Zn, 1.499–1.505 Å; Al, 1.560 Å; Sb, 1.718–1.781 Å for BaZn_1.5Al_0.5Sb₂; and Ba, 2.473–2.489 Å; Eu, 2.489 Å; Zn, 1.494–1.499 Å; Al, 1.551 Å; Sb, 1.721–1.784 Å for Ba_0.75Eu_0.25Zn_1.5Al_0.5Sb₂. The basis sets included 6s, 6p, 5d, and 4f orbitals for Ba; 6s, 6p, and 5d orbitals for Eu; 4s, 4p, and 3d orbitals for Zn; 3s, 3p and 3d orbitals for Al; and 5s, 5p, 5d, and 4f orbitals for Sb. The Löwdin downfolding technique was applied for the Ba 6p, Eu 6p, Al 3d, Sb 5d, and 4f orbitals.³² The self-consistent charge density was obtained using 832 irreducible k-points in the Brillouin zone for all three models.³³

2.4. Thermoelectric property measurements

The density (ρ) of the pellet used for physical property measurements was determined using the geometric method on the dense, regularly shaped hot-pressed pellets. The relative density of all samples was confirmed to be higher than 90% of the theoretical density calculated from the SXRD data, indicating successful densification.

The electrical conductivity σ and the Seebeck coefficient S were simultaneously measured between room temperature and 873 K under a He atmosphere by a ZEM-3 instrument system (ULVAC-RIKO Inc., Yokohama, Japan). To ensure data accuracy, the instrument's calibration is periodically verified against a constantan standard reference material, with measured values consistently falling within a < 7% error tolerance. Furthermore, before each measurement, the system performs an automated V–I plot linearity check to ensure a high-quality ohmic contact between the probes and the sample. Temperature-dependent measurements were only carried out when the correlation coefficient R² exceeded 0.9999.

Finally, thermal diffusivity (D) was measured from room temperature to 873 K under an N₂ atmosphere using a Netzsch LFA 457 HyperFlash instrument (NetzschGerätebau GmbH, Selb, Germany). The accuracy of this instrument is confirmed by standard reference materials (Pyroceram 9606, Poco Graphite) to be within ±3% for diffusivity and ±5% for heat capacity. Then, the total thermal conductivity κ_tot was calculated using the equation: κ_tot = DCpρ (D = thermal diffusivity, Cp = heat capacity, and ρ = density).³⁴ This work exploited the Dulong–Petit value (3R/atom, R = gas constant) for Cp. The n and carrier mobility µ were determined using a van der Pauw four-point probe configuration (HMS-3000, Ecopia) at room temperature with an applied magnetic field of 0.55 T.³⁵

2.4. Machine learning

The publicly available experimentally synthesized TE materials (ESTM)¹⁴ dataset was combined with our in-house AM₂M′₂ (A = alkaline-earth/rare-earth-metals, M = transition metals, M′ = pnictogens) Zintl phase dataset for training the three ML models. A composition-based feature vector (CBFV) generation method was used to generate feature vectors from the chemical compositions presented in our dataset.^36,37 A total of 154 feature vectors were extracted from the materials agnostic platform for informatics and exploration (MAGPIE) descriptors.³⁸ Furthermore, our customized dataset was divided into training and testing sets, where 80% of the data was assigned for training purposes and 20% was reserved for testing. Then, the input data was scaled and normalized using scikit-learn's built-in standardScaler class and normalize function, respectively. The XGBRegressor ensemble-learning algorithm was used to predict the three TE properties including σ, S, and κ_tot values, and the feature importance function of the algorithm was employed to compute the significance of the individual features using an importance score.³⁹ The hyperparameters, generated using the Bayesian optimization technique from the specialized Python package optuna for the three independent regression models are listed in SI Table S3. Finally, the five-fold cross-validation method was used to assess the accuracy of the models using the coefficient of determination (R²) as the evaluation metric.³⁶ Further details regarding the ML procedures are described in the SI.

3. Results and discussion

3.1. Crystal structure analysis

Five new Zintl phase solid-solutions belonging to the Ba_1−xEu_xZn_2−yAl_ySb₂ (0 ≤ x ≤ 0.04(2); 0.07(2) ≤ y ≤ 0.39(2)) system have been successfully synthesized by the molten Pb-flux method for the TE materials application purpose. The initial assessment of these synthesized products, including phase analysis and purity check, was processed by indexing the collected PXRD peaks with the simulated PXRD peaks based on the SXRD refinement result for Ba_0.97(1)Eu_0.03Zn_1.83(1)Al_0.17Sb₂. As shown in Fig. 1(a)–(f), nearly all the collected peaks nicely matched the simulated peaks, except a small peak at 2θ = 36° descended from a trace of Pb-metal flux. Then, the definitive crystal structure analysis, including phase identification, substituent solid solubility, atomic site occupancies, unit cell parameters, and atomic coordinates inside the unit cell, was determined through full structural refinement of SXRD data, which is the most powerful technique for this purpose. The results proved that all five title compounds adopted the orthorhombic BaCu₂S₂-type phase (Pnma, Z = 4, Pearson code oP20), where the five independent atomic sites included one Ba or Ba/Eu mixed site, two Zn/Al mixed sites, and two Sb sites.

The SXRD refinement provides direct, quantitative proof of substituent incorporation and is the foundation of our structure–property analysis. It confirmed mixed occupancy on both the cation and anion sites. Interestingly, the Al-substituent was successfully introduced to both Zn sites with a slightly higher site-preference to the Zn2 site. The refinement results of SXRD data, atomic coordinates, and some selected bond distances are presented in Tables 1–3.

The additional structural assessments, including the phase purity and the lattice parameters, were performed through Rietveld refinement of the PXRD data. As displayed Table S1 and Fig. S2, all five title compounds were proven to be single-phase products.

A detailed structural description of the BaCu₂S₂-type five isotypic title compounds can be found in our previous reports.¹⁵ Therefore, we will just briefly illustrate the overall crystal structure using a representative quinary compound, Ba_0.97(1)Eu_0.03Zn_1.83(1)Al_0.17Sb₂. The overall structure shown in Fig. 2(a) can be described as a combination of (1) the 3D “cage-shaped” ³_∞[(Zn/Al)Sb3/4Sb1/4] anionic frameworks built by the 1D vertex-sharing anionic chains (Fig. 2(b)) and (2) the Ba/Eu mixed cationic sites filling the central voids of these cage-shaped anionic frameworks. In particular, an individual coordination environment of the mixed-cationic site can be seen as surrounded by nine Zn and seven Sb, eventually resulting in forming the 3D cage-shaped coordinates (Fig. 2(c)). Moreover, all five title compounds satisfied the aforementioned r₊/r₋ ratio criterion, which was proven to be a critical factor in terms of the phase selectivity between the BaCu₂S₂-type and the CaAl₂Si₂-type phases: r₊/r₋ = 1.04 and 1.03 for the two quaternary and two quinary title compounds, respectively. Therefore, the observed BaCu₂S₂-type phase should readily be taken by all title compounds.

Fig. 2 (a) Overall crystal structure of Ba_0.97(1)Eu_0.04Zn_1.83(1)Al_0.17Sb₂ illustrated by a combination of ball-and-stick and polyhedral representations. The tetrahedral [(Zn/Al)Sb₄] moiety is highlighted in light pink polyhedra. (b) Two-vertex–sharing 1D anionic chains building up the 3D infinite ³_∞[(Zn/Al)Sb3/4Sb1/4] framework and (c) the Ba/Eu mixed-cationic site surrounded by 16 anionic elements are also illustrated. Atomic labels are provided, and a unit cell is outlined in a black line. Color codes: Ba/Eu mixed-site, gray; Zn/Al mixed-site, pink; Sb, blue.

Finally, in addition to diffraction analysis, the micro-scale homogeneity and elemental composition of the synthesized crystals were confirmed using SEM coupled with EDS. The results, provided in the SI (Fig. S1 and S3–S7), show a uniform distribution of all constituent elements, corroborating the successful formation of a homogeneous solid solution at the microscopic level.

3.2. Electronic structure analysis

An influence of the Eu and Al co-substitution on electronic structures of the Ba_1−xEu_xZn_2−yAl_ySb₂ system was investigated based on the DFT calculation results performed using three model structures: (1) two hypothetical compounds with idealized compositions of BaZn_1.5Al_0.5Sb₂ and Ba_0.75Eu_0.25Zn_1.5Al_0.5Sb₂, and (2) one ternary reference compound BaZn₂Sb₂. The structural details of these models were described in the earlier Experimental section, and DOS and COHP curves of BaZn₂Sb₂ are provided in SI Fig. S9(a) and S9(b). By calculating the properties of perfect, idealized structures, we create a theoretical baseline that shows the “pure” effect of substitution on its own. We then compare this ideal result to our real-world experiments to help separate and understand all the complex factors involved.

First, the total and partial DOS (TDOS and PDOS) curves of BaZn_1.5Al_0.5Sb₂ and Ba_0.75Eu_0.25Zn_1.5Al_0.5Sb₂ models showed quite a strong orbital mixing throughout the entire energy range, as shown in Fig. 3(a) and (b). In particular, we observed shifts of the Fermi levels (E_F) toward a relatively higher energy region, where a steep DOS level increase was observed in both the Al-substituted models, as compared to that of the reference BaZn₂Sb₂, which was located near a pseudogap (SI Fig. S9(a)). This kind of E_F shift in two Al-containing title compounds should be attributed to the addition of extra electrons via the n-type Al-substituent, and the location of these shifted E_F corresponded to the valence electron count (VEC) of 16.5. However, the experimentally obtained VECs, corresponding to the SXRD refinement results of Ba_0.97(1)Eu_0.03Zn_1.83(1)Al_0.17Sb₂ and BaZn_1.61(2)Al_0.39Sb₂, turned out to be 16.17 and 16.39, respectively, which were smaller than 16.5. As a result, the DOS level at E_F of these two refined compositions also decreased compared to that of the ternary reference. In other words, the successfully synthesized title compounds tended to keep the Al substituent lower than 0.5 in order to maintain the low DOS levels for their structural stability in the BaZn₂Sb₂-type phase.

Fig. 3 TDOS, PDOS, and COHP curves for BaZn_1.5Al_0.5Sb₂ ((a) and (c)) and Ba_0.75Eu_0.25Zn_1.5Al_0.5Sb₂ ((b) and (d)). In the COHP diagram, the “+” and “−” regions indicate bonding and antibonding interactions, respectively. The dashed and dotted lines correspond to the two idealized compositions and the two SXRD refined compositions, respectively. The VECs for the corresponding compositions are also provided. Color codes: TDOS, bold black outline; Ba PDOS, gray region; Eu PDOS, purple region; Zn PDOS, pink region; Al PDOS, yellow region; Sb PDOS, blue region.

However, this theoretical prediction of the E_F shifting towards the conduction band appears to directly contradict our experimental transport measurements (discussed in Section 3.3), where all synthesized samples remain p-type conductors. This apparent discrepancy is, in fact, a central scientific finding of this work. The ideal DFT model in Fig. 3 correctly shows the intended electronic effect of Al substitution in a perfect, defect-free lattice. However, the experimental reality is governed by the real-world chemistry of the Zintl phase. As we will demonstrate with our Hall effect data in the following section, the final carrier type is a net balance between several competing factors: (1) the baseline p-type character of the parent compound from incomplete electron transfer, (2) the hole compensation from the Al-donor, and (3) an increase in holes caused by the higher electronegativity of Eu. This chemical model, rather than one based on thermodynamic self-compensation via native vacancies,⁴⁰ provides a direct, evidence-based explanation for why the samples remain p-type.

The COHP curves analysis can also rationalize this kind of reduced VECs from the experimentally obtained compositions. As shown in Fig. 3(d), the strong antibonding character of the Zn/Al–Sb COHP curve at E_F in the Ba_0.75Eu_0.25Zn_1.5Al_0.5Sb₂ model (16.5 VEC) significantly decreased at E_F of the experimentally refined Ba_0.97(1)Eu_0.03Zn_1.83(1)Al_0.17Sb₂ (16.17 VEC) due to the smaller amount of Al-substituent. In addition, E_F in BaZn_1.61(2)Al_0.39Sb₂ (16.39 VEC) was also slightly moved away from the shifted E_F to reduce the antibonding character of the Zn/Al–Sb COHP curve, as shown in Fig. 3(c). Furthermore, a strong bonding character of the Ba(Ba/Eu)-Zn/Al COHP curve successfully compensated for the remaining antibonding characters of the Zn/Al–Sb COHP curve in each model. Therefore, based on these comprehensive DOS and COHP curve analyses, we can claim that although some additional electrons were introduced to the title system by the n-type Al substitution, these title compounds found a way to stabilize their crystal structures by minimizing the occupation of the high DOS level and the strong antibonding states of the Zn/Al–Sb COHP curves.

Last, the changes in the distribution of paired electron density due to the co-substitution in the title system were thoroughly studied by a series of ELF calculations.⁴¹ The resultant ELF diagrams of the ternary BaZn₂Sb₂ and the quinary Ba_0.75Eu_0.25Zn_1.5Al_0.5Sb₂ are displayed in Fig. 4(a), (b), and SI Fig. S10. Since the electronegativity of the substituting Al is lower than Zn (Al = 1.47, Zn = 1.66, Sb = 1.82 in the Allred-Rochow scale),⁴² a polarity of the Al–Sb bond, which was newly formed in the two Al-containing compounds, was expected to increase compared to that of the Zn–Sb bond in the ternary compound. This kind of polarity change clearly resulted in the distortion of isosurfaces as shown in Fig. 4(b), in which the size of the isosurfaces surrounding Al2 shrank, while those surrounding Sb1 and Sb2 expanded considerably. Moreover, the size of the isosurface around Eu also significantly decreased in the quinary model compared to that around Ba in the ternary model. Therefore, we conclude that the Al-substitution can cause the changes in polarity of chemical bonding within the anionic frameworks, and with the addition of Eu-substitution causing atomic disordering, the currently executed co-substitution can influence the carrier mobility in the title system. In other words, this increased bond polarity leads to greater electron localization, which hinders charge carrier movement. This is one of the microscopic mechanisms responsible for the observed reduction in µ in the Al-substituted compounds. This type of correlation between the polarity modification and the resultant mobility change in Zintl phase TE materials has been reported in several previous articles.^15,43 We will further discuss the electrical transport properties in the subsequent thermoelectric properties measurements section.

Fig. 4 Two ELF diagrams for the unit cells of (a) BaZn₂Sb₂ and (b) Ba_0.75Eu_0.25Zn_1.5Al_0.5Sb₂. The 2D contour maps of the (0 1 0) sliced plane and the superimposed 3D isosurfaces are illustrated together. Each unit cell is outlined by white lines. A color scheme from blue to red for the 2D contour map indicates ELF values between 0.0 and 1.0, and the yellow 3D isosurfaces show the area having ELF values of 0.61.

3.3. Thermoelectric properties measurements

A series of TE properties, including σ, S, and κ, were measured for the five title compounds in the title Ba_1−xEu_xZn_2−yAl_ySb₂ system over the temperature range between 300 and 773 K. The n and µ of the selected samples were also measured for an in-depth interpretation of TE transport properties as the composition varied.

First, the temperature-dependent σ of the three quaternary and two quinary compounds displayed lower magnitudes than that of the ternary BaZn₂Sb₂, as shown in Fig. 5(a). This phenomenon, which was both intended and expected, should be attributed to the trivalent Al substitution for a divalent Zn, in which an extra electron was donated to the title system, resulting in a decrease of n and consequently a decrease of σ. In addition, the Al substitution also enhanced carrier scattering and bond polarity between anionic elements, eventually reducing µ due to the increased atomic disordering. In particular, boosting bond polarity is in agreement with the ELF result discussed earlier. On the other hand, an additional Eu substitution for Ba in the two quinary compounds slightly increased the overall σ values compared to those of the three Al-substituted quaternary compounds. This result can also partly be attributed to the relatively higher electronegativity of Eu than Ba (Eu = 1.01, Ba = 0.97 in the Allred-Rochow scale).⁴² Since Eu was relatively reluctant to donate electrons to the anionic framework compared to Ba, there existed somewhat increased n in the anionic frameworks, eventually enhancing σ magnitudes.

Fig. 5 Temperature-dependent (a) electrical conductivity σ and (b) Seebeck coefficient S of the five title compounds in the Ba_1−xEu_xZn_2−yAl_ySb₂ (0 ≤ x ≤ 0.04(2), 0.07(2) ≤ y ≤ 0.39(2)) system measured between 300 and 773 K.

As provided in Table 4, the results of the Hall measurement clearly indicated that as Al was introduced to the parental BaZn₂Sb₂, both n and µ of BaZn_1.84(2)Al_0.16Sb₂ decreased, causing a substantial reduction in σ value. However, a further addition of Eu to this quaternary compound somewhat improved n and µ, eventually resulting in a moderately increased σ in Ba_0.97(1)Eu_0.93Zn_1.83(1)Al_0.17Sb₂. This observation was also consistent with the experimentally measured σ values shown in Fig. 5(a). Moreover, the three quaternary compounds displayed the semiconducting characteristics, increasing σ with temperature, whereas the two Eu and Al co-substituted quinary compound varied with the relative Eu and Al content. To be more precise, Ba_0.96(2)Eu_0.04Zn_1.93(2)Al_0.07Sb₂ with the lower Al content showed a gradual decreasing σ trend with temperature, followed by a quite subtle increase, which was a typical behavior of the heavily doped semiconductor. On the other hand, the other quinary compound Ba_0.97(1)Eu_0.03Zn_1.83(1)Al_0.17Sb₂, with the relatively higher Al content, showed a typical semiconducting characteristic. Lastly, a substantial rise of σ observed at higher temperatures could be attributed to the intrinsic carrier excitation resulting in the increase of n. The observed maximum σ values were 117, 68, 100, 98, and 126 S cm⁻¹ for Ba_1.84(2)Zn_1.84(2)Al_0.16Sb₂, BaZn_1.75(2)Al_0.25Sb₂, BaZn_1.61(2)Al_0.39Sb₂, Ba_0.96(2)Eu_0.04Zn_1.93(2)Al_0.07Sb₂, and Ba_0.97(1)Eu_0.03Zn_1.83(1)Al_0.17Sb₂, respectively.

Table 4 Carrier concentration n and mobility µ of the selected compounds in the Ba_1−xEu_xZn_2−yAl_ySb₂ (0 ≤ x ≤ 0.04(1); 0 ≤ y ≤ 0.17(1)) system

Compound	Carrier concentration n (×10^19 carrier cm^−3)	Mobility µ (cm^2 V^−1s^−1)
BaZn2Sb2 (ref. 44)	3.2	54.0
Ba0.96(1)Eu0.04Zn2Sb2 (ref. 20)	6.6	30.1
BaZn1.84(2)Al0.16Sb2	2.2	5.9
Ba0.97(1)Eu0.03Zn1.83(1)Al0.17Sb2	3.4	15.2

---

Second, all the temperature-dependent S of the title compounds showed positive values and confirmed the p-type characteristics as displayed in Fig. 5(b). This persistence of p-type behavior is the crucial experimental confirmation of the chemical interplay model outlined in Section 3.2. The system does not undergo a p-to-n transition because the electron donation from Al, while significant, is simply insufficient to overcome the combined p-type influence of the parent compound's intrinsic hole concentration and the electronegativity-driven effect of Eu. This demonstrates a predictable chemical balance, not a pinning of the E_F by acceptor like defects as suggested previously by Pomrehn et al.⁴⁰ Overall, the S values increased with temperature up to a specific critical temperature, and then slightly decreased. This observed decrease in S was presumably a result of intrinsic carrier excitation across the bandgap at high temperatures, and it was consistent with the corresponding drastic rise in σ values as shown in Fig. 5(a). Furthermore, we have observed increased S with increasing Al content in most samples. This type of phenomenon, explained by an overall decrease in background hole n and enhanced carrier scattering upon the introduction of Al into the system, can be understood using the following equation: , where k_B = Boltzmann constant, h = Planck's constant, e = electron charge, and m* is the hole's effective mass.² S is directly proportional to m*, which in turn is inversely related to µ, and inversely proportional to n. Therefore, the observed decrease in the n and µ (Table 4) as a result of Al inclusion nicely justified the increase in the magnitude of S. Lastly, the quinary compounds displayed relatively smaller S values compared to those of the two quaternary compounds, and these results should be attributed to the relatively higher n and µ in the two Eu-containing compounds. The observed maximum S values were 205, 246, 153, 171, and 175 µV K⁻¹ for Ba_1.84(2)Zn_1.84(2)Al_0.16Sb₂, BaZn_1.75(2)Al_0.25Sb₂, BaZn_1.61(2)Al_0.39Sb₂, Ba_0.96(2)Eu_0.04Zn_1.93(2)Al_0.07Sb₂, and Ba_0.97(1)Eu_0.03Zn_1.83(1)Al_0.17Sb₂, respectively. In addition, the power factor (PF), which was evaluated using a simple equation of PF = S²σ, showed an increase with temperature in all cases, but the values were relatively smaller than those of the parental compound, as shown in SI Fig. S11(a). An observed maximum PF was 3.85 µW cm⁻¹K⁻² at 760 K for the quaternary BaZn_1.61(2)Al_0.39Sb₂. While the PF of the substituted samples is lower than that of the parent compound, this trade-off is a necessary consequence of the strategy to optimize n and, more importantly, to drastically reduce κ_tot.

Third, the temperature-dependent κ_tot of the four title compounds displayed significantly lower κ_tot values than the ternary BaZn₂Sb₂ as shown in Fig. 6(a). This result should be attributed to the introduction of Eu and/or Al substituents, which lead to enhanced cationic and/or anionic site disorder and in turn were responsible for increased overall phonon scattering. Our co-substitution strategy is, therefore, exceptionally effective at suppressing thermal transport.

Fig. 6 Temperature-dependent (a) total thermal conductivity κ_tot and (b) figure-of-merit ZT of the four title compounds in the Ba_1−xEu_xZn_2−yAl_ySb₂ (0 ≤ x ≤ 0.03(1), 0.16(2) ≤ y ≤ 0.39(2)) system measured between 300 and 773 K.

Furthermore, κ_tot decreased gradually as the Al content increased due to a substantial reduction in the lattice component κ_latt of the κ_tot as shown in SI Fig. S11(b). Interestingly, the κ_tot of the quinary Ba_0.97(1)Eu_0.03Zn_1.83(1)Al_0.17Sb₂ with both cationic and anionic site disordering was slightly larger than those of the two quaternaries BaZn_1.75(2)Al_0.25Sb₂ and BaZn_1.61(2)Al_0.39Sb₂ with higher levels of Al substitutions. The observed discrepancy should be attributed to a relatively higher hole n of the quinary compound, which increased the electrical component κ_elec of κ_tot. Overall, we have observed ultra-low κ_tot values in the title Ba_1−xEu_xZn_2−yAl_ySb₂ system, and the minimum κ_tot values were 0.60, 0.24, 0.26, and 0.31 W m⁻¹K⁻¹ for Ba_1.84(2)Zn_1.84(2)Al_0.16Sb₂, BaZn_1.75(2)Al_0.25Sb₂, BaZn_1.61(2)Al_0.39Sb₂, and Ba_0.97(1)Eu_0.03Zn_1.83(1)Al_0.17Sb₂, respectively. In general, the κ_tot value can be expressed as a summation of the electrical term (κ_elec) and the lattice term (κ_latt): κ_tot = κ_elec + κ_latt, and the κ_elec value can be directly calculated by the Wiedemann–Franz law: κ_elec = LσT (L = Lorenz number).⁴⁵ Based on these correlations, we calculated the κ_elec values of 0.04, 0.01, 0.03, and 0.06 W m⁻¹K⁻¹ at 298 K for Ba_1.84(2)Zn_1.84(2)Al_0.16Sb₂, BaZn_1.75(2)Al_0.25Sb₂, BaZn_1.61(2)Al_0.39Sb₂, and Ba_0.97(1)Eu_0.03Zn_1.83(1)Al_0.17Sb₂, respectively, and the resultant κ_latt values were plotted in SI Fig. S11(b). Therefore, among the title compounds, we can claim that the contribution from κ_elec appeared to be slight, and κ_latt was dominant for the overall κ_tot.

Finally, the temperature-dependent TE figure of merit ZT values are plotted in Fig. 6(b). Three Al-substituted quaternary compounds possessed relatively smaller ZT values than the ternary and quinary compounds at lower temperatures. However, these values continued to increase with temperature, and the Al-richest quaternary BaZn_1.61(2)Al_0.39Sb₂ displayed the maximum ZT of 0.49 at 755 K based on its moderately high σ and ultra-low κ_tot at high temperatures. In contrast, the Eu and Al co-substituted quinary Ba_0.97(1)Eu_0.03Zn_1.83(1)Al_0.17Sb₂ showed slightly smaller ZT values than the ternary BaZn₂Sb₂ in lower temperatures. However, within the 550–600 K temperature range, it eventually exceeded the ternary compound and reached a peak ZT value of 0.50 at 653 K. While the ZT enhancement is modest, the results demonstrate a successful proof-of-concept and reveal the synergistic effect of our co-substitution strategy.

3.4. Machine learning validation

The independent and comprehensive AI-assisted ML predictions were performed for the five title compounds as well as the reference ternary compound to cross-validate and elucidate the experimentally obtained TE properties and ZT values. We exploited our own customized dataset containing 5798 experimentally measured, temperature-dependent TE properties and compositions across 965 unique materials. Crucially, this dataset is enriched with 755 data points from 85 chemically similar AM₂M′₂ Zintl compounds, providing the model with deep domain-specific experience in the exact chemical space of this study. The elemental count within the dataset is visualized in a periodic table heat map presented in Fig. 7. In addition, SI Fig. S12–S14 provide additional visualizations to complement the explanation of the ML prediction process, including a pie chart illustrating the proportion of TE material families in the dataset, a figure depicting the correlation between TE property values and temperature intervals for TE families, and histograms showing the distribution of TE property values within the dataset, respectively. Based on our customized dataset and the 154-dimensional feature vectors generated by the MAGPIE library, three independent ML models for σ, S, and κ_tot predictions were initially trained. By transforming data to a logarithmic scale, models can better capture relationships, stabilize variance, and improve performance, especially when dealing with skewed target variables as shown in SI Fig. S13 and S14.⁴⁶

Fig. 7 Visualization of elemental frequency in the customized dataset, shown as a heat map on a periodic table.

The three independent ML models for targeting TE properties were initially created using the XGBRegressor algorithm, and the performance of each model for predicting the log₁₀ σ, S, and log₁₀ κ_tot was evaluated by the various scatter plots as shown in Fig. 8(a)–(f). These models exhibited exceptional predictive performance on the training dataset, as evidenced by the close alignment of data points with the red dashed line representing perfect prediction. In particular, the predictions for S and log₁₀ κ_tot demonstrated remarkable reliability to the training data, which indicated the capacity of models for accurately capturing underlying trends. As applied to the unseen test dataset, the models also displayed robust generalization capabilities although a somewhat greater degree of scatter was observed relative to the training set. Interestingly, the predictions for S and log₁₀ κ_tot remained consistent and reliable across the test dataset, while predictions for log₁₀ κ_tot exhibited the most considerable deviations indicating significant dispersion from the ideal prediction line. In addition, several pronounced outliers were observed in the predictions for log₁₀ κ_tot in both the training and test datasets, which also implied potential limitations of these models to predict this property accurately.

Fig. 8 Scattered plot for the predictions of TE properties log₁₀ σ ((a) and (b)), S ((c) and (d)), and log₁₀ κ_tot ((e) and (f)) using the training and testing dataset at different temperatures. The red dashed line represents perfect predictions.

The predictive performance of the ML models was rigorously evaluated using a 5-fold cross-validation methodology, and the results are provided in SI Fig. S15.³⁶ The model for the prediction of κ_tot demonstrated the strongest generalization capability, achieving an average coefficient of determination (R²) of 0.8672 (±0.0253) and a Mean Absolute Error (MAE) of 0.0863 (±0.0079) on the test sets. In addition, the model for S also showed good predictive power with an average test R² of 0.7985 (±0.0544) and an MAE of 49.9162 (±6.3718). On the other hand, the model for σ was found to be the most challenging to predict, and the result showed a test R² of 0.7300 (±0.0481) and an MAE of 0.2704 (±0.0223). As a result of this cross-validation, we can claim that all three models fit the training data exceptionally well (R² > 0.95), which strongly indicates that the underlying patterns were successfully learned. Over-fitting in ML refers to the phenomenon where a model learns the training data too well, even including noise as well as irrelevant details. As a result, the trained model shows a poor performance on new and unseen data.⁴⁷ In this study, a performance gap was observed between the training and test sets for all three models, which was indicative of some degree of overfitting, and it was most pronounced for the σ. We attributed this type of overfitting to the larger variance in the σ dataset, containing a mixture of metallic, semiconducting, and heavily doped semiconducting materials. However, despite some degree of overfitting, overall these models demonstrated strong predictive power on unseen data, particularly for κ_tot (R² = 0.8672), confirming their ability to generalize. To rigorously evaluate the performance and generalization capability of the trained XGBRegressor models, a residual error analysis was also performed and visualized in SI Fig. S16(a)–(f). The results proved that error distributions, while centered at zero, were significantly wider for the test set than the training set, which indicated a degree of model overfitting and reduced generalization performance. The more detailed residual error analysis and the overall “error analysis” for these three predictive models were presented in appropriately named sections in the SI.

Feature importance in an XGBRegressor provides an importance score for each feature, indicating its contribution to the predictions of the model.^47,48 A higher score suggests a more influential feature in predicting the target property of the model. On the other hand, the feature engineering test verifies the effectiveness of engineered features by evaluating their impact on the performance metrics of a model. The resultant 20 most important features are presented in SI Fig. S17(a)–(c), and the performance of a model as a function of reducing the number of features (i.e., feature engineering) for log₁₀ σ, S, and log₁₀ κ_tot is also displayed in SI Fig. S18(a)–(c), respectively. More details regarding the feature importance and feature engineering can be found in appropriately named sections in the SI.

The resultant ML predicted σ, S, and κ_tot are displayed in Fig. 9(a)–(c). In addition, the calculated ZT values based on those ML-predicted TE values are also shown in Fig. 9(d). First, in Fig. 9(a), the ML predicted σ of the quaternary and quinary compounds was lower than that of the ternary BaZn₂Sb₂ at all temperature regions. The σ magnitudes increased with temperature for the three quaternary title compounds, indicating a typical semiconducting behavior. On the other hand, it fluctuates around initial values with increasing temperature in the case of the two quinary title compounds. Furthermore, the σ of the Al-substituted three quaternary compounds was observed to be relatively smaller than those of the Eu and Al co-substituted quinary compounds. These observations were mainly consistent with our experimental results, and the rationale behind them, in terms of net hole concentration, carrier mobility, and electronegativity differences between elements, had already been explained in the previous section. Second, the predicted S showed all positive values, indicating holes as the dominant carriers, and these generally increased with increasing temperature, as shown in Fig. 9(b). The S values were observed to be the largest for the ternary compound, and these values decreased upon the introduction of cationic and anionic substituents. In particular, the quaternary compounds possessed a relatively higher S value than the quinary compounds, which was generally consistent with the experimental observations. Third, the predicted magnitudes of the κ_tot were relatively small as usually observed in many Zinlt phases, and these values decreased with increasing temperatures, as shown in Fig. 9(c). Moreover, the magnitudes of κ_tot were smaller in the quinary compounds than in the quaternary compounds, and this type of reduction should be attributed to the additional disorder upon the Eu introduction. However, the ML predicted κ_tot was relatively larger than our experimental observation, where the reduction of the experimental κ_tot came from compositional disorder at the cationic and/or anionic sites. The relaxed structure of undoped ternary BaZn₂Sb₂ was devoid of any local atomic disorder. Thus, the κ_tot of this reference compound was solely determined by its complex crystal structure, phonon-hole scattering, and phonon–phonon scattering.⁴⁹ Furthermore, the ternary BaZn₂Sb₂ tended to have a larger phonon mean free path compared to its Eu and Al co-substituted variants. On the other hand, multiple cationic and anionic substitutions with different atomic mass in the quaternary and quinary compounds induced mass disorder scattering, enhanced electron-phonon coupling, and enabled manipulation of phonon vibrational modes.⁵⁰ Therefore, this kind of complicated relationship between substituents and structural distortions caused additional scattering mechanisms to come into play, which eventually resulted in experimentally obtained ultra-low κ_tot values as observed in our experiments. However, our ML model used to predict κ_tot was not designed to incorporate these complex electronic and physical environments during the predictions, resulting in a relative overestimation for κ_tot. Fourth, the calculated TE figure of merit ZT, based upon the ML predicted σ, S, and κ_tot value, increased with increasing temperatures as shown in Fig. 9(d). In particular, the ZT values showed a clear distinction between the quaternary and quinary compounds, and the latter displayed slightly higher ZT values than the former. The maximum ZT value of 0.43 was predicted for the quinary Ba_0.97Eu_0.03Zn_1.83Al_0.17Sb₂ at 750 K, and this value was quite close to the experimental observations.

Fig. 9 ML predicted (a) σ, (b) S, (c) κ_tot, and (d) ZT for the five title compounds over the temperature range of 300–773 K.

A direct comparison between the experimental and the ML predicted ZT values was made for the four title compounds, and the resultant relative error percentages were also illustrated in Fig. 10. The experimentally measured ZT values at the highest measurement temperatures for the four title compounds of BaZn_1.84Al_0.16Sb₂, BaZn_1.75Al_0.25Sb₂, BaZn_1.61Al_0.39Sb₂, and Ba_0.97Eu_0.03Zn_1.83Al_0.17Sb₂ were 0.47, 0.43, 0.49, and 0.44, respectively, while our ML predicted corresponding ZT values were 0.41, 0.39, 0.38, and 0.43. Notably, this comparison revealed a strong predictive accuracy for our ML models with relative errors falling within an acceptable 3–23% range,^51–53 and these results signified their ability to capture the complex, non-linear relationships between composition and the resulting TE properties.

Fig. 10 Comparison of the experimental and the ML-predicted ZT values for the four title compounds at 760 K. The dashed-line plot (right y-axis) illustrates the corresponding relative error percentage between the experimental and the ML-predicted values.

The most significant result from the ML analysis, however, is the correct prediction that all Al-substituted compounds remain p-type materials (Fig. 9(b)). This finding contradicts the idealized DFT calculations but aligns perfectly with our experimental data. This is not a coincidence, but a powerful, data-driven validation of our central hypothesis. The ML model does not operate on simple theoretical rules; it learns the complex, non-linear patterns of real-world chemistry directly from the thousands of experimental data points it was trained on. This dataset, rich in Zintl phases, has implicitly taught the model the nuanced rules of their chemical physics: that a baseline p-type character from incomplete charge transfer is common, and that the electronegativity of substituents plays a critical role. Therefore, the model's prediction of p-type behavior is a direct confirmation that the final electronic properties in this system are governed by a net chemical balance of competing influences.

Ultimately, this work demonstrates a complete and successful research cycle. We proposed a substitution strategy and observed a p-type outcome that contradicted idealized models. We then formulated a new hypothesis based on the competing chemical influences of the constituent elements within the Zintl phase and used a data-driven ML model as an independent arbiter. The model, by learning from the collective experience of thousands of real experiments, not only predicted the final performance with reasonable accuracy but also validated the subtle chemical physics governing the system. Although the experimentally obtained maximum ZT value near 0.50 is modest, these results prove the viability of our combined experimental-computational pathway. The synergistic effect between experiment and ML demonstrated here is a powerful paradigm for understanding complex materials and accelerating the search for new advanced thermoelectrics.

4. Conclusion

This work investigated a synergistic co-substitution strategy in the BaZn₂Sb₂ Zintl phase, leading to the successful synthesis of five new solid solutions in the Ba_1−xEu_xZn_2−yAl_ySb₂ system, all crystallizing in the orthorhombic BaCu₂S₂-type structure. The overall crystal structure was described as an assembly of the complex 3D ³_∞[(Zn/Al)Sb3/4Sb1/4] anionic framework and the cationic elements filling the cage-shaped voids within the framework. The observed phase selectivity for the BaCu₂S₂-type phase over the CaAl₂Si₂-type phase was nicely elucidated by our r₊/r₋ radius ratio criterion. A series of DFT calculations proved that the title compounds found a way to maintain their structure type energetically favorable by properly adjusting the amount of substituents, and the ELF diagram successfully visualized the modified bond polarity due to the atomic substitutions. The measured temperature-dependent κ_tot values of all compounds were ultra-low, which nearly compensated for the relatively low PF values, and as a result, the maximum ZT of 0.50 at 653 K and 0.49 at 755 K were observed for Ba_0.97Eu_0.03Zn_1.83Al_0.17Sb₂ and BaZn_1.61Al_0.39Sb₂, respectively. The AI-assisted ML models trained on a customized dataset showcased strong predictive performance and generalization capabilities for overall trends in σ and S, while somewhat over-predicted the κ_tot values. The ML-predicted ZT values aligned closely with experimental results with relative errors of 3–23% and confirmed the capability of these models to capture complex, non-linear relationships.

However, the central finding of this study is the elucidation of the competing chemical influences that govern the electronic properties of this system. We discovered a critical discrepancy: while DFT calculations on an ideal, defect-free lattice correctly predicted an n-type shift from Al substitution, our experimental transport measurements confirmed that all samples remain p-type. We have demonstrated that this is, most likely, not due to a defect-based self-compensation mechanism, but rather the predictable outcome of a chemical balance. The final carrier concentration is a net result of three competing factors: (1) the intrinsic p-type character of the parent BaZn₂Sb₂ phase, (2) the hole compensation from the n-type Al-donor, and (3) an increase in holes caused by the higher electronegativity of the Eu-substituent. Our Hall effect data provides direct, quantitative proof of this interplay. Crucially, this chemical model was independently validated by our custom-trained ML model, which, by learning from a vast dataset of real-world Zintl phases, correctly predicted the persistent p-type behavior, a confirmation of the underlying chemical physics, not just defect formation or simple electron-counting rules.

Ultimately, this study demonstrates a complete and powerful research cycle: a targeted synthesis strategy led to an experimental observation that challenged idealized models, which was explained by a new hypothesis based on competing chemical influences, and finally validated by an independent, data-driven ML model. This synergy of ML-driven validation and interpretation with complex experimental results is a powerful paradigm for accelerating the discovery and optimization of advanced TE materials.

Author contributions

Yunjeong Lee: formal analysis, investigation, visualization; Aziz Ahmed: software, writing – original draft; Jihyun Lee: data curation, formal analysis, investigation; Kang Min Ok: data curation, investigation, resources; Tae-Soo You: conceptualization, formal analysis, funding acquisition, project administration, writing – original draft, writing – review & editing.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

CCDC 2470327–2470330 and 2471602 contain the supplementary crystallographic data for this paper. CCDC 2470327 for Ba_0.97(1)Eu_0.03Zn_1.83(1)Al_0.17Sb₂, CCDC 2470328 for Ba_0.96(2)Eu_0.04Zn_1.93(2)Al_0.07Sb₂, CCDC 2470329 for BaZn_1.84(2)Al_0.16Sb₂, CCDC 2470330 for BaZn_1.61(2)Al_0.39Sb₂, and CCDC 2471602 for BaZn_1.75(2)Al_0.25Sb₂.^54a–e

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information: additional experimental details for ML; additional ML-related results and discussion; Tables S1–S3; Fig. S1–S18; references. See DOI: https://doi.org/10.1039/d5ta08713a.

Acknowledgements

This work was financially supported by the Basic Science Research Program and the Nano-Material Technology Development Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Science and ICT (RS-2024-00337629 and RS-2022-NR068194).

References

1. F. Li, M. Ruan, B. Jabar, C. Liang, Y. Chen, D. Ao, Z. Zheng, P. Fan and W. Liu, Nano Energy, 2021, 88, 106273.
2. G. J. Snyder and E. S. Toberer, Nat. Mater., 2008, 7, 105–114.
3. Y.-X. Zhang, Q.-Y. Huang, X. Yan, C.-Y. Wang, T.-Y. Yang, Z.-Y. Wang, Y.-C. Shi, Q. Shan, J. Feng and Z.-H. Ge, Nat. Commun., 2024, 15, 2736.
4. L. Wang, W. Zhang, S. Y. Back, N. Kawamoto, D. H. Nguyen and T. Mori, Nat. Commun., 2024, 15, 6800.
5. Z. Zhang, Y. Yan, X. Li, X. Wang, J. Li, C. Chen, F. Cao, J. Sui, X. Lin, X. Liu, G. Xie and Q. Zhang, Adv. Energy Mater., 2020, 10, 2001229.
6. Z. Zhang, X. Wang, Y. Liu, C. Chen, H. Yao, L. Yin, X. Li, S. Li, F. Zhang, F. Bai, J. Sui, B. Yu, F. Cao, X. Liu, J. Mao, G. Xie and Q. Zhang, J. Materiomics, 2019, 5, 583–589.
7. X. Shi, S. Song, G. Gao and Z. Ren, Science, 2024, 384, 757–762.
8. S. Fang, J. Li, K. Zou, H. Shuai, L. Xu, W. Deng, G. Zou, H. Hou and X. Ji, Chem. Eng. J., 2022, 433, 133841.
9. M. Jiang, Y. Fu, Q. Zhang, Z. Hu, A. Huang, S. Wang, L. Wang and W. Jiang, Natl. Sci. Rev., 2023, 10, nwad095.
10. Z. Zhang, J. Li, H. Yao, Q. Wang, L. Yin, K. Liu, X. Ma, M. Yuan, R. Wang, S. Duan, X. Bao, J. Cheng, X. Wang, X. Li, J. Shuai, J. Sui, X. Lin, X. Tan, X. Liu, J. Mao, G. Xie and Q. Zhang, Acta Mater., 2024, 268, 119777.
11. K.-F. Liu and S.-Q. Xia, J. Solid State Chem., 2019, 270, 252–264.
12. D. Morgan and R. Jacobs, Annu. Rev. Mater. Res., 2020, 50, 71–103.
13. X. Wang, Y. Sheng, J. Ning, J. Xi, L. Xi, D. Qiu, J. Yang and X. Ke, J. Phys. Chem. Lett., 2023, 14, 1808–1822.
14. G. S. Na and H. Chang, npj Comput. Mater., 2022, 8, 214.
15. J. Jeong, D. Shim, P. Yox, M.-H. Choi, K. M. Ok, K. Kovnir, G. J. Miller and T.-S. You, Chem. Mater., 2023, 35, 3985–3997.
16. F. Gascoin, S. Ottensmann, D. Stark, S. M. Haïle and G. J. Snyder, Adv. Funct. Mater., 2005, 15, 1860–1864.
17. X.-J. Wang, M.-B. Tang, H.-H. Chen, X.-X. Yang, J.-T. Zhao, U. Burkhardt and Y. Grin, Appl. Phys. Lett., 2009, 94, 092106.
18. J. Jeong, D. Shim, M.-H. Choi, K. M. Ok and T.-S. You, Bull. Korean Chem. Soc., 2024, 45, 165–170.
19. J. Jeong, D. Shim, M.-H. Choi, Z. Yunxiu, D.-H. Kim, K. M. Ok and T.-S. You, J. Alloys Compd., 2024, 1002, 17527.
20. D. Shim, J. Lee, A. Ahmed, J. H. Pi, M.-H. Choi, K. M. Ok, K. H. Lee and T.-S. You, Molecules, 2025, 30, 310.
21. X.-J. Wang, M.-B. Tang, J.-T. Zhao, H.-H. Chen and X.-X. Yang, Appl. Phys. Lett., 2007, 90, 232107.
22. Rigaku Corporation, PDXL2 (Version 2.8.4.0), Rigaku, Tokyo, Japan, 2011.
23. Bruker Nano, Bruker APEX3 (Version 2019.1-0), Bruker AXS Inc., Madison, WI, USA, 2006.
24. A. Bruker, Bruker SAINT (Version 8.40A), Bruker AXS Inc., Madison, WI, USA, 2002.
25. S. M. Sheldrick, SADABS (Version 2016/2), University of Göttingen, Göttingen, Germany, 2003.
26. J. W. Shin, K. Eom and D. Moon, J. Synchrotron Radiat., 2016, 23, 369–373.
27. Z. Otwinowski and W. Minor, Methods Enzymol., 1997, 276, 307–326.
28. O. K. Andersen and O. Jepsen, Phys. Rev. Lett., 1984, 53, 2571–2574.
29. W. R. Lambrecht and O. K. Andersen, Phys. Rev. B:Condens. Matter Mater. Phys., 1986, 34, 2439–2449.
30. O. Jepsen, and A. Burkhardt, The Program TB-LMTO-ASA (Version 4.7), Max-Plank-Institut fur Festkorperforschung: Shuttgart, Germany, 1999.
31. O. K. Andersen, O. Jepsen and D. Glötzel, in Highlights of Condensed Matter Theory, ed. F. Bassani, F. Fumi and M. Tosi, Elsevier North Holland, New York, USA, 1985, 65–72.
32. O. Jepsen, O. Andersen and Z. Physik, Phys. B Condens. Matter, 1995, 97, 35–47.
33. P. E. Blöchl, O. Jepsen and O. K. Andersen, Phys. Rev. B:Condens. Matter Mater. Phys., 1994, 49, 16223–16233.
34. K. A. Borup, J. De Boor, H. Wang, F. Drymiotis, F. Gascoin, X. Shi, L. Chen, M. I. Fedorov, E. Muller and B. B. Iversen, Energy Environ. Sci., 2015, 8, 423–435.
35. L. J. van der Pauw, Philips Res. Rep., 1958, 13, 174–182.
36. N. Parse, N. J. Recatala-Gomez, R. Zhu, A. K. Y. Low, K. Hippalgaonkar, T. Mato, Y. Katsura and S. Pinitsoontorn, Adv. Theory Simul., 2024, 7, 2400308.
37. R. J. Murdock, S. K. Kauwe, A. Y.-T. Wang and T. D. Sparks, Integr. Mater. Manuf. Innov., 2020, 9, 221–227.
38. L. Ward, A. Agrawal, A. Choudhary and C. Wolverton, npj Comput. Mater., 2016, 2, 16028.
39. R. Arora, K. Kumar and S. Dixit, Asian J. Civ. Eng., 2024, 25, 1517–1530.
40. G. S. Pomrehn, A. Zevalkink, W. G. Zeier, A. van de Walle and G. J. Snyder, Angew. Chem., Int. Ed., 2014, 53, 3422–3426.
41. A. Savin, R. Nesper, S. Wengert and T. E. Fassler, Angew. Chem., Int. Ed. Engl., 1997, 36, 1808–1832.
42. J. Emsley, The Elements, Oxford University Press, New York 1989.
43. K. Guo, Q. Cao and J. Zhao, J. Rare Earths, 2013, 31, 1029–1038.
44. R. Yan, W. Lv, K. Wang, K. Guo, X. Yang, J. Luo and J. Zhao, J. Mater. Chem. A, 2016, 4, 12119–12125.
45. H.-S. Kim, Z. M. Gibbs, Y. Tang, H. Wang and G. J. Snyder, APL Mater., 2015, 3, 041506.
46. K. E. Vipin and P. Padhan, J. Mater. Chem. C, 2024, 12, 7415–7425.
47. X. Zhong, B. Gallagher, S. Liu, B. Kailkhura, A. Hiszpanski and T. Y.-J. Han, npj Comput. Mater., 2022, 8, 204.
48. F. Oviedo, J. L. Ferres, T. Buonassisi and K. T. Butler, Acc. Mater. Res., 2022, 3, 597–607.
49. E. S. Toberer, A. F. May and G. J. Snyder, Chem. Mater., 2010, 22, 624–634.
50. P. Sharma, P. Singh and G. Balasubramanian, Carbon, 2024, 223, 119015.
51. Z. Guo, P. R. Chowdhury, Z. Han, Y. Sun, D. Feng, G. Lin and X. Ruan, npj Comput. Mater., 2023, 9, 95.
52. N. K. Barua, A. Golabek, A. O. Oliynyk and H. Kleinke, J. Mater. Chem. C, 2023, 11, 11643–11652.
53. S. Kim, M. Alizamir, S. Heddam, S. W. Chang, I.-M. Chung, O. Kisi and C. Kulls, Sci. Rep., 2025, 15, 35523.
54. (a) CCDC 2470327: Experimental Crystal Structure Determination, 2026,  DOI:10.25505/fiz.icsd.cc2nxkz4; (b) CCDC 2470328: Experimental Crystal Structure Determination, 2026,  DOI:10.25505/fiz.icsd.cc2nxl06; (c) CCDC 2470329: Experimental Crystal Structure Determination, 2026,  DOI:10.25505/fiz.icsd.cc2nxl17; (d) CCDC 2470330: Experimental Crystal Structure Determination, 2026,  DOI:10.25505/fiz.icsd.cc2nxl28; (e) CCDC 2471602: Experimental Crystal Structure Determination, 2026,  DOI:10.25505/fiz.icsd.cc2nyx3n.

---