Thermoelectric Properties of Pnma and Rocksalt SnS and SnSe

: Thermoelectric materials convert waste heat to electricity and are part of the package of technologies needed to limit global warming. The tin chalcogenides SnS and SnSe are promising candidate thermoelectrics, with orthorhombic SnSe showing some of the highest ﬁgures of merit ZT reported to date. As for other Group IV chalcogenides, SnS and SnSe can form rocksalt phases under certain conditions, but the thermoelectric properties of these phases are largely unexplored. We have applied a fully ab initio modelling protocol to compare the ZT of the orthorhombic and rocksalt phases of SnS and SnSe. Electronic structures from hybrid density-functional theory were used to calculate the three electrical transport properties, including approximate models for the electron relaxation times, and lattice dynamics calculations were performed to model the phonon spectra and lattice thermal conductivities. We obtained good estimates of the ZT of the well-studied orthorhombic phases. The rocksalt phases were predicted to show larger electrical conductivities and similar Seebeck coefﬁcients to the orthorhombic phases, resulting in higher thermoelectric power factors, but these were offset by larger thermal conductivities. These results therefore motivate further investigation of the recently discovered “ π -cubic” phases of SnS and SnSe, which are based on distorted rocksalt supercells, to establish their thermoelectric performance.


Introduction
The growing worldwide demand for energy and the urgent need to reduce anthropogenic greenhouse gas emissions to limit global warming are among the most important challenges in contemporary science. Meeting ambitious climate goals requires a package of technologies including clean energy sources, energy storage systems, and means to improve the efficiency of energy-intensive processes. It is estimated that around 60% of the energy used globally is currently wasted as heat [1], which has led to significant interest in thermoelectric (TE) power.
Thermoelectric generators (TEGs) harness the Seebeck effect in a TE material to extract electrical energy from a temperature gradient [2]. TEGs are solid-state devices with no moving parts and, as such, have diverse applications from powering remote sensing devices, to recovering energy from exhaust gases in combustion engines, to repurposing decommissioned oil rigs as geothermal power plants [2,3].
The performance of a TE material is typically expressed by the dimensionless figure of merit ZT [1,2]: where S is the Seebeck coefficient, σ is the electrical conductivity, and κ el and κ latt are the electronic and lattice (phonon) components of the thermal conductivity κ. S, σ, and κ el depend on the electronic structure of the material and are interdependent through the carrier concentration n, such that the best balance is typically found in heavily-doped initio protocol for predicting the ZT yields results for the orthorhombic phases that are in good agreement with experiments. We predict that the rocksalt phases would show higher electrical conductivities and comparable Seebeck coefficients to the orthorhombic phases, resulting in higher thermoelectric power factors S 2 σ, but that these would generally be offset by larger electrical and lattice thermal conductivities. Our results suggest that if RS SnS can be stabilised at a geometry close to its equilibrium lattice constant it may show a larger ZT than the Pnma phase, and this strongly motivates further investigation of the thermoelectric properties of the rocksalt-derived π phases of both chalcogenides.

