Intrinsic properties and dopability effects on the thermoelectric performance of binary Sn chalcogenides from first principles

High-performance thermoelectric (TE) materials rely on semiconductors with suitable intrinsic properties for which carrier concentrations can be controlled and optimized. To demonstrate the insights that can be gained in computational analysis when both intrinsic properties and dopability are considered in tandem, we combine the prediction of TE quality factor (intrinsic properties) with first-principles simulations of native defects and carrier concentrations for the binary Sn chalcogenides SnS, SnSe, and SnTe. The computational predictions are compared to a comprehensive data set of previously reported TE figures-of-merit for each material, for both p-type and n-type carriers. The combined analysis reveals that dopability limits constrain the TE performance of each Sn chalcogenide in a distinct way. In SnS, TE performance for both p-type and n-type carriers is hindered by low carrier concentrations, and improved performance is possible only if higher carrier concentrations can be achieved by suitable extrinsic dopants. For SnSe, the p-type performance of the Cmcm phase appears to have reached its theoretical potential, while improvements in n-type performance may be possible through tuning of electron carrier concentrations in the Pnma phase. Meanwhile, assessment of the defect chemistry of SnTe reveals that p-type TE performance is limited by, and n-type performance is not possible due to, the material’s degenerate p-type nature. This analysis highlights the benefits of accounting for both intrinsic and extrinsic properties in a computation-guided search, an approach that can be applied across diverse sets of semiconductor materials for TE applications.

The discovery of new thermoelectric (TE) semiconductors is most commonly based on predicting or measuring the figure of merit zT (Madsen, 2006; Zhu et al., 2015; Fang et al., 2017; Li et al., 2018; Chami et al., 2020; Gan et al., 2021; Fu et al., 2015; Zhao et al., 2014). While zT is effective for evaluating a material’s TE performance, it is a composite parameter that arises from intrinsic material properties such as carrier effective masses and thermal conductivity, as well as extrinsic properties such as temperature and carrier concentration (Snyder and Toberer, 2008; Gorai et al., 2017). Consequently, measurements or predictions of zT alone may mask key performance bottlenecks. If a semiconductor is shown to have low zT, it is not always obvious if improvements are needed in the material itself or if the carrier concentration should be further optimized by doping. Additionally, semiconductors that exhibit high carrier mobilities and low thermal conductivity can still be poor TEs if they cannot be doped to the desired carrier concentration, as in n-type GeTe (Samanta et al., 2019) and Sr 3 GaSb 3 (Zevalkink et al., 2012). For accelerated TE materials discovery and optimization, the roles of intrinsic material properties and dopability should be considered in tandem and their effects disambiguated.

One well-known and technologically-relevant class of TE materials that exhibit a diversity of behaviors in intrinsic properties and dopability are the binary IV-VI chalcogenides. These compounds exhibit crystal structures including orthorhombic (GeS and GeSe—space group Pnma, SnS and SnSe—space groups Pnma and Cmcm), rhombohedral (GeTe—space group R3m), and rock salt (SnTe, PbS, PbSe, and PbTe—space group Fm3m). As this family of materials already exhibits anharmonicity-induced ultra-low thermal conductivities (Lee et al., 2014), opportunities for improvement lie in optimizing electronic properties and carrier concentrations. Experimentally it is observed that PbTe (Wu et al., 2015) achieves a zT of 2.3 and GeTe (Li et al., 2018) achieves a zT of 2.4 both at hole concentrations ∼ 1 0 20 cm −3, while SnSe (Zhao et al., 2014) shows higher performance (zT = 2.6) at a lower hole concentration ∼ 1 0 19 cm −3. Experimental work on binary IV-VI chalcogenides has historically focused on p-type compositions, yielding state-of-the-art performances, e.g. PbTe (Pei et al., 2011a; Wu et al., 2015; Hu et al., 2016; Biswas et al., 2012), GeTe (Li et al., 2016; Pei et al., 2011b; Samanta and Biswas, 2017; Adamczyk et al., 2022), SnSe (Zhao et al., 2014; 2015). More recently, their n-type performances have been characterized as well, such as SnSe (zT = 2.2, Duong et al., 2016) and PbTe (zT = 1.8, Fu et al., 2017). Detailed analysis of dopability using both experiment and computation (Goyal et al., 2017; Male et al., 2019) established the route for n-type doping of PbTe. These studies invite questions about whether optimal performance has been achieved, or if any of the IV-VI binaries have potential for TE improvement by doping optimization. Assessing this possibility requires analysis of both intrinsic properties and dopability, to identify the most effective routes to increase zT.

