Modulating phonon dynamics: tailoring lattice vibrations to enhance thermoelectric efficiency in Mg_3(Sb, Bi)_2 alloy

Domain: physics
Published: 2025-11-24
Editorial outcome: minor_revision

Abstract
Heat in crystalline materials is transported by phonons from lattice vibrations, and lattice thermal conductivity critically determines thermoelectric performance. Different from conventional approach that reduce thermal conductivity via extrinsic additives sacrificing electrical transport, here, we demonstrate a notable advancement in the n-type Mg3Sb1.5Bi0.5 by modulating phonon dynamics through lattice softening and simultaneously suppressing the phonon mean free path in a more localized manner while remaining compositionally invariant. Originating from Mg vacancies and derivative defects, elevated internal strain degrades bonding rigidity and localize phonons at the lattice-constant level, yielding an ultra-low thermal conductivity of 0.3 Wm⁻¹K⁻¹, close to the theoretical minimum. This intrinsic strategy, combined with electron concentration optimization, yields a ZTmax of 2.06 and an extraordinary ZTave of 1.58, exceeding state-of-the-art n-type materials. Furthermore, a single-leg generator and two-pair module deliver conversion efficiencies of 12.5% (ΔT=440 K) and 7.4% (ΔT=300 K), respectively, highlighting exceptional potential for waste heat recovery.

Introduction
More than half of the energy from fossil fuels is lost as unusable heat during combustion, especially waste heat in the range of near room temperature to below 500°C, for which there is currently no effective method to capture and convert this low-grade energy into useful power1,2. Thermoelectric (TE) generators could directly convert heat and electricity via the Seebeck effect, thereby efficiently utilizing waste heat and promoting sustainable energy development3–5. The crucial factor in enhancing the conversion efficiency (η) of TE generators relies on achieving a high dimensionless figure of merit, ZT=S2σΤ/κtot, where S, σ, Τ and κtot stand for the Seebeck coefficient, electrical conductivity, temperature, and total thermal conductivity, containing electronic (κele), lattice (κlat) and bipolar (κbip) thermal conductivity, respectively. Realizing a high ZT hinges on a high power factor (S2σ) and low thermal conductivity, yet these parameters are coupled together and varied by carrier concentration and mobility6,7.

It is, in general, difficult to enhance any individual TE parameter without compromising the others8,9. Among S, σ, and κ, the κlat is regarded as the only parameter that can be tuned independently10,11. Thermodynamically, the κlat is basically dominated by the phonon propagation dynamics, quantized lattice vibrations that govern heat transport. Usually, traditional strategies for κlat regulation predominantly rely on introducing nanostructures, including secondary phases, interfaces, or grain boundaries on the tens to hundreds of nanometer scale to disrupt vibrational coherence, which reduces the phonon mean free path (l) and lowers the κlat12,13. While these structures effectively scatter mid-to-low-frequency phonons, their comparable dimensions to carrier mean free paths inevitably deteriorate carrier mobility. According to the phonon gas model14,15, the κlat could be expressed as16:

κlat = 1/3 C v l (1)

where C, v represent heat capacity and the sound velocity, respectively. Apart from tailoring lattice vibrations through phonon mean free path regulation, the modulation of sound velocity is often overlooked. To minimize κlat, two synergistic strategies appear. As illustrated in Fig. 1b12,17–22, it is crucial to simultaneously lessen the product of the mean free path and sound velocity (υ×l) in a localized manner, to effectively decouple thermal and electronic transport. When phonon-phonon interactions are the primary scattering mechanism, the phonon gas model formula could also be written as follows23,24:

κlat = (6π2)2/3 M v3 / V2 34π2γ2 T (2)

where V, γ, and M symbolize the atomic volume, Grüneisen parameter, and mean atomic weight, respectively. It can be seen that κlat is more sensitive to lattice bonding stiffness. In other words, reducing the elastic properties, comprising average sound velocity, Shear modulus, Young’s modulus, or increasing the Grüneisen parameter and lowering the Debye temperature, weakens the bonding stiffness, which causes lattice softening and loosens up phonon propagation, and minimizes the effect on electron transport due to the defect-dominated charge scattering. This dual strategy, centered on tailoring lattice vibrations locally, reduces the adverse correlation between κlat and carrier mobility, which is a cornerstone for optimizing TE performance.

Bi2Te3 is currently the most widely used TE material25–28, but its scarce tellurium and poor mechanical properties, especially for its n-type counterpart, have long plagued the industry29,30. Ever since Tamaki et al. achieved a high ZT value of 1.5 in n-type Mg3(Sb,Bi)2 by adding excess Mg content to lower the formation energy of Mg vacancies, Mg3(Sb,Bi)2 has attracted widespread attention for its high performance, low cost and plasticity31–42. Owing to the complex hierarchical structure, Mg3(Sb,Bi)2 naturally exhibits a low intrinsic lattice thermal conductivity (κlat ≈ 1 Wm−1K−1 at room temperature)43,44, though it is still sizably larger than the theoretical limit. Due to the high vapor pressure of Mg, many strategies have emerged to refine the electrical transport properties. For example, increasing the sintering temperature, synthesizing single crystals and annealing in the Mg atmosphere to boost grain size45–48, or doping with elements, such as Nd, Mn, and Cu to alter electron scattering mechanisms have all been employed42,49–53. These efforts primarily focus on electronic conduction. However, the critical dimension of phonon dynamics modulation through lattice vibrational tailoring remains insufficiently addressed.

In this work, in contrast to the conventional approach that depends on substantial fluctuations in chemical makeup by extrinsic additives to reduce κlat at the price of electrical properties, we propose an atomic-scale localized defect engineering strategy: constructing sub-nanometer point defects, such as vacancies with associated local strain fields within the lattice, prioritizing atomic-scale phonon scattering while strategically suppressing electronic degradation. Vacancy-induced strain fields are especially pronounced compared to those from extrinsic dopants. Taking n-type Mg3(Sb,Bi)2-based material as a prototype, its TE properties are known to be highly sensitive to Mg vacancies. However, conventional approaches through crude reduction of initial Mg content fail to precisely control vacancy concentration, with excessive vacancies triggering conductivity type inversion to p-type conduction. To address this, we exploit the characteristic of high saturated vapor pressure of Mg at elevated temperatures, dynamically tuning vacancy concentration via prolonged sintering time. After carrier concentration tuning by MnTe, intrinsic lattice softening amplification coupled with selective phonon mean free path restriction by leveraging conventional sintering techniques, synergistically suppresses phonon propagation. Originating from Mg vacancies and their derivative defects, elevated internal strain degrades bonding rigidity and defects suppress phonon mean free path, resulting in an ultra-low lattice thermal conductivity. The schematic of the main mechanism is shown in Fig. 1a. Although a certain number of dislocations and nanopores are likely to form due to vacancy clustering during prolonged sintering to affect the electrical performance, their impact may further mitigate by the carrier concentration and density of states effective mass decrease. As a result, through fine-tuning the Mg vacancy, the carrier concentration and mobility remain at high levels, maintaining the power factors in the x=0.01 (20 min) sample, and the product of phonon velocity and phonon mean free path decreases to an extremely low value below previous samples31,51,53–56, close to the theoretical minimum, achieving a decoupling of electrical and thermal transport. Consequently, one of the highest ZTmax of 2.06, notable average ZTave of 1.04 (300–500 K) and extraordinary ZTave = 1.58 (300–723 K) are obtained, showing its competitiveness to any other n-type TE materials near room temperature, and superiority in the mid-temperature region, as shown in Fig. 1d57–65. Most importantly, we achieve a high conversion efficiency of 12.5% (ΔT=440 K) and 7.4% (ΔT=300 K) for a single TE-leg generator and a two-pair module in Fig. 1e21,22,29,51,56,58,60,66–71, respectively. Compared to both commercial Bi2Te3 modules or other TE systems, the initial conversion efficiency obtained in this work exhibits strong competitiveness, highlighting great potential for waste heat recovery.

Results
Modulating thermal properties by tailoring lattice vibrations, including lattice softening and suppressed phonon mean free path

Firstly, the doped MnTe in Mg3.2Sb1.5Bi0.5 is to optimize the carrier concentration. As shown in Fig. S1, the electrical conductivity increases from 219 Scm−1 for x=0.005 to 481 Scm−1 for x=0.02 sample. The increase in electrical conductivity is primarily due to the rise in carrier concentration following MnTe doping. In addition, with the enhancement of MnTe content, the carrier-scattering mechanism shifts from being dominated by grain boundary scattering to mixed scattering before 400 K, which is beneficial for optimizing TE performance in the room-temperature range. Ultimately, a peak ZT of 1.73 at 723 K is achieved in the x=0.01 sample. So, the subsequent discussion will focus on the x=0.01 composition.