Computational Modelling
Calculations were performed using pseudopotential plane wave density functional (DFT) as implemented in the Vienna Ab initio Simulation Package (VASP) code [47].
Electron exchange and correlation were described using the PBEsol generalised gradient approximation (GGA) functional [48] with the DFT-D3 dispersion correction [49]. The ion cores were modelled using projector-augmented wave (PAW) pseudopotentials [50,51] with the Sn 5s, 5p, and 4d and the S 3s/3p and Se 4s/4p electrons treated as valence states. The valence wavefunctions were described using a plane wave basis with a 600 eV kinetic energy cutoff. The electronic Brillouin zones were sampled using Γ-centred Monkhorst-Pack k-point meshes [52] with 4 × 8 × 8 subdivisions for the Pnma structures and 10 × 10 × 10/6 × 6 × 6 subdivisions for the rocksalt primitive/conventional cells. The electronic wavefunctions were optimised to a tolerance of 10 −8 eV on the total energy.
Optimised equilibrium structures of Pnma and RS SnS and SnSe were taken from our previous study [45]. Since the equilibrium structure of RS SnS possesses phonon instabilities [44,45], we also investigated a compressed structure with a 5% smaller lattice constant (14.3% smaller volume) for which the imaginary mode becomes real.
Lattice dynamics calculations were performed using the Phonopy package [53]. Secondorder (harmonic) force constants, computed using the supercell finite-differences approach with 1 × 6 × 6 expansions of the Pnma unit cells and 3 × 3 × 3 expansions of the rocksalt conventional cells, were taken from our previous work [45]. A transformation matrix was applied to convert from the rocksalt conventional to the corresponding primitive cell during post-processing. Atom-projected phonon density of states (PDoS) curves were computed by interpolating the phonon frequencies onto regular Γ-centred q-point grids with 8 × 48 × 48 and 32 × 32 × 32 subdivisions for the Pnma and RS phases, respectively, and using the linear tetrahedron method for Brillouin zone integration. Phonon dispersions were computed by evaluating the frequencies at strings of q-points passing through the high-symmetry points in the Pnma and RS Brillouin zones.
Electronic transport calculations were performed using the AMSET code [54]. Initial electronic structure calculations were performed using the HSE 06 hybrid functional [55] to obtain accurate bandgaps. Uniform band structures and sets of Kohn-Sham wavefunction coefficients were then computed using PBEsol + D3 and denser 8 × 16 × 16 and 20 × 20 × 20 k-point meshes for the Pnma and RS structures, respectively, and the bandgaps increased to the HSE 06 values using appropriate scissors operators. AMSET computes electronic relaxation times by summing scattering rates from four different processes, namely acoustic deformation potential (ADP), piezoelectric (PIE), polar optical phonon (POP), and ionised impurity (IMP) scattering. For ADP scattering, deformation potentials were computed by performing a series of single-point energy calculations on deformed structures generated using AMSET with HSE 06, and elastic constants were computed using PBEsol + D3 and the finite-differences routines in VASP. PIE scattering is not relevant to any of the materials considered here as the Pnma and RS (Fm3m) space groups are centrosymmetric and the piezoelectric moduli vanish. For POP scattering, high-frequency and static dielectric constants ε ∞ and ε s = ε ∞ + ε ionic were determined using the density-functional perturbation theory (DFPT) and finite-difference routines in VASP [56], and the POP frequencies ω po were determined as a weighted average of the phonon frequencies at q = Γ, obtained using Phonopy, with the infrared (IR) activities computed using Born effective charges Z * from DFPT [57]. As described in the text, for the RS phases, we also computed ε ∞ , ε s , and Z * using HSE 06, with the finite-field method used to determine ε ∞ and ε s [58,59]. Finally, for IMP scattering, the required ε s were determined as for the POP scattering. For all of these calculations, we switched to Sn PAW pseudopotentials with the Sn 4d electrons in the core. We found this had no significant impact on the calculated electronic structures or bandgaps, but by reducing the number of valence electrons per formula unit by a factor of two it made the HSE 06 calculations significantly less demanding.
Thermal conductivity calculations were performed using the Phono3py code [60]. The thermal conductivities of Pnma SnS and SnSe were taken from another of our previous studies, which used a very similar technical setup [61]. A comparison of the optimised lattice parameters computed for these two structures in this work and [61] is provided as Supplemental Materials. Third-order (anharmonic) force constants for the RS models were computed in 2 × 2 × 2 expansions of the conventional cell and combined with the second-order force constants from the larger 3 × 3 × 3 expansions. The force constants were calculated using PBEsol + D3 and with additional support grids with 8 × the number of points as the standard charge density grids to ensure accurate forces. The thermal conductivities were then computed from modal properties evaluated on 28 × 28 × 28 sampling meshes.

Structure and Lattice Dynamics
Representative structures of the orthorhombic Pnma and cubic rocksalt (RS) structures of the tin monochalcogenides are shown in Figure 1. The RS structure has an octahedral bonding environment in which each Sn(II) cation has six nearest-neighbour chalcogen ions with equal Sn-Ch bond lengths. The Pnma structure can be thought of as a distortion of the RS structure where alternate layers are misaligned along one direction to produce a pseudo-2D layered structure. Within each layer, the Sn(II) cations adopt a distorted local geometry with three bonds to neighbouring chalcogen ions and a stereochemically active Sn 5s lone pair that projects into the interlayer spacing to facilitate a dispersive (van der Waals) interaction between layers. The optimised lattice parameters of the five structures obtained with PBEsol + D3 are collected in Table 1. The Pnma lattice constants are a good match for the 295 K neutronscattering measurements in [28], namely a = 11.143, b = 3.971, and c = 4.336 Å for SnS (0.15-3.1% smaller) and a = 11.501, b = 4.153, and c = 4.445 Å for SnSe (0.7-2.5% smaller). These discrepancies may in part be due to the fact that the DFT calculations are "athermal", i.e., the optimised structures are those at 0 K without corrections for the vibrational zeropoint energy, whereas the experimental measurements at finite temperature may include a small amount of thermal expansion. It is interesting to note that the mismatch between the measured and optimised lattice parameters is anisotropic and is consistently smallest along the b axes and largest along the c axes. As seen in molecular dynamics simulations, the interatomic distances along this direction are substantially affected by the incipient Pnma → Cmcm phase transition, which, given its second-order nature, leads to structural changes that begin to manifest even at room temperature [63]. Nonetheless, we still consider the discrepancies with experimental measurements to be acceptable, and we therefore do not consider this to be too big an issue.
While we would not necessarily expect a good match between optimised bulk lattice constants and those of epitaxial thin films, the optimised lattice constants for equilibrium and compressed RS SnS are 1.5 smaller than the value of 5.8 Å measured for RS SnS grown epitaxially on NaCl in [38], and the optimised lattice constant for RS SnSe is only 1.3% smaller than the 5.99 Å reported in [37]. The calculated phonon dispersion and atom-projected density of states (PDoS) curves of the five structures are shown in Figure 2. The Pnma structure has eight atoms in the unit cell, resulting in 24 branches at each phonon wavevector q, whereas the highersymmetry rocksalt structure has a smaller two-atom primitive cell and only six branches in the dispersion. The DoS curves show distinct low-and high-frequency components corresponding respectively to the motion of the Sn and chalcogen atoms, resulting in the appearance of a "phonon bandgap" in the Pnma dispersions. The separation is more pronounced in the sulphides than in the selenides, which reflects the smaller mass difference and more covalent bonding in the selenides. The upper limit of the frequency spectrum is 25 cm −1 higher in the Pnma phases of both chalcogenides than in the corresponding equilibrium RS structures, although compression of RS SnS leads to a stiffening of the optic modes, which is most noticeable around q = Γ.
As found in previous studies [44,45], the equilibrium structure of RS SnS shows a prominent imaginary mode at q = X, which hardens and becomes real under compression, whereas RS SnSe is predicted to be dynamically stable. (We note that the dispersion of compressed SnS shows a small imaginary mode along the Γ-W path, which we ascribe to an interpolation artefact as a result of the finite supercell size used to compute the force constants.) The contrasting stability of the equilibrium rocksalt SnS and SnSe structures can be explained using the revised lone pair model proposed by Walsh et al. [64] The interaction between the Sn 5s and chalcogen p orbitals produces high-energy antibonding orbitals that can interact with the Sn 5p and indirectly mediate a coupling between the Sn 5s and 5p states. This interaction is forbidden by symmetry in the octahedral coordination coordination environment in the RS phase, but is supported by the distorted local geometry in the Pnma structure. The lower-energy Sn 3p orbitals are more closely matched in energy to the Sn 5s, which results in a strong orbital interaction and a stronger propensity to distort away from the symmetric RS structure. Conversely, the poorer energy match between the Sn 5s and Se 4p orbitals results in a weaker preference for the distorted Pnma structure, so RS SnSe is both dynamically stable and considerably closer in energy to the orthorhombic phase [45]. The same structural chemistry is seen in the bismuth chalcogenides, where Bi 2 S 3 preferentially adopts a low-symmetry orthorhombic Pnma phase while Bi 2 Se 3 and Bi 2 Te 3 are more stable in a rhombohedral R3m phase with the Bi cations in an octahedral bonding environment [13]. On the other hand, compression shortens the Sn-S bond lengths and introduces a barrier to the structural distortion, which causes the imaginary mode to harden and become real.