Amongst a variety of methods that predict TE performance (Madsen and Singh, 2006; Madsen et al., 2018; Ganose et al., 2021), the thermoelectric quality factor approach (Chasmar and Stratton, 1959) (Figure 1) obtains zT in terms of two independent variables, the quality factor β and the reduced Fermi level η. Since β is a measure of intrinsic properties and η depends on doping and temperature, the quality factor approach cleanly separates zT into its intrinsic and extrinsic constituents. The development of semi-empirical models (Yan et al., 2015; Miller et al., 2017) to evaluate carrier mobilities and thermal conductivity (Figure 1A) has made high-throughput assessment of β possible (Gorai et al., 2015; Ohno et al., 2018; Male et al., 2019; Qu et al., 2020). The dopability of a material can be evaluated through first-principles analysis of defect formation energies (Freysoldt et al., 2014) to establish the limits of achievable n-type and p-type carrier concentration and η (Figure 1B). As shown in Figure 1C), the union of intrinsic properties and dopability yield the achievable zT for a given material. The peak values observed in plots of zT vs. carrier concentration typically occur at carrier concentrations between 10 19 –10 20 cm −3, making a defect engineering approach essential. In materials where compensating defects are prevalent, the peak zT may not be accessible, leading to wasted time and effort. Compensating defects refer to low energy native defects in a semiconductor that counteract attempts to dope the material; for example it is not possible to dope CuInTe 2 n-type due to the prevalence of compensating copper vacancies V Cu that counteract the doping (Adamczyk et al. (2020)). Additionally, the more ideal the intrinsic properties the less stringent the demands on carrier concentration: the larger the quality factor β, the lower the carrier concentration at which zT peaks and the higher the peak value.

To highlight the effects of assessing both intrinsic properties and dopability in computational analysis, in this work we apply the TE quality factor approach to the prediction of zT in well-known Sn-based IV-VI binary compounds SnS, SnSe, and SnTe. Amongst the IV-VI binaries, we focus on Sn-based compounds due to their varying crystal structure (Pnma, Cmcm, and Fm3m) as well as recent demonstrations of n-type performance (Duong et al., 2016; Cai J et al., 2020; Gu et al., 2021; Zhang et al., 2015; Cha et al., 2019; Hu et al., 2019; Zhi-Cheng et al., 2019). We first apply the TE quality factor formalism for both p-type and n-type carriers by estimating β for each compound to predict the dependence of zT on carrier concentration. Our results are compared to a comprehensive literature search of reported zT values for each compound, which reveals opportunities for improvement for all three materials via carrier concentration refinement. We next carry out analysis of achievable carrier concentrations for each, and find that native defects and dopability restrict the performance of all three materials differently. Even with the limitations imposed by carrier concentration analysis, some of the Sn chalcogenides show opportunities for improvement if good extrinsic dopants can be found. The application of the quality factor approach highlights the advantages of distinguishing between the roles of intrinsic properties and dopability, as an avenue for accelerated optimization of thermoelectric materials.

2.1 DFT simulations—General approach

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 10:5s 2 5p 2, S: [Ne]:3s 2 3p 2, Se: [Ar]4d 10:4s 2 4p 4, Te: [Kr]4d 10:5s 2 5p 4. 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.

2.2 Thermoelectric quality factor β

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.

2.3 Dopability analysis

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.

2.3.1 Phase stability

The thermodynamic stability of each binary Sn chalcogenide was determined against all competing phases, including elemental phases (Sn, S, Se, and Te) and binaries SnS 2 and SnSe 2, whose total energies were obtained via geometry optimization as described above. The range of chemical potentials Δ μ for which the parent compound is stable was obtained from the total energies, by requiring that the compound of interest be stable with respect to decomposition to all possible competing phases.

2.3.2 Native defect formation energies and carrier concentrations

