Anharmonicity in the high-temperature Cmcm phase of SnSe: soft modes and three-phonon interactions

The layered semiconductor SnSe is one of the highest-performing thermoelectric materials known. We demonstrate, through a first-principles lattice-dynamics study, that the high-temperature Cmcm phase is a dynamic average over lower-symmetry minima separated by very small energetic barriers. Compared to the low-temperature Pnma phase, the Cmcm phase displays a phonon softening and enhanced three-phonon scattering, leading to an anharmonic damping of the low-frequency modes and hence the thermal transport. We develop a renormalisation scheme to quantify the effect of the soft modes on the calculated properties, and confirm that the anharmonicity is an inherent feature of the Cmcm phase. These results suggest a design concept for thermal insulators and thermoelectric materials, based on displacive instabilities, and highlight the power of lattice-dynamics calculations for materials characterization.

( + ) ⁄ , where and are the Seebeck coefficient and electrical conductivity, and and are the lattice and electronic thermal conductivities, respectively. The key to developing high-performance thermoelectric materials is to reduce the thermal conductivity while maintaining a high thermopower 2 . Widespread application requires a above 2 at the target operating temperature. [3] Historically, the lead chalcogenides set the benchmark for high thermoelectric performance due to their unique combination of favourable electrical properties and stronglyanharmonic lattice dynamics. [4][5][6][7][8][9][10][11] However, bulk SnSe was recently shown to be a very promising high-temperature thermoelectric, with a score surpassing the record set by nanostructured PbTe. [12] This has been ascribed to an ultralow lattice thermal conductivity, arising from its pseudolayered structure and strong high-temperature anharmonicity. [13,14] Compared to PbTe, SnSe achieves superior thermoelectric performance without the need for doping or other material modifications, [15] which can be detrimental to electrical properties. It is therefore important to elucidate the microscopic origin of its high performance, in order that the same ideas may be incorporated into design strategies for improving thermoelectrics.
SnSe displays a second-order phase transition from Pnma to Cmcm phases between 700-800 K. [13,14] Modelling studies on the thermal transport of the low-temperature phase have been carried out, [16] while combined theoretical and computational studies of both phases [14,17] have connected the low lattice thermal conductivity to the anharmonicity arising from the displacive phase transition. However, there has not yet been a detailed comparison of the intrinsic lattice dynamics of the two phases, in part due to the difficulties associated with modelling phonon instabilities.
Previous theoretical studies have shown that firstprinciples lattice-dynamics calculations can provide deep insight into the material physics underlying thermoelectric performance. [18][19][20][21][22][23][24][25] In this work, we have employed firstprinciples lattice-dynamics calculations to model both phases of SnSe, with a particular focus on characterising the lattice dynamics and thermal transport in the high-temperature phase, and on developing a practical approach to the theoretical challenge of modelling soft-mode instabilities.
Our lattice-dynamics calculations were performed using the Phonopy [26,27] and Phono3py [28] packages, with the VASP code [29] as the force calculator. We used the PBEsol functional, [30] which we have previously found to give a very good description of the lattice dynamics in a range of semiconductors. [31] Projector augmented-wave (PAW) pseudopotentials [32,33] treating the Sn 4d, 5s and 5p and the Se 4s and 4p electrons as valence were used to model the ion cores. The electronic wavefunctions were expanded in a plane-wave basis with a kinetic-energy cutoff of 500 eV. An 8×4×8 Monkhorst-Pack k-point mesh was used to sample the Brillouin zone, correspondingly reduced for the supercell force calculations. [34] The second-order (harmonic) force constants were computed using a carefully-chosen 6×1×6 supercell expansion for both phases and a displacement step of 10 -2 Å. The third-order force constants were computed using 3×1×3 expansions and a 3×10 -2 Å step. Full technical details of our simulations, along with supercell-size convergence tests for the harmonic-phonon calculations, can be found in the supporting information. [41]). The 2D potential-energy surface along the two corresponding normal-mode coordinates ( 1 / 2 ) is shown in plot (e). The left-and right-hand snapshots in (c)/(d) are taken through the ab and ac planes, respectively, and the Sn and Se atoms in the models are coloured silver and green. The images were prepared using VESTA (Ref. [35]).
The optimised equilibrium structures of the two phases are shown in Figs. 1c and 1d. The calculated lattice parameters of the Pnma and Cmcm phases are = 4.367, = 11.433 and = 4.150 Å, and = 4.217 Å, = 11.525 and = 4.204 Å, respectively. These values are consistent with experimental characterization (Pnma: a = 4.44 Å, b = 11.49 Å, c = 4.135 Å; Cmcm: a, c = 4.31 Å, b = 11.70 Å), and with other theoretical studies. [13,16] The low-temperature phase is a distortion of the more symmetric high-temperature structure, leading to a primitive unit cell with twice the volume, and both structures can be thought of as distortions to the rocksalt structure of the Pb analogue PbSe.
Comparing the phonon dispersion and density-of-states (DoS) curves of the two phases (Figs. 1a/1b), the upperfrequency part of the Cmcm DoS displays a noticeable red shift with respect to the Pnma phase. From a numerical integration,[41] around 15 % of the states in the lowtemperature phase lie between 3-4 THz, compared to ~30 % in the high-temperature phase. Moreover, whereas the dispersion of the Pnma structure shows real frequencies across the Brillouin zone (Fig. 1a), the dispersion of the hightemperature phase displays prominent imaginary modes along the line between Γ and Y (Fig. 1b), with the primary soft modes at the two symmetry points. [14,17] The imaginary modes make an almost negligible contribution to the overall phonon DoS, suggesting them to be confined to a small reciprocal-space volume. This was confirmed by inspecting the frequencies at the commensurate -points in the various supercells used for convergence testing, as well as in a set of 2× ×2 expansions with a high density of commensurate points along the Γ-Y line. [41] From an analysis of the phonon eigenvectors,[41] both modes correspond to symmetry-breaking displacive instabilities. The optic mode at Γ corresponds to motion of the four (two) Sn and Se atoms in the conventional (primitive) cell in opposing directions parallel to the c axis. The acoustic mode at Y is a similar motion along the a axis, but with the atoms in the two bonded layers moving in opposite directions. From mapping out the 2D potential-energy surface spanned by the two mode eigenvectors (Fig. 1e), it can be seen that the Y mode represents the primary distortion to the lowtemperature phase: the minimum along the Γ-point mode (λ 1 in Fig. 1b/ 1 in Fig. 1e) is another saddle point on this potential-energy surface, whereas the global energy minimum lies along the Y mode (λ 2 / 2 ). It can also be seen from a comparison of the Pnma and Cmcm structures in Fig.  1c/d that the symmetry breaking induced by this mode maps between the high-temperature to the low-temperature phase.  Fig. 1b as a function of the corresponding normal-mode coordinates =1,2 . The markers show the calculated points, the solid lines are fits to a 20-power polynomial, and the black lines inside the potentials show the eigenvalues obtained by solving a 1D Schrödinger equation. Plots (b) and (d) show the effective renormalized frequencies of the two modes (̃= 1,2 ) as a function of temperature, calculated as the harmonic frequencies which reproduce the contributions of the modes to the vibrational partition function. Plot (e) compares phonon dispersions and densities of states calculated without renormalization (black) and with the imaginary modes renormalized to the calculated 0 (blue), 300 (red) and 800 K (orange) frequencies.
Each imaginary mode forms an anharmonic double-well potential (Figs. 2a/2c). The potential along the Γ-point mode is steeper, accounting for its more negative harmonic frequency, but, as seen in Fig. 1e, the potential along the Y mode possesses a deeper minimum. The respective minima are both shallow, at 5 and 11 meV below the average structure. These findings indicate that the Cmcm phase is effectively an average structure, in a manner analogous to the high-temperature phases of some oxides [36] and oxide/halide perovskites (e.g. CsSnI3). [25,37,38] To confirm that the Cmcm phase is unlikely to become stable under lattice expansion (or contraction) at finite temperature, we computed phonon dispersions and DoS curves over a range of expansions and contractions about the athermal equilibrium volume,[41] which indicated that the imaginary modes soften further under lattice expansion, and also persist under moderate applied pressure.
Soft modes present a challenge to lattice-dynamics calculations, as techniques for treating this anharmonicity (e.g. self-consistent phonon theory [39,40]) are often impractically expensive. Whereas in real systems phonon instabilities are localized to regions of infinitesimal volume in reciprocal space (i.e. points, lines or planes), and thus do not contribute to thermodynamic properties such as the free energy, [36] the finite Brillouin-zone sampling and interpolation used in lattice-dynamics calculations "spreads" the imaginary modes over a finite volume.
In order to investigate the effect of the imaginary modes in the high-temperature phase on our calculations, we implemented a simple renormalization scheme to calculate approximate effective harmonic frequencies. The potential energy along each mode as a function of the respective normal-mode coordinate, =1,2 , is fitted to a 20-power polynomial and 1D Schrödinger equations for the potentials are solved to obtain the eigenvalues inside the anharmonic double wells (Figs. 2a/2c). Effective (real) harmonic frequencies for the imaginary modes are then calculated as a function of temperature to reproduce the contribution of the modes to the anharmonic vibrational partition function (Figs. 2b/2d). These effective frequencies are used to adjust the harmonic force constants, by back transforming the dynamical matrices at the commensurate -points in the supercell expansion after adjusting the frequencies of the imaginary modes. The corrected force-constant matrices are then used as input for further post processing. A detailed overview of the renormalization scheme and its implementation is given in the supporting information. [41] The Γ and Y points in the Cmcm primitive cell are commensurate with our chosen 6×1×6 supercell expansion, but the expansion has no commensurate points along the Γ-Y line. The dispersion appears to be well reproduced with this expansion (c.f. Fig. 1b), however, and it was thus employed for renormalization at the symmetry points. The renormalized 0 K frequencies of the imaginary modes at Γ and Y are 0.54 and 0.45 THz, respectively. The frequencies initially decrease with temperature, as the system explores higher energy levels within the minima, and reaches a minimum close to where the available thermal energy is comparable to the barrier height (approx. 40 and 80 K, respectively, for the two modes). Above these temperatures, the renormalized frequencies increase monotonically as the system accesses the steeper parts of the potential. The renormalized frequencies at 300 K are 0.35 (Γ) and 0.23 THz (Y), and increase to 0.45 (Γ) and 0.30 THz (Y) at 800 K (the temperature above which the Cmcm phase is observed crystallographically [13,14]). Fig. 2e shows the impact of renormalizing the imaginary modes at 0, 300 and 1000 K on the phonon dispersion and DoS. The effect on the orthogonal phonon branches is minimal, and the correction is strongly localized to the lines along which the imaginary modes occur. The latter is primarily a result of the large supercell used for calculating the force constants, which we consider a prerequisite for this renormalization scheme. The effect of the renormalization on the phonon DoS is minimal, with the most notable differences occurring up to ~1 THz (Fig. 2e, inset). (as defined in Ref. [28]), overlaid on the phonon densities of states. The error bars show the standard deviation on the average in each histogram bin. Plots (d) -(e) show the calculated lattice thermal conductivity ( ) along each of the three crystallographic axes, together with the isotropic average. The Cmcm-phase data shown in (b) and (e) was modelled without renormalization of the imaginary modes, while that shown in (c) and (f) was modelled with the second-order (harmonic) force constants corrected using the 800 K renormalized frequencies.
The renormalization was also found to have a negligible impact on thermodynamic functions calculated within the harmonic approximation,[41] with a maximum difference in the high-temperature (1000 K) vibrational Helmholtz free energy of < 0.3 kJ mol -1 per SnSe formula unit (~0.3 %; c.f. = 8.314 kJ mol -1 at this temperature).
To investigate the effect of the phase transition on the thermal transport, we modelled the phonon lifetimes within the single-mode relaxation-time approximation, which were used to solve the Boltzmann transport equations. The lifetimes are calculated as the phonon self-energy within many-body perturbation theory, using the three-phonon interaction strengths ϕ ′ ′′ , obtained from third-order force constants ϕ along with an expression for the conservation of energy. This method is implemented in the Phono3py code, and a detailed overview may be found in Ref. [28]. In these calculations, we did not attempt to correct the third-order force constants for the imaginary modes in the Cmcm phase, but instead performed post-processing for the hightemperature structure using renormalized second-order force constants.
Neglecting quasi-harmonic effects from volume expansion, ϕ ′ ′′ is temperature independent, and a histogram showing the average phonon-phonon interaction strengths, , across the phonon DoS ( Fig. 3a-3c) provides a convenient means to compare the inherent "anharmonicity" of the two phases. This comparison clearly shows that the low-frequency modes in the Cmcm phase experience a substantially stronger interaction with other phonon modes compared to the lower-frequency branches in the Pnma phase. Correcting the second-order force constants reduces the interaction strength up to ~0.5 THz, but the average up to ~3 THz remains consistently higher than in the Pnma phase. These calculations therefore indicate that the lattice dynamics of the high-temperature phase are inherently more anharmonic.
The calculated 300 K lattice thermal conductivity ( ) of the Pnma phase ( Fig. 3d) is 1.43, 0.52 and 1.88 W m -1 K -1 along the , and axes, respectively. These values are in reasonable agreement with measurements of ~0.7 and ~0.45 W m -1 K -1 along the two short and one long axes, [13] and the predicted axial anisotropy is consistent with other studies. [13,16] For the Cmcm phase ( Fig. 3e), we calculated an isotropic average of 0.33 W m -1 K -1 at 800 K compared to the measured value of ~0.25 W m -1 K -1 . [13] Given the neglect of volume expansion, higher-order anharmonicity and other potential scattering processes in the present calculations, this is again reasonable agreement. As in the Pnma phase, the calculations predict an axial anisotropy in , with significantly reduced transport along the long axis. With the second-order force constants corrected, we calculated an average of 0.34 (an increase of 2.4 %), indicating that renormalization has only a small quantitative impact. Analysis of the modal contributions to the thermal transport at 300 and 800 K[41] shows that the low-frequency modes account for the bulk of the heat transport in both phases, and that the larger phonon-phonon interaction strengths in the Cmcm phase significantly reduce the lifetimes of these modes, thus explaining its reduced .
These findings are consistent with the mechanism proposed in Refs. [14,17]. Our calculations show that the high-temperature phase exhibits a substantial phonon softening compared to the Pnma phase, and that enhanced phonon-phonon interactions in the Cmcm phase lead to an anharmonic damping of the low-frequency modes and hence the thermal transport. Moreover, our renormalisation scheme provides a means to explore the physics of the soft modes and their effect on calculated properties, such as the thermodynamic free energy and the thermal transport, from first principles, without requiring experimental measurements of the phonon frequencies to fit to. [14] In this regard, our scheme could easily be combined with the quasi-harmonic approximation, which would allow us to explore the effect of thermal expansion at finite temperature. This would be expected to further dampen the thermal conductivity, [24] and detailed lattice-dynamics calculations could allow the relative contributions from changes in phonon lifetimes and group velocities to be discerned. [14] In summary, our calculations demonstrate that the low high-temperature lattice thermal conductivity of SnSe is due to anharmonic damping of the low-frequency phonon modes. We have developed a simple renormalization scheme to quantify the impact of the phonon instabilities in the hightemperature phase on properties calculated using secondorder force constants, and hence shown the enhanced anharmonicity to be an inherent property of this system. This scheme may form a practical basis for studying other important classes of system with displacive instabilities, e.g. halide perovskites. From a materials-design perspective, similar anharmonic phonon dampening may occur in other systems at the boundary of a phase transition, and so this could serve as a selection criterion for identifying materials with ultra-low thermal conductivity. [14] In these materials, the poor thermal transport is a bulk property, and so the potential negative impact on electrical properties of modifications such as doping and nanostructuring may be avoided. Understanding this phenomenon may thus provide a robust design strategy for developing thermal insulators and high-performance thermoelectric materials.
We gratefully acknowledge helpful discussions with Dr J. M. Frost. JMS is supported by an EPSRC programme grant (grant no. EP/K004956/1). LAB is currently supported by the JSPS (grant no. 26.04792). JB acknowledges support from the EPSRC (grant no. EP/K016288/1). Calculations were carried out on the ARCHER supercomputer, accessed through membership of the UK HPC Materials Chemistry Consortium, which is funded by the EPSRC (grant no. EP/L000202). We also made use of the Bath University HPC cluster, which is maintained by the Bath University Computing Services, and the SiSu supercomputer at the IT Center for Science (CSC), Finland, via the Partnership for Advanced Computing in Europe (PRACE) project no. 13DECI0317/IsoSwitch.

