Kondo-like phonon scattering in thermoelectric clathrates

Crystalline solids are generally known as excellent heat conductors, amorphous materials or glasses as thermal insulators. It has thus come as a surprise that certain crystal structures defy this paradigm. A prominent example are type-I clathrates and other materials with guest-host structures. They sustain low-energy Einstein-like modes in their phonon spectra, but are also prone to various types of disorder and phonon-electron scattering and thus the mechanism responsible for their ultralow thermal conductivities has remained elusive. Our thermodynamic and transport measurements on various clathrate single crystal series and their comparison with ab initio simulations reveal an all phononic Kondo effect as origin. This insight devises design strategies to further suppress the thermal conductivity of clathrates and other related materials classes, with relevance for thermoelectric waste heat recovery and, more generally, phononic applications. It may also trigger theoretical work on strong correlation effects in phonon systems.

T he increasing world wide energy consumption and the associated climate change call for enhancing the overall energy efficiency of technological energy conversion processes. Thermoelectrics [1][2][3][4][5][6] can convert waste heat into electricity and could thus contribute to such an efficiency increase. The perfect thermoelectric material combines low phonon thermal conductivity with high electrical conductivity and high thermopower. Materials classes with cage-like structures such as the clathrates 7,8 or skutterudites 4,9 show surprisingly low phonon thermal conductivities κ ph , even in perfect single crystals. Type-I clathrates G 8 H 46 , that show promising thermoelectric figures of merit 5,10,11 , consist of host atoms (H) that encapsulate guest atoms (G) in a framework of oversized cages (Fig. 1a). Recent inelastic neutron and X-ray scattering studies uncovered energetically low-lying optical phonon modes which were attributed to a so-called rattling motion of the guest atoms in the cages and are thus referred to as Einstein modes [12][13][14] . Similar observations were also made for skutterudites 9 and clathrate hydrates 15 . A severe flattening of the acoustic phonon branches at energies near the optical modes was observed [12][13][14] , and attributed to a finite coupling between the guest atoms and the host cages. Early thermal conductivity studies on type-I clathrates 16,17 , inspired by investigations on glasses 18 , tried to capture this guest-host coupling by including terms for resonance scatteringoriginally proposed to describe the resonance interaction between phonons and non-paramagnetic defects 19 , and for scattering from tunneling states 18 into phenomenological treatments 20 . Whereas empirical multi-parameter fits including these terms can indeed model thermal conductivity data of several type-I clathrates below about 50 K (refs. 16,17 ), inconsistencies with other results have been pointed out more recently. Notably, such modeling can neither explain the exceptionally long phonon lifetimes 12,21 nor the large thermal conductivity differences between structurally very similar n-and p-type versions of these materials 22 . The recent achievement of ab initio calculations of the thermal conductivity of a type-I clathrate, based on intrinsic phonon-phonon Umklapp scattering processes 23 , represents a major step forward. It showed that, unlike in the resonance scattering picture, the phonon lifetimes are reduced in a wide frequency range. Nevertheless, also DFT calculations can explain experimental thermal conductivity data only at low temperatures, but strongly undershoot them at high temperatures.
Here we propose an alternative mechanism, a Kondo-like phonon-phonon interaction, that can explain the difference as the disentanglement of Kondo-coupled acoustic and rattling phonon modes above the phonon Kondo temperature. In general terms, the Kondo effect describes non-commutative scattering of an extended wave from a localized entity with internal degree of freedom. It has been studied in many different settings, including magnetic, orbital, charge, and local vibrational degrees of freedom (see Supplementary Note 8 for more details). However, to the best of our knowledge, it has never been studied in an all phononic context.

