Metavalent alloying and vacancy engineering enable state-of-the-art cubic GeSe thermoelectrics

Domain: earth_science
Published: 2025-04-01
Editorial outcome: major_revision

Abstract
Conventional alloying strategies often require high alloying concentrations, leading to impurity phases and additional phase transition that limit the figure of merit of thermoelectric materials. Here, we introduce metavalent alloying and vacancy engineering as transformative strategies to facilitate the orthorhombic-to-cubic phase transition, in which we stabilize pure cubic GeSe under ambient conditions with just 10% alloying concentration using Sb2Te3 as an effective alloying agent. Compared to the covalently bonded orthorhombic phase, the metavalently bonded cubic GeSe features lower cation vacancy formation energy, reduced band gap, enhanced band degeneracy, weaker chemical bonding, stronger lattice anharmonicity, and multiple phonon scattering centers. These properties synergistically improve the power factor and suppress the lattice thermal conductivity. Subsequent trace Pb doping further reduces the lattice thermal conductivity, achieving an unprecedented ZT of 1.38 at 723 K in cubic (Ge0.95Pb0.05Se)0.9(Sb2Te3)0.1, along with a remarkable energy conversion efficiency of 6.13% under a 430 K temperature difference. These results advance the practical application of GeSe-based alloys for medium-temperature power generation and provide critical insights into the orthorhombic-to-cubic phase transition mechanism in chalcogenides.

Introduction
Thermoelectric energy conversion technology is crucial for addressing the global energy crisis by efficiently converting low-grade thermal energy into valuable electrical energy 1–5. The energy conversion efficiency (η) of thermoelectric devices is governed by the dimensionless figure of merit, ZT = S2σT/κ, where σ, S, κ, and T represent the electrical conductivity, Seebeck coefficient, thermal conductivity (comprising both carrier κe and lattice κL contributions), and absolute temperature, respectively 6–8. Enhancing thermoelectric performance typically involves band engineering 9,10 and phonon engineering 11–13. Band engineering focuses on optimizing the electronic band structure, such as through band convergence 14, band anisotropy 15, or the introduction of resonant energy levels 16,17, to improve the power factor (S2σ). In contrast, phonon engineering aims to reduce the κL by lowering the sound velocity (vm) or introducing multi-scale defects to shorten the phonon mean free path 18,19.

Recent advancements have brought layered orthorhombic SnSe into the spotlight due to its unparalleled ZT values 20,21. Similarly, GeSe, another IV-VI group semiconductor, has emerged as a promising candidate for medium-temperature thermoelectric applications 22,23. Density functional theory (DFT) calculations predict that optimally doped orthorhombic GeSe could achieve a peak ZT of 2.5 along the b-axis at 800 K 22. However, the limited solubility of dopants in orthorhombic GeSe restricts the effective carrier concentration (nH), leading to a low ZT of approximately 0.22 24. In contrast, transitioning to rhombohedral or cubic GeSe has demonstrated significantly enhanced ZT, driven by its unique metavalent bonding (MVB) mechanism 25–27. Unlike the conventional covalent bonding (CVB) in orthorhombic GeSe where two electrons are shared 28,29, resulting in a bond order of 1, MVB involves the sharing of approximately one electron between neighboring atoms, yielding a bond order of ~1/2. This bonding type, distinct from covalent, ionic, and metallic bonds, endows MVB solids with exceptional properties for carrier and phonon transport 28,30. Electrically, compared to its covalently bonded orthorhombic counterpart, metavalently bonded rhombohedral or cubic GeSe exhibits a reduced band gap (Eg), lower Ge vacancy formation energy, narrower band edges, and greater band degeneracy (Nv). Thermally, it demonstrates a lower vm and stronger lattice anharmonicity. These factors collectively contribute to the enhanced ZT value observed in metavalently bonded GeSe 31.

Current research on high-performance p-type GeSe primarily focuses on its rhombohedral phase 32. Significant advancements, such as Sb-Cd co-doping 33 and alloying with compounds like InTe 34, CdTe 27, MnCdTe 35, AgInTe2 36, LiBiTe2 37, and Ag(Sb, Bi)(Se, Te) 26,38–41, have successfully transformed orthorhombic GeSe into the rhombohedral structure, achieving ZT values ranging from 0.90 to 1.36 26,27,33–41. However, cubic GeSe, which is unaffected by Peierls distortion, offers larger band degeneracy, lower vm and greater lattice anharmonicity 25, theoretically surpassing the thermoelectric performance of the rhombohedral phase. Additionally, the phase transition from orthorhombic (or rhombohedral) to cubic at elevated temperatures causes a significant change in the coefficient of thermal expansion, leading to internal stresses that can damage the materials or disrupt the material-electrode interfaces. This, in turn, compromises the long-term reliability of thermoelectric devices. Therefore, facilitating the orthorhombic-to-cubic phase transition and stabilizing the cubic phase across the entire temperature range hold both theoretical and practical importance.

Understanding the mechanisms driving the orthorhombic-to-cubic phase transition in chalcogenides is therefore critical. In these materials, the orthorhombic structure exhibits strong Peierls distortion, characterized by high electron sharing (ES) between neighboring atoms, while the cubic structure, free from such distortion, displays minimal ES due to the absence of s- and p-band hybridization 28. By chemically or physically suppressing Peierls distortion (effectively reducing the ES capacity between neighboring atoms), the bonding can be transitioned from CVB to MVB, enabling the progression of GeSe from the orthorhombic to the rhombohedral, and ultimately, to the cubic phase 28,29. Strategic alloying of GeSe with compounds such as InSnTe2 and AgSnTe2 has been explored to facilitate this phase transition 42,43. However, limited solubility in these alloys has hindered the complete elimination of residual orthorhombic phases and the prevention of secondary phase precipitation, presenting challenges to achieving fully cubic GeSe.

Previous studies have investigated the critical role of chemical bonding mechanism in governing dopant solubility. High miscibility is achieved when both the host material and the dopant share the same MVB mechanism, as demonstrated in systems like SnTe-AgSbTe2 44 and GeTe-PbSe-PbS 45. In contrast, mixing metavalent compounds with covalent solids often leads to phase separation, as observed in the GeTe-GeSe-GeS system 46. Consequently, MVB compounds, such as IV-VI, V2VI3, and I-V-VI2, exhibit excellent solubility in GeSe, minimizing residual orthorhombic or precipitated phases (Fig. 1a). However, alloying cubic GeSe with IV-VI or I-V-VI2 compounds typically requires extremely high alloying concentrations, which can negatively impact μH and, consequently, the ZT value. For example, stabilizing cubic GeSe through I-V-VI2 alloying necessitates at a minimum of 30% content, significantly lowering the ZT compared to rhombohedral GeSe 38,43,47. Quantum chemical calculations and enthalpy of formation analyzes reveal that cation vacancies destabilize covalently bonded orthorhombic GeSe 48, but are essential for mitigating unfavorable anti-bonding Ge-Se interactions in metavalently bonded rhombohedral or cubic GeSe 49. This suggests that the introduction of cation vacancies is critical for stabilizing high-symmetry phases (Fig. 1a). Furthermore, DFT calculations and X-ray diffraction (XRD) studies indicate that the concentration of Ge vacancies in cubic GeTe is significantly higher than in its rhombohedral counterpart 50. These findings imply that alloying with V2VI3 compounds, which offer high solubility and additional cation vacancies, can accelerate the orthorhombic-to-cubic phase transition in chalcogenides, potentially bypassing the intermediate rhombohedral phase.

Results
According to the Hume-Rothery rules 51, the atomic size mismatch between Ge and Bi is significantly larger compared to that between Ge and Sb, resulting in a higher solubility of Sb2Te3 in GeSe than Bi2Te3 (Supplementary Fig. 1). Based on this principle, we introduced a small amount (>7.5%) of metavalent Sb2Te3 as an alloying agent to accelerate the stabilization of cubic GeSe (Fig. 1b). Compared to orthorhombic GeSe, cubic GeSe exhibits optimized nH, enhanced μH, and a higher density-of-states effective mass (m*), leading to a significant improvement in the S2σ. Additionally, the weakened chemical bonding and enhanced lattice anharmonicity in cubic (GeSe)1−x(Sb2Te3)x, combined with multiple phonon scattering mechanisms, effectively reduce the κL. To further minimize κL, trace amounts of Pb was doped into (GeSe)0.9(Sb2Te3)0.1, achieving a record-breaking ZT = 1.38 at 723 K in cubic (Ge0.95Pb0.05Se)0.9(Sb2Te3)0.1. The single-leg device demonstrated an impressive η of 6.13% under a ΔT of 430 K (Fig. 1c, d). These results not only represent a significant breakthrough in the ZT of cubic GeSe but also provide valuable insights into the key factors driving the orthorhombic-to-cubic phase transition in chalcogenides.

To determine whether the integration of non-stoichiometric metavalent Sb2Te3 can accelerate the stabilization of cubic GeSe, powder X-ray diffraction (XRD) was performed on the full series of (Ge1−yPbySe)1−x(Sb2Te3)x samples (x=0~15%, y=3~7%). As shown in Fig. 2a, the room-temperature XRD patterns reveal significant structural changes in GeSe with increasing Sb2Te3 content (x). The initial sample with x=0 exhibited diffraction peaks characteristic of the typical orthorhombic structure (Pnma). After alloying with Sb2Te3, the intensity of the orthorhombic phase peaks between 30° and 33° decreased significantly, while the intensity at 31.3° increased, indicating the formation of a cubic phase (Fm-3m). This suggests the coexistence of cubic and orthorhombic phases in the sample with x=0.05. When the Sb2Te3 content increased to x≥0.075, a predominantly cubic structure was consistently observed in all samples at room temperature. Notably, further addition of trace Pb to the x=0.10 sample had minimal impact on the phase structure, as shown in Supplementary Fig. 2.

To quantify the impact of Sb2Te3 alloying on the crystal structure of GeSe, structural refinement was conducted using the Fullprof Suite 52, allowing precise determination of the lattice parameters for the (Ge1−yPbySe)1−x(Sb2Te3)x samples (Supplementary Table 1 and Supplementary Figs. 3 and 4). Specifically, the x=0 sample exhibited a purely orthorhombic phase. As x increased, the proportion of the orthorhombic phase decreased significantly, while the fraction of the cubic phase correspondingly increased. For instance, in the x=0.05 sample, the ratio of cubic to orthorhombic phases was quantified as 55.84% to 44.16%. Interestingly, while the room-temperature XRD pattern of the x=0.075 sample suggested no detectable orthorhombic phase, structural refinement revealed a residual orthorhombic fraction of 6.36%. Complete transformation to a 100% cubic phase was achieved at x=0.10, emphasizing the critical requirement for Sb2Te3 content to exceed 7.5% to stabilize cubic GeSe at ambient conditions.

As shown in Fig. 2b and Supplementary Fig. 5, the lattice parameters (a, b, and c) of both orthorhombic and cubic GeSe increase significantly with the rising values of x or y. This lattice expansion is attributed to the larger ionic radii of Sb (145 pm), Pb (175 pm), and Te (221 pm) compared to those of Ge (73 pm) and Se (198 pm), respectively 53. More importantly, as shown in Fig. 2c and Supplementary Fig. 6, the bond lengths in orthorhombic GeSe generally range from 2.52 to 2.61 Å, whereas those in cubic GeSe extend from 2.84 to 2.89 Å. Based on the phenomenological bond length-bond strength relationship, shorter bonds typically correspond to stronger interactions, while longer bonds tend to indicate weaker interactions. Thus, the elongation of bond lengths in the cubic structure signifies a weakening of chemical bonding. This reduction in bond strength leads to weaker atomic interactions, which in turn reduces the vm and significantly lowers the κL.

To assess the effect of Sb2Te3 alloying on phase stability at high temperatures, temperature-dependent XRD measurements were conducted on the x=0.075 and x=0.10 samples (Fig. 2d, e). Upon heating, trace orthorhombic phases were observed in both samples, consistent with thermodynamic predictions that the metastable cubic phase transitions to the more stable orthorhombic structure under thermal stress. Notably, the orthorhombic phase appears above 623 K in the x=0.075 sample and above 673 K in the x=0.10 sample. Furthermore, the intensity of the orthorhombic peaks is significantly weaker in the x=0.10 sample, indicating better retention of the cubic structure at higher temperatures. While increasing the Sb2Te3 content could potentially suppress the formation of the orthorhombic phase at elevated temperatures, it would likely compromise the ZT value.

To validate whether metavalent Sb2Te3 alloying with additional cation vacancies can accelerate the formation of cubic GeSe, we compared the required alloying concentrations with conventional approaches (Fig. 1b). Traditional alloying to achieve a fully cubic structure typically demands a high alloying content, such as over 40% for AgSnTe2 43 and at least 30% for AgBiTe2 26, AgSbTe2 38, and AgBiSe2 47. In contrast, our study demonstrates that orthorhombic GeSe can successfully transition to a pure cubic phase with Sb2Te3 alloying concentrations as low as 10% or less. To further clarify the influence of different alloying types, we quantitatively evaluated the configurational entropy (ΔS) required for various cubic GeSe-based alloys using Boltzmann theory 54: ΔS = R Σni=1 xi ln xi, where R represents the gas constant, n denotes the number of element types, and xi specifies the mole fraction of the i-th element (Fig. 1b). Notably, alloying with AgSnTe2, AgBiTe2, and AgSbTe2 requires ΔS values exceeding the typical threshold for high-entropy alloys (ΔS=1.50R) 55. In the case of Sb2Te3 alloying, the cubic structure becomes mostly stable at ΔS=0.98R and fully stable at ΔS=1.21R. Compared to existing literature, the observed alloying content and ΔS in this study are significantly reduced, highlighting the importance of understanding the key factors governing the orthorhombic-to-cubic phase transition in chalcogenides. This reduced requirement underscore the efficiency and practicality of using Sb2Te3 as an alloying agent to stabilize cubic GeSe.

The driving force for the formation of cubic (GeSe)1−x(Sb2Te3)x alloys was investigated through the theoretical calculations of Ge vacancy formation energy (EV). As shown in Fig. 2f, we simulated three different crystal structures, namely Ge31Se32 (Pnma), Ge26Se27 (Fm-3m), and Ge22Sb4Se21Te6 (Fm-3m), and calculated EV using DFT. Notably, EV decreased significantly from 0.99 eV for binary orthorhombic Ge31Se32 to 0.35 eV for binary cubic Ge26Se27. For cubic structure, Ge22Sb4Se21Te6, EV was further reduced to 0.07 eV. Given that cubic GeSe typically require extreme conditions, such as high temperatures (>900 K) or high pressures (>8 GPa) 56,57, our results confirm that metavalent Sb2Te3 alloying can accelerate the formation of the cubic phase under ambient conditions by intentionally introducing cation vacancies. This highlights the role of Sb2Te3 in reducing vacancy formation energy, enabling the stabilization of cubic GeSe without the need for extreme synthesis conditions.