Electronic Structure and Transport Properties
The general features of the electronic band dispersion and density of states (DoS) curves of the Pnma and rocksalt phases of SnS and SnSe were discussed in our previous work [45], and we provide the calculated HSE 06 band structures of the five structures examined in this work as Supplemental Materials. The two Pnma structures are indirect gap semiconductors with the valence band maximum and conduction band minimum (VBM/CBM) close to the k = Z and Γ wavevectors, respectively. Using the HSE 06 hybrid functional, we obtained bandgaps of 0.93 eV and 0.82 eV for SnS and SnSe, which agree well with experimental measurements of 1.06 ± 0.15 eV and 0.86 eV [19,65]. The smallest direct gaps are considerably larger at 1.41 eV and 1.27 eV, respectively. The three rocksalt phases have a direct bandgap at k = L, as in the rocksalt-structured Pb chalcogenides [66]. We calculated a gap of 0.65 eV for both the equilibrium and compressed RS SnS structures and a gap of 0.6 eV for RS SnSe. All three RS structures were also predicted to have an indirect gap ∼50-130 meV smaller than the direct gaps, which is not seen along the calculated dispersion path, although the differences between the two are much smaller than in the Pnma phases.
Using semi-classical Boltzmann transport theory, the electrical transport properties can be determined from the spectral conductivity Σ and the nth-order moments of the generalised transport coefficients L n , defined as [54,67]: kj , ν kj , and τ kj are the energies, group velocities, and relaxation times of electrons with wavevector k and band index j, e is the elementary charge, and f 0 ( , T) is the Fermi-Dirac distribution given by: where F is the Fermi energy and k B is the Boltzmann constant. The three electrical properties in Equation (1) can then be obtained as: where the explicit temperature dependence of the transport properties and the L n has been omitted for brevity. The electron relaxation times τ kj in Equation (2) are computed as the sum of four scattering processes, namely acoustic deformation potential (ADP), piezoelectric (PIE), polar optical phonon (POP), and ionised impurity (IMP) scattering: 1 Full details of how the four scattering rates (inverse lifetimes) in Equation (8) are calculated can be found in [54]. The σ, S, and κ el are 3 × 3 tensors for which the three diagonal elements correspond to transport along the Cartesian x, y, and z directions. To more easily compare the three compounds, we also computed the scalar averages σ, S, and κ el : We first investigated the calculated power factor S 2 σ (PF), i.e., the numerator in Equation (1), and electronic thermal conductivity κ el of the Pnma and RS phases of SnSe as a function of temperature from 200-1000 K and at hole carrier concentrations n from 10 15 -10 20 cm −3 (see the Supplemental Materials).
For undoped SnSe, typical n are on the order of 2 × 10 17 cm −3 , and concentrations of up to 4 × 10 19 cm −3 can be achieved with hole doping [20,68]. Our initial survey of the calculated properties clearly shows that n on the order of 10 18 cm −3 and above are required to obtain a reasonable power factor, and, therefore, that doping is required to obtain good thermoelectric performance. At larger carrier concentrations, we predict the RS phase to support an up to 4 × larger PF than the Pnma phase, but this would be partially negated by a higher κ el , which appears in the denominator of Equation (1).
At carrier concentrations approaching the largest n = 10 20 cm −3 in our calculations, the Pnma phase shows a significant increase in the electronic thermal conductivity. Furthermore, at low n the κ el increases with temperature while at higher n, the reverse trend is observed. According to the Wiedemann-Franz law, in metals and degenerate semiconductors κ el is proportional to the electrical conductivity σ according to: where L is the Lorentz number and the explicit T dependence of κ el and σ has again been omitted. These observations therefore suggest that at around n 10 19 cm −3 Pnma SnSe becomes degenerated and the conductivity switches from a semiconductor-to metallic-like temperature dependence. This has been observed experimentally in hole-doped SnSe [20]. RS SnSe shows similar behaviour to the Pnma phase up to ∼500 K, whereas at higher temperatures, the κ el peaks at carrier concentrations around 10 18 -10 19 cm −3 and then falls as n is increased further. A similar analysis for the Pnma and equilibrium/compressed rocksalt phases of SnS yielded similar qualitative conclusions.
To compare the five structures, Figure 3 shows the three electrical properties σ, S, and κ el as a function of carrier concentration at a fixed T = 800 K and as a function of temperature at a fixed n = 2.15 × 10 19 cm −3 .
Comparing the PFs S 2 σ to the electrical conductivities and Seebeck coefficients shows that the dependence on n is primarily governed by the large increase in σ with carrier concentration. The two Pnma phases show a larger maximum Seebeck coefficient than the three RS structures, but the S peaks at a smaller n and falls at larger carrier concentrations, whereas the Seebeck coefficients of the RS phases peak at larger n. Together with their generally larger maximum σ, this results in the RS structures showing larger PFs than the Pnma phases at higher carrier concentrations. Interestingly, our calculations predict the Seebeck coefficients of the two Pnma phases to be negative at low hole carrier concentrations and to become positive for n 10 16 -10 17 cm −3 .
The κ el of the RS structures is predicted to be large for n up to ∼10 19 cm −3 , but falls sharply at higher carrier concentrations. With reference to the substantial increase in the conductivity of the RS phases with n, and given that we would expect an increasingly degenerate behaviour at higher carrier concentrations, this is in apparent violation of the Wiedemann-Franz law (Equation (12)). To investigate further, we compared the κ el obtained for the five structures using Equations (7) and (12) as a function of temperature and carrier concentration (see the Supplemental Materials). We found that the Wiedemann-Franz model predicted a 20-40% increase in the κ el of the Pnma structures and a reduction of up to 80% in the κ el of the RS phases at higher temperatures and carrier concentrations. We note, however, that the Lorentz number in Equation (12) can vary between systems and with temperature, and calculations for hole-doped SnSe, based on a multi-band model, predict a value ∼20% smaller at 300 K and ∼25% smaller at ∼800 K [69]. This would account for a large part of the increased κ el of the Pnma phases predicted with this method.
In spite of these issues, as we discuss in Section 3.4 below, we found that the two different models for κ el ultimately had a relatively small impact on the predicted maximum thermoelectric figures of merit calculated using Equation (1).
With a fixed n = 2.15 × 10 19 cm −3 , the conductivities of the Pnma phases show a metallic-like reduction with temperature as noted previously. The σ of the RS phases show a similar fall at low-to-moderate T, but increase again at high temperatures. The result is that whereas the two RS SnS models are predicted to have significantly lower σ than Pnma SnS at 200 K, they have ∼2 × higher conductivities at 1000 K. At this carrier concentration, the Seebeck coefficients of the Pnma phases increase with temperature and converge toward very similar values at T = 1000 K. The S of the three RS phases, on the other hand, all peak around 600 K and fall sharply above this temperature, such that they are significantly smaller than those of the Pnma phases at 1000 K. This sharp reduction in the Seebeck coefficient at higher temperature results in the RS phases showing smaller S than the corresponding Pnma phases at 1000 K and is clearly reflected in a reduction in the power factors. Finally, a comparison of the κ el as a function of temperature again highlights the much larger electrical thermal conductivities of the RS structures compared to the Pnma phases. Taken together, these results indicate that RS SnSe would require a larger carrier concentration (heavier doping) for its high-temperature thermoelectric performance to be competitive with the Pnma phase. Whereas the electrical properties suggest RS SnSe may show better performance at lower temperatures, this depends on any differences in the lattice thermal conductivity κ latt , which we address in Section 3.3 below. Similar qualitative conclusions can be made for SnS.
Whereas the cubic symmetry of the rocksalt structure means that the three diagonal components of the conductivity, Seebeck, and electrical thermal conductivity tensors are equivalent and equal to the average, the anisotropic bonding in the orthorhombic Pnma structure results in a significant directional dependence of all three quantities. Figure 4 shows the xx, yy, and zz components of the power factors and κ el of Pnma SnSe as a function of n and T, as in Figure 3, together with the isotropic averages for the Pnma and RS phases (cf. Equations (9)-(11)). In the Pnma structure, the xx, yy, and zz components of the tensors correspond to transport along the crystallographic a, b, and c directions, respectively (cf. Figure 1). Similar analyses comparing the directional anisotropy in the σ, S, PF, and κ el of the Pnma phases of both SnS and SnSe to the average values of the Pnma and RS phases shown in Figure 3 are provided as Supplemental Materials. We found that the Pnma structures have distinct "easy" and "hard" axes for transport, corresponding respectively to the b axes along which the bonding is strongest and the layered a axes. This leads to a large anisotropy in the electrical conductivity both as a function of carrier concentration and temperature. On the other hand, the Seebeck coefficients are strongly anisotropic at small n, but converge at carrier concentrations above ∼10 18 cm −3 . The anisotropy in the σ with temperature and the similar S at larger n are both consistent with experimental measurements [20]. The directional dependencies are reflected in the power factors, although since appreciable PFs are only obtained at larger carrier concentrations, the PFs largely follow the anisotropy in the conductivity. However, the electronic thermal conductivity also reflects the anisotropy in σ, as would be expected from Equation (12), and this may serve to limit the directional variation in the thermoelectric figures of merit.
Finally, it is of interest to examine the relative impact of the acoustic deformation potential, piezoelectric, polar optic phonon, and ionised impurity scattering mechanisms on the electron relaxation times τ kj (cf. Equation (8)). Figure 5 shows the calculated scattering rates (i.e., the inverse lifetimes τ −1 kj ) as a function of the electronic state energy up to 0.4 eV below the Fermi energy = F for n = 2.15 × 10 19 cm −3 and T = 800 K. For the Pnma phase, our calculations predict that POP scattering is dominant, with contributions from ADP and IMP scattering in states close to the Fermi energy. For RS SnSe, on the other hand, POP scattering dominates close to F , whereas at lower energies, the POP and ADP scattering have similar rates, and IMP scattering does not seem to play a significant role. Piezoelectric scattering is not relevant to either of the Pnma or RS structures because they are both centrosymmetric and the piezoelectric moduli therefore vanish. A similar analysis of the calculated scattering rates of the three SnS structures (see the Supplemental Materials) leads to similar qualitative conclusions.   (11). The isotropic averages of the rocksalt phase, for which the three diagonal components are equal, are also shown for comparison. As in Figure 3, both properties are shown as a function of carrier concentration n for a fixed T = 800 K (a,c) and as a function of temperature for a fixed n = 2.15 × 10 19 cm −3 (b,d). We noted during our calculations that the high-frequency and static dielectric constants of the three rocksalt phases computed with PBEsol + D3 were considerably larger than those of the Pnma phases, which we tentatively attributed to the near-metallic electronic structures obtained for the equilibrium RS SnS and SnSe structures with this method [45]. Since the ε ∞ and ε s are used to calculate the POP and IMP scattering rates, we also computed transport properties with the dielectric properties of these structures calculated using HSE 06 (see the Supplemental Materials). In general, we found that the HSE 06 dielectric constants led to an increase in the conductivities and electronic thermal conductivities while leaving the Seebeck coefficients largely unchanged. While at some carrier concentrations and temperatures, the calculated σ was increased by a factor of two, κ el was increased by a similar factor, so we would not expect a significant overall effect on the calculated figures of merit (cf. Equation (1)). We also note that using the different dielectric constants did not affect the reduction in κ el at large n predicted using Equation (7).
The transport properties predicted from these calculations are generally a reasonable match to the experiments. As one of the current flagship chalcogenide thermoelectrics, the transport properties of Pnma SnSe have been studied extensively. The doped samples in [20] achieved carrier concentrations around 4 × 10 19 cm −3 . This study reported Seebeck coefficients of 160 µV K −1 and 300 µV K −1 at 300 K and 773 K, with little variation along the three crystallographic axes, and we obtained similar values of ∼150 µV K −1 and 250 µV K −1 at 300 K and 780 K with n = 4.64 × 10 19 cm −3 . The measured conductivities of 1486 S cm −1 and 148 S cm −1 along the bonding direction at the two temperatures are on the same order of magnitude as the 2178/3006 S cm −1 and 641/918 S cm −1 calculated for n = 3.16 × 10 19 /4.64 × 10 19 cm −3 . With the higher carrier concentration, we calculated PFs of 56.8 µW cm −1 K −2 and 53.9 µW cm −1 K −2 at 300 K and 780 K. The former compares reasonably well to the measured value of 40 µW cm −1 K −2 at 300 K, whereas the latter is nearly 4 × larger than the measured 14 µW cm −1 K −2 at 773 K. The measurement in [68] on polycrystalline samples with n 10 19 cm −3 yielded conductivities in the range of 140-180 S cm −1 and 120-130 S cm −1 at 423 K and 783 K, a Seebeck coefficient of 342 µV K −1 at 673 K, and PFs on the order of 10 µW cm −1 K −2 from 473-783 K. At the same carrier concentration, we obtained comparable averaged σ of 287 S cm −1 and 122 S cm −1 at 420 K and 780 K, a similar S of 360 µV K −1 at 680 K, and power factors between 17.2 µW cm −1 K −2 and 24.7 µW cm −1 K −2 . SnS has been comparatively less well studied, but the experiments on polycrystalline samples in [31] with n 2 × 10 18 cm −3 suggest a maximum Seebeck coefficient of 500 µV K −1 at 600 K, at which temperature the conductivity is ∼3 S cm −1 and the power factor is 0.75 µW cm −1 K −2 . With n = 2.15 × 10 18 cm −3 , we predicted an averaged S and σ of 433 µV K −1 and 33 S cm −1 , giving a PF of 6.2 µW cm −1 K −2 .
We therefore conclude that the protocol in [54] gives reasonable predictions of the Seebeck coefficients in these systems, but tends to overestimate the electrical conductivity. The fact that the S are well reproduced suggests the calculated electronic structures are reasonable, so we might attribute the overestimation of the σ either to approximations in the models used to calculate the electronic relaxation times in Equation (2) or to errors in calculating the requisite material properties. However, it is also worth noting that our calculations are based on perfect bulk crystals, whereas real materials are likely to contain defects such as grain boundaries that could limit the electrical transport. Our model also does not consider the effect of temperature on the structure through, e.g., thermal expansion at finite temperature. This could in principle be addressed by using the quasiharmonic approximation to model the thermal expansion, as in our previous study [45], or, more simply, by fixing the cell volumes to experimentally measured values at appropriate temperatures. Using either method would, however, require calculations on multiple structures, which would significantly increase the computational workload. In any case, we would hope that any fundamental issues in the transport calculations would have a similar effect on all five systems and, therefore, that the results would remain comparable, which was the main focus of this study.