Results
Steps beyond the state-of-the-art. To demonstrate the all phononic Kondo effect, we had to go beyond the state-of-the-art in thermoelectric clathrates' research in several respects. Firstly, to overcome the bottleneck of large uncertainties and systematic errors of standard thermal conductivity measurements at elevated temperatures, we developed a high-precision implementation ( Fig. 2a and "Methods" section) of the 3ω technique 24 . A second key ingredient for our work are sample series of single crystals, which allow us to pin down the dominating instrinsic phonon scattering mechanisms by employing an empirical model 20 , modified to be in line with the most recent findings 21,23 on the phonon band structure and phonon dynamics in type-I clathrates. Thirdly, we compare both the thermal conductivity data and thermodynamic data to ab initio lattice dynamics calculations, thereby providing compelling evidence for the phonon Kondo effect.
Phonon thermal conductivities. We start by giving an overview of the phonon thermal conductivities κ ph of all type-I clathrate single crystals studied here (Fig. 2b). Although type-I clathrates are often referred to as model systems for the enigmatic phononglass electron-crystal concept 1 , the most prominent feature of κ ph (T) is a large crystal-like maximum. It is narrow and occurs at a surprisingly low temperature of only about 10 K. As will be further detailed below, this narrow low-temperature maximum is a strong indication for the phonon transport not being controlled by the natural energy scale k B Θ D of the Debye temperature Θ D but by a much smaller energy scale k B Θ E .
To visualize the effect of a reduced Debye temperature, we use a simple phenomenological lattice thermal conductivity model: A modified version of the standard Callaway model 20 , where we replaced the Debye temperature Θ D by a variable temperature Θ E (Supplementary Note 1). The left panel of Fig. 3b shows temperature dependences of normalized phonon thermal conductivities calculated within this model for a series of different characteristic energies k B Θ E , keeping fixed scattering rates τ À1 D for defect scattering, τ À1 B for boundary scattering, and τ À1 phÀel for phonon-electron scattering. Resonance scattering is of minor importance here ( Supplementary Fig. 4). With decreasing Θ E the broad maximum sharpens and is shifted to lower temperatures.   Fig. 3 Comparison of spin and phonon Kondo effect. a Schematic dispersion relations for electrons in heavy electron systems (left) and phonons in an analogously defined "heavy" phonon system (right). The blue and red curves represent the non-interacting dispersive and localized entities, respectively, the violet curves the hybridized interacting states. As function of temperature, the system evolves from non-interacting well above the Kondo temperature to interacting well below it. For simplicity, phonon branches above ω E are neglected. This is justified in real clathrates by the presence of multiple Einstein modes, resulting in multiple anticrossings 13,14 . The Debye model (blue line, right) assumes ω = v s q where v s is the sound velocity. The new dispersion relation (violet, right) is characterized by the group velocity v g = ∂ω/∂q. It equals the sound velocity only at low wave vectors and frequencies.
b Temperature-dependent phonon thermal conductivity, normalized to its maximum, calculated using a modified Callaway model (Supplementary Note 1) for various Einstein temperatures Θ E (left) and corresponding data for the Ge-based type-I clathrate BCGG1.0 (Supplementary Table 1) and elemental Ge, electron irradiated and annealed at 77 K to similar defect densities 64 (right). c Specific heat and thermal expansion phonon Kondo anomaly, obtained by subtracting the experimental from the theoretical specific heat and thermal expansion curves of Ba 8 Table 1) and single crystalline Ge, irradiated 25 to a similar defect concentration. It is clear that the clathrate is governed by a strongly reduced energy scale. Both materials have similar Debye temperatures (Supplementary Note 2), but the clathrate has, in addition, a low-lying Einstein temperature Θ E , that apparently controls the behavior.
Dominating role of Umklapp scattering. Next, we prove, with two clathrate single crystal series, that the strong Umklapp scattering due to the new low-energy scale indeed dominates the lattice thermal conductivity in a wide temperature range. In the first such series, Ba 8 Supplementary Table 1), an increasing Ga content x is accompanied by a decreasing content y of host atom vacancies □. This is supported by the increase of the lattice parameter a (Fig. 4a) and the Hall mobility μ H = R H /ρ = 1/(neρ) (combination of the two panels in Fig. 4c) with x, as well as by the charge neutrality of this substitution (both the electronic γ term of the specific heat and the charge carrier concentration n are essentially independent of x, see Fig. 4b, c bottom), as discussed in the Supplementary Note 3. Because the mass difference between Ga and Ge is very small, Ga represents a much weaker scattering potential 26 in Ba 8 Cu 4.8 Ge 41.2−x −y □ y Ga x than a vacancy. Thus, we expect a decrease of defect scattering with x. The low-temperature phonon thermal conductivity κ ph (T) of all samples in this series shows a maximum at low temperatures that is systematically enhanced at nearly constant temperature with increasing x (Fig. 4d left). This trend is indeed well reproduced by a decrease of the point-defect scattering rate τ À1 D with x, as seen by the good agreement of simulation and data for this case (full and dashed lines in Fig. 4d left). A much less satisfying agreement is observed if the boundary scattering rate τ À1 B is allowed to change with x instead (Fig. 4d  right). The phonon-electron scattering rate τ À1 phÀel can be assumed to be x independent because of the above discussed charge neutrality. Above about 50 K, however, the κ ph (T) data for different x merge ( Fig. 4d left), indicating that scattering from point defects has become negligible. This can be rationalized by comparing the (Debye) phonon wavelength λ = 2πv s / ω, where v s is the sound velocity (related to β, Fig. 4b), with the size of the scattering center. A vacancy with the distortion surrounding it was estimated to have a diameter d ≈ 5 Å (ref. 27 ). Strong point-defect scattering of Rayleigh type, with an ω 4 dependence, occurs only if λ is at least an order of magnitude larger than d, corresponding to phonon energies ħω of 2.5 meV (30 K) and below. At much larger energies (and temperatures), defect scattering should be weak and frequency independent 28 .
The second sample series is the prototypical clathrate Ba 8 Ga 16−x Ge 30+x that has been much investigated in the past. It is an ideal system to study the importance of phonon-electron  to the x = 0 data (red full line) and simulations (dashed lines) using the parameters of the x = 0 fit except for the defect scattering rate τ À1 D (left) and the boundary scattering rate τ À1 B (right). The respective scattering rate was decreased to the indicated percentage in direction of the arrow. The speed of sound of each sample was determined from β, the Einstein temperature Θ E from bell-shaped contributions to C p /T 3 versus logT. e Phonon thermal conductivity vs temperature for two different Ba 8 Ga 16−x Ge 30+x samples from the literature 30,31 , with lines corresponding to our fit to the data from ref. 30 and a simulation using essentially the same parameters except for a strongly enhanced phonon-electron scattering rate τ À1 phÀel for the data from ref. 31 scattering because it has a low and essentially constant amount of point defects (as seen from the above, shown in ref. 26 , and explained and demonstrated in the Supplementary Note 9 and Supplementary Fig. 5, respectively, Ga does not act as strong point defect in Ge clathrates) but a charge carrier concentration that varies strongly with x (ref. 29 ). Indeed, different Ba 8 Ga 16−x Ge 30 + x single crystals reported in the literature 30,31 show severely different κ ph (T) at low temperatures. Figure 4e replots two extreme cases. The red line is our lowtemperature fit (T < Θ E /2) with the modified Callaway model discussed in Supplementary Note 1 to the data of ref. 30 . It takes Umklapp scattering τ À1 , and phonon-electron scattering τ À1 phÀel into account. The latter is calculated from the reported materials properties 30 as discussed in Supplementary Note 1. Interestingly, the κ ph (T) data of a different Ba 8 Ga 16−x Ge 30+x single crystal 31 can be very well reproduced by strongly increasing τ À1 phÀel and by only slightly adjusting τ À1 U (dashed line in Fig. 4e). Above about 50 K, the differences in κ ph caused by a different τ À1 phÀel vanish. This observation is in line with the fact that phonon-electron scattering can only occur for phonons with a wave vector q smaller than twice the Fermi wave vector k F . For the crystal of ref. 30 we estimate this to hold below about 140 K (Supplementary Note 2). At higher temperatures phonon-electron scattering is unlikely to be relevant.
A suppression of large low-temperature differences in κ ph (T) at higher temperatures has also been seen in other type-I clathrate single crystal series 32,33 , but precise high-temperature data on these are not available to date.
Taking both our clathrate series together we have managed to rule out the influence of defect and phonon-electron scattering on κ ph (T) of various type-I clathrates above 50 K. In addition, boundary scattering cannot contribute significantly in single crystals at these temperatures. Thus, at elevated temperatures, κ ph is dominated by intrinsic phonon-phonon (Umklapp) scattering. This allows us to pin down its microscopic origin, as shown in what follows.
Universal scaling. Remarkably, a broad range of clathrates, including even a gas hydrate 15,[34][35][36] , shows a universal scaling of the room-temperature phonon thermal conductivity with the product of sound velocity and Einstein temperature of the lowestlying rattling mode(s), κ ph ∝ v s Θ E (Fig. 5). The simple kinetic gas relation κ ph ¼ c v v 2 s τ=3 predicts κ ph to depend on the square of the sound velocity. In the Debye model the sound velocity is proportional to the Debye temperature and thus κ ph / v 2 s / Θ 2 D is expected for simple Debye solids. The modified scaling Þ shows that the above discussed energy renormalization is universal in clathrates.
Analogy with spin Kondo effect. A similar energy renormalization is seen in heavy fermion metals, where the (spin) Kondo effect rescales the Fermi temperature T F to the (spin) Kondo temperature T K (refs. 37,38 ). Figure 3a illustrates this analogy with schematic dispersion relations for electrons and phonons. In a band picture for heavy electron systems (left) a broad conduction band (blue) hybridizes with a flat, essentially non-dispersing 4f band of (renormalized) energy ϵ 4f (red). The hybridized bands (violet) are extremely flat near ϵ 4f , corresponding to quasiparticles with strongly renormalized effective masses. In the phonon case the strongly dispersing acoustic phonon mode, approximated here by the linear dispersion of the Debye model, takes the role of the broad conduction band, and the flat Einstein-like rattling mode at ħω E = k B Θ E corresponds to the narrow 4f band. The resulting hybridized band is again extremely flat in large portions of the Brillouin zone, giving rise to "heavy" phonons with extremely low group velocity v g at finite wave vectors. Even though the resulting T = 0 dispersion relation may look similar to the one obtained from ab initio lattice dynamics simulations 23 , there is an important difference: Being a strong correlation phenomenon, the phonon Kondo effect has a characteristic temperature dependence. Well above the Kondo temperature, the interacting states (violet) go over to the noninteracting ones (red and blue), an effect referred to as crossover from infrared slavery to asymptotic freedom 38 , which is absent in the lattice dynamics simulations (see also "Discussion" section).
The fact that the observed universal scaling of κ ph still contains v s (to linear power) is attributed to the absence of a Fermi level in phonon systems. Whereas the electrical resistivity in heavy fermion metals is dominated by the heavy electrons at the Fermi wavevector k F , phonons of all wave vectors, including longwavelength ones propagating with v s , contribute to κ ph .
Finally, we compare specific heat, thermal expansion, and thermal conductivity data for the prototypical clathrate Ba 8 Ga 16 Ge 30 with ab initio lattice dynamics calculations (see Methods) and reveal key characteristics of the Kondo effect. For the specific heat, the agreement between experiment and theory is excellent at low and high temperatures, but a distinct deviation is seen in between (Supplementary Fig. 1b). This difference gives rise to an anomaly near Θ E (Fig. 3c). It releases an entropy of order R ln2 per rattling atom at the 6d site. This is the behavior expected for any Kondo effect involving a 2-fold degenerate localized entity: Its degeneracy is lifted by the Kondo interaction, giving rise to the above entropy release. The entropy reaches 0.65 ln2 at 40 K. This temperature is considered as good estimate of the Kondo temperature 39 , which is thus found to be comparable to Θ E . The experimental thermal expansion is sizeably smaller than the theoretical prediction below 150 K, but agrees well with it above this temperature (Supplementary Fig. 2b). The difference between the theoretical and experimental curve closely resembles the specific heat anomaly (Fig. 3c). These discrepancies in specific heat and thermal expansion translate into a corresponding discrepancy in the mode-averaged Grüneisen parameter: Ab initio calculations predict an upturn at low temperatures, an effect usually associated with enhanced phonon anharmonicity at low frequencies, that is however absent in the experimental curve ( Supplementary Fig. 2c). Finally, the experimental phonon thermal conductivity is well described by the ab initio calculations below 50 K, but severely overshoots them at higher temperatures ( Supplementary Fig. 3b). The difference, plotted as inverse to represent a thermal resistivity (Fig. 3d), increases with decreasing temperature, with a slow (about −ln T) dependence, the hallmark of incoherent Kondo scattering above T K in the original spin Kondo effect 37,38 .