In addition to theoretical calculations, cation vacancy concentration can be quantitatively evaluated through XRD refinement by analyzing the atomic occupancy of cation sites 52 (Fig. 2g). Ideally, the atomic occupancy of cation sites should be 100%; however, with increasing Sb2Te3 content (x), this occupancy gradually deviates from 100%, indicating a significant increase in cation vacancy concentration. For instance, the cation site occupancy decreased sharply from 100% in the x=0 sample to 96.8% in the x=0.10 sample and further dropped to 93.2% in the x=0.15 sample. This trend demonstrates a direct correlation between the increasing Sb2Te3 alloying content and the generation of cation vacancies. In contrast, the addition of Pb to the x=0.10 sample increases the cation site occupancy, thereby reducing the concentration of cation vacancies, as shown in Supplementary Table 2.

The solubility of dopants, in addition to the cation vacancy, plays a pivotal role in facilitating the rapid formation of cubic GeSe-based materials. The powder XRD patterns of the (Ge1−yPbySe)1−x(Sb2Te3)x series (Fig. 2a and Supplementary Fig. 2), alongside scanning electron microscopy (SEM) combined with energy-dispersive X-ray spectroscopy (EDS) (Supplementary Figs. 7–12), confirmed the absence of any impurity phases. Additionally, X-ray photoelectron spectroscopy (XPS) measurements confirmed the valence states of Ge, Pb, Sb, Se, and Te as 2+, 2+, 3+, 2−, and 2−, respectively (Supplementary Fig. 13). To further investigate the chemical composition and distribution in the (GeSe)0.9(Sb2Te3)0.1 sample, transmission electron microscopy (TEM), known for its high chemical sensitivity and spatial resolution, was employed. Figure 3a, obtained by high-angle annular dark-field (HAADF) imaging, reveals a smooth surface with no evidence of secondary phases. EDS mapping further supports this observation, showing a uniform distribution of the matrix elements Ge and Se, as well as the dopants Sb and Te, throughout the sample. These findings, combined with XRD and SEM results, confirm the effective incorporation of metavalent Sb2Te3 into the GeSe matrix. A similar absence of impurity phases was observed in other metavalent Ag-V-VI2 alloyed GeSe systems 38,39. In contrast, when GeSe was alloyed with compounds possessing different bonding types (e.g., covalent, ionic, or metallic bonding), such as InSe 58, InTe 34, CdTe 27, MnCdTe 35, AgInTe2 36, and AgSnTe2 43, secondary phases like In2Se3, GeIn2Se4, MnCdSe2, CdSe, and Ag2Te were formed. This highlights the superior solubility of MVB compounds, such as Sb2Te3, in GeSe, which promotes the full stabilization of the cubic phase while effectively preventing the formation of undesired secondary phases. Thus, alloying with high-solubility, non-stoichiometric MVB compounds, combined with the deliberate introduction of cation vacancies, is identified as a critical factor in accelerating the orthorhombic-to-cubic phase transition in chalcogenide materials. This approach highlights the synergy between chemical design and defect engineering in stabilizing the cubic phase under ambient conditions.

The high-resolution TEM (HRTEM) image in Fig. 3b further confirms the cubic structure of (GeSe)0.9(Sb2Te3)0.1. A magnified view of the white boxed region in Fig. 3b reveals a well-defined atomic arrangement, where alternating Ge and Se atoms occupy the lattice positions. Due to the higher atomic number of Se (34) compared to Ge (32) 53, Se anions appear with brighter contrast, as illustrated in Fig. 3c. The observed atomic structure aligns perfectly with the atomic-level crystal structure model of cubic GeSe. To further verify the crystal structure, a fast Fourier transform (FFT) analysis was performed to identify the diffraction pattern. As shown in Fig. 3d, the FFT pattern corresponds to a cubic GeSe phase projected along the [010] direction. This analysis confirms that Sb2Te3 alloying successfully induces the formation of cubic GeSe.

Typically, rhombohedral and cubic GeSe exhibit lower Ge vacancy formation energy compared to the orthorhombic phase, creating a more favorable environment for Ge precipitation and resulting in Ge vacancy defects 35,48. Specifically, in rhombohedral GeSe systems, Ge vacancy planes are frequently observed, primarily due to larger atoms substituting smaller Ge atoms, leading to a pinning effect on the cation vacancy migration process 34. In this work, HRTEM analysis was conducted to explore structural defects in cubic GeSe. As shown in Fig. 3e, a high density of stripe-like defects was observed in (GeSe)0.9(Sb2Te3)0.1, consistent with previously reported characteristics of Ge vacancy planes 34. Figure 3f provides a magnified view of the boxed region in Fig. 3e, clearly showing weaker contrast for vacancy planes compared to lattice positions occupied by atoms. To further characterize these defects, Fig. 3g illustrates the line intensity distribution corresponding to Se layers, Ge layers, and Ge vacancy stripe layers. A significant intensity reduction (~35~40%) is observed for the vacancy layers, directly visualizing the presence of Ge vacancy defects. These Ge vacancy defects effectively scatter mid-frequency phonons, thereby suppressing the κL, contributing to the improved thermoelectric performance of the material.

To investigate the effect of Sb2Te3 alloying on the electrical transport properties, we used DFT to calculate the electronic band structures of orthorhombic GeSe, cubic GeSe, and Sb2Te3-alloyed cubic GeSe. Figure 4a–c show band structures of the orthorhombic GeSe, cubic GeSe, and cubic (GeSe)0.9(Sb2Te3)0.1, and Fig. 4d-f show corresponding partial density-of-states (PDOS). Orthorhombic GeSe exhibits an indirect Eg of 0.79 eV (Fig. 4a), where the relatively large Eg leads to low nH and limits the potential for achieving a high S2σ (Fig. 4d). In contrast, pure cubic GeSe shows a direct Eg of 0.05 eV along the L-symmetry direction (Fig. 4b), making it highly susceptible to minority carrier excitation and, consequently, to bipolar conduction effects (Fig. 4e). Upon Sb2Te3 alloying, the Eg widens significantly to 0.33 eV (Fig. 4c), effectively mitigating the detrimental impact of bipolar conduction (Fig. 4f). This Eg optimization highlights the critical role of Sb2Te3 alloying in improving the electronic properties of cubic GeSe and enhancing its thermoelectric performance.

Additionally, in orthorhombic GeSe, the valence band maximum (VBM) is located along the Γ-Y direction, exhibiting a valley degeneracy (Nv) of 2 (Supplementary Fig. 14a). By increasing the doping level and shifting the Fermi level deeper into the valence band, a second maximum can emerge at the Γ point, increasing Nv to 3. However, due to the inherently low solubility of dopants in covalently bonded GeSe, achieving such high nH experimentally remains challenging 24. In contrast, cubic GeSe, with its higher symmetry Brillouin zone, has its VBM located at the L point, resulting in a high Nv of 4 (Supplementary Fig. 14b). At higher doping levels, the Σ band, which exhibits an Nv of 12, can also contribute to conduction, leading to a total Nv of 16 in cubic GeSe. This behavior is similar to that observed in cubic PbTe and GeTe 59,60. However, in binary cubic GeSe, the significant energy separation (ΔE=0.39 eV) between the first and second valence band maxima limits the contribution of deeper valence bands to electrical transport. In this context, Sb2Te3 alloying plays a crucial role in reducing ΔE, facilitating the convergence of multiple valence bands. This promotes substantial orbital degeneracy and enhances the m*, thereby improving the thermoelectric performance of cubic GeSe.

To elucidate the impact of non-stoichiometric metavalent Sb2Te3 alloying on the electrical transport properties of all (Ge1−yPbySe)1−x(Sb2Te3)x samples, Fig. 5a presents the room-temperature nH and μH. Initially, pure orthorhombic GeSe exhibited an extremely low nH of 1.2 × 1016 cm−3, attributed to its relatively wide Eg and sparse Ge vacancies. Furthermore, the limited solubility of dopants in orthorhombic GeSe hindered the effective tuning of the Fermi level required to optimize nH. Metavalent Sb2Te3 alloying rapidly converted the CVB-dominated orthorhombic GeSe into MVB-connected cubic GeSe by deliberately introducing additional cation vacancies. This structural transition led to a dramatic increase in room-temperature nH, skyrocketing from 1.2 × 1016 cm−3 in binary orthorhombic GeSe to 2.7 × 1020 cm−3 and 5.6 × 1020 cm−3 in the x=0.075 and x=0.15 samples, respectively. This significant increase arose from the phase transition from the orthorhombic to the cubic structure, which considerably narrowed the Eg. Moreover, unlike materials dominated by strong CVB interactions with few cation vacancies, MVB materials inherently generate cation vacancies to eliminate energetically unfavorable anti-bonding interactions 31, leading to a four-order-of-magnitude increase in nH. Additionally, the larger atomic radius of Pb compared to Ge suppresses the formation of Ge vacancies, resulting in a slight reduction in nH in Pb-containing samples.

According to the Kane band model 61, the effective mass of the band edge carriers (m*), which is inversely proportional to the Eg, plays a critical role in determining electrical transport properties. In orthorhombic GeSe, the larger Eg results in flatter valence band edges, which are associated with a higher m*. This leads to a significantly lower μH at room temperature, measured at 1.2 cm2 V−1 s−1 35. Alloying with Sb2Te3 induces a dual transformation effect. For one thing, Sb2Te3 alloying enhances the crystal symmetry and significantly reduces Eg, leading to a lower m* and a rapid increase in μH. For another, the introduction of Sb2Te3 and Pb introduces a substantial number of substitutional point defects and cation vacancies, which enhance the point defect scattering on the carriers. These opposing effects resulted in an initial sharp increase in room-temperature μH for the cubic x=0.10 sample, reaching 7.8 cm2 V−1 s−1, followed by a gradual decline to 2.3 cm2 V−1 s−1 in the cubic x=0.15 sample. This trend reflects the balance between improved band structure and increased point defect scattering in the alloyed system. Figure 5b presents the temperature-dependent σ of the (Ge1−yPbySe)1−x(Sb2Te3)x samples. Initially, σ increases significantly with the addition of Sb2Te3, followed by a gradual decrease as the Sb2Te3 content continues to rise. The initial increase in σ is directly linked to the substantial enhancements in both nH and μH. However, despite the continued rise in nH, the subsequent decrease in σ is primarily attributed to the reduction in μH. For example, at room temperature, σ increases sharply from 0.3 S m−1 for the orthorhombic x=0 sample to an impressive 3.3 × 104 S m−1 for the cubic x=0.075 sample, before slightly decreasing to 2.0 × 104 S m−1 for the cubic x=0.15 sample. Additionally, pristine orthorhombic GeSe and low Sb2Te3-alloyed samples (0≤x<0.10) show an increase in σ with temperature, indicative of semiconductor-like transport dominated by carrier grain boundary scattering. In high Sb2Te3-alloyed samples (x≥0.10), the high nH reduces grain boundary scattering, leading to a decrease in σ near room temperature, reflecting metal-like behavior. However, the introduction of Pb slightly lowers nH, reinstating thermally activated σ and semiconductor-like behavior.

In contrast to the σ, the S initially undergoes a significant decrease, followed by an increase with increasing Sb2Te3 and Pb content, as shown in Fig. 5c. Specifically, at room temperature, S drops sharply from 775 μV K−1 in pristine GeSe to 130 μV K−1 for the x=0.05 sample, before increasing to 232 μV K−1 for the x=0.10, y=0.07 sample. The initial decrease in S for the Sb2Te3-alloyed samples is primarily attributed to the dramatic increase in nH, which rises by four orders of magnitude compared to the original orthorhombic GeSe. This significant increase in nH reduces the thermoelectric voltage generated per unit temperature gradient, leading to the observed reduction in S. Moreover, the S of pristine orthorhombic GeSe reaches its peak around 473 K before decreasing at higher temperatures, indicating the onset of a bipolar effect. This effect is intensified by its ultra-low nH, despite the large Eg, and severely diminishes S and consequently the ZT at elevated temperatures. However, Sb2Te3 and Pb alloying facilitate the formation of cubic GeSe, which significantly increases nH and shifts the onset temperature of the bipolar effect to higher values.

In addition to the nH, the S is closely tied to the m*, a key parameter influencing thermoelectric performance. Using the single parabolic band (SPB) model 62, we calculated the relationship between S and nH at room temperature, highlighting the significant impact of Sb2Te3 alloying on m* (Fig. 5d). As shown in Supplementary Fig. 14, orthorhombic GeSe exhibits a relatively low m* of 0.72 m0, attributed to its limited Nv=2. However, the incorporation of non-stoichiometric metavalent Sb2Te3 effectively mitigates Peierls distortion, facilitating the transition of GeSe to the cubic structure. This structural evolution leads to a substantial increase in Nv=4. Moreover, electronic band structure calculations reveal that Sb2Te3 alloying promotes multivalley band convergence, further enhancing m*. As a result, cubic (Ge0.95Pb0.05Se)0.9(Sb2Te3)0.1 achieves an exceptional m* of 4.82 m0, the highest reported value for rhombohedral and cubic GeSe-based alloys (Fig. 5e) 27,34,35,38,39,42,43,58. This significant increase in m* explains the observed improvements in thermoelectric performance in the Sb2Te3-alloyed samples.

Weighted mobility (μW) analysis provides deeper insights into the electronic band structure and scattering mechanisms within the material, offering a clearer evaluation of electron transport properties by considering the interplay between μH and m* 63. Figure 5f illustrates the variation of μW over the entire temperature range, calculated using the method described in the literature. Notably, the cubic (GeSe)0.9(Sb2Te3)0.1 sample exhibits the highest μW across the entire temperature range. Compared to pristine orthorhombic GeSe, the significant increase in μW highlights the effective manipulation of the electronic band structure facilitated by Sb2Te3 alloying. While Pb substitution leads to a slight reduction in μW, attributed to enhanced alloy scattering, particularly at ambient temperatures. The μW of cubic (Ge0.95Pb0.05Se)0.9(Sb2Te3)0.1 consistently remains above 40 cm2 V−1 s−1 at all temperatures. This robust performance underscores the tremendous potential of the material for further enhancing its ZT, demonstrating its viability for high-performance thermoelectric applications.

