Lattice dynamics of Pnma Sn(S_1-xSe_x) solid solutions: energetics, phonon spectra and thermal transport

Alloying is widely used as a means to fine-tune the properties of thermoelectric materials by reducing the lattice thermal conductivity. However, the effects of compositional variation on the lattice dynamics of alloy systems are not well understood, due in part to the difficulty of building realistic first-principles models of structurally-complex solid solutions. This work builds on our previous study of Sn n (S 1- x Se x ) m solid solutions [Gunn et al. , Chem. Mater. 31 , 10 , 3672, 2019 ] to explore the lattice dynamics of the Pnma Sn(S 1- x Se x ) system, which has been widely studied for potential thermoelectric applications. We find that the vibrational internal energy and entropy have a large quantitative impact on the mixing free energy and are likely to be particularly important in alloy systems with competing phases. The thermodynamically-averaged phonon dispersions and density of states curves show that alloying preserves the structure of the low-frequency bands of modes associated with the Sn sublattice but broadens the high-frequency chalcogen bands into a near-continuous spectrum at the 50/50 mixed composition. This results in a general reduction in the phonon mode group velocities and an increase in the number of energy-conserving scattering channels for heat-carrying low-frequency modes, which is consistent with the decrease in thermal conductivity observed in experimental measurements. Finally, we discuss some of the limitations of our first-principles modelling approach and propose methods to address these in future studies.

A c c e p t e d M a n u s c r i p t

