# Extended anharmonic collapse of phonon dispersions in SnS and SnSe

Sep 4, 2020

### Diffraction/phase transition

The Pnma–Cmcm structural distortion and our single-crystal neutron diffraction results on SnS are shown in Fig. 1. Pnma notation is used throughout to describe the (hkl) peak indices, corresponding to conventional orthorhombic cells in Fig. 1a, e. Owing to the doubling of the primitive cell volume in Pnma, superstructure Bragg peaks occur at zone-boundary Y points of the primitive Cmcm cell, for example at (201) and (102) in the low-temperature phase, and disappear above TC (~880 K) (Fig. 1b, f). Panels 1c, d, g show that the intensities of superstructure Bragg peaks (102) and (201) continuously decrease on warming until full suppression at TC and above. Conversely, (002) and (101) are not superlattice peaks and their intensities increase with temperature. The continuous evolution of the peak intensities and lattice parameters (Fig. 1h) are supportive of a second-order (or weakly first-order) transition scenario. A strong anisotropy is also observed, with the c-axis contracting, while a and b exhibit a normal expansion on heating, see Fig. 1h and prior studies11,13,24. The thermal contraction along c in the Pnma phase was partially explained by a quasiharmonic (QH) analysis of the free-energy, yet the stronger lattice parameter variations close to TC pointed to strong intrinsic anharmonicity in this regime25. Above TC, the change in lattice parameters becomes weaker. From our complementary heat capacity measurements (Supplementary Fig. 18), the phase transition in SnS occurs at 875 K which is consistent with our neutron data and prior results24.

### Phonon dispersions

A drastic change in phonon dispersions is observed across the structural phase transition in both SnS (Fig. 2) and SnSe (Supplementary Figs. 12). High energy resolution was required to accurately map low energy modes in the dynamical susceptibility  χ”(Q,E), and, therefore, we used the cold neutron chopper spectrometer (CNCS) at the spallation neutron source (SNS) and the cold triple-axis spectrometer (CTAX) at the high flux isotope reactor (HFIR) (details in “Methods”). Additional measurements with thermal neutrons were also performed with a thermal triple-axis spectrometer (HB-3) at HFIR and the wide angular-range chopper spectrometer (ARCS) at SNS. The temperature evolution of INS phonon measurements in SnS is shown and compared with simulations in Fig. 2. Panels a–d show the Pnma phase and e–h the Cmcm phase. Calculated phonon dispersions for the two phases are shown in panels a and e with polarization effects being clearly evident in b and f where the simulated χ”(Q,E) gives intensity only for the c-polarized TAc and TOc phonons. Experimental χ”(Q,E) maps for these dispersions are seen in the remaining panels showing the TAc and TOc modes, for phonon wavevectors along [100]. As seen in panels c and d, a drastic and extended softening (decreasing energy) of both branches occurs on warming in the Pnma phase. Above TC, the re-emergence of a single TAc branch is observed in the Cmcm phase, as seen in panels g, h. Our measurements do directly reveal the condensation of a soft phonon mode at the transition (at H = 2 in panel h for Cmcm and H = 1 in panel d for Pnma which are both superstructure peaks), but also show how the instability leads to an extreme softening affecting extended portions of TA branches in both phases, as well as the folded-back lowest-energy TO branch in the distorted Pnma phase. In the Cmcm phase, it is a zone boundary TAc mode that condenses at Q’s corresponding to the Pnma super-structure Bragg peaks (e.g., (102) and (201)). The soft-mode re-emerges in Pnma as the back-folded zone-center (Γ) TOc mode. This striking evolution of the phonons was observed in multiple Brillouin zones and on measurements from different instruments, as shown in Supplementary Fig. 6. An almost identical behavior is observed in SnSe as seen in Supplementary Figs. 12 but with energies shifted slightly lower relative to SnS due to the mass of selenium being higher than sulfur.