APPENDIX: Data-access statement
In addition to the data in the electronic supporting information,[41] the optimised structures and raw/processed data from the phonon calculations are available online, free of charge, from https://github.com/WMD-group/Phonons. This repository also includes the codes used to implement our approximate renormalization scheme. *a.walsh@bath.ac.uk

Ab initio lattice-dynamics calculation protocol
First-principles calculations were carried out within the pseudopotential plane-wave density-functional theory formalism, as implemented in the Vienna ab initio Simulation Package (VASP) code. [1] We used projector augmentedwave (PAW) pseudopotentials, [2,3] treating the Sn 5s, 5p and 4d and the Se 4s and 4p electrons as valence states, in conjunction with the PBEsol exchange-correlation functional [4] and a plane-wave basis with a kinetic-energy cutoff of 500 eV. An 8×4×8 Monkhorst-pack -point mesh [5] was used to sample the first Brillouin zone of the conventional Pnma and Cmcm cells, which was correspondingly reduced for the supercell-phonon calculations. The PAW projection was performed in reciprocal space, and non-spherical contributions to the gradient corrections inside the PAW spheres were taken into account. The electronic-structure calculations were performed within the scalar relativistic approximation, excluding spin-orbit coupling.
During geometry optimizations, a tolerance of 10 -8 eV was applied during the electronic minimisation, and the ion positions and lattice parameters were optimized until the magnitude of the forces on the ions was below 10 -2 eV Å -1 .
Lattice-dynamics calculations were performed with the Phonopy [6,7] and Phono3py [8] packages, which were used to obtain sets of second-and third-order force-constant matrices, respectively, via the supercell finite-displacement method. [9] The second-order force constants were calculated using 6×1×6 expansions of the conventional Pnma and Cmcm cells, containing 288 atoms, with a displacement step size of 10 -2 Å. We found that these cells were sufficiently large to converge the shape of the phonon density of states (DoS; see Section 2, below). We also performed additional calculations on 2×6×2 and 6×6×2 supercells to obtain higher-quality dispersions along the Γ-Y and S-Γ segments of the Cmcm band dispersion, respectively (see Fig. 1b in the text). Due to the unfavourable scaling of the number of inequivalent two-atom displacements with supercell size, the third-order force constants were calculated from 3×1×3 supercell expansions containing 72 atoms using a step size of 3×10 -2 Å. During the post processing, the phonon DoS curves were constructed by evaluating the phonon frequencies on a uniform 48×48×48 Γ-centred -point grid, while the phonon lifetimes used to model the thermal transport were sampled on 16×16×16 Γ-centred grid.
The phonon DoS and thermal conductivity of the Cmcm phase were calculated with respect to the conventional cell, while the dispersion was calculated in the primitive basis. The transformation from the conventional to primitive cell was performed using the transformation matrix: The dispersions of the Pnma and Cmcm phases (the latter with respect to the primitive cell) were constructed by evaluating the phonon frequencies along a path connecting the high-symmetry points in the respective Brillouin zones.
The reduced -vectors of the symmetry points used to generate the Pnma and Cmcm dispersions are listed in Tables S1 and S2, respectively.