Introduction
Around two thirds of the energy used globally is wasted as low-grade heat from sources including industry, refrigeration and transportation [1]. Thermoelectric (TE) power provides a means to recover electrical energy from temperature gradients and thereby enable significant efficiency gains in these energy-intensive processes, making it an important complement to primary renewable energy sources such as photovoltaics (PV). In particular, in the drive to meet ambitious climate change targets, emissions from internal combustion engines are becoming a huge area of concern, and employing thermoelectric generation as part of a mild hybrid powertrain could provide an interim solution while electric vehicle technology and infrastructure matures to allow for more widespread adoption.
The efficiency of a thermoelectric material is commonly expressed by the dimensionless figure of merit [2]: The Seebeck coefficient , electrical conductivity and electronic thermal conductivity are interdependent el such that the balance of a high power factor and low is generally found in doped semiconductors. At 2 el low to moderate temperatures, the lattice thermal conductivity makes the most significant contribution to latt the denominator, and thus minimising has proven to be a facile means of improving TE performance. In latt some cases the electronic and thermal transport can be almost completely decoupled [1], as in the "phonon glass, electron crystal" concept [3]. All four parameters in Eq. 1 are implicit functions of temperature, requiring that materials are either optimised to produce a high peak at a target operating temperature range or tuned to display a high across a wide range of temperatures.
Current flagship TE materials include PbTe, SnSe and Bi 2 Te 3 , all three of which display favourable electrical properties and intrinsically low thermal conductivity due to a combination of heavy elements and strongly anharmonic lattice dynamics [4][5][6]. However, PbTe and Bi 2 Te 3 are not suitable for widespread adoption due to the low abundance of Te, and there are also concerns with SnSe due to the environmental toxicity of Se. There has therefore been significant effort devoted to alternative systems including the more earth-abundant SnS [2] and metal oxides [7][8][9][10].
Alloying is commonly used as a means to enhance thermoelectric performance, as a suitable choice of components can maintain or improve a favourable electronic structure while reducing by introducing latt variation in atomic masses and chemical bond strength to promote stronger phonon scattering [11]. Pb(S,Se,Te), Sn(S,Se) and (Bi,Sb) 2 (Se,Te) 3 alloys have all been studied as thermoelectrics and the alloying shown to improve [12][13][14][15][16][17]. Due to the record-breaking high-temperature of SnS, the Te-free Sn(S,Se) alloys are of particular interest, and experimental studies have demonstrated a reduction in the thermal conductivity of mixed compositions [12,14].
Thermoelectricity is unique in that all four of the terms in the figure of merit are amenable to calculation using first-principles electronic-structure methods such as density-functional theory [18]. There have been a number of theoretical studies on existing and potential new single-component bulk thermoelectrics [6,[19][20][21][22][23], and there have recently been efforts to use high-throughput modelling to screen large numbers of compounds in a bid to identify new candidate thermoelectrics [24,25].  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 A c c e p t e d M a n u s c r i p t -Page 3 -Studying multicomponent alloy systems with first-principles methods is challenging, however, as building a realistic model requires calculations on a large number of configurations that for complex systems with low-symmetry unit cells can easily approach the cost of a screening study [26]. Ektarawong and Alling used a first-principles cluster expansion method to study the Pnma Sn(S 1-x Se x ) system and confirmed the stability of the mixed phases at finite temperature [27]. We recently carried out a study of four Sn n (S 1-x Se x ) m alloys [26], including the Pnma and rocksalt monosulphides, by performing first-principles calculations on the full sets of symmetry-inequivalent structures formed by substituting chalcogen atoms in small supercells of the parent structures [28]. In this study we also confirmed the stability of the mixed Pnma Sn(S 1-x Se x ) alloys and further predicted the effect of alloying on the structural, electrical and optical properties [26].
While many large-scale modelling studies make use of (semi-)local DFT methods to evaluate and compare energetics and electrical properties [29][30][31][32], modelling lattice dynamics and in particular thermal transport can be much more resource intensive, and there are thus comparatively fewer examples of screening studies focussing on these properties [33]. As a result, despite the widespread use of alloying to tune the performance of thermoelectric materials, there has been little theoretical investigation into how composition affects the dynamics and related properties. In this work, the model developed in Ref. [26] is extended to include explicit evaluation of the lattice dynamics of the Pnma Sn(S 1-x Se x ) solid solution. We quantify the effect of the phonon free energy on the mixing energetics, establish the effect of alloying on the phonon spectrum and frequency dispersion, and explore some of the implications for the thermal transport. Finally, we discuss the broader applicability of this modelling technique and possible future enhancements.

Computational modelling
The starting point for our calculations was the set of optimised structures in the Pnma Sn(S 1-x Se x ) solid solution model from our previous study [26]. In our previous work, the Transformer code [34] was used to enumerate the symmetry-independent structures formed by successively substituting the chalcogen atoms in a 32-atom 2 1 2 expansion of the Pnma SnS/SnSe structure, yielding a total of 2,446 unique structures × × across seventeen compositions from = 0 (SnS) to = 1 (SnSe). In the present study, lattice dynamics The accuracy of the phonon frequencies obtained using the finite-displacement method depends on the range of the real-space interatomic force constants (IFCs) evaluated with the chosen supercell expansion. Taking into account the symmetry of individual structures in the alloy model, where present, between 8 and 192 accurate calculations were required to evaluate the IFCs for each structure, generally representing a considerably larger computational workload than the initial geometry optimisations. For this reason, the force constants were calculated using the 32-atom alloy supercell structures, i.e. without further expansion.
Phonon density of states (DoS) and atom-projected DoS (PDoS) curves were evaluated using Fourier interpolation to obtain frequencies on uniform -centered grids of phonon wavevectors with 24 16 24 Γ × × subdivisions, with a Gaussian smearing of width = 0.032 THz (FWHM 2.5 cm -1 ). The vibrational contributions t o the thermodynamic free energy were evaluated using denser sampling meshes with 36 24 36 × × subdivisions. Phonon band dispersions were evaluated using the band-unfolding scheme implemented in Phonopy [36] to project the vibrational modes calculated for the alloy supercells back onto the parent Pnma structure. The dispersion curves were evaluated across a series of segments passing through all the highsymmetry wavevectors in the Pnma Brillouin zone. Calculations of the lattice thermal conductivity of the SnS and SnSe endpoints, and evaluation of the two-phonon joint density of states of the alloy models, were performed using the Phono3py package. [37]  A c c e p t e d M a n u s c r i p t -Page 4 -All calculations were performed using pseudopotential plane-wave density-functional theory (DFT) as implemented in the Vienna Ab initio Simulation Package (VASP) code [38]. As in our previous study [26], a planewave cutoff of 600 eV was employed alongside Monhkhorst-Pack k-point sampling grids [39] with 4 4 4 × × subdivisions. The PBEsol exchange-correlation functional with the DFT-D3 dispersion correction was used with projector augmented-wave (PAW) pseudopotentials [40,41] including the Sn 5s, 5p and 4d and the outermost 3d/3p and 4s/4p electrons on S and Se in the valence shell. The precision of the charge-density grids was set automatically to avoid aliasing errors, and a support grid with 8 as many points was used to evaluate the × forces. The PAW projection was performed in reciprocal space and nonspherical contributions to the gradient corrections inside the PAW spheres were accounted for.

a. Mixing thermodynamics
The solid solution model for a given composition comprises a set of optimised structures with Se total energy and degeneracy , which are used to form the thermodynamic partition function at a temperature , , according to [28]: where is the Boltzmann constant. is then used to obtain the constant-volume (Helmholtz) free

From
, the free energy of mixing can be obtained as: where and are the free energies of the SnS and SnSe endmembers respectively, which, since , are simply equal to the corresponding . = 1 For a given composition and temperature the occurrence probability of individual structures in the model are calculated from the partition function as: A c c e p t e d M a n u s c r i p t -Page 5 - The can then be used to calculate the thermodynamic average of a general property using: where are the properties of the individual structures in the model.
This model can be improved upon by using the phonon frequencies obtained from the lattice-dynamics calculations to calculate for each structure the temperature-dependent Helmholtz energy including contributions from vibrations: where is the vibrational Helmholtz energy and and are the vibrational internal energy and entropy.

vib vib vib
The superscript is used for clarity to explicitly identify quantities calculated based on the lattice dynamics. vib As for the free energy of the solid solution, is calculable from the vibrational partition function vib ( ) vib ( ) using the bridge relation: Here the phonon frequencies are indexed by a phonon wavevector and band index and is the reduced ℏ Planck constant. The product in Eq. 8 runs over all bands on a grid of wavevectors sampling the phonon Brillouin zone, which appears as a normalisation constant in Eq. 9. The expansion in Eq. 9 shows that is a vib sum of the harmonic zero-point energies and a temperature-dependent term from the population of 1 2 ℏ higher mode energy levels at finite temperature.
We note that this model assumes the distribution of the different supercell configurations is in thermodynamic equilibrium at the formation temperature [28], and that the chosen alloy supercell is F sufficient to sample the range of local configurations that may be present. In real alloy systems, there may be medium-or long-range deviations in local structure and/or composition that are not probed by this method, for example due to the flexible oxidation state and lone-pair activity of Sn [42] .  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 A c c e p t e d M a n u s c r i p t -Page 6 - obtained using the lattice energies (a) and the Helmholtz free energies (b) including contributions from the lattice dynamics. Each plot shows for a range of alloy formation temperatures from 300 -1100 mix ( Se ) F K, which are denoted by line colours from blue (low ) to orange (high ).

F F
By substituting the in Eq. 2 with (Eq. 7), we can build the partition function for the solid solution and evaluate the mixing free energies including the effect of the lattice dynamics ( Fig. 1). Interestingly, we find that including the vibrational free energy terms destabilises the solid solution, decreasing for the 50/50 mix composition ( = 0.5) by a factor of two from -1.72 to -0.85 kJ mol -1 per atom at 900 K compared to using the Se lattice energies. This difference corresponds to 0.1 , which is in the typical range of vibrational entropy ~× B differences observed by van de Walle and Ceder [43]. Se Se much larger difference of 0.12 kJ mol -1 atom -1 is obtained when including the free energy, indicating a larger deviation from ideality, which is clearly evident when comparing the two sets of mixing profiles in Fig. 1.
Of particular note are the large deviations from the general trend in Fig. 1a at = 0.125 and 0.375, Se the reasons for which are not clear. One possible explanation is that the 32-atom supercell does not adequately capture the full range of local structures present in the bulk alloy, although in our previous study we found that this was at least sufficient to converge the mixing energies calculated from the [26]. Another possibility is that the deviation is an artefact of the limited real-space range of the calculated force constants and its impact on the interpolation of the phonon frequencies to a dense grid of wavevectors to evaluate . Although we vib would naively expect a similar level of accuracy across the nine compositions, we found that the proportion of imaginary frequencies in the wavevector sampling meshes -which are excluded from the partition function in Eq. 8 -ranged from 2.06 0.91 % for the = 0.875 composition to 2.81 0.6 % for = 0.375, so this is a ± Se ± Se possibility. As an additional test we also recalculated the mixing energy based on the obtained with a denser   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 A c c e p t e d M a n u s c r i p t -Page 7 - In both cases, the mixing is still energetically favoured at all compositions and the qualitative stability of the alloy is thus unaffected. However, the substantial difference in the mixing energy may be important in systems with competing phases. In our previous work [26,42], we predicted based on lattice energies that the energy difference between competing Pnma and rocksalt Sn(S 1-x Se x ) monochalcogenide solid solutions varies from 10.5 and 2.6 kJ mol -1 per F.U. between the SnS and SnSe endpoints [26]. We also found that the vibrational free energy stabilises the rocksalt phase of SnS relative to the Pnma phase at higher temperatures, reducing the energy difference to 4.5 kJ mol -1 per F.U. at 900 K [42]. It is therefore conceivable that differences in the free energy may bring the rocksalt phase closer to the convex hull at Se-rich compositions. To test this would require an additional set of lattice dynamics calculations on the rocksalt solid solution models, which we have not attempted here.
Including the vibrational free energy in the partition function also leads to a marked change in the distribution of the occurrence probabilities of individual structures. For the 50/50 composition ( = 0.5) and a Se formation temperature of 900 K, the calculated using the lattice energies range from 8.8 10 -5 to 4.2 F × 10 -3 , whereas including the effect of lattice dynamics produces a significantly wider spread of 2. 3 10 -5 -× × 5.2 10 -2 . As shown in Fig. 2, the discrepancy is such that 20 % of structures with the highest account for > × 70 % of the partition function formed using the free energies, compared to just under 30 % with the lattice energies.  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 A c c e p t e d M a n u s c r i p t -Page 8 - The vibrational free energy is derived from the phonon frequencies, which in this system are much more sensitive to the structure than the lattice energy. Again taking the = 0.5 model as an example, the largest Se difference in between two structures is 3.8 meV atom -1 (0.36 kJ mol -1 atom -1 ), whereas the difference in the zero point energies is 0.06 kJ mol -1 atom -1 and the difference in at 900 K is 2.06 kJ mol -1 atom -1 . Given the vib number of examples in the literature showing that differences in vibrational free energy -in particular the vibrational entropy -are a key driver of temperature-induced phase transitions [42,[44][45][46] and are important in the determining the stability of alloy systems [43,47,48] -this finding is not surprising, but is nonetheless noteworthy.
In principle, the large change in may lead to a significant change in the averaged properties, although this is likely to be mitigated if the variation in properties among the structures is small. To test this, we took the electronic bandgaps of the Pnma structures obtained in our previous study and calculated the averaged bandgap as a function of composition using the 900 K occurrence probabilities calculated using the and Finally, we note that in this analysis we have used the Helmholtz rather than the Gibbs free energy as in our previous study. In the absence of contributions to the free energy from the lattice dynamics, the effect of finite pressure can be modelled by replacing the internal energies with the enthalpies = [49]. This may be done approximately by adding a correction to each or by relaxing the structures + under applied pressure to evaluate the explicitly. However, thermal expansion at finite temperature has a considerable effect on the phonon spectra and the vibrational contributions to the free energy, which should be accounted for by calculating for each structure within the quasi-harmonic approximation (QHA) [35]. However, to do so for all of the structures in the solid solution model would require an order of magnitude more calculations and is therefore not feasible using first-principles techniques.  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 A c c e p t e d M a n u s c r i p t -Page 9 -

b. Phonon spectra
For a given composition and formation temperature , a thermodynamically averaged phonon density F of states (DoS) can be obtained by summing the DoS curves of individual structures weighted by the ( ) corresponding occurrence probabilities (Eq. 5/Eq. 6). Obtaining an averaged phonon dispersion is more ( F ) involved, as the vibrational modes calculated in the alloy supercells must first be "unfolded" onto the parent structure. In the present study we employed the method proposed by Allen et al. [36], which is implemented in recent versions of the Phonopy lattice-dynamics code [35] and accessible through the Python API. Reference structures for the unfolding were generated for each composition by linear interpolation of the lattice vectors and atomic positions between the SnS and SnSe structures.   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 A c c e p t e d M a n u s c r i p t -Page 10 -  Figs. 4c-4e were averaged assuming a 900 K formation temperature (there is only one unique structure for each endpoint, so no averaging is required to generate the spectra in Figs. 4a/4b). A full set of phonon spectra for the nine compositions can be found in Figs. S2-S10.
The spectra of the endpoints each comprise two bands of modes primarily associated with motion of the Sn and chalcogen atoms. At intermediate compositions, the structure of the low-frequency Sn peak is largely unchanged, while the higher-frequency chalcogen feature splits into a mid-frequency Se and high-frequency S component. To a large extent, the Sn sublattice in the mixed phases retains the parent Pnma structure, and this is reflected in the low-frequency part of the unfolded dispersions retaining significant detail. On the other hand, substituting the chalcogen atoms spreads out the mid-and high-frequency parts of the dispersion, leading to almost continuous bands in the 50/50 composition with = 0.5.

Se
The broadening of the frequency spectrum at intermediate compositions may have significant implications for the thermal transport, as it could both decrease the mode group velocities and increase the density of energy-conserving phonon scattering events and thereby suppress the mode lifetimes. These points will be returned to in the following subsection.
It can also be seen from the unfolded dispersions and DoS curves that some of the structures in the mixed compositions display imaginary modes. As noted in the previous section, the force constants were calculated to a limited real-space range, and the unfolded dispersions and DoS curves are derived by interpolating the zone-centre ( -point) modes of the alloy supercells. None of the structures in the solid-Γ solution model show imaginary modes at the zone centre, which implies that these imaginary modes may indeed be interpolation artefacts. However, evaluating the IFCs in an expanded cell would not be practical given the number of structures in the solid-solution model.
Finally, we also examined the spread in the averaged phonon DoS curves by overlaying the weighted standard deviations (Figs. S11-S19). As suggested by the dispersions of the intermediate compositions in Fig. 4, the calculations predict the spread in the DoS, in particular the chalcogen bands, to increase substantially towards the 50/50 composition. We also generated a series of DoS meshes using the larger 36 24 36 × × sampling mesh used to evaluate the vibrational free energies and the tetrahedron method for Brillouin zone integration, which yielded similar overall shapes and predicted spreads (Figs. S20-S28).

c. Implications for thermal transport
Within the single-mode relaxation time approximation (RTA) solution to the Boltzmann transport equation, the macroscopic lattice thermal conductivity tensor can be calculated as a sum of macroscopic latt contributions from individual phonon modes according to [37]: where are the modal heat capacities, are the group velocities, are the mode lifetimes, and the  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 A c c e p t e d M a n u s c r i p t -Page 11 -to integrate over the Brillouin zone. The product is the phonon mean-free path , which appears in a λ frequently-used variant of Eq. 10. and are calculated within the harmonic approximation as: where are the mode frequencies, equivalent to the that appear in Eqs. 8 and 9, are the corresponding mode eigenvectors, and is the dynamical matrix for the phonon wavevector . are the phonon lifetimes ( ) and are given by the inverse of the phonon linewidths : Γ are calculated as the imaginary part of the phonon self-energy, which can be computed perturbatively to Γ third order from: are the phonon occupation numbers from the Bose-Einstein distribution: are the three-phonon interaction strengths calculated from: Φ ′ ′′ Page 11 of 23 AUTHOR SUBMITTED MANUSCRIPT -JPENERGY-100181.R1 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 A c c e p t e d M a n u s c r i p t -Page 12 - The indices label atoms with mass amd , and label the Cartesian directions. The third-order IFCs Φ are calculated using the finite-displacement method: The indices label crystallographic unit cells, and the functions and enforce conservation of energy and ∆ crystal momentum, respectively.
Due to the larger number of pairwise atomic displacements that need to be considered, calculating the third-order IFCs in Eq. 16 to compute the phonon linewidths is typically 1-2 orders of magnitude more computationally expensive than obtaining the second-order force constants for the harmonic phonon calculation, and represents the bulk of the computational workload when modelling within the RTA. latt We have previously modelled the thermal conductivity of Pnma SnS and SnSe using the RTA [6,50], based on which we obtained room-temperature (300 K) averaged values of 0.74 and 1.28 W m -1 K -1 respectively. These are in reasonably good agreement with experimental measurements of 1.2 and 0.4-0.7 W m -1 K -1 [2,51]. The higher thermal conductivity of the selenide predicted by these calculations is possibly unexpected given the general decrease in with the formula mass [37], but given the difficulty of preparing high-quality bulk single latt crystals for measurements it is not clear whether this is a real phenomenon or an issue with the calculations [52,53]. , the modal heat capacities , the group velocity norms and the phonon lifetimes . In both Tr( ) | λ | systems, the span roughly five orders of magnitude. The heat capacities are a shallow function of frequency and vary by 10 -5 eV K -1 over the 10 THz range of the SnS phonon spectrum. The in both systems ~λ vary over two orders of magnitude from approx. 10 2 -10 4 ms -1 and are generally higher in the sulphide than the selenide, reflecting the stronger chemical bonding. The mode lifetimes vary by around two orders of magnitude from 0.1 -10 ps in SnS and 0.5 -50 ps in SnSe. In both compounds, there is a marked difference in the spectrum of for the lower-and higher-frequency modes in the phonon spectrum, with the latter being generally shorter lived.  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 A c c e p t e d M a n u s c r i p t -Page 13 -  [50] and [6] respectively. A c c e p t e d M a n u s c r i p t -Page 14 -  Due to the tensor product in Eq. 10 the are proportional to and the spectrum of group λ ⊗ λ 2 velocities largely determines the overall spectrum of . are directly proportional to the to the lifetimes, and the variation in with frequency thus imposes additional structure. Comparing SnS and SnSe, it can be seen that the balance of smaller but longer in the selenide results in similar predicted overall thermal λ conductivities. The larger averaged thermal conductivity of the SnSe endpoint seems in particular to arise from a group of low-frequency modes ( < 1.5 THz) with group velocities comparable to SnS but with significantly longer lifetimes.
are calculated within the harmonic approximation, making it straightforward to compute an λ averaged frequency spectrum for three of the mixed compositions in the Pnma solid-solution model (Fig. 6). Due to the large number of structures, the of each were generated on a small sampling mesh with 5 3 5 λ × × subdivisions. Tests on the SnS and SnSe endpoints using an equivalent 9 3 10 sampling mesh yielded 300 K × × thermal conductivities to within 5 % of the values obtained with the larger 16 16 16 meshes used in the × × previous studies [6,50]. We note that we are comparing data for the mixed compositions to previous calculations performed without a dispersion correction and using larger supercell expansions to calculate the force constants. In general, the main impact of the exchange-correlation functional on the lattice dynamics is through differences in the predicted equilibrium volume [54]. The optimised volumes of SnS and SnSe obtained in Refs. [6] and [50] are similar to those of the two endpoints in our solid-solution model 46.9/45.8 and 51.8/50.7 Å 3 per F.U.
( respectively), so we do not expect this to be big issue. However, the larger supercell expansion used in the endpoint calculations will in principle make the calculated more accurate by improving the accuracy of the interpolated dispersion and hence the derivatives (c. f. Eq. 12). This may in particular explain why = ∂ ∂ some of the lower-frequency modes in the mixed compositions are predicted to reach higher than in either | | of the endpoints, which is otherwise unexpected.  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 A c c e p t e d M a n u s c r i p t -Page 15 -Due to the significantly higher cost of calculating the third-order IFCs, explicitly computing the lifetimes of even a subset of the solid solution models is impractical. However, it is possible to obtain some qualitative insight into how the alloying may affect the heat transport from the (harmonic) phonon DoS. As noted above, the phonon lifetimes in Eq. 10 are calculated from the phonon linewidths (Eq. 13), which are in turn computed Γ as a sum of contributions from three-phonon scattering processes with pairs of modes and (Eq. 14) that ′ ′′ depend on the three-phonon interaction strength , the mode occupation numbers and the difference |Φ ′ ′′ | 2 in the phonon frequencies. An approximate linewidth can thus be written as the product of an averaged Γ phonon interaction strength and a joint density of states (JDoS) counting the number of energy-2 ( , ) conserving scattering channels as a function of frequency: [37] Γ = 18 The interaction strength is defined as: where is the number of atoms in the primitive cell and is the number of bands at each phonon 3 wavevector. The JDoS is a sum of two functions corresponding to collision (Type 1 -two phonons in, 2 ( , ) one out) and decay processes (Type 2 -one phonon in, two out): and are defined as follows: Page 15 of 23 AUTHOR SUBMITTED MANUSCRIPT -JPENERGY-100181.R1 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 A c c e p t e d M a n u s c r i p t -Page 16 - strengths . As in Fig. 5 the heat maps are colour coded by from red (small ) to yellow (large ). Note that the y-axes in subplots (c) and (d) are on a linear scale while the other subplots use a logarithmic scale. As in Fig. 5 the data for SnS and SnSe was generated from the calculations in Refs. [50] and [6] respectively.
are the phonon occupation numbers at the calculation temperature (Eq. 15) and the functions and ∆ enforce conservation of energy and crystal momentum, respectively (c.f. Eqs. 14/16). Finally, we further define a function and component parts and averaged over wavevectors as follows: Page 16 of 23 AUTHOR SUBMITTED MANUSCRIPT -JPENERGY-100181.R1 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 A c c e p t e d M a n u s c r i p t -Page 17 - The JDoS functions are calculable from the harmonic phonon frequencies and, if the can be assumed to be similar in the intermediate phases and endpoints, could provide an indication of the expected composition dependence of the phonon lifetimes. Figure 7 compares the spectrum of phonon linewidths of the SnS and SnSe endpoints against the JDoS functions and the . The distribution of reflects the trends in the lifetimes shown Γ in Fig. 5: in both spectra, the higher-frequency bands of modes show broader linewidths (shorter lifetimes) than the lower-frequency modes, and the sulphide shows broader linewidths than the selenide. Comparison of the linewidth spectra against shows that variation in the density of energy-conserving scattering pathways 2 ( ) with frequency is largely responsible for the overall structure, with peaks in overlapping with increases 2 ( ) in . Collision (Type 1) and decay (Type 2) events are the dominant scattering processes at low and high Γ frequencies, respectively, while both processes are active at intermediate frequencies.
On the other hand, differences in the magnitude of the phonon interaction strengths are clearly responsible for the marked difference in the linewidths of the low-and high-frequency bands of modes, and the lower in the selenide counteract the increase in the JDoS, resulting in narrower linewidths (longer lifetimes) than in the sulphide. distinct peaks due to the well-defined bands of modes in the respective phonon spectra, the JDoS functions of the mixed compositions are much broader, implying stronger scattering across a broad spectrum of mid-and high-frequency modes. However, the calculations predict the JDoS of all three mixed compositions to be on a similar scale to the endpoints, which suggests that the broader range of phonon frequencies in the mixed compositions does not necessarily lead to as large an enhancement in the number of scattering channels as might be expected. This is, however, consistent with the relatively small absolute variation in observed latt among different alloy compositions in experiments [12,14].
The analysis in Fig. 5 suggests that the majority of the heat transport is through low-frequency modes with large group velocities and long lifetimes. The calculated JDoS functions for the mixed phases suggest an overall increase in the density of both collision and decay pathways for modes up to 2 THz compared to the ẽ ndpoints. For SnS this can be ascribed to chalcogen substitution introducing low-lying optic modes to participate in collision processes, whereas for SnSe the dominant effect is most likely the broadening of the mid-frequency Se phonon bands creating additional scattering channels. In both cases, the small increase in the number of scattering channels would be expected to result in some suppression of heat transport in the alloy, consistent with experimental findings. However, based on the analysis in Fig. 7 we would expect differences in the interaction strengths with composition to have a prominent effect on the heat transport and, as it is not feasible to investigate this further, we therefore treat this finding with caution.  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 A c c e p t e d M a n u s c r i p t -Page 18 - (Type 2 -) events separately. The data in Fig. 7 for the SnS and SnSe endpoints is included for comparison (2) 2 and was generated using the previous calculations in Refs. [6] and [50].  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 A c c e p t e d M a n u s c r i p t -Page 19 -

Conclusions
This work has presented a detailed investigation of the lattice dynamics in Pnma Sn(S 1-x Se x ) solid solutions and the impact of the dynamics on the energetics and thermal transport.
Our calculations show that the phonon free energy makes a significant contribution to the thermodynamic partition function, in this case leading to a reduction in the mixing energy and a marked difference in the distribution of occurrence probabilities among individual configurations in the alloy model. Here the qualitative stability of the mixed phases remains unchanged, and the similar electrical properties of the individual structures in the mixed compositions results in negligible changes to the averaged bandgaps. However, given the size of the effects predicted by the calculations, it is very likely that the lattice dynamics could have a significant qualitative impact in other systems, particularly those with energetically-similar competing phases, as has been observed in other work [43,47]. In the case of the Sn n (S 1-x Se x ) m , it is possible that differences in lattice dynamics may bring the metastable rocksalt phase closer to the Sn(S 1-x Se x ) convex hull at intermediate compositions and/or stabilise the sesquisulphide Sn 2 (S 1-x Se x ) 3 relative to the competing mono-and di-chalcogenide phases, but we defer exploration of this to future work.
Comparing the unfolded phonon dispersions and density of states over the range of compositions shows a clear pattern whereby the well-defined bands of high-and mid-frequency S/Se modes in the SnS and SnSe endpoints, respectively, are spread over a wider range of frequencies while the band structure of the lowfrequency Sn modes is largely retained. This results in some lowering of the group velocities and produces additional energy-conserving scattering channels in the joint density of states, both of which could in principle reduce the thermal conductivity in the mixed phases. However, comparison of the endpoints suggests that differences in the three-phonon interaction strengths make a significant contribution to the difference in heat transport between them, and it is not practical at present to investigate how these are affected by alloying.
This leads us to close by noting several challenges to this type of study that should be addressed in future work. The practical restriction of calculating the harmonic force constants in the 32-atom alloy supercells is likely restrict the accuracy of calculated phonon spectra and derived properties such as the free energies and group velocities. In general, this issue may be further exacerbated by potentially needing to reduce the planewave cutoff and k-point sampling in alloy calculations to minimise the computational cost. Secondly, more accurate free energies could be obtained by taking into account the volume expansion at finite temperature using the quasi-harmonic approximation, but to perform QHA calculations for the full set of structures in a complex alloy such as Pnma SnS is not currently practical using first-principles methods. Finally, as noted above, our analysis of the thermal conductivity of the SnS and SnSe endpoints indicates that differences in the phonon interaction strengths are likely to be a significant contributor to the differences between the compounds, but to quantify this requires computing the third-order force constants which is again not practical even for a small subset of the structures.
Several routes to address these challenges can be envisaged. The cost of the first-principles calculations could be reduced by using alternatives to plane-wave DFT, for example periodic codes using efficient localorbital or LCAO approaches [55][56][57]. Alternatively and/or in addition to this, a sampling approach could be employed to obtain the second-and/or third-order force constants using a smaller set of explicit calculations. The relatively predictable set of chemical environments in an alloy should also make them amenable to modelling with force fields, which could be parameterised against a smaller subset of first-principles calculations in the spirit of the aforementioned sampling approach.