Figure 5g elucidates the relationship between S2σ and nH. The addition of Sb2Te3 causes a rapid increase in nH, bringing it closer to the optimal threshold. Notably, the actual nH of the x=0.10 sample aligns closely with the ideal value. This synergistic optimization of nH and μW significantly enhances S2σ across the entire temperature range, as depicted in Fig. 5h. Specifically, the cubic x=0.10 sample achieves a peak S2σ of 15.12 μW cm−1 K−2 at 723 K, which is approximately 32 times higher than that of pristine orthorhombic GeSe (0.46 μW cm−1 K−2). Although the introduction of trace Pb slightly reduces nH and μW, the cubic (Ge0.95Pb0.05Se)0.9(Sb2Te3)0.1 sample still demonstrates a highly competitive maximum S2σ of 13.52 μW cm−1 K−2. This exceptional performance makes it one of the most efficient GeSe-based thermoelectric materials reported to date 26,27,33–43,58,64, as illustrated in Fig. 5i, underscoring its significant potential for high-performance thermoelectric applications.

Non-stoichiometric metavalent Sb2Te3 alloying has a profound impact on the thermal transport properties, as shown in Fig. 6. Evidently, the κ decreases rapidly with increasing Sb2Te3 content near room temperature, and this trend becomes even more pronounced with subsequent Pb alloying (Fig. 6a). Although Pb alloying slightly reduces S2σ, it intensifies phonon scattering through point defects, further lowering κL, and also reduces κe due to a decrease in σ. Together, these effects result in a significant reduction in total κ. Specifically, the room-temperature κ gradually decreases from 1.76 W m−1 K−1 in pristine orthorhombic GeSe to 0.86 W m−1 K−1 for the cubic x=0.10 sample and further to 0.65 W m−1 K−1 for the cubic x=0.10, y=0.05 sample. This substantial reduction in κ demonstrates the effectiveness of Sb2Te3 and Pb alloying in suppressing thermal transport, which is crucial for enhancing the thermoelectric performance of GeSe-based materials.

The κL was calculated by subtracting the electronic contribution (κe) from the κ. κe was determined using the Wiedemann-Franz law 65, κe = LσT, where L is the Lorenz number estimated using the formula L=1.5 + exp(−|S|/116) (Supplementary Fig. 15). Figure 6b visually illustrates the effect of Sb2Te3 and Pb alloying on κL. Notably, near room temperature, κL decreases progressively with increasing Sb2Te3 and Pb content. Specifically, room-temperature κL significantly decreases from 1.71 W m−1 K−1 for pristine orthorhombic GeSe to 0.72 W m−1 K−1 for cubic x=0.10 sample, and further to 0.61 W m−1 K−1 for cubic x=0.10, y=0.05 sample. Remarkably, for the (Ge0.95Pb0.05Se)0.9(Sb2Te3)0.1 sample, the κL reaches a minimum of 0.42 W m−1 K−1 at 723 K, approaching the amorphous limit of GeSe 24. Figures 6c, d compare the κL values obtained in this study with literature-reported data, revealing that the current κL values are among the lowest for GeSe-based thermoelectric materials. These findings highlight the exceptional ability of Sb2Te3 and Pb alloying to suppress κL and improve thermoelectric performance 26,27,33–43,58,64. Furthermore, pristine orthorhombic GeSe exhibits a κL ~T−1.21 trend, aligning with acoustic phonon scattering (κL ~T−1.0). Alloying with Sb2Te3 and Pb reduces the power law exponent. Specifically, the cubic (Ge0.95Pb0.05Se)0.9(Sb2Te3)0.1 sample demonstrates a κL ~T−0.45 trend, consistent with dominant alloy scattering (κL ~T−0.5).

To explore the mechanisms contributing to the reduction in κL toward the theoretical minimum for GeSe, the following factors can be identified. Firstly, according to the phonon transport model, vm is a critical determinant of κL. vm is influenced by the interatomic force constant (F), which reflects the chemical bond strength, and the atomic mass (M) according to vm ∝ F/M 31. In orthorhombic GeSe, the strong CVB, where a pair of electrons is shared between neighboring atoms, results in rigid bonds and a higher vm. In contrast, cubic GeSe with MVB, characterized by single-electron sharing between neighboring atoms, exhibits weaker chemical bonds, leading to a lower vm. This reduction in bond strength is corroborated by the elongated bond lengths observed in XRD refinements (Fig. 2c) and the reduced bulk/Young’s modulus and Debye temperature listed in Supplementary Table 3, all indicating weaker bonds in cubic GeSe. Secondly, the replacement of Ge with heavier atoms like Sb and Pb, and Se with Te, further reduces vm. For example, the room-temperature vm decreases significantly from 2242 m s−1 in pristine orthorhombic GeSe to 2139 m s−1 in the cubic x=0.10 sample and further to 2073 m s−1 in the cubic x=0.10, y=0.05 sample (Fig. 6e). These mechanisms demonstrate how Sb2Te3 and Pb alloying effectively weaken the bonds and reduce the vm, thereby suppressing κL and enhancing thermoelectric performance.

Additionally, lattice anharmonicity plays a critical role in controlling intrinsic phonon scattering by the Umklapp process (U-process), which significantly affects κL 31. U-process scattering has been identified as the dominant phonon scattering mechanism in most thermoelectric materials 66. In contrast to the behavior of covalently bonded orthorhombic GeSe, the highly delocalized p-electrons in metavalently bonded cubic GeSe exhibit long-range interactions along the <100> direction. These interactions induce significant electron-transverse optical phonon coupling, which softens optical phonon modes and enhances lattice anharmonicity. The softened optical phonons strongly couple with heat-carrying acoustic phonons through anharmonic scattering, reducing the group velocity of acoustic phonons and expanding the phase space for Umklapp scattering processes 66. To quantify this effect, we calculated the Grüneisen parameter (γ), an indicator of lattice anharmonicity 31. Higher γ values indicate greater anharmonicity, which, in turn, results in lower κL. As shown in Fig. 6f, at room temperature, γ increases from 1.36 in orthorhombic GeSe to 1.67 in the cubic x=0.10 sample, and then gradually rises to 1.71 in the cubic x=0.10, y=0.07 sample. Additionally, we observe a slight decrease in γ (and a corresponding increase in vm) with higher Sb2Te3 contents (x>0.05), likely due to the influence of Sb2Te3 alloying on the phonon dispersion relationship. In summary, this increase in anharmonicity, driven by MVB and the softening of optical phonons, is a key factor in the further reduction of κL.

Additionally, Sb2Te3 and Pb co-alloying introduces substitution and vacancy point defects, disrupting the lattice periodicity and inducing mass and strain field fluctuations, which enhance phonon scattering. Therefore, a point defect scattering model was employed to simulate the impact of mass and strain contrast on κL 67,68, as illustrated in Fig. 6g. Specifically, we considered three types of point defects: Sb/Ge substitution, Te/Se substitution, and Ge vacancies. The simulation revealed that even with these combined scattering mechanisms, the experimentally observed κL in (GeSe)1−x(Sb2Te3)x is significantly lower than the theoretical predictions. This discrepancy suggests the presence of additional phonon scattering mechanisms beyond those considered in the point defect model. To further investigate the impact of multiscale microstructures on heat transport in the (Ge0.95Pb0.05Se)0.9(Sb2Te3)0.1 sample, the Debye-Callaway model was applied to calculate the spectral lattice thermal conductivity (κs) 69 (Fig. 6h). This model accounted for various phonon scattering mechanisms, including Umklapp and normal processes (U+N), grain boundaries, substitutions, Ge vacancies, and Ge vacancy layers (detailed calculation parameters are provided in Supplementary Table 4 of the Supporting Information). The area between the two κs curves corresponds to the reduction in κL caused by each specific phonon scattering mechanism: grain boundaries primarily scatter low-frequency phonons, substitutions and Ge vacancies primarily impede the high-frequency phonon transport, and planar defects significantly reduce κL by scattering mid-frequency phonons. In summary, the combination of weakened chemical bonding, strong lattice anharmonicity, and the presence of multiple phonon-scattering centers collectively reduces κL, bringing it closer to the theoretical minimum. These factors highlight the synergy between material design and multiscale defect engineering in optimizing thermoelectric performance. Besides, stabilizing cubic GeSe through non-stoichiometric metavalent Sb2Te3 alloying enables a synergistic optimization of electrical and thermal transport properties. This is clearly demonstrated by the correlation between the μW and the reciprocal of the κL (1/κL) at 300 K (Fig. 6i) 27,35,42 and the temperature-dependent quality factor (βB) (Supplementary Fig. 16) for (Ge1−yPbySe)1−x(Sb2Te3)x samples.

As shown in Fig. 7a and Supplementary Table 5, the cubic (Ge0.95Pb0.05Se)0.9(Sb2Te3)0.1 sample exhibits an optimized nH, enhanced μH, increased m*, and reduced κL, collectively leading to a maximum ZT (ZTmax) of 1.38 at 723 K, representing a remarkable 27-fold improvement compared to pristine orthorhombic GeSe, which shows a ZTmax of 0.05 at 673 K. This ZTmax is the highest reported to date for GeSe-based alloys, spanning orthorhombic, rhombohedral, and cubic structures (Figs. 1c, 7b) 24,26,33–38,40–43,64. Moreover, as illustrated in Fig. 7c, the ZTave values for samples with x=0.10 across the y series consistently exceed 0.60, highlighting their excellent thermoelectric performance over a wide temperature range. These results underscore the significant potential of Sb2Te3-alloyed cubic GeSe for high-performance thermoelectric applications. Although trace amounts of Pb are incorporated in this work, its usage is reduced by an order of magnitude compared to traditional PbTe-based alloys. Moreover, there is an inherent trade-off between material toxicity and achieving high ZT values, as superior performance is essential for thermoelectric applications.

The mechanical properties of thermoelectric materials have become increasingly important as they determine their processability, operational stability, and cost-effectiveness in thermoelectric devices. Sb2Te3 and Pb alloying significantly enhance critical mechanical parameters, such as fracture toughness (KIC), Vickers hardness (HV), and compressive strength (σbc). Pristine orthorhombic GeSe exhibits moderate mechanical properties, with KIC of 0.36 MPa m1/2, HV of 1.1 GPa, and σbc of 255 MPa (Fig. 7d–f). The introduction of Sb2Te3 leads to substantial improvements: for the cubic x=0.10 sample, KIC increases to 0.88 MPa m1/2, HV rises to 1.8 GPa, and σbc reaches a peak of 495 MPa. Further substitution with Pb results in slight additional enhancements, with KIC reaching 0.90 MPa m1/2, HV improving to 2.1 GPa, and σbc climbing to 555 MPa. Since metavalent bonding is inherently weaker and associated with a high concentration of cation vacancies, cubic GeSe is expected to exhibit lower HV and σbc than orthorhombic GeSe. However, alloying with Sb2Te3 and Pb introduces significant solid-solution strengthening, enhancing HV and σbc beyond those of the original orthorhombic GeSe. Although solid-solution strengthening can sometimes reduce toughness, MVB acts as an intermediate state between localized covalent and delocalized metallic bonding, enabling materials to absorb and distribute strain more effectively under stress. This mechanism mitigates brittle fracture and enhances toughness, explaining the superior KIC of Sb2Te3 and Pb co-alloyed cubic GeSe compared to its orthorhombic counterpart. This enhanced mechanical performance makes Sb2Te3- and Pb-alloyed cubic GeSe more suitable for practical thermoelectric applications, offering both improved durability and processability.

After achieving enhanced thermoelectric and mechanical performance, we fabricated a single-leg thermoelectric device using the optimized cubic (Ge0.95Pb0.05Se)0.9(Sb2Te3)0.1 sample and evaluated its performance with a custom-designed system. Figure 7g illustrates the current (I) and output voltage (V) characteristics of the single-leg device under various ΔT, with the cold-side temperature (Tc) maintained at 285 K. The internal resistance (Rin) of the device was calculated from the slope of the V-I curves, and the open-circuit voltage (VOC) increased significantly as the hot-side temperature (Th) rose from 523 to 723 K. Specifically, the VOC increased from 49.6 to 98.9 V, while Rin decreased from 1.71 to 1.39 Ω. Concurrently, the V decreased linearly with increasing I. Figures 7h, i show the output power (P) and η as functions of I under different ΔT. Both P and η initially increased with I, reaching their maximum values when the external load resistance matched Rin, and then decreased. Ultimately, the single-leg device achieved a maximum power density of 0.21 W cm−2 and a peak η of 6.13% under a ΔT of 430 K. The experimentally measured η remains significantly lower than the theoretical value of 8.36%, which can be attributed to several factors, including inconsistencies in the soldering process, heat losses within the testing apparatus, non-ideal characteristics of the thermoelectric model, dimensional inaccuracies, and measurement errors in the module 70. These factors highlight the need for further optimization of device fabrication and measurement techniques to fully realize the theoretical performance potential of the materials. Furthermore, the development of high-performance n-type GeSe-based alloys to complement their p-type counterparts is critical for fabricating efficient p-n legs in thermoelectric devices. This advancement represents a key area for future research and is essential for facilitating the industrial-scale application of GeSe-based materials.

In conclusion, we synthesized a series of (Ge1−yPbySe)1−x(Sb2Te3)x samples to explore the strategies for accelerating the orthorhombic-to-cubic phase transition in chalcogenides. Our findings reveal that only 10% metavalent Sb2Te3 alloying is sufficient to stabilize cubic GeSe under ambient conditions without the formation of detectable impurity phases. This stabilization is attributed to the high solubility of Sb2Te3 and its role in introducing additional cation vacancies. Compared to its covalently bonded orthorhombic counterpart, the metavalently bonded cubic GeSe exhibits a lower Ge vacancy formation energy, a narrower Eg, and higher valley degeneracy. The Sb2Te3 alloying further facilitates multivalley band convergence, enhancing orbital degeneracy, and optimizing nH, improve μH, and increase m*. These improvements lead to a significant enhancement in S2σ. Moreover, the weak chemical bonding, strong lattice anharmonicity, and multiple phonon-scattering centers in the cubic (Ge0.95Pb0.05Se)0.9(Sb2Te3)0.1 effectively suppress κL to 0.42 W m−1 K−1 at 723 K, approaching the amorphous limit of GeSe. As a result, cubic (Ge0.95Pb0.05Se)0.9(Sb2Te3)0.1 achieves an impressive ZT of 1.38 at 723 K and a peak η of 6.13% under a ΔT of 430 K in a single-leg thermoelectric device. These findings highlight the critical role of chemical bonding and vacancy engineering in facilitating the orthorhombic-to-cubic phase transition in chalcogenides and provide valuable insights for designing next-generation thermoelectric materials.

