Structural Dynamics and Thermal Transport in Bismuth Chalcogenide Alloys
We present a detailed study of the structural dynamics, energetic and dynamical stability and thermal transport of the bismuth chalcogenides Bi2S3, Bi2Se3 and Bi2Te3 and their alloys. The differing activities of the Bi cation lone pairs in Bi2S3, Bi2Se3 and Bi2Te3 lead to competition between orthorhombic Pnma and rhombohedral R3 ̅ m phases, with the latter favored by heavier chalcogen atoms, while the reported non-ambient phases of Bi2Se3 and Bi2Te3 show phonon instabilities under ambient conditions. The Pnma structure has weaker chemical bonding and stronger phonon anharmonicity than the R3 ̅ m phase, resulting in an intrinsically lower lattice thermal conductivity. A thermodynamic model of Bi2(Se1-xSx)3 indicates that the R3 ̅ m structure is energetically favored only for low S content, but the stability window can potentially be extended at lower formation temperatures. Bi2(Se1-xTex)3 in the R3 ̅ m phase shows substantial deviation from ideal solid-solution behavior, due to a strong energetic preference for the Se and Te atoms to occupy the interior and exterior sites, respectively, in the constituent quintuple layers. Strain-field fluctuations induced by inhomogeneities in the chemical bonding are shown to play a significant role in determining the heat transport in the alloy systems, and chalcogen disorder away from the preferred symmetric structure is found to be an important factor in the reduced thermal conductivity of Bi2SeTe2 compared to the Bi2Te3 endpoint. The microscopic insight from these modelling studies provides new insight into the bismuth chalcogenides and their alloys, which may inform ongoing research to optimize the thermoelectric performance of these and related materials.
INTRODUCTION<!>𝑍𝑇 =<!>METHODS<!>RESULTS AND DISCUSSION<!>Thermodynamics of Binary Solid Solutions. Based on the relative energies of the Pnma<!>𝑛=1
<p>Rising concerns over climate change and dwindling non-renewable resources have led to increased emphasis on developing technologies to improve energy sustainability and to meet growing energy demand. Since more than 60 % of the energy used worldwide is wasted as heat, 1 thermoelectric generators (TEGs), which can extract electrical energy from a temperature gradient, are currently the subject of intensive research. 2 TEGs are easily scalable to a range of sizes for applications including in the automotive, manufacturing and nuclear industries, and can thus potentially improve energy sustainability across a wide range of sectors. 2 The performance of a thermoelectric material in a TEG is typically expressed by the dimensionless figure of merit 𝑍𝑇: 3</p><!><p>𝑆 2 𝜎𝑇 𝜅 latt + 𝜅 el (Eq. 1)</p><p>The Seebeck coefficient 𝑆, electrical conductivity 𝜎 and electronic thermal conductivity 𝜅 el depend on the electronic structure of the material and can be fine-tuned by chemical doping. All three are interdependent through the carrier concentration 𝑛 and must be balanced to optimize the 𝑍𝑇 at a target operating temperature. The lattice thermal conductivity 𝜅 latt is the transport through phonons, which depends on the structure and bonding in the material and can in principle be optimized independently of the electrical properties. [3][4][5] The requirement for a favorable electronic structure and low 𝜅 latt has led to a number of high-performance TEs based on "heavy" chalcogenides with strongly-anharmonic lattice dynamics and intrinsically low thermal conductivity. PbTe is generally considered as a standard for hightemperature thermoelectric applications and can reach a 𝑍𝑇 of 2.2 at 915 K with hierarchical nanostructuring. 6 Its high performance has been attributed to anharmonic lattice dynamics 7 and to a convergence of the electronic bands at elevated temperature leading to enhanced band degeneracy and a larger Seebeck coefficient. 8,9 More recently, SnSe was shown to have a very high bulk 𝑍𝑇 of up to 2.6 at 923 K, due to an ultra-low 𝜅 latt associated with the high-temperature phase transition. [10][11][12][13][14] With a 𝑍𝑇 of ~1.0 between 300-400 K, Bi2Te3 is the commercial standard for room-temperature TEGs at present, due to its favorable electronic properties and inherently low 𝜅 latt , 2,15 the latter of which has been attributed to weak chemical bonding and highly anharmonic lattice dynamics. 15,16 However, the low terrestrial abundance of Te (ca. 1 ppb) limits large-scale adoption. 2 Therefore, there has been substantial investigation into variants of Bi2Te3 in the same tetradymite family 17 (e.g. Bi2Se2Te) and using strategies such as alloying to reduce or eliminate Te. [18][19][20][21][22] Doping and alloying are common design strategies for enhancing the 𝑍𝑇 of thermoelectrics, as they can both improve the Seebeck coefficient by promoting band convergence and reduce the 𝜅 latt by introducing point defects that scatter phonons through mass-and strain-field fluctuations. 3 A number of recent studies on bismuth chalcogenide alloy systems, including Bi2(S,Se)3, Bi2(S,Te)3 and Bi2(Se,Te)3, have reported a 𝑍𝑇 of 0.5-1.3 (Table 1). [19][20][21][22][23][24] Table 1 Reported thermoelectric figures of merit 𝑍𝑇 of some bismuth chalcogenide alloys.</p><p>𝑍𝑇 𝑇 [K] a Bi2(STe2) 19 0.8 573 Bi2(S2Se) 19 0.8 573 Bi2(Se0.3Te2.7) 20,23 0.54/1.27 RT Bi2(Se0.5Te2.5) 21 1.28 RT Bi2(Se0.6Te2.4) 24 0.79 RT Bi2(Se1.5Te1.5) 22 0.67 473 a RT -room temperature.</p><p>Liu et al. recently carried out a systematic study of the Bi2(S,Se,Te)3 ternary phase space and established several compositional trends in the thermoelectric performance. 19 The thermal conductivities of mixed compositions along the Bi2(Se1-xTex)3 tie line were found to be lower than the closest endpoints and to reach a minimum for Bi2Se2Te. Bi2STe2 and Bi2S2Te were found to show lower thermal conductivities than either of Bi2S3 and Bi2Te3, and a global minimum in the 𝜅 latt was identified at Bi2Te1.2Se1.8. The reduced 𝜅 latt of the mixed compositions was attributed to enhanced phonon scattering in the alloys. In contrast, measurements of the electrical properties showed a strong dependence on the Te content, and the weighted mobility was optimized towards Bi2Te3.</p><p>Thermoelectrics are attractive candidates for modelling studies because the terms in Eq. 1</p><p>are amenable to calculation. The electrical properties can be predicted using semi-classical Boltzmann transport theory, 25,26 and it has also recently become possible to use first-principles modelling to investigate the structural dynamics and lattice thermal conductivity. [27][28][29][30] The latter calculations have played an important role in understanding the intrinsic anharmonicity in materials with low 𝜅 latt , 11,14,[31][32][33][34] including several flagship thermoelectrics, 11,12,14,33 and have also been used to identify and characterize the thermal transport in novel candidate TEs. 34 There are several reports in the literature of modelling studies on the thermoelectric properties of bulk bismuth chalcogenides [35][36][37][38][39] and nanostructures (e.g. single layers), [40][41][42][43] including a handful of studies that investigate the structural dynamics and thermal conductivity. 36,[40][41][42][43] Modelling studies on the alloys are more scarce, owing to the higher complexity of the calculations, and tend to be restricted to compositions with well-defined structures. 37,41,42 Previous theoretical studies of alloy systems using first-principles density functional theory (DFT) and statistical thermodynamics have demonstrated the capacity systematically to explore a range of compositions between the endpoints of a two-component alloy system. 44,45 The structure, properties and energetic and dynamical stability of different compositions can be evaluated in silico, providing a means to predict the electronic structure and thermal transport of the alloy and to guide experimental design towards optimizing the properties for thermoelectric applications. [44][45][46] The high computational resource requirements for modelling the lattice dynamics and thermal transport, however, pose significant challenges to studying alloy systems this way, which has resulted in few investigations in this direction. 42,[44][45][46][47][48] In this work, we use first-principles DFT calculations to perform a comprehensive study of the lattice dynamics, energetic and dynamical stability and thermal transport of all of the reported phases of Bi2S3, Bi2Se3 and Bi2Te3, and we employ methodology developed in previous studies to investigate the thermodynamics and properties of Bi2(Se,S)3 and Bi2(Se,Te)3 alloys. The differing activities of the Bi cation lone pairs in Bi2S3, Bi2Se3 and Bi2Te3 lead to competition between R3 ̅ m and Pnma phases with very different thermal transport properties, which we attribute quantitatively to differences in the phonon group velocities and lifetimes. Whereas the Bi2(Se,S)3 alloy in the Pnma phase forms a homogenous solid solution, the Bi2(Se,S)3 and Bi2(Se,Te)3 alloys in the R3 ̅ m phase show strong structural preferences for certain chalcogen arrangements. This leads to significant deviation from ideal mixing behavior, which can also be rationalized in terms of electrostatics. Finally, we use a variety of models to investigate the thermal conductivity in the alloys, and we show that the strain field introduced by inhomogeneous chemical bonding is likely to have a significant impact on the transport. The detailed microscopic insight from these studies highlights the benefits of this theoretical modelling approach for studying alloy systems, providing the potential for its extension to a wider range of solid-state materials for thermoelectric applications.</p><!><p>First-principles calculations were carried out using pseudopotential plane-wave DFT as implemented in the Vienna ab initio Simulation Package (VASP) code. 49,50 Initial structures of the Pnma, R3 ̅ m, and P42/nmc phases of the bismuth chalcogenides were obtained from the Materials Project 51 database and fully optimized using a plane-wave cutoff of 600 eV and the Γ-centred</p><p>Monkhorst-Pack 𝐤-point sampling meshes 52 listed in Table 2. These parameters were determined to converge the total energy to < 1 meV atom -1 and the pressure to < 1 kbar (0.1 GPa). The Materials Project has two entries for R3 ̅ m Bi2Te3, the second of which has a different ordering of the atoms within the quintuple layers, and which we denote R3 ̅ m*. The PBEsol GGA functional 53 with the DFT-D3 dispersion correction 54 (i.e. PBEsol+D3) was used to model exchange and correlation effects. PAW pseudopotentials 50,55 were used to model the core electrons, with the valence electrons of S, Se, and Te and the valence and semi-core d electrons of Bi treated as valence states. Tolerances of 10 -8 eV and 10 -2 eV Å -1 were applied to total energy and forces during wavefunction and geometry optimisation, respectively.  The size of the charge-density grids was set automatically to avoid aliasing errors, the PAW projections were performed in reciprocal space, and non-spherical contributions to the gradient corrections inside the PAW spheres were accounted for.</p><p>Lattice dynamics calculations were performed on the optimized structures to obtain the 2 ndand 3 rd -order interatomic force constants (IFCs) using the supercell finite-displacement approach implemented in the Phonopy and Phono3py packages. 27,28 The IFCs were determined using the supercells shown in Table 2 and displacement distances of 10 -2 and 3 × 10 -2 Å for the harmonic and anharmonic calculations respectively. Force calculations were performed using VASP. The 𝐤point meshes in Table 2 were reduced proportionally for the chosen supercells, and additional charge-density grids with 8 × the number points as the standard grids were used to ensure accurate forces. Phonon density of states (DoS) curves were obtained using Fourier interpolation to evaluate frequencies on uniform 48 ×48 × 48 Γ-centered 𝐪-point meshes and the liner tetrahedron method for Brillouin Zone integration. The phonon contributions to the constant-volume (Helmholtz) free energy were also evaluated using these sampling meshes. Phonon dispersions were computed by interpolating phonon frequencies along strings of high symmetry 𝐪-points, with the band paths determined using the SeeK-path code. 56 Thermal conductivity calculations were performed within the single-mode relaxation-time approximation (RTA). 28 The modal properties were sampled on the regular Γ-centered 𝐪-point meshes in Table 2. These were sufficient to converge the calculated lattice thermal conductivities to within 5 %, except for Pnma Bi2Se3, which due to erratic convergence was converged to within ~10 %. Finally, we also performed phonon and thermal-conductivity calculations on two selected structures in the Bi2(Se1-xTex) alloy, with the 2 nd -order and 3 rd -order IFCs determined in supercells matching those used for the endpoint structures.</p><!><p>Chemical and Lattice Stability of the Bismuth Chalcogenides. we identified four reported crystal phases across the three chalcogenides Bi2S3, Bi2Se3 and Bi2Te3 (Figure 1). Bi2S3 is reported in an orthorhombic Pnma phase with 20 atoms in the primitive cell. 19 Bi2Se3 is reported in an equivalent Pnma phase and rhombohedral R3 ̅ m and tetragonal P42/nmc phases with five and 40 atoms in the respective primitive cells. 19,[58][59][60] Finally, two R3 ̅ m phases of Bi2Te3 are reported, both with five atoms in the primitive cell; 19,61 the first of these is isostructural with the R3 ̅ m phase of Bi2Se3, while the second, which we denote R3 ̅ m*, has different ordering of atoms within the quintuple layers along the crystallographic c axis. S, Se and Te atoms are colored yellow, green and olive, respectively. These images were prepared using the VESTA software. 62 Since the Pnma and R3 ̅ m phases are common to more than one chalcogenide, we also performed calculations on Bi2Te3 in the Pnma structure and Bi2S3 in the R3 ̅ m structure. In total, we therefore investigated eight structures: two of Bi2S3, and three each of Bi2Se3 and Bi2Te3.</p><p>The optimized lattice parameters and relative lattice energies ∆𝑈 latt are shown in Table 3.</p><p>The lattice parameters are in good agreement with available measurements, 19,[59][60][61] suggesting that PBEsol+D3 provides a good description of the structures. With the exception of P42/nmc Bi2Se3 and R3 ̅ m* Bi2Te3, the calculated and experimental values are within 4 %. For these two structures, a larger deviation of ~10% is observed, which can be explained by the fact that these phases are prepared under high-pressure and high-temperature (i.e. non-ambient) synthesis conditions. 60,61 For Bi2S3, the calculations predict the Pnma phase to be more stable than the R3 ̅ m phase by 2.9 kJ mol -1 atom -1 . For Bi2Se3, the stability is reversed, with the R3 ̅ m being more stable than the Pnma phase by 1.8 kJ mol -1 atom -1 . The P42/nmc phase is predicted to be much higher in energy than either of these, which we again ascribe to it being prepared under non-ambient conditions. The calculations again predict the R3 ̅ m phase of Bi2Te3 to be more stable than the Pnma phase, this time by a larger energy difference of 5.6 kJ mol -1 atom -1 . This is small enough that the Pnma phase may be accessible, but there are no experimental reports of its synthesis. As for the non-ambient P42/nmc phase of Bi2Se3, the R3 ̅ m* phase of Bi2Te3 is predicted to be considerably higher in energy than both the R3 ̅ m and Pnma phases. This can be understood by the different atomic ordering in the quintuple layers: in the R3 ̅ m phase, layers of Bi 3+ are sandwiched by oppositelycharged layers of Te 2layers, which is electrostatically favorable, whereas the R3 ̅ m* phase has adjacent layers of Te atoms that depending on the electron density distribution may result in less favorable electrostatic interactions and/or some level of electrostatic repulsion.</p><p>The energetic preference of the heavier chalcogenides to adopt the R3 ̅ m structure can be rationalized by the stereochemical activity of the Bi 6s lone pair. Bi 3+ experiences an inert-pair effect due to the relativistic contraction of the 6s orbital. Interactions between the Bi 6s and anion p orbitals produces a high-energy antibonding orbitals that can interact with the Bi 6p, thus indirectly mediating Bi 6s/6p hybridization. 63 Since this s/p mixing is forbidden by symmetry in an octahedral bonding environment, a structural distortion is required that results in the lone pair occupying a coordination site. The distortion is only favored when the stabilization due to the anion-mediated s/p mixing offsets the loss of electrostatic stabilization due to the reduced coordination number. The orbital interaction between the Bi 6s and S 3p orbitals is the strongest among the three bismuth chalcogenides, as the relative energy difference is the smallest. Bi2S3 therefore preferentially adopts the low-symmetry orthorhombic Pnma phase with an active lone pair, and remains stable even under pressure. 64,65 Table 3 Optimized lattice parameters, relative lattice energies ∆𝑈 latt and dynamical stabilities of the eight bismuth chalcogenide phases considered in this work.  19,[59][60][61] shown in parentheses for comparison. b Lattice energies reported relative to the most stable phase of each chalcogenide.</p><p>On the other hand, the poorer energy match between the Bi 6s and Se 4p/Te 5p leads to weaker orbital interactions that are insufficient to stabilize an active lone pair, resulting in "quenching". 66 Bi2Se3 and Bi2Te3 therefore preferentially adopt the higher-symmetry R3 ̅ m phase. relatively high energies, we therefore conclude that these phases are highly unlikely to be formed under ambient conditions. The dynamical stabilities of the eight structures, assessed from the presence or absence of imaginary modes in the calculated phonon dispersions, are collected alongside the relative energies in Table 3.</p><p>We can further use the phonon calculations to compute the vibrational contributions to the temperature-dependent constant-volume (Helmholtz) free energy 𝐴: where 𝐴 vib = 𝑈 vib (𝑇) − 𝑇𝑆 vib (𝑇) is the contribution to the free energy from the phonon modes and 𝑈 vib and 𝑆 vib are the phonon internal energy and entropy respectively. When two phases are close in energy, differences in phonon spectra and 𝐴 vib can lead to qualitative changes in stability at 0 K and/or to temperature-induced phase transitions. 67 Bi2S3 and the R3 ̅ m phase for Bi2Se3 and Bi2Te3 (c.f. Table 3).</p><p>Figure 3 compares the energy differences ∆𝐴 between the Pnma and R3 ̅ m phases of Bi2S3, Bi2Se3 and Bi2Te3. In all three, we find that the Pnma phase is lowered in energy relative to the competing R3 ̅ m phase with increasing temperature, indicating that the former is stabilized by differences in the vibrational free energy. For Bi2S3, this means that the energy difference between the ground-state Pnma and the R3 ̅ m phase increases sharply with temperature. For Bi2Se3 and Bi2Te3, on the other hand, the free-energy differences between the Pnma and R3 ̅ m phases fall with temperature. For Bi2Se3, the ∆𝐴 falls from ~1.8 kJ mol -1 atom -1 at 0 K to 0.75 kJ mol -1 atom -1 at 1000 K, whereas for Bi2Te3 the differences range from ~5.5 kJ mol -1 atom -1 at 0 K to ~3.9 kJ mol -1 atom -1 at 1000 K. (The melting temperatures of Bi2Se3 and Bi2Te3 are around 983 and 858 K, respectively.) While this analysis further highlights the chemical trend in the energetic preference for the R3 ̅ m structure in the heavier chalcogenides, the calculations do not suggest any thermal phase transitions and predict the same qualitative stability ordering as the lattice energies (c.f.</p><p>Table 3). Nonetheless, the Pnma phase of Bi2Se3 is close in energy to the R3 ̅ m phase, and it may therefore be possible to stabilize it by alloying with Bi2S3 or by epitaxial growth on a suitable substrate. We return to the former point in a following subsection.</p><p>Lattice Thermal Conductivity of Pnma and R𝟑 ̅ m Bi2S3, Bi2Se3 and Bi2Te3. Using the single-mode relaxation-time approximation, the lattice thermal conductivity 𝜅 latt is computed as a sum over contributions 𝜿 𝜆 from individual phonon modes 𝜆 according to:</p><p>where 𝐶 𝜆 are the modal heat capacities, 𝒗 𝜆 are the model group velocities, 𝜏 𝜆 are the phonon lifetimes, 𝑁 is the number of wavevectors included in the summation, and 𝑉 0 is the unit cell volume. 𝐶 𝜆 and 𝒗 𝜆 are calculated within the harmonic approximation as:</p><p>𝒗 𝜆 = 𝜕𝜔 𝜆 𝜕𝐪 𝜆 (Eq. 5)</p><p>Where 𝜔 𝜆 are the (angular) phonon frequencies and 𝑘 B is the Boltzmann constant. The 𝜏 𝜆 are determined as the inverse of the phonon linewidths Γ 𝜆 according to:</p><p>The Γ 𝜆 are computed as the imaginary part of the phonon self-energies from a perturbative expression using the harmonic frequencies and eigenvectors and the third-order force constantsfull details of this method are given in Ref. 28 . (𝜅 𝑥𝑥 + 𝜅 𝑦𝑦 + 𝜅 𝑧𝑧 ). We calculate a 𝜅 ave of 0.9 for Pnma Bi2S3, which compares favorably to the experimental value of 0.87 W m -1 K -1 . 68 Our predicted 𝜅 ave for R3 ̅ m Bi2Se3 and Bi2Te3 are 1.82 and 0.88 W m -1 K -1 . The value for Bi2Se3 again compares favorably to the measured value of 1.7 W m -1 K -1 quoted in Ref. 15. On the other hand, the latter significantly underestimates the values of 1.28 and 1.37 W m -1 K -1 measured for Bi2Te3. 15,69 The reasons for this discrepancy are not clear. Bi2Se3/Bi2Te3 (c.f. Table 2), and that in the R3 ̅ m phases 𝜅 𝑥𝑥 = 𝜅 𝑦𝑦 by symmetry.</p><p>Several trends are apparent from the calculations. For a given structure type (i.e. Pnma or R3 ̅ m), the 𝜅 latt decreases with increasing chalcogen mass as expected. More interestingly, for a given chalcogen the average 𝜅 latt of the Pnma structure is consistently smaller than that of the R3 ̅ m phase, by 60-70 %. This is most noteworthy for Bi2Se3, which has been reported in both structure types and for which our calculations predict the Pnma and R3 ̅ m phases to be energetically competitive. Finally, we find that for both structure types the 𝜅 latt is predicted to be anisotropic. In the Pnma phases, the conductivity is largest along the short axes (the a axis for Bi2S3 and the b axes for Bi2Se3 and Bi2Te3; c.f. Table 3). The 𝜅 latt is consistently around 80 % smaller along one of the two orthogonal axes, and is 52, 46 and 37 % smaller along the other axis in Bi2S3, Bi2Se3 and Bi2Te3 respectively.</p><p>The anisotropy in the thermal transport decreases for the heavier chalcogens, but all three Pnma systems retain distinct "easy" and "hard" axes. For the R3 ̅ m, the transport along the x and y directions corresponds to transport in the ab plane, and is 4.4, 2.4 and 1.6 × larger than the transport in the cross-plane z direction (c axis) in Bi2S3, Bi2Se3 and Bi2Te3 respectively. As for the Pnma systems, the anisotropy therefore decreases with increasing chalcogen mass, but even for Bi2Te3 the cross-plane transport is 40 % smaller than the in-plane transport.</p><p>To better understand the origin of the trends in the 𝜅 latt , we use the constant relaxationtime approximation (CRTA) model developed in our previous work. 46,70 We write the 𝜿 latt as the product of a temperature-dependent harmonic function 𝜿 𝜏 CRTA ⁄ and weighted average lifetime 𝜏 CRTA defined according to:</p><p>where the terms are defined as in Eq. 3 above. As for 𝜿 latt the 𝜿 𝜏 CRTA ⁄ are tensors, and the 𝜏 CRTA may be taken as a scalar, like the 𝜏 𝜆 , or as a tensor to provide some flexibility to account for anisotropy in the averaged mode lifetimes. 46,70     phases the lifetimes increase by a factor of two from Bi2S3 to Bi2Se3, but are reduced by ~25 % in Bi2Te3 compared to Bi2Se3. An increase in lifetime counteracting a reduced harmonic term was also observed when comparing the thermal conductivity of SnS and SnSe. 46 This method of analyzing the 𝜿 latt and comparing between materials is relatively new, but we have previously performed similar calculations on SnS/SnSe and CoSb3, 46,70 all of which are widely-studied thermoelectric materials. For SnSe, we calculated a room-temperature 𝜅 ave of 1.58 W m -1 K -1 together with a 𝜅 𝜏 CRTA ⁄ of 0.372 W m -1 K -1 ps -1 and a 𝜏 CRTA of 4.23 ps. 46 The 𝜅 ave falls between the Pnma and R3 ̅ m phases of Bi2Se3. The harmonic term is larger than in both materials, and the higher 𝜅 ave of R3 ̅ m Bi2Se3 is thus due to its longer average lifetimes while the lower thermal conductivity of Pnma Bi2Se3 is due to both the smaller harmonic term and shorter 𝜏 CRTA . For CoSb3, we obtained 𝜅 ave = 9.98 W m -1 K -1 with 𝜅 𝜏 CRTA ⁄ = 0.273 W m -1 K -1 ps -1 and 𝜏 CRTA = 36.6 ps. The harmonic term is smaller than in both sulfides, comparable to the R3 ̅ m phase of Bi2Se3, and larger than both tellurides. However, the lifetime is much larger than all six chalcogenide materials investigated here, by a factor of 6-17 ×, and this is the main reason for its much larger thermal conductivity.</p><!><p>and R3 ̅ m phases of Bi2S2, Bi2Se3 and Bi2Te3 discussed above, we would expect Bi2(Se1-xSx)3 to potentially form in the Pnma or R3 ̅ m structures, whereas we would only expect Bi2(Se1-xTex)3 to form in the R3 ̅ m structure. This is largely borne out in experimental reports, although Bi2Te2S has been reported in the R3 ̅ m structure. 19 We therefore investigated the mixing thermodynamics of these three alloy compositions following the approach adopted in our previous work on the Snn(S1-xSex)m system. 44,71 Systematically substituting the chalcogen atoms in the 20-atom Pnma Bi2S3 primitive cell yields 1,072 unique structures in 13 compositions between the Bi2S3 and Bi2Se3 endpoints, from which we analyzed every second composition (560 structures). Applying the same procedure to a 2 × 2 × 1 supercell expansion of the R3 ̅ m primitive cell yields 378 unique structures in 13 compositions between the Bi2S3/Bi2Se3 and Bi2Se3/Bi2Te3 endpoints (i.e. 756 structures in total), from which we analyzed all 13 compositions.</p><p>For a given composition, the solid-solution model thus comprises 𝑛 𝑥 symmetryinequivalent chalcogen arrangements with total energy 𝐸 𝑛 , which we equate to the lattice internal energy 𝑈 𝑛 latt , and an associated degeneracy 𝑔 𝑛 . These are used to build the thermodynamic partition function 𝑍(𝑥; 𝑇):</p><!><p>(Eq. 8)</p><p>The Helmholtz energy 𝐴(𝑥; 𝑇) can be computed from the 𝑍(𝑥; 𝑇) using the bridge relation from (Eq. 10)</p><p>where 𝐴(𝑥 = 0; 𝑇) and 𝐴(𝑥 = 1; 𝑇) are the mixing Helmholtz energies of the endpoints which, given the 𝑔 𝑛 are unity by construction, are equal to the corresponding 𝑈 𝑛 latt .</p><p>Given a 𝑍(𝑥; 𝑇), the occurrence probabilities 𝑃 𝑛 of the individual structures in each composition can be computed as:</p><p>(Eq. 11)</p><p>The 𝑃 𝑛 can the be used to compute the thermodynamic average 𝑋 ̅ of a physical property 𝑋 according to: 𝑋 ̅ (𝑥; 𝑇) = ∑ 𝑃 𝑛 (𝑥; 𝑇) × 𝑋 𝑛 𝑛 𝑥 𝑛=1 (Eq. 12)</p><p>In particular, the 𝑈 ̅ latt in Eq. 9 can be computed as a weighted average of the 𝑈 𝑛 latt , and subtracted from the Helmholtz energy to obtain the −𝑇𝑆 conf term. One can then separate the mixing energies 𝐴 mix into contributions from internal energy and entropy (𝑈 mix /−𝑇𝑆 mix ). It is worth noting that the structure of the alloy is generally assumed to be set during its formation, so the averaged properties are predicted using 𝑃 𝑛 computed for a set formation temperature 𝑇 F . 71 Consideration of the contributions to the 𝐴 mix from the 𝑈 mix and −𝑇𝑆 mix terms (Figures S9-S11)</p><p>indicates that these compositions are associated with a significantly reduced (i.e. more negative/favorable) enthalpy and a decrease in entropy, both of which are again most pronounced at lower formation temperatures.</p><p>Taken together, these results indicate a strong energetic preference for particular configurations. This can be understood by the fact that there are two distinct chalcogen environments in the R3 ̅ m structure, which we refer to as Chal (1) and Chal (2) . The single Chal (1) position is in the interior of the layers, where the chalcogen atom forms three short and three long bonds to six Bi atoms in a distorted octahedral environment. The two Chal (2) positions are at the edges of the quintuple layers, where the chalcogens form three covalent bonds to Bi atoms and project a lone pair of electrons into the interlayer spacing. It is thus more favorable for smaller and more electronegative anions to occupy the Chal (1) sites, which benefit from stronger electrostatic interactions with the Bi cations, and for larger and more polarizable anions to occupy the Chal (2) sites, which benefit from more space to accommodate sterically-depending lone pairs.</p><p>Occupying the Chal (2) sites with more polarizable anions may also facilitate stronger interlayer dispersion forces. In the Bi2(Se1-xSx)3 alloy, the preference would be to have S in the Chal (1) sites and Se in the Chal (2) sites, giving the favored 𝑥 = 1/3 stoichiometry. In Bi2(Se1-xTex)3, the preference would be for Se in the Chal (1) sites and Te in the Chal (2) sites, favoring the 𝑥 = 2/3 stoichiometry. That the favored compositions appear to have a stronger influence on the thermodynamics of the Bi2(Se1-xTex)3 alloy may be explained by the larger electronegativity difference between Se and Te and the fact that Te is the largest and most polarizable of the three chalcogens. We also note that this structural preference argument would explain why Bi2Te2S is stable in the R3 ̅ m phase. 19 As expected following these arguments, the distribution of 𝑃 𝑛 for structures in the R3 ̅ m alloy models, particularly Bi2(Se1-xTex)3, are heavily skewed towards particular chalcogen arrangements in the alloy supercells (Figure 7). For Bi2Se2S and Bi2SeTe2 at a formation temperature 𝑇 F = 800 K, the most probable chalcogen arrangements have 𝑃 𝑛 = 0.19 and 0.42, respectively, and the three structures (of 21) with the largest 𝑃 𝑛 collectively account for 41 and 76 % of the cumulative probability. Even for the 𝑥 = 1/2 solid solution of Bi2(Se1-xTex)3, where the number of configurations and hence the configurational entropy are a maximum, the most probable four structures with the largest 𝑃 𝑛 (4.8 % of the total) account for 68 % of the cumulative probability. This is in stark contrast to the Pnma Bi2(Se1-xSx)3 alloy, where for the 𝑥 = 1/3, 1/2 and 2/3 compositions the difference between the smallest and largest 𝑃 𝑛 are 6.3 × 10 -2 , 1.9 × 10 -2 and 5.1 × 10 -2 , respectively, indicating a much more homogenous distribution of configurations. It is interesting to note, however, that the calculations do predict a distribution over chalcogen arrangements in both Bi2Se2S and Bi2SeTe2, at least for high-temperature syntheses, which suggests these "line" compounds might contain significant concentrations of defects where chalcogen atoms swap from their ideal crystallographic sites (i.e. antisite defects).  Despite the skewed distribution of the 𝑃 𝑛 , the predicted cell volume as a function of composition 𝑉(𝑥) at 𝑇 F = 800 K (Figure 8) is close to linear in the Se content, as would be expected from Vegard's law. This implies that there is a small variation in the volumes 𝑉 𝑛 of individual structures.</p><p>Based on the mixing energies in Figure 6 we would expect Bi2(Se1-xSx)3 to form stable solid solutions in both the Pnma and R3 ̅ m phases. Since we calculate the absolute free energies of each of the alloy compositions, we can also use our model to predict which of the two are most stable for a given composition and formation temperature (Figure 9). We find that at most temperatures and compositions the Pnma phase is the most energetically stable, with the R3 ̅ m phase favored up to around x = 1/2 (25-50 % S) below ~550 K, above which it is the most stable phase only very close to the SnSe endpoint (x = 0). Experimentally, Bi2(Se1-xSx)3 is reported to form a Pnma solid solution with compositions from Bi2S3 (x = 1) to the Bi40Se32S28 (x = 0.47), 17 Pnma phases with composition Bi2Se2S and Bi2SeS2 (x = 0.33/0.67), 19 and an R3 ̅ m solid solution from Bi2Se3 to Bi40Se56S4 (𝑥 = 0.07). 17 This is remarkably consistent with our predictions for synthesis temperatures above ~600 K.</p><p>Based on our results, alloying Bi2Se3 with a small amount of Bi2S3 would constitute a reliable route to stabilizing the Pnma phase, whereas a low-temperature synthesis route may allow an R3 ̅ m alloy with up to ~50 % S content to be formed. The former is notable as the Pnma phase is predicted to have a significantly lower thermal conductivity than the R3 ̅ m phase, which would be beneficial for the thermoelectric properties. On the other hand, the latter is potentially also interesting as R3 ̅ m Bi2Se3 and Bi2Te3 are both topological insulators. 17  We also note that while alloys with low S content may be unstable with respect to phase separation into Pnma Bi2S3 and R3 ̅ m Bi2Se3 at lower temperature, the very different structures mean this would likely require significant atomic rearrangement, and alloys formed at high temperature and quenched may therefore remain kinetically stable.</p><p>Effect of Alloying on Thermal Transport. Finally, we investigate the thermal transport of the three alloy systems. In principle, one can calculate phonon spectra and thermal conductivities of individual structures in the alloy models and form a weighted average as per Eq.</p><p>12. However, given that large supercells are likely to be required for accurate results, this "brute force" approach is not feasible. 45,46 The reduction in 𝜅 latt due to alloying and atomic substitution can be understood conceptually as arising from mass disorder and fluctuations in the strain field introduced by inhomogeneous chemical bonding. 3 The former contribution can be modelled approximately using an isotope-scattering model, 73 and we previously found that this approach provided a remarkably good prediction of the reduced 𝜅 latt in the Pnma Sn(S0.1875Se0.8125) alloy. 46 A significant effect from strain-field fluctuations would be expected when the endpoints in the alloy show significantly different bond strengths. Since the bond strengths are captured by the force constants used to compute the harmonic phonon spectra and thermal conductivities, we can obtain some insight into this by recalculating the 𝜅 latt of the endpoint structures with the masses adjusted to the alloyed element. Taking Pnma Bi2(Se1-xSx)3 as an example, if setting the mass of S to Se in the Bi2S3 structure gives a good reproduction of the 𝜅 latt of Bi2Se3, and vice versa, it would imply that the chemical bonding (force constants) in both phases are similar, and the primary difference in the 𝜅 latt is due to the different chalcogen masses. In this case, we might then expect an isotopescattering model to give a good reproduction of the thermal conductivity of the alloy. To investigate further, we approximate the 𝜅 latt of the three alloys as a function of composition using the method employed in our previous work on the Li(No,Mn,Co)O2 (NMC) alloy system. 74 We calculate the 𝜅 ave of both endpoint structures as the average mass 𝑚 ave of the chalcogen is varied to the other endpoint:</p><p>where the sum runs over the 𝑛 𝑖 species present at the disordered site with fractional occupation 𝑎 𝑖 and mass 𝑚 𝑖 . We also use the isotope-scattering model in Ref. 73 to model the effect of mass variance, which is parameterized by a mass variance 𝑚 var computed as:</p><p>For each alloy, we obtain two predictions of 𝜅 ave as a function of composition (one for each endpoint), which we take as upper and lower bounds for predictions.  46 Although the predictions in the previous section indicate that Bi2(Se0.55S0.45)   would adopt the Pnma structure, the steep decrease in 𝜅 latt with S content for 0 < 𝑥 < 0.1 suggests that alloying with a small amount of S may be beneficial to reducing the 𝜅 latt of R3 ̅ m Bi2Se3.</p><p>Similarly, although Bi2(Se1-xSx)3 with 𝑥 = 0.1 may preferentially adopt the R3 ̅ m structure, compositions up to around 𝑥 = 0.3 are predicted to yield a similar or lower 𝜅 latt than Pnma Bi2Se3, which is itself significantly lower than the thermal conductivity of the R3 ̅ m phase.</p><p>Despite the approximations in this model, several of the predictions are nonetheless a reasonable qualitative match to the experimental measurements in Ref. 19 . The measurements on the Bi2(S,Se)3 system indicate that Bi2S2Se and Bi2SSe2, which both adopt the Pnma phase, have a lower 𝜅 latt than Bi2S3. Similarly, the measurements on Bi2(Se,Te)3 predict a minimum in the 𝜅 latt for Bi2Se0.3Te2.7, a reduced 𝜅 latt for Bi2SeTe2, and an increased 𝜅 latt for Bi2Se2Te compared to Bi2Te3. The interpolation predicts the minimum to occur at Bi2Se0.15Te0.85, which is between the two compositions measured to have 𝜅 latt lower than Bi2Te3, and the model also predicts an increase in 𝜅 latt above the Bi2Te3 endpoint at Bi2Se2Te (𝑥 = 1/3). We therefore conclude that this model can make useful qualitative, if not quantitative, predictions.</p><p>As noted above, with a 𝑇 F of 800 K we find that the Bi2SeTe2 composition is heavily skewed towards a small number of structures, which makes it feasible to compute an approximate 𝜅 latt from explicit calculations on these structures. We select the two most probable structures, with 𝑃 𝑛 = 0.42 and 0.2, and compute an averaged 𝜅 latt as a function of temperature from a partial weighted sum according to (c.f. Eq. 12): where the sum runs over a subset 𝑛 𝑥 ′ = 2 of the 𝑛 𝑥 = 21 structures in the model and we have omitted the dependence on the composition 𝑥 for brevity . Figure 11 shows the predicted temperature dependence of the 𝜅 latt of Bi2SeTe3 obtained for 𝑇 F = 800 K. We also compare this to two models using the Bi2Te3 force constants, viz. one where we change the mass of the Te atom at the interior Chal (1) site from 𝑚 Te to 𝑚 Se , and a second where we use an average mass and isotope-scattering model to approximate a random distribution of Se and Te over both chalcogen positions. The latter third model is equivalent to the analysis in Figure 10. Surprisingly, and in contrast to experiments, we find that the partial weighted sum predicts a 𝜅 latt higher than Bi2Te3. The most probable of the two structures with 𝑃 𝑛 = 0.42 corresponds to an ordered phase where the Chal (1) and Chal (2) sites are fully occupied by Se and Te atoms. The calculated in-plane and out-of-plane 𝜅 latt of this configuration at 𝑇 = 300 K is 2.55 and 1.37 W m -show stronger vdW interactions in the Chal (2) position, so this structural arrangement of Bi2SeTe2 may in fact optimize the chemical bonding strength compared to either endpoint.</p><p>The theoretical study on monolayer Bi2SeTe2 with the same structural ordering in Ref. 48 reports a 𝜅 latt of 1.96 W m -1 K -1 , which is 23 % lower than our calculated in-plane value. The similar theoretical studies in Ref. 42 report room-temperature thermal conductivities of 3.92 and 1.35 W m -1 K -1 , depending on whether or not spin-orbit coupling (SOC) is included in the calculations, and our value falls between these two extremes. (We note in passing, however, that the large effect of SOC suggested by the latter study is somewhat surprising, especially given the small effect on the phonon dispersion.)</p><p>The second most probable structure (𝑃 𝑛 = 0.2), on the other hand, has a much smaller 𝜅 ave of 0.62 W m -1 K -1 , which is lower than Bi2Te3. Since the other 19 of the 21 structures in the alloy model would also show some level of chalcogen disorder compared to the ideal arrangement, we would expect these to have a 𝜅 latt closer to the second structure, which would in turn bring the complete weighted sum closer to the ~10 % smaller 𝜅 latt relative to Bi2Te3 measured in Ref. 19 .</p><p>The second model for the 𝜅 latt of Bi2SeTe2, using the Bi2Te3 structure and force constants but with the mass of the Chal (1) ion adjusted from 𝑚 Te to 𝑚 Se , gives remarkably similar results to the partial weighted sum. This model incorporates the increased Chal (1) mass but neglects the enhanced bond strength, and thus gives a result intermediate between the two explicit calculations in the partial weighted sum. Finally, the third model, assuming a random distribution of Se and Te over the three chalcogen sites, predicts a 300 K 𝜅 ave of 0.59 W m -1 K -1 , which is a ~30 % reduction compared to Bi2Te3. This larger than the 10 % observed experimentally, and emphasizes the impact of the inhomogeneous distribution of chalcogen atoms predicted in the previous section. shows remarkable agreement with experiment, and suggests that the stability window over which the R3 ̅ m alloy forms could possibly be extended using low-temperature synthesis routes.</p><p>Finally, analysis of the 𝜅 latt of the alloys using approximate models suggests that the reduction in the thermal conductivity of R3 ̅ m Bi2(Se1-xSx)3 is likely to be driven predominantly by mass disorder, whereas strain-field fluctuations due to inhomogeneous chemical bonding play a significant role in the reduced 𝜅 latt of Pnma Bi2(Se1-xSx)3 and R3 ̅ m Bi2(Se1-xTex)3. Nonetheless, we find that using an isotope-scattering model to predict the effect of mass disorder allows for useful qualitative predictions. A more detailed analysis of the thermal transport in Bi2SeTe2 indicates that disordered occupation of the chalcogen atoms away from the preferred sites, favored by configurational entropy, is likely to be an important factor in the reduced 𝜅 latt of this compound compared to Bi2Te3 observed in experiments.</p><p>Overall, this study highlights the utility of first-principles modelling in predicting and understanding the stability, structural dynamics and heat transport in thermoelectric materials, for predicting the structure and energetics of moderately-complex alloys, and for establishing the impact of alloying on the thermal conductivity. Given that alloying is an accepted strategy for optimizing the performance of thermoelectric materials, we hope that this work and similar studies will provide useful guidance to ongoing work on the bismuth chalcogenides and related systems.</p>