We used the standard supercell approach (Lany and Zunger, 2008) to determine the formation energies of native point defects to establish dopability. The formation energy Δ H D, q of point defect D in charge state q is given by where E D, q and E H represent total energies of a defect supercell and a host supercell without defects, respectively; μ i is the chemical potential of element i added (n i > 0) or removed (n i < 0) from the host to create the defect; E F is the Fermi energy varying from the VBM to CBM; and E corr refers to the finite size corrections within the supercell approach. The corrections included in evaluating Δ H D, q were (i) potential alignment corrections and (ii) image charge corrections for charged defects as described by Lany and Zunger (2008). To evaluate image charge corrections, the electronic and ionic dielectric constant was calculated using density functional perturbation theory (DFPT) as implemented in VASP (Kresse and Furthmüller, 1996).

For defect calculations of SnS, SnSe, and SnTe we built supercells of 72, 72, and 64 atoms, respectively, and relaxed the structure with HSE06 to calculate the total energies of the supercells. The supercells were relaxed using Brillouin zone sampling with a Γ-centered 2 × 2 × 2 k-point grid.

Defect concentration, carrier concentration, and equilibrium Fermi energy were determined by assuming equilibrium defect concentrations and imposing charge neutrality. The charge neutrality condition is given by where q is the charge state of the defects, n and p are the free electron and hole concentrations, and C D, q is the defect concentrations, respectively. C D, q is obtained by where N is the concentration of the lattice sites, k B is the Boltzmann constant, and H D, q is the defect formation energy.

In Eq. 2, n and p can be obtained by where D C (ϵ), D V (ϵ) and f (ϵ) are the conduction band density of states, valence band density of states, and the Fermi-Dirac distribution function, respectively. For a non-degenerate semiconductor the carrier concentration can be analytically simplified by the parabolic band approximation as where m DOS * represents DOS effective masses for electrons and holes in the two expressions, respectively. In general, we used Eq. 5 to calculate intrinsic carrier concentration using the calculated DOS effective mass for electrons and holes. Eq. 4 was used for SnTe and p-type SnS, since we found degenerate doping for these cases (see Section 3.3).

3.1 TE quality factor assessment

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 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 σ 0 the “intrinsic” electrical conductivity (Chasmar and Stratton, 1959; Yan et al., 2015). The electronic conductivity σ = σ 0 exp(η), where η = −Δ/ k b T and Δ is the positive distance between the Fermi level and the band edge. The quantity σ 0 is determined by the intrinsic mobility μ 0 and density of states (DOS) effective mass m DOS *. Putting everything together, the expression for β (T) is The intrinsic mobility μ 0 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 μ 0 and κ L at T = 300 K are available (Yan et al., 2015; Miller et al., 2017). The model for μ 0 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 zT s (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).

