First-principles density functional theory (DFT) calculations were used for estimating both the thermoelectric quality factor β and the dopability of each compound. In all cases, these simulations were performed using the Vienna Ab-initio Simulation Package (VASP, Kresse and Furthmüller (1996)), within the projector augmented wave (PAW) formalism (Blöchl (1994)). We used pseudopotentials with core/valence electrons as follows: Sn: [Kr]4d¹⁰:5s²5p², S: [Ne]:3s²3p², Se: [Ar]4d¹⁰:4s²4p⁴, Te: [Kr]4d¹⁰:5s²5p⁴. As described below, both the Perdew-Burke-Eznerhof (DFT-PBE, Perdew et al. (1996)) and the Heyd-Scuseria-Eznerhof (DFT-HSE, Krukau et al. (2006)) description of the exchange correlation potential were used. A plane wave energy cutoff of 400 eV and, for bulk unit cells, a Γ-centred k-point mesh of 5 × 5 × 5 were used to perform the relaxation. The total energy and force convergence criteria used for all geometry optimizations were 10^–5 eV and 0.01 eV/Å, respectively.

The thermoelectric quality factor was obtained using DFT-PBE, with spin-orbit coupling included for PbTe and SnTe. Here, we used DFT-PBE rather than DFT-HSE since dense k-meshes are required to converge computed parameters. Relaxed structures obtained as described above were used for subsequent electronic band structure and bulk modulus calculations.

The semi-empirical model to predict each material’s intrinsic mobility (Yan et al., 2015; Miller et al., 2017) requires density of states effective mass m DOS * and band degeneracy N _b . From these, the band effective mass m b * is obtained by m DOS * = N b 2 / 3 m b * , under the assumption of symmetry equivalent and isotropic bands. The DOS effective masses and band degeneracies themselves were obtained using dense k-mesh calculations carried out on unit cells. For SnTe and PbTe, a kpoint mesh of 20 × 20 × 20 was used for the 2-atom primitive cell. Since the VBM and CBM of both SnS and SnSe occur between high-symmetry corners of the Brillouin zone, a denser k-mesh of 28 × 28 × 28 was needed to converge m DOS * and N _b for the 8 atom unit cell (see Supplementary Figure S1). The DOS effective mass is extracted from a 100 meV energy window from the DOS at the VBM and CBM under the parabolic band approximation.

The semi-empirical model used to predict lattice thermal conductivity (Yan et al., 2015; Miller et al., 2017) requires the bulk modulus as an input. This was obtained by fitting the Birch–Murnaghan (Birch, 1952) equation of state to a set of total energies computed at 5 different volumes near the minimum volume of the relaxed structure.

We note that both models predict intrinsic properties assuming pure material without any defects, while quantities like mobility and lattice thermal conductivity can be affected if defects or impurities are present in the system at high concentration.

For all aspects of dopability assessment (except m DOS * as described above) including phase stability, band gaps, and defect formation energies, we used DFT-HSE. Optimized lattice constants and internal coordinates were obtained form HSE06 (Krukau et al., 2006) with an exchange mixing of α = 0.2. A comparison of DFT-predicted and experimental lattice constants and band gaps are shown in Supplementary Table S1.

To begin, we applied the quality factor approach to SnS, SnSe, and SnTe. The quality factor β was evaluated for each material in both the p-type and n-type regime, in order to determine the attainable zT and its dependence on carrier concentration. The approach is summarized here for completeness; additional details on the quality factor framework itself can be found in Chasmar and Stratton (1959) and on the semi-empirical models to obtain needed parameters in Yan et al. (2015); Chasmar and Stratton (1959).

Within the relaxation time approximation, the solutions to the Boltzmann transport equation offer an expression for zT based on material descriptor β, as z T = α 2 σ T κ L + κ e = u β 1 + v β (6) Here, α, σ, κ _L , and κ _e are the Seebeck coefficient, electrical conductivity, lattice thermal conductivity, and electrical thermal conductivity, respectively. The functions u and v are entirely dependent on the reduced Fermi level η and charge carrier scattering mechanism ν, while β is an intrinsic material parameter (Chasmar and Stratton (1959)). From Figure 1C, it is evident that both high β and optimized η should be concurrently achieved to maximize zT. The precise shape of the curve and the optimal carrier concentration depends on the dominant scattering mechanism (e.g. polar optical phonon scattering, impurity scattering, etc).