The extremely soft nature of the TAc branch in the SnS Cmcm phase near TC can be appreciated by considering data at 930 K (Fig. 2g). The branch rises from 0  meV at (1,0,1) and peaks at only 1.36 meV at (1.5,0,1) before decreasing to 0.90 meV at (2,0,1). This softening is extended across the entire dispersion as opposed to being localized at zone centers or boundaries. As T is increased above TC, the TAc branch stiffens back up along its entirety (Supplementary Fig. 6q–t). The extended nature of this ultra-low energy TA branch is captured well by the temperature dependent effective potential (TDEP) method26 (see “Methods”) for the Cmcm calculation (panels e, f). The TDEP method directly incorporates anharmonic renormalization by extracting effective force constants through fitting the potential energy surface of temperature-dependent ab initio molecular dynamics (AIMD). The renormalized calculations show the continuous collapse of the TOc and TAc modes with temperature in the Pnma phase in Supplementary Fig. 4. The TOc softening occurs over a large temperature range, starting far below TC. At 150 K the TOc mode at Γ has an energy of 5.3 meV and it is softened by 25.5% (to 3.9 meV) at 740 K, still 140 K below the transition.

We note that Fig. 2 panels c and d show systematic temperature variations in the intensities of phonon branches. This evolution of χ(Q,E) reflects thermal changes in both Bragg peak intensities and anharmonic phonon eigenvectors (Supplementary Figs. 45). At 150 K the TAc branch (<4 meV) shows the greatest intensity closest to (102), whereas at 860 K intensity is highest near (002). In Supplementary Fig. 6 we also observe drastic changes in the relative intensities of the TOc mode at (103) and (203) (panels a–e) and (002) and (102) (panels f–j). The thermal evolution of the TOc mode eigenvector obtained from our anharmonic phonon simulations is shown in Supplementary Fig. 27.

### Soft mode behavior

We now proceed to quantify the energy, E0, and linewidth (ΓLW, full width at half maximum) of the soft mode. The phonon linewidths are inversely proportional to phonon lifetimes (τ) while frequencies impact κlat by altering group velocities (vg = kω), Bose–Einstein occupation factor (n(ET)) and the phase space for scattering. Extracting linewidths and energies from INS requires careful consideration of the instrument-specific resolution function, R(QE), which we calculated for each instrument, Q-point and energy, and then convolved with fitting functions (details in “Methods”). In addition to INS data, we also collected high-resolution Raman spectra on single crystals of SnS and SnSe to probe zone-center optical modes (see Supplementary Note 4 for details), including the soft-mode in the Pnma phase. Phonon lineshapes were successfully modeled as damped harmonic oscillators (DHO), see Supplementary equation 5, which approximates to a Lorentzian lineshape in the usual case where the phonon energy E0 is much larger than the damping ΓLW. A strong departure from Lorentzian behavior is observed in the χ(E) spectra (Fig. 3a, b) close to TC, where the phonons are overdamped with ΓLW > E0. In this regime, the traditional phonon picture breaks down as modes no longer have well-defined oscillations.

The soft mode behavior is clearly observed in the Raman spectra in Fig. 3a where the TOc (Ag) mode softens and broadens upon approaching TC from below. The soft-mode disappears above TC in the Raman spectra, since it re-emerges as a zone-boundary Y mode in the Cmcm phase. The lowest energy zone center optical mode in the Cmcm phase (seen in Fig. 2e, f at  ≈ 4.5 meV) is Raman inactive from symmetry considerations27 and as such we did not observe any peaks above TC. Fits to INS data are shown in Fig. 3b, revealing the pronounced overdamping around TC. In Supplementary Fig. 21, the DHO profiles extracted from fitting the INS data show the intrinsic lineshapes, corrected for instrumental resolution. Figure. 3c shows the temperature evolution of E(TOc), which is clearly continuous up to the condensation at TC, and can be fitted with a power law, consistent with a second-order transition. Remarkably good agreement is found between INS measurements on different instruments and the Raman data. The reason why E0 does not exactly reach zero is because of strong broadening in the DHO profile at TC, a common feature of soft-mode transitions28. The soft mode behavior can be seen on both sides of TC in the INS data (CTAX), since no modes are symmetry-forbidden in INS, unlike Raman. In both phases the mode (TOc for Pnma and zone-boundary TAc for Cmcm) softens upon approaching TC. This behavior is characteristic of soft mode transitions in general, as detailed in seminal studies29.