Supercell-size convergence for the finite-displacement calculations
As noted in the text and in Section 1 above, we carefully converged the supercell sizes used in the finitedisplacement calculations on both the Pnma and Cmcm phases.     A more quantitative way to assess the supercell-size convergence is to compare the vibrational constant-volume (Helmholtz) free energies, , which are computed by summation over a grid of -points sampling the phonon Brillouin zone according to (eq. S1): is the Boltzmann constant, is the temperature, and the product runs over the phonon modes with frequencies . Cmcm phases, respectively, which corresponds to < 1 %. This is also an order of magnitude smaller than at 1000 K (8.314 kJ mol -1 ), and is at least comparable to the error one would expect from using an approximate exchangecorrelation functional, suggesting that for this system the free energies are relatively insensitive to the set of forceconstant matrices.
One potential issue with ×1× expansions is that while the Γ and Y points in the Cmcm primitive cell, where the principal soft modes occur, are both commensurate, these supercells contain no commensurate points along the line between them. To check the convergence of the Γ-Y segment of the dispersion, we therefore performed a systematic set of calculations on 2× ×2 expansions, with = 1, 2, 4, 6 and 8 (Fig. S3).
As noted in the text, these tests confirmed that the soft modes lie along the line between Γ and Y, rather than being localised to the symmetry points. However, we found that the smaller expansions, including = 1, reproduced the shape of the dispersion of the soft modes along this segment obtained from the larger calculations. The main effect of larger values of is seen in the dispersion of the highest-energy optic branch, introducing oscillations which may be due to our not including non-analytical corrections for LO/TO splitting in these calculations. However, the effect on the phonon DoS is subtle, confirming that including longer-range interactions along the a and b directions are more important for reproducing the distribution of phonon frequencies. Based on these considerations, we selected 6×1×6 expansions for our "production" calculations, as a balance between computational cost and accuracy.  directions, and is not seen in supercells with larger expansions along these directions.