Non-dimensional parameter β involves quantities such as carrier effective mass and “intrinsic” mobility, band degeneracy, and lattice thermal conductivity. Rigorously β is given by ( k B / e ) 2 ( σ 0 T / κ L ) where k _B is the Boltzmann constant, e the electron charge, T the temperature, and σ ₀ the “intrinsic” electrical conductivity (Chasmar and Stratton, 1959; Yan et al., 2015). The electronic conductivity σ = σ ₀ exp(η), where η = −Δ/k _b T and Δ is the positive distance between the Fermi level and the band edge. The quantity σ ₀ is determined by the intrinsic mobility μ ₀ and density of states (DOS) effective mass m DOS * . Putting everything together, the expression for β(T) is β T = 5.76 × 1 0 − 6 μ 0 m DOS ∗ / m e 3 / 2 κ L T 5 / 2 (7) The intrinsic mobility μ ₀ and lattice thermal conductivity κ _L themselves exhibit a temperature dependence, that itself is sensitive to the scattering regime.

To facilitate high-throughput evaluations of β, semi-empirical models that predict μ ₀ and κ _L at T = 300 K are available (Yan et al., 2015; Miller et al., 2017). The model for μ ₀ accounts for acoustic and optical electron-phonon scattering, and depends on the bulk modulus (B) and carrier band effective mass ( m b * ) . Typically, this model predicts mobilities that lie within a factor of ten of experiment across a diverse set of materials with mobilities spanning several orders of magnitude (Yan et al., 2015). The model for κ _L includes contributions from acoustic and optical phonon modes and depends on the material density and speed of sound. This model has been extended to better capture anharmonicity, resulting in accurate predictions of κ _L at 300 K for a range of materials (Miller et al., 2017). We also note that the semi empirical models used here do not explicitly account for anisotropy, but estimate average values, which may contribute to uncertainties in computed properties.

Table 1 reports the value of β we obtained for SnS, SnSe, and SnTe. The classic TE material PbTe is also shown for reference. For each material, the value of β is shown for the principal crystal structure(s) in which it is observed. SnS and SnSe are stable in layered orthorhombic Pnma crystal structure, a distorted version of rock salt (Figure 2A), at room temperature. Although both SnS and SnSe undergo a phase transition around 950 K and 750–800 K from Pnma to Cmcm crystal structure (Figure 2B), in Table 1 we considered only the Cmcm phase for SnSe as the transition temperature for SnS is quite high. In the Cmcm structure, the crystallographic orientation of Sn and Se atoms is altered so that Sn cations are bonded to five equivalent anions and form a mixture of corner and edge-sharing square pyramids. On the other hand, SnTe is stable in the cubic rock salt structure, where each Sn is bonded with six Te atoms forming edge and corner-sharing symmetric octahedral bonds as shown in Figure 2C.

The predicted β and associated parameters that we obtained are summarized in Table 1. The values reported here are updated from those in Yan et al. (2015) due to the inclusion of spin-orbit coupling and the improved thermal conductivity model (Miller et al., 2017). For each material we show the quality factor normalized by the corresponding value in p-type PbTe, for both p-type (β _p/β _PbTe,p) and n-type (β _n/β _PbTe,p) carriers.

Focusing on p-type performance, from Table 1, both SnSe (Cmcm) and SnTe have higher β _p than PbTe. The large β _p for SnSe is expected; experimentally, p-type SnSe exhibits the highest reported zT amongst all single crystal materials (Zhao et al., 2014). On the other hand, SnTe does show reasonable zT experimentally, but larger values are only achieved when alloyed (Al Rahal Al Orabi et al., 2016; Doi et al., 2019; Ma et al., 2019; Tang et al., 2018; Tan et al., 2015). The Cmcm phase of SnSe shows higher β _p than the Pnma phase, consistent with the experimentally reported higher zTs (Zhao et al., 2014; Qin et al., 2018; Zhao et al., 2015; Chere et al., 2016). The larger β _p for Cmcm arises from the higher hole mobility and comparable lattice thermal conductivity to the Pnma phase; the higher hole mobility in turn arises from a lower band effective mass. Similarly, we obtained good β _p for SnTe due to very low valence band effective mass resulting in high mobility. Although we predicted β _p to be low in SnS, recent experiments have shown that a zT of 1.6 is achievable in p-type SnS when alloyed with Se (He et al., 2019).

For n-type performance indicator β _n, we predict the highest performance again for SnSe, but now in the low symmetry Pnma phase due to high band degeneracy in the conduction band and high mobility which also matches previous computational analysis on SnSe (Kutorasinski et al., 2015). Experiments also support this prediction, as almost all n-type reports of high zT in SnSe are in the Pnma phase (Zhang et al., 2015; Duong et al., 2016; Chang et al., 2016; Shang et al., 2019; Cha et al., 2019; Cai J et al., 2020; Gu et al., 2021). Meanwhile the n-type ratio β _n/β _PbTe,p is 0.73 for SnTe, and is only 0.27 for SnS. There are no reports of n-type SnTe unless heavily alloyed with PbTe (Pang et al., 2021), and very few reports of n-type SnS (Zhi-Cheng et al., 2019; Hu et al., 2019).