Lattice Thermal Conductivity
To compute the κ latt in Equation (1), we used the single-mode relaxation time approximation (RTA) model [60]. The macroscopic thermal conductivity κ latt was computed as the sum of contributions from individual phonon modes with wavevector q and band index j as: V is the unit cell volume, and N is the number of wavevectors included in the summation, which is equivalent to the number of unit cells in the crystal. C qj are the modal heat capacities calculated according to: (14) ν qj are the mode group velocities, which are the derivatives of the frequencies ω qj with respect to the wavevector: (15) τ qj are the phonon lifetimes, calculated as the inverse of the phonon linewidths Γ qj : The calculation of the linewidths requires the harmonic phonon frequencies ω qj and eigenvectors W qj together with the third-order force constants, and the method is documented in detail in [60]. We previously demonstrated that the RTA model predicts lattice thermal conductivities for Pnma SnS and SnSe in reasonable agreement with experiments [61], albeit with a tendency to overestimate as observed with the calculated electrical conductivities in Section 3.2. Figure 6 compares the temperature dependence of the κ latt of the five models examined in this study, and values at T = 800 K are listed in Table 2. As for the electrical transport properties, we considered the principal xx, yy, and zz components of the κ latt tensors of the two Pnma phases, corresponding to heat transport along the a, b, and c axes, respectively (cf. Figure 1), together with the scalar average calculated as: where the explicit temperature dependence has been omitted for brevity. For the rocksalt phases, we only considered the average, as in these systems, the three diagonal components are equivalent by symmetry and equal to the average. Table 2. Calculated thermal conductivities of the five structures examined in this work at T = 800 K. Each row lists the three diagonal components κ xx , κ yy , and κ zz of the κ latt tensor together with the diagonal average κ ave = 1 3 κ xx + κ yy + κ zz and its decomposition into harmonic and lifetime components, κ/τ CRTA ave and τ CRTA , according to Equation (18). The data for Pnma SnS and SnSe is from [61] with the κ xx , κ yy , and κ zz relabelled to be consistent with the orientation of the unit cells in this work.  For the Pnma phases, the principal κ xx , κ yy , and κ zz components of the κ latt tensor, corresponding to transport along the a, b, and c axes, respectively, are shown together with the diagonal average κ ave = 1 3 κ xx + κ yy + κ zz . For the rocksalt phases, the three diagonal components and the average are equivalent by symmetry, so we only show κ ave . The data for Pnma SnS and SnSe is from [61], but the three principal components have been relabelled to match the orientation of the unit cells in this study.
The thermal conductivity shows a steep temperature dependence, which is particularly apparent for the rocksalt phases. As with the electrical transport properties, the Pnma phases show the largest κ latt along the strongly bonded b axes with much lower thermal conductivities along the layered a directions. In Pnma SnS and SnSe, the ratios of κ yy and κ xx are 2.35 and 2.3, respectively. Our calculations predict the thermal conductivity of rocksalt SnSe to be ∼2-5 × larger than along any of the axes in the corresponding Pnma structure. On the other hand, the κ latt of the equilibrium structure of RS SnS is predicted to be ∼25% smaller than the average for Pnma SnS. The thermal conductivity in the rocksalt-structured Pb chalcogenides is strongly volume dependent [66], as the bond distances are straightforwardly related to the lattice constant and cell volume such that compression and expansion of the lattice strengthens and weakens the chemical bonding, respectively. The same phenomenon is evident here in the ∼5 × larger predicted κ latt of the compressed rocksalt SnS structure compared to the equilibrium structure. Noting that the equilibrium RS SnS structure is dynamically unstable, the higher predicted κ latt of the compressed RS SnS and RS SnSe structures compared to the corresponding Pnma phases would serve to further negate the higher power factors predicted by the transport calculations in Section 3.2.
We note in passing that our results for RS SnSe are inconsistent with the modelling study in [46], which modelled the thermal conductivity using a similar approach to the present study but predicted κ latt to be smaller than the average for the Pnma phase. The study noted that a small displacement of the Se atoms along the <111> direction (i.e., a rhombohedral distortion) was applied to remove imaginary modes in the phonon dispersion, and given that the imaginary mode in equilibrium RS SnS appears to result in a very low lattice thermal conductivity this may explain the discrepancy. Based on a previous comparison of how different DFT functionals predict the quasi-harmonic lattice dynamics of the Pb chalcogenides, it is possible that the imaginary modes are an artefact of the tendency of the PBE functional to overestimate the unit cell volume [70], but the lattice constant of 6.06 A is similar to our predicted value of 5.99 Å. It is possible that RS SnSe possesses an incipient structural distortion that manifests under thermal expansion at elevated temperature-such a phenomenon has been proposed to occur in PbTe, but is heavily disputed [71][72][73][74]. However, we also note that the calculations in [46] used a smaller plane wave cutoff and a smaller supercell for calculating the second-order force constants of RS SnSe, which could also lead to differences in the calculated phonon dispersion and thermal conductivity.
Conceptually, it is useful to discuss the differences in the thermal transport between the five systems in terms of differences in the group velocities ν qj and the lifetimes τ qj . This distinction can be made using the constant relaxation time approximation (CRTA) model developed in [61,75]. In this model, κ latt are written as the product of a harmonic term and a weighted average lifetime τ CRTA : We note that both the harmonic and lifetime terms are implicitly temperature dependent, the latter through the temperature dependence of the C qj . Comparing τ CRTA and the harmonic term in the summand then allows the difference in κ latt to be attributed quantitatively to differences in the harmonic term (mainly the group velocities) and the phonon lifetimes. Figure 7 compares the averaged κ latt as a function of temperature for the five models to the corresponding averaged harmonic function κ/τ CRTA and lifetime τ CRTA , and values of the harmonic and lifetime terms at T = 800 K are also included in Table 2 for comparison.
In general, the harmonic κ/τ CRTA function saturates rapidly with temperature as C qj tend to the Dulong-Petit limit, and the sharp decrease in τ CRTA with temperature primarily determines the temperature dependence of κ latt . As discussed in our previous study [61], we found that the lower predicted κ latt of Pnma SnSe results from a balance of a smaller harmonic term and a longer lifetime. This CRTA analysis also provides insight into the anomalously small κ latt computed for the equilibrium RS SnS structure-for this structure, κ/τ CRTA is predicted to be larger than for the compressed phase, while τ CRTA is predicted to be significantly smaller than for all the other systems including both Pnma phases. Since κ/τ CRTA primarily reflects the differences in group velocities, which are given by the slope of the phonon dispersion, we attributed the former anomaly to the distortion of the band structure due to the presence of the imaginary mode (cf. Figure 2). The dynamical instability may also have the effect of increasing the third-order force constants, which could explain the extremely short τ CRTA . On the other hand, both the compressed RS SnS and RS SnSe structures show larger κ/τ CRTA and shorter τ CRTA than the corresponding Pnma phases, highlighting a significant difference in the lattice dynamics of the two structure types.