Imaginary-mode mapping and renormalization scheme
The potential-energy surfaces along the two imaginary modes in the equilibrium Cmcm structure (Figs. 1e and 2a/c in the text) were mapped by performing single-point total-energy calculations on a series of structures of the conventional cell with atomic displacements along the mode eigenvectors over a range of amplitudes "frozen in". Given the mode eigenvectors of a set of phonon modes , together with the corresponding normal-mode coordinates (amplitudes) , the displacement , of the th atom in the th unit cell can be calculated according to: where , is the component of on atom , are the atomic masses, is the number of atoms in the supercell used to model the displacement, is the phonon wavevector and are the positions of the atoms. We note that absorbs the time dependence of the position.
As noted in the text, the mode at in the Brillouin zone of the Cmcm primitive cell folds to Γ in the reciprocal space of the conventional cell. In the special case of mapping Γ-point modes, Eq. S2 can be simplified to: In the present study, the potential-energy surfaces were mapped for values of 1 = ±30 and 2 = ±45 amu Each polynomial was evaluated on a 1D grid of points, and the resulting potential substituted into the 1D SE and solved for the eigenvalues and eigenvectors by means of a Fourier transform followed by matrix diagonalization using the EISPACK routines. [11] This procedure and its advantages compared to solutions in real space are explained in detail in Ref. [12], and we have verified its efficiency against more conventional "shooting method" approaches.
The eigenvalues, , are used to determine the partition function, , via: where is the Boltzmann constant and is the temperature. The number of terms included is chosen so that the addition of the final term changes by less than 10 -6 . By setting this expression for equal to the harmonic partition function: where ̃ is the effective harmonic frequency, we can derive that: which provides an effective renormalized harmonic frequency for the mode at a target temperature that reproduces its contribution to the thermodynamic partition function (Figs. 2b/2d in the text).
The extent of the potential in and the grid density were both carefully converged, with a basis of 512 grid points and 1 = ±25 and 2 = ±35 amu 1 2 Å, in both cases extending to ~1 eV in energy above the average structure, found to be sufficient to obtain renormalized frequencies converged to < 10 -2 THz.
These effective frequencies were used to adjust the force constants using the Python API exposed by the Phonopy code. [6,7] The original force/displacement sets and corresponding force constants were used to calculate the phonon frequencies and eigenvectors in the conventional cell at the -points commensurate with our chosen supercell expansions. The frequencies of the imaginary modes at the Γ point were then set to the calculated effective (real) ones, and the corrected dynamical matrices back-transformed to a new set of force constants to be used in subsequent postprocessing steps.
A set of scripts for setting up and post-processing the displacement-mapping calculations, source code and example input files for the 1D SE solver, and a script for patching the force constants using the Phonopy API, are available as additional data (see the Appendix in the main text).