The TOc linewidth shows a drastic increase at the transition (Fig. 3d), reaching a value 9.4 ×  higher than at 295 K, corresponding to a drastically increased scattering rate. Importantly, this increase is far larger than the linear temperature dependence expected within low-order perturbation theory (n(ET) T at high T). We observe a similar behavior in SnSe INS and Raman spectra (Supplementary Figs. 13). This divergence behavior near a second order (continuous) phase transition can be ascribed to critical behavior, similar to the divergence in the heat capacity (Supplementary Fig. 18). Further, a strong polarization dependence of the softening is revealed from the fitting of Raman data (Supplementary Fig. 20). Thus, the anharmonicity depends anisotropically on atomic motions through the bonding directionality. ({A}_{{rm{g}}}^{1}) is a-polarized (out-of-plane) and exhibits significantly less softening than the ({A}_{{rm{g}}}^{0}) and B3g modes polarized in the bc plane. This agrees well with the findings of Liu et al.22 except here we focus on the lowest energy modes in the SnSe/SnS compounds whereas Liu et al.22 only reported modes above 50 cm−1.

We find a clear departure from mean field theory (MFT) over a large temperature range in both SnS and SnSe. MFT predicts E(T) = AT − TCα, where A is a constant and (alpha =frac{1}{2}) based upon an expansion of the free energy28,30. However, fitting the CTAX and Raman SnS (SnSe) data for the soft mode gives α = 0.23 and 0.24 (0.23), respectively. Furthermore, the c position of Sn atoms in the unit cell (offset from Cmcm mirror plane) can be identified as the order parameter of the Pnma–Cmcm transition. Fitting this order parameter (data taken from ref. 13) gives α = 0.24, in remarkable agreement with the soft phonon frequency. This exponent could potentially reflect the rather 2D nature of bonding, as well as long-range interactions from resonant bonding in SnS and SnSe (refs. 31,32).

We now compare the experimental TOc evolution with simulations combining AIMD with TDEP26 to incorporate anharmonic renormalization. As a baseline, phonon dispersions of the Pnma phase at 0 K were simulated within the harmonic approximation (Fig. 2a). Further, AIMD was performed at 600 K and combined with TDEP to extract renormalized second-order force-constants (Φ(2)) and phonon dispersions. A linear interpolation between the harmonic and renormalized Φ(2) was subsequently used to approximate the temperature dependence of Φ(2) in the Pnma phase between 0 K and 600 K and extrapolated to the phase transition. The resulting temperature dependence of the TOc mode could be fitted with the same power law as above, resulting in a transition temperature TC = 624 K and α = 0.28, the latter of which is in reasonable agreement with experimental results. The underestimation of TC is consistent with recent anharmonic simulations on SnSe (ref. 18). Simulation temperatures were scaled to match the experimental TC = 880 K, therefore 600 K is scaled to 846 K when comparing with experiments. The resultant soft mode energies are plotted in Fig. 3c.

### Anharmonic renormalization and thermal conductivity

Dramatic shifts in frequency and linewidth are observed for the soft mode. However, this giant anharmonicity, which is reflected in the changing Φ(2) in Fig. 4d, is not limited to the soft mode. Supplementary Figs. 2225 also show considerable phonon energy shifts and broadenings for different branches (s) and wavevectors (q). Extracting phonon energies allows for the calculation of group velocities (vg), which is a crucial factor in the lattice thermal conductivity:

$${{boldsymbol{kappa }}}_{{rm{lat}}}=sum _{q,s} {C}_{{rm{v}},q,s} {v}_{{rm{g}},q,s}otimes {v}_{{rm{g}},q,s} {tau }_{q,s},$$

(1)