We next used the computed values of β to predict the carrier concentration dependence of zT for each material. Figures 3, 4 show the carrier concentration dependence for p-type and n-type materials respectively, which are compared to previous experimentally reported values we obtained from literature. Since the predicted values of β in Table 1 are valid for 300 K, to make each curve we used Eq. 7 to estimate the corresponding value for β at 900 K for SnS (Pnma) and SnTe (rock salt) and 800 K for SnSe (Pnma and Cmcm) phases, which are more typical of the temperatures at which zT has been reported. The predicted zT at different temperatures are also shown in Supplementary Figures S2, S3.

To extrapolate β to these higher temperatures, we assumed that the ratio of the mobility μ ₀ to the lattice thermal conductivity κ _L remains constant. In Figures 3, 4, dots indicate reported zT for doped materials, while asterisks are used to denote prior measurements for alloyed systems. Doping typically refers to the controlled addition of a small quantity of impurities to adjust electrical properties, whereas alloying is the formation of a solid solution; here the ‘doped’ and ‘alloyed’ labels are assigned based on the original authors’ classification.

Predicted zT in the p-type Sn chalcogenides is shown in Figure 3. The red line shows the dependence of zT on carrier concentration for the case of polar optical phonon limited carrier mobility, while the blue line corresponds to impurity scattering limited mobility. Both of these regimes are shown since, although polar optical phonon scattering has been predicted to be the dominant bulk scattering mechanism in most binary chalcogenides (Ganose et al., 2021), impurity scattering can also dominate when doping levels are substantial. Each line also contains a shaded region around it, corresponding to a factor of ten uncertainty in the semi-empirical model for carrier mobility.

For both Pnma and Cmcm SnSe, the computational predictions are in good agreement with experiment (Figures 3B,C). This comparison also shows that the predicted high zT has been achieved experimentally (Zhao et al., 2014; Qin et al., 2018; Zhao et al., 2015; Chere et al., 2016) for p-type SnSe. On the other hand for SnS and SnTe, there is more discrepancy between the model predictions and experiment. These discrepancies likely arise 1) from uncertainties in the model and computed parameters, and 2) because many of the experimentally reported values correspond to alloys (especially for SnTe). In spite of these differences, from Figure 3A it appears that higher zT may be possible for SnS if the hole carrier concentration could be increased beyond typical values (Zhou et al., 2017; Tan et al., 2014; He et al., 2019; Cai B et al., 2020) of 10¹⁸–10¹⁹ cm⁻³. For p-type SnTe Figure 3D suggests the opposite: that even higher zT could be achieved if hole concentrations could be reduced below typical values (Al Rahal Al Orabi et al., 2016; Doi et al., 2019; Ma et al., 2019; Tang et al., 2018; Tan et al., 2015) of ∼ 1 0 20 cm⁻³).

Similarly, predicted zT for n-type Sn chalcogenides is shown in Figure 4. Compared to p-type counterparts, there are fewer experimental measurements of n-type SnS and SnSe (Pnma), and no n-type performance has been reported for either SnSe (Cmcm) or SnTe (unless heavily alloyed with PbTe Pang et al. (2021)). For SnS (Figure 4A), only two experimental measurements (Zhi-Cheng et al., 2019; Hu et al., 2019) are available, showing zT = 0.12 and zT = 0.17 (the dots in Figure 4A overlap). Here, the predictions indicate that higher zT could be achieved by increasing the electron carrier concentration. For SnSe (Pnma), reasonable n-type zTs have been reported. Within the uncertainty of the computational predictions, slight improvements may be possible with further carrier concentration tuning (Figure 4B). The predicted n-type performance of the Cmcm phase of SnSe is quite low compared to the Pnma phase (Figure 4C). For SnTe, although the inability to achieve n-type material is well-known it is interesting to note that the model predicts that for electron carrier concentrations on the order of ∼ 1 0 19 cm⁻³, zT ≈ 1.5 could be achievable (Figure 4D).

