Thermodynamics, Electronic Structure, and Vibrational Properties of Sn n (S1-x Se x ) m Solid Solutions for Energy Applications.

The tin sulfides and selenides have a range of applications spanning photovoltaics and thermoelectrics to photocatalysts and photodetectors. However, significant challenges remain to widespread use, including electrical and chemical incompatibilities between SnS and device contact materials and the environmental toxicity of selenium. Solid solutions of isostructural sulfide and selenide phases could provide scope for optimizing physical properties against sustainability requirements, but this has not been comprehensively explored. This work presents a detailed modeling study of the Pnma and rocksalt Sn(S1-x Se x ), Sn(S1-x Se x )2, and Sn2(S1-x Se x )3 solid solutions. All four show an energetically favorable and homogenous mixing at all compositions, but rocksalt Sn(S1-x Se x ) and Sn2(S1-x Se x )3 are predicted to be metastable and accessible only under certain synthesis conditions. Alloying leads to a predictable variation of the bandgap, density of states, and optical properties with composition, allowing SnS2 to be "tuned down" to the ideal Shockley-Queisser bandgap of 1.34 eV. The impact of forming the solid solutions on the lattice dynamics is also investigated, providing insight into the enhanced performance of Sn(S1-x Se x ) solid solutions for thermoelectric applications. These results demonstrate that alloying affords facile and precise control over the electronic, optical, and vibrational properties, allowing material performance for optoelectronic applications to be optimized alongside a variety of practical considerations.

Introduction

The current drive toward energy solutions to reduce our dependence on fossil fuels has prompted wide-ranging research to identify new, improved functional materials for cleaner energy generation and storage. Photovoltaics (PV) is an established technology for primary energy generation, but identifying cost-effective and environmentally benign materials for efficient light-to-electricity conversion remains an active research area.[1] Thermoelectric materials are another highly promising development area,[2,3] allowing waste heat to be recovered as electricity from, for example, industrial processes and combustion engines, and recent breakthroughs are steadily pushing toward higher-efficiency thermoelectric generators.[3] Hydrogen has long been seen as a clean alternative fuel for combustion engines, and identifying efficient photocatalysts to generate H2 by splitting water remains an important research goal.[4,5]

The tin sulfides and selenides are unique in having broad applications across the energy materials domain. Pnma SnS has a large optical absorption coefficient and a bandgap well matched to the solar spectrum and as such has been widely studied as a PV absorber.[6−8] The structurally analogous Pnma SnSe was recently shown to have a record-breaking thermoelectric figure of merit,[9] with hole doping enabling efficient energy harvesting over a wide range of operating temperatures.[10] SnS2 and SnSe2 have potential applications as photocatalysts for water splitting[11] and as photodetectors,[12] respectively, with the added advantage that the pseudo-two-dimensional (2D) layered structure can be formed into a variety of high-surface-area nanomaterials.[13] The tin selenides have also shown potential for PV applications[11,14−16] and as components in lithium-ion batteries[11,14] and supercapacitors.[17]

These materials are not without their drawbacks, however. The environmental toxicity of selenium limits the widespread use of tin selenides, but while Pnma SnS has been studied as an alternative thermoelectric,[18−20] the greatly reduced performance indicates that the complete substitution of selenium by sulfur cannot sustain the requisite material properties. Another major issue is that although SnS has been intensively researched as a cost-effective PV absorber, SnS-based devices have so far fallen short of the efficiencies obtained with alternative materials,[6,7,21] which may be due to several factors including the propensity of SnS to form phase impurities, for example, by excess oxidation of Sn during synthesis, a mismatch in energy levels (band offsets) between SnS and the contact materials in devices, and chemical incompatibility between SnS and the contact materials leading to the formation of impurities through detrimental reactions.[21,22]

Due to the ability of Sn to adopt both Sn(II) and Sn(IV) oxidation states[23] and the structural flexibility of the Sn lone pair, the tin sulfides and selenides form several isostructural stable and metastable phases. Five monosulfide/monoselenide phases have been identified or proposed, viz., orthorhombic Pnma and Cmcm(24−26) and cubic rocksalt,[27] zincblende,[28] and P213 (“π-cubic”),[29,30] alongside the di- and sesqui-sulfide/selenide phases SnS2/SnSe2 and Sn2S3/Sn2Se3.[31−33]

The similar sulfide and selenide phase spaces and the similar covalent radii of S and Se (100 vs 115 pm)[34] make it a natural step to consider Sn(S1–Se) solid solutions as a route to fine tuning the material properties and balancing the performance against sustainability requirements. A small number of studies have explored this in the context of thermoelectric materials[35−37] and have demonstrated that substituting a small amount of the Se in Pnma SnSe can enhance the performance compared to both the SnS and SnSe endpoints. However, there has been very little research into solid solutions of the other phases nor more generally into using this strategy to adapt the material properties for other applications.

In this work, we expand our recent theoretical studies of the SnS system[23,38,39] and make use of high-throughput first-principles computational modeling and statistical thermodynamics to explore the stability and properties of Pnma and rocksalt Sn(S1–Se), Sn(S1–Se)2, and Sn2(S1–Se)3 solid solutions. Our results show that forming the solid solutions is facile and enables the structural, electronic, and vibrational properties of the four systems to be systematically tuned between the endpoints, providing a route to optimizing the properties for energy applications. We also address the practical challenges associated with accurate, high-throughput calculations on large collections of structures, providing a general direction for future studies on other solid solutions.

Experimental Section

The starting point for our calculations is a set of optimized supercells built from the primitive cells of Pnma and rocksalt SnS, SnS2, and Sn2S3 and the corresponding isostructural selenides. The structures of the sulfide endpoints[23,39] are shown in Figure .

Figure 1

Structures of Pnma and rocksalt SnS, SnS2, and Sn2S3 with the sulfur and tin atoms shown in yellow and blue, respectively. The images were prepared using the VESTA software.[55]

Solid-solution models were built using the open-source Transformer library[40] by enumerating all of the symmetry-inequivalent structures formed by successively substituting the chalcogen atoms in the parent supercells together with the associated degeneracies. We found that taking supercells with 20–32 atoms yielded a manageable number of unique configurations for high-throughput density functional theory (DFT) calculations (Table ).