Tailoring lattice vibrations at the local scale is essential to decouple electrical-thermal performance by weakening the adverse correlation between κlat and carrier mobility. Unlike conventional extrinsic doping methods that rely on substantial fluctuations in chemical composition, we modulate phonon transport locally by creating sub‑nanometer point defects, such as vacancies with associated strain fields, that weaken bond stiffness and enhance phonon scattering, while strategically minimizing adverse effects on electron transport. However, excess Mg vacancies are believed to be detrimental to the performance of n-type Mg3(Sb, Bi)2-based materials. It is critical to reconcile this competing effect and tune the vacancy concentration to an optimal range in order to balance electrical and thermal transport. We adopt an initial composition with slight excess Mg and MnTe compound to ensure robust n-type conductivity with an adequate carrier concentration. Exploiting the characteristic of high saturated vapor pressure of Mg at elevated temperatures, Mg vacancies are dynamically tuning by extending the sintering time to achieve balanced TE transport, which is further discussed later. To further demonstrate the necessity of such fine Mg vacancy regulation, we also prepare samples with different initial Mg contents to investigate the effect of starting composition as shown in Fig. S2. These results show that the conventional approaches through crude reduction of initial Mg content fail to precisely control vacancy concentration. Mg3(Sb,Bi)2-based alloys belong to the α-La2O3-type structure. As Fig. S3 shows, the X-ray diffraction (XRD) patterns of (Mg3.2Sb1.5Bi0.5)1-x(MnTe)x samples at room temperature are well indexed to the Mg3Sb2 structure (JCPDS#03-0375). There is a single-phase structure without any impurity phases within the detection limits of XRD. We have also performed SEM analysis on two representative samples, x=0.01 (10 min) and x=0.01 (20 min) shown in Figs. S4 and S5. In both samples, the elemental mapping shows that Mg, Sb, Bi, Mn and Te are homogeneously distributed throughout samples. No obvious segregated MnTe compound is observed. To elucidate phonon transport behavior, the thermal performance parameters are presented in Fig. 2. As shown in Fig. 2a, the κtot near the room temperature gradually decreases as the sintering holding time is extended. A relatively low κtot of 0.6 Wm−1K−1 at 623 K is acquired for the x=0.01 (20 min) sample. The variations in κtot are associated with fluctuations in κele, shown in Fig. S6b, where κele is calculated by the Wiedemann–Franz law, κele = LσT, and L is the Lorenz number estimated by the empirical formula72 L=1.5+exp(–|S|/116)×10−8 V2K−2 and plotted in Fig. S6a. The κtot-κele, displayed in Fig. 2b, gradually decreases with rising temperature, while it experiences a slight increase at high temperatures due to the bipolar effect. It is interesting that a lower value of κlat is acquired after extending the sintering holding time, dropping from 0.88 Wm–1K–1 for the x=0.01 (10 min) sample to 0.57 Wm–1K–1 for the x=0.01 (20 min) sample at room temperature, a reduction of over 30%. Notably, the lowest value decreases to 0.3 Wm–1K–1 at 623 K, which is close to the theoretical minimum predicted by the diffusion limit73. The calculation of the theoretical minimum is provided in the Supporting Information. The reason for the diminishment of lattice thermal conductivity is primarily attributed to the phonon mode softening and enhanced phonon scattering, which is further discussed later.

Elastic properties are commonly used to gauge the inter-atomic bonding strength in a crystal lattice. To further scrutinize the reasons behind the reduction in thermal conductivity, we undertake investigation of the elastic properties to estimate the chemical bond strength. The calculation method for the elastic properties is provided in the Supporting Information. As we can see from Fig. 2e, at room temperature, the longitudinal (νl), transverse (νt) and mean (νs) sound velocities decrease obviously when the sintering holding time is extended. Meanwhile, elastic parameters including the Shear and Young’s modulus, have also noticeably decreased, as shown in Table S1. The lessened elastic properties suggest a weakening of the chemical bonding stiffness in the Mg3.2Sb1.5Bi0.5 material, which commonly engenders lattice softening and impedes phonon propagation. Besides elastic properties, the Debye temperature can also serve as an indicator of atomic bonding strength. As shown in Fig. 2f, the Debye temperature steadily decreases with extended sintering holding time, which aligns with the conclusions derived from the elastic properties. According to the ball-and-spring model by Heremans74, weaker inter-atomic bonding strength is associated with stronger anharmonicity, and as shown in Table S1, the sample’s anharmonicity also increases. To gain deeper insights into the reduction in elastic properties, we calculate the internal strain of powder samples via XRD. Figure 2h presents the measurement of internal strain according to the Williamson–Hall function75, with the calculation details provided in the Supporting Information, and ε represents the magnitude of the internal strain. It is obvious that there exists a significant increase in internal strain with prolonged sintering holding time, which primarily originates from the stress field induced by the higher density of VMg and their derivative defects, with further discussion later. Typically, defect-induced internal strain fields usually modify phonon frequencies, which in turn lead to lattice softening. Therefore, the lattice softening induced by internal strain leads to a reduction in thermal conductivity. The peak positions of the Raman spectrum for samples at room temperature gradually shift to the lower frequency area, shown in Fig. 2g, indicating the lattice softening as well. By the way, the phonon mean free path of samples decreases significantly with the extending sintering holding time displayed in Fig. 2f, especially for the x=0.01 (20 min) sample, which is reduced to 0.57 nm and close to the lattice parameters of the unit cell. This indicates that phonon transport has been nearly suppressed to its limit through scattering, contributing to the obtained extremely low lattice thermal conductivity.

To further quantitatively elucidate the effect in the reduction of κlat, we employ the Debye–Callaway model to theoretically analyse76:

κlat = kB / 2π2v (kB T / ℏ)3 ∫0θD/T τ(x) x4ex / (ex-1)2 dx (3)

where x=ℏω/kBT, kB, ℏ, τ are the reduced phonon frequency, Boltzmann constant, reduced Plank constant, and phonon relaxation time, respectively. More details and corresponding parameters can be seen in Supporting Information and Table S2. This model primarily considers phonon-phonon Umklapp process (U), grain boundary scattering (B), point defect scattering, dislocations scattering, strain scattering (DS), nanoprecipitates (NP) and nanopores (Pore). As shown in Fig. 2c, the experimental data align well with the theoretical calculations, confirming that lattice softening and defect scattering are the primary reasons for the reduction in lattice thermal conductivity. It can be also observed that the bipolar effect becomes more pronounced with the extending sintering holding time, because of the decrease in carrier concentration. The frequency-dependent κs for the x=0.01 (20 min) sample at room temperature is plotted in Fig. 2d. The reduction in lattice thermal conductivity originates from Mg vacancies and their derived defects, which act at different levels. The isolated Mg vacancies contribute at the atomic scale through point-defect scattering. And during prolonged sintering, Mg vacancies tend to aggregate, giving rise to vacancy clusters, dislocations, nanopores77, that act at larger length scales and further enhance phonon scattering. In this sense, the internal strain generated by Mg vacancies and their derivative defects leads to lattice softening, while these defects suppress the phonon mean free path. Together, they jointly constitute the main origin of the lattice thermal conductivity reduction. Thus, the obtained κlat is lower than those of most reported values, both at room and the peak ZT temperature, as shown in Fig. 2i51,53–56,78,79.

We further explore more microstructure analysis to gain insights into phonon transport behavior. A certain of Mg vacancies’ derived defects, such as dislocations and nanopores, form due to vacancy clustering during prolonged sintering. Some defects, including obvious dislocations and NP, have found through the low-magnification transmission electron microscopy images plotted in Fig. 3a–c. A certain of nanoscale pores forms from the STEM images in Fig. S7, which have the ability to scatter low- and medium frequency phonons. From the EDS elemental mapping in Fig. 3d, there is a significant Bi-rich aggregation near the grain boundaries, which is also confirmed through line scanning in Fig. S8. The same phenomenon is also observed in other regions of the same sample shown in Fig. S9. These NP have a diameter of nearly tens of nanometers, which often results in the formation of stacking faults and induces internal stress. The high density of lattice distortions, stacking fault and NP can also be noticed in the high-resolution transmission electron microscopy images within the grain and near the grain boundary shown in Fig. 3e, f, which is beneficial to impede phonon propagation. Figure 4g exhibits the strain maps by the geometric phase analysis (GPA) along different directions, inducing large strain fluctuation near the grains, which regularly triggers phonon mode softening. Furthermore, compared with the x=0.01 (10 min) sample, the dislocation density increases by nearly 50% in x=0.01 (20 min) sample from the EBSD analysis seen in Fig. 3h, i. All in all, widespread stress fluctuations and numerous defects are responsible for the lattice softening and phonon scattering, as depicted in the primary phonon scattering mechanisms diagram of Fig. 1a, which results in the obtained ultra-low κlat, and is consistent with the analysis above. Although the newly introduced dislocation defects also adversely affect electron transport, the reduced carrier concentration and density of states effective mass can mitigate this impact, as discussed later.

In addition, from the SEM images shown in Fig. S10, the number of micron-scale pores increases significantly from x=0.01 (20 min) to x=0.01 (30 min) samples as the sintering time further increases. The gradual decrease in sample density listed in the Table S3 with increasing sintering time also indirectly indicates a corresponding rise in pores. In the geometric-scattering regime, large pores scatter primarily long-wavelength phonons, which contribute partly to reduction of thermal resistance. Therefore, upon further extension to 30 min, the nanoscale pores tend to coarsen and agglomerate into micron-scale pores, resulting in a slight enhancement in the phonon mean free path and κlat compared to the x=0.01 (20 min) sample shown in Fig. 2b, f.

High electrical transport maintenance and ZT value enhancement

Besides thermal transport, maintaining high electrical performance is also crucial. As shown in Fig. 4a, with extending the sintering holding time, the carrier-scattering mechanism shifts toward grain boundary scattering. And the electrical conductivity deteriorates severely when the sintering holding time extends to 30 min. From the grain size distribution in EBSD images shown in Fig. S11, there is no significant change in the average grain size. So, the carrier-scattering mechanism alteration is attributed to the formation of Mg vacancies during prolonged holding time, which is also supported by the gradual decrease in carrier concentration, as shown in Fig. 4b. From the XRD refinement results shown in Figs. S12–S14 and Table S4, the Mg site occupancy decreased markedly, further confirming an increase in Mg vacancy concentration. Therefore, sustaining high electrical conductivity necessitates an appropriate sintering holding time, further evidencing the existence of a critical Mg vacancy concentration and the imperative for dynamically regulating vacancy concentration to decouple phonon-electron interactions in Mg3(Sb, Bi)2-based materials. Additionally, with the increase of Mg vacancy concentration, the carrier mobility does not deteriorate due to the reduction in carrier concentration, which mitigate the adverse impact of dislocation defects on carrier scattering.

The temperature-dependent Seebeck coefficient of samples is shown in Fig. 4c. The Seebeck coefficient gradually increases with rising temperature, and the negative values denote a typical n-type conducting behavior. The Seebeck coefficient shows no significant change in x=0.01 (10 min) and x=0.01 (20 min) samples, whereas x=0.01 (30 min) sample exhibits a marked increase, primarily due to the reduction in carrier concentration. The curve between carrier concentration and Seebeck coefficient is plotted in Fig. 4d to further consider the electronic transport properties on the basis of single parabolic band model. In samples with extended sintering holding time, the density of states effective mass decreases partially. We have performed EPMA compositional analysis on samples with different holding times. As we can see from Figs. S15–S17 and Table S5, all elements are homogeneously distributed, and no MnTe compound is observed. By the way, a slight decrease in Mg content is detected as the holding time enhancement, indicating an increased concentration of Mg vacancies. This phenomenon is consistent with the decrease of electron concentration. The observed decrease of the density of states effective mass can be related to the formation of Mg vacancies. A lower electron concentration shifts the Fermi level downward, away from the regions of higher density of states, thereby reducing the effective mass. Some similar phenomena have also been reported in previous studies54. The reduction in density‑of‑states effective mass can also mitigate the effect of carrier mobility by defect‑induced electron scattering.

As a result, one of the highest ZT value of 2.06 has been achieved at 623 K, as shown in Fig. 4e. Interestingly, the ZT curve quickly reaches its peak and then gradually levels off, which is very beneficial behavior for obtaining a high average ZT value. This temperature dependence is mainly attributed to the fact that as the temperature increases past the optimum performance, i.e., peak ZT, the enhanced carrier mobility scattering slightly weakens the improvement in electrical transport, while the lattice thermal conductivity in this work, already achieving values close to its theoretical limit, exhibits a relatively small decrease with temperature. Excitingly, the ZTave reaches 1.04 and 1.58 within the temperature range of 300–500 and 300–723 K shown in Fig. 4f, respectively. In comparison with previous reports on Mg3(Sb,Bi)2, or any other known n-type TEs, the ZTave in this work is at the state-of-the-art level and higher than that of commercial Bi2Te3. This indicates that the n-type Mg3(Sb, Bi)2 obtained in this work exhibits significant advantages for waste heat recovery across both the near-room-temperature and mid-temperature ranges. In addition, we reproduce multiple samples of x=0.01 (20 min) and repeatedly test them, and the consistency of the properties is excellent, as plotted in Figs. S18 and S19.

Design, fabrication and assessment of TE module

Both a single-leg and a two-pair TE device have been designed to better demonstrate the TE performance across different temperature ranges80–82. Figure 5a illustrates the structure of the single-leg device. To better harness the excellent TE performance and fully exploit the temperature difference, the finite element simulation is carried out. As shown in Fig. 5b, only the appropriate height and cross-sectional area can yield the maximum TE conversion efficiency. Thus, the dimensions of 3.4×3.4×6.5 mm3 have been adopted to fabricate the single-leg device. The contact resistance between the Mg3(Sb,Bi)2-based leg and the barrier layer is examined using the four-probe method, shown in Fig. S20. The low contact resistance around 10 μΩcm2, plotted in Fig. 5c, indicates a good contact. The slight difference in contact resistances at both ends of the TE leg can be attributed to minor variations in local thermal-mechanical conditions during sintering, but it remains within the typical experimental uncertainty. In addition, the microscopic morphology and element distribution of the interface have been shown in Fig. S21. The elemental distribution is uniform, and there is no apparent elemental diffusion or voids at the interface, consistent with the result of low contact resistance.

As shown in Fig. 5d for the single-leg device, the open-circuit voltage slowly improves with the rising temperature difference. The internal resistance, obtained by the slope of the V–I curve shown in Fig. S22b, increases gradually, primarily due to the decrease in electrical conductivity as the temperature rises. In addition, as the current increases, the output heat flow displayed in Fig. S22a gradually rises due to the enhanced Joule heating and Peltier effect. When the external resistance matches the internal resistance of the device, the output power and conversion efficiency reach their maximum, as shown in Fig. 5e. Ultimately, the maximum conversion efficiency of the single-leg device realizes a conversion efficiency of 12.5% at the 440 K temperature difference, competitive with the highest performance previous reports, as plotted in Fig. 5f60,78,83–88. Furthermore, coupled with MgAgSb, the prepared two-pair module after geometrical design, shown in Fig. S23, achieves a high conversion efficiency of 7.4% (ΔT=300 K). It is indicated that this work has great potential for waste heat recovery, whether in the near-room-temperature or mid-temperature range.

Methods
Sample synthesis
High-purity raw Mg (4N), Bi (4N), Sb (4N) and MnTe (3N) were weighed according to the ratio of (Mg3.2Sb1.5Bi0.5)1−x(MnTe)x (x=0, 0.005, 0.01, 0.015, 0.02). The elements were loaded into the ball milling jar inside an Ar-filled glove box with an O2 concentration below 0.1 ppm and then ball-milled (SPEX-SamplePrep 8000 Mixer Mill) for 5 h. The obtained powders were loaded into a φ10 mm graphite die inside the glove box. It was then immediately sintered using SPS (SPS-1080 System, SPS SYNTEX INC) at 700 °C under 60 MPa for 10, 20, and 30 min.

Measurement and characterization
The electrical conductivity and Seebeck coefficient of all samples were simultaneously measured using ZEM-3 under the He atmosphere. The thermal conductivity was estimated through κtot = DρCp, where thermal diffusivity (D) was measured by Netzsch LFA467, the density of samples (ρ) was calculated by the Archimedes method, and the specific heat capacity (Cp) was determined by the formula Cp = 3NR (1+1.3×10−4T−4×103T−2)/M89. The room temperature electron carrier concentration (n) was measured by the physical properties measurement system, Quantum Design based on the formula n=1/RH e and μ=σRH. The samples’ phase structure was identified by the XRD (SmartLab 3, Rigaku) with Cu Kα radiation. The microstructure of samples was characterized using electron backscatter diffraction (EBSD, JSM-7001F, JEOL Inc.), scanning electron microscope (SEM, Hitachi SU4800) fitted with an energy-dispersive spectroscopy (EDS, Horiba EMAX Evolution X-Max) and electron probe microanalysis (EPMA, JXA-8500F). The Raman spectroscopy was performed with the 532 nm excitation laser (inVia, Renishaw). The finer microstructure was characterized by scanning transmission electron microscopy (STEM, JEM-ARM-200F-B). The room temperature sound velocity was measured by the Model UVM-2 (Ultrasonic Engineering Co., Ltd).

TE module fabrication and stimulation method
To fabricate the single-leg generator, the sandwiched structure of (Mg3.2Sb1.5Bi0.5)0.99(MnTe)0.01 powder and Fe alloy (304 stainless steel) powder was employed, using the one-step sintering process. The obtained ingot was cut into dice with dimensions of 3.4×3.4×6.5 mm3. The two-pair TE module was fabricated on the basis of n-type (Mg3.2Sb1.5Bi0.5)0.99(MnTe)0.01 leg with the height of 5 mm and cross-section of 12.25 mm2, and p-type MgAgSb leg with the height of 5 mm and cross-section of 10.89 mm2. The MgAgSb leg preparation process can be seen from previous work22. The TE module conversion efficiency (η) was measured by Mini-PEM, (ADVANCE RIKO, Japan)90. The details of the TE module test and numerical modeling are provided in the Supporting Information.

The numerical simulation of power generation was performed with COMSOL Multiphysics software, version 6.1. The Heat Transfer Module and the Electric Currents Module were coupled through the TE Effect multiphysics interface to simulate the coupled temperature and electrical potential fields. A three-dimensional model was built in the software interface to represent TE legs with a similar geometry and dimension. The temperature-dependent electrical conductivity, Seebeck coefficient, and thermal conductivity were taken from actual measurement results. Herein, the electrical and thermal contact resistances between interfaces were not considered in our simulation model. The material domains were discretized using a physics-controlled mesh, with the element size set to the default (normal) setting, and a fixed temperature difference was applied between the hot and cold ends.

Statistical analysis
The electronic and thermal transport measurements are independent of sample dimensions, and the uncertainties mainly reflect the instrumental deviations. In particular, the systematic errors are approximately 3% for the Seebeck coefficient and 5% for the electrical conductivity. For the total thermal conductivity, the uncertainty is about 7% calculated from dκtot/κtot = √ (dρ/ρ)2 + (dCp/Cp)2 + (dD/D)2 (1% for the density ρ, 5% for the specific heat Cp and 5% for the thermal diffusion D). Consequently, the propagated error for the calculated ZT values is estimated to be approximately 10%. The computational formula involved is d(ZT)/ZT = √ 2×(dS/S)2 + (dσ/σ)2 + (dκ/κ)2.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this study, the authors aim to precisely control the concentration of Mg vacancies in Mg₃(Sb,Bi)₂ alloys to enhance the\nthermoelectric properties of the material. They begin with a slightly Mg-rich Mg₃.₂Sb₁.₅Bi₀.₅ alloy as the initial material, and\nsubsequently introduce MnTe to modulate the carrier concentration. By extending the sintering time, they dynamically adjust\nthe Mg vacancy concentration, ultimately achieving an excellent ZT value (ZTmax = 2.06). The average ZT value across the\ntemperature range of 300-723 K is also significantly high at 1.58. The authors believe that factors such as internal strain and\ndefect proliferation lead to the ultra-low lattice thermal conductivity (0.3 W m⁻¹K⁻¹), which reduces interatomic bond stiffness\nand confines phonon propagation to the scale of the lattice parameter.\nThis paper is likely to be of interest to researchers in the field of thermoelectric materials and devices. However, substantial\nimprovements are necessary before publication. Below are some questions that need to be clarified.\n1.) All equations in the main text and Supporting Information should be numbered, and the source of each equation should\nbe cited. In Line 165 of the Supporting Information, the expression for (tao_DS)^(-1) contains “gamma” and “r”. What do these\nsymbols represent? Is “gamma” the Grüneisen parameter?\n2.) On page 3, line 56 of the manuscript, the authors list an equation, but its origin is unclear. Ref. 22 does not appear to\ncontain this formula.\n3.) How was the Lorenz number, used for calculating the electronic thermal conductivity and plotted in Fig.S4a, derived?\n4.) In Fig. 3b, are the defects indicated by the arrows nanopores?\n5.) The authors have employed the Debye-Callaway model for theoretical analysis of their results. This model primarily\nconsiders phonon-phonon Umklapp processes (U), grain boundary scattering (B), point defect scattering (PD), dislocation\nscattering (DC), and strain scattering (DS).\nIn this study, what specific point defects are considered for point defect scattering? Does it refer to mass fluctuations and\nstrain field fluctuations caused by Bi atoms substituting Sb atoms? If so, how are Mg vacancies incorporated into this model?\nThe authors should clearly demonstrate how to calculate the mass fluctuations scattering parameter (Gamma_M) and strain\nfield fluctuations scattering parameter (Gamma_S).\n6.) What physical mechanisms cause strain scattering (DS) in the Debye-Callaway model?\n7.) What is the unit for N_D (Dislocation density) in Table S2? Is it cm⁻³ or cm⁻²?\n8.) In the Abstract, the authors state, “Bonding rigidity degradation and phonon mean free path suppression are dominantly\ngoverned by elevated internal strain, achieving lattice-constant-level phonon localization.”\nThis statement seems to imply that an increase in Mg vacancy concentration leads to increased internal strain, reduced\nbonding rigidity, and suppression of the phonon mean free path, suggesting that the decrease in lattice thermal conductivity\nis primarily due to the increase in Mg vacancy concentration.\nOn the other hand, the authors state, “Except for the vacancy-induced strain field decreasing κ_lat, a certain amount of\ndislocations and nanopores form due to vacancy clustering during prolonged sintering. Numerous defects, including obvious\ndislocations and nanoprecipitates, have been found through the low-magnification transmission electron microscopy images\nplotted in Figure 3(a)-(c).” (Line 221, page 10)\nThis suggests that the decrease in lattice thermal conductivity is not only due to the increase in Mg vacancy concentration,\npoint defect scattering (PD), dislocation scattering (DC), and strain scattering (DS); nanopores and nanoprecipitates, which\nincrease with sintering time, also contribute. Can the Debye-Callaway model account for the contributions of nanopores and\nnanoprecipitates? Does the theoretical analysis based on the Debye-Callaway model indicate that the increase in Mg\nvacancy concentration is the primary cause of the decrease in lattice thermal conductivity?\n9.) On Line 182, page 15 of the Supporting Information, \"k_B\" should be the Boltzmann constant.\n10.) In Fig. 5c, how was the contact resistance measured? What does the \"distance\" on the X-axis represent?\n11.) The authors need to provide a detailed explanation of how the COMSOL calculations were performed, including the\nversion used, the specific modules employed, and any other relevant details.", "reviewer_2": "The author precisely controls the Mg vacancy concentration in the n-type Mg3Sb1.5Bi0.5 system by regulating the holding\ntime of SPS. modulating phonon dynamics through lattice softening and simultaneously suppressing the phonon mean free\npath in a more localized manner while remaining compositionally invariant. This strategy combined with electron\nconcentration optimization, yields one of the highest ZTmax of 2.06, accompanied by an extraordinary average ZTave of\n1.58 (300-723 K), which exceeds other state-of-the-art n-type thermoelectric materials. In general, this is valuable research\nwork, but the following issues need to be addressed before formal publication.\n1, Existing studies have shown that Te element doping is an effective means to convert Mg-based systems into n-type.\nIn this paper, the authors mixed MnTe with raw materials via ball milling, followed by sintering and forming through SPS. Then,\nafter the above experimental treatments, in what form does MnTe exist in the matrix? Does it form a second phase or is it\ncompletely doped into the lattice? If it is a second phase, will the doping amount of MnTe increase with the extension of\nholding time, thereby exerting a significant impact on the performance? In addition, since the added content of MnTe is\nextremely low, which may be below the detection limit of XRD, it is suggested to supplement other characterization methods\nsuch as scanning electron microscopy for further analysis.\n2, In Figure S3, the sample with x = 0.01 shows an obvious diffraction peak of the second phase around 25°, while this\ndiffraction peak disappears when the holding time reaches 20 minutes. How do the authors explain this phenomenon?\n3, Prolonged holding time leads to significant changes in effective mass, which the authors attribute to possible internal\nstrain generated during the sintering process. However, in terms of rationality, it may be more appropriate to explain the\nchanges in effective mass by the generation of Mg vacancies or the dissolution of MnTe into the matrix. It is suggested that\nthe authors add EPMA compositional characterization to analyze the actual composition of samples with different holding\ntimes.", "reviewer_3": "In this work, the authors tailored lattice vibrations via dynamically regulating Mg vacancy concentration, pronounced internal\nstrain, and defect proliferation in n-type Mg3Sb1.5Bi0.5-based alloy, which softens interatomic bond stiffness and confines\nphonon propagation to distances comparable to the lattice parameter. As a result, one of the highest ZTmax of 2.06, ZTave\nof 1.04 (300-500 K) and unprecedented 1.58 (300-723 K) are obtained. Also, the superior TE materials alongside low\ncontact resistivity enable an excellent conversion efficiency of 12.5% (ΔT = 440 K) and 7.4% (ΔT = 300 K) for a single TE-leg\ngenerator and a two-pair module. This work is suggested for publication after addressing some comments/questions below.\n1. Figure 2, the figure caption should be in alphabetical order. For example, “(e) Sound velocities and (g) Raman spectrum\nfor the (Mg3.2Sb1.5Bi0.5)1-x(MnTe)x samples at room temperature. (f) The comparison of phonon mean free path, Debye\ntemperature”\n2. TEM characterization. It is recommended to reformat the scale bar and label the type of defect/dislocation\ncorrespondingly. For example, in Figure 3g, the scale bar of GPA should be -10% to 10%. Figure S6, the HRTEM\ncharacterization and line profile of atomic intensity would be preferable for the examination of vacancy and strain variation.\n3. The measurement uncertainty or error analysis should be provided to ensure the robustness and reproducibility of TE\nproperties (Figures in the main text and supplementary info) for the reader.\n4. Figure 4c, it is suggested to format the title of figures (e.g., position, font, size) in the revised submission.\n5. It is impressive to benchmark the ZT in Figure 4f, while the relevant refs are required to be cited appropriately.\n6. Device section. The geometry optimization of TE legs and contact characterization are crucial for efficient power\ngeneration (Figures 5 and S17). Some refs are suggested to be included in the device section (doi.org/10.1038/s41565-023-\n01457-5; doi.org/10.1126/sciadv.adf5701; http://doi.org/10.1016/j.joule.2024.08.009). Also, detailed modelling and contact\ncharacterizations are required to be supplemented in the methods.\n7. Figure 5c, the low contact resistivity achieved Mg3(Sb, Bi)2-based leg is impressive. Could the authors elaborate more on\nthe difference between contact resistivities obtained in the hot- (7.1 microohm square cm) and cold-side (10.2 microohm square\ncm)? Besides, in light of the minimization of contact resistivity and element oxidation, what are the selection criteria of\nbarrier/coating materials for Mg3(Sb, Bi)2-based legs?"}, "author_response": "We thank the Reviewers for their constructive comments which have helped us improve\nthe paper.\nResponse to Reviewer #1:\nIn this study, the authors aim to precisely control the concentration of Mg vacancies in\nMg (Sb,Bi) alloys to enhance the thermoelectric properties of the material. They begin\n3 2\nwith a slightly Mg-rich Mg Sb Bi alloy as the initial material, and subsequently\n3.2 1.5 0.5\nintroduce MnTe to modulate the carrier concentration. By extending the sintering time,\nthey dynamically adjust the Mg vacancy concentration, ultimately achieving an\nexcellent ZT value (ZT = 2.06). The average ZT value across the temperature range\nmax\nof 300-723 K is also significantly high at 1.58. The authors believe that factors such as\ninternal strain and defect proliferation lead to the ultra-low lattice thermal conductivity\n(0.3 W m-1K-1), which reduces interatomic bond stiffness and confines phonon\npropagation to the scale of the lattice parameter.\nThis paper is likely to be of interest to researchers in the field of thermoelectric\nmaterials and devices. However, substantial improvements are necessary before\npublication. Below are some questions that need to be clarified.\n1. All equations in the main text and Supporting Information should be numbered, and\nthe source of each equation should be cited. In Line 165 of the Supporting Information,\nthe expression for (tao_DS)^(-1) contains “gamma” and “r”. What do these symbols\nrepresent? Is “gamma” the Grüneisen parameter?\nResponse: We thank the reviewer for her/his kind remind. We have numbered all\nequations in both the main text and Supporting Information and cited their sources\naccordingly. In the expression of , the symbols “gamma” and “r” represent the\n−1\nGrüneisen parameter and the Poiss𝜏𝜏o𝐷𝐷n𝐷𝐷 ratio, respectively. For consistency and clarity,\nwe have replaced “r” with “υ ”. In addition, we have added corresponding notes in the\np\nSupporting Information to clarify these symbols (see page 18, middle).\n2. On page 3, line 56 of the manuscript, the authors list an equation, but its origin is\nunclear. Ref. 22 does not appear to contain this formula.\nreference for the equation (2) listed on page 3 is incorrect in the manuscript. We have\ncorrected this in the revised version and cite the appropriate references: (see page 3, up)\n1. Toberer ES, Zevalkink A, Snyder GJ. Phonon engineering through crystal chemistry.\nJ. Mater. Chem. 21, 15843 (2011).\n2. Hanus R, et al. Lattice Softening Significantly Reduces Thermal Conductivity and\nLeads to High Thermoelectric Efficiency. Adv. Mater. 31, e1900108 (2019).\n3. How was the Lorenz number, used for calculating the electronic thermal conductivity\nand plotted in Fig.S4a, derived?\nResponse: We thank the reviewer for her/his question. The Lorenz number used for\ncalculating the electronic thermal conductivity in Fig. S4a is obtained based on the\nempirical relation L = 1.5 + exp(–|S|/116) × 10–8 V2 K–2 between the Seebeck coefficient\nand the Lorenz number (APL MATERIALS 3, 041506 (2015)), which is widely adopted\nin thermoelectric research. This relation provides a more realistic estimation of the\nLorenz number compared with assuming a constant Sommerfeld value. We have added\nsome description about Lorenz number in the revised manuscript (see page 8, up).\n4. In Fig. 3b, are the defects indicated by the arrows nanopores?\nResponse: We thank the reviewer for her/his kind remind. In Fig. 3b, the defects\nindicated by the yellow arrows correspond to dislocations, while those marked by the\nred arrows represent nanoprecipitates. These defects arise during prolonged sintering\nand contribute to phonon scattering through dislocation (DC) and strain (DS)\nmechanisms, as well as nanoprecipitate scattering. Nanopores are also formed at longer\nholding times, and they are presented separately in Figure S7. We have also added some\ndescription in the figure caption in the revised manuscript (see page 10, down).\n5. The authors have employed the Debye-Callaway model for theoretical analysis of\ntheir results. This model primarily considers phonon-phonon Umklapp processes (U),\ngrain boundary scattering (B), point defect scattering (PD), dislocation scattering (DC),\nand strain scattering (DS).\nDoes it refer to mass fluctuations and strain field fluctuations caused by Bi atoms\nsubstituting Sb atoms? If so, how are Mg vacancies incorporated into this model?\nThe authors should clearly demonstrate how to calculate the mass fluctuations scattering\nparameter (Gamma_M) and strain field fluctuations scattering parameter (Gamma_S).\nResponse: We thank the reviewer for this insightful question. In the Debye-Callaway\nframework, the point-defect (PD) scattering term indeed includes both mass\nfluctuations and strain-field fluctuations arising from Bi atoms substituting Sb atoms.\nIn addition, Mg vacancies can be incorporated into the PD channel by treating them as\nmass-zero and radius-zero atoms on the Mg sublattice. The mass fluctuations ( ) and\nstrain-field fluctuations ( ) can be calculated by the following formulas, 𝛤𝛤𝑀𝑀\n𝛤𝛤𝐷𝐷\nPD 4\n−1 𝑉𝑉𝜔𝜔\n𝜏𝜏 = 3(𝛤𝛤𝑀𝑀 +𝛤𝛤𝐷𝐷)\n4𝜋𝜋𝜈𝜈\n2\n∑\n𝑛𝑛\n𝑖𝑖=1 𝑐𝑐𝑖𝑖\n𝑀𝑀𝑖𝑖\n Γ𝑀𝑀,𝑖𝑖 𝑘𝑘 2\n𝑀𝑀‾ 𝑘𝑘 𝑀𝑀𝑖𝑖\n𝛤𝛤𝑀𝑀 = 𝑛𝑛 ,𝛤𝛤𝑀𝑀,𝑖𝑖 = ∑𝑘𝑘 𝑓𝑓𝑖𝑖 1− \n∑𝑖𝑖=1 𝑐𝑐𝑖𝑖 𝑀𝑀𝑖𝑖\n𝑛𝑛\n∑𝑖𝑖=1 𝑐𝑐𝑖𝑖𝑀𝑀𝑖𝑖\n𝑀𝑀‾ = 𝑛𝑛\nwhere c is number of atoms on site i, ∑𝑖𝑖=1 𝑐𝑐a𝑖𝑖nd are the mass and fractional\ni\n𝑘𝑘 𝑘𝑘\noccupancy of the k-th atom, is the ave𝑀𝑀ra 𝑖𝑖 ge mass𝑓𝑓 𝑖𝑖 at site i, and is the average\nmass of the compound. 𝑀𝑀𝑖𝑖 𝑀𝑀‾\n2\n∑\n𝑛𝑛\n𝑖𝑖=1 𝑐𝑐𝑖𝑖\n𝑀𝑀‾𝑖𝑖\n 𝜀𝜀𝑖𝑖Γ𝑅𝑅,𝑖𝑖 𝑘𝑘 2\n𝑀𝑀‾ 𝑘𝑘 𝑅𝑅𝑖𝑖\n𝛤𝛤𝐷𝐷 = 𝑛𝑛 ,𝛤𝛤𝑅𝑅,𝑖𝑖 = ∑𝑘𝑘 𝑓𝑓𝑖𝑖 1− \n∑𝑖𝑖=1 𝑐𝑐𝑖𝑖 𝑅𝑅𝑖𝑖\n2\n2\n1+𝜈𝜈𝑝𝑝\n𝜀𝜀𝑖𝑖 = 6.4𝛾𝛾 \n9 1−𝜈𝜈𝑝𝑝\nHere, and are the radius and fractional occupancy of the k-th atom at site i,\n𝑘𝑘 𝑘𝑘\nis the a𝑅𝑅ve𝑖𝑖rage ra𝑓𝑓d𝑖𝑖ius at site i, is a coefficient determined by the Grüneisen paramet𝑅𝑅er 𝑖𝑖\nγ and the Poisson ratio . 𝜀𝜀𝑖𝑖\n𝜈𝜈𝑝𝑝\nobtained, and the recalculated , are 0.1117 and 0.3961 for the x = 0.01 (10 min)\nsample, and are 0.1120 and 0.4𝛤𝛤2𝑀𝑀27 𝛤𝛤f𝐷𝐷or the x = 0.01 (20 min) sample. In addition, Mg\nvacancies can aggregate during prolonged sintering, forming vacancy clusters,\ndislocations, and nanopores (Nature Communications 8, 13828 (2017)). In this way, the\ninfluence of Bi/Sb substitutional disorder and isolated Mg vacancies is consistently\nconsidered. Isolated point defects are included in the PD channel through and ,\nwhile vacancy clustering and extended defect formation are incorporat𝛤𝛤e𝑀𝑀d throu𝛤𝛤g𝐷𝐷h\ndifferent scattering channels. We have updated the Figure 2(c-d) and added some\ndescription in the Supporting Information (see page 17).\nFigure. (a) Experimental and calculated data of the samples after considering main\nphonon scattering, (b) the frequency-dependent κ for the x = 0.01 (20 min) sample at\ns\nroom temperature on the basis of the Debye-Callaway model.\n6. What physical mechanisms cause strain scattering (DS) in the Debye-Callaway\nmodel?\nResponse: We thank the reviewer for the insightful question. In the Debye-Callaway\nmodel, strain scattering (DS) refers to phonon scattering by the elastic strain fields of\ndislocations. These long-range strain fields distort the local lattice potential and alter\nphonon frequencies, which induce phonon scattering and lattice softening. The isolated\nMg vacancies can aggregate to form dislocations (Nature Communications 8, 13828\n(2017)). And the derivative defects of Mg vacancies are represented through different\nscattering channels. The DS contribution is distinct from the point-defect channel which\nalready accounts for the mass and size fluctuations of Bi/Sb substitutions and isolated\nand extended-scale effects of Mg vacancies are consistently considered. We have also\nadded some description in the revised Supporting Information (see page 18).\n7. What is the unit for N_D (Dislocation density) in Table S2? Is it cm-3 or cm-2?\nResponse: Yes, the dislocation density (N ) indicates the dislocation line density,\nD\nwhich is defined as the total length of dislocation lines per unit volume. Therefore, the\nunit for dislocation density (N ) is cm-2.\nD\n8. In the Abstract, the authors state, “Bonding rigidity degradation and phonon mean\nfree path suppression are dominantly governed by elevated internal strain, achieving\nlattice-constant-level phonon localization.”\nThis statement seems to imply that an increase in Mg vacancy concentration leads to\nincreased internal strain, reduced bonding rigidity, and suppression of the phonon mean\nfree path, suggesting that the decrease in lattice thermal conductivity is primarily due\nto the increase in Mg vacancy concentration.\nOn the other hand, the authors state, “Except for the vacancy-induced strain field decreasing κ ,\nlat\na certain amount of dislocations and nanopores form due to vacancy\nclustering during prolonged sintering. Numerous defects, including obvious\ndislocations and nanoprecipitates, have been found through the low-magnification\ntransmission electron microscopy images plotted in Figure 3(a)-(c).” (Line 221, page\n10)\nThis suggests that the decrease in lattice thermal conductivity is not only due to the\nincrease in Mg vacancy concentration, point defect scattering (PD), dislocation\nscattering (DC), and strain scattering (DS); nanopores and nanoprecipitates, which\nincrease with sintering time, also contribute. Can the Debye-Callaway model account\nfor the contributions of nanopores and nanoprecipitates? Does the theoretical analysis\nbased on the Debye-Callaway model indicate that the increase in Mg\nvacancy concentration is the primary cause of the decrease in lattice thermal conductivity?\nResponse: We thank the reviewer for her/his kind remind. We agree that our initial\n5\nexpression in the abstract is not sufficiently precise. The reduction in lattice thermal\nconductivity originates from Mg vacancies and their derived defects, which act at\ndifferent levels. The isolated Mg vacancies contribute at the atomic scale through point-\ndefect scattering. And during prolonged sintering, Mg vacancies tend to aggregate,\ngiving rise to vacancy clusters, dislocations, nanopores (Nature Communications 8,\n13828 (2017)), that act at larger length scales and further enhance phonon scattering.\nIn this sense, the internal strain generated by Mg vacancies and their derivative defects\nleads to lattice softening, while these defects suppress the phonon mean free path.\nTogether, they jointly constitute the main origin of the lattice thermal conductivity\nreduction. To provide a more comprehensive description, we have now explicitly\nincluded nanoprecipitate (τ -1 ) and nanopore (τ -1 ) scattering within the Debye-\nNP Pore\nCallaway framework. Their scattering rates are added following the formulas,\n-1 -1\n4\n-1\nτ -1= ν 2πR2 + πR2 ΔD D 2 ωR ν 4 N\nNP 9 P\n   ( ⁄ ) ( ⁄ )  \nwhere N is the number density of precipitates, R is the average radius for the\np\nprecipitates, D is the matrix density, and ΔD is the density difference between matrix\nand precipitates.\nν\n−1\n𝜏𝜏Pore = 3\n𝜋𝜋\n −1𝑙𝑙\nwhere l is the size of pores and p is the poro6si𝑝𝑝ty. These terms are incorporated into the\ntotal relaxation rate according to Matthiessen’s rule,\nτ -1\nNP\n−1 −1 −1 −1 −1 −1 −1\n𝜏𝜏𝑡𝑡𝑡𝑡𝑡𝑡 = 𝜏𝜏𝑈𝑈 + 𝜏𝜏𝑃𝑃𝐷𝐷 + 𝜏𝜏𝐵𝐵 + 𝜏𝜏𝐷𝐷𝐷𝐷 + 𝜏𝜏𝐷𝐷𝐷𝐷 + +𝜏𝜏Pore\nWe have updated the Figure 2(c-d) in the revised manuscript. This result confirms that\nMg vacancies and their derivative defects are the dominant cause of the lattice thermal\nconductivity reduction, with nanoprecipitates providing secondary contributions. We\nhave also added some description in the revised manuscript (see page 1, middle; page\n9, down) and Supporting Information (see page 19, up).\nphonon scattering, (b) the frequency-dependent κ for the x = 0.01 (20 min) sample at\ns\nroom temperature on the basis of the Debye-Callaway model.\n9. On Line 182, page 15 of the Supporting Information, \"k \" should be the Boltzmann\nB\nconstant.\nResponse: We thank the reviewer for her/his kind remind. We have corrected this note\nin the supporting information. (see page 19, down)\n10. In Fig. 5c, how was the contact resistance measured? What does the \"distance\" on\nthe X-axis represent?\nResponse: We thank the reviewer for her/his question. The contact resistance is\nobtained based on the four-probe method. As we can see the schematic image, when a\ncurrent is driven between A and B, the voltage drop is recorded between probes C and\nD. Probe C remains stationary, while probe D is shifted stepwise from one side into\nanother. Using the applied current (I ) and the measured voltage (V ), the resistance\nAB CD\nis obtained as R = V /I . As probe D traverses the thermoelectric leg, the\nCD CD AB\nrelationship can be expressed as follows.\nwhere R , R and R are resis𝑅𝑅ta 𝐷𝐷 n 𝐷𝐷 c=es 𝑅𝑅o 𝐷𝐷 f 1 co+n𝑅𝑅ta 𝑚𝑚 ct +lay𝑅𝑅e 𝐷𝐷 r 2 1, contact layer 2 and material,\nC1 C2 m\nrespectively. R and R were obtained from the intercept of the resistances versus\nC1 C2\nposition plot. The contact electrical resistivity (ρ ) is obtained from formula ρ = R ∙ ,\nc c C\nwhere S is the cross-sectional area of the thermoelectric leg. Accordingly, the “distanc𝑆𝑆e”\non the X-axis represents the distance between probes C and D. We have also added\nsome description in the revised Supporting Information (see page 11).\nfour-probe method.\n11. The authors need to provide a detailed explanation of how the COMSOL\ncalculations were performed, including the version used, the specific modules\nemployed, and any other relevant details.\nResponse: We thank the reviewer for her/his kind remind. The numerical simulation\nof power generation was performed with COMSOL Multiphysics software, version 6.1.\nThe Heat Transfer Module and the Electric Currents Module were coupled through the\nThermoelectric Effect multiphysics interface to simulate the coupled temperature and\nelectrical potential fields. A three-dimensional model is built in the software interface\nto represent thermoelectric legs with a similar geometry and dimension. The\ntemperature-dependent electrical conductivity, Seebeck coefficient, and thermal\nconductivity were taken from actual measurement results. Herein, the electrical and\nthermal contact resistances between interfaces were not considered in our simulation\nmodel. The material domains are discretized using a physics-controlled mesh, with the\nelement size set to the default (normal) setting, and a fixed temperature difference is\napplied between the hot and cold ends. We have also added some description in the\nrevised manuscript (see page 17, middle).\nThe author precisely controls the Mg vacancy concentration in the n-type Mg Sb Bi\n3 1.5 0.5\nsystem by regulating the holding time of SPS. modulating phonon dynamics through\nlattice softening and simultaneously suppressing the phonon mean free path in a more\nlocalized manner while remaining compositionally invariant. This strategy combined\nwith electron concentration optimization, yields one of the highest ZT of 2.06,\nmax\naccompanied by an extraordinary average ZT of 1.58 (300-723 K), which exceeds\nave\nother state-of-the-art n-type thermoelectric materials. In general, this is valuable\nresearch work, but the following issues need to be addressed before formal publication.\n1. Existing studies have shown that Te element doping is an effective means to convert\nMg-based systems into n-type. In this paper, the authors mixed MnTe with raw\nmaterials via ball milling, followed by sintering and forming through SPS. Then, after\nthe above experimental treatments, in what form does MnTe exist in the matrix? Does\nit form a second phase or is it completely doped into the lattice? If it is a second phase,\nwill the doping amount of MnTe increase with the extension of holding time, thereby\nexerting a significant impact on the performance? In addition, since the added content\nof MnTe is extremely low, which may be below the detection limit of XRD, it is\nsuggested to supplement other characterization methods such as scanning electron\nmicroscopy for further analysis.\nResponse: We thank the reviewer for her/his insightful question. Yes, due to the\ndetection limit of XRD, it needs more evidence to verify the MnTe form in the matrix.\nWe have performed SEM analysis on two representative samples, x = 0.01 (10 min)\nand x = 0.01 (20 min). In both samples, the elemental mapping shows that Mg, Sb, Bi,\nMn and Te are homogeneously distributed throughout samples. No obvious segregated\nMnTe compound is observed. And from the STEM analysis in Figure 3 in manuscript\nabout x = 0.01 (20 min) sample, the result also confirms that Mn and Te are uniformly\ndistributed without detectable MnTe compound. These results indicate that the added\nMnTe does not remain as a distinct second phase but is incorporated homogeneously\ninto the matrix during the ball milling or SPS process. Additionally, after prolonged\n9\nholding, no MnTe compound was detected, suggesting that the holding time used in\nthis work do not lead to noticeable second phase formation. We have also added some\ndescription in the revised manuscript (see page 7, down).\nFigure. The EDS elemental mapping of x = 0.01 (10 min) sample.\n-\nFigure. The EDS elemental mapping of x = 0.01 (20 min) sample.\n2. In Figure S3, the sample with x = 0.01 shows an obvious diffraction peak of the\nsecond phase around 25°, while this diffraction peak disappears when the holding time\nreaches 20 minutes. How do the authors explain this phenomenon?\nResponse: We thank the reviewer for her/his kind remind. We have carefully\nreexamined the diffraction peak around 25° in the x=0.01 sample. By comparing with\nthe standard diffraction data, we find that this peak corresponds to the carbon (PDF#26-\n10\n1076). Since this phase is not expected from the thermodynamic stability of the\nMg (Sb,Bi) system, we believe that the peak might be due to contamination during\n3 2\nmeasurement. To verify this, we repeatedly test the XRD on the same sample, and the\npeak at 25° completely disappeared, confirming that it is indeed caused by\nmeasurement contamination. We have revised Figure S3 in the Supporting Information\n(see page 3, down).\nFigure. The comparison of the XRD pattern of the x=0.01 sample with the standard\ndiffraction peaks of the carbon phase.\nFigure. The XRD patterns of (Mg\n3.2\nSb\n1.5\nBi\n0.5\n)\n1-x\n(MnTe)\nx\n(x = 0, 0.005, 0.01, 0.015, 0.02)\nsamples at room temperature.\n3. Prolonged holding time leads to significant changes in effective mass, which the\nauthors attribute to possible internal strain generated during the sintering process.\neffective mass by the generation of Mg vacancies or the dissolution of MnTe into the\nmatrix. It is suggested that the authors add EPMA compositional characterization to\nanalyze the actual composition of samples with different holding times.\nResponse: We thank the reviewer for her/his constructive suggestion. We have\nperformed EPMA compositional analysis on samples with different holding times. As\nwe can see, all elements are homogeneously distributed, and no MnTe compound is\nobserved. By the way, a slight decrease in Mg content is detected as the holding time\nenhancement, indicating an increased concentration of Mg vacancies. Although EPMA\nresults may contain certain uncertainties due to surface condition of the samples and\ninstrument precision, the overall trend of decreasing Mg content is consistently\nobserved. This phenomenon is consistent with the decrease of electron concentration.\nYes, the observed decrease of the density of states effective mass can be more\nrationalized by the formation of Mg vacancies. A lower electron concentration shifts\nfermi level downward, away from the regions of higher density of states, thereby\nreducing the effective mass. We have also added some description in the revised\nmanuscript (see page 13, middle).\nFigure. The EPMA elemental mapping of x = 0.01 (10 min) sample at room temperature.\nFigure. The EPMA elemental mapping of x = 0.01 (30 min) sample at room temperature.\nTable. EPMA mapping analysis results in the (Mg Sb Bi ) (MnTe) samples at\n3.2 1.5 0.5 0.99 0.01\ndifferent sintering condition.\nAtom% Mg Sb Bi Mn Te\nx=0.01 (10 min) 63.2 26.6 8.78 0.20 0.198\nx=0.01 (20 min) 61.5 30.2 7.97 0.20 0.160\nx=0.01 (30 min) 60.9 30.6 8.10 0.20 0.182\nIn this work, the authors tailored lattice vibrations via dynamically regulating Mg\nvacancy concentration, pronounced internal strain, and defect proliferation in n-type\nMg Sb Bi -based alloy, which softens interatomic bond stiffness and confines\n3 1.5 0.5\nphonon propagation to distances comparable to the lattice parameter. As a result, one\nof the highest ZT of 2.06, ZT of 1.04 (300-500 K) and unprecedented 1.58 (300-\nmax ave\nenable an excellent conversion efficiency of 12.5% (ΔT = 440 K) and 7.4% (ΔT = 300\nK) for a single TE-leg generator and a two-pair module. This work is suggested for\npublication after addressing some comments/questions below.\n1. Figure 2, the figure caption should be in alphabetical order. For example, “(e) Sound\nvelocities and (g) Raman spectrum for the (Mg Sb Bi ) (MnTe) samples at room\n3.2 1.5 0.5 1-x x\ntemperature. (f) The comparison of phonon mean free path, Debye temperature”\nResponse: We thank the reviewer for her/his kind remind. We have corrected the\ncaption of Figure 2 so that the subfigures are now described strictly in alphabetical order\nin the manuscript (see page 6, middle).\n2. TEM characterization. It is recommended to reformat the scale bar and label the type\nof defect/dislocation correspondingly. For example, in Figure 3g, the scale bar of GPA\nshould be -10% to 10%. Figure S6, the HRTEM characterization and line profile of\natomic intensity would be preferable for the examination of vacancy and strain variation.\nResponse: We thank the reviewer for her/his constructive suggestion. We have\nreformatted the scale bars of the TEM/GPA images, adjusted the GPA color scale to -\n10% to +10%, and labeled the type of dislocations directly in the figures. Yes, the line\nprofile of atomic intensity would be helpful for vacancy and strain analysis. However,\ndue to the imaging conditions of HRTEM, it is not good to provide reliable quantitative\ninformation in this time. In addition, the variations in carrier concentration and the\nWilliamson-Hall strain analysis consistently support our discussion on vacancy\nformation and strain effects.\n14\n3. The measurement uncertainty or error analysis should be provided to ensure the\nrobustness and reproducibility of TE properties (Figures in the main text and\nsupplementary info) for the reader.\nResponse: We thank the reviewer for her/his kind remind. We agree that providing\nmeasurement uncertainty is essential for demonstrating the robustness and\nreproducibility. The electronic and thermal transport measurements are independent of\nsample dimensions, and the uncertainties mainly reflect the instrumental deviations. In\nparticular, the systematic errors are approximately 3% for the Seebeck coefficient and\n5% for the electrical conductivity. For the total thermal conductivity, the uncertainty is\nabout 7% calculated from (1% for the density ρ, 5%\ndd tot d d Dp d D\n2 2 2\nfor the specific heat C a d n t d o t 5 = %  f ( o r d ) t h e + t h ( e D D p ) + ( D D ) ) . Consequently, the\np\npropagated error for the calculated ZT values is estimated to be approximately 10%.\nThe computational formula involved is . For\nd(ZT) d D d d d tot\n2 2 2\nclarity and readability of the figures, we d Z i Z s p l = a y  t h ( e 2 × e r r o D D r ) b + a r s ( i d n ) t h + e Z ( T d t o p t l o ) t , as it\nintegrates the uncertainties of all three properties and provides a clear indication of\noverall measurement reliability without overloading the individual transport-property\nplots with excessive visual information. We have also added some description in the\nrevised manuscript (see page 17, down).\n4. Figure 4c, it is suggested to format the title of figures (e.g., position, font, size) in the\nrevised submission.\nResponse: We thank the reviewer for her/his kind remind. We have reformatted the\nposition, font, and size of Figure 4c to ensure consistency and clarity (see page 12, up).\n5. It is impressive to benchmark the ZT in Figure 4f, while the relevant refs are required\nto be cited appropriately.\nResponse: We thank the reviewer for her/his constructive suggestion. We have added\nthe references in Figure 4f to ensure proper attribution (see page 12, up).\n6. Device section. The geometry optimization of TE legs and contact characterization\nare crucial for efficient power generation (Figures 5 and S17). Some refs are suggested\n15\nto be included in the device section (doi.org/10.1038/s41565-023-01457-5;\ndoi.org/10.1126/sciadv.adf5701; http://doi.org/10.1016/j.joule.2024.08.009). Also,\ndetailed modelling and contact characterizations are required to be supplemented in the\nmethods.\nResponse: We thank the reviewer for her/his kind remind. We have added the\nreferences in the device section. As for the geometry optimization modelling, the\nnumerical simulation of power generation was performed with COMSOL Multiphysics\nsoftware, version 6.1. The Heat Transfer Module and the Electric Currents Module were\ncoupled through the Thermoelectric Effect multiphysics interface to simulate the\ncoupled temperature and electrical potential fields. A three-dimensional model was\nbuilt in the software interface to represent thermoelectric legs with a similar geometry\nand dimension. The temperature-dependent electrical conductivity, Seebeck coefficient,\nand thermal conductivity were taken from actual measurement results. Herein, the\nelectrical and thermal contact resistances between interfaces were not considered in our\nsimulation model. The material domains were discretized using a physics-controlled\nmesh, with the element size set to the default (normal) setting, and a fixed temperature\ndifference was applied between the hot and cold ends.\nAs for the contact characterizations, the contact resistance is obtained based on the\nfour-probe method. As we can see the schematic image, when a current is driven\nbetween A and B, the voltage drop is recorded between probes C and D. Probe C\nremains stationary, while probe D is shifted stepwise from one side into another. Using\nthe applied current (I ) and the measured voltage (V ), the resistance is obtained as\nAB CD\nR = V /I . As probe D traverses the thermoelectric leg, the relationship can be\nCD CD AB\nexpressed as follows.\nwhere R , R and R are resis𝑅𝑅ta 𝐷𝐷 n 𝐷𝐷 c=es 𝑅𝑅o 𝐷𝐷 f 1 co+n𝑅𝑅ta 𝑚𝑚 ct +lay𝑅𝑅e 𝐷𝐷 r 2 1, contact layer 2 and mater(ia1l),\nC1 C2 m\nrespectively. R and R were obtained from the intercept of the resistances versus\nC1 C2\nposition plot. The contact electrical resistivity (ρ ) is obtained from formula ρ = R ∙ ,\nc c C\nwhere S is the cross-sectional area of the thermoelectric leg. We have also added som𝑆𝑆e\n16\ndescription in the revised manuscript (see page 17, middle) and Supporting Information\n(see page 11).\nFigure. Schematic diagram of a measurement of the contact electrical resistance by a\nfour-probe method.\n7. Figure 5c, the low contact resistivity achieved Mg (Sb, Bi) -based leg is impressive.\n3 2\nCould the authors elaborate more on the difference between contact resistivities\nobtained in the hot- (7.1 microohm square cm) and cold-side (10.2 microohm square\ncm)? Besides, in light of the minimization of contact resistivity and element oxidation,\nwhat are the selection criteria of barrier/coating materials for Mg (Sb, Bi) -based legs?\n3 2\nResponse: We thank the reviewer for her/his insightful question. Yes, there is a slight\ndifference between the contact resistances at both ends of the thermoelectric leg. The\npossible reasons for this may include subtle variations in local thermal and mechanical\nconditions during sintering, as the temperature gradient and pressure distribution may\nnot be fully identical, and the thermal stresses experienced at both ends could be slightly\ndifferent. However, the observed difference remains within the typical experimental\nuncertainty. Regarding the choice of barrier/coating materials, the main principles\ninclude: (i) chemical stability to suppress interfacial reactions and prevent the formation\nof resistive intermediate layers, (ii) work function matching to ensure ohmic contact\nand low interfacial resistivity, (iii) thermal expansion coefficient compatibility to\nminimize thermal stress during operation, (iv) low diffusivity to inhibit interfacial\nelement interdiffusion and (v) mechanical robustness to maintain long-term interfacial\nreliability under thermal cycling. According to Yin et al. (Acta Materialia 198, 25\n(2020)), the 304 stainless steel basically satisfies these requirements. In addition, the\n17\n304 stainless steel is also relatively inexpensive and readily available, and it has already\nbeen widely employed as a barrier material in Mg (Sb,Bi) -based thermoelectrics.\n3 2\nIn our work, by employing 304 SS as the barrier layer, we also achieved a low interfacial\nresistivity, further confirming its suitability.\n18"}, {"round": 2, "reviewer_comments": {"reviewer_1": "Thank you for addressing my concerns. I am satisfied with the clarifications provided, and no further questions remain.", "reviewer_2": "The author has adequately addressed my concerns; therefore, I recommend acceptance.", "reviewer_3": "Tha authors well addressed previous concerns with supplemented data and discussion, and this work is suggested for\npublication as it is."}, "author_response": ""}]

Decision letter