3.2 zT vs carrier concentration: Computation and experiments

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 μ 0 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 18 –10 19 cm −3. 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 −3).

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 zT s 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 −3, 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 20 –10 21 cm −3 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 −3. SnSe in the Cmcm phase and SnTe have high β and so, the maximum TE potential can be achieved around lower carrier concentrations of 10 19 cm −3 and 10 18 cm −3 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 19 –10 20 cm −3 (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 19 –10 20 cm −3. However, for SnSe (Pnma) and SnTe we find that a carrier concentration of 10 19 cm −3 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 2 /(Vs) while the experimental hole mobilities are 80, 170, and 650 cm 2 /(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

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).

3.3 Dopability predictions

3.3.1 Phase stability

The first step for dopability analysis is to establish the phase stability limits of SnS (Pnma), SnSe (Pnma), and SnTe (Fm3m). Thermodynamic stability is set by an allowable range of chemical potentials Δ μ i (T, P) that give the availability of each element i, with more (or less) negative values indicating relative scarcity (or abundance) in the thermodynamic environment. Outside of this range, the compound becomes unstable with respect to the formation of competing boundary phases. Beyond determining stability limits, by controlling chemical potentials during growth it is possible to tune defect concentrations and gain better control over both carrier type and quantity, dopability, and transport (Ohno et al., 2018; Male et al., 2019; Ortiz et al., 2019).

The chemical potential stability regions of SnS, SnSe and SnTe are indicated by the teal lines in Figure 5, as a function of the chemical potential Δ μ Sn. For binary Sn compounds, fixing the chemical potential Δμ Sn correspondingly fixes the chemical potential of the remaining element. For SnS for instance where Δ H f (SnS) is the enthalpy of formation of SnS. The stability region for each compound is therefore shown as a one-dimensional line in Figure 5 that spans the allowed range of Sn-poor to Sn-rich conditions. The line is bounded on both ends by “invariant points,” fixed points that correspond to two-phase equilibria between the compound of interest and a second, competing phase. The competing boundary phases that limit the stability of SnS are elemental Sn (under Sn-rich conditions) and SnS 2 (under S-rich conditions). Elemental Sn and SnSe 2 limit the stability of SnSe for Sn-rich and Se-rich conditions, respectively. For SnTe, the stability region is bounded by Sn and Te for Sn-rich and Te-rich conditions respectively, as there are no other competing phases in the Sn-Te binary phase diagram.

We found formation energies from DFT-HSE for SnS, SnSe, and SnTe to be Δ H f = −0.43, −0.44, and −0.28 eV per atom, while the reported experimental formation energies are Δ H f = −0.57, −0.47, and −0.32 eV per atom (Stevanović et al., 2012; Castanet et al., 1972). The difference between DFT-HSE and experimental reports is within 0.14 eV/atom for all three compounds, which is within the expected range of DFT uncertainties (Stevanović et al., 2012; Kim et al., 2016). Considering the competing phases of all three compounds the chemical potential boundaries are −0.43 eV ≤ Δ μ Sn ≤ 0 eV, −0.60 eV ≤ Δ μ Sn ≤ 0 eV and −0.55 eV ≤ Δ μ Sn ≤ 0 eV for SnS, SnSe and SnTe respectively ranging from Sn-poor to Sn-rich growth conditions (Figure 5). This accessible range of chemical potentials from Sn-poor to Sn-rich growth conditions leads to variations in defect and carrier concentrations and establishes the dopability window for each compound, as described below.

3.3.2 Approach to dopability prediction

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.

3.3.3 Defect chemistry of SnTe

Formation energies of native point defects in SnTe in both Sn-rich and Sn-poor growth conditions are shown in Figures 6A,B. The native defects include vacancies V Sn and V Te and antisite defects Sn Te and Te Sn, all of which were considered in varying charge states. Following convention, we show only the minimum energy charge state of each defect in Figure 6. Interstitial defects were originally considered for all three Sn chalcogenides at the PBE level, but were found to have higher formation energies compared to vacancies and antisites (Supplementary Figure S4); therefore we did not include them in the DFT-HSE results for dopability analysis.

Under both Sn-rich (Figure 6A) and Sn-poor (Figure 6B) growth conditions, Sn vacancies (V Sn) act as the dominant acceptor defect with ultra-low formation energy. The next lowest energy defects are the two antisites, which flip their relative energies between Sn-rich and Sn-poor conditions. The only donor defect, V Te, shows high formation energy across the full stability region. The ultra-low formation energy of acceptor V Sn for both extremes of growth conditions indicates that V Sn will be the dominant defect in SnTe across the full stability region. This result is consistent with the degenerate p-type behavior that is experimentally reported in SnTe (Doi et al., 2019; Ma et al., 2019; Tang et al., 2018; Tan et al., 2015). It is also consistent with the experimental observation that it is nearly impossible to synthesize stoichiometric SnTe (Zhou et al., 2014), with materials always showing cation deficiency. As also suggested in Wang et al. (2014), we believe that the reason for the ultra-low formation energy of V Sn is related to the high band edge position of the VBM of the Fm3m phase of SnTe (on an absolute scale, in comparison to the subsequent compounds to be discussed below). For example, we found the formation energy of V Sn and other acceptor defects to be comparatively higher in a hypothetical Pnma phase of SnTe and more similar to the Pnma phases of SnS and SnSe (see Supplementary Figure S5).

For SnTe, the previous analysis from Figure 3D indicated the possibility of improved zT if p-type carrier concentrations could be reduced. However, Figure 6 shows that this is not possible and provides an explanation for the large p-type carrier concentrations around 10 20 cm −3 in Figure 3D. For both Sn-rich and Sn-poor cases, we found that the equilibrium Fermi energy E F e q resides inside the valence band: from charge neutrality constraints, the calculated E F e q ’s are 0.13 and 0.55 eV below the VBM at 400 K for Sn-rich to Sn-poor condition respectively. The corresponding hole and electron concentrations for these two extremes are shown as a function of temperature from 400 K to 900 K in Figure 6C. Under the Sn-rich conditions most amenable to n-type material, the hole carrier concentration is around 10 19 cm −3 and the electron carrier concentration ranges from 1.22 × 10 14 to 6.24 × 10 16 cm −3. At Sn-poor growth conditions the hole carrier concentration increases to ∼ 1 0 20 cm −3, with negligible small electron carrier concentrations. Therefore, despite that the p-type and n-type TE performance of SnTe is predicted to exceed that of PbTe according to material descriptor β, defect analysis shows that the reduction of hole carriers (or realization of n-type material) is not feasible. The only reported n-type SnTe was found to occur upon substantial alloying with PbTe, as a result of lowering the position of the valence band edge (Pang et al., 2021). Therefore, combined analysis shows that improved p-type and good n-type TE performance will be challenging, if not impossible, in SnTe even in highly alloyed compositions.

Although the defect analysis suggests that tuning carrier concentrations in bulk material is challenging, other strategies for engineering SnTe may still be possible. For example, alloying with PbTe (Pang et al., 2021) and CaTe (Al Rahal Al Orabi et al., 2016) may help push the VBM lower, increasing V Sn formation energies and reducing carrier concentrations. Alternatively, recent approaches to controlling disorder have been suggested as routes to TE performance optimization (Roychowdhury et al., 2021; Agne et al., 2021). In these cases, the use of disorder-induced carrier localization could be used for TE improvement via tuning of carrier mobilities in order to increase Seebeck coefficient while maintaining electrical conductivity. This strategy might also be helpful to reduce lattice thermal conductivity through nanostructuring as shown in Ref. Roychowdhury et al. (2021). Finally, considerations of alternative scattering mechanisms or more sophisticated treatments of scattering, such as the four coefficient approach (Crawford et al., 2021), have been found to yield different optimized values for parameters such as carrier concentrations.

3.3.4 Defect chemistry of SnS

For SnS, it is of interest to determine whether it is possible to increase carrier concentrations to improve TE performance, particularly for n-type material. The formation energies of native point defects and their applicable charge states under Sn-poor and Sn-rich growth conditions from pure DFT are shown in Figures 7A,B (the empirically corrected diagrams are shown in the Supplementary Figure S6A). The equilibrium Fermi energy E F e q at 400 K is shown by the dashed line. Distinct from SnTe, now in both growth conditions V Sn and V S act as the main acceptor and donor defects respectively. These two defects compete, and their relative formation energies at each invariant point determine the position of E F e q. In Sn-poor conditions, V Sn has a relatively low, and V S a relatively high, formation energy giving rise to more p-type conditions. In Sn-rich conditions, V Sn has a relatively high, and V S a relatively low, formation energy giving rise to more n-type conditions. The equilibrium Fermi energy at 400 K varies between 0.61 and 0.68 eV above the VBM within the chemical potential range. Overall, since the cation vacancy V Sn no longer dominates as the governing defect as in SnTe, a larger dopability window for both p-type and n-type material is present.

However, SnS also has a wider band gap compared to SnTe (1.31 eV vs. 0.18 eV), making it more difficult to achieve high carrier concentration. To assess the maximum p-type and n-type carrier concentration achievable in SnS, we calculated both the intrinsic and the ideal carrier concentration in Sn-poor and Sn-rich growth conditions. The intrinsic carrier concentration refers to the carrier concentration that can be achieved without the introduction of extrinsic dopants. It arises from the range of accessible equilibrium Fermi levels E F e q throughout chemical potential space. On the other hand, the ideal carrier concentration describes the maximum achievable p-type or n-type carrier concentration throughout chemical potential space. As indicated in Figure 7B, it is found by assuming the presence of an ideal extrinsic dopant that pushes the Fermi level to where the lowest energy compensating defect (here, an acceptor in n-type material) has formation energy of zero. Although in practice it is difficult or impossible to identify an ideal extrinsic dopant, the ideal dopability serves as an upperbound to the achievable carrier concentration (Qu et al., 2021).

The predicted intrinsic and ideal hole carrier concentrations are shown in Figure 7C. Light colors indicate intrinsic and dark colors ideal hole concentrations at Sn-poor conditions. The intrinsic hole carrier concentration at 923 K from pure DFT is predicted to be 2.8 × 10 17 cm −3, and increases to 9.8 × 10 17 cm −3 when empirical corrections to band gap and chemical potentials are applied. These values match fairly closely with reported experimental measurements for undoped material which range from 4 × 10 17 –10 18 cm −3 at similar temperatures (Tan et al., 2014; Abraham et al., 1978). As an upperbound, with an ideal p-type dopant, hole carrier concentrations on the order of 10 20 cm −3 would be achievable. Reported carrier concentrations for doped or alloyed SnS range between 2 × 10 18 to 3 × 10 19 cm −3 (Zhou et al., 2017; Tan et al., 2014; He et al., 2019; Cai B et al., 2020). Our combined analysis using β and dopability predictions suggest that there is opportunity to achieve higher p-type zT through improving hole carrier concentration by extrinsic doping.

To assess achievable n-type performance, predicted intrinsic and ideal electron carrier concentrations are shown in Figure 7D. Compared to hole carrier concentrations in Figure 7C, achievable electron carrier concentrations are generally lower. For intrinsic material again at 923 K, we predict carrier concentrations of 1.4 × 10 16 cm −3 for pure DFT and 3.6 × 10 16 cm −3 for empirically-corrected DFT, about an order of magnitude lower than the corresponding p-type values. The predicted ideal electron concentration suggests that with an ideal n-type dopant, electron carrier concentrations on the order of 10 19 cm −3 might be achievable suggesting that SnS is n-type dopable. To date, n-type carrier concentrations in extrinsically doped or alloyed SnS lie around ∼ 1 0 17 cm −3 (Zhi-Cheng et al., 2019; Hu et al., 2019). Similar to p-type SnS, the combined intrinsic property and dopability analysis here suggests that there is an opportunity for improvement in n-type zT if suitable extrinsic dopants can be found.

3.3.5 Defect chemistry of SnSe

Lastly, focusing on n-type performance, we analyze the dopability of SnSe in the Pnma phase to evaluate whether it is possible to tune carrier concentration to improve TE performance. Only the defect chemistry at Sn-rich growth condition is shown in Figure 8A from pure DFT (the empircally corrected diagrams for both Sn-rich and Sn-poor condition are shown in the Supplementary Figure S6B). Similar to SnS, V Sn and V Se act as the primary acceptor and donor defects respectively, and their relative energies determine the position of E F e q. In Sn-rich conditions, the equilibrium Fermi energy is located at 0.59 eV above the VBM at 400 K. Moreover, SnSe has a smaller band gap than SnS, making it easier to achieve a higher carrier concentration.

Predicted intrinsic and ideal electron carrier concentrations are shown in Figure 8B. For intrinsic material at 800 K, we predict carrier concentrations of 9.6 × 10 15 cm −3 for pure DFT and 6.11 × 10 16 cm −3 for empirically-corrected DFT. Here the order of magnitude difference arises from the ∼0.2 eV band gap difference between DFT-HSE and experiment. The predicted ideal electron concentration suggests that electron carrier concentrations ranging between 10 19 to 10 20 cm −3 might be achievable from pure and empirically-corrected DFT, respectively. To date, n-type carrier concentrations in extrinsically doped or alloyed SnSe lie on the order of ∼ 1 0 19 cm −3 (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) which suggests that slight improvements in n-type zT might be possible in SnSe.

In this work, we studied the TE quality factor approach to identify improvement opportunities in TE performance in the well-known Sn-based IV-VI binary compounds SnS, SnSe, and SnTe focusing on intrinsic properties and dopability. First, we predicted zT using the TE quality factor formalism through estimations of parameter β. The predicted value of β for n and p-type versions of all three materials results in distinct target carrier concentrations for each material, with SnS requiring larger, and SnTe requiring lower, carrier concentrations. Although comparison of these predictions to a comprehensive literature search revealed opportunities for improvement in all three Sn chalcogenides, the subsequent dopability analysis showed that the defect chemistry of each material constrains its TE performance in a distinct way. In SnS, improved TE performance in both p-type and n-type material can be achieved if suitable extrinsic dopants can be found. For SnSe, the p-type performance in the Cmcm phase appears to have reached its potential, but improved n-type TE performance can be achieved in the Pnma phase through proper tuning of electron carrier concentration. For SnTe both p-type and n-type performance is constrained by the degenerate p-type carrier concentrations arising from the abundance of low energy V Sn defects. Through this work, we aim to show the importance of carefully considering both intrinsic and extrinsic properties in a computation-guided search, as a general strategy to be utilized across diverse arrays of materials for TE application. It is of future interest to carry out comprehensive quantitative comparisons between computationally-predicted and experimentally-measured performance in diverse families of materials as demonstrated here for binary chalcogenides as a way to benchmark the accuracy attainable and typical uncertainties in first-principles approaches to TE materials.