Table 1

Summary of the Solid Solutions Studied in This Work with the Supercell Sizes, Number of Configurations, k-Point Sampling, and Phonon Supercells Used for Finite-Displacement Lattice-Dynamics Calculations

system	supercell (# atoms)	# structures (unique)	k-point samplinga	phonon supercell (# atoms)
Pnma Sn(S1–xSex)	2 × 1 × 2 (32)	65 536 (2446)	4 × 4 × 4	2 × 1 × 2 (128)
rocksalt Sn(S1–xSex)	2 × 4 × 2 (32)	65 536 (652)	4 × 2 × 4	2 × 1 × 2 (128)
Sn(S1–xSex)2	2 × 2 × 2 (24)	65 536 (1056)	4 × 4 × 2	2 × 2 × 2 (192)
Sn2(S1–xSex)3	1 × 1 × 1 (20)	4096 (1072)	4 × 8 × 3	2 × 2 × 2 (160)

The k-point sampling is specified as the number of subdivisions along each reciprocal lattice vector in a Monkhorst–Pack k-point grid.[56]

Each unique structure was fully volume relaxed. To speed up this process, structures with higher S and Se content were optimized starting from the lattice parameters of the respective sulfide and selenide endpoints. The total energies after optimization were combined with the degeneracies to construct the thermodynamic partition functions and to obtain the free energies of mixing and statistical averages of material properties.

The DFT calculations were performed within the pseudopotential plane-wave DFT formalism implemented in the Vienna Ab initio Simulation Package (VASP) code.[41] Based on previous work on the tin sulfide system,[23] we used the PBEsol generalized-gradient approximation (GGA) exchange-correlation functional[42] with the DFT-D3 dispersion correction[43] (i.e., PBEsol-D3). The basis set was defined with a 600 eV kinetic-energy cutoff, and the k-point sampling listed in Table was used to integrate the electronic Brillouin zone. The core electrons were treated with projector augmented-wave (PAW) pseudopotentials including the Sn 5s, 5p and 4d, S 3s and 3p, and Se 4s and 4p electrons in the valence shell.[44,45] Tolerances of 10–8 eV and 10–2 eV Å–1 were applied to the total energy and forces during wavefunction minimization and geometry optimization, respectively. The precision of the charge-density grids was set automatically to avoid aliasing errors, the PAW projection was performed in reciprocal space, and nonspherical contributions to the gradient corrections inside the PAW spheres were accounted for.

As detailed in the Results and Discussion, after benchmarking the accuracy and efficiency of several techniques for modeling the electronic structure, we opted to perform electronic-structure calculations using the strongly constrained and appropriately normed (SCAN) meta-GGA functional[46] on the full set of structures in the four solid-solution models and to assess the accuracy of a subset using “single-shot” (non-self-consistent)[47] hybrid calculations with the HSE 06 functional.[48] For both, the plane-wave cutoff was reduced to 400 eV, which we found sped up the calculations without a significant loss of accuracy. The origin of the k-point meshes was shifted to = (0, 0, 0) (Γ) to include a larger number of irreducible k-points in the Brillouin-zone sampling, and the integration was performed using the Blöchl-corrected tetrahedron method.[49] Optical properties were evaluated from frequency-dependent dielectric functions obtained within the independent-particle approximation,[50] with equal numbers of occupied and virtual states being included in the calculations to converge the sum over empty states.

Lattice-dynamics calculations were performed using the Phonopy package[51,52] with the VASP code as the force calculator. The supercell expansions used to evaluate the force-constant matrices with the Parlinski–Li–Kawazoe method[53] are listed in Table . For the single-point force calculations, the k-point sampling was reduced according to the chosen supercell expansions, and the same tolerance of 10–8 eV on the total energy used for the other calculations was applied during the wavefunction optimization. Phonon calculations on the sulfide endpoints for comparison were taken from our previous work,[23] and analogous calculations were performed on the selenide endpoints. Phonon density of states (phonon DoS) curves were computed by interpolating the phonon frequencies onto uniform Γ-centered q-point meshes with between 24 × 24 × 24 and 48 × 48 × 48 subdivisions and applying a Gaussian broadening with σ = 0.15 THz (approx. 5 cm–1). As discussed below, it was not feasible to perform lattice-dynamics calculations on all of the structures, so selected compositions were sampled using a structural-fingerprinting technique to assess the structural similarity and thereby identify structurally diverse chalcogen arrangements.

Results and Discussion

Mixing Energies

To assess the thermodynamic stability of the solid solutions, we calculated for each system the Gibbs free energy of mixing Gmix(xSe, T) as a function of Se fraction xSe and temperature T.

For a given composition, the Gibbs energy G(xSe, T) is calculated from the thermodynamic partition function Z(xSe, T)where the sum runs over N symmetry-inequivalent structures with total energy E, degeneracy g and kB is Boltzmann’s constant. Since the structures are fully volume relaxed (i.e., p = 0), the enthalpy = U + pV is equal to the internal energy U, and G (xSe, T) can be calculated via the bridge relationFurther details can be found in ref (54). We note that the effect of finite pressure could be accounted for by adding a correction pV to each E, and that the total energies do not include vibrational contributions to the free energy[23] as to do so would be infeasible for the ∼5000 configurations in this study.

Gmix(xSe, T) is then calculated by comparing the free energy of the mixed phase to the phase-separated endpoints according towhere G(xSe = 0) and G(xSe = 1) are the Gibbs energies of the sulfide and selenide endmembers, respectively, which, since g is unity, are simply equal to the calculated DFT total energies E.

Figure shows the calculated mixing energies for each of the four systems as a function of composition and formation temperature. We predict favorable mixing at all compositions, with Gmix for the 50/50 mixtures ranging from −1.63 kJ mol–1 atom–1 for rocksalt Sn(S1–Se) to −2.24 kJ mol–1 atom–1 for Sn(S1–Se)2 at a formation temperature of 900 K. Comparing the enthalpy and entropy terms (see Figure S1.1) shows a mixing enthalpy Hmix of ∼0.49–0.70 kJ mol–1 atom–1, which we ascribe to the strain induced by the small size difference between S and Se and a larger TSmix term of 2.18–2.93 kJ mol–1 atom–1 from the configurational entropy.

Figure 2

Per-atom mixing free energies Gmix as a function of Se fraction xSe for Pnma (a) and rocksalt (b), Sn(S1–Se), Sn(S1–Se)2 (c), and Sn2(S1–Se)3 (d) solid solutions. The line colors denote formation temperatures from 300 K (blue) to 1100 K (orange).

A recent study employed an alternative first-principles cluster expansion technique to examine the stability of Pnma Sn(S1–Se) solid solutions.[57] This study reported formation energies ranging from 0.24 to 0.48 kJ mol–1 atom–1 for structures with the 50/50 composition and predicted homogenous mixing with calculated mixing energies of −0.96, −1.93, and −2.89 kJ mol–1 atom–1 at 300, 450, and 600 K. Our calculations predict an averaged mixing enthalpy of 0.48–0.59 kJ mol–1 atom–1 between 300 and 1100 K and mixing energies between −0.25 and −0.99 kJ mol–1 atom–1 from 300 to 600 K, which we consider sufficiently similar given the differences in methodology and computational setup.

To check the convergence of our results with respect to the supercell size, we performed additional calculations on the Pnma Sn(S1–Se) and Sn(S1–Se)2 solid solutions using smaller 1 × 1 × 2 and 2 × 2 × 1 expansions with 12 and 16 atoms, respectively. We obtained differences of up to 0.04 and 0.26 kJ mol–1 atom–1 in the values of Hmix and Gmix calculated at 900 K, which suggests that our results are reasonably well converged with respect to the cell size.

In our recent computational study of the energetic and dynamical stability of the tin sulfides,[23] we established that Pnma SnS was energetically more stable than the rocksalt phase, and that the sesquisulfide Sn2S3 was slightly above the energetic convex hull but ultimately stabilized with respect to decomposition by its larger vibrational entropy. The calculated energies of rocksalt SnS and SnSe place them 10.45 and 2.6 kJ mol–1 per SnS formula unit (F.U.), respectively, above the corresponding Pnma phases on the convex hull (see Figure S1.2). The difference falls almost linearly with increasing Se content, indicating that the rocksalt structure becomes more energetically accessible toward the selenide endmember but remains metastable with respect to the competing Pnma phase.

Sn2S3 is predicted based on lattice energy to be unstable to decomposition into SnS2 and Pnma SnS but only by 0.03 kJ mol–1 per F.U., which is in the range of differences in vibrational zero point energy.[23] On the other hand, the hypothetical selenide Sn2Se3 is somewhat less stable, with disproportion into SnSe and SnSe2 favored by a much larger 2.4 kJ mol–1 per F.U. (Figure S1.2). The disproportion of the mixed composition Sn2(S0.5Se0.5)3 into Sn(S0.5Se0.5) and Sn(S0.5Se0.5)2 is favored by 0.9 kJ mol–1 at 300 K, and higher temperatures further destabilize the mixed phase, with the energy difference rising to ∼1.8 kJ mol–1 at 900 K. This suggests that it may be possible to prepare Sn2(S1–Se)3 solid solutions with small Se content, whereas the higher Se content would lead to phase separation. We note that an alternative structure has been proposed for Sn2Se3 which was found to be 1.7 meV per atom (0.8 kJ mol–1 per F.U.) lower in energy than the Sn2S3 structure examined here,[33] although based on the present results this would still be unstable to phase separation.

These results demonstrate that the formation of Pnma Sn(S1–Se) and Sn(S1–Se)2 solid solutions should be facile. Rocksalt Sn(S1–Se) solid solutions are predicted to be metastable, although it may be possible to obtain them by epitaxial growth on a suitable substrate,[23,27] and the calculations predict that solid solutions with high Se content should be closer to the energetic convex hull. Sn2(S1–Se)3 solid solutions with low Se content should also be accessible. Finally, it is also worth noting that the mixing energy from alloying could potentially be used to stabilize a bulk phase and minimize or prevent unfavorable reactions with contact materials, such as have been identified as a major contributing factor in the poor performance of SnS PV devices.[21]

Structural Properties

The thermodynamic average X̅ (xSe) of a general physical property X at a given composition can be formed by weighting the properties of each structure in the solid-solution model by the occurrence probabilities P obtained from the partition function asThe averaged volume for a specific composition, for example, can be calculated according toOur calculations predict that the volume expansion of all four solid solutions is close to linear with Se content (Figure ). The relative expansion between the sulfide and selenide endpoints mirrors the chalcogen content, falling in the order of Sn(S1–Se)2 > Sn2(S1–Se)3 > Pnma Sn(S1–Se) ≈ rocksalt Sn(S1–Se).

Figure 3

Volume expansion in solid solutions of Pnma and rocksalt Sn(S1–Se), Sn(S1–Se)2, and Sn2(S1–Se)3 as a function of Se fraction xSe for a 900 K formation temperature. The markers show the calculated data points, and the solid lines are fits to the model function described in the text.

To quantify the trend, we fitted the (xSe) data to the model functionwhere V(xSe = 0) and V(xSe = 1) are the volume per formula unit of the sulfide and selenide endpoints, respectively, and b is a “bowing parameter” that captures the deviation from linearity. The fitted b values are small, ranging from 0.49% of the cell volume of the Sn2Se3 endpoint to 1.12% of the volume of Pnma SnS (see Table S2.1). Comparing the results for the Pnma Sn(S1–Se) and Sn(S1–Se)2 alloys against those obtained using the smaller supercell expansions (see the previous section) showed a difference in the averaged volume of <0.2 Å3 per formula unit, indicating that the structural properties are well converged with respect to the supercell size.

A comparison of the averaged lattice parameters (see Tables S2.2–S2.5) indicates that all four solid solutions maintain the parent lattice structure and suggests homogenous mixing. This was confirmed by calculating averaged pair-distribution functions (PDFs; gAB(r)) as a function of composition for each of the four phases. gAB(r) expresses the probability of finding an atom of species B between r and r + Δr from a reference atom of species A relative to the number density ρB = NB/VgAB(r) shows peaks at preferred interatomic distances and can be evaluated for pairs of a given atom type (i.e., A = B), for pairs of two different atoms (A ≠ B), or across all atomic pairs in the structure.

The all-atom and Sn–S/Sn–Se partial PDFs for five compositions in the Sn2(S1–Se)3 solid solution are shown in Figure . g(r) shows a series of sharp peaks at the various interatomic distances in the Sn2(S1–Se)3 structure, and the features corresponding to the tin–chalcogen distances can be assigned with reference to gSn–S(r) and gSn–Se(r). At intermediate compositions, the peaks split into two components that can be assigned to shorter Sn–S and longer Sn–Se distances, indicating that the chalcogen substitution leads to small local variations in the Sn coordination environment. Despite this, there is a smooth progression between the PDFs of the sulfide and selenide endpoints, further supporting a homogenous distribution of the chalcogen ions over lattice sites. Simulated PDFs of the Pnma and rocksalt Sn(S1–Se) and Sn(S1–Se)2 compositions show similar phenomena (see Figures S2.1–S2.3).

Figure 4

Calculated total (black) and Sn–S (orange)/Sn–Se (green) partial pair-distribution functions g(r) for Sn2(S1–Se)3 solid solutions with compositions xSe = 0 (a), 0.25 (b), 0.5 (c), 0.75 (d), and 1 (e). The partial correlation functions have been rescaled so that they represent the area under the total g(r) from all of the Sn–S and Sn–Se bond distances within a coordination sphere of radius r. The histograms of interatomic distances were generated using a bin width Δr = 0.01 Å and broadened with Gaussian functions of width σ = 0.05 Å.

The ability to fine-tune the cell volume and lattice parameters by varying the composition is noteworthy because it could be used to optimize the lattice match to other materials in a device structure. One example where this would be of value is in the recently reported use of SnS2 in place of CdS as a buffer layer in Cu(In1–Ga)Se2 and Cu2ZnSn(S1–Se)4 solar cells.[58,59] The homogenous distribution of the chalcogen ions is also noteworthy as the microscopic interfaces introduced by preferential clustering could negatively affect the electronic/thermal transport and/or act as recombination centers for electrons and holes. Furthermore, the ability to accommodate the lattice strain induced by the cation substitution suggests that these materials may be relatively defect tolerant, with further implications for doping, for example, to improve electrical properties for thermoelectric applications.[10,18]

Electronic Structure and Optical Properties

Several of the endpoints in the tin sulfide and selenide families have optoelectronic properties that make them well suited to important technological applications. The monosulfide and selenide SnS/SnSe have both been extensively studied as PV absorbers,[6−8,11,14−16] whereas SnS2 has recently been highlighted as a potential photocatalyst for water splitting,[11] and SnSe2 has been used in high-performance photodetectors.[12] The electronic structure of SnSe also plays a key role in its application as a high-performance thermoelectric.[9,10] However, save for a small number of experimental studies on tin sulfide/selenide thermoelectric materials,[35−37] relatively little attention has been given to the possibility of tuning the electrical properties by alloying. It is therefore of interest to investigate how forming solid solutions affects the electronic structure and optical responses of these systems.

High-throughput modeling studies must necessarily balance the accuracy of calculated properties against the computational cost of obtaining them. Modeling the electronic structure and related properties presents a particular challenge, as lower-level quantum-mechanical methods such as semilocal DFT invariably underestimate the size of the electronic bandgap[60,61] with a consequent effect on the dielectric function and derived properties such as the optical absorption coefficient. The PBEsol functional used for geometry optimization was found to underestimate the bandgaps of the endpoint structures by ∼0.5–1 eV with respect to experimental measurements. More seriously, although more accurate electronic-structure methods showed Sn2Se3 to be a narrow-gap semiconductor, PBEsol predicted a qualitatively incorrect near-metallic electronic structure. As discussed in detail in Section 3 of the Supporting Information, we compared the predicted bandgaps obtained with a range of electronic-structure methods against available experimental measurements and found that the SCAN meta-GGA functional struck a good balance between accuracy and cost, allowing us to evaluate the bandgap, density of states, and optical properties of the ∼5000 structures in our models with a reasonable accuracy.

Figure shows the thermodynamically averaged bandgaps of the four systems as a function of composition, with tabulated values in Tables S3.9–S3.12. The bandgaps fall roughly in the order of Sn(S1–Se)2 > Pnma Sn(S1–Se) > rocksalt Sn(S1–Se) > Sn2(S1–Se)3. Strikingly, for the Sn(S1–Se)2 solid solution (Figure c), the composition dependence is close to linear and suggests that the gap can be tuned through ∼0.9 eV by varying the ratio of S and Se. The calculated SCAN bandgaps of the SnS2 and SnSe2 endpoints are both underestimated at 1.6 and 0.7 eV, respectively, compared to the experimental values of 2.28 and 0.9 eV,[22,62] suggesting that a wider range of ∼1.38 eV would be achievable in practice. For a single-junction solar cell, the Shockley–Queisser limit points to a bandgap of 1.34 eV for the theoretical maximum conversion efficiency of 33.7%.[63] Presuming a simple linear variation of the gap with composition and taking the experimental gaps as the endpoints, our calculations suggest that this could be achieved by a solid solution with (2.28 – 1.34)/(2.28 – 0.9) ≈ 70% Se content, i.e., Sn(S0.3Se0.7)2. Taking into account the slight curvature of the graph in Figure c, our results predict that a lower Se content of about 60% would be optimum, i.e., Sn(S0.4Se0.6)2. The possibility of “tuning down” SnS2 for PV applications as an alternative to Pnma SnS is an interesting one, as the dichalcogenide system has several practical advantages: the 2D van der Waals structure should allow for a better control over the crystal growth and alignment, and the Sn(IV) oxidation state would avoid the need to carefully control the chalcogen chemical potential during the synthesis to prevent the introduction of phase impurities by “over oxidation” of the metal. We also note that the calculated mixing energies predict that the ∼1:2 SnS2/SnSe2 solid solution should be readily formed (cf. Figure ).

Figure 5

Composition dependence of the bandgap Eg of Pnma (a) and rocksalt Sn(S1–Se) (b), Sn(S1–Se)2 (c), and Sn2(S1–Se)3 (d) solid solutions. On each plot, the markers show the weighted averages calculated for a 900 K formation temperature and the shaded regions show the range of values within ±1 weighted standard deviation. The four subplots span the same 1 eV energy range to allow the direct comparison of the range and spread of the gaps but note the different starting values.

The calculations predict a much smaller variation in the bandgap of the monochalcogenide Sn(S1–Se) solid solutions (Figure a/b). The SCAN calculations give a range of ∼0.2 eV between the bandgaps of the Pnma SnS and SnSe endpoints of 0.85 and 0.66 eV, respectively (Figure a), compared to the experimental values of 1.06 and 0.95 eV.[22,62] As for Sn(S1–Se)2, the dependence is approximately linear with a comparably narrow spread, and the linear trend is consistent with experimental measurements on Sn(S1–Se) solid solutions.[35] The bandgap of Pnma SnS is below the Shockley–Queisser limit for optimum efficiency, and both the present calculations and experimental measurements on the endpoints suggest that alloying with SnSe will not improve upon this. On the other hand, the smaller bandgap could lead to improved electrical properties for thermoelectric applications, which is supported by experiments demonstrating larger carrier concentrations and higher mobilities toward the SnSe endpoint, leading to an improved electrical conductivity and a better thermoelectric performance.[35,36]

The rocksalt Sn(S1–Se) solid solutions show markedly different behavior (Figure b). The SnS and SnSe endpoints are predicted to have very similar bandgaps, and the calculations suggest a small reduction of ∼0.1 eV for compositions close to the Sn(S0.7Se0.3) midpoint. This behavior, together with the narrow spread, suggests that the rocksalt solid solutions would be tolerant to local variations in composition, although given its predicted metastability, the presence of Pnma impurities is likely to be a much larger issue.

Although the experimental bandgaps of Pnma SnS and Sn2S3 are very similar at 1.06 and 1.10 eV,[22] the predicted bandgap of the Sn2Se3 endpoint is considerably smaller than that of Pnma SnSe, suggesting a higher degree of compositional tunability in this system (Figure d). The SCAN calculations predict Sn2S3 to be a narrow-gap semiconductor with a 0.63 eV bandgap, whereas Sn2Se3 is predicted to be a near-metallic system with a gap of 0.19 eV, suggesting a highly tunable high electrical conductivity that, in conjunction with the “pseudo 1D” structure (see Figure ), would, in principle, make Sn2(S1–Se)3 solid solutions attractive for electronics applications. However, as noted above, stabilizing this solid solution may be challenging, and for this reason, we do not envisage these materials finding widespread applications.

To estimate the error on the SCAN bandgaps, the gaps of a subset of the structures were evaluated with more accurate single-shot (non-self-consistent) HSE 06 calculations using the SCAN wavefunctions (see Section 3 of the Supporting Information). Calculations on the endpoint structures indicated that these non-self-consistent hybrid calculations typically yield bandgaps within a few percent of the self-consistent results at an order of magnitude lower computational cost and compare favorably to available experimental measurements. The SCAN calculations underestimate the bandgaps of the Pnma Sn(S1–Se), Sn(S1–Se)2, and Sn2(S1–Se)3 models by 0.2–0.3, 0.4–0.6, and 0.2–0.3 eV (approx. 20, 30, and 45%), respectively, with larger relative errors for the smaller-gap Se-rich compositions (Tables S3.13, S3.15, and S3.16). Rocksalt Sn(S1–Se) is a notable exception, with smaller errors ranging from −21 to 35 meV and an average 0.5% overestimate of the gap (Table S3.14). The fact that the errors for each system appear to behave consistently lends a degree of confidence to the trends in Figure . In particular, HSE 06 calculations on selected Sn(S0.5Se0.5)2 and Sn(S0.25Se0.75)2 structures give bandgaps around 1.4 and 1.2 eV, respectively (see Table S3.15), in support of the 70% Se content predicted to give a target bandgap close to 1.34 eV, and the 0.93/0.37 eV HSE 06 gaps of Sn2S3 and Sn2Se3 support the qualitative conclusion above regarding the potential for tuning the electrical conductivity.

We also computed averaged bandgaps for the Pnma Sn(S1–Se) and Sn(S1–Se)2 solid solutions using smaller supercell expansions, as outlined above. For the monochalcogenide solid solution, we observed a large deviation of around 100 meV for all compositions, which suggests that this property is more sensitive to the chosen supercell expansion than the structural properties. However, the deviation also occurs for both endpoints, suggesting that the error is constant across the set of structures, and the relative changes in bandgaps are reasonably well converged. The two expansions tested for the dichalcogenide solid solution show smaller absolute differences of up to 46 meV. Crucially, for both systems, the trends shown in Figure are reproduced with both the supercell expansions tested, and so we would not expect our conclusions on the compositional tunability to be affected by our choice of cell size.

To establish the origin of the trends in Figure , the effect of composition on the band edges was investigated by calculating averaged electronic density of states (DoS) curves referenced to the Sn 1s core levels (Figures and S3.12–S3.14). The valance band edges in Sn(S1–Se)2 are composed predominantly of chalcogen states,[22] whereas the conduction band has contributions from both Sn and S/Se. Varying the composition raises the valence-band maximum (VBM) and lowers the conduction-band minimum (CBM) relative to the average Sn 1s core level, which together produces the narrowing of the bandgap seen in Figure c.

Figure 6

Electronic density of states (DoS) of the valence and conduction bands of Sn(S1–Se)2 solid solutions with Se fractions xSe = 0 (a), 0.25 (b), 0.5 (c), 0.75 (d), and 1 (e). The DoS curves are drawn as stacked area plots showing the projections onto Sn (blue), S (red), and Se (orange) atomic states. The energies are referenced to the average Sn 1s core-level eigenvalues, and the zero is set to the valence-band maximum of the SnS2 endpoint (xSe = 0; subplot (a)).

For the two monosulfide Sn(S1–Se) solid solutions, the substitution of sulfur with selenium lowers the energy of both band edges but shifts the CBM more prominently, which for the Pnma phase results in a reduction in the gap with increasing Se content, as shown in Figure a (Figures S3.12/S3.13). The large reduction in bandgap on going from Sn2S3 to Sn2Se3 arises primarily from a pronounced lowering of the CBM, whereas the VBM shows a subtle upshift in energy (Figure S3.14).

These analyses show that forming solid solutions provides a means to tune both the band-edge positions and the bandgap. It has previously been shown that poor band alignment to the contact materials in SnS-based PV devices is a likely contributor to the low efficiency realized experimentally,[21,64] and indeed forming a solid solution to improve the energy level alignment between SnS and the Zn(S1–O) buffer layer led to a significant improvement in the efficiency.[6] The present results demonstrate that alloying could be used to adjust the band offsets at interfaces for better device performance,[65] perhaps balancing this against conversion efficiency if required, or to optimize the electronic structure for photocatalytic reactions.[5]

For thermoelectric applications, the changes in the density of states at the band edges also suggest that alloying could be used to promote so-called “band convergence”, where increasing the density of states at the band edges allows the Seebeck coefficient to be increased without degrading the electrical conductivity. This is a proven mechanism for stabilizing the performance of thermoelectric devices over a wider range of operating temperatures.[3,66] This is in addition to the potential enhancements in carrier concentrations, mobilities, and electrical conductivity that can result from lowering the bandgap, albeit often at the expense of a reduced Seebeck coefficient, as well as the potential for higher defect tolerance in the mixed phases.

The electronic structure defines the energy-dependent complex dielectric function ε(E) = εRe(E) + iεIm(E), which in turn determines optical properties including the wavelength-dependent absorption coefficient α(λ), a key parameter for PV and photocatalysis applications. Using the SCAN electronic-structure calculations on the four solid solutions, we modeled the composition dependence of the dielectric functions,[50] allowing us to predict the variation in the wavelength-dependent refractive index n(λ), the extinction coefficient κ(λ), and the absorption coefficient α(λ) (Figures and S3.15–S3.18).

Figure 7

Calculated wavelength-dependent absorption coefficient α(λ) for Pnma (a) and rocksalt Sn(S1–Se) (b), Sn(S1–Se)2 (c) and Sn2(S1–Se)3 (d) solid solutions. The lines are color coded by Se fraction xSe from blue (dark; xSe = 0) to cyan (light; xSe = 1).

The optical properties of all four solid solutions are predicted to vary smoothly with composition, which again suggests that alloying should enable controlled tuning between the sulfide and selenide endpoints. There is a consistent red shift in the absorption profile with increasing Se content, as expected given the typical bandgap reductions in Figure , together with an enhancement of α(λ) at infrared and high-energy UV wavelengths and a decrease in the visible and mid-to-near UV. It should be noted that the underestimation of the bandgap with SCAN would lead to a red-shifted absorption profile, and thus it is likely that the enhanced infrared absorption would translate to an enhanced visible absorption in practice. On the other hand, the independent-particle approximation used to simulate the dielectric function tends to overestimate transition energies and blue shift the simulated absorption profile,[67] so there may be a partial cancellation of errors.

As for the bandgaps, we recalculated the DoS and optical properties of a subset of the structures in each solid solution using non-self-consistent HSE 06 (see Figures S3.19–S3.38). Comparing the SCAN and HSE 06 DoS curves shows a rigid shift of the band edges and some scaling of the energies of deeper valence and conduction states with relatively small changes to the feature intensity. The shifts of the band edges produce a uniform shift and a reduction in the intensity of the low-energy dielectric functions, confirming both that the absorption intensity shifts toward shorter visible wavelengths with more accurate electronic-structure techniques, and that the qualitative changes in band positions and optical properties with composition are reasonably well reproduced by the SCAN calculations.

Perhaps, the most important result in Figure is that the Pnma Sn(S1–Se) and Sn(S1–Se)2 solid solutions are predicted to have a comparable α, which suggests that the exemplary optical absorption properties of SnS should be similar in Sn(S1–Se)2, allowing it to be applied in a thin-film configuration. Indeed, experiments have highlighted the strong optical absorption of very thin SnSe2 flakes for photodetection.[12] For the more established SnS, the calculations suggest that alloying with SnSe would enhance the absorption in the longer-wavelength part of the visible spectrum, although this would be counteracted by a reduction in the absorption at shorter wavelengths, suggesting little overall benefit for PV applications.

Lattice Dynamics

Designing an optimized thermoelectric material entails balancing the electrical conductivity σ and Seebeck coefficient S against a low lattice and electronic thermal conductivity κlatt + κel, which is usually expressed by the dimensionless figure of merit ZT(3)κel is typically negligible in semiconducting systems, so a major focus of the thermoelectric research has been reducing κlatt while maintaining the power factor S2σ, leading to strategies such as exploiting phonon anharmonicity,[68−70] nanostructuring,[71] and performing atomic substitutions to increase the atomic mass variance.[3,72]

For a commercially viable thermoelectric device, materials with a ZT score above 2 are considered to be the benchmark. The lead chalcogenides, in particular, PbTe, have been widely researched due to their intrinsically low lattice thermal conductivity,[73] arising from strong phonon anharmonicity,[68,74] together with their favorable electrical properties.[66] However, although a ZT > 2 can be achieved with nanostructured PbTe,[71] the environmental toxicity of Pb and the rarity of Te are both major barriers to widespread application. More recently, bulk SnSe was shown to have a very high ZT of 2.6 at 923 K,[9] due, in part, to the strong phonon anharmonicity and very low lattice thermal conductivity.[69,70] Although the headline efficiency was only observed over a narrow temperature range, it was subsequently found that doping allows good thermoelectric performance to be achieved over a wide range of operating temperatures by optimizing the electronic structure.[10]

Alloying is a natural approach to optimizing thermoelectric performance, since the larger variations in atomic mass and chemical bond strength generally promote stronger phonon scattering and reduce the lattice thermal conductivity.[72] There are various reports of the thermoelectric properties of Sn(S1–Se) solid solutions in the literature, with ZT scores ranging from 0.64 to 1.67 in Ag-doped samples.[35−37] In particular, several studies have shown that solid solutions with 80% Se improve the thermoelectric efficiency relative to either of the endpoints by reducing the thermal conductivity and increasing the electrical conductivity at the expense of a slight reduction in the Seebeck coefficient.[35,36] More recent studies of Ag-doped Sn(S1–Se) solid solutions achieved a ZT of 1.67 at 823 K with 85% Se and attributed the high thermoelectric performance to the formation of nanoscale phonon-scattering point defects.[37]

As with accurate electronic-structure calculations, including lattice dynamics in high-throughput modeling studies is a technical challenge. To obtain an accurate set of force constants for evaluating the phonon frequencies using a finite-displacement approach requires up to 6 N accurate force calculations, depending on the crystal symmetry, on supercell expansions containing ∼100–200 atoms. Taking the system sizes and phonon supercells listed in Table , this “worst-case scenario” ranges from 192 calculations on 128-atom cells to 144 calculations on 192-atom cells, which is not practical for our full set of ∼5000 structures.

To obtain instead a qualitative idea of how the phonon spectrum changes with composition, we adopted the “structural fingerprint” method proposed by Valle and Organov as a measure of structural similarity.[75,76] This method uses a set of fingerprinting functions FAB(r) to characterize a structure in terms of its spectrum of interatomic distances, which is straightforwardly related to the pair-distribution functionA unique fingerprint for a given structure is obtained by concatenating the FAB(r) for all unique pairs of atoms (A, B), computed out to a radius rmax that captures all pair distances within a single crystallographic unit cell. The similarity between two structures can then be assessed by calculating the (vector) distance between their respective fingerprints, calculated for a common set of atomic pairs and rmax, for example using the cosine formulaWe note that the definition used in ref (75) includes an additional factor of 1/2 to adjust the distances to the range [0, 1]. Given the P for the structures in a given set, we can define a weighted average fingerprint vector (“center of mass”) asBy computing the cosine distances of structures in the set to F̅, we can then identify “most average” and “least average” structures, the phonon spectra of which we take as an approximate measure of the average and spread.

We performed this analysis for the 50/50 compositions of each of the four solid solutions, using a histogram bin width and Gaussian broadening of 0.05 and 0.02 Å, respectively, as suggested in ref (76). For all four systems, we found that the majority of the structures fell close to the center of mass (see Figures S4.1/S4.2). In the Pnma and rocksalt Sn(S1–Se) and Sn(S1–Se)2 solid solutions, the distances range from 0 to approx. 0.3, with >90% of the structures within a range of 0.1. The Sn2(S1–Se)3 solid solution is an exception, with structures spread over a wider range of distances from the center of mass and outliers at further extremes. This increased structural diversity is likely due to its more complex crystal structure and the presence of different chalcogen sites (cf. Figure ) as well as possibly to the relatively small number of atoms in the single cell used to build the solid-solution model.

Figure compares the simulated phonon spectra of Pnma SnS, SnSe, and representative average structures selected from the Sn(S0.5Se0.5) mixed phase. The phonon DoS curves of both endpoints consist of high- and low-frequency peaks corresponding predominantly to modes involving motion of the Sn and chalcogen atoms, respectively. The phonon modes in SnS span a larger range of frequencies than in the selenide, as expected given the larger difference in atomic mass and stronger chemical bonding.[77]

Figure 8

Phonon density of states (phonon DoS, g(v)) curves for Pnma Sn(S1–Se) solid solutions with Se fractions xSe = 0 (a), 0.5 (b), and 1 (c). The curves are drawn as stacked area plots showing the projections of phonon modes onto Sn (blue), S (red), and Se (orange) atoms. For the xSe = 0.5 composition in (b), the main curve shows the phonon spectrum of the most average structure, and the dashed black lines show ± the difference to the spectrum of the least average structure to give an indication of the expected spread (see text).

Comparing the endpoints to the mixed phase suggests that the DoS varies somewhat systematically with composition, although a comparison of the spectrum to a linear combination of the endpoints shows some degree of deviation (Figure S4.3). Taking the difference between the spectrum of the “most” and “least” average structures as a measure of the anticipated spread in the phonon frequencies shows that, for this compound, the variation is likely to be relatively small, which in turn suggests that, as for the bandgap, predictable tuning of the lattice dynamics through the composition should be possible.

Although a direct calculation of the lattice thermal conductivity is not feasible given the computational cost of doing so, the change in the phonon DoS with composition does provide a qualitative basis to explore the likely effect of alloying on heat transport in the solid solutions. The macroscopic thermal conductivity can be written as a sum of contributions from individual phonon modes λ according to[78]where Cλ is the model heat capacity, νλ is the mode group velocity, τλ is the phonon lifetime, V0 is the unit-cell volume, and N is the number of phonon wavevectors included in the summation. We note that the product of νλ and τλ gives the phonon mean-free path, which appears in a widely used alternative expression.

In semiconductors, the dominant contributors to finite phonon lifetimes are three-phonon scattering events, either collisions (two phonons in, one out) or decay (one in, two out). The lifetimes are inversely related to the phonon linewidths Γλ, which can be approximated as[78]Here, the quantity P–λ is the averaged three-phonon interaction strength and captures the physical coupling between interacting modes, and N2(qλ, ωλ) is a two-phonon density of states that gives the number of allowed energy- and momentum-conserving scattering events for a mode with wavevector qλ and frequency ωλ.

Spreading the phonon DoS over a larger range of frequencies, as in the 50/50 solid solution, should increase the number of energy-conserving three-phonon interactions without a large effect on the coupling strength, thereby resulting in an overall decrease of the thermal conductivity compared to either of the endpoints. The high-temperature (825 K) thermal conductivity of SnS and SnSe has been measured at ∼0.5 and 0.2–0.3 W m–1 K–1, respectively.[9,18] Measurements on Sn(S1–Se) solid solutions broadly show a reduction in the thermal conductivity with increasing Se content up to xSe = 0.8, with measured thermal conductivities ranging from 0.1 to 0.6 W m–1 K–1 depending on the sample preparation and doping.[35−37] However, none of these three studies provides measurements for pure SnSe synthesized under the same conditions, and the wide variability in the experimental results makes it difficult to ascertain whether or not the alloys can show lower thermal conductivities than the selenide endpoint. Nonetheless, the softening of the phonon spectra with composition suggested by our analyses is consistent with a reduction in the lattice thermal conductivity, tentatively supporting the experimental observations that alloying optimizes both the thermal transport and electrical properties for thermoelectric performance.[35,36]

A similar analysis of the rocksalt Sn(S1–Se), Sn(S1–Se)2, and Sn2(S1–Se)3 solid solutions (see Figures S4.4–S4.9) suggests similar trends, i.e., a systematic change in the phonon DoS between the sulfide and selenide endpoints with composition, albeit with notable deviations from linearity. We would, therefore, expect alloying to similarly allow for the manipulation of the thermal conductivity of the other sulfide/selenide phases.

As noted above, the lattice dynamics also contribute to the thermodynamic free energy. Within the harmonic approximation employed here, the Helmholtz free energy A(xSe, T) for a given composition, including lattice dynamics contributions, can be writtenwhere Ulatt is the lattice energy, Sconfig is the configurational entropy, and Uvib and Svib are the contributions to the internal energy and entropy due to the occupation of phonon modes. Evaluating Uvib and Svib quantitatively would require a statistical average over the full set of structures, which is not feasible. However, we can compare the vibrational free energies calculated for the endpoints and our selected xSe = 0.5 structures to obtain an approximate magnitude of the contribution to the mixing energies. Across the four systems examined in this work, we calculate contributions to Gmix(xSe, T) ranging from −0.75 to +0.03 kJ mol–1 atom–1 at 900 K which, with reference to Figure , are quite significant. For all but one of the xSe = 0.5 structures examined, the vibrational contribution decreases Gmix, suggesting that differences in lattice dynamics would further stabilize the mixed phases. However, since our fingerprinting technique is designed to select configurations with structural, and not necessarily energetic diversity, we would need to perform a full set of phonon calculations to confirm this.

Finally, we note that the relatively common use of alloying to optimize thermoelectric performance indicates that the exploration of methods for analyzing more quantitatively the effect of alloying on the lattice dynamics is an important future development area. A potentially promising direction here would be to use the energetics calculations on the mixed phases to parameterize an empirical force field to perform these calculations, which could then allow for fully thermodynamically averaged phonon spectra and/or lattice thermal conductivities to be modeled with moderate computational cost, although this is beyond the scope of the present work.

Conclusions

In summary, we have carried out a detailed modeling study of four tin sulfide/selenide solid solutions, viz., Pnma and rocksalt Sn(S1–Se), Sn(S1–Se)2, and Sn2(S1–Se)3, with a view to using alloying to optimize the physical properties for energy applications.

The calculated phase diagrams indicate facile formation with mixing strongly favored by configurational entropy. We predict Pnma Sn(S1–Se) and Sn(S1–Se)2 solid solutions to be accessible across the full range of compositions through typical high-temperature syntheses, whereas rocksalt Sn(S1–Se) and Sn2(S1–Se)3 solid solutions with high and low selenium content, respectively, may also be accessible under certain synthesis conditions. Structural analyses of the solid-solution models indicate a homogenous distribution of the chalcogen atoms over lattice sites, resulting in a predictable variation of the volume and lattice parameters with composition that should help with lattice matching to reduce the interfacial strain in device structures.

Our calculations predict a systematic composition dependence of the bandgap, band-edge positions, and dielectric properties in the Pnma Sn(S1–Se) and Sn(S1–Se)2 solid solutions, providing a means to fine-tune the optoelectronic properties for their proven photovoltaic and photocatalysis applications. We also predict that the Sn(S1–Se)2 bandgap can be tuned down to the theoretical optimum for PV while maintaining a similar absorption coefficient to Pnma SnS. To the best of the authors’ knowledge, Sn(S1–Se)2 solid solutions have not been widely investigated for photovoltaic applications, but the present work suggests that they could be a superior alternative to SnS: the crystal growth direction could be better controlled, given the 2D structure, and the full oxidation of Sn would reduce the need to tightly control the chalcogen chemical potential during synthesis. After the years of optimization that have been spent on achieving the ∼5% photovoltaic efficiency of SnS, we believe that this result merits further investigation. For the Pnma Sn(S1–Se) solid solution, on the other hand, increasing the Se content decreases the bandgap from what is already below the optimum for photovoltaic performance while not significantly improving the optical absorption coefficient.

Our calculations also show that alloying can allow the band levels to be adjusted for a better electrical match to other device components, potentially allowing this to be balanced with other properties if required. In a similar vein, the entropic stabilization obtained by alloying may reduce or prevent undesirable chemical reactions at interfaces with contact materials.

Our calculations support the experimental observations that Pnma Sn(S1–Se) solid solutions with high Se content can improve the thermoelectric properties over the champion SnSe material, through a mixture of enhancing the electrical conductivity by reducing the bandgap and reducing the thermal conductivity by enhancing the phonon scattering. Similar trends in the phonon spectra of the other Sn(S1–Se) phases with composition suggest that the thermal transport in these systems could also be controlled by alloying.

The predicted metastability of the majority of the Sn2(S1–Se)3 and all of the rocksalt Sn(S1–Se) compositions implies that preparing and working with these materials in bulk would be challenging. On the other hand, both are at some compositions close enough to the convex hull that they may be of concern for phase purity. We consider it unlikely that either of these will find widespread applications in the near future.

Finally, from a methodological standpoint, our study shows that first-principles modeling of the structural and electrical properties of relatively complex solid solutions is within the reach of modern quantum-chemical techniques and can provide valuable information to guide experimental syntheses. Moreover, given the importance of alloying in the optimization of thermoelectric materials, future developments toward quantitative modeling of the composition dependence of the lattice dynamics and thermal transport will make this a valuable addition to the growing capabilities of computational materials design.