The origin of the lattice thermal conductivity enhancement at the ferroelectric phase transition in GeTe

The proximity to structural phase transitions in IV-VI thermoelectric materials is one of the main reasons for their large phonon anharmonicity and intrinsically low lattice thermal conductivity $\kappa$. However, the $\kappa$ of GeTe increases at the ferroelectric phase transition near $700$ K. Using first-principles calculations with the temperature dependent effective potential method, we show that this rise in $\kappa$ is the consequence of negative thermal expansion in the rhombohedral phase and increase in the phonon lifetimes in the high-symmetry phase. Negative thermal expansion increases phonon group velocities, which counteracts enhanced anharmonicity of phonon modes and boosts $\kappa$ close to the phase transition in the rhombohedral phase. A drastic decrease in the anharmonic force constants in the cubic phase increases the phonon lifetimes and $\kappa$. Strong anharmonicity near the phase transition induces non-Lorentzian shapes of the phonon power spectra. To account for these effects, we implement a novel method of calculating $\kappa$ based on the Green-Kubo approach and find that the Boltzmann transport equation underestimates $\kappa$ near the phase transition. Our findings elucidate the influence of structural phase transitions on $\kappa$ and provide guidance for design of better thermoelectric materials.

The proximity to structural phase transitions in IV-VI thermoelectric materials is one of the main reasons for their large phonon anharmonicity and intrinsically low lattice thermal conductivity κ. However, the κ of GeTe increases at the ferroelectric phase transition near 700 K. Using firstprinciples calculations with the temperature dependent effective potential method, we show that this rise in κ is the consequence of negative thermal expansion in the rhombohedral phase and increase in the phonon lifetimes in the high-symmetry phase. Negative thermal expansion increases phonon group velocities, which counteracts enhanced anharmonicity of phonon modes and boosts κ close to the phase transition in the rhombohedral phase. A drastic decrease in the anharmonic force constants in the cubic phase increases the phonon lifetimes and κ. Strong anharmonicity near the phase transition induces non-Lorentzian shapes of the phonon power spectra. To account for these effects, we implement a novel method of calculating κ based on the Green-Kubo approach and find that the Boltzmann transport equation underestimates κ near the phase transition. Our findings elucidate the influence of structural phase transitions on κ and provide guidance for design of better thermoelectric materials.
Recent computational work has predicted that driving IV-VI materials closer to the ferroelectric phase transition via strain or alloying can lead to a drastically lower lattice thermal conductivity [22,23]. Under the assumption of the displacive phase transition, it was found that coupling between soft transverse optical (TO) modes and the heat carrying acoustic modes is the main reason for the κ reduction. At the displacive transition, the frequency of the soft TO mode collapses, becoming effectively zero. Since scattering rates are inversely proportional to phonon frequencies, the lifetimes of the acoustic modes that couple to soft TO modes decrease dramatically, leading to a considerable κ reduction [22,23]. Surprisingly, experimental studies have shown that the lattice thermal conductivity increases at the ferroelectric phase transition in GeTe [25][26][27][28][29][30][31][32]. This is at odds with measurements in some other materials going through ferroelectric phase transitions, where a significant decrease in κ is observed [33]. The reason for the anomalous behaviour of κ at the phase transition in GeTe remains unknown. Understanding the microscopic origin of the κ increase at the ferroelectric phase transition in GeTe may lead to design of improved thermoelectric materials.
Here we study how driving GeTe near the ferroelectric phase transition via temperature affects its lattice thermal conductivity, making no assumptions about the nature of the phase transition. Unlike the previous work [22,23], we calculate interatomic force constants at different temperatures using the state-of-the art, temperature dependent effective potentials (TDEP) method [34][35][36]. We find that the increase of the lattice thermal conductivity of GeTe at the phase transition in the rhombohedral phase comes from negative thermal expansion that enhances the phonon group velocities. In the cubic phase the phonon lifetimes increase, leading to even more substantial increase of κ. Large anharmonicity of phonon modes minimizes the phonon lifetimes at the phase transition in the rhombohedral phase and leads to non-Lorentzian power spectra of phonon modes. We implement a new method of calculating κ that includes these non-Lorentzian lineshapes of the phonon power spectra near the phase transition. This approach further increases κ at the phase transition, which can be attributed to further softening of the phonon frequencies due to phonon-phonon interaction.

