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. The IFCs were determined using the supercells shown in Table 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 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. Thermal conductivity calculations were performed within the single-mode relaxation-time approximation (RTA). The modal properties were sampled on the regular Γ-centered 𝐪-point meshes in Table . 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.