Microscopic mechanism of unusual lattice thermal transport in TlInTe$_2$

We investigate the microscopic mechanism of ultralow lattice thermal conductivity ($\kappa_l$) of TlInTe$_2$ and its weak temperature dependence using a unified theory of lattice heat transport that considers contributions arising from the particle-like propagation as well as wave-like tunneling of phonons. While we use the Peierls-Boltzmann transport equation (PBTE) to calculate the particle-like contributions ($\kappa_l$(PBTE)), we explicitly calculate the off-diagonal (OD) components of the heat-flux operator within a first-principles density functional theory framework to determine the contributions ($\kappa_l$(OD)) arising from the wave-like tunneling of phonons. At each temperature, T, we anharmonically renormalize the phonon frequencies using the self-consistent phonon theory including quartic anharmonicity, and utilize them to calculate $\kappa_l$(PBTE) and $\kappa_l$(OD). With the combined inclusion of $\kappa_l$(PBTE), $\kappa_l$(OD), and additional grain-boundary scatterings, our calculations successfully reproduce the experimental results. Our analysis shows that large quartic anharmonicity of TlInTe$_2$ (a) strongly hardens the low-energy phonon branches, (b) diminishes the three-phonon scattering processes at finite T, and (c) recovers the weaker than T$^{-1}$ decay of the measured $\kappa_l$.