Results and discussion.
TO mode softening at the phase transition. GeTe is a ferroelectric material, exhibiting a spontaneous polarization below 600 − 700 K [37][38][39][40] (the critical temperature strongly depends on the free charge carrier concentration). This occurs due to a slight offset of the Te sublattice along one of the body diagonals of the rocksalt structure. At temperatures higher than 600 − 700 K, GeTe transforms to the rocksalt structure, losing its ferroelectric nature [41][42][43][44]. Below 600−700 K, the GeTe structure can be described by the following set of lattice vectors: where a is the lattice constant, b = 2(1 − cos θ)/3, c = (1 + 2 cos θ)/3, and θ is the angle between the primitive lattice vectors. The atomic positions in this structure are taken to be: Ge (0.0,0.0,0.0) and Te (0.5 + µ, 0.5 + µ, 0.5 + µ) in reduced coordinates. If the phase transition is displacive, as assumed in Refs. [22,23], the angle θ becomes 60 • at the phase transition, while the interatomic displacement parameter µ becomes zero [41,42]. In this type of the phase transition, the harmonic frequency of the soft mode also collapses to zero [45]. Here we define the harmonic frequency of the phonon mode as the square root of the eigenvalue of the dynamical matrix for that phonon mode at a certain temperature. In the order-disorder phase transition, both the harmonic frequency of the soft mode and the local interatomic displacement are non-zero [43,44,46]. It is still under debate which type of the phase transition occurs in GeTe [41][42][43][44][45].
We have calculated the phonon dispersion of GeTe at different temperatures using the TDEP method. This method allows us to calculate harmonic frequencies of phonon modes at different temperatures. Fig. 1 shows the harmonic frequencies of the two TO modes at the Brillouin zone centre in GeTe. The phase transition temperature obtained in our calculations is T C ≈ 634 K (see Methods). The harmonic frequencies of both phonon modes are strongly renormalized by anharmonic interaction and lattice thermal expansion. The harmonic frequencies of these phonon modes soften drastically at the phase transition, but they do not become zero, indicating that the phase transition in GeTe might not be of the displacive type.
The observed phonon frequency in experiments is not the harmonic frequency, but what we will call in the rest of the paper the anharmonic frequency i.e. the peak of the phonon mode power spectrum. Blue squares in Fig.  1 represent the computed anharmonic frequencies of both a The lower transverse optical (TO) mode (TO1) and b the higher TO mode (TO2) at the zone centre. Red circles and blue squares represent the calculated harmonic and anharmonic phonon frequencies, while other symbols represent the experimental results from Refs. [43,45,47]. The harmonic frequency is the square root of the eigenvalue of the dynamical matrix for q = 0 at a particular temperature, while the anharmonic frequency is the peak of the phonon mode power spectrum. The black vertical line represents the phase boundary between the rhombohedral and rocksalt structures in our calculations (≈ 634 K).
TO modes at the zone centre, which do fall to zero at the phase transition, in contrast to the harmonic frequencies.
Our results thus suggest that the observation of phonon mode softening is not a conclusive proof of the displacive type of the phase transition, as previously argued in the case of GeTe [45].
Our computed anharmonic TO frequencies at the zone centre agree fairly well with those measured in experiments, see Fig. 1. This agreement highlights the accuracy of the TDEP method even for the challenging cases of materials undergoing structural phase transitions. We note that the critical temperature in Ref. [45] is around 600 K, in Ref. [43] approximately 700 K and in Ref. [47] 650 K. This difference in the calculated and measured critical temperature is expected since the critical tem-perature strongly depends on the number of free charge carriers [48,49].
Phonon spectral function. Harmonic frequencies are a valid description of lattice dynamics only in the absence of phonon-phonon interaction. In the inelastic neutron scattering experiments that measure phonon spectral functions, harmonic phonons would produce zero linewidth signals, revealing infinitely long lived quasiparticles. However, in real materials phonons interact with each other, thus broadening phonon spectral functions, with linewidths inversely proportional to phonon lifetimes. This effect can be described using the concept of phonon self-energy that quantifies the strength of phonon-phonon interaction [50][51][52].
The probability of an incoming neutron to interact with a phonon system acquiring/losing energy Ω and momentum q is proportional to the spectral function [50][51][52]: where u q,s is the phonon displacement operator, ω q,s is the harmonic frequency of the phonon mode with wave vector q and phonon branch s, β is 1/k B T with k B being the Boltzmann constant and T the temperature, and is the reduced Planck constant. ∆ q,s and Γ q,s are the real and imaginary part of the phonon-self energy, respectively. In a weakly anharmonic material, this expression can be reduced to a Lorentzian with a halfwidth of Γ q,s and the position of the peak at ω q,s + ∆ q,s . In this case the phonon lifetime can be calculated as τ q,s = 1 2Γ q,s . However, in materials with strong anharmonicity, phonon spectral functions can exhibit exotic behaviour with satellite peaks, shoulders etc. [18,20,[53][54][55][56].
First we show the phonon spectral function of GeTe near the phase transition (631 K) along a high symmetry path, see Fig. 2. The black lines represent the harmonic frequencies ω q,s calculated at 631 K and the cyan dashed lines are the harmonic frequencies at 300 K. The harmonic frequencies of optical modes soften from 300 to 631 K. The acoustic mode frequencies soften as well in the F -Γ and the in plane (Γ -X) directions. This is the consequence of the overall softening of the second order force constants with temperature. On the other hand, the acoustic modes along the Γ -Z direction (the direction along the trigonal axis) stiffen because of the negative thermal expansion, as we will discuss in the next section. The spectral function shows further softening of the phonon frequencies due to anharmonic phononphonon interaction. As expected, the TO phonon at Γ can not be clearly resolved in the graph of the phonon spectral function.  Figure 3 shows the unusual features of the spectral function for the second transverse optical mode (TO2) of GeTe at the zone centre for several different temperatures. A non-Lorentzian behavior of the phonon spectral function is evident even at 300 K, very far from the phase transition. The broadening of the power spectrum is large, revealing the short lifetime of this phonon mode. The distortion of the power spectrum is stronger at 625 K, with a very large shift of the peak of the spectral function compared to the harmonic frequency. For temperatures near the phase transition, the power spectrum peaks around 0 THz as expected at the phase transition (see Fig. 1). At temperatures higher than the phase transition temperature, the peak of the power spectrum is at non-zero frequencies, but is still strongly renormalized compared to the harmonic frequency.
Non-Lorentzian shapes of the phonon power spectra of the TO modes in GeTe at the phase transition are the consequence of the coupling of these phonon modes to the entire phonon bath, rather than coupling to specific phonons. We test this by calculating the power spectrum of the TO2 phonon mode disregarding the coupling to specific phonon branches. We find that the change in the power spectrum does not substantially vary depending which phonon branch we disregard. A similar behaviour can be seen in the spectral function of the TO1 mode (see Supplementary Note 1).
Lattice thermal conductivity in the Boltzmann transport approach. Next we calculate the lattice thermal conductivity of GeTe for a range of temperatures including both rhombohedral and rocksalt phases (see Fig. 4), combining the TDEP method with the Boltzmann transport approach. Overall, the lattice thermal conductivity is inversely proportional to temperature, as a result of the linear dependence of phonon populations with temperature above the Debye temperature ( 200 K for GeTe). The calculated κ deviates from the 1/T law near the phase transition, where there is a large κ increase at the phase transition and in the cubic phase. At high temperatures, the κ in the cubic phase regains the 1/T dependence.
There is an anisotropy in the lattice thermal conductivity in the rhombohedral phase ( Fig. 4), as a consequence of the van der Waals gaps formed due to the Te sublattice offset. The direction perpendicular to the van der Waals gaps (i.e. parallel to the trigonal [111] axis) has weaker bonding, leading to the lower phonon group velocities and κ in that direction. Anisotropy of the lattice thermal conductivity disappears close to the phase transition and is not present in the rocksalt phase.
The agreement between the computed and experimental κ values is very good in the whole temperature range of the rhombohedral phase, especially if we consider the average κ (dashed line in Fig. 4). Almost all experiments show an increase in the lattice thermal conductivity in the vicinity of the phase transition [25][26][27][28][29][30][31][32], similarly to our results. The discrepancy between our results and experiments increases as we get closer to the phase transition and particularly in the cubic phase. Our results consistently overestimate κ compared to experiments. In the rhombohedral phase we would expect scattering from lattice imperfections, such as ferroelectric domain walls [57][58][59][60][61], to further reduce the computed lattice thermal conductivity, but this is not the case in the cubic phase. The discrepancies could also stem from ommiting the higher order terms in the Taylor expansion of the interatomic forces (we include the second and third order terms only). However, to the lowest approximation, the fourth order anharmonic terms would only affect the real part of the self-energy [51] and not imaginary part, which would mean that the phonon lifetimes should remain unchanged. Additionally, since the force constants in TDEP are obtained through a fitting procedure, the fourth order force constants should be smaller than the third order ones. We checked the influence of the fourth order anharmonicity on the real part of the self-energy and found that it is small even for highly anharmonic soft modes at the phase transition. Therefore, higher order anharmonicity is not the reason for the differences in the calculated and experimental κ values.
Experimental values of lattice thermal conductivity are usually extracted from the total thermal conductivity measurements using the Wiedemann-Franz law to eliminate the electronic contribution to the thermal conductivity, whose validity at structural phase transitions is not well understood. Additionally, most references use the single parabolic band Kane model to extract the Lorenz factor from measurements of the Seebeck coefficient, which is not appropriate in GeTe due to the intrinsically complicated Fermi surface [31]. Such lattice thermal conductivity values can differ widely near structural phase transitions, and sometimes an increase in the total thermal conductivity is assigned to the electronic contribution. Here we show that the increase in the thermal conductivity of GeTe at the phase transition, at least partially, comes from the lattice thermal conductivity. The difference between our theoretical and experimental results in the cubic phase might be due to an inaccurate estimation of the electronic contribution to the total thermal conductivity in experiments. Measuring the thermal conductivity of GeTe in an applied magnetic field (to exclude the electronic thermal conductivity) would test our predictions of the increased lattice thermal conductivity at the phase transition.
To understand the anomalous behaviour of κ near the phase transition, we calculate the average phonon lifetimes and group velocities at different temperatures, Fig. 3. For example, the average values of the phonon lifetimes in the vicinity of the phonon frequency ω 0 are given as:τ where the sum goes over all phonon modes λ, and σ is the smearing parameter taken to be σ where ω D is the Debye frequency and N is the number of ω 0 frequencies.
The phonon group velocities of GeTe are mostly independent of temperature, except very close to the phase transition, see Fig. 5 a. In this temperature region (600-− 675 K), there is an increase in the phonon group velocities across most of the frequency range and most noticeably for phonons between 1 and 3 THz. This is the frequency region that contributes most to the thermal conductivity (see Supplementary Note 2). We thus conclude that the anomalous increase of the thermal conductivity at the phase transition is partially due to this rise in the phonon group velocities. We find that the increase in group velocities originates from the lattice contraction near the critical temperature [41,42,62]. This can be understood from the observation that the increase of phonon group velocities happens most prominently along out of the plane direction, where negative thermal expansion occurs [62].
Phonon lifetimes in weakly anharmonic materials usually follow the 1/T law, similar to thermal conductivity. This is the case in our calculations in the rhombohedral phase far from the phase transition. At the phase transition, however, the phonon lifetimes of acoustic modes decrease more than expected from the 1/T scaling. This is a signature of stronger anharmonicity of acoustic phonon modes closer to the phase transition in the rhombohedral phase. Optical modes have more complicated behaviour. While soft transverse optical modes near the zone centre have much lower lifetimes at the phase transition, this is not the case for transverse optical modes in the rest of the Brillouin zone. Even more intriguing is the behaviour of longitudinal optical modes. At low temperatures, they usually have frequency independent lifetimes, but at the phase transition they decrease exponentially with frequency. In the cubic phase there is a substantial increase in the phonon lifetimes with respect to the rhombohedral phase(see Fig. 5 b). The phonon lifetimes at 637 K are larger in most of the frequency range compared to the temperatures closest to the phase transition in the rhombohedral phase (625 K and 631 K). Interestingly, the phonon lifetimes at 675 K are larger than the phonon lifetimes at 300 K rescaled by temperature (i.e. by 675 K/300 K), revealing stronger intrinsic anharmonicity of the rhombohedral phase. Third order force constants are much stronger in the rhombohedral phase, even for very similar temperatures and structures (631 K vs 637 K), see Supplementary Note 3.
In conclusion, the increase of the lattice thermal conductivity near the phase transition can be attributed to two phenomena, depending on the structure of GeTe. In the rhombohedral phase, the increase in the thermal conductivity is due to negative thermal expansion that causes phonon group velocities to increase, increasing phonon mean free paths. On the other hand, in the cubic phase, phonon lifetimes increase dramatically due to lower intrinsic anharmonicity of this phase.
Lattice thermal conductivity using the Green-Kubo method. The non-Lorentzian behaviour of the phonon spectral function raises the question whether the Boltzmann transport equation employed in the calculation of the κ values in Fig. 4 is valid close to the phase transition. It is assumed in the derivation of the Boltzmann equation that phonons are well defined quasiparticles with unique frequencies and lifetimes. This implies that their spectral weights are the Lorentzian functions centred at the harmonic frequencies and the widths equal to the phonon lifetimes. However, evidently this does not hold close to the phase transition, see Fig. 3.
We include the effect of non-Lorentzian lineshapes on lattice thermal conductivity following the Green-Kubo approach derived in Refs. [63,64]. In this approach, the heat current in the Cartesian direction j is defined as [63,65,66]: where v j q,s is the group velocity of the phonon mode with wave vector q and branch s, V is the volume of the unit cell, N is the number of q points, and p q,s (t) and u q,s (t) are the Fourier transforms of the phonon momentum and position operators. The lattice thermal conductivity tensor is obtained by employing the Green-Kubo relation [67]: The spectral theorem is then used to relate the heat current autocorrelation function in the equation above to the one particle retarded Green function [50,64,[68][69][70]. The final expression for the thermal conductivity in this approach is [63,64]: where q,s is: v q,s,s is the generalized phonon group velocity [65]: where X q,s is the eigenvector of the phonon with wave vector q and branch s, Φ a,b ( R) is the force constant between atoms with masses m a and m b and the vector distance R.
We can separate Eq. 6 into two parts. The first part is diagonal (s = s ). In the limit of small anharmonicity (∆ q,s = 0 and Γ q,s ω q,s ), this part reduces to the standard solution of the Boltzmann equation in the single relaxation time approximation. The second, non-diagonal part (s = s ), can be reduced in the limit of small anharmonicity to the expressions similar to the ones given in Refs. [71,72]. The non-diagonal contribution to the lattice thermal conductivity will become prominent only if there is a substantial overlap in the spectral functions of two phonon modes with the same wave vector. This is only true in the case of strong anharmonicity or when spectral functions broaden due to disorder. We have implemented the expression given by Eq. 6 with the TDEP method. Figure 6 a shows the difference between our results obtained using Eq. 6 and BTE (Fig.  4). We can see that the BTE underestimates the thermal conductivity in the whole temperature range. Additionally, we can see that the underestimation is not large, around 10% even at high temperatures. Overall, the difference scales linearly with temperature. We can also see that the difference is largest at the phase transition, which is expected considering large deviations from the Lorentzian shape of the phonon spectral functions in this region (see Fig. 2). The contribution of the non-diagonal part of the lattice thermal conductivity is comparable to the overall enhancement of the diagonal part due to non-Lorentzian shapes of the phonon spectral functions.
To understand the reason for the increased difference in κ at the phase transition obtained by the standard BTE method and using the Green-Kubo relation (Eq. 6), we show the phonon lifetimes of GeTe calculated using the two methods in Fig. 6 b. We define the phonon lifetimes in the Green-Kubo method as: where c q,s is the harmonic heat capacity of the phonon mode ( q, s). The increase of the phonon lifetimes in the Green-Kubo method is visible in the whole Brillouin zone. It is, however, most prominent in the region of soft phonon modes, where the phonon lifetimes increase by a factor of 100. The increase in phonon lifetimes mostly comes from the shifts of the peaks of the spectral functions, which increases the heat capacity of the phonon modes compared to the harmonic heat capacity.
Discussion. Here we highlight the advantages of the Green-Kubo method we implemented here over the standard approaches for computing κ in strongly anharmonic materials. Unlike the Boltzmann transport equation, the Green-Kubo method accounts for non-Lorentzian shapes of phonon spectral functions. It is possible to include these effects also by running long molecular dynamics (MD) simulations. Compared to MD, the Green-Kubo method presented here is faster and easier to converge, which is particularly important if these methods are combined with first principles calculations. The results obtained using this approach are easier to interpret and analyse compared to traditional MD approaches. Unlike MD simulations, the Green-Kubo method uses the Bose-Einstein statistics for phonons. Additionally, this method gives the non-diagonal part of lattice thermal conductivity, accounting explicitly for the whole phonon power spectra. The method of Refs. [71,72] includes only the values of the phonon self-energy at the harmonic frequency in the evaluation of the non-diagonal contribution to κ and is not applicable to strongly anharmonic materials.
In conclusion, we have performed a detailed first principles study of the lattice thermal conductivity κ of GeTe close to the ferroelectric phase transition. The harmonic frequencies of the soft modes, although dramatically softened, do not become zero at the phase transition. On the other hand, strong anharmonicity causes the spectral functions of the soft modes to collapse and effectively peak at zero frequency. Strong anharmonicity minimizes the acoustic phonon modes lifetimes at the phase transition. However, we calculate an increase in the lattice thermal conductivity at the phase transition, in agreement with experiments. In the rhombohedral phase, this effect is due to negative thermal expansion that increases phonon group velocities. In the cubic phase, the increase in κ is primarily driven by increased phonon lifetimes due to smaller anharmonicity of phonon modes compared to the rhombohedral phase. We implement a novel approach to compute lattice thermal conductivity that includes the observed non-Lorentzian power spectra of phonon modes. Using the new approach, we find that the calculated κ increases even further at the phase transition, which is the consequence of larger phonon populations due to softening of the phonon modes caused by phonon-phonon interaction.

Computational methods.
We calculate the lattice thermal conductivity of GeTe from first principles using density functional theory (DFT) and the Boltzmann transport equation (BTE). In this approach, the lattice thermal conductivity tensor is given as [63]: where i and j are the Cartesian directions, N is the number of q points, V is the unit cell volume, c q,s is the phonon mode heat capacity and v i q,s is the group velocity of the phonon mode ( q, s) in the direction i. The relaxation time of the same mode is τ q,s = 1/2Γ q,s , where the imaginary part of the phonon self-energy due to threephonon scattering is given as [50]: Here λ is a short hand notation for ( q, s) and phonon momentum is conserved in the three-phonon processes above. The three phonon matrix element Φ λλ λ is calculated using: where m i is the mass of atom i, iα λ is the component α of the eigenvector for mode λ and atom i, and r i is the vector associated with atom i. ω λ is the harmonic frequency of the phonon mode λ, and Φ αβγ ijk is the third order interatomic force constant.
To get temperature dependent interatomic force constants, we use the temperature dependent effective potential (TDEP) method [34][35][36]. This approach employs a fitting procedure of the second and third order force constants to DFT forces on forces sampled along a molecular dynamics (MD) trajectory. To perform MD simulations, we developed a very accurate interatomic potential based on the Gaussian Approximation Potentials scheme [73,74] (see Supplementary Note 4). We calculated the structural parameters of GeTe at different temperatures using MD. We ran NVT ensemble on these structures to obtain atomic configurations. After 50 ps of equilibration we sampled 24 configurations on the 300 ps long trajectory. We used a 512 atom supercell to converge phonon properties with respect to the second and third order force constants cutoff (12 and 8 Å, respectively). We extracted selected configurations and carried out density functional theory calculations on them to obtain forces.
DFT calculations were performed using the ABINIT software package [75,76]. We use a generalized gradient approximation with the Perdew-Burke-Ernzerhof parametrization (GGA-PBE) [77] for the exchangecorrelation functional and the Hartwigsen-Goedecker-Hutter (HGH) pseudopotentials [78]. Wave functions are represented in a plane wave basis set with the cutoff of 16 Ha, and the Γ point is used for sampling of electronic states. Sampling of phonon states is carried out using a 30 × 30 × 30 q-point grid.

Code availability.
The code that implements the temperature effective potential (TDEP) method is available from Olle Hellman upon reasonable request. The additional data processing scripts and codes are available from the corresponding authors. GAP software is available for non-commercial use from www.libatoms.org.

Data availability.
The authors declare that the data supporting the present work is available from the corresponding authors upon reasonable request.  figure 1 shows the spectral function of the zone centre second transverse optical (TO2) phonon mode at 631 K, computed including and excluding transverse acoustic-TO2 mode coupling. When we neglect the coupling of the TO2 mode with transverse acoustic (TA) modes in our calculation, the peak of the spectral function shifts closer to the harmonic frequency. A similar behaviour is observed if we disregard the coupling of the TO2 zone centre mode with TA1 and TO2 phonon branches. This illustrates that although coupling to TA modes is strong, it is not the sole source of the exotic behaviour of the soft TO2 phonon power spectrum. In Supplementary figure 1 we show the types of coupling that bring the phonon power spectrum closest to the Lorentzian shape.

Acknowledgements.
The spectral function of the zone centre first transverse optical mode (TO1) in GeTe is illustrated in Supplementary figure 2 a for several temperatures. We can see a secondary peak in the TO1 spectral function at 625 K, indicating that anharmonicity of this mode is even larger than that of the TO2 mode. This is somewhat to be expected since the harmonic frequency of the TO1 mode is lower than that of the TO2 mode, leading to larger coupling to the rest of the modes. At 631 K the spectral function of this mode has a peak at almost 0 THz. In the rocksalt phase, the TO1 and TO2 modes are degenerate and have the same lineshapes. In Supplementary figure 2 b we show the spectral function of the zone centre TO1 mode at 631 K with full anharmonicity, excluding the coupling of TO1 mode with transverse acoustic modes and excluding the coupling to the TA2 mode and itself. We can see that the Lorentzian shape of the spectral function is not regained after excluding different types of interaction and thus we conclude that the non-Lorentzian shape of the TO1 mode power spectrum and the softening of the TO1 frequency is due to coupling to the entire phonon bath, rather than a particular phonon branch. To analyze the contribution of specific phonon modes to the total lattice thermal conductivity, we calculated the spectral thermal conductivity at different temperatures by convolving the total lattice thermal conductivity with a Gaussian of an appropriate width (see the main part for more details). We have found that acoustic modes give the dominant contribution to the lattice thermal conductivity, as one would expect (see Supplementary Fig. 3). At temperatures near the phase transition, the overall contribution of transverse modes (acoustic and optical) to the lattice thermal conductivity of GeTe diminishes. At the phase transition, the largest contribution to κ comes from the phonon modes in the frequency range between 1 and 3 THz. This is the frequency region that shows the most prominent enhancement of phonon group velocities due to negative thermal expansion in the rhombohedral phase.  figure 4 shows the average phonon lifetimes of GeTe at different temperatures scaled by T /300 K, where T is the temperature. This enables us to see whether the decrease in the phonon lifetimes near the phase transition is due to phonon populations or increased anharmonicity of the material. Phonon lifetimes scale inversely with temperature for the temperatures far from the phase transition temperature (T C ≈ 634 K), which indicates that anharmonicity is not increased for those temperatures. However, close to the phase transition (631 K) in the rhombohedral phase, we can see a dip in the phonon lifetimes for most of the frequency range, revealing increased anharmonicity near the ferroelectric phase transition in the rhombohedral phase.
In the cubic phase, however, we see an increase in the scaled phonon lifetimes compared to the rhombohedral phase. Further from the phase transition (675 K), the scaled lifetimes are larger than the scaled phonon lifetimes at 300 K, which we attribute to lower intrinsic anharmonicity of the cubic phase.
To investigate in detail the reason for this behaviour of phonon lifetimes, we checked the change in anharmonic force constants with temperature (see Supplementary Fig. 5). To compare anharmonic force constants, we define a norm of the anharmonic force constant matrix as: N AF C (i, j, k) = α,β,γ |Φ αβγ ijk |, where i, j, k denote atoms in triplet and α, β, γ are the Cartesian directions. Surprisingly, most of the anharmonic force constants decrease with temperature, in both rhombohedral and cubic phases. There is a sudden drop in the norm of the anharmonic force constants at the phase transition. The anharmonic force constants are drastically smaller in the cubic phase, explaining higher scaled phonon lifetimes in this phase. However, the temperature dependence of anharmonic force constants does not explain the increased anharmonicity of the transverse acoustic modes at the phase transition in the rhombohedral phase (see Supplementary  Figure 4). It is the scattering phase space (SPS) that increases at the phase transition appreciably, which leads to higher anharmonicity in the rhombohedral phase. We show this in Supplementary Figure 6 by calculating the scattering phase space for phonons at different temperatures. We do this by setting the matrix element Φ λλ λ in Eq. (11) of the main part to 1 and calculating the imaginary part of self-energy at the harmonic frequency for each phonon mode. The results are then scaled by temperature to minimize the effect of phonon population on the calculated value of the scattering phase space. We can notice that for the phonons in the frequency re- gion around 1 THz the SPS increases dramatically near the phase transition, which explains the lower contribution of transverse optical modes to total κ at the phase transition compared to 300 K (see Supplementary Figure 3) and prominent dip in the scaled phonon lifetimes for this frequency region (see Supplementary Figure 4). Phonons in the frequency region around 2 THz have a smaller SPS, which again is consistent with the results presented in Supplementary Figure 4 (the scaled phonon lifetimes at the phase transition and 300 K are comparable in this frequency region). However, this can not be the reason for the increased lattice thermal conductivity at the phase transition, because this effect is noticeable only due to the temperature scaling of phonon lifetimes and SPS and does not exist if one takes phonon populations into account (see Fig. 5 b of the main part). Finally, the available scattering phase space of longitudinal optical (the highest frequency) phonons increases dramatically with temperature, which explains their unusual frequency dependence.
In the cubic phase the SPS of the majority of phonons (except LO phonons) increases linearly with temperature, which balances out the decrease in the anharmonic force constants and leads to almost constant phonon lifetimes with temperature. The decrease in κ in the cubic phase between 700 K and 800 K comes primarily from a decrease in phonon group velocities.

Supplementary note 4: Gaussian Approximation Potential
To accurately sample atomic displacements at a certain temperature for the temperature dependent effective potential (TDEP) method, one must perform molecular dynamics (MD) simulations. Most commonly, one would run ab initio MD to achieve this. However, ab initio MD is extremely computationally expensive and the accuracy of such calculations is typically reduced in order to speed up the computation of forces and energies. Another important drawback of ab initio MD is that it does not scale linearly with the number of atoms, so carrying out these calculations on large supercells needed to reach convergence in GeTe would be very expensive.
To overcome these difficulties, we developed an accurate interatomic potential based on the Gaussian Approximation Potential (GAP) framework using ab initio calculated forces for different atomic configurations in different GeTe supercells. The details of the fitting of this interatomic potential will be the focus of the future publication. To check the validity of this interatomic potential for sampling atomic configurations through MD simulations, we compared the phonon band structures calculated using forces from DFT and GAP (see Supplementary Fig. 6 a). The comparison between these two methods is very good.
To further check the appropriateness of GAP for MD simulations, we compared the errors from GAP and DFT. We calculate forces on a set of structures and atomic configurations using DFT with fully converged parameters. Then we calculate the same forces using GAP and find the differences with respect to DFT. We then bin the errors and present the results in Fig. 6 b. We highlight that the structures taken in this study are not in the training set used to obtain the GAP potential which makes the low value of errors in forces even more remarkable.
As we previously said, to actually run ab initio MD simulations, one would have to reduce the accuracy of DFT calculations. To mimic this, we calculated forces on 216 atom supercells of GeTe using a 2x2x2 k-point grid and the energy cutoff of 16 Ha for plane waves. We consider this calculation as high accuracy. To perform low accuracy calculations, we chose an 1x1x1 k-point grid and the energy cutoff of 10 Ha for the plane wave expansion of electronic wave functions. We then take the differences in the forces in these two calculations and bin them in a similar manner to the previous case with the GAP potential (see Supplementary Fig. 7 b). We can see that the errors from the GAP potential are comparable to the errors obtained from the reduced accuracy first principles calculations, which gives us confidence in the sampling method we used. b The force error probability distribution function for GAP and DFT with reduced accuracy (see text for more details).