Stability of SnSe 1 − x S x solid solutions revealed by first-principles cluster expansion

The configurational thermodynamics of a pseudo-binary alloy SnSe 1 − x S x in the Pnma phase is studied using first-principles cluster-expansion method in combination with canonical Monte Carlo simulations. We find that, despite the alloy having a tendency toward a phase decomposition into SnSe and SnS at 0 K, the two constituent binaries readily mix with each other to form random SnSe 1 − x S x solid solutions even at a temperature below room temperature. The obtained isostructural phase diagram of SnSe 1 − x S x reveals that the alloy is thermodynamically stable as a single-phase random solid solution over a whole composition range above 200 K. These findings provide a fundamental understanding on the alloying behavior of SnSe 1 − x S x and bring clarity to the debated clustering tendency in this alloy system.

S Supplementary material for this article is available online (Some figures may appear in colour only in the online journal)

Letter
Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. dependent on the mixing thermodynamics of SnSe 1−x S x alloy, whether it displays at a given temperature a tendency of ordering to form ordered structures or a tendency of clustering toward phase segregation or a tendency of mixing to form a single-phase solid solution.
Notwithstanding several reports on successful syntheses of SnSe 1−x S x solid solutions [4,5,[11][12][13][14][15], the recent study on the alloying behavior of SnSe 1−x S x at the atomistic level, based on the quantitative analysis of scanning tunnelling microscopy measurements and density functional theory calculations, suggested that SnSe 1−x S x alloy has a tendency toward local phase segregation into SnSe and SnS, rather than forming a disordered solid solution [16]. This gives rise to an ambiguity in the mixing behavior of SnSe 1−x S x alloy. In order to clarify the issue, the thermodynamic of mixing between SnSe and SnS needs to be further investigated and comprehended, so that the alloy can be optimally designed and fully functional in its applications.
In the present work, we aim at revealing the alloying behavior of SnSe 1−x S x in the Pnma structure as a function of temperature using the cluster-expansion (CE) formalism [17]. Since the alloying process in SnSe 1−x S x is particularly involved with the chalcogen sublattice, containing both Se and S atoms, we exclude the Sn atoms, residing at the group-IV sublattice sites, from the cluster expansion, and consider them as spectator atoms. Consequently, the mixing energy ∆E mix (σ) of SnSe 1−x S x of a given atomic arrangement (σ) on the lattice with SnSe and SnS contents ( x SnSe and x SnS = 1 − x SnSe ) can be written as where the energy of SnSe 1−x S x alloy of a given σ, E(σ), is expanded into a sum over correlation functions ξ Although the expansion, expressed in equation (2), is mathematically complete in the limit of inclusion of all possible f, it has to be truncated for practical purposes. Here, we use the MIT Ab initio Phase Stability code [18], as implemented in the Alloy-Theoretic Automated Toolkit (ATAT) [19], to truncate the expansion and to determine the ECIs in such a way that equation (2) returns the energies E(σ) of SnSe 1−x S x as close to those obtained from first-principles calcul ations as possible for σ, included in the expansion. The first-principles total energies are, on the other hand, calculated from the density functional theory (DFT), where the projector augmented wave method [20] as implemented in the Vienna ab initio simulation package [21,22] is used and the exchangecorrelation effects are modeled using the generalized gradient approximation [23]. The valence electron configurations, used for pseudopotentials, are 5s 2 5p 2 , 4s 2 4p 4 , and 3s 2 3p 4 for Sn, Se, and S, respectively. Since the standard DFT calculations cannot accurately describe the van der Waals interactions, we use the correction, proposed by Grimme (DFT-D3) [24], to account for the interactions. To minimize the energy, the internal atomic coordinates, volume, and cell shape of all included σ are fully relaxed. We also assure the calculated energies are converged within an acc uracy of 1 meV/atom with respect to the plane-wave energy cutoff and the number of Monkhorst-Pack k-points grids for the Brillouin zone integration [25].
To investigate the configurational thermodynamics of SnSe 1−x S x alloy, a database of different configurations σ of SnSe 1−x S x alloy has to firstly be established. To this end, we employ an algorithm developed by Hart and Forcade [26] to generate a set of 3002 configurations of SnSe 1−x S x alloy with up to 24 atoms in the primitive supercell, equivalent to 3 primitive orthorhombic unitcells. We then single out the first few hundreds of σ, calculate their total energy using first-principles approach, and include them in the cluster expansion to determine ECIs. For this particular case, the ECIs are considered up to 3-site interactions. The initial ECIs, extracted from the first expansion, are employed to predict the total energy of all generated σ through equation (2). We note that those initial ECIs may not do the prediction accurately. Thus, to improve the predictive power of cluster expansion, we use the CE-predicted energies as a guideline to single out a few more tens or hundreds of σ, not included in the first expansion, calculate their total energy, and include them, together with those from the first expansion, to redetermine the ECIs in the next cluster expansion. This procedure can be iterated, until the cluster expansion of desired quality is reached. The final expansion includes 476 σ and employ 40 ECIs. That is, apart from the 0-site and 1-site interactions, the ECIs are composed of 23 2-site interactions and 15 3-site interactions. The final expansion fits the 476 input σ with the cross-validation score of 0.195 meV/f.u. For comparison to the cluster expansion approach, we also calculate through first-principles approach ∆E mix of completely random solid solutions of SnSe 1−x S x with the composition x = 0.25, 0.50, and 0.75, modeled within 64 atom supercells by using the special quasirandom structure (SQS) method [27]. The ground-state diagram at 0 K of SnSe 1−x S x alloy is displayed in figure 2. Our calculations reveal that ∆E mix of SnSe 1−x S x alloy, evaluated with respect to SnSe and SnS, is positive for all considered σ, including those modeled by the SQS method. This implies a tendency toward a phase decomposition of SnSe 1−x S x into SnSe and SnS under thermodynamic equilibrium. As far as we are aware, no ordered structure of SnSe 1−x S x alloy has been reported so far in the literature, which is in line with our prediction. The next step is to investigate the configurational thermodynamics of SnSe 1−x S x alloy. To do so, we utilize the obtained ECIs in canonical Monte Carlo (MC) simulations using the Easy Monte Carlo Code (EMC2) [28], as implemented in the ATAT [19]. In this work, the simulation boxes of 16 × 15 × 6 orthorhombic primitive unitcells (5760 chalcogen atoms) are used. The simulations are performed at fixed compositions x, where 0 x 1 and ∆x = 0.05. At each composition, the SnSe 1−x S x alloy is cooled from 3000 to 25 K, using simulated annealing ∆T = 25 K and at each temperature, the simulations include 15 000 MC steps for equilibrating the system and then 12 000 more steps for obtaining the proper average of ∆E mix and configurational specific heat (C V ) as a function of temperature and alloy composition. The configurational thermodynamics of SnSe 1−x S x alloy is then evaluated through the Gibbs free energy of mixing, ∆G mix , given by where the mixing entropy ∆S mix can be obtained from thermodynamic integration of C V : The term ∆S MF mix stands for the mixing entropy per formula unit of the ideally random solid solution, stable in the limit of T → ∞, and can be obtained from the mean-field approach to be For this particular case, we assume that ∆S mix (x, T = 3000 K) ≈∆S MF mix (x), and thus the thermodynamic integration in equation (4) is performed from this high temperature downwards to the temperature of interest. The mixing free energies ∆G mix of a completely random solid solution of SnSe 1−x S x modeled with the SQS method are, on the other hand, obtained from the mean-field approximation of the mixing entropy, given by equation (5).
The resulting ∆G mix (x, T) curves of SnSe 1−x S x alloy, calculated by using the MC approach and the mean-field approximation, are displayed in figure 3(a). As can be seen from figure 3(a), the ∆G mix , derived from the MC approach (open circle), has a positive curvature for a whole composition already at T = 300 K and above. This clearly suggests a formation of continuous solid solution of SnSe 1−x S x alloy. Interestingly, the results, calculated within the mean-field approximation (shaded squares), are close to those, derived from the MC approach, which is typically the case only at temperature considerably above any phase transitions. By applying common tangent construction to the ∆G mix curves at different fixed temperature, we outline an isostructural T-x phase diagram of SnSe 1−x S x alloy, as depicted by figure 3(b). The phase diagram exhibits a miscibility gap at T 200 K, in which a continuous solid solution of SnSe 1−x S x , as denoted by α phase and thermodynamically stable at T 200 K, decomposes into a mixture of two phases of SnSe 1−x S x of different compositions x, given by α and α phases. Although SnSe 1−x S x alloy displays a tendency toward a phase decomposition into SnSe and SnS at T = 0 K, as demonstrated by positive ∆E mix at all compositions, see figures 2 and 3(a), SnSe readily mixes with SnS to form random solid solutions of SnSe 1−x S x at T > 0 K. This can be attributed to the increasingly strong contribution of the configurational entropy −T∆S mix to the Gibbs free energy of mixing ∆G mix as the temperature increases, and also to weak ECIs in the system. The latter gives rise to a driving force for SnSe 1−x S x to form a random solid solution, stable in the limit V/T → 0, where V is the strongest interaction in the system [29]. We note that the magnitude of the ECIs, derived from our cluster expansion of SnSe 1−x S x , is smaller than 1 meV/f.u. (see also supplementary figure S1 (stacks.iop.org/JPhysCM/30/29LT01/ mmedia)), and thus the criterion V/T → 0, characterizing atomic configuration of the random alloy, can be fulfilled at elevated temperature for this par ticular alloy. We note that, in the present work, the vibrational contributions to ∆G mix of SnSe 1−x S x alloy, induced by the lattice vibrations, are neglected. As discussed above, SnSe 1−x S x alloy has been predicted to be thermodynamically stable as a single-phase random solid solution over a whole composition range at T 200 K. This temperature is way below the melting temper atures of both SnSe and SnS, i.e. 1134 K and 1155 K, respectively [15]. Moreover, the vibrational contributions are typically of minor importance for isostructural alloys, as compared to the configurational parts [30]. As a consequence, the influence of lattice vibrations in this case should be very tiny at such a low temperature and does not have a significant impact on the phase stability of SnSe 1−x S x alloy, predicted in the present work (figure 3).
Our prediction is indeed in line with the experimental syntheses of a single-phase solid solution of SnSe 1−x S x alloy, generally performed at high temperature (T > 800 K) [4, 5,  11 -15]. This is, however, in contrast to the recent quantitative analysis of scanning tunnelling microscopy measurements of SnSe 1−x S x alloy, suggesting that the alloy has a tendency toward local phase segregation into SnSe and SnS, rather than forming a random solid solution [16]. Even though SnSe 1−x S x alloy is predicted in accordance with our calculations to display a clustering tendency to form SnSe and SnS as T → 0 K, it seems unlikely that such a situation could be reached in experiment due to a lack of atomic diffusion at low temperature.
Concerning the positive mixing energies of the random solid solution of SnSe 1−x S x , modeled in [16], we would like to point out that those energies were carried out only at T = 0 K. This is in fact identical to our calculations of ∆G mix (x, T = 0 K), as shown in figure 3(a). However, we have just demonstrated above that the random SnSe 1−x S x alloy is favored, if the influence of temperature on the mixing thermodynamics of SnSe 1−x S x alloy is taken into consideration. It is thus questionable whether SnSe 1−x S x alloy really exhibits a tendency toward a phase segregation into SnSe and SnS, as reported in [16]. This, however, requires further experimental elaborations.
The calculated lattice parameters of SnSe (SnS) are a = 11.571(11.277) Å, b = 4.180(4.008) Å, and c = 4.486(4.345) Å. These values are in excellent agreement with the experiments [4,31], in which they differ by less than 1%, as compared to the corresponding experimental values. The lattice parameters of the SnSe 1−x S x alloy, modeled using the SQS method, are found to slightly deviate from the Vegard's law by less than 1.8% (see supplementary figure S2). We note that the small deviation of the lattice parameters of SnSe 1−x S x alloy from the Vegard's law is also observed in experiment [5]. On the other hand, the band gaps of SnSe and SnS are estimated in this work to be 0.6 and 0.75 eV, respectively. The calculated band gaps are, however, underestimated by approximately 30% with respect to the experimentally measured values, i.e. 0.9 and 1.08 eV for SnSe and SnS, respectively [15]. This can be attributed to the used of the generalized gradient approximation to account for the exchange-correlation effects [32]. The band gap of the SnSe 1−x S x alloy is found to increase with the sulfur content, together with a positive deviation from the Vegard's law up to 6% (see supplementary figure S3). As a compliment, the calculated electronic density of states, showing the estimated bandgap, of SnSe, SnS and SnSe 1−x S x alloys, modeled using the SQS method, is also provided in the supplementary material (see figures S4 and S5).
In conclusion, we investigate the configurational thermodynamics of SnSe 1−x S x in the Pnma phase by using the first-principles cluster expansion in combination with the MC simulations as well as the mean-field approximation. The two approaches reveal that, although SnSe 1−x S x alloy displays a tendency toward a phase decomposition into SnSe and SnS at T = 0 K, the two constituent binaries readily mix with each other to form random solid solutions of SnSe 1−x S x at T > 200 K. This can be explained through the increasingly strong contrib ution of the configurational entropy to the mixing Gibbs free energy as the temperature increases, and the weak effective interactions in this particular alloy system. These findings are not only in line with the high-temper ature syntheses of a single-phase solid solution of SnSe 1−x S x , reported in the literature [4,5,[11][12][13][14][15], but also put in question the report of clustering in [16], or at least a thermodynamic explanation for such clustering.  All of the calculations are carried out using supercomputer resources provided by the Swedish National Infrastructure for Computing (SNIC) performed at the National Supercomputer Centre (NSC) and the High Performance Computing Center North (HPC2N).