INTRODUCTION
Crystalline semiconductors with ultralow lattice thermal conductivity (κ l ) are important for effective utilization and management of thermal energy in highperformance thermoelectrics [1][2][3][4], photovoltaics [5][6][7][8], thermal barrier coatings [9], and thermal data storage devices [10]. Although it is quite natural for compounds with complex crystal structures having large unit cells [11,12] and heavy atoms to exhibit low-κ l , relatively simple crystalline materials having small unit cells and even with light atoms [13] could also possess ultralow-κ l due to the presence of rattler atoms [14], strong lattice anharmonicity [15][16][17][18][19], or bonding heterogeneity [20][21][22][23]. Investigation of the microscopic mechanism behind the ultralow-κ l that often approaches the glassy limit in ordered compounds is not only fundamentally interesting, but it also helps to unravel the complex correlation between the crystal structure, bonding, and anharmonic lattice dynamics. Results of such investigations provide new criteria to find hitherto unknown low-κ l materials as well as paves the way to engineer the heat transport properties in already known compounds. Despite significant research efforts, a comprehensive theoretical understanding of the mechanism behind extremely poor heat transport in ultralow-κ l materials, which approach the limit of their theoretical minimum (κ min ) has remained challenging [24][25][26]. This is in part due to the fact that many of these materials are often so highly anharmonic that a harmonic description of the phonon frequencies fails to describe the lattice dynamics of the compounds correctly. Hence, a finite temperature treatment of the phonon modes becomes necessary to * koushik.pal.physics@gmail.com † c-wolverton@northwestern.edu account for the renormalization of the phonon frequencies arising from the temperature induced anharmonic effects. In some cases, the mean free paths of the phonon modes approach the smallest atomic distance in the crystal, leading to a breakdown of the conventional particlelike description of phonons towards a glass-like thermal conductivity [27]. However, recent theoretical and computational developments [28][29][30][31][32][33] have enabled the treatment of phonons at finite temperature considering anharmonicity arising from high-order phonon-phonon interactions, and examination beyond the particle-like description by including the contributions arising from the wave-like tunneling of phonons [32,33]. In recent theoretical studies, it was shown that the ultralow κ l in many well known family of crystalline solids such as clathrates [30], double-halide perovskites [34], tetrahedrites [35], Tl 3 VSe 4 [36] can only be explained successfully if higherorder anharmonic phonon-phonon interactions are taken into account in the description of their lattice dynamics and phonon transport [29,37].
In this work, we develop a microscopic understanding and uncover the physical principles underlying the unusual lattice thermal transport properties of TlInTe 2 which exhibits ultralow-κ l that is close to its theoretical minimum [38,39] and shows weak (milder than T −1 decay) temperature dependence in κ l . TlInTe 2 represents a class of structurally similar ABX 2 (A=Tl 1+ , In 1+ ; B=Tl 3+ , In 3+ , Ga 3+ ; X=Se 2− , Te 2− ) compounds, many of which are shown to exhibit ultralow-κ l due to the presence of rattler cation at the A-site [38][39][40][41]. The rattler cations manifest nearly dispersion-less optical phonon branches at low-energies which effectively scatter the heat-carrying phonons by creating additional phonon scattering channels [18,42,43]. Since phonon scattering rates crucially depend on the phonon frequencies [44] and the rattling phonon modes are quite sensitive to temperature, an accurate treatment of the thermal transport properties in the above compounds requires the treat- ment of phonons at finite temperature including anharmonic effects [29,45]. However, a comprehensive theoretical understanding of the lattice thermal transport including the effects of temperature-induced anharmonic renormalizations to the phonon frequencies is missing in this family of compounds. Here, we use a unified theory of lattice thermal conductivity that considers the contributions arising from the particle-like propagation as well as wave-like tunneling of phonons. We utilize the phonon frequencies that are renormalized at finite temperatures including the effects of quartic anharmonicity to calculate both these contributions to κ l . The particle-like contributions (κ l (PBTE)) which are obtained after solving the Peierls-Boltzmann transport equation (PBTE) only account for the diagonal components of the heat-flux operator [32,[46][47][48]: J α = −κ αβ l ∇T β , where ∇T is the temperature gradient and α, β are the Cartesian coordinates. We also calculate the off-diagonal (OD) contributions, κ l (OD), associated with the OD components of the heat-flux operator, which are not present in the PBTE formalism. κ l (OD) is related to the wave-like tunneling of phonons, which is the heat carried by the coupled vibrational eigenstates as a result of the loss of coherence between them [32,33,46,49]. In a recent theoretical work, Simoncelli et el. [32] combined both the pictures (i.e., particle-like and wave-like) of phonon transport in a unified theory of lattice heat transport, where the total κ l is given by κ αβ l (tot) = κ αβ l (PBTE) + κ αβ l (OD), which successfully describes the κ l of anharmonic crystals, harmonic glasses as well as complex compounds such as tetrahedrites [32,35,36]. Our calculated κ l within the PBTE formalism in conjunction with the OD contributions using the renormal-ized phonon frequencies, and additional grain-boundary scatterings successfully reproduce two sets of available experimental results [38,39] of TlInTe 2 . Our analysis reveals that TlInTe 2 possesses large quartic anharmonicity that (a) strongly hardens the low-energy rattling and other optical phonon branches with temperature, (b) diminishes the three-phonon scattering rates at finite temperatures, and (c) recovers the magnitude as well as the correct T-dependence of κ l that shows weaker than T −1 decay found in experimental measurements.

Crystal structure and anharmonicity
TlInTe 2 has a chain-like body-centered tetragonal (space group: I4/mcm) crystal structure ( Fig. 1a) with eight atoms in the primitive unit cell. Within the unit cell, In 3+ cations are covalently bonded to four Te atoms, forming the InTe 4 tetrahedra which share their edges and extend like chains along the crystallographic c-axis. The empty spaces between these chains are bridged by the Tl atoms which stabilize the structure through electron transfer from Tl 1+ cation to the [InTe 2 ] −1 anion sublattices. Analysis of the second-order interatomic force constants (IFCs) reveals that In 3+ are strongly bonded to the lattice, whereas the Tl 1+ cations are only loosely connected (see Supplementary Figure 2). In the crystal structure, Tl atoms are coordinated by eight Te atoms in a square anti-prismatic coordination environment (Fig.  2b) forming an oversized cage (also known as Thompson cube) inside which the Tl atoms rattle due to their weaker  Fig. 1c and Fig. 1d, respectively. As will be discussed later, the R LO branch is strongly affected by the temperature, which suppresses the phonon-scattering rates at high temperature, giving rise to a milder T −0.86 decay of the calculated κ l that is closer to the temperature dependence of T −0.88 and T −0.82 found in the experimentally measured κ l of TlInTe 2 in two sets of experiments [38,39].
Since the presence of strong anharmonicity is an important characteristic of many ultralow-κ l compounds, it is necessary to asses its strength in TlInTe 2 . Anharmonicity of the phonon modes is estimated by the mode Gruniesen parameters (γ λ = − dlnω λ dlnV ) that quantify the change of the phonon frequencies (ω λ ) with respect to the change in the unit cell volume (V ), where λ is a composite index that combines the wave vector (q) and phonon branch index (s). While for weakly anharmonic solids γ λ 's are close to 1, for highly anharmonic materials γ λ 's become much larger than 1. Some examples of compounds that possess large γ λ 's (hence, strong anharmonicity) and ultralow-κ l are AgBiSe 2 [16] and SnSe [15,50]. Previous studies [39][40][41]51] have shown that the phonon modes of this ABX 2 family of materials exhibit large γ λ 's. Thus, in the presence of strong anharmonicity, the harmonic description of the phonon frequencies of the compounds in this family becomes inappropriate at finite temperatures as the anharmonicity induces multiphonon interactions, giving rise to shifts in the phonon frequencies and broadening of the phonon states. As a consequence, the phonon frequencies of these anharmonic solids are expected to have strong renormalization effects at finite temperatures, which can crucially alter both the magnitude and T-dependence of the κ l that deviates from the ideal T −1 behavior found in weakly anharmonic solids. Indeed, the weak temperature dependence of the experimentally measured [38,39] κ l of TlInTe 2 , which decays as T −0.88 and T −0.82 , signifies the prevalence of severe anharmonicity in the crystal structure of TlInTe 2 .

Temperature-induced anharmonic phonon renormalization
We use the self-consistent phonon theory (SCPH) [29,31,52,53] to renormalize the phonon frequencies of TlInTe 2 at finite temperatures including anharmonic effects [28] which are treated as the phonon self-energies [31,54]. Within the SCPH theory, the anharmonically renormalized phonon frequency is determined from the pole of the many-body Green's function. Considering only the first-order contribution to the phonon self-energy arising from the quartic anharmonicity, the SCPH equation [29] is written as Ω 2 where ω λ is the harmonic phonon frequency at T = 0 K and Ω λ is the anharmonically renormalized phonon frequency at finite T. The quantity I λ is defined as: , n, and V (4) (λ, −λ, λ 1 , −λ 1 ) are the number of sampled wave vectors, the reduced Planck constant, the temperature-dependent phonon population, and the reciprocal representation of the fourth-order IFCs, respec- 3. Heat transport due to particle-like propagation of phonons. (a) Particle-like contributions to the lattice thermal conductivity (κ l ) of TlInTe2 calculated utilizing the harmonic (i.e., HA) and the anharmonically renormalized (i.e., SCPH) phonon frequencies using the PBTE. κ l and κ ⊥ l are the components of κ l , which are parallel and perpendicular to the chain directions in the crystal structure of TlInTe2, respectively. (b) The averages of the κ l and κ ⊥ l components are compared against the experimentally measured values which are denoted with Exp1 and Exp2 taken from Ref. [38] and Ref. [39], respectively. (c) Three-phonon scattering rates, and (f) phonon group velocities rates obtained using the harmonic (i.e., HA) and anharmonically renormalized (i.e., SCPH) phonon frequencies at T = 300 K. Mode contributions to κ l and κ ⊥ l are analyzed through (d) the cumulative plots and (e) their first-order derivatives with respect to the anharmonically renormalized phonon frequency at T = 300 K. The solid black line in (c) assumes the scattering rates of the phonon modes to be twice their frequencies according to the Cahill-Watson-Pohl model [24].
tively. The temperature dependence of the SCPH equation is contained in the phonon population term that obeys the Bose-Einstein statistics. Since Ω λ and I λ are inter-dependent on each other, the SCPH equation is solved self-consistently until the convergence in Ω λ is achieved.
We present the anharmonically renormalized phonon dispersions of TlInTe 2 in Fig. 2a in the temperature range between 0-700 K. The harmonic phonon dispersion (i.e., T= 0 K) of TlInTe 2 exhibits two groups of low-frequency optical phonon branches with very small dispersions, that are characteristic of the rattler atoms in the crystal structure. The lowest-energy rattling phonon branch (denoted as R LO ) arises from the longitudinal vibration of Tl 1+ ions along the direction of the InTe 4 chains (Fig. 1c) inside the hollow Thompson cube. The second group of rattling phonon branches (R T O ) appears above R LO , where the Tl 1+ cations vibrate along the transverse direction (i.e., perpendicular to the InTe 4 chains, Fig. 1d). Both R LO and R T O phonons are highly localized at T = 0 K as revealed by analysis based on the phonon participation ratio [42,55], which does not change as they are anharmonically renormalized at finite temperatures (see see Supplementary Figure 6 and Supplementary Note 4). The atom-resolved phonon density of states at T=0 and T = 300 K are shown in Fig.  2b, which clearly show the two peaks associated with R LO and R T O phonons, which gradually convolute as the temperature increases. It is seen from Fig. 2a that only the phonons up to 100 cm −1 are strongly hardened while the phonons above 100 cm −1 show very weak hardening or softening. The most drastic temperatureinduced changes are observed in the frequency hardening of R LO and R T O phonons (Fig. 2c). It is interesting to note that the frequency of R LO phonon is smaller than that of R T O due to the chain-like crystal structure and weak bonding of Tl atoms which vibrate slowly but with much larger amplitudes along the c-direction in the empty space within the lattice.
Particle-like contributions to κ l We start our analysis of the thermal conductivity by examining the calculated κ l obtained from solving the PBTE using the harmonic (i.e., T= 0 K) and anharmonically renormalized phonon frequencies (at finite T) obtained using the SCPH method. Hereafter, we denote these results by κ l (HA), and κ l (SCPH), respectively, as shown in Fig. 3a. It is seen that κ l (HA) changes significantly when the renormalized phonon frequencies are incorporated. For example, κ ⊥ l (HA) and κ l (HA) increase by 55 % and 38 %, respectively, at 300 K when the renormalized phonon frequencies are used to solve the PBTE. Here, and ⊥ symbols indicate components of κ l parallel and perpendicular to the chain direction (i.e., the c-axis) in the crystal structure of TlInTe 2 , respectively. We compare the average κ l calculated under both HA and SCPH methods with that of available experimental measurements [38,39] of κ l in TlInTe 2 . From Fig.  3b, we see that the average κ l (HA) is significantly underestimated in magnitude compared to the two sets of available experiments [38,39]. While κ l (HA) decays as T −1 according to well-known behavior of the lattice thermal conductivity found in weakly anharmonic solids, the two experimentally measured κ l data decay as T −0.82 and T −0.88 . The deviation from T −1 behavior and the weaker temperature dependence signify the presence of a strong higher-order anharmonicity in TlInTe 2 , and thus necessitate the anharmonic renormalization of the harmonic phonon frequencies. The calculated average of κ l within the SCPH method increases in magnitude with respect to κ l (HA), but still does not agree well with either sets of experimentally measured values of κ l . However, the effect of anharmonic renormalization significantly improves the temperature dependence, which varies as T −0.86 (Fig.  3b), bringing it closer to the experimental observations. Next, we analyze how the anharmonic renormalization affects the two key ingredients that enter into the PBTE (see Supplementary Note 2), namely the three-phonon scattering rates (Fig. 3c) and the phonon group velocities (Fig. 3f). According to the phonon gas model, where the phonons behave like particles, the maximum phonon scattering rate of a phonon mode is assumed to be twice its frequency [24], which is denoted with a black line in Fig. 3c. The HA scattering rates are strongly suppressed by the anharmonic effects in the frequency region below 100 cm −1 which enhances the phonon lifetimes, leading to a large increase in κ l (SCPH) compared to κ l (HA). The SCPH scattering rates are well below the solid black line, indicating the dominant particle-like nature of the phonons in TlInTe 2 that rules out the existence of a hopping channel of instantaneously localized vibrations as shown in Ref. [51]. The phonon group velocities are weakly hardened below 100 cm −1 (Fig. 3f) due to the anharmonic effects at finite temperatures and its effect on κ l (SCPH) is less significant than that of the scattering rates. To examine the phonon mode specific contributions to κ l , we show the cumulative plot of κ l (SCPH), and its derivative κ l (SCPH) with respect to renormalized phonon frequency in Fig. 3d and Fig. 3e, respectively. The cumulative plot for κ ⊥ l changes rapidly up to 60 cm −1 and reaches a plateau above it, showing that only the acoustic and low-energy optical phonons mainly contribute to it. This is also clear from Fig. 3e, where large peaks are present mainly up to 60 cm −1 . On the other hand, examining κ l in both Fig. 3e and Fig. 3f, reveals that κ l has large contributions coming not only from the acoustic and low-energy optical phonons (< 80 cm −1 ) but also from the high-energy optical phonons (∼ 140-180 cm −1 ).

Off-diagonal contributions to κ l
To understand the difference between the measured κ l and calculated κ l (within the SCPH method) of TlInTe 2 , we recognize that κ l obtained after solving the PBTE i.e., κ l (PBTE) only accounts for the diagonal components of the heat-flux operator [32,47,48]. We calculate the offdiagonal (OD) contributions i.e., κ l (OD) using the renormalized phonon frequencies and obtain the total [32] κ l as: κ l (tot) = κ l (PBTE) + κ l (OD) (see Supplementary Note 1). It was shown that while κ l (OD) is negligible in compounds like Si and diamond, while it is very important in CsPbBr 3 and tetrahedrites [32,35,36,46]. One key quantity that enters into the expressions for κ l (PBTE) and κ l (OD) is the generalized phonon group velocity operator (see its expression in SI), V s1,s2 , where s 1 and s 2 are phonon branch indices. While κ l (PBTE) utilizes only the diagonal components (i.e., s 1 = s 2 ) of V s1,s2 , the κ l (OD) term uses the off-diagonal (i.e., s 1 = s 2 ) components of V s1,s2 . Fig. 4a shows the κ ⊥ l and κ l components of κ l (OD), and their average as a function of temperature. Although on an absolute scale these values are small, for low-κ l compounds, these values are quite significant. For example at 300K, κ ⊥ l (OD) and κ l (OD) account for 14 % and 10 % contributions to that of κ ⊥ l (SCPH) and κ l (SCPH) values, respectively. However, with increasing temperature, the relative OD contributions also increase. For example, the contributions of κ ⊥ l (OD) and κ l (OD) increase to 26 % and 15 %, respectively at 700 K. When the average κ l (OD) is added on top of the average κ l (SCPH) terms, the resulting κ l agrees very well both in magnitude with one set of experiments, Exp1 [38] and in temperature dependence that decays as T −0.80 (Fig. 4b). The off-diagonal (OD) contributions to κ l (i.e., κ l , κ ⊥ l , and their average) calculated using the anharmonically renormalized phonon frequencies at each temperature. (b) Total κ l (calculated within SCPH method) including the OD contributions reproduces the Exp1 [38] well. Additional grain-boundary scatterings on top of SCPH+OD contributions to κ l reproduce Exp2 [39]. Grain sizes ranging from 40 to 1000 nm are considered here. Three-dimensional visualizations of the mode specific contributions to κ l (OD) as a function of anharmonically renormalized phonon frequencies (i.e., Ωs) at T = 300 K for the (c) κ ⊥ l and (d) κ l components. These quantities are plotted per cross-sectional area of the Brillouin zone to clearly highlight the phonon frequencies that primarily contribute to the κ l (OD). A color scale corresponding to the variation of colors in Figs. 4c and 4d is shown, which has the same unit as indicated in the z-axes of both figures.

Effects of grain boundaries
Having reproduced the first set of experiments (i.e., Exp1 [38]), we notice that the measured κ l of TlInTe 2 in an another experiment (i.e., Exp2 [39]) is lower than the former. Since the measurements were performed in the polycrystalline samples of the compound in which the presence of grains boundaries are generally common, the lower κ l in Exp2 [39] most likely originates from the additional scatterings of the heat carrying phonons due to the grain boundaries. We introduced grain-boundary scatterings on top of the three-phonon and isotope scatterings and calculated κ l to see if this additional scattering mechanism can explain the second experiment (i.e., Exp2 [39]). Assuming that the boundary scatterings are predominantly diffuse, the grain-boundary scattering rates are given by τ −1 gb,λ = 2|v λ | L , where v λ and L are the phonon group velocity and the averaged grain size, respectively.
The total scatterings rates of the phonons are obtained by applying the Matthiessen's rule, which are then used to calculate κ l (see Supplementary Note 2 for details). We show in Fig. 4b the effect of grain boundary scatterings on κ l as a function of grain size. It can be seen that κ l is significantly modified in magnitudes particularly in low temperature (< 550 K). Our calculations show that while the phonon-scatterings due to the grains with an average size of 40 nm reproduce Exp2 [39] quite well below T= 550 K, larger grain sizes give rise to better agreement with Exp2 above 550 K (Fig. 4b).

DISCUSSION
After successfully explaining the experimental results, we now closely examine the mode-specific contributions to the κ ⊥ l (OD) and κ l (OD), which are shown in Fig. 4c and Fig. 4d, respectively, as three-dimensional plots av-eraged over the planar area in the Brillouin zone. We see that phonons with very similar frequencies with s 1 = s 2 (i.e., near the diagonal in Fig. 4c) below 100 cm −1 primarily contribute to κ ⊥ l (OD). On the other hand, quasidegenerate phonons in the frequency ranges of 20-100 cm −1 and near 180 cm −1 mainly contribute to κ l (OD). It is interesting to note that the anisotropy (i.e., κ l /κ ⊥ l ) of the PBTE calculated contributions (κ l (SCPH)) is much stronger than the anisotropy present in the κ l (OD) contributions. For example, κ l /κ ⊥ l = 2.6 for the PBTE results while for the OD terms, κ l /κ ⊥ l =1.8 at T = 300 K. While the anisotropy in the phonon dispersion generally results in the anisotropy in κ l obtained from PBTE as the group velocity in the PBTE is related to the single phonon mode, the much weaker anisotropy in κ l (OD) stems from the fact that the OD heat transport takes place between coupled phonon eigenstates where the group velocity is associated with two different (quasi degenerate) frequencies. Thus, κ l (OD) has a weak dependence on the slopes of the phonon branches, which partially counteracts the anisotropy in κ l (PBTE). It is interesting to see that while κ ⊥ l (OD) increases with temperature, κ l (OD) decreases with temperature. Our analysis shows that this opposite temperature-dependence arises from the renormalization of the phonon frequencies at finite temperatures. We show in the Supplementary Figure  7 that when κ ⊥ l (OD) and κ l (OD) are calculated using the harmonic (i.e., T= 0 K) phonon frequencies, both the components increase with temperature. However, we note that in all cases, the temperature-dependence of both components of κ l (OD) are quite weak.
A previous attempt to understand the lattice thermal conductivity of TlInTe 2 was performed by Wu et al. [51] where the two-channel model [56] of thermal conductivity was used to explain the experimentally measured κ l . The two-channel model has been invoked in cases where the particle-like description of phonons becomes insufficient to describe the observed κ l of materials. Thus, a second channel of lattice heat conduction is introduced to compensate for the κ l , which is attributed to arise from the localized hopping among the uncorrelated phonon oscillators [56]. In the two-channel model used by Wu et al. [51], the phonon frequencies were treated within the harmonic approximation (T= 0 K) with no effects of temperature, which were then used to solve the particlelike contribution to κ l using the PBTE. Hence, the calculated [51] κ l underestimated the measured κ l values of TlInTe 2 [38,39]. To account for this difference, the contribution to κ l coming from the second channel was calculated using the Cahill-Watson-Pohl (CWP) model [24] which estimates the minimum κ l achievable in disordered solids. We note that our calculated κ l within the unified theory that utilizes the anharmonically renormalized phonon frequencies successfully reproduces the experiments without requiring to invoke the second channel of the two-channel model [56] of κ l . Also, our calculated total κ l shows better agreement with Exp1 [38] than the κ l calculated using the two-channel model by Wu et al. [51] (see Supplementary Figure 4).
We note that we did not consider the effect the lattice thermal expansion in our calculations. In general, the thermal expansion will soften the phonon frequencies which give rise to the reduced phonon lifetimes and hence a lower-κ l at high temperatures. However, this reduction is partially counterbalanced by the increase in the κ l due to the anharmonic heat flux [46] at high temperature, which is also absent in the current formalism. We did not calculate the thermal expansion coefficient and the anharmonic heat flux of TlInTe 2 as their calculations are computationally very expensive for non-cubic lattice. Nonetheless, the interplay of the thermal expansion and anhrmonic heat flux in TlInTe 2 and any ultralow-κ l materials in general is worth exploring in future studies.
In summary, we have investigated the microscopic origin and the underlying physical principles governing the extremely low and weakly temperature-dependent lattice thermal conductivity in the chain-like crystalline semiconductor TlInTe 2 using a unified theory of lattice heat transport that combines both the particle-like propagation and wave-like tunneling of phonons. To treat the strong anharmonicity present in TlInTe 2 , we have applied the SCPH theory to anharmonically renormalize the phonon frequencies at finite temperatures considering the quartic anharmonicity. Our calculated κ l using the PBTE coupled with the SCPH method (i.e., κ l (SCPH)) and the off-diagonal contributions (i.e., κ l (OD)) arising from the wave-like tunneling of phonons successfully reproduce both the magnitude and temperature dependence (milder than T −1 decay) of the measured κ l in one set of experiments. Adding of additional grain-boundary scatterings on top of κ l (SCPH)+ κ l (OD), our calculated κ l reproduces well the second set of experiments. Our work thus highlights the important roles of (i) temperature induced renormalization of the harmonic phonon frequencies, particularly the low-energy rattling and other optical phonons, by the anharmonic effects, and (ii) the OD contributions to κ l in the heat-flux operator to correctly explain the origin of unusual lattice thermal transport of strongly anharmonic solids. The detailed microscopic understanding of the heat transfer mechanism obtained in this work will provide guidance towards the rational design and discovery of hitherto unknown ultralowκ l compounds for various energy applications.

Density functional theory calculations
We perform all density functional theory (DFT) calculations using the Vienna Ab-initio Simulation Pack-age (VASP) [57,58] employing the projector-augmented wave (PAW) [59] potentials of Tl (5d 10 6s 2 6p 1 ), In (4d 10 5s 2 5p 1 ) and Te (5s 2 5p 4 ) and utilized the PBEsol [60] parametrization of the generalized gradient approximation (GGA) [61] to the exchange-correlation energy functional. We use a kinetic energy cut-off of 520 eV, and Γ-centered k-point mesh of 12 ×12×12 to sample the Brillouin zone. The fully optimized lattice constants (a=8.406Å, c=7.134Å) agree very well (absolute error < 1 %) with the experimentally reported values (a=8.478 A, c=7.185Å) [62]. We choose the high-symmetry k-path in the Brillouin zone while potting the phonon dispersion following the convention of Setyawan et al. [63]. Phonon dispersions are calculated using the second-order interatomic force constants (IFCs) using Phonopy [64]. Since the calculated phonon scattering rates strongly depends on the phonon frequencies which can be quite sensitive to the supercell size, we have performed convergence tests (see Supplementary Figure 1), and adopted 2×2×2 supercell for the calculations of the second-order IFCs.

Thermal conductivity calculations
To renormalize the phonon frequencies at finite temperature, and to calculate the lattice thermal conductivity (κ l ) using the PBTE, an accurate estimation of anharmonic IFCs, namely the third-order and fourth-order IFCs are required. We obtain these anharmonic IFCs using the compressive sensing lattice dynamics (CSLD) method [28,29,65,66] using 2×2×2 supercell and 6×6×6 Γ-centered k-point mesh. While constructing the third and fourth order IFCs, a cut-off radius (r c ) is chosen, above which all three-body and four-body atomic interactions are discarded, respectively. Although with the increasing order of the IFCs, the atomic interactions become very short ranged, the value of r c can be quite critical in correctly calculating κ l [15]. We have performed convergence tests of κ l against r c (see Supplementary  Figure 3) which are chosen to be 7.56Å and 4Å for the third and fourth-order IFCs, respectively. The r c value for the third-order IFCs is carefully examined to give good convergence in the calculated κ l . While we have calculated the particle-like contributions to κ l using the PBTE considering the three-phonon, isotope and grain-boundary scatterings (see Supplementary Note 2 for details), the off-diagonal contributions to κ l has been calculated by explicitly evaluating the off-diagonal terms of the heat-flux operator (see Supplementary Note 3 for details). Both these terms have been calculated by utilizing the anharmonically renormalized phonon frequencies obtained at finite temperatures.

Self-consistent phonon calculations
We solved the SCPH equations self-consistently until the phonon frequencies are converged to a given small threshold (e.g., 10 −3 cm −1 ). using a relatively dense 6× 6×6 mesh of q-points. The renormalized phonon frequencies and eigenvectors are then used to obtain the renormalized harmonic IFCs through the inverse Fourier transformation. These renormalized IFCs are utilized to calculate phonon dispersions at finite temperatures. We note that in this study we did not take into account the second-order correction to ω λ due to the cubic anharmonic term because their contributions have been found less significant than the phonon frequency hardening by the quartic anharmonicity in some of the low-κ l systems such as PbTe [67] and clathrates [30]. However, we note that there are complex compounds like the tetrahedrites where the role of the cubic anharmonicity is found to be quite significant [35]. We calculate the particle-like contributions, κ l (PBTE), by iteratively solving the PBTE using the ShengBTE code [44]. We use 12× 12×12 qpoint mesh to obtain κ l with good convergence (see Figs. S3c-d). We present two components of the κ l in the results section: κ l and κ ⊥ l , which are parallel and perpendicular to the InTe 4 chain direction (i.e., along the c-axis in Fig. 1a) in the crystal structure of TlInTe 2 , respectively. While κ l is the zz-component of the κ l -tensor, κ ⊥ l is the directional average of the xx and yy-components of the κ l -tensor. The average κ l is defined as the directional average of xx, yy and zz components of the κ l -tensor throughout the manuscript. The T −α fitting of the experimentally measured and calculated κ l data of and (c) Quest high-performance computing facility at Northwestern University which is jointly supported by the Office of the Provost, the Office for Research, and Northwestern University Information Technology.