Thermodynamics, Electronic Structure, and Vibrational Properties of Snn(S1–xSex)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–xSex), Sn(S1–xSex)2, and Sn2(S1-xSex)3 solid solutions. All four show an energetically favorable and homogenous mixing at all compositions, but rocksalt Sn(S1–xSex) and Sn2(S1–xSex)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–xSex) 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.

.2 Calculated free energies as a function of Se fraction Se for the conversion of Pnma Sn(S1-xSex) solid solutions to the rocksalt phase (a) and for the conversion of Sn2(S1-xSex)3 solid solutions to Pnma Sn(S1-xSex) and Sn(S1-xSex)2 phases with equivalent composition (b). We note that we do not consider the possibility of decomposition into phases with different composition (e.g. Sn2(S0.33Se0.66)3 → SnS + SnSe2), since the mixing energies in Fig. S1.1/ Fig. 2 in the text predict that both the Sn(S1-xSex) and Sn(S1-xSex)2 products should readily form solid solutions. As in Fig. S1.1, the line colours denote synthesis temperatures from 300 K (blue) to 1100 K (orange). 3 ] Pnma Sn(S1-xSex) 183.2 202.9 2.05 Rocksalt Sn(S1-xSex) 46.59 51.72 -0.293 Sn(S1-xSex)2 65.40 73.97 0.595 Sn2(S1-xSex) 3 444.0 500.3 -2.45 Table S2. 1 Fit parameters for the variation of the cell volume as a function of Se fraction Se to the model function in Eq. 6 in the text. All parameters are given with respect to the primitive unit cell used to build the supercells when forming the solid solution; note that for rocksalt Sn(S1-xSex) this is the two-atom primitive rather than the eight-atom conventional cell.  Table S2. 2 Thermodynamically averaged lattice parameters for Pnma Sn(S1-xSex) solid solutions as a function of Se fraction Se . For each of the averaged unit-cell parameters ̅, ̅ , ̅ , ̅, ̅ and ̅ , we also report the weighted standard deviation as a measure of the expected spread.

Structural properties
-  Table S2.5 Thermodynamically averaged lattice parameters for Sn2(S1-xSex)3 solid solutions as a function of Se fraction Se . For each of the averaged unit-cell parameters ̅, ̅ , ̅ , ̅, ̅ and ̅ , we also report the weighted standard deviation as a measure of the expected spread.   -Page 12 -