where Cv is the specific heat at constant volume. We extracted vg for the TAc branch at Q = (1.5,0,1), (1.8,0,1) from the CTAX data over a large temperature range, as shown in Fig. 4a. From 295 to 864 K, a drastic reduction of 35% and 77% is observed in vg at q = (0.2, 0, 0) and (0.5, 0, 0), respectively. By contrast, vg for the LOc mode along [001] increases for small q as shown in Fig. 4b. We note that this complex and anisotropic vg behavior is qualitatively reproduced in the simulations. However, based on our simulations, ({bar{v}}_{{rm{g}}}) averaged over all branches and throughout the Brillouin zone remains almost constant: ({bar{v}}_{{rm{g}}}=) 948 ms−1 at 0 K and 973 ms−1 at 600 K. This anisotropic temperature dependence of the group velocities requires carefully investigating the effects of renormalization on κlat, which we now discuss.

Phonon renormalization with temperature is currently neglected in the majority of first-principles thermal conductivity simulations. In a widely embraced approach pioneered by Broido et al.33, the phonon Boltzmann transport equation (BTE) is solved using first-principles simulations of second and third-order force-constants, Φ(2) and Φ(3) based on low-order perturbation theory33,34,35,36. The scattering rates (Γ) in the BTE are usually obtained from three-phonon scattering processes (sometimes extended to higher order), but energies and group velocities are typically from the harmonic dispersions. While this formalism has been successful in analyzing many materials and has led to important predictions that were later verified, such as ultra-high thermal conductivity in boron arsenide37,38,39, the range of anharmonic strength over which it remains valid is poorly established.

Here, we show that anharmonic phonon renormalization, achieved through using the TDEP method combined with temperature-dependent AIMD, can have large effects on thermal conductivity. The anharmonic renormalization in SnS leads to a strong reduction in the calculated κlat of 39% just by changing Φ(2) (0 K harmonic →  AIMD-TDEP at 600 K), whilst keeping Φ(3) unchanged. Group velocities directly depend on the phonon frequencies, which can be obtained from either the harmonic or renormalized phonon models. In addition, scattering rates Γs for different phonon branches s also depend on energies (Eq. (2))40, although more indirectly, through the scattering phase space and the Bose–Einstein occupation factor (Eq. (3)).

$${Gamma }_{s}={frac{hslash pi }{16}}sum_{{s}_{1}{s}_{2}}{big|{Phi }_{s,{s}_{1},{s}_{2}}big|}^{2}big[({n}_{{s}_{1}}+{n}_{{s}_{2}}+1) delta ({omega }_{s}-{omega }_{{s}_{1}}-{omega }_{{s}_{2}})\ + 2({n}_{{s}_{1}}-{n}_{{s}_{2}}) delta ({omega }_{s}-{omega }_{{s}_{1}}+{omega }_{{s}_{2}})big]$$

(2)

$$n({omega }_{q,s},{T}_{{rm{B}}E})={big(exp (hslash {omega }_{q,s}/{k}_{{rm{B}}}{T}_{{rm{B}}E})-1big)}^{-1}$$

(3)

The lattice thermal conductivity was calculated with almaBTE41 as well as TDEP (see “Methods”). Our κlat calculations agree well with our direction-dependent measurements of κtot on SnS crystals (details in “Methods”), shown in Fig. 4c. We note that the κlat obtained with TDEP agrees better with experimental values than the κlat obtained with almaBTE, reflecting the importance of renormalization effects. Our measured data are in good agreement with previous single-crystal measurements of SnS42, in particular showing a similar anisotropy: κb > κc > κa (Supplementary Fig. 16). The electronic component κel = κtot − κlat is expected to only contribute a few percent of the total42. Our calculations also agree with trends in previous simulations35,36, yielding κb,lat > κc,lat > κa,lat, but the absolute values differ significantly (Supplementary Fig. 15). We first consider the effect of unit cell expansion and shape at the quasiharmonic (QHA) level. Comparing the first two lines in Table 1, we see that κlat at occupation temperature TBE = 300 K is strongly suppressed when using Φ(2) calculated with an expanded unit cell in QHA (see “Methods”) compared to results from Φ(2) evaluated on the theoretical relaxed cell. In particular, the expanded cell, representative of T = 862 K, results in 62% suppression along a. Further, the anisotropy between κb,lat and κc,lat decreases, reflecting the nearly tetragonal cell shape near TC (b = 4.07 Å, c = 4.05 Å24) and capturing some of the experimental trend.