Methods
Synthesis procedure
All the GeSe-based samples were prepared using a systematic process involving melting, quenching, ball milling, and spark plasma sintering (SPS). High-purity element chunks of Ge, Se, Sb, Te, and Pb (≥99.999%, Beijing Yipin Chuancheng Technology Co., Ltd.) were weighed according to the stoichiometric compositions (GeSe)1−x(Sb2Te3)x (x=0.05, 0.075, 0.10, 0.125, 0.15) and (Ge1−yPbySe)0.9(Sb2Te3)0.1 (y=0.03, 0.05, 0.07). The weighted materials were placed into quartz tubes with a 20 mm diameter and sealed under a vacuum of 10−3 Pa. The encapsulated mixtures were heated in a muffle furnace to 1323 K for 8 hours, held at this temperature for 10 hours, and then cooled to 1123 K, followed by rapid quenching in water. The resulting ingots were pulverized using a ball mill (MSK-SFM-3, MTI Corporation) at 1000 rpm for 25 minutes under vacuum. The resulting fine powders were then consolidated using SPS at 723 K for 5 minutes under a uniaxial pressure of 50 MPa, producing high-density cylindrical samples with a 20 mm diameter, suitable for microstructural analysis and thermoelectric performance evaluation.

Phase and microstructure characterization
The phase purities and crystal structures of the samples were analyzed at room temperature using XRD (Cu-Kα radiation, SmartLab, Rigaku®, Japan). Lattice parameters for the various crystal structures were determined via Rietveld refinement using the Fullprof Suite software 52. The uncertainty in the XRD structural refinement was quantified as ±0.002 Å. The chemical valence states were confirmed using XPS (Thermo ESCALAB 250, Japan) with monochromatic Al Kα radiation. Microstructural characterizations were conducted using SEM (Hitachi® SU-70, Japan) equipped with EDS, along with high-resolution TEM (HRTEM, FEI Titan Themis G2 300).