Electronic structure and optical properties
As noted in the text, accurate electronic-structure calculations pose a technical challenge to high-throughput modelling studies, since density-functional theory (DFT) with semi-local generalised-gradient approximation (GGA) functionals tend to underestimate the electronic bandgap 1,2 and incorrectly predict related properties such as the dielectric function and optical-absorption coefficient. In practice, this deficiency is usually mitigated with so-called hybrid functionals that incorporate a fraction of the non-local exact exchange energy and yield more accurate electronic structures, but at a much higher cost and less favourable scaling with system size.
Since it is not practical to perform hybrid calculations on the 1000s of structures in our four solid solution models, we examined the possibility of selecting a smaller subset of representative low-energy structures. However, plotting the distribution of the occurrence probabilities over the structures in each system/composition (Fig. S3.1) showed a small difference of just an order of magnitude between largest and smallest occurrence probabilities, which we do not consider sufficient to justify excluding a significant fraction of the structures from the thermodynamic averages.
In a bid to balance accuracy and computational cost, we tested the predictive capabilities of five electronicstructure methods on the primitive cells of the sulphide and selenide endpoints of each system (i.e. Pnma and rocksalt SnS/SnSe, SnS2/SnSe2 and Sn2S3/Sn2Se3). Taking structures optimised with PBEsol+D3, we compared the band gaps obtained with the PBEsol GGA, 3 the SCAN 4 and mBJ 5,6 meta-GGA functionals and the HSE 06 hybrid functional. 7 We note in passing that the empirical DFT-D3 dispersion correction 8 does not correct the electronic structure, and thus the band energies from "bare" PBEsol and PBEsol-D3 calculations on a given structure are equivalent. We also tested the G0W0 method, 9,10 which gives a non-self-consistent perturbative correction to the energy levels obtained with another method, in our case either PBEsol or SCAN. In principle, hybrid functionals can also be used in a similar non-selfconsistent manner, i.e. by taking self-consistent wavefunctions obtained with a less demanding functional and recalculating the band energies with the hybrid. Although not common, this approach has been shown to be suitable for high-throughput or screening calculations, 11,12 so we included it in our testing.
Γ-centred Monkhorst-Pack meshes with 8×4×8, 8×8×8, 8×8×4 and 4×8×3 subdivisions were employed in the calculations on the primitive cells of Pnma and rocksalt SnS/SnSe, SnS2/SnSe2 and Sn2S3/Sn2Se3 respectively. For the HSE 06 calculations, the plane-wave cutoff was reduced to 300 eV, a reduced FFT grid was used to evaluate the Hartree-Fock exchange (the PRECFOCK = Fast setting in VASP 13 ), and the exact exchange was evaluated at a subset of the kpoints in the sampling meshes (EVENONLY = .TRUE.). From extensive testing, we found that these substantially reduced the computational requirements of the calculations without significant loss of accuracy. Similarly, we found that in the G0W0 calculations the cutoff for the basis used to represent the response functions could be reduced to 200 eV and the exact-exchange contribution could be evaluated with a reduced FFT grid; however, attempting to reduce the k-point mesh led to a notable loss of accuracy, so we did not do so. We found that the convergence of the bandgap with respect to the number of virtual states (empty bands) was somewhat slow, so a total of 528, 480, 504 and 240 bands were included in the Pnma and rocksalt SnS/SnSe, SnS2/SnSe2 and Sn2S3/Sn2Se3 calculations, respectively, of which the quasiparticle corrections were computed for the first 48, 24, 24 and 120 bands, respectively, covering the valence and low-lying conduction bands.
SnS2 is a wide-bandgap semiconductor with an extrapolated 0 K indirect gap of 2.28 eV. 14 The PBEsol gap of SnSe2 has direct and indirect gaps of 1.7 and 0.9 eV, respectively. 15 The Although mBJ was designed to predict accurate bandgaps, and typically produces results close to the HSE 06 values at a much smaller cost, it does not appear to reproduce the shape of the band dispersion as well; a comparison of the calculated SnS2 DoS between functionals shows that the smaller PBEsol and SCAN gaps result primarily from a uniform shift of the valence and conduction bands compared to the HSE 06 calculations (Fig. S3.11), whereas mBJ predicts a markedly different shape to the conduction-band edge.
-Page 14 -Based on all these considerations, we opted to perform a full set of electronic-structure calculations using SCAN, and to perform non-self-consistent calculations on a subset of structures to obtain an estimate of the error.
Finally, to simulate optical properties we used the SCAN electronic-structure calculations to obtain the complex dielectric function ( ) = Re ( ) + Im ( ) within the independent-particle approximation using the linear-optics routines in the Vienna Ab initio Simulation Package code. 16 The calculated wavelength-dependent refractive index ( ), extinction coefficient ( ) and absorption coefficient ( ) were then calculated as: where | ( )| = √ Re 2 + Im 2 in Eq. 3.1 and 3.2 is the complex modulus.       Table S3.14 Comparison of the bandgaps g of five structures in the rocksalt Sn(S1-xSex) solid model solution obtained with the SCAN meta-GGA and from a single-shot non-self-consistent HSE 06 hybrid calculation based on the SCAN wavefunctions (see Section 3c in the text). For each structure, the first three columns give the Se fraction Se , spacegroup (S. G.) and structure number in the composition/spacegroup set, the fourth and fifth columns give the calculated bandgaps, and the last two columns give the absolute and percentage differences.  Table S3. 15 Comparison of the bandgaps g of five structures in the Sn(S1-xSex)2 solid model solution obtained with the SCAN meta-GGA and from a single-shot non-self-consistent HSE 06 hybrid calculation based on the SCAN wavefunctions (see Section 3c in the text). For each structure, the first three columns give the Se fraction Se , spacegroup (S. G.) and structure number in the composition/spacegroup set, the fourth and fifth columns give the calculated bandgaps, and the last two columns give the absolute and percentage differences.  Table S3. 16 Comparison of the bandgaps g of five structures in the Sn2(S1-xSex)3 solid model solution obtained with the SCAN meta-GGA and from a single-shot non-self-consistent HSE 06 hybrid calculation based on the SCAN wavefunctions (see Section 3c in the text). For each structure, the first three columns give the Se fraction Se , spacegroup (S. G.) and structure number in the composition/spacegroup set, the fourth and fifth columns give the calculated bandgaps, and the last two columns give the absolute and percentage differences.                  Sn(S1-xSex)2 (c) and Sn2(S1-xSex)3 solid solution models with composition Se = 0.5. The structures are ordered by occurrence probability to emphasise the similarity of the most likely atomic arrangements in most of the systems.