Effect of renormalisation on thermodynamic functions
To quantify the effect of renormalizing the imaginary modes in the equilibrium Cmcm structure on the thermodynamic functions derived within the harmonic approximation, we compared the temperature dependence of four The results are shown in Fig. S8, and the calculated vibrational zero-point energies ( ) and 1000 K free energies ( ,1000K ) are collected in Table S5.   Fig. 1b in the text (̃= 1,2 ) that were used are listed. The and ,1000K values for the temperature-dependent renormalization are the same as those for the 0 and 1000 K constant renormalized frequencies, respectively.
The data in Table S5 suggests that the renormalization has a minimal impact on the free energy, leading to differences in the zero-point energy on the order of 10 -3 kJ mol -1 per SnSe formula unit, and to differences in the hightemperature free energy of < 0.3 kJ mol -1 (~0.3 %). This is an order of magnitude smaller than at 1000 K (8.314 kJ mol -1 ). From the comparison in Fig. S8, this difference occurs predominantly as a result of an increase in , although there is also a small increase in , which is roughly the same for all four renormalized curves.

Detailed analysis of the thermal transport in the Pnma and Cmcm phases
To investigate the origin of the reduction in lattice thermal conductivity ( ) between the low-temperature Pnma and high-temperature Cmcm phases, we performed a detailed analysis of the modal contributions to the roomtemperature (300 K) and 800 K lattice thermal conductivities of the two phases.
A comparison of the cumulative lattice thermal conductivity as a function of phonon frequency for the Pnma phase at 300 K against the phonon DoS (Fig. S9a) indicates that ~70 % of the heat transport is due to the lower-frequency modes (with frequencies < ~3 THz). This rises to around 80 % in the Cmcm phase (Fig. S9b), which can be partly explained by the red shift of some of the high-frequency modes visible in Fig. S4 (see Section 3, above). Renormalizing the soft modes has little effect on the shape of the DoS nor the cumulative thermal conductivity distribution (Fig. S9c). A similar comparison at 800 K (Fig. S13) yields a practically identical analysis.
Within the single-mode relaxation-time approximation, the macroscopic thermal-conductivity tensor is obtained as a summation of modal contributions according to: where is the number of unit cells in the crystal (equivalently, the number of reciprocal-space grid points used to sample the lifetimes), 0 is the unit-cell volume, are the modal heat capacities, , are the mode group velocities, and are the mode lifetimes. The group velocities and modal heat capacities are calculable within the harmonic approximation according to: The procedure for calculating lifetimes followed in this work is explained in detail in Ref. [8]. The lifetimes are derived from three-phonon interaction strengths, ϕ ′ ′′ , together with and an expression for the conservation of energy. ϕ ′ ′′ are calculated according to:  where are the phonon occupation numbers, and the delta functions enforce conservation of energy. Assuming threephonon processes to be the dominant scattering effect, Γ ( ) are related to the phonon lifetimes according to: where 2Γ ( ) are the full-width at half-maxima of the phonon lines.
To analyse the thermal transport in more detail, we compared the spread of the (isotropically-averaged) modal contributions to the thermal conductivity, , as a function of frequency against the spreads of the group velocity norms, | , |, lifetimes, , and average three-phonon interaction strengths, , . Following Ref. [8], , are defined here as: where is the number of atoms in the unit cell, and 3 is thus the number of phonon bands at each -point.
Figs. S10 and S11 give a breakdown of the modal contributions to the room-temperature thermal conductivity of the Pnma and Cmcm phases. In the Pnma phase, the dominant contribution is clearly from long-lived low-frequency modes, while higher-frequency phonons with larger group velocities account for a secondary contribution. The lifetimes of the low-frequency modes are significantly reduced in the high-temperature phase, and mid-frequency phonons with large group velocities make a proportionately higher contribution to the bulk thermal transport. A corresponding analysis of the 800 K thermal conductivity (Figs. S14, S15) shows a very similar pattern, but with a marked reduction in the phonon lifetimes of both phases, particularly of the low-frequency modes.
Finally, we also analysed the contributions of modes with different mean-free paths (MFPs) to the roomtemperature thermal conductivity ( Fig. S12; the MFP is calculated as MFP = | , | ). We found that phonons with MFPs below ~1nm make a very small contribution to the overall heat transport, while the longest phonon MFP observed in either phase was ~1 μm. The distribution of cumulative lattice thermal conductivity as a function of increasing phonon MFP is skewed more towards longer paths in the low-temperature Pnma phase than in the high-temperature Cmcm system, although in both around 50 % of the thermal conductivity is through phonons with MFPs < 10-20 nm. Again, renormalization of the soft modes in the high-temperature phase makes little difference to the spread of the /MFPs, and has a negligible effect on the shape of the thermal-conductivity distribution. At 800 K (Fig. S16), the distributions maintain a similar overall shape, but are skewed towards phonons with shorter MFPs, as would be expected given the reduction in lifetimes evident from comparing Figs. S10/11 and S14/15.
As discussed in the text, these analyses indicate that the comparative reduction in the lifetimes of the lowfrequency phonon modes in the Cmcm phase is primarily responsible for its lower lattice thermal conductivity. This can be ascribed to the higher average three-phonon interaction strength experienced by modes in the lower-frequency part of the DoS in the high-temperature phase (Fig. 3 in the text), which belies an inherently higher anharmonicity.