Thermoelectric Figure of Merit
We now combine the calculated electrical properties from Section 3.2 with the lattice thermal conductivities from Section 3.3 to calculate the thermoelectric figure of merit ZT as a function of carrier concentration and temperature using Equation (1) (Figure 8).
As observed in Section 3.2, these plots clearly show that achieving high figures of merit for all five systems requires high carrier concentrations to optimise the electrical transport. The analysis also shows that ZT is generally maximised at higher temperatures, which can be attributed mainly to the sharp reduction in κ latt with temperature highlighted in Section 3.3. For the Pnma phases, the optimum ZT appears to be obtained for carrier concentrations in the region of 10 19 cm −3 and close to the maximum T = 1000 K examined in our calculations. This optimum n balances the Seebeck coefficient, conductivity, and electrical thermal conductivity, while the optimum T is because at these carrier concentrations, higher temperatures maximise S, compensating the reduced σ, while also minimising κ latt . For the RS phases, the optimum ZT scores were predicted to occur at the maximum carrier concentration of 10 20 cm −3 we investigated, but at lower temperatures. In these systems, unlike the Pnma phases, the Seebeck coefficient remains high at larger n, but κ el increases significantly at higher temperatures. Since the RS phases generally have higher lattice thermal conductivities than the Pnma phases, the optimum temperature is determined mainly as a balance between the two contributions to the thermal conductivity.
We note in passing that combining the predicted electrical transport properties with the calculated lattice thermal conductivity implicitly assumes that κ latt is independent of n, i.e., that the required doping levels do not introduce impurities at a level that would materially disrupt the heat transport through phonons.
From each of the plots in Figure 8, we extracted a maximum ZT, which we summarise together with the corresponding carrier concentration/temperature and the properties from Equation (1) in Table 3.
For the two Pnma phases, we predict maximum averaged ZT values of 1.75 and 2.81, respectively. The latter compares favourably to the ZT of 3.1 for polycrystalline SnSe reported in [68], whereas our prediction for SnS is considerably higher than the value of 1.1 reported for polycrystalline SnS [31,33]. For both systems we predict the optimum ZT is obtained for n = 4.64 × 10 19 cm −3 and T = 1000 K. The carrier concentrations are comparable to those achieved in the studies on SnSe in [20,68], but are larger than the 2 × 10 18 cm −3 and 2 × 10 19 cm −3 obtained for SnS in [31,32]. Studies on both materials also showed the best ZT at elevated temperature, which is again consistent with our predictions. As for the constituent properties, there is significant anisotropy in the figure of merit, with the ZT along the layering a and strongly-bonded b axes calculated to be 1.96 and 1.24 for SnS and 3.15 and 1.85 for SnSe. This is consistent with the single-crystal measurements in [19,20], although the maximum reported ZT scores of ∼2 and 2.6 along the b axis are smaller than the maximum predicted in our calculations. We note however that 1000 K is above the phase transition temperatures of both Pnma SnS and SnSe, which may lead to a discrepancy between our predictions and measurements. It is also of note that the optimal figure of merit in polycrystalline SnSe reported in [68] was obtained by purification to remove tin oxides, and this or similar issues could explain why the experimental studies on SnS have so far reported considerably lower ZT scores than our predicted maximum [31][32][33].  Table 3. Predicted maximum thermoelectric figures of merit ZT for the five systems examined in this work, extracted from Figure 8. For each system, we list the carrier concentration n and temperature T at which the maximum ZT was obtained together with the properties from Equation (1), namely the electrical conductivity σ, the Seebeck coefficient S, the power factor S 2 σ (PF), the electronic and lattice contributions to the thermal conductivity κ el /κ latt , and the total thermal conductivity κ tot . For the equilibrium and compressed rocksalt SnS structures, our calculations predict maximum ZT scores of 2.99 and 1.48, respectively. These are both obtained for our maximum calculated carrier concentration of 10 20 cm −3 , but at lower temperatures of 720 K and 800 K compared to Pnma SnS. Both RS structures are predicted to have 3-4 × higher electrical conductivities than the Pnma phase and comparable Seebeck coefficients, resulting in much larger power factors. This, plus the low predicted κ latt , leads to the large predicted ZT of the equilibrium RS SnS structure. On the other hand, the high κ latt of compressed RS SnS, when combined with its high κ el , results in a total thermal conductivity ∼4 × higher than the Pnma phase and a smaller ZT. For RS SnSe, we predict a maximum ZT of 2.6 with n = 10 20 cm −3 at T = 800 K. In this case, σ is around 3 × higher than Pnma SnSe and the Seebeck coefficient is again comparable, but the higher κ el and κ latt results in an overall lower figure of merit. Since 10 20 cm −3 is the largest of the n considered in the transport calculations in Section 3.2, it is possible that the ZT of the RS phases could be enhanced further with higher doping levels. However, since experiments on the Pnma phases have so far only reported carrier concentrations on the order of 10 19 cm −3 , such high doping levels may not be achievable in practice.
Given the deviation in the predicted κ el from the Wiedemann-Franz law, we performed a similar analysis to determine the maximum ZT with the electronic thermal conductivity computed using Equation (12) (see the Supplemental Materials). As expected, given that this predicted a larger κ el for the Pnma phases, the average ZT scores of Pnma SnS and SnSe decrease to 1.5 and 2.13, respectively. However, despite the marked reduction in the κ el of the RS phases at some carrier concentrations, the maximum figures of merit are largely unchanged. The ZT of the equilibrium and compressed RS SnS structures decrease and increase to 2.81 and 1.56, respectively, while the maximum predicted ZT of RS SnSe decreases to 2.33. Given this, plus the fact that the ZT of Pnma SnSe obtained using the κ el from Equation (7) is in better agreement with experiments, we consider the results in Table 3 to be robust to this anomaly.