Discussion
To understand these results, it is important to assess the limitations of the used ab initio calculations (see "Methods" section). The phonon density of states (DOS) and the specific heat are calculated in the harmonic approximation, the thermal expansion in the quasi-harmonic approximation, and the thermal conductivity from anharmonic interatomic force constants 23 . The latter thus takes phonon-phonon interactions to lowest order into account. However, temperature dependences resulting from strong phonon correlation effects, most notably the Kondo disentanglement of acoustic and rattling modes above the phonon Kondo temperature, cannot be captured by these simulations. This is why the comparison of experimental data and such calculations can be used to extract the characteristic temperature dependences due to the correlations (somewhat like non-f reference materials are used as background to reveal Kondo physics in heavy electron compoundsunfortunately, empty type-I clathrates without the rattling atoms, that would represent such background here, do not exist). Specifically, the temperaturedependent specific heat is overestimated by the calculations as, unlike in the Kondo effect, no entropy is released by Bose populating temperature-independent anticrossings (circles in Fig. 3c). The thermal expansion calculations do not contain the non-trivial temperature dependence of the anharmonicities due to the phonon Kondo interaction, which leaves a corresponding imprint in the thermal expansion difference (stars in Fig. 3c). Finally, the strong Umklapp scattering that is captured by the thermal conductivity calculations does not contain the weakening above T K characteristic of an asymptotically free theory (Fig. 3d). Further support, independent of any theoretical modeling, comes from recent inelastic neutron scattering experiments. They reveal unexpected changes of the optical dynamical structure factors as function of temperature 21 that could be an indication of such a disentanglement of interacting acoustic and rattling modes at high temperatures. All these observations provide strong evidence for a previously unknown correlation effectan all phononic Kondo effectas microscopic origin of the peculiar thermodynamic and thermal transport behavior of clathrates.
We do not want to conclude without proposing a possible microscopic realization of the phonon Kondo effect in clathrates. The rattling motion of the guest atoms at the crystallographic 6d site within a soft plane 12,31,40 in the tetracaidecahedra (Fig. 1b) is known to have the lowest frequency 23 . Within this plane, due to the four-fold symmetry of the potential well, rattling occurs preferentially along two perpendicular soft directions 22,41,42 , which we refer to as 1 and 2. The two corresponding rattling modes (e 1 and e 2 in Fig. 1b), which are degenerate in energy for symmetry reasons, are thus identified as the pseudospins of the Kondo model (they represent spin-up and spin-down states in the spin Kondo effect). They can be approximated as simple Einstein oscillators, with the first excited state corresponding to the Einstein temperature Θ E observed in inelastic scattering experiments [12][13][14] . The two polarizations of the transverse acoustic phonons represent the corresponding degree of freedom of the itinerant species (e T1 and e T2 in Fig. 1b, thus representing the conduction electrons in the spin Kondo effect). The hybridization of the local (rattling) and extended (acoustic) phonon modes has been observed in inelastic neutron and X-ray scattering experiments [12][13][14][15]34 . A spin flip process in the spin Kondo effect corresponds to a scattering process with polarization change in the phonon Kondo effect. It can be visualized as follows: Assume the guest atom rattles in mode e 1 and a propagating phonon of polarization e T2 "hops" onto the cage, creating a distortion of the cage that resembles the effect of mode e 2 . This additional cage distortion will facilitate a change of rattling direction from e 1 to e 2 , and an accompanying change in polarization of the outgoing acoustic wave from e T2 to e T1 . Such an interaction of two modes with "opposite" polarization is analogous to an antiferromagnetic exchange interaction in the spin Kondo effect. Further details as well as a suggestion for a mapping of these ingredients onto a Kondo-type Hamiltonian are given in the Supplementary Note 8. We hope that our discovery will be taken up by a broad community of theorists, to develop both ab initio and model-based approaches that can describe the unusual observed phonon properties.
Our discovery of Kondo physics in an all phononic system is not only of fundamental interest (see discussion of non-canonical Kondo physics in other settings in Supplementary Note 8), but also has practical implications. First, it gives a concise description of how to tailor rattling materials for thermoelectric applications at elevated temperatures, which are most relevant for waste heatrecovery applications: Materials with the lowest possible Einstein and thus phonon Kondo temperatures should be found. Microscopically, this might be achieved by changing the mass and size of the guest atoms, but also by structural disorder on the guest site 43,44 and/or tailored guest-host charge transfer 43 .
An equally striking consequence is that Cahill's definition of the minimum thermal conductivity κ min ph (ref. 45 ), that is being referred to so extensively, needs to be reconsidered. With the energy renormalization also κ min ph is drastically reduced (Supplementary Note 1). This implies that for the materials in Fig. 5 there is still room for further improvement as their phonon thermal conductivities all lie well above their respective new κ min ph value (dashed red line in Fig. 5). Further lowering the κ ph of a given material could, for instance, be realized by nanostructuring or the introduction of dense dislocation arrays 46 . Long-wavelength phonons that are only weakly affected by the phonon Kondo effect could, at high temperatures, be effectively scattered by nanostructures that are small compared to their mean free path. Indeed, a few clathrates and skutterudites appear to violate the original Cahill κ min ph limit, though this has not been explicitly recognized in these works (Supplementary Note 2). We expect our results to trigger more systematic efforts along these lines.
The occurrence of low-lying Einstein-like phonon modes that interact with acoustic phonons is not limited to the class of clathrates and related cage-like materials. For the recently discovered family of Cu 12−x M x Sb 4 S 13 (M = transition metal) tetrahedrites, optical phonon branches involving out-of-plane vibrations of the three-fold coordinated Cu ions were predicted by ab initio calculations and were suggested as origin of the low thermal conductivities 47 . Interestingly, κ ph of this material fits perfectly into our universal scaling plot (Fig. 5). Other candidate materials are, e.g., PbTe (ref. 48 ), Bi 2 Te 3 (ref. 49 ), BiCu(Se,Te)O (ref. 50 ), Cu 3 SbSe 3 (ref. 51 ), and rattling-induced superconductors such as the β-pyrochlore oxides 52 , dodecaborides 53 , or VAl 10 (ref. 54 ), and possibly even amorphous and glassy materials 55 . Detailed investigations, such as presented here, are needed to test whether also in these materials the phonon Kondo effect is at work. Phonon Kondo systems transfer heat largely via low-frequency phonons of long mean free paths. As such they may be seen as promising intrinsic thermocrystals, for applications such as heat waveguides or thermal diodes in the emerging field of phononics 56 .