To investigate intrinsic renormalization effects, Φ(3) were combined with Φ(2) from either harmonic 0 K or 600 K TDEP calculations. These Φ(2) are then used to calculate κlat (almaBTE) at different TBE with fixed Φ(3). The strong reduction in κlat upon changing Φ(2) is found to be dominated by changes in phonon lifetimes rather than group velocities, reflecting changes in the phase space from the pronounced phonon renormalization. To establish this, we calculate the change in κlat from changing Φ(2) calculated at 0 K (harmonic) to those obtained at 600 K (AIMD-TDEP), whilst fixing TBE = 300 K (Table 1, first versus third row and Fig. 4f blue curve versus red curve). This drastically reduces the direction-averaged ({bar{kappa }}_{{rm{lat}}}) from 2.25 to 1.37 Wm−1 K−1. This reduction of 39% in ({bar{kappa }}_{{rm{lat}}}) results purely from phonon renormalization effects as Φ(3) were kept fixed. Further, we separate the effects of renormalization on vg and τ in Eq. (1), which both depend on Φ(2). A change from vg,0K to vg,600K, whilst keeping τ0K fixed, has a relatively small effect on thermal conductivity, as seen in the cumulative κlat plot in Fig. 4f. In fact, this actually acts to slightly increase κlat by 2.3%, which arises from a mixture of increasing and decreasing vg along different branches upon warming, as seen in Fig. 4a, b. Conversely, a change from τ0K to τ600K, while fixing vg = vg,0K, leads to a large reduction of 38% in κlat (dotted brown curve), which is almost equivalent to the 39% reduction found when using τ600K with vg,600K (red curve). Thus, the suppression in κlat from Φ(2) renormalization occurs mainly through an increase in scattering rates (linewidths), shown in Fig. 4e, especially pronounced for low-energy TO and TA modes below 5 meV. In Supplementary Table 1, we benchmark the calculated linewidths against those obtained from fitting Raman spectra and find good agreement. The linewidths increase across the entire spectrum when thermally renormalizing Φ(2) while keeping Φ(3) fixed (purple circles versus red crosses), illustrating the importance of anharmonic renormalization of dispersions.

Using TDEP, we computed κlat with renormalized force-constants (Φ(2) and Φ(3)) from AIMD at 300 K and 600 K, with fixed TBE = 300 K to isolate renormalization effects. Renormalizing both Φ(2) and Φ(3) suppresses κlat from 2.10 Wm−1 K−1 to 1.77 Wm−1 K−1, (a 16% decrease). Including only the effect of renormalized group velocities would actually increase κlat from 2.10 Wm−1 K−1 to 2.24 Wm−1 K−1, owing to the increased average group velocity (from 1041 ms−1 to 1117 ms−1). On the other hand, changing the phonon lifetimes from 300 K TDEP to 600 K TDEP results in a thermal conductivity suppression from 2.10 Wm−1 K−1 to 1.59  Wm−1 K−1, a 24% suppression. Thus, our TDEP simulations further support our conclusion that the change in scattering rates is the dominant renormalization effect, compared to the more minor effect of group velocity renormalization.

### Conclusion

Our results reveal drastically temperature-dependent and anisotropic lattice dynamics of SnS and SnSe. Spectacular softening is observed for extended regions of the TAc and TOc dispersions along [100] in Pnma below TC, and an extremely soft TAc branch persists in the Cmcm phase. However, renormalization effects were not limited to these dispersions and instead were found to be significant throughout the Brillouin zone leading to strong renormalization and pronounced anisotropic changes in group velocities. While not included in typical calculations based on low-order perturbation theory, similarly strong anharmonic renormalization effects are expected to occur in many materials near structural transitions, including halide perovskites, oxide ferroelectrics, or thermoelectrics near instabilities (e.g., PbTe, PbSe, GeTe, Cu2−xSe) besides SnS and SnSe. This work provides impetus for future first-principles κ simulations to consider renormalization effects and indeed to re-examine older work within this context. Furthermore, our results provide additional motivation to investigate materials exhibiting lattice instabilities in the search for ultra-low κ materials.