Discussion
In this study, we have compared the thermoelectric performance of the well-characterised Pnma phases of SnS and SnSe to the less well studied rocksalt structures.
By combining the protocol for modelling electronic transport properties in [54], including approximate models for the electron scattering rates, with the relaxation time approximation for the lattice thermal conductivity in [60], we were able to make a fully ab initio prediction of the thermoelectric figures of merit ZT. Comparison with the large body of experimental studies on SnS and SnSe suggested that these calculations gave very good estimates of the Seebeck coefficient, but tended to overestimate the electrical conductivity and the thermal conductivity. While this may occur for a variety of reasons, the cancellation of errors results in a good overall prediction of ZT, and we therefore consider this approach to be suitable for use in a predictive capacity.
Our calculations predict that, provided that carrier concentrations on the order of 10 20 cm −3 can be obtained by doping, the rocksalt phases may show larger electrical conductivities and comparable Seebeck coefficients to the orthorhombic phases, but the resulting enhanced power factors are likely to be negated by higher thermal conductivities. Since the phonon contributions to the thermal conductivity are generally large in these systems, there may be scope to further improve ZT by using structural modifications to suppress the heat transport through the lattice [1,9]. Based on our results, we would conclude that rocksalt SnSe is unlikely to show superior thermoelectric performance to the Pnma phase, and whereas some improvement on the relatively poorer ZT of Pnma SnS might be possible in an alternative RS structure, it would require stabilising it in a geometry close to the equilibrium structure. Given that there is experimental evidence to indicate that RS SnS can be grown epitaxially [37,38], this could plausibly be achieved in a thin film device using a suitable contact material. It is also possible that the RS phase could be stabilised with suitable doping or atomic substitution-e.g., partial replacement of Sn with Pb, since the analogous PbS and PbSe adopt the RS structure-but we are not aware of any experimental data that would support this suggestion. On the other hand, it is worth noting that maintaining an RS structure with an incipient dynamical instability under the high temperatures in working thermoelectric generators may be difficult to achieve in practice.
However, both SnS and SnSe have been prepared in π-cubic phases that are closely related to the rocksalt structure [42,43,45], but with distortions that have been predicted to reduce the lattice thermal conductivity [76], and our results therefore motivate an in-depth investigation to establish the thermoelectric performance of these systems.