Methods
Synthesis and structural characterization. As starting material for the single crystal growth of BCGGx, two cylindrical rods with the same nominal composition Ba 8 Cu 4.8 Ge 41.2−x−y □ y Ga x (BCGGx) were prepared for each sample with x = 0.0, 0.2, 0.5, and 1.0 from high-purity elements using a high-frequency induction furnace. One rod with 7 mm in diameter and 60 mm in length served as the feed rod, the other one with the same diameter and 20 mm in length as the seed for the crystal growth. The single crystal growth was performed in a 4-mirror furnace equipped with 1000 W halogen lamps. The pulling speed of the rod was 3-5 mm h −1 . Both rods rotated in opposite direction (speed:~8 revolutions min −1 ) to ensure efficient mixing of the liquid and a uniform temperature distribution in the molten zone. A pressure of 1.5 bar Ar was used during the crystal growth.
X-ray powder diffraction data on BCGGx were collected using a HUBER-Guinier image plate system (Cu K α 1 , 8°≤ 2θ ≤ 100°). The lattice parameters (Fig. 4a) were obtained from least squares fits to indexed 2θ values employing Ge (a Ge = 0.5657906 nm) as internal standard.
Thermal conductivity. Commonly used laser flash methods measure the thermal diffusivity and thus need to be combined with specific heat and density measurements to calculate the thermal conductivity, which typically reduces the accuracy of this technique. By contrast, the 3ω method is an ac technique for direct thermal conductivity measurements. During a 3ω experiment the sample is heated locally and thus, in contrast to steady-state heat-flow experiments, errors due to radiation at elevated temperatures are reduced to a negligibly low level 24 . Furthermore, this method is insensitive to geometrical errors. This is because the only geometrical parameter entering is the length of a line heater. As it is usually prepared by means of photo or electron-beam lithography and sputtering, its length is very well defined (see below). In fact, the error of our 3ω thermal conductivity data, which we estimate to be below 5%, is dominated by uncertainties in the heater resistance and its temperature dependence.
For our studies, the narrow metal line serving as both the heater and the thermometer had a width of 20 μm and a length of 1 mm, with an uncertainly of 1 μm. To avoid electrical contact between heater and sample a thin layer of SiO 2 was first deposited on the polished sample surfaces by chemical vapor deposition. Then, a 4 nm thick titanium sticking layer and the 64 nm thick gold film were sputtered in an Ardenne LS 320 S sputter system. The heater structures were made by standard optical lithography techniques using a Karl Suess MJB4 mask aligner.
The metal line was heated by an oscillating current at a circular frequency ω, which thus leads to a 2ω temperature oscillation of both the heater and the sample. Due to the linear temperature dependence of the metallic heater, the 2ω temperature oscillation translates into a 3ω voltage oscillation, which is detected using a lock-in amplifier (7265, Signal Recovery). Applying the 3ω method 24 to bulk geometry, the measured in-phase temperature oscillation of the heater/ thermometer line is expected to be linear in logarithmic frequency f as long as the thermal penetration depth is large compared to the heater half width b and at least five times smaller than the sample thickness t (see boundaries indicated in Fig. 2a). The thermal conductivity κ of the material can then be extracted from the slope of the in-phase temperature oscillation ΔT vs log f.
Prior to the thermal 3ω voltage detection, the first harmonic and all related higher harmonics are subtracted from the signal using a carefully gain and phase calibrated active filter 24 based on a technique that allows to adjust the magnitude and phase of a reference signal as a function of frequency. In this way the main error sources of a 3ω experiment, spurious 3ω signals arising from harmonic distortions of the involved amplifiers, can be largely eliminated. Using an ultraprecise 100 Ω Z-foil resistor (temperature coefficient ±0.05 ppm°C −1 , tolerance ±50 ppm; Vishay) as a sample we found the background signal to be negligibly small within the entire frequency range (0.5 Hz to 50 kHz) of the experiment. Using the same resistor, the accuracy of our voltage controlled AC constant current source was found to be ±0.1% within the whole frequency range and for all tested excitations (100 μA to 10 mA). By further measurements on different resistors, load dependences were ruled out.
3ω measurements were done in the temperature range between 80 and 330 K. With this setup, we reproduced the thermal conductivity data on single crystalline Ba 8 Ga 16 Ge 30 of Sales et al. 41 within the error bar, estimated to only 3% in that work. This precision, which is remarkable for a steady-state experiment, was reached with special radiation shields and specific sample geometries 41,57 .
Below 100 K the data were completed by additional steady-state heat-flow experiments. The phonon thermal conductivities of all investigated materials were calculated by subtracting an electronic contribution, determined using the Wiedemann-Franz law with a constant Lorenz number of L 0 = 2.44 · 10 −8 WΩK −2 and electrical resistivity data measured on the same samples.
Specific heat. The specific heat was measured with a relaxation-type method using the 4 He specific heat option of a Physical Property Measurement System (PPMS) from Quantum Design. The addenda was measured separately prior to each sample measurement.
To study the phonon contribution to C p , first the electronic contribution was determined. At low temperatures, the specific heat can be approximated by C p /T = γ + βT 2 , where the Sommerfeld coefficient γ represents the electron contribution and the β parameter quantifies the Debye-like phonon contribution. For BCGGx, the data below 3.5 K are very well described by such linear fits (not shown).
Rattling modes, originating from localized oscillations of the guest atoms, can be revealed by analyzing C p /T 3 vs logT. Within such a representation, rattling modes appear as bell-shaped contributions on top of a Debye-like phonon background (not shown).
Hall effect and electrical resistivity. The electrical resistivity and the Hall coefficient were determined by a standard 6-wire technique using the horizontal rotator option of a Physical Property Measurement System (PPMS) from Quantum Design. Temperature-dependent Hall effect measurements were performed in a magnetic field of 9 T. The Hall resistivity was confirmed to be linear in fields up to 9 T at all temperatures down to 2.5 K. The Hall coefficient was analyzed within a simple one-band model, R H = 1/(ne).
Thermal expansion. Measurements of the coefficient of thermal expansion α L (T) = l −1 dl/dT were carried out by using a high-resolution capacitive dilatometer 58 , which enables the detection of length changes Δl ≥ 10 −2 Å. Relative length changes were measured along a principle axis of cubic Ba 8 Ga 16 Ge 30 . α L (T) is obtained by numerical differentiation of the Δl(T)/l data with respect to temperature (Supplementary Note 7).
Ab initio calculations. Ab initio density functional theory (DFT) simulations were conducted using the Vienna Ab initio Simulation Package (VASP) 59 , applying the projector augmented wave method 60 and the generalized gradient approximation (GGA) as proposed by Perdew, Burke, and Ernzerhof (PBE) 61 .
In a first step, a fully ordered, cubic 54 atom unit cell of Ba 8 Ga 16 Ge 30 with Ba at the Wyckoff sites 2a and 6d, Ge at 6c and 24k, and Ga at 16i was investigated, using a 5 × 5 × 5 k-point mesh and a plane wave cutoff of 500 eV. For a fixed lattice constant of 10.74 Å, corresponding to the experimental value, the atomic positions were relaxed. To test the impact of a more realistic Ga/Ge distribution 22 , we have in addition performed selected calculations for a disordered structure (see Supplementary Notes 9 and 11, and Supplementary Fig. 5). After reaching convergence (residual forces of less than 10 −3 eV Å −1 ), symmetry non-equivalent displacements were introduced into the relaxed structure of both the ordered and the disordered (see Supplementary Table 2) unit cell. For displacements of 0.02 Å, the restoring forces were determined, again by using the VASP code with the above described settings. From the obtained forces, the dynamical matrix was extracted and the ALAMODE code 23,62 was used to determine the phonon DOS and the specific heat ( Supplementary Fig. 1, main parts and insets) in the harmonic approximation.
The thermal expansion coefficient for both the ordered and disordered structure was determined within the quasi-harmonic approximation. For this purpose, the ground state structures and energies were determined for a series of unit cells with lattice parameters, corresponding to both decreased and increased cell volumes. For each of these volumes the dynamical matrix was then obtained in the same way as described above for the experimental volume. By using the phonopy code 63 the free energy was determined at the given volume and as a function of temperature. The free energy as a function of volume for a given temperature was then fitted to the Vinet equation of state, again using the phonopy code. Finally, from the free energy minima at a given temperature, the corresponding equilibrium volume at this temperature can be extracted, which then allows to access the thermal expansion coefficient. The ab initio thermal expansion curves are slightly rescaled to match the high-temperature data ( Supplementary Fig. 2).
The mode-specific Grüneisen parameter was obtained for both the ordered and the disordered model structures (Supplementary Fig. 5e). Again, the lattice parameter of the equilibrium structure was fixed to the experimental value of 10.74 Å and the dynamical matrix was determined as discussed above. To evaluate the behavior of the vibrational frequencies with respect to a volume change, the dynamical matrix was also determined for increased and decreased lattice parameter (±0.5%). Relating the frequency and volume change of a phonon mode allows then to extract the mode-specific Grüneisen parameter. This procedure is also available within the phonopy code 63 .
The lattice thermal conductivity of the fully ordered structure was calculated from ab initio results for the anharmonic force constants and by applying the relaxation time approximation, as implemented in the ALAMODE code 23,62 (see Supplementary Notes 10 and 11 for details). Our calculations reproduce the temperature dependence of the lattice thermal conductivity obtained by Tadano et al. 23 for the same structure model accurately, and the absolute value of the lattice thermal conductivity within a factor of 1.1 (which is most likely due to slightly different lattice constant used in both calculations). Expanding the range of anharmonic interactions from nearest to next-nearest neighbor shells leaves the temperature dependence largely unaffected (for Fig. 3d these results are used). The same is expected if the anharmonic force constants are allowed to vary (within ranges compatible with experiments) with temperature (Supplementary Note 11). Finally, the temperature dependence of the calculated thermal conductivity is shown to be robust against disorder within the Ga-Ge framework (for details see Supplementary Notes 9 and 11 and Supplementary Fig. 5f).

Data availability
The data sets generated and/or analyzed during the current study are available from the corresponding author on reasonable request.