Thermoelectric performance measurement
The temperature-dependent σ and S were simultaneously measured using a ZEM-3 (Ulvac-Riko®, Japan) under a protective helium atmosphere, with the measuring error of 10% for σ and 7% for S. In this vein, the σ was measured using four-probe method and the S was measured by static direct-current methods. Thermal diffusivity (D) was determined by the laser flash method (LFA 467, Netzsch®, Germany), with the measure error of 3%. The sample density (d) was evaluated using the Archimedes method 71, while the specific heat capacity (Cp) was estimated using the Dulong-Petit law 72. The κ was then calculated using the formula κ=DdCp. The κe was computed using the Wiedemann-Franz law 73, κe = LσT, where the Lorenz number (L) was estimated using the formula L=1.5 + exp(−|S|/116) 74. Subsequently, the κL was derived by subtracting κe from the κ. The maximum measurement temperature for all samples was 723 K. The Hall coefficient (RH) was measured at room temperature using the physical properties measurement system (PPMS-9, Quantum Design®, America). The uncertainty of RH measurements was ±5%. nH and μH were calculated via the equation nH=1/eRH and μH=σRH, respectively. The longitudinal (vl) and transverse (vt) components of sound velocity were measured using an ultrasonic pulse receiver (MS-100, Ritec Laboratories, France).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "We have carefully reviewed the manuscript titled \"Metavalent Alloying and Vacancy Engineering Enable State-of-the-Art Cubic GeSe Thermoelectric\" by Haoran Luo et al. The manuscript offers important insights into the stabilization of cubic GeSe through Sb2Te3 and Pb alloying, along with its exceptional thermoelectric properties. While the study is promising and contributes significantly to the understanding of thermoelectric materials, we have identified several aspects that need further clarification and expansion to strengthen the manuscript before publication to Nature Communications. Below, we outline our detailed feedback and suggestions for improvement, concentrating on areas such as experimental validation, data interpretation, and an in-depth discussion of the material’s characteristics.\n• The introduction lacks a clear explanation of the significance of the orthorhombic-to-cubic phase transition.\n• Could the authors provide additional experimental evidence to support the phase stabilization mechanism of cubic GeSe?\n• A more thorough explanation of the impact of Pb alloying on the thermoelectric properties is needed.\n• How does the increase in valley degeneracy affect the overall thermoelectric performance of the material?\n• The manuscript would benefit from a clearer description of the computational methods used in determining the electronic structure.\n• The manuscript does not sufficiently address the potential challenges of scalability for industrial applications.\n• Can the authors provide more details about the processing methods used to fabricate the samples?\n• The discussion on lattice anharmonicity could be expanded to better explain its role in reducing κL.\n• A direct comparison between the thermoelectric properties of the cubic and orthorhombic GeSe phases would strengthen the manuscript.\n• The temperature-dependent behavior of the material should be analyzed more thoroughly to capture its full potential.\n• More clarification is needed on the experimental setup for measuring the thermoelectric properties.\n• The impact of Sb2Te3 content on the overall material stability at high temperatures could be discussed in more detail.\n• The authors should clarify the role of point defects and vacancies in enhancing thermoelectric performance.\n• A deeper discussion of the potential environmental impact of using Pb in the alloying process would be valuable.\n• Can the authors explain the computational techniques used for the electronic structure calculations and their accuracy in predicting thermoelectric properties?\n• How did the choice of computational models affect the electronic band structure predictions, and what limitations might exist?\n• Could the authors provide more details on the pseudopotentials and exchange-correlation functional used for the electronic structure calculations?\n• What are the potential sources of error in the electronic structure predictions, and how can they be mitigated in future studies?", "reviewer_2": "Title: Metavalent Alloying and Vacancy Engineering Enable State-of-the-Art Cubic GeSe Thermoelectrics\nThis study provides valuable insights into the effects of meta valent alloying on the thermoelectric of GeSe materials. The work well written and can be accepted after addressing given comments.\nComment 1:\nThere is an error in the caption for Figure S3, where the composition is incorrectly written as x = 0.005 instead of x = 0.05.\nCorrect this inconsistency throughout the manuscript to reflect the actual data.\nComment 2:\nIn Figure 2c, the bond lengths in orthorhombic GeSe are showed to range from 2.52 to 2.61 Å, whereas those in cubic GeSe are longer, extending from 2.84 to 2.89 Å. Based on the bond length-bond strength relationship, shorter bonds generally indicate stronger interactions and higher stability, while longer bonds suggest weaker interactions. Consequently, the orthorhombic phase appears more stable than the cubic phase due to stronger bonding. Please clarify how the x = 0.10 sample, with weaker bonding in the cubic phase, demonstrates higher thermal stability compared to the x = 0.075 sample.\nComment 3:\nIn Figure 2e, it is evident that the x = 0.01 sample reverts to the orthorhombic phase (peaks) at temperatures exceeding 673 K. Given this observation, the authors cannot use the word \"excellent thermal stability.\" Revise the statement.\nComment 4:\nIn Figure 1b, the evaluated ΔS values for x = 7.5% and 10% are shown as approximately 0.2 and 0.25 R, respectively. However, on page 9, line 178, the authors state ΔS values of 0.98 and 1.21 for these compositions. This discrepancy needs to be resolved. Please clarify which set of values is correct and provide a consistent explanation for the difference in reported data.", "reviewer_3": "The submitted manuscript, \"Metavalent Alloying and Vacancy Engineering Enable State-of-the-Art Cubic GeSe Thermoelectrics,” offers a compelling exploration of superior thermoelectric transport in cubic p-type GeSe systems, which could be of interest to the materials research community. The authors shed new light on the Sb2Te3 metavalent alloying and vacancy engineering as a strategy for orthorhombic to cubic phase transitions in GeSe with enhanced TE performance. The manuscript is well-written and includes extensive analysis.\nWhile the manuscript shows promise, some aspects require revisions to enhance its impact. Before considering publication in Nature Communications, I recommend the authors address the following comments.\n1. Phase formation\n• The authors must report lattice parameters and bond lengths with proper error bars in the graphs (Fig 2 (b), (c), S5 and S6).\n• Authors discussed the change in site occupancy in alloyed samples, evaluated from Rietveld analysis. Did authors observe such variations due to Pb incorporation as well? Please list essential refined parameters (R factors, site occupancy, Caglioti parameters, Biso etc.) with standard deviations in supplementary information to support the arguments.\n• Do authors have any experimental data other than EDS, to support the existence of Sb2Te3 as a minor matrix? Furthermore, did the authors experimentally verify the ionic states of all the elements, especially Sb and Te using techniques such as X-ray photoelectron spectroscopy? Please justify.\n2. Thermoelectric transport\n• Temperature dependency of electrical conductivity (σ) found to be different for different x values (Fig. 5(b)), for instance, Pb doped samples (y = 0.03, 0.05, 0.07 with x=0.1) show a strict increase in σ with temperature as compared to other samples. Please give a brief regarding the temperature dependency of σ.\n• Interestingly, the temperature dependency of weighted mobility of alloyed/doped samples are different than the pristine sample. Please describe this in terms of carrier scattering with temperature. Please mention the equation used in estimating weighted mobility in supplementary details.\n• In general, the heat capacity Cp values estimated by Dulong - Petit’s law will underestimate the thermal conductivity values at high temperatures. I suggest the authors refer to an empirical equation for Cp from the article- Wang, H., Gibbs, Z. M., Takagiwa, Y., & Snyder, G. J. (2014). Tuning bands of PbSe for better thermoelectric efficiency. Energy & Environmental Science, 7(2), 804-811. Please compare the values calculated from both the methods and use only the appropriate one in this case. Justify it in reply letter.\n• In Figure 6f, the Grüneisen parameter shows an increase with x value up to an optimum, beyond which it decreases. It would be helpful to the readers if the authors can provide a brief regarding this.\n• Even though the authors plotted the variation of μW versus 1/κL, it would be better to mention the values of the thermoelectric quality factor as well.\n• Please mention the uncertainty in the electrical measurements. Propagate the error in estimating ZT values.\n3. Mechanical Properties\n• Can substantial improvement in mechanical properties of alloyed/doped samples are also connected to the improved microstructure? Please provide a discussion on this.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "Response to reviewers’ comments\nReviewer #1\nWe have carefully reviewed the manuscript titled \"Metavalent Alloying and Vacancy Engineering\nEnable State-of-the-Art Cubic GeSe Thermoelectric\" by Haoran Luo et al. The manuscript offers\nimportant insights into the stabilization of cubic GeSe through Sb Te and Pb alloying, along with its\n2 3\nexceptional thermoelectric properties. While the study is promising and contributes significantly to the\nunderstanding of thermoelectric materials, we have identified several aspects that need further\nclarification and expansion to strengthen the manuscript before publication to Nature Communications.\nBelow, we outline our detailed feedback and suggestions for improvement, concentrating on areas such\nas experimental validation, data interpretation, and an in-depth discussion of the material’s\ncharacteristics.\nAutor reply: We sincerely appreciate the positive feedback and constructive suggestions, which are\ninvaluable for improving the quality of our manuscript.\nComment 1: The introduction lacks a clear explanation of the significance of the orthorhombic-to-\ncubic phase transition.\nAutor reply: Thank you for this valuable comment. Orthorhombic GeSe, characterized by covalent\nbonding, is theoretically expected to exhibit high thermoelectric performance (ZT) due to its favorable\nelectronic structure and inherently low lattice thermal conductivity. However, its practical ZT value is\nlimited to approximately 0.2, mainly due to its high band gap, large Ge vacancy formation energy, flat\nband edge, and low band degeneracy. In contrast, cubic GeSe, which features unique metavalent\nbonding, exhibits superior carrier and phonon transport properties. Electrically, cubic GeSe has a\nreduced band gap, lower Ge vacancy formation energy, narrower band edge, and higher band\ndegeneracy. These factors enhance carrier concentration, carrier mobility, and the density-of-states\neffective mass, ultimately improving the power factor. Thermally, cubic GeSe benefits from greater\nlattice anharmonicity and a lower sound velocity, which reduce lattice thermal conductivity and further\nenhance ZT. Additionally, the orthorhombic-to-cubic phase transition occurs at elevated temperatures\nand is accompanied by a significant change in the coefficient of thermal expansion. This phase\ntransition can induce internal stresses, potentially damaging the material or disrupting the interface\nbetween the material and electrodes, thereby compromising the long-term reliability of thermoelectric\ndevices. Promoting the orthorhombic-to-cubic phase transition and stabilizing the cubic phase across\nall temperatures is therefore of significant theoretical and practical importance. This approach offers a\ndual benefit: enhancing the ZT and improving the durability of GeSe-based thermoelectric devices.\nThese details have been included on Pages 3 and 4 of the revised manuscript.\nComment 2: Could the authors provide additional experimental evidence to support the phase\nstabilization mechanism of cubic GeSe?\nAutor reply: The phase stabilization mechanism of cubic GeSe is supported by extensive theoretical\nstudies, experimental data, and findings from our previous research. The orthorhombic-to-cubic phase\ntransition in GeSe is driven by a shift in chemical bonding, transitioning from covalent bonds with\nsignificant Peierls distortion in the orthorhombic phase to metavalent bonds with minimal or no Peierls\ndistortion in the rhombohedral or cubic phases. The key to promoting this transition is suppressing\nPeierls distortion through chemical or physical pressure, which alters the bonding mechanism of GeSe.\nTo accelerate the orthorhombic-to-cubic phase transition, we propose two guiding strategies: (1)\nalloying with metavalently bonded compounds and (2) introducing additional cationic vacancies. The\nsolubility of dopants, governed by the bonding mechanisms, ensures high miscibility when both the\nhost material and dopant exhibit metavalent bonding (MVB). Conversely, mixing metavalent\ncompounds with covalent or ionic solids typically causes phase separation. Thus, alloying GeSe with\nmetavalently bonded compounds is critical to stabilize the rhombohedral or cubic phase. For instance,\nalloying with InSe does not induce a phase transformation in GeSe, while InTe results in a coexistence\nof orthorhombic and rhombohedral phases but often introduces impurity phases (Figure R1a and b).\nNotably, only metavalently bonded compounds, such as IV-VI, V VI , and I-V-VI compounds, due to\n2 3 2\ntheir high solubility, can stabilize pure rhombohedral or cubic phases. However, alloying with these\ncompounds often requires high concentrations, which can reduce carrier mobility and ZT. For example,\nstabilizing cubic GeSe requires more than 30% AgBiSe , AgBiTe , or AgSbTe , and over 70% SnTe\n2 2 2\nor PbTe (unpublished data, Figure R1c and d). Quantum chemical calculations and enthalpy of\nformation analysis suggest that introducing additional cationic vacancies reduces unfavorable anti-\nbonding Ge–Se interactions, accelerating the formation of cubic GeSe. Among metavalently bonded\nsolids, only V VI compounds, specifically Sb Te and Bi Te , effectively stabilize cubic GeSe at low\n2 3 2 3 2 3\nconcentrations.\nFigure R1. Room-temperature powder X-ray diffraction (XRD) patterns of (a) GeSe(InSe) 1, (b)\nu\nGeSe(InTe) 1, (c) (GeSe) (SnTe) and (d) (GeSe) (PbTe) .\nv 1-x x 1-x x\nOur experimental results confirm these predictions. Alloying with just over 7.5% Sb Te stabilizes\n2 3\ncubic GeSe at room temperature, while 10% Bi Te also yields the cubic phase (Figure S1). However,\n2 3\nthe large atomic size mismatch between Bi and Ge leads to BiSe precipitation, making Sb Te the\n2 3\npreferred alloying agent to avoid secondary phase formation. These experimental findings strongly\nsupport the proposed mechanism for cubic GeSe stabilization.\nComment 3: A more thorough explanation of the impact of Pb alloying on the thermoelectric\nproperties is needed.\nAutor reply: Pb substitution for Ge is an isoelectronic replacement and, to a first approximation, does\nnot significantly alter the carrier concentration. However, due to the larger atomic radius of Pb\ncompared to Ge, Pb alloying suppresses the formation of Ge vacancies, resulting in a slight reduction\nin carrier concentration. Additionally, the increased scattering of carriers by point defects caused by\nPb alloying leads to a moderate decrease in carrier mobility. These combined effects slightly reduce\nthe electrical power factor. For example, at 723 K, the power factor decreases from 15.12 μW cm−1\nK−2 for the x = 0.10 sample to 13.52 μW cm−1 K−2 for the x = 0.10, y = 0.05 sample.\nIn contrast, Pb alloying significantly reduces the total thermal conductivity of GeSe, thereby\nenhancing the ZT value. At room temperature, the total thermal conductivity decreases from 0.86 W\nm−1 K−1 for the x = 0.10 sample to 0.65 W m−1 K−1 for the x = 0.10, y = 0.05 sample. This reduction\noccurs through two primary mechanisms. First, since electronic thermal conductivity is directly\nproportional to electrical conductivity, the reduced carrier concentration and mobility from Pb\nsubstitution lower the electronic thermal conductivity. Second, the large atomic mass and radius\ndifference between Pb and Ge creates substantial mass and strain field fluctuations, which enhance\nphonon scattering by point defects and significantly reduce lattice thermal conductivity.\nWe have included this information on Pages 18 and 22 of the revised manuscript.\nComment 4: How does the increase in valley degeneracy affect the overall thermoelectric\nperformance of the material?\nAutor reply: The increase in valley degeneracy has a significant impact on the electrical transport\nproperties. The Seebeck coefficient (S) is inversely related to the carrier concentration (n ) but\nH\npositively correlated with the density-of-states effective mass (m*). Therefore, enhancing m* is crucial\nfor achieving a high S. In this study, the formation of cubic GeSe leads to a significant increase in n .\nH\nTo counteract the negative effect of high n on S, a substantial increase in m* is essential. The m* is\nH\ndetermined by the band valley degeneracy (N ) and the band effective mass (m *) according to the\nV b\nequation m* = (N )2/3m *. While a heavy m * contributes to a higher m*, it often results in a reduction\nV b b\nin carrier mobility (μ ). In contrast, increasing N can effectively enhance m* without significantly\nH V\naffecting μ , provided valley scattering is negligible. Therefore, increasing N is critical for enhancing\nH V\nm*, which, in turn, improves S and ultimately the ZT.\nComment 5: The manuscript would benefit from a clearer description of the computational methods\nused in determining the electronic structure.\nAutor reply: Density-functional theory (DFT) calculations were performed using the all electron\nprojected augmented wave (PAW) method, as implemented in the Vienna Ab initio Simulation Package\n(VASP) 2-7. The generalized gradient approximation (GGA) with the fully relativistic Perdew-Burke-\nErnzerhof (PBE) functional was employed to treat the exchange correlations 8. The valence wave\nfunctions were expanded in a plan-wave basis with a cut-off energy of 500 eV. All atoms were allowed\nto relax in their geometric optimizations until the Hellmann–Feynman force is less than 1×10–3 eV Å–\n1. The convergence criterion and the Monkhorst-Pack k-mesh adopted for ionic relaxation are 1×10–7\neV and 0.06π per Å, respectively. A denser 0.02π Å–1 Monkhorst-Pack k-mesh was adopted for\ncalculating density-of-state (DOS), and a line-mode k-path based on Brillouin path features indicated\nby the AFLOW framework was adopted for calculating band structures 9, 10. The spin-orbital coupling\nwas considered for Sb Te -doped GeSe because Te is a heavy element. The LOBSTER package was\n2 3\nused to analyze bonding characteristics based on the contracted Slater-type orbitals 11. The harmonic\nlattice dynamics were calculated using the frozen-phonon method based on the force constant matrix\nobtained from explicitly displacing atoms, as implemented in the Phonopy package 12, during which a\n2×2×2 supercell for orthogonal GeSe and a 3×3×3 supercell for cubic GeSe were constructed to reflect\nlong-range interactions.\nWe have added the “electronic band structure calculation” in the experimental section of the\nrevised Supporting Information.\nComment 6: The manuscript does not sufficiently address the potential challenges of scalability for\nindustrial applications.\nAutor reply: Many chalcogenide compounds, such as GeS, GeSe, SnS, and SnSe from the IV-VI\ngroup, adopt low-symmetry orthorhombic structures under standard conditions. As discussed in\nComment 1, understanding the orthorhombic-to-cubic phase transition is critical for optimizing\nthermoelectric performance and ensuring the longevity of thermoelectric devices. This study highlights\nmetavalent alloying and vacancy engineering as effective strategies to promote the orthorhombic-to-\ncubic phase transition. For example, incorporating 10% Sb Te into GeSe stabilizes the cubic phase\n2 3\nacross all temperatures, preventing phase transitions at elevated temperatures and significantly\nimproving the long-term reliability of thermoelectric devices. Additionally, co-alloying with Sb Te\n2 3\nand Pb enhances both electrical and thermal transport properties, achieving a record-high ZT value of\n1.38 for GeSe-based thermoelectrics. This corresponds to an energy conversion efficiency of 6.13%\nunder a temperature gradient of 430 K, demonstrating the potential of these materials for industrial\napplications. However, further research is needed to develop high-performance n-type GeSe-based\nalloys to complement the p-type material and enable efficient p-n leg fabrication for thermoelectric\ndevices. Addressing this challenge is essential for realizing the full industrial potential of GeSe-based\nthermoelectrics.\nThis information has been included on Page 28 of the revised manuscript.\nComment 7: Can the authors provide more details about the processing methods used to fabricate the\nsamples?\nAutor reply: All GeSe-based samples were prepared through a systematic process involving melting,\nquenching, ball milling, and spark plasma sintering (SPS). High-purity elements of Ge, Se, Sb, Te,\nand\nPb (99.999% purity, supplied by Beijing Yipin Chuancheng Technology Co., Ltd.) were weighed\naccording to the stoichiometric compositions of (GeSe) (Sb Te ) (x = 0.05, 0.075, 0.10, 0.125, 0.15)\n1−x 2 3 x\nand (Ge Pb Se) (Sb Te ) (y = 0.03, 0.05, 0.07). The weighted materials were loaded into quartz\n1−y y 0.9 2 3 0.1\ntubes (20 mm diameter) and sealed under a vacuum of 10−3 Pa. The sealed tubes were heated to 1323\nK in a muffle furnace and held for 10 hours, followed by cooling to 1123 K. The samples were then\nquenched in water to rapidly solidify the material. The resulting ingots were pulverized into fine\npowders using a ball mill (MSK-SFM-3, MTI Corporation) operated at 1000 rpm for 25 minutes under\nvacuum. The powders were then consolidated using SPS at 723 K for 5 minutes under a uniaxial\npressure of 50 MPa, yielding high-density cylindrical samples with a 20 mm diameter. These samples\nwere prepared for subsequent microstructural analysis and thermoelectric performance evaluations.\nDetails of this synthesis procedure have been expanded in the experimental section of the revised\nsupporting information to ensure a thorough understanding of the experimental methods.\nComment 8: The discussion on lattice anharmonicity could be expanded to better explain its role in\nreducing κ .\nL\nAutor reply: In thermoelectric materials, lattice thermal conductivity is primarily determined by\nphonon transport. Even in a perfect crystal without structural defects, lattice thermal conductivity\nremains finite due to anharmonic phonon scattering from Umklapp processes. In the spring-ball model,\nwhich describes the crystal lattice and its vibrations, a harmonic system has a displacement-\nindependent spring constant, meaning phonon interactions do not affect net heat transport. In contrast,\nan anharmonic system exhibits a non-linear force-displacement relationship, where the spring constant\nchanges with atomic displacement. This anharmonicity enables phonon-phonon interactions, altering\nphonon propagation velocity and reducing the net phonon flow, thereby lowering lattice thermal\nconductivity. The degree of anharmonicity is commonly quantified by the Grüneisen parameter (γ).\nFor processes dominated by acoustic phonon scattering, the lattice thermal conductivity can be\nexpressed as13:\nMv3  1 \n  m   (R1)\nL TV2/3N1/3\nApparently, a larger Grüneisen parameter (γ) corresponds to a lower lattice thermal conductivity.\nA previous study compared the lattice thermal conductivity of materials with metavalent and\ncovalent bonding 14. The results showed that metavalently bonded materials exhibit long-range\ninteractions along the <100> direction, in contrast to the shorter-range interactions in covalently\nbonded materials. These long-range interactions primarily arise from p-electrons, while the\npolarization of s-electrons remains short-range, typically limited to the fourth-nearest neighbor. In this\nstudy, orthorhombic GeSe is characterized by covalent bonding, with the presence of Ge s-states at the\nvalence band edge indicating s-p hybridization. This results in entirely localized electrons. In contrast,\ncubic GeSe features metavalent bonding, dominated by p-p orbital interactions with weak s-p\nhybridization. This bonding nature combines localized and delocalized electronic states. The highly\ndelocalized p-electrons in cubic GeSe enable long-range interactions, which induce significant\nelectron–transverse optical phonon coupling. This coupling softens the optical phonon modes and\nenhances lattice anharmonicity. The softened optical phonons strongly couple with heat-carrying\nacoustic phonons through anharmonic scattering, reducing the group velocity of acoustic phonons and\nexpanding the phase space for Umklapp scattering processes. As a result, cubic GeSe with metavalent\nbonding exhibits significantly lower lattice thermal conductivity compared to its orthorhombic\ncounterpart.\nThis information has been included on Page 24 of the revised manuscript.\nComment 9: A direct comparison between the thermoelectric properties of the cubic and orthorhombic\nGeSe phases would strengthen the manuscript.\nAutor reply: A direct comparison of key thermoelectric properties, including room-temperature carrier\nconcentration (n ), carrier mobility (μ ), effective mass (m*), electrical conductivity (σ), Seebeck\nH H\ncoefficient (S), total thermal conductivity (κ), lattice thermal conductivity (κ ), and the maximum ZT\nL\n(ZT ) value, has been presented in Table R1 (Supplementary Table 5). This comparison highlights\nmax\nthat cubic GeSe, compared to orthorhombic GeSe, exhibits optimized carrier concentration, higher\ncarrier mobility, increased density-of-states effective mass, and lower lattice thermal conductivity.\nThese improvements collectively contribute to enhanced thermoelectric performance.\nThis information has been included on Page 27 of the revised manuscript.\nTable R1 (Supplementary Table S5). The comparison of room-temperature phase structures, carrier\nconcentration (n ), carrier mobility (μ ), effective mass (m*), electrical conductivity (σ), Seebeck\nH H\ncoefficient (S), total thermal conductivity (κ), lattice thermal conductivity (κ ), and maximum ZT\nL\n(ZT ) for representative GeSe-based samples (Orth. is the abbreviation of orthorhombic and Mixture\nmax\nrepresents the mixed phase of orthorhombic and cubic).\nn μ m* σ S κ κ\nH H L\nPhase Sample ZT\nmax\n(cm−3) (cm2V−1s−1) (m ) (Sm−1) (μVK−1) (Wm−1K−1) (Wm−1K−1)\nMixture x = 0.05 2.1×1020 4.02 2.65 1.36×104 130 1.26 1.18 0.48\nx = 0.10 3.0×1020 5.34 4.42 2.55×104 158 0.86 0.72 1.14\nCubic\nx = 0.10,\n1.57×1020 4.54 4.82 8.51×103 213 0.65 0.61 1.38\ny = 0.05\nComment 10: The temperature-dependent behavior of the material should be analyzed more\nthoroughly to capture its full potential.\nAutor reply: We first examine the temperature-dependent Seebeck coefficient. As shown in Figure\n5c, the Seebeck coefficient of pristine orthorhombic GeSe initially increases to a maximum value\nbefore decreasing at higher temperatures, indicating a bipolar effect around 473 K. This behavior arises\nfrom the ultra-low carrier concentration in orthorhombic GeSe, despite its large band gap. The onset\nof bipolar conduction reduces the Seebeck coefficient and diminishes thermoelectric performance at\nelevated temperatures, limiting its suitability for mid-temperature applications. In contrast, alloying\nwith Sb Te and Pb promotes the formation of cubic GeSe, reducing the band gap and the formation\n2 3\nenergy of Ge vacancies. This enhances the carrier concentration and suppresses the bipolar effect, even\nwith the narrowed band gap, shifting the onset of the bipolar effect to higher temperatures.\nRegarding electrical conductivity (Figure 5b), pristine orthorhombic GeSe and low Sb Te -\n2 3\nalloyed samples (0 ≤ x < 0.10) exhibit increasing electrical conductivity with temperature, consistent\nwith semiconductor-like transport behavior. This thermally activated behavior is attributed to carrier\ngrain boundary scattering, as reported in the literature 15-17. For high Sb Te -alloyed samples (x ≥ 0.10),\n2 3\nthe increased carrier concentration reduces grain boundary scattering, leading to metal-like transport\nbehavior and a decrease in electrical conductivity at lower temperatures. The addition of Pb further\nreduces the carrier concentration, reintroducing thermally activated electrical conductivity. Future\nstudies focusing on grain size optimization, such as single-crystal growth, could mitigate grain\nboundary scattering and improve carrier mobility, electrical transport, and ZT values.\nIn terms of lattice thermal conductivity (κ ) of orthorhombic GeSe (Figure 6b), a T−1.21\nL\ndependence is observed, consistent with acoustic phonon scattering (κ ~ T−1.0). Alloying with Sb Te\nL 2 3\nand Pb decreases the power law exponent, indicating a change in the phonon scattering mechanism.\nSpecifically, the cubic (Ge Pb Se) (Sb Te ) sample exhibits a κ ~ T−0.45 dependence,\n0.95 0.05 0.9 2 3 0.1 L\nconsistent with alloy scattering (κ ~ T−0.5). This shift highlights the role of point defects in thermal\nL\ntransport.\nThese details have been included on Pages 19 and 23 of the revised manuscript to provide a\ncomprehensive understanding of the temperature-dependent electrical and thermal behavior of the\nGeSe-based samples.\nComment 11: More clarification is needed on the experimental setup for measuring the thermoelectric\nproperties.\nAutor reply: In this work, electrical and thermal transport properties were measured using commercial\nequipment. Electrical conductivity and the Seebeck coefficient were measured using a ZEM-3 (Ulvac-\nRiko®, Japan) under a helium atmosphere, with measurement errors of 10% for electrical conductivity\nand 7% for the Seebeck coefficient. Electrical conductivity was measured using the four-probe method,\nand the Seebeck coefficient was determined using the static direct-current method. The measurement\nmechanism is shown in Figure R2a. The thermal diffusivity was measured using the laser flash method\n(LFA 467, Netzsch®, Germany), with a measurement error of 3%. The measurement mechanism is\nillustrated in Figure R2b. Sample density (ρ) was determined using the Archimedes method, and the\nspecific heat capacity (C ) was estimated using the Dulong-Petit law. The total thermal conductivity\np\n(κ) was calculated using the formula κ = DρC . Electronic thermal conductivity (κ ) was calculated\np e\nusing the Wiedemann-Franz law, κ = LσT, where the Lorenz number (L) was estimated using L = 1.5\ne\n+ exp(−|S|/116). Lattice thermal conductivity (κ ) was derived by subtracting κ from κ.\nL e\nAll\nmeasurements were conducted up to a maximum temperature of 723 K. The Hall coefficient (R ) was\nH\nmeasured at room temperature using a physical properties measurement system (PPMS-9, Quantum\nDesign®, USA), with a measurement uncertainty of ±5%. Carrier concentration (n ) and carrier\nH\nmobility (μ ) were calculated using the equations n = 1/eR and μ = σR , respectively.\nH H H H H\nFigure R2. (a) ZEM-3 (Ulvac-Riko®, Japan) test principle schematic diagram. (b) Schematic diagram\nof LFA 467 (Netzsch®, Germany) test system.\nThe “thermoelectric performance measurement” section is included in the supporting information.\nAdditionally, we have supplemented details on the measurement errors and mechanisms in the revised\nversion.\nComment 12: The impact of Sb Te content on the overall material stability at high temperatures could\n2 3\nbe discussed in more detail.\nAutor reply: The Sb Te content plays a crucial role in the formation and stabilization of the cubic\n2 3\nGeSe phase. As shown in Figure 2a, a gradual phase transition from the orthorhombic to cubic phase\noccurs with increasing Sb Te content. Specifically, the cubic phase is formed upon alloying with 7.5\n2 3\nat% Sb Te , indicating that even a small amount of Sb Te can stabilize the cubic phase at room\n2 3 2 3\ntemperature. However, in-situ XRD measurements reveal additional diffraction peaks near 32.03° at\ngroup Pnma), suggesting a thermally driven phase transition from the cubic to the orthorhombic\nstructure.\nIn contrast, increasing the Sb Te content to x = 0.10 significantly enhances the phase stability of\n2 3\ncubic GeSe. As shown in Figure 2e, in-situ XRD measurements of the x = 0.10 sample indicate a\nstable, pure cubic phase across the temperature range of 300 to 623 K, with only a trace amount of the\northorhombic phase appearing above 673 K. This represents an increase of 50 K in the temperature at\nwhich the orthorhombic phase reappears. Additionally, the intensity of the orthorhombic phase peaks\nin the x = 0.10 sample is significantly weaker than in the x = 0.075 sample. These results highlight the\ncritical role of Sb Te in stabilizing the cubic GeSe phase. While increasing Sb Te content may\n2 3 2 3\nprevent the formation of the orthorhombic phase at high temperatures, it may also reduce the ZT value.\nThis information has been incorporated on Page 9 of the revised manuscript.\nComment 13: The authors should clarify the role of point defects and vacancies in enhancing\nthermoelectric performance.\nAutor reply: This can be explained in two aspects. First, orthorhombic GeSe has a large band gap and\na high formation energy for Ge vacancies, resulting in an ultra-low carrier concentration of 1.2 × 1016\ncm−3, which limits its potential to achieve higher ZT values. Increasing the carrier concentration is\ntherefore crucial for optimizing electrical conductivity and improving ZT. Aliovalent doping is a\ncommon strategy to enhance carrier concentration, such as introducing Cu, Na, or Ag into the Ge site\nin p-type GeSe. However, the low solubility of dopants in orthorhombic GeSe limits the effectiveness\nof this approach. Alternatively, creating Ge vacancies offers a viable pathway, as these vacancies act\nas acceptors that contribute holes to the matrix.\nDensity Functional Theory (DFT) calculations show that high-symmetry phases like\nrhombohedral and cubic GeSe have lower formation energies for Ge vacancies compared to the\northorhombic phase. Crystal Orbital Hamilton Population (COHP) analysis further indicates that Ge\nvacancies stabilize these high-symmetry phases. Alloying GeSe with Sb Te facilitates the formation\n2 3\nof Ge vacancies and promotes the cubic phase. The cubic phase further reduces the formation energy\nof Ge vacancies, effectively increasing the carrier concentration from 1.2 × 1016 cm−3 for orthorhombic\nGeSe to 3.0 × 1020 cm−3 for cubic (GeSe) (Sb Te ) . However, the larger atomic size of Pb compared\n0.9 2 3 0.1\nto Ge slightly increases the formation energy of Ge vacancies, leading to a minor reduction in carrier\nconcentration in Pb-containing samples. Nonetheless, the significantly increased carrier concentration\n(>1020 cm−3) contributes to higher electrical conductivity and an improved power factor.\nAdditionally, Sb Te and Pb co-alloying introduces substitution and vacancy point defects,\n2 3\ndisrupting lattice periodicity and inducing mass and strain field fluctuations, which enhance phonon\nscattering. This significantly reduces lattice thermal conductivity. The combined optimization of\nelectrical and thermal transport properties due to point defects and vacancies supports high\nthermoelectric performance.\nThis information has been included on Pages 18 and 24 of the revised manuscript.\nComment 14: A deeper discussion of the potential environmental impact of using Pb in the alloying\nprocess would be valuable.\nAutor reply: Pb, as a heavy metal, is widely recognized for its environmental hazards, making the\ndevelopment of eco-friendly, Pb-free materials a priority in advanced functional material research.\nHowever, the environmental impact of Pb depends significantly on its concentration in the material.\nIn this study, only a trace amount of Pb (5%) is incorporated, which plays a critical role in enhancing\nthe thermoelectric performance of cubic GeSe. Unlike conventional Pb-chalcogenides with much\nhigher Pb content, the minimal Pb concentration used here ensures a negligible impact on\nenvironmental pollution. Additionally, the Pb atoms are securely stabilized within the crystal structure\nduring the operation of the thermoelectric material, further minimizing the risk of environmental\nrelease. These factors indicate that the trace amount of Pb introduced in this study has a negligible\nenvironmental footprint, particularly compared to widely studied Pb-chalcogenides. Furthermore,\nwhile toxicity and high ZT values often represent a trade-off, thermoelectric performance remains a\nprimary focus. Nonetheless, we are committed to future research aimed at identifying non-toxic\nalternatives to replace Pb, ensuring high thermoelectric performance and maximum ZT without\ncompromising environmental safety.\nThis information has been incorporated on Page 27 of the revised manuscript.\nComment 15: Can the authors explain the computational techniques used for the electronic structure\ncalculations and their accuracy in predicting thermoelectric properties?\nAutor reply: The details of the DFT calculations are provided in the preceding response. Regarding\nthe accuracy of predicting thermoelectric (electrical) properties based on electronic structure\ncalculations, DFT results can be integrated into linearized Boltzmann transport equations (BTE) to\ndetermine electronic transport and, consequently, the fundamental electrical properties. Specifically,\nthe BTE is solved for small perturbations in the distribution function caused by thermal or electric\nfields. From this, expressions for the Seebeck coefficient, electrical conductivity, and electronic\nthermal conductivity can be derived, which are then combined to evaluate the ZT. By incorporating\nscattering mechanisms, such as impurity scattering and phonon-electron scattering, alongside the band\nstructure, the linearized BTE provides a comprehensive framework for predicting thermoelectric\nproperties. This approach is widely used in thermoelectric research and has proven effective in\nunderstanding and optimizing material performance.\nComment 16: How did the choice of computational models affect the electronic band structure\npredictions, and what limitations might exist?\nAutor reply: The choice of pseudopotentials in VASP plays a crucial role in the calculation of band\nstructures, as pseudopotentials simplify the treatment of atomic cores while preserving the essential\nphysics of valence electrons. Among the commonly used pseudopotentials in VASP, norm-conserving\npseudopotentials ensure that pseudo-wavefunctions match all-electron wavefunctions outside a cutoff\nradius. However, they require significant computational resources and time. Ultrasoft pseudopotentials\nallow greater flexibility near the core and reduce the number of plane waves needed for a given\naccuracy but may inaccurately estimate band gaps due to their inherent limitations. Projector\naugmented-wave (PAW) potentials combine the advantages of both norm-conserving and ultrasoft\npseudopotentials, maintaining accuracy for core states while optimizing computational efficiency. In\nthis study, DFT calculations were performed using the generalized gradient approximation (GGA) with\nthe fully relativistic Perdew-Burke-Ernzerhof (PBE) functional in the PAW formalism to handle\nexchange-correlation effects. While GGA-PBE pseudopotentials are known to overestimate lattice\nparameters and underestimate band gaps, we addressed this issue by conducting ionic optimization\nbased on GGA-PBEsol (a version of PBE revised for solids). This approach relaxes both lattice\nparameters and atomic positions, reducing discrepancies between experimental results and DFT\npredictions.\nComment 17: Could the authors provide more details on the pseudopotentials and exchange-\ncorrelation functional used for the electronic structure calculations?\nAutor reply: We utilized the all-electron Projector Augmented-Wave (PAW) method to calculate the\nelectronic structure efficiently and accurately. This method explicitly treats both core and valence\nelectrons, combining the advantages of all-electron and pseudopotential approaches. The generalized\ngradient approximation (GGA) with the fully relativistic Perdew-Burke-Ernzerhof (PBE) functional\nwas used to handle exchange-correlation effects. Compared to the local density approximation (LDA),\nGGA-PBE incorporates not only the local electron density but also its gradient, providing a more\naccurate description of exchange-correlation interactions. This is particularly beneficial for systems\nwith significant electron density variations, such as narrow-gap semiconductors like GeSe.\nComment 18: What are the potential sources of error in the electronic structure predictions, and how\ncan they be mitigated in future studies?\nAutor reply: The potential sources of error in electronic structure predictions and their mitigation\nstrategies are outlined as follows:\n1. Inherent limitations of pseudopotentials: Pseudopotentials are commonly used in DFT\ncalculations to reduce computational costs by simplifying the treatment of core electrons. However,\nthey approximate the full atomic potential, which may lead to inaccuracies in describing electron-core\ninteractions.\nPotential solutions: Develop advanced pseudopotentials that more accurately describe core\nelectrons and minimize transferability issues. Adopt hybrid functionals to improve accuracy or\ncomplement DFT calculations with quantum-chemistry methods.\n2. Supercell method: Simulating doped thermoelectric materials requires constructing supercells\nwhere host atoms are replaced by dopants. However, these supercells may not fully represent the actual\ndoping level, and multiple substitution sites can introduce variability.\nPotential solutions: With advances in computational power, larger supercells can better\napproximate real doping levels. Two approaches are commonly considered: One is the virtual crystal\napproximation which treats the different atoms in the supercell as a virtual single atom with averaged\nproperties (such as the average atomic number, mass, or other relevant properties) to represent the\nentire mixture of species. But this method is physically incorrect. The other is the generation of quasi-\nrandom structure based on evolutionary algorithms, which is still timely constant and also lacks\nphysical explanation. The development of AI might help this method to be more precise and efficient.\n3. Eigenvalues vs finite temperature: DFT calculations typically represent ground-state\neigenvalues, whereas thermoelectric materials operate at high temperatures. This discrepancy means\nthat effects such as thermal excitation, temperature-dependent phase transformations, and lattice\ndynamics are not accounted for.\nPotential solutions: Perform molecular dynamics simulations at finite temperatures under a\ngiven ensemble to predict temperature-dependent structures for DFT input. This approach is currently\nthe standard for incorporating temperature effects. In the future, the development of temperature-\ninclusive fundamental theories could revolutionize DFT calculations for thermoelectric properties.\nReviewer #2\nTitle: Metavalent Alloying and Vacancy Engineering Enable State-of-the-Art Cubic GeSe\nThermoelectrics. This study provides valuable insights into the effects of meta valent alloying on the\nthermoelectric of GeSe materials. The work well written and can be accepted after addressing given\ncomments.\nAutor reply: We sincerely appreciate the positive feedback and constructive suggestions, which have\nsignificantly contributed to improving the quality of our manuscript.\nComment 1: There is an error in the caption for Figure S3, where the composition is incorrectly written\nas x = 0.005 instead of x = 0.05. Correct this inconsistency throughout the manuscript to reflect the\nactual data.\nAutor reply: We sincerely apologize for the oversight. The inconsistency has been addressed, and the\nrevised caption for Figure S3 now accurately reflects the correct composition of x = 0.05 throughout.\nComment 2: In Figure 2c, the bond lengths in orthorhombic GeSe are showed to range from 2.52 to\n2.61 Å, whereas those in cubic GeSe are longer, extending from 2.84 to 2.89 Å. Based on the bond\nlength-bond strength relationship, shorter bonds generally indicate stronger interactions and higher\nstability, while longer bonds suggest weaker interactions. Consequently, the orthorhombic phase\nappears more stable than the cubic phase due to stronger bonding. Please clarify how the x = 0.10\nsample, with weaker bonding in the cubic phase, demonstrates higher thermal stability compared to\nthe x = 0.075 sample.\nAutor reply: While shorter bonds are often associated with stronger interactions and higher thermal\nstability, thermal stability depends on multiple factors beyond bond length, including the bonding\nmechanism and structural context. The bond length criterion is most relevant when comparing\ncompounds with similar bonding types and structures. For instance, many metavalently bonded\nmaterials, such as PbTe and SnTe, exhibit long bonds and relatively weak chemical interactions, yet\nmaintain a stable cubic phase across a wide temperature range. This stability arises from their unique\nbonding mechanism, where adjacent atoms share half an electron pair, forming a half-filled sigma bond\nwith a bond order of approximately 1/2.\nIn this study, increasing the Sb Te content transitions the bonding mechanism in GeSe from\n2 3\ncovalent to metavalent bonding, stabilizing the cubic phase. Rietveld refinement shows that the x =\n0.075 sample has shorter bond lengths than the x = 0.10 sample, yet both retain the cubic phase at room\ntemperature. Upon heating, both samples develop partial orthorhombic phases, as observed in in-situ\nXRD (Figures S2d and 2e). However, the orthorhombic phase emerges at 623 K for the x = 0.075\nsample, whereas it appears only at 673 K for the x = 0.10 sample. Furthermore, the intensity of the\northorhombic phase peaks is significantly weaker in the x = 0.10 sample, indicating that it retains the\ncubic structure more effectively at higher temperatures. These findings demonstrate that increasing the\nSb Te content enhances the stability of the cubic phase and reduces the formation of the orthorhombic\n2 3\nphase at elevated temperatures. We acknowledge that the term “thermal stability” may not have fully\ncaptured this behavior and have revised the manuscript to more precisely reflect these results. We are\ngrateful for the reviewer’s constructive feedback, which has been instrumental in improving the clarity\nand accuracy of our discussion.\nThe corresponding revisions have been made on Page 9 of the revised manuscript.\nComment 3: In Figure 2e, it is evident that the x = 0.10 sample reverts to the orthorhombic phase\n(peaks) at temperatures exceeding 673 K. Given this observation, the authors cannot use the word\n\"excellent thermal stability.\" Revise the statement.\nAutor reply: As noted, a discernible orthorhombic phase appears in the x = 0.10 sample at\ntemperatures exceeding 673 K. While increasing the Sb Te content could potentially suppress this\n2 3\nphase transition, it may also negatively impact the ZT value. Therefore, we have revised the text to\navoid describing the thermal stability as “excellent” and have updated the wording accordingly on\nPage 9 of the revised manuscript.\nComment 4: In Figure 1b, the evaluated ΔS values for x = 7.5% and 10% are shown as approximately\n0.2 and 0.25 R, respectively. However, on page 9, line 178, the authors state ΔS values of 0.98 and\n1.21 for these compositions. This discrepancy needs to be resolved. Please clarify which set of values\nis correct and provide a consistent explanation for the difference in reported data.\nAutor reply: We sincerely apologize for any confusion caused. After thoroughly reviewing both\nFigure 1b and the main text, we have confirmed that the data are consistent. The apparent discrepancy\narises from the presentation of two distinct sets of values in Figure R3 (Figure 1b in the revised\nmanuscript). The bar chart corresponds to the alloying content (plotted on the left vertical axis), while\nthe red line represents the configurational entropy (plotted on the right vertical axis). To ensure clarity,\nwe have added explanatory notes to the caption of Figure 1b in the revised manuscript, providing a\nclearer understanding for reviewers and readers.\nFigure R3 (Figure 1b in the revised manuscript). Alloying content (x) and configurational entropy\n(ΔS) required for obtaining cubic GeSe studied herein in comparison with literature data (the bar chart\ncorresponds to x and the red line represents ΔS).\nReviewer #3\nThe submitted manuscript, \"Metavalent Alloying and Vacancy Engineering Enable State-of-the-Art\nCubic GeSe Thermoelectrics.” offers a compelling exploration of superior thermoelectric transport in\ncubic p-type GeSe systems, which could be of interest to the materials research community. The\nauthors shed new light on the Sb Te metavalent alloying and vacancy engineering as a strategy for\n2 3\northorhombic to cubic phase transitions in GeSe with enhanced TE performance. The manuscript is\nwell-written and includes extensive analysis. While the manuscript shows promise, some aspects\nrequire revisions to enhance its impact. Before considering publication in Nature Communications, I\nrecommend the authors address the following comments.\nAutor reply: We sincerely appreciate the positive feedback and constructive suggestions, which have\nsignificantly improved the quality of our manuscript.\n1. Phase formation\nComment 1: The authors must report lattice parameters and bond lengths with proper error bars in the\ngraphs (Fig 2 (b), (c), S5 and S6).\nAutor reply: In response to your comments, we have updated the graphs (Figures 2b, 2c, S5, and S6)\nto include appropriate error bars (Figures R4–R6). The uncertainty in the X-ray diffraction (XRD)\nstructural refinement, performed using the Fullprof Suite, has been quantified as ±0.002 Å. This level\nof precision is also detailed in the “Phase and Microstructure Characterization” section of the\nexperimental details in the revised supporting information.\nFigure R4 (Figure 2b and c in the revised manuscript). Calculated (b) lattice parameters a, b, c and\n(c) bond distances.\nFigure R5 (Figure S5 in the revised Supporting Information). Calculated lattice parameters a, b, c of\n(Ge Pb Se) (Sb Te ) samples.\n1−y y 0.9 2 3 0.1\nFigure R6 (Figure S6 in the revised Supporting Information). Calculated bond distance of\n(Ge Pb Se) (Sb Te ) samples.\n1−y y 0.9 2 3 0.1\nComment 2: Authors discussed the change in site occupancy in alloyed samples, evaluated from\nRietveld analysis. Did authors observe such variations due to Pb incorporation as well? Please list\nessential refined parameters (R factors, site occupancy, Caglioti parameters, Biso etc.) with standard\ndeviations in supplementary information to support the arguments.\nAutor reply: The XRD data were refined using Rietveld analysis with the Fullprof software. The key\nrefined parameters, including R factors, site occupancy, Caglioti parameters, and Biso, are provided in\nTable R2 of the revised supporting information. Sb\n3\nalloying reduces the cation site occupancy,\nthereby increasing the concentration of cation vacancies. In contrast, the addition of Pb to the x = 0.10\nsample increases the cation site occupancy, reducing the vacancy concentration, as shown in Table R3.\nSpecifically, the cation site occupancy increases from 96.8% in the x = 0.10 sample to 99.1% in the x\n= 0.10, y = 0.07 sample.\nThe corresponding information has been revised and included on Page 11 of the revised\nmanuscript.\nTable R2 (Table S1 in the revised supporting information). The essential refined parameters of R\nfactors: R-pattern (R ), R-weighted pattern (R ), R-expected (R )18, Caglioti parameters (U, V, W),\np wp exp\nand Biso of (Ge Pb Se) (Sb Te ) samples.\n1−y y 1−x 2 3 x\nCaglioti parameters Biso\nSamples Phase R R R\np wp exp\nU V W Ge/Sb/Pb Se/Te\nx = 0 Orthorhombic 7.37% 9.38% 1.84% 0.21 −0.04 0.02 3.59 1.70\nOrthorhombic 0.01 −0.03 0.04 0.42 1.36\nx = 0.05 4.57% 6.10% 1.75%\nCubic 0.81 −0.12 0.01 0.13 1.21\nOrthorhombic 0.08 −0.07 0.07 0.16 0.48\nx = 0.075 5.78% 7.50% 1.89%\nCubic 1.23 −0.22 0.01 1.13 0.20\nx = 0.10 Cubic 4.80% 6.31% 1.86% 0.28 −0.05 0.01 1.49 0.12\nx = 0.125 Cubic 5.54% 7.06% 2.01% 0.12 −0.03 0.02 2.56 0.42\nx = 0.15 Cubic 5.00% 6.71% 1.91% 0.65 −0.11 0.02 0.05 0.07\nx = 0.10,\nCubic 5.23% 6.74% 1.78% 0.15 −0.03 0.01 1.12 0.02\ny = 0.03\nx = 0.10,\nCubic 5.62% 7.42% 1.73% 0.11 −0.01 0.01 1.62 0.23\ny = 0.05\nx = 0.10,\nCubic 5.31% 6.76% 1.72% 0.13 0.02 0.01 1.47 0.31\ny = 0.07\nTable R3 (Table S2 in the revised supporting information). The cation site occupancy of\n(Ge Pb Se) (Sb Te ) samples according to the XRD Rietveld refinement.\n1−y y 0.9 2 3 0.1\nCations Anions Cation site\nSamples\nGe Sb Pb Se Te occupancy\n(GeSe) (Sb Te ) 0.3882 0.0956 0.0000 0.3748 0.1242 0.9675\n0.9 2 3 0.1\n(Ge Pb Se) (Sb Te ) 0.3796 0.0936 0.0166 0.3852 0.1248 0.9796\n0.97 0.03 0.9 2 3 0.1\n(Ge Pb Se) (Sb Te ) 0.3840 0.0907 0.0202 0.3754 0.1242 0.9898\n0.95 0.05 0.9 2 3 0.1\n(Ge Pb Se) (Sb Te ) 0.3780 0.0908 0.0266 0.3792 0.1244 0.9908\n0.93 0.07 0.9 2 3 0.1\nComment 3: Do authors have any experimental data other than EDS, to support the existence of\nSb Te as a minor matrix? Furthermore, did the authors experimentally verify the ionic states of all the\n2 3\nelements, especially Sb and Te using techniques such as X-ray photoelectron spectroscopy? Please\njustify.\nAutor reply: To address the first question, we carefully considered the reviewer’s concern regarding\nthe potential presence of Sb Te secondary phases within the GeSe matrix. Energy-dispersive X-ray\n2 3\nspectroscopy (EDS) is widely recognized as a reliable method for analyzing phase compositions and\nelemental distributions. Accordingly, SEM-EDS mapping, TEM-EDS mapping, and their quantitative\nanalyses were performed, as shown in Figures S7–S12 and Figure 3a. The results demonstrate that\nall elements, including the matrix constituents (Ge and Se) and the dopants (Sb and Te), are uniformly\ndistributed within the measurement error. This confirms the successful formation of a single-phase\ncubic GeSe structure in the samples. Additionally, XRD analysis further supports these findings. As\nshown in Figures 2a and S2, no diffraction peaks corresponding to Sb Te secondary phases were\n2 3\ndetected in any of the samples. To further validate this conclusion, Rietveld refinement of the XRD\ndata was conducted, which also confirmed the presence of a pure cubic GeSe phase. These results\ncollectively indicate that Sb Te is fully incorporated into the GeSe matrix without forming secondary\n2 3\nphases.\nIn response to the reviewer’s second suggestion, X-ray photoelectron spectroscopy (XPS) was\nused to analyze the ionic states of all elements in the (Ge Pb Se) (Sb Te ) sample, as shown in\n0.95 0.05 0.9 2 3 0.1\nFigure R7. The analysis revealed that the valence states of Ge, Pb, Sb, Se, and Te are 2+, 2+, 3+, 2−,\nand 2−, respectively (Figure S13).\nFigure R7 (Figure S13 in the revised Supporting Information). XPS spectra of the scanning of (a) Ge,\n(b) Pb, (c) Sb, (d) Se, and (f) for the (Ge Pb Se) (Sb Te ) sample.\n0.95 0.05 0.9 2 3 0.1\n2. Thermoelectric transport\nComment 1: Temperature dependency of electrical conductivity (σ) found to be different for different\nx values (Fig. 5(b)), for instance, Pb doped samples (y = 0.03, 0.05, 0.07 with x = 0.1) show a strict\nincrease in σ with temperature as compared to other samples. Please give a brief regarding the\ntemperature dependency of σ.\nAutor reply: In pristine orthorhombic GeSe and low Sb Te -alloyed samples (0 ≤ x < 0.10),\n2 3\nthe\nelectrical conductivity increases with temperature, exhibiting typical semiconductor transport behavior.\nThis is primarily due to the large band gap and high Ge vacancy formation energy, which result in an\nultralow carrier concentration. The thermally activated behavior is attributed to carrier grain boundary\nscattering, as reported in the literature 15-17. Alloying GeSe with Sb Te induces a phase transition from\n2 3\nthe covalently bonded orthorhombic phase to the metavalently bonded cubic phase, significantly\nincreasing the carrier concentration. This higher carrier concentration reduces the effect of grain\nboundary scattering, which explains why, in high Sb Te -alloyed samples (x ≥ 0.10), the electrical\n2 3\nconductivity decreases with temperature at low temperatures, reflecting metal-like transport behavior.\nIntroducing Pb into the system slightly reduces the carrier concentration, shifting the electrical\ntransport behavior back toward semiconductor-like conduction. This illustrates the intricate interplay\nbetween chemical bonding, phase structure, carrier concentration, and electrical transport mechanisms\nin determining electronic transport properties. To further enhance carrier mobility and minimize the\neffects of grain boundary scattering, optimizing grain size—particularly through single-crystal\ngrowth—is recommended for future studies. This approach is expected to improve electrical transport\nproperties and contribute to higher ZT values.\nA more detailed discussion of these electrical transport behaviors is provided on Page 19 of the\nrevised manuscript.\nComment 2: Interestingly, the temperature dependency of weighted mobility of alloyed/doped\nsamples are different than the pristine sample. Please describe this in terms of carrier scattering with\ntemperature. Please mention the equation used in estimating weighted mobility in supplementary\ndetails.\nAutor reply: As mentioned in response to Comment 1, we addressed the potential carrier scattering\nmechanisms influencing these observations.\nIn pristine orthorhombic GeSe, the weighted mobility increases with temperature, reflecting\ntypical semiconductor transport behavior primarily driven by carrier grain boundary scattering. In\ncontrast, Sb Te -alloyed samples show a decrease in weighted mobility with temperature due to the\n2 3\nsignificantly enhanced carrier concentration, indicative of metal-like transport behavior. For instance,\nthe x = 0.10 sample exhibits a mobility-temperature relationship of μ ~ T−0.45, suggesting that point\nW\ndefect scattering is the dominant mechanism. With the addition of Pb, the carrier concentration\ndecreases slightly, shifting the electrical transport properties back toward semiconductor-like behavior\nat room temperature, where grain boundary scattering again becomes the primary mechanism.\nAdditionally, the formula used to calculate weighted mobility, as referenced in 19 is detailed in\nthe experimental section in the revised supporting information.\n  S  3 S \n3\n exp 2 \n 331 cm2 mΩ cm T   2    k B /e   2 k B /e  (R2)\n    \nW Vs   300   S    S \n \n1exp5 1 1exp5 1\n    k B /e    k B /e   \nWe have added the “weighted mobility and quality factor calculation model” in the experimental\nsection of the revised supporting information.\nComment 3: In general, the heat capacity C values estimated by Dulong - Petit’s law will\np\nunderestimate the thermal conductivity values at high temperatures. I suggest the authors refer to an\nempirical equation for C from the article- Wang, H., Gibbs, Z. M., Takagiwa, Y., & Snyder, G. J.\np\n(2014). Tuning bands of PbSe for better thermoelectric efficiency. Energy & Environmental Science,\n7(2), 804-811. Please compare the values calculated from both the methods and use only the\nappropriate one in this case. Justify it in reply letter.\nAutor reply: We appreciate the reviewer’s perspective. Generally, C increases with temperature,\np\nwhich can lead to an underestimation of thermal conductivity at high temperatures. While measuring\nthe temperature-dependent C curve would provide more accurate thermal conductivity values, such\np\nmeasurements are both complex and challenging in terms of precision. To the best of our knowledge,\nDulong-Petit’s law is widely employed to estimate the C of GeSe in most reported studies. For\np\nconsistency and comparability with existing literature, using C estimated by Dulong-Petit’s law is\np\nmore appropriate (Table R4). We thank the reviewer for their understanding.\nTo address the reviewer’s concern, we also calculated C using the empirical equation C /k per\np p B\natom = 3.07 + 4.7 × 10–4(T/K - 300) 20, with the results presented in Table R5. Additionally, we\ncompared the total thermal conductivity of the (Ge Pb Se) (Sb Te ) sample obtained using\n0.05 0.05 0.9 2 3 0.1\ndifferent methods for C calculation, along with comparisons of lattice thermal conductivity and ZT\np\nvalues (Figure R8). The results show that C values calculated using the empirical equation are slightly\np\nhigher than those estimated by Dulong-Petit’s law, especially at elevated temperatures. However, the\nresulting differences in total thermal conductivity and lattice thermal conductivity are minimal,\namounting to only 0.06 W m−1 K−1. This leads to a minor decrease in the ZT value by 0.12 compared\nto the value reported in our manuscript, which is well within the measurement uncertainties.\nTherefore, for consistency and to ensure better comparability with previously published studies,\nwe believe that using C estimated by Dulong-Petit’s law remains the most appropriate approach. We\np\nappreciate the reviewer’s understanding.\nFigure R8. (a) Comparison of the total thermal conductivity of (Ge Pb Se) (Sb Te ) sample\n0.05 0.05 0.9 2 3 0.1\nobtained using different methods for specific heat calculation, along with the corresponding\ncomparisons of (b) lattice thermal conductivity and (c) ZT values.\nTable R4. The heat capacity C values estimated by Dulong - Petit’s law. The unit of C is J g−1K−1.\np p\n(Ge Pb Se) (Sb Te )\n1-y y 1-x 2 3 x\nx=0.10 x=0.10 x=0.10\nx=0 x=0.05 x=0.075 x=0.10 x=0.125 x=0.15\ny=0.03 y=0.05 y=0.07\n0.329 0.305 0.296 0.288 0.280 0.274 0.283 0.279 0.276\nTable R5. The heat capacity C values estimated by empirical equation. The unit of C is J g−1K−1.\np p\n(Ge Pb Se) (Sb Te )\n1-y y 1-x 2 3 x\nx=0.10 x=0.10 x=0.10\nT (K) x=0 x=0.05 x=0.075 x=0.10 x=0.125 x=0.15\ny=0.03 y=0.05 y=0.07\n300 0.337 0.313 0.303 0.295 0.287 0.281 0.290 0.286 0.283\n323 0.338 0.314 0.304 0.296 0.288 0.282 0.291 0.287 0.284\n373 0.341 0.317 0.307 0.298 0.291 0.284 0.293 0.289 0.286\n423 0.343 0.319 0.309 0.301 0.293 0.286 0.295 0.292 0.288\n473 0.346 0.321 0.311 0.303 0.295 0.288 0.297 0.294 0.290\n523 0.348 0.324 0.314 0.305 0.297 0.290 0.300 0.296 0.293\n573 0.351 0.326 0.316 0.307 0.299 0.292 0.302 0.298 0.295\n623 0.353 0.329 0.318 0.310 0.302 0.295 0.304 0.300 0.297\n673 0.356 0.331 0.321 0.312 0.304 0.297 0.306 0.303 0.299\n723 0.359 0.333 0.323 0.314 0.306 0.299 0.308 0.305 0.301\nComment 4: In Figure 6f, the Grüneisen parameter shows an increase with x value up to an optimum,\nbeyond which it decreases. It would be helpful to the readers if the authors can provide a brief regarding\nthis.\nAutor reply: We also observed an irregular variation in the Grüneisen parameter with increasing\nSb Te content (Figure 6f). This trend is reflected in the measured sound velocity (Figure 6e), which\n2 3\ndecreases initially with increasing Sb Te content up to x = 0.05, before gradually increasing thereafter.\n2 3\nPrevious studies have shown that metavalently bonded materials exhibit longer-range interactions\ncompared to covalently bonded counterparts, leading to optical phonon softening and strong\nanharmonic scattering 14. This suggests that materials with metavalent bonding inherently possess\nhigher intrinsic anharmonicity. In this study, alloying GeSe with Sb Te induces a transition in the\n2 3\nchemical bonding mechanism from covalent to metavalent bonding. The initial increase in the\nGrüneisen parameter is attributed to this transition, which enhances lattice anharmonicity.\nHowever, the subsequent decrease in the Grüneisen parameter at higher Sb Te content indicates\n2 3\nthe involvement of more complex mechanisms, potentially related to the impact of Sb Te alloying on\n2 3\nthe phonon dispersion relationship. To better understand phonon transport behavior in this system, we\ninitially planned to calculate the phonon dispersion curves for (GeSe) (Sb Te ) . Unfortunately, we\n1−x 2 3 x\nencountered significant computational challenges due to the high resource demands, which exceeded\nour current capabilities. We kindly request the reviewers’ understanding regarding this limitation.\nIn summary, the variation in the Grüneisen parameter likely arises from a combination of factors,\nincluding the transition in bonding mechanism and changes in phonon dispersion induced by Sb Te\n2 3\nalloying. Further investigations, particularly focusing on phonon dispersion, would provide deeper\ninsights into these phenomena. Additional details have been included on Page 24 of the revised\nmanuscript.\nComment 5: Even though the authors plotted the variation of μ versus 1/κ , it would be better to\nW L\nmention the values of the thermoelectric quality factor as well.\nAutor reply: The quality factor (β ) characterizes the balance between electron and phonon transport,\nB\nindependent of the optimization of charge carrier concentration. The β is determined by the ratio of\nB\nμ to κ , as described by the following equation 19:\nW L\n3\nk  2 8e2m k T 2 \n  B e B W (R3)\n \nB  e  3h3 \nL\nThe calculated β factor as a function of temperature is shown in the Supporting Information,\nB\nFigure R9 (Figure S16). At room temperature, the β increase from 3.68 × 10−4 for pristine GeSe to\nB\n6.63 × 10−2 for (GeSe) (Sb Te ) . Notably, β shows a more pronounced increase at elevated\n0.9 2 3 0.1 B\ntemperatures. For example, at 673 K, the β of pristine GeSe is 0.28, which rises to 0.50 for the x =\nB\n0.10 sample and further increases to 0.65 with 5 at% Pb doping. These findings highlight the\nsynergistic optimization of electrical and thermal transport properties in cubic GeSe.\nFigure R9 (Figure S16 in the revised Supporting Information). Quality factor β\nB\nas a function of\ntemperature of all samples.\nComment 6: Please mention the uncertainty in the electrical measurements. Propagate the error in\nestimating ZT values.\nAutor reply: To ensure rigor, the electrical and thermal transport properties in this study were\nmeasured using commercial equipment. Electrical resistivity and the Seebeck coefficient were\nmeasured on a ZEM-3 (Ulvac-Riko®, Japan) under a protective helium atmosphere, with measurement\nerrors of 10% and 7%, respectively. The thermal diffusivity was determined using the laser flash\nmethod (LFA 467, Netzsch®, Germany), with a measurement error of 3%. The Hall coefficient (R )\nH\nwas measured at room temperature using the physical properties measurement system (PPMS-9,\nQuantum Design®, USA). Carrier concentration (n ) and mobility (μ ) were calculated using the\nH H\nequations n = 1/eR and μ = σR , respectively, with an uncertainty of ±5% for R . As a result, the\nH H H H H\ntotal uncertainty in the calculated ZT values is approximately 20%. These measurement errors have\nbeen detailed in the Supporting Information under the “Thermoelectric Performance Measurement”\nsection. For clarity and visual aesthetics, error bars were not included in the figures. We appreciate the\nreviewer’s understanding of this approach.\nReferences\n1. Duan B, et al. The Role of Cation Vacancies in GeSe: Stabilizing High-Symmetric Phase Structure and\nEnhancing Thermoelectric Performance. Adv Energy Sustain Res 3, 2200124 (2022).\n2. Kresse G, Hafner J. Ab initio molecular-dynamics simulation of the liquid-metal-amorphous-\nsemiconductor transition in germanium. Phys Rev B 49, 14251-14269 (1994).\n3. Kresse G, Hafner J. Ab initio molecular dynamics for liquid metals. Phys Rev B 47, 558-561 (1993).\n4. Kresse G, Furthmüller J. Efficiency of ab-initio total energy calculations for metals and semiconductors\nusing a plane-wave basis set. Comp Mater Sci 6, 15-50 (1996).\n5. Kresse G, Hafner J. Norm-conserving and ultrasoft pseudopotentials for first-row and transition\nelements. J Phys-Condens Mat 6, 8245-8257 (1994).\n6. Kresse G, Furthmüller J. Efficient iterative schemes for ab initio total-energy calculations using a plane-\nwave basis set. Phys Rev B 54, 11169-11186 (1996).\n7. Kresse G, Joubert D. From ultrasoft pseudopotentials to the projector augmented-wave method. Phys\nRev B 59, 1758-1775 (1999).\n8. Perdew JP, Burke K, Ernzerhof M. Generalized Gradient Approximation Made Simple. Phys Rev Lett\n77, 3865-3868 (1996).\n9. Setyawan W, Curtarolo S. High-throughput electronic band structure calculations: Challenges and tools.\nComp Mater Sci 49, 299-312 (2010).\n10. Tran F, Blaha P. Accurate Band Gaps of Semiconductors and Insulators with a Semilocal Exchange-\nCorrelation Potential. Phys Rev Lett 102, 226401 (2009).\n11. Dronskowski R, Blöchl PE. Crystal orbital Hamilton populations (COHP): energy-resolved\nvisualization of chemical bonding in solids based on density-functional calculations. J Phys Chem 97,\n8617-8624 (1993).\n12. Togo A, Chaput L, Tadano T, Tanaka I. Implementation strategies in phonopy and phono3py. J Phys-\nCondens Mat 35, 353001 (2023).\n13. Toberer ES, Zevalkink A, Snyder GJ. Phonon engineering through crystal chemistry. J Mater Chem 21,\n15843-15852 (2011).\n14. Lee S, Esfarjani K, Luo T, Zhou J, Tian Z, Chen G. Resonant bonding leads to low lattice thermal\nconductivity. Nat Commun 5, 3525 (2014).\n15. Zhang X, et al. Thermoelectric properties of GeSe. J Materiomics 2, 331-337 (2016).\n16. Martin J, Wang L, Chen L, Nolas GS. Enhanced Seebeck coefficient through energy-barrier scattering\nin PbTe nanocomposites. Phys Rev B 79, 115311 (2009).\n17. Hu C, Xia K, Fu C, Zhao X, Zhu T. Carrier grain boundary scattering in thermoelectric materials. Energy\nEnviron Sci 15, 1406-1422 (2022).\n18. Toby BH. R factors in Rietveld analysis: How good is good enough? Powder Diffr 21, 67-70 (2006).\n19. Snyder GJ, Snyder AH, Wood M, Gurunathan R, Snyder BH, Niu C. Weighted Mobility. Adv Mater 32,\n2001537 (2020).\n20. Wang H, Gibbs ZM, Takagiwa Y, Snyder GJ. Tuning bands of PbSe for better thermoelectric efficiency.\nEnergy Environ Sci 7, 804-811 (2014)."}, {"round": 2, "reviewer_comments": {"reviewer_1": "After thoroughly reviewing and evaluating the revised manuscript, I recommend the acceptance of the research paper titled \"Metavalent Alloying and Vacancy Engineering Enable State-of-the-Art Cubic GeSe Thermoelectrics\" for publication.", "reviewer_2": "I have reviewed the responses to the comments and the revised manuscript. The author has addressed all the comments thoroughly and implemented the changes in the revised version. I don't have any further comment on the given paper.", "reviewer_3": "The manuscript quality has increased considerably after the revisions. The authors have answered two sets of comments (1. Phase formation and 2. Thermoeletric transport) out of three satisfactorily. But the answer to the 3rd set of comment is not included in the revised manuscript. Please find and address the comment before accepting for publication.\nThe comment is as follows.\n3. Mechanical Properties\n• Can substantial improvement in mechanical properties of alloyed/doped samples are also connected to the improved microstructure? Please provide a discussion on this."}, "author_response": "Autor reply: Thank you for your thoughtful feedback and for acknowledging the improvements in our manuscript. We sincerely appreciate your insightful comments and apologize for the oversight regarding the third set of comments. We will promptly address the remaining issues and ensure that the revised manuscript fully incorporates your suggestions.\n3. Mechanical Properties\nComment 1: Can substantial improvement in mechanical properties of alloyed/doped samples are also connected to the improved microstructure? Please provide a discussion on this.\nAutor reply: Thank you for this insightful question. Mechanical properties are fundamentally governed by both chemical bonding and microstructure. Pristine orthorhombic GeSe is characterized by covalent bonding, where adjacent atoms are held together by a shared electron pair, resulting in a bond order of one. In contrast, rhombohedral and cubic GeSe exhibit metavalent bonding (MVB), in which adjacent atoms share only half an electron pair (one electron) to form a σ-bond. Consequently, the bond order in MVB is 0.5, indicating weaker chemical bonds compared to covalent bonding. Additionally, covalent GeSe has a very low concentration of Germanium (Ge) vacancies, whereas metavalent GeSe contains a significantly higher concentration of Ge vacancies. These factors contribute to the theoretically lower hardness and strength of metavalent GeSe compared to covalent GeSe. However, the introduction of Sb Te and Pb co-alloying in cubic GeSe induces substantial solid-solution strengthening, leading to Vickers hardness (HV) and compressive strength (σ ) that surpass those of the original orthorhombic GeSe.\nIt is important to note that while solid-solution strengthening enhances hardness and strength, it can sometimes reduce toughness. However, the improved fracture toughness (K ) observed in Sb Te and Pb co-alloyed cubic GeSe is primarily a result of its unique chemical bonding rather than microstructural changes. MVB represents an intermediate state between covalent and metallic bonding, balancing localized and delocalized electronic interactions. This partial electron delocalization in metavalent materials enables greater atomic displacements, enhancing plasticity. Under external stress, metavalent materials effectively absorb and redistribute strain, preventing brittle fractures and exhibiting superior toughness. This explains why Sb Te and Pb co-alloyed cubic GeSe demonstrate significantly enhanced fracture toughness (K ) compared to pristine orthorhombic GeSe.\nThe revised discussions can be found on Page 20 of the revised manuscript."}, {"round": 3, "reviewer_comments": {"reviewer_3": "Authors have addressed all the comments satisfactorily. I recommend acceptance of the manuscript for publication in Nature Communications."}, "author_response": ""}]

Decision letter