We note that the required carrier concentration to achieve maximum zT depends on the β of each material, which in turn arises from the material’s electronic structure. Starting with p-type compounds, SnS has a comparatively low β, therefore to reach its maximum zT high hole carrier concentration on the order of 10²⁰–10²¹ cm⁻³ is required. For p-type SnSe in the Pnma phase, β is higher than p-type SnS, so the required hole carrier concentration is on the order of ∼ 1 0 20 cm⁻³. SnSe in the Cmcm phase and SnTe have high β and so, the maximum TE potential can be achieved around lower carrier concentrations of 10¹⁹ cm⁻³ and 10¹⁸ cm⁻³ respectively for SnSe (Cmcm) and SnTe. Our prediction for p-type SnSe matches previous computational analysis that good p-type performance might be achievable when the carrier concentration is on the order of 10¹⁹–10²⁰ cm⁻³ (Shi and Kioupakis, 2015). In the n-type cases, we again observe similar behavior. For SnS and SnSe (Cmcm) phase β is low, so the required carrier concentration is between 5 × 10¹⁹–10²⁰ cm⁻³. However, for SnSe (Pnma) and SnTe we find that a carrier concentration of 10¹⁹ cm⁻³ would be sufficient for achieving the maximum TE potential. Therefore, our analysis reflects that if a material has high β then the required carrier concentration e.g. doping level would be lower than a material with low β.

However, as shown by the shaded curves in Figures 3, 4, there are uncertainties in the semi-empirical models used to estimate mobilities and lattice thermal conductivity, and therefore quality factor β. The evaluated mobilities introduce a particularly large uncertainty, as they may vary by an order of magnitude from experimentally measured ones. For instance, the predicted hole mobilities for SnS, SnSe, and SnTe are around 7, 15, and 1,000 cm²/(Vs) while the experimental hole mobilities are 80, 170, and 650 cm²/(Vs) at 300 K (Al Rahal Al Orabi et al., 2016; Ganose et al., 2021). The discrepancies arise from the empirical nature of the scattering models and the neglect of the role of defects themselves when estimating mobility. One effective approach to overcome the uncertainty of the predicted mobility would be to obtain the mobility from readily available experimental data instead, using the concept of the “weighted mobility” μ _w (Snyder et al., 2020). The weighted mobility can be obtained from measurements of the Seebeck coefficient and electrical resistivity, and is related to the drift mobility by μ w = μ ( m DOS * m e ) 3 / 2 . These experimentally-obtained values can instead be used to estimate the quality factor at specific temperatures of interest, so that Eq. 7 becomes β T = 5.76 × 1 0 − 6 μ w κ L T 5 / 2 . (8)

From here, target carrier concentrations can again be identified for maximum thermoelectric performance.

In summary, by optimizing carrier concentration improvements may be possible in p-type SnS and SnTe, and n-type SnS, SnSe (Pnma), and SnTe. Meanwhile the p-type performance observed in SnSe appears to be close to its peak predicted value. However, realizing improved performance relies on the ability to tune carrier concentrations, which in turn depends on native defects and dopability. To assess the possibility of achieving the target carrier concentrations, in the next section we analyze the dopability of SnS, SnTe, and SnSe (Pnma).

We next analyzed the defect chemistry and dopability of all three compounds to determine if the optimal carrier concentrations we identified in Section 3.2 are accessible via control over defects. The dopability window of a compound is established from the formation energies of native defects at each of the two invariant points of the stability region. For the binary Sn chalcogenides, Sn-poor conditions are most favorable for acceptor Sn vacancies V_Sn to form and therefore represent conditions that promote larger p-type carrier concentrations. In contrast, Sn-rich (anion-poor) conditions are most favorable for anion donor vacancies to form, and therefore most suitable when trying to achieve large n-type carrier concentration. By investigating defect formation at these two extremes, the complete accessible range of carrier concentrations under equilibrium can be established for each compound.

Defect formation energies ΔH _D,q are typically shown as a function of semiconductor Fermi level E _F , which varies from the valence band maximum (VBM) to conduction band minimum (CBM). Accurate chemical potentials and band gaps are prerequisites to obtaining accurate defect formation energies. In addition to discrepancies between experiment and DFT-HSE in formation energies, there are also discrepancies between DFT-HSE and experimental band gaps. The band gaps obtained using DFT-HSE for SnS, SnSe, and SnTe are 1.31, 1.09, and 0.18 eV while the experimental band gaps are 1.2 eV (Albers et al., 1960; Zhou et al., 2017), 0.9 eV (Lefebvre et al., 1998), and 0.18 eV (Rogers, 1968; Dimmock et al., 1966), showing differences that lie within the expected range of HSE band gap uncertainties.(see Supplementary Table S1).

Carrier concentrations are exponential functions of defect formation energies, so small errors in defect formation energies can translate to large errors in predicted carrier concentrations. In the following, predicted carrier concentrations are shown for two cases: 1) the pure DFT-HSE prediction (labeled “pure DFT”), and 2) the DFT-HSE prediction with empirical corrections for chemical potentials and band gaps (labeled “DFT + Exp”). For the latter, we adjusted the band gap by applying equal shifts to both band edge positions to match the experimental band gap. We further shifted the calculated total energies of each Sn binary to match the experimental formation energy. The comparison of the two sets of results allows us to estimate the degree of uncertainty in computed carrier concentrations.