Impact of temperature and mode polarization on the acoustic phonon range in complex crystalline phases: A case study on intermetallic clathrates

The low and weakly temperature-varying lattice thermal conductivity, κL(T), in crystals with a complex unit cell such as type-I clathrates is assumed to originate from a reduced momentum and energy space available for propagative lattice vibrations, which is caused by the occurrence of low-energy optical phonon modes. In the context of ab initio self-consistent phonon (SCP) theory, it has been shown that the cubic and quartic anharmonic interactions result in a temperature-induced energy renormalization of these low-lying optical branches which contributes to the anomalous behavior of κL(T) in structurally ordered type-I clathrates [T. Tadano and S. Tsuneyuki, Phys. Rev. Lett. 120, 105901 (2018)]. By means of inelastic neutron scattering, we provide evidence for this energy renormalization in temperature, which has been resolved for transversely and longitudinally polarized phonons in the single crystal type-I clathrate Ba7.81Ge40.67Au5.33. By mapping the neutron intensity in the momentum space, we demonstrate the coherent character of the low-lying optical phonons. The overall phonon spectrum and dynamical structure factors are satisfactorily reproduced by ab initio harmonic calculations using density functional theory with the meta-GGA SCAN functional and a fully ordered structure. However, a polarization-dependent cutoff energy with opposing temperature shifts for longitudinal and transverse acoustic dispersions is experimentally observed which is not reproduced by the simulations. Anharmonicity affects the energies of the low-lying optical phonons in the transverse polarization, which compares quantitatively well with available results from SCP theory, whereas differences are observed for the longitudinal polarization.

(Received 13 July 2020; accepted 25 November 2020; published 8 January 2021) The low and weakly temperature-varying lattice thermal conductivity, κ L (T), in crystals with a complex unit cell such as type-I clathrates is assumed to originate from a reduced momentum and energy space available for propagative lattice vibrations, which is caused by the occurrence of low-energy optical phonon modes. In the context of ab initio self-consistent phonon (SCP) theory, it has been shown that the cubic and quartic anharmonic interactions result in a temperature-induced energy renormalization of these low-lying optical branches which contributes to the anomalous behavior of κ L (T) in structurally ordered type-I clathrates [T. Tadano and S. Tsuneyuki, Phys. Rev. Lett. 120, 105901 (2018)]. By means of inelastic neutron scattering, we provide evidence for this energy renormalization in temperature, which has been resolved for transversely and longitudinally polarized phonons in the single crystal type-I clathrate Ba 7.81 Ge 40.67 Au 5. 33 . By mapping the neutron intensity in the momentum space, we demonstrate the coherent character of the low-lying optical phonons. The overall phonon spectrum and dynamical structure factors are satisfactorily reproduced by ab initio harmonic calculations using density functional theory with the meta-GGA SCAN functional and a fully ordered structure. However, a polarization-dependent cutoff energy with opposing temperature shifts for longitudinal and transverse acoustic dispersions is experimentally observed which is not reproduced by the simulations. Anharmonicity affects the energies of the low-lying optical phonons in the transverse polarization, which compares quantitatively well with available results from SCP theory, whereas differences are observed for the longitudinal polarization. DOI: 10.1103/PhysRevResearch.3.013021

I. INTRODUCTION
Tailoring the lattice thermal conductivity, κ L , of energyefficient semiconductors is a common materials issue in many applications such as for thermoelectric [1] and photovoltaic [2,3] conversion, phase change memories [4], and battery electrodes [5,6]. In the search for low κ L , the main strategy is the use of "complexity" at multiple length scales, from structural complexity within the crystal unit cell, to disorder, short-range order, and nanostructuring [7][8][9]. Crystals with * To whom correspondence should be addressed: stephane.pailhes@univ-lyon1.fr a high structural complexity and chemical bonding inhomogeneity [10], such as tetrahedrites [11] or type-I clathrates [12], often have a very low and almost temperature independent κ L of ∼0.5-2 Wm −1 K −1 in the 50-500 K range. The current understanding is that the heat conduction is mostly conveyed by well-defined acoustic phonons, which exist only in a limited range of the energy and momentum phase space, delimited by a continuum of nondispersive optical phonon bands [13][14][15][16][17][18]. The onset of this continuum at low energy, labeled E 1 , defines the upper energy limit of the acoustic regime such that it has been associated with a phononic low-pass acoustic filter [15] or a modified Debye energy [16,19]. E 1 can be changed by varying the chemical composition [13] or the structural topology [20]. Using the Boltzmann transport equation for phonons, the acoustic contribution of κ ac L is given by a cumulative spectral integral [13,15,19]: where κ ac,ν L (ω) is the mode thermal conductivity, ρ ν (ω) is the density of states per mode (DOS), ν is the longitudinal/transverse polarization index and ω q stands for the phonon relation dispersions. In addition, any variation in E 1 changes the whole phonon-phonon scattering phase space, thus impacting the acoustic phonon lifetimes entering into κ ν L (ω) [17,21]. Explaining the combined relationship between the complex crystal structure and related defects, the acoustic phonon properties, and the nature of E 1 , remains a fundamental challenge.
Type-I clathrates contain 46 framework atoms of mostly group 14 elements which arrange in a 3D covalent host network of face-sharing polyhedral cages that encapsulate alkali or alkali earth guest cations [12]. The structure is usually described using the cubic space group Pm3n (group 223) with a lattice parameter of about 1 nm. The cutoff energy E 1 is defined as the center of a distribution of optical phonon modes related to the dynamics of guest atoms located at the 6d-Wyckoff sites in tetrakaidecahedral (5 12 6 2 ) host cages. These modes lead to well-defined peaks in the phonon DOS [15,17,18,21] which results in a large deviation from the Debye-like T 3 temperature dependence of the lattice specific heat (C p ) at low temperature [25,26]. In literature, the temperature at the maximum of the T 3 -normalized heat capacity is commonly referred to as an Einstein temperature and corresponds well to E 1 in the phonon spectrum which we recall as the upper limit of the integral in Eq. (1). Recently, a phenomenological universal relation has been revealed in type-I clathrates between κ L and the product of the average sound velocity and E 1 [19]. The nature of the low-lying guest optical phonons with energies around E 1 is particularly intriguing. These modes are characterized by a very low phonon participation ratio ( 0.1), which is interpreted either as a signature of localization [27,28] or an effect of mode confinement [13,14,17]. In this latter case, the phonon is viewed as a Bloch state confined to a relatively small atomic pattern, the Ba(6d) atoms in this case, within the large complex unit cell whose periodic repetitions result in a special character. Moreover, the flatness of their dispersion and the concomitant high DOS provides a large momentum-and energy-conserving phase space for three-phonon scattering processes involving the acoustic modes [17,21]. The Ba 7.81 Ge 40.67 Au 5.33 structure [29][30][31], which will be dealt with in this paper, is shown in Fig. 1(a). Besides the structural complexity, the dative Au-Ba bonding of the Ausubstitution in Ba 7.81 Ge 40.67 Au 5.33 results in an off-centering of Ba atoms in the tetrakaidecahedral cages [13,32,33]. No correlation among the off-centering sites has been observed experimentally, indicating that they can be viewed as point defects. In the type-I structure Ba 8 Ge 40 Au 6 , molecular dynamics simulations performed at 300 K on a 2 × 2 × 2 supercell with independent random substitutions of Au atoms found no evidence for a correlated Au/Ge defect structure, thus no ordering among the off-centered Ba positions (see the Supplementary Material in Ref. [13]).
All ab initio phonon and κ L calculations reported for type-I clathrates have been carried out on the fully ordered model structure, with centered Ba atoms and full site occupation of all sites. The experimental phonon spectrum in type-I clathrates is qualitatively reproduced quite well by harmonic The cubic unit cell (space group Pm3n) contains two tetrakaidecahedral (5 12 6 2 ) and five dodecahedral (5 12 ) host cages formed by Ge atoms (light gray) with guest Ba (green) atoms encapsulated inside. One to three Au atoms (gold) substitute Ge atoms at the Wyckoff site 6c, which results in a slight distortion of the tetrakaidecahedron and an off-centering of the Ba atoms inside [13]. (b) The lattice thermal conductivity, κ L , for BGA (black circles) [13] is compared to different experimental measurements of κ L for Ba 8 Ge 30 Ga 16 (BGG), including those by Avila et al. [22], Sales et al. [23], May et al. [24], and the theoretical calculations of Tadano and Tsuneyuki [43]. Black and blue dashed lines show the deviation of κ L from 1/T for both BGA and BGG, respectively, at higher temperatures.
ab initio calculations [13][14][15]. Quantitatively, ab initio lattice dynamics studies of Ge clathrates showed decreased acoustic and low-lying optical mode energies in comparison to the experimental data (up to 40%), a discrepancy which has been recently overcome by the use of the meta-generalized-gradient approximation (meta-GGA) functional SCAN (strongly constrained and appropriately normed) for the exchange and correlation energy [34]. Ab initio simulations in perturbation theory, when limited to three-phonon processes, predicts a 1/T temperature dependence for acoustic phonon lifetimes and κ L (T ), which drastically fails to reproduce the experimental ∼T −0.25 in Ba 7.81 Ge 40.67 Au 5.33 [13] [see Fig. 1(b)], and also in other complex crystals like the ∼T −0.1 dependence in the quasicrystal approximant o-Al 13 Co 4 [35]. For the latter case, molecular dynamics simulations on an ordered model structure yield a T −0.5 dependence, whereas the inclusion of random disorder results in much closer agreement with experiment. This points not only to the importance of disorder, but also to either higher order anharmonicity or the effect of phonon energy renormalization including the polarization mixing of phonon eigenvectors, which are included in the molecular dynamics simulations but not in most ab initio based calculations.
A significant improvement to the ab initio approach has been achieved by self-consistent phonon (SCP) theory, which nonperturbatively treats the effects of anharmonicity [36][37][38][39][40][41][42]. An SCP study on an ordered model of Ba 8 Ge 30 Ga 16 reveals that the quartic anharmonicity leads to a softening of E 1 (T ) upon cooling, surpassing the usual hardening effect due to thermal expansion. Although this softening accounts for less than 10% of E 1 (T ) in the 0-300 K range, i.e., about 1 meV ∼ 12 K, it leads to a closer κ L (T ) matching [43]. This high sensitivity of κ L (T ) to E 1 (T ) results in the peculiar vibrational nature of the optical modes at E 1 (T ). It should be mentioned here that the Ba 8 Ge 30 Ga 16 SCP calculations assume that phonon polarization vectors are not affected by anharmonicity [43], which is a common approximation in complex crystals in order to limit the computational cost [38,41].
In this paper, we provide substantial experimental evidence of the importance of the anharmonic processes (resulting from cubic and quartic terms) for the type-I clathrate Ba 7.81 Ge 40.67 Au 5.33 by investigating the propagation direction, and the polarization and temperature dependencies of the cutoff energy E ν 1 (T ) by inelastic neutron scattering on a high-quality single crystal. By probing the mode symmetry in momentum space at different energies, we show that the distribution of the neutron intensity related to the low-lying optical vibrations is structured in the momentum space within a Brillouin zone, and from one Brillouin zone to another. This distribution of intensity in momentum as well as the overall phonon energies are satisfactorily reproduced by ab initio harmonic calculations using the (meta-GGA) functional SCAN done on a fully ordered structure. We confirm the agreement between the experimental results and simulations for phonon dispersions obtained along different high-symmetry directions and with both transverse and longitudinal polarizations.
However, some differences remain, especially in the region of the phase space where acoustic and optical phonons are hybridized. Experimentally, a polarization dependence of the acoustic-optical phonon coupling is observed such that the origin and the value of E 1 in the longitudinal and transverse acoustic (LA, TA) polarizations are different and exhibit opposite temperature dependencies, the effect of optical mode hardening with increasing temperature being found only for the lowest mode in the transverse polarization. This difference, which is not observed in our ab initio harmonic simulations using a fully ordered model, indicates either an effect of the particular defect cage structure caused by the Au substitutions or a more subtle anharmonic effect involving phonon polarization.
Furthermore, and in light of the recent SCP method calculations on Ba 8 Ge 30 Ga 16 [43], this experimental case allows us to quantitatively compare clathrate anharmonicity as found by both experimental and theoretical approaches. The rate of change in energy of E 1 (T ) in Ba 7.81 Ge 40.67 Au 5.33 was properly measured in the temperature range 100-550 K. After  [50]. The numerical derivative of a/a is the linear thermal expansion coefficient, α L . The volumetric thermal expansion coefficient, the subtraction of the thermal expansion contribution from this rate of change, the anharmonic contribution in the thermal shift of E 1 (T ) has been experimentally determined and directly compared to that of the quartic and cubic terms of the SCP method calculations for Ba 8 Ge 30 Ga 16 , providing experimental validation of the SCP method calculations for determining anharmonic effects in clathrates.

II. METHODS
Inelastic neutron scattering (INS) measurements were performed on the same high-quality single crystal of the type-I clathrate Ba 7.81 Ge 40.67 Au 5.33 as in our previous work [13], whose temperature dependence of the lattice thermal conductivity, κ L (T ), is shown in Fig. 1(b). The structural study, the chemical disorder caused by Au substitutions, and the thermal characterizations were reported in Ref. [13]. In addition, we have precisely measured its lattice thermal expansion, as depicted in Fig. 2, by means of neutron Larmor diffraction on the triple-axis spectrometer TRISP at the Heinz Maier-Leibnitz Zentrum (FRM-II, Germany) (see Appendix A 3). The INS intensity was recorded over a wide range of the momentum and energy phase space at 150, 300, and 530 K on the cold-neutron time-of-flight (TOF) spectrometer IN5 at the Institut Laue-Langevin (ILL, France). Details on the experimental settings, on the TOF instrumental resolution in momentum and in energy, and on the integration parameters used to produce the experimental phonon dispersions from fits of 1D-energy cuts at constant momentum (raw data and fits are shown in Appendix F), the high resolution Generalized Vibrational Density of States (GVDOS) in Fig. 3, and the mappings shown in Fig. 5, are given in Appendix A 1. The temperature dependence of the low-lying optical bands was further investigated on the cold-neutron triple-axis spectrometer IN12 at the ILL (see Appendix A 2). Preliminary experiments were also conducted on the triple-axis spectrometer 1T at the Laboratoire Léon Brillouin (LLB, France). In all INS experiments, the single crystal was mounted in a cryofurnace and aligned in the ([110]; [001]) scattering plane such that wave vectors of the form Q = 2π a (ζ , ζ , ξ ), with a = 10.7987(1) [13], were accessible.
Of particular importance in this work is the polarization term, which appears in the coherent one-phonon scattering function, S(Q, ω), which in turn is proportional to the double differential inelastic neutron cross section (see Appendix A). In the case of a coherent, one-phonon scattering process by a phonon of branch i, with energy ω q,i and polarization vector ξ i ω q,i , the neutron scattering function is written as [44] S ph (Q, ω) = n(ω) where n(ω) = 1 1−exp(−hω/k B T ) and comes from the detailed balance factor, and Q = q + G is the scattering vector given by the nearest reciprocal lattice vector G and the phonon wave vector q. Usually, phonons are measured in a Brillouin zone far from the -point (around Bragg peaks with high Miller indices) such that |G| |q| and Q ∼ G. In this work, we mainly discuss measurements performed around the Bragg peaks (006) and (222) whose moduli are much higher than π/a. The function F i D (Q) in Eq. (2) is called the dynamical structure factor (DSF) and is defined as [44] where b j , r j , M j , and W j (Q) are the coherent scattering length, fractional coordinates, mass, and Debye-Waller factor of the jth element, respectively. This expression closely relates to the nuclear structure factor which determines the Bragg peak intensity. It includes an additional term, the scalar product {Q · ξ j i (Q)} ∼ {G · ξ j i (Q)}, which contains the phonon polarization and can thus be used to distinguish longitudinal and transversal phonon modes by choosing the appropriate combination of phonon polarization wave vector and reciprocal lattice vector. The polarization vectors of longitudinal and transversal phonons are parallel and perpendicular to the phonon wave vector q, respectively.
The lattice dynamics were simulated using a fully ordered model of a type-I clathrate with the Ba 8 Ge 40 Au 6 composition in which the gold atoms occupy all the 6c host sites. The periodic density functional theory (DFT) code VASP [45][46][47] was used for structure optimization as well as for the determination of the harmonic force constants. While the projector augmented wave method was applied for describing the ionic cores, the meta-GGA functional SCAN was used to account for exchange and correlation [34,48]. The SCAN functional has been recently proven to reproduce the phonon spectrum in type-I clathrates with a much higher accuracy. The unit cell was relaxed to the ground state using a k-point mesh (5 × 5 × 5) centered at the zone center ( ) and a convergence criterion of residual forces of less than 10 −4 eV/Å using a plane wave energy cutoff of 500 eV. The lattice parameter of the optimized structures obtained by the SCAN functional is of 10.78 Å, very close to the experimental value of 10.7987(1) Å [13]. The Hellmann-Feynman forces were then calculated after introducing symmetrically non-equivalent displacements of ± 0.03 Å in the relaxed unit cell. These forces are given as an input to the Phonopy code [49] for the calculation of the dynamical matrix.

III. THE GRÜNEISEN PARAMETER
Before discussing anharmonicity at the phononic level, we first address the Grüneisen parameter, γ , which is a material constant that gives an idea about the amount of anharmonicity that exists in a material. This macroscopic property provides context to our experimental findings in this paper.
The mode specific Grüneisen parameter, γ i , for a phonon mode E i at molar volume V is defined as As a first approximation of the experimental Grüneisen parameter, however, we assume an averaged and temperature-dependent Grüneisen parameter for all modes, which, in the quasiharmonic approximation, depends on the volumetric thermal expansion coefficient α V , the Bulk modulus B, the molar volume V , and the specific heat at constant volume C V : The temperature dependence of C V (T ) was deduced from the measurement of C p (T ) as detailed in Appendix C. As discussed in the introduction, C p (T ) in type-I clathrates is dominated by the contribution of the optical phonon branches and mostly by the lowest energy guest modes at E 1 such that the γ extracted from Eq. (4) is mainly specific to the contribution of the low-lying guest modes. The change in lattice spacing with temperature in our single crystal of Ba 7.81 Ge 40.67 Au 5.33 was experimentally measured by neutron Larmor diffraction on the most intense Bragg peak (006), and is plotted as black circles in Fig. 2(a). The result is compared to the miniature capacitance dilatometer measurement made by Falmbigl et al. [50] of a type-I clathrate with a similar chemical composition. Figure 2(a) shows the consistency between these two different experimental methods. The temperature dependence of the linear thermal expansion coefficient, α L (T), is then the numerical derivative of a/a, and the volumetric expansion The temperature-dependent molar volume for Ba 7.81 Ge 40.67 Au 5.33 has been taken from the conversion of the coefficient of lattice expansion data in Fig. 2(a) to the experimental lattice parameter, and since we find no change in sound velocity within 300 ± 150 K (see Fig. 3), a temperature-independent B is assumed. For Ba 7.81 Ge 40.67 Au 5.33 , we find B = 65.60 GPa from our measurement of the phononic sound velocities. (More details are given in Appendix B.) The Grüneisen parameter can then be experimentally deduced from Eq. (4). It is found to be temperature independent in the range of interest for this study, with a value of γ = 1.38. A similar method was used by Ikeda et al. [19] for Ba 8 Ge 30 Ga 16 in which γ 300K = 1.67 was observed. In literature, the Grüneisen parameter of type-I clathrates, obtained by various methods and for different chemical compositions, is typically found to be in the range of 1.2-2.0 [16,33,[50][51][52][53][54][55][56][57][58].

A. Polarization dependence of optical branches
The LA and TA phonon dispersions obtained from the experimental mappings (see Appendix F) of the phonon energy at 150, 300, and 530 K covering several Brillouin zones are shown in Figs. 3(a) and 3(b). Measurements have been performed around the most intense Bragg peak (006) in   The experimental phonon spectra are compared to the simulated phonon spectrum obtained by DFT calculations using the SCAN functional (see Methods section), as shown in Figs. 3(a) and 3(b) (gray lines are the simulated phonon dispersions). A good agreement is observed on the whole spectrum especially on the transverse and longitudinal acoustic branches which are well reproduced. That corresponds to a significant improvement in the theoretical approach in comparison to the simulations performed with the PBE functional for which acoustic phonon energies are strongly underestimated [13,15].
For all temperatures and both polarizations, the experimental phonon spectra exhibit an acoustic regime at low energy, which contains well-defined phonon peaks whose dispersions are delimited by low-lying optical bands. The energy at which the acoustic dispersions are interrupted is higher for the longitudinal polarization, 6.5 meV than for the transverse polarization, 4.8 meV. As emphasized by the use of the two sets of experimental data in Figs. 3(a) and 3(b), this is consistent across [001] and [111] longitudinal polarizations, and for [110] and [111] transverse polarizations. The cutoff effect of the acoustic branches by the low-lying optical branches is also seen in the simulated phonon spectrum. The simulation perfectly reproduces the TA dispersion, while for the longitudinal polarization, the computed acoustic branch is interrupted at 4.8 meV and not at 6.5 meV as observed experimentally.
Referring again to the Ba 7.81 Ge 40.67 Au 5.33 crystal structure in Fig. 1(a), the motions of Ba(6d) atoms in the soft plane of the large tetrakaidecahedral cages dominate the optical band centered at 4.8 meV (E Ba ) [59], and hybridized vibrations of Au(6c)-Ba(6d) atoms are thought to dominate the optical band centered at 6.5 meV (E AuBa ) [13,29]. This mode assignment is confirmed by the plot of the partial phonon density of states (pDOS) shown in Fig. 4. Indeed, only the Ba(6d) atoms contribute in the energy range around 4.8 meV while the optical band centered at 6.5 meV contains contributions of the Ba(6d) and the framework Au(6c) atoms. For energies higher than E Ba,AuBa up to a cutoff energy of around 35 meV, the phonon spectrum consists of several broad distributions of optical bands, such as those in the range E Ba -15 meV in Fig. 3(a). The section in momentum at 3.5 meV, in Fig. 5(a), is a cut through the acoustic branches. The Brillouin zones in which the zone center corresponds to an intense Bragg peak result in a strong dynamical structure factor of the acoustic modes and therefore contain two well-defined rings of high intensity. The intensity along the rings is not homogeneously distributed in momentum space as it is weighted by the polarization factor [described by Eq. (3)] such that when the phonon wave vector (q) is aligned/perpendicular to the Bragg wave vector (G), the longitudinal/transverse polarization is selected. Thus, looking at the two rings surrounding the zone center G 006 , the outer ring which is intense along the [110] direction (perpendicular to G 006 ) corresponds to TA phonons, and, reciprocally, the inner ring with maximum intensity along [001] corresponds to LA phonons. The intensity maxima along the rings follow the polarization factor and are rotated by 45 • between those surrounding the zone centers 006 and 222 . This intensity distribution is reproduced on the simulated map, shown in Fig. 6(a), and carries the signature of the coherent character of the acoustic modes.
The section in momentum at a fixed energy equal to E Ba , shown in Fig. 5(b), reveals the intensity distribution of the lowest optical band which cuts the TA dispersion. It also contains the contributions of the LA modes which form the inner rings closest to the zone centers. The intensity of the optical band at E Ba shows a distinct Q-dependence in momentum space with intensity maxima at the zone boundaries that can be associated with the Bragg intensity at a zone center, for instance, (1.5 1.5 6) and (1.5 1.5 3.5) with Bragg peaks (006) and (222), respectively. The intensities along these maxima follow the transverse polarization and are turned by 45 • between the zone centers 006 and 222 as for the acoustic modes. This intensity distribution of the outer ring related to the optical mode at E Ba is also reproduced in the simulated map shown in Fig. 6(b). It confirms the coherent nature of the optical modes contained in this band which thus cannot be associated with vibrations of localized and independent Einstein oscillators, in agreement with previous phonon studies on clathrates [14] and other cage compounds [60]. Indeed, an isolated Einstein oscillator-type behavior would result in a Q 2 -dependent intensity distribution. Note that there are also intensity maxima at zone centers such as (116) and (114), which is not seen in the simulated map shown in Fig. 6(b) simply because the intensity is in the low range of the color scale.
Moving again in energy to E AuBa in Fig. 5(c), one sees the maxima of intensity related to the distribution of the optical band which cuts off the LA dispersion, in addition to the ones related to E Ba which are pinned at the zone boundary. The intensity distribution to E AuBa around 006 is more spread in momentum space than that at E Ba which is also observed in the simulated map shown in Fig. 6(c). Looking carefully at the experimental map, one can distinguish minima along the [110] direction following the longitudinal polarization factor which are less obvious to see in the simulated map. The optical vibrations at E AuBa are much less coherent than at E Ba and exhibit the trend of being polarized longitudinally.  Fig. 3, and the INS measurements of Ref. [59]. E Ba,AuBa and E 3,4 correspond to the peaks labeled in Fig. 7(e), and thermal expansion to the fit TE labeled in Fig. 7(f). Theoretical values from the Self-Consistent Phonon method, which includes the cubic and quartic contributions (SCPB), have been reported [43] and SCPB-based thermal expansion is also given, all in 1 × 10 −4 meV/K.
Thus, the energy dispersions and momentum distributions shown in Figs. 3 and 5 demonstrate that the low-lying optical bands interact with the acoustic dispersion that is present in a given polarization. It appears that the TA dispersion couples largely with the 4.8 meV band (E Ba ), while the LA dispersion couples mainly with the 6.5 meV band (E AuBa ). Last, Fig. 5(d) is a experimental data integration in the interval 7-12 meV, which represents host optical band energies. A similar momentum map obtained by integrating the simulated data in the same interval is shown in Fig. 6(d). The intensity distribution of these modes on the experimental and simulated maps reveal the coherence of these host bands, which display mostly longitudinally polarized intensity.

B. Temperature dependence of optical branches
Focusing now on the temperature dependence, and in contrast to the acoustic phonon energies for which no change of their energies in temperature is observed, we clearly see sizable shifts with temperature of the optical band energies in Figs. 3(a) and 3(b). For the transverse polarization, the energy E Ba increases upon heating, while the energies of the other optical bands, including E AuBa , follow the opposite trend. As a consequence, the temperature-dependent changes of the energy range for acoustic phonons is opposite for TA and LA phonons as it is directly related to the temperature dependence of the optical cutoff bands E Ba,AuBa [see Eq. (1)]. This is better seen by comparing the GVDOS between 150 K and 530 K, depicted in Fig. 3(c). In the low-energy range below 10 meV, the GVDOS exhibits mostly three peaks at energies E Ba , E AuBa , and E 3 . Only E Ba follows a hardening shift upon heating, which has been similarly reported in type-I clathrates of different chemical compositions [59,[61][62][63][64]. The overall structure of the measured GVDOS is reproduced well by the phonon DOS obtained from ab initio simulations using the meta-GGA SCAN functional, as shown in Fig. 3(c). Some differences appear in the energy range between 5 and 7 meV, which might be linked to those observed for the LA dispersions in Figs.

3(a) and 3(b).
We then more systematically investigated these temperature dependencies by use of a triple-axis spectrometer between 200 and 500 K. Performing energy scans at constant wave-vector Q = (113), where the acoustic phonon intensity is expected to be very low, allows for a more selective study of the optical branches, as seen in Figs. 7(a)-7(d). E AuBa is not visible in this polarization. However, following the GVDOS peaks in Fig. 3(c), we find a similar E 3 = 7.5 meV peak, and a higher energy peak at E 4 = 9.8 meV. The energy fits of E Ba and E 3,4 , along with E Ba,AuBa and E 3 from the GVDOS, are then plotted in Fig. 7(e). All bands display linear trends in the overall temperature range, with only E Ba hardening with a rate of 7.0 × 10 −4 meV/K. For comparison, the 9.3 × 10 −4 meV/K temperature dependence of E Ba for Ba 8 Ge 30 Ga 16 , measured by inelastic neutron scattering [59], has been included as well in Fig. 7(e). The rates of softening of the other optical bands in Ba 7.81 Ge 40.67 Au 5.33 are summarized in Table I. 013021-7

V. DISCUSSION
From Eq. (1), it is evident that the reproduction of the experimentally observed κ L (T) can be achieved only if the temperature dependence of the low-lying phonon modes is correctly accounted for. Therefore, we now seek direct comparison of the experimentally observed temperature dependence of the E Ba band to results from SCP theory. The combined cubic and quartic anharmonicity terms give the isochoric contribution of the thermal change of phonon energy, ( ∂E ∂T ) V , while experimental measurements, such as the ones presented above, are usually performed at constant pressure and give access to the isobaric thermal variation, ( ∂E ∂T ) P . To the first order, the thermodynamic relation between those quantities is given by where the far right term corresponds to the mode-specific thermal expansion which we have experimentally quantified in Ba 7.81 Ge 40.67 Au 5.33 and Ba 8 Ge 30 Ga 16 through their volumetric thermal expansion coefficients, α V (T), average Grüneisen parameters, γ , and energies of the mode under investigation (with the help of α V and inelastic neutron scattering measurements in Refs. [19,59] for Ba 8 Ge 30 Ga 16 ). While experimentally we cannot further separate the cubic from the quartic anharmonic term, we discuss a conceptual first approximation in Appendix E. In Table I  Conversely, for the SCP simulations, the SCP calculation which includes the quartic and cubic contributions (SCPB) with a rate of 14.8 × 10 −4 meV/K can be used along with γ SCPB (T ) in order to find the equivalent isobaric rate of 11.2 × 10 −4 meV/K. (More details are given for SCPB in Appendix D.) We therefore find close isochoric matching between SCPB and experimentally deduced anharmonicity in Ba 8 Ge 30 Ga 16 , allowing us to experimentally validate the SCPB method in clathrates.

VI. CONCLUSIONS AND PERSPECTIVES
In summary, we confirm the good agreement between our measurements and the ab initio harmonic calculations using the meta-GGA SCAN functional of the overall phonon spectrum along different directions and for the transverse and longitudinal polarizations on the type-I clathrate Ba 7.81 Ge 40.67 Au 5.33 . However, experimentally, the TA and LA branches are delimited by two optical phonon bands of different nature which is not reproduced by our simulations. While the former hybridizes with transverse optical vibrations centered at E Ba = 4.8 meV associated with the coherent guest motions of Ba(6d) in the soft plane of the tetrakaidecahedral cages, the latter is interrupted by the longitudinal optical band centered at E AuBa = 6.5 meV related to the coherent hybridized motions of the substituted Au(6c) host and the Ba(6d) guest atoms. The TA dispersion is perfectly reproduced by our ab initio simulation in the whole Brillouin zone. On the other hand, a difference appears on the longitudinal branch in the region of the phase space where optical and acoustic modes hybridize such that the LA branch is predicted to be cut off at the same energy as the transverse one.
Upon cooling, E AuBa increases, following a rate of change guided by thermal expansion, while E Ba decreases. The experimental isochoric rate of change of E Ba extracted in Ba 8 Ge 30 Ga 16 , which quantifies the amount of cubic and quartic anharmonicity, is in good agreement with the SCPB simulations in Ba 8 Ge 30 Ga 16 . In Ba 7.81 Ge 40.67 Au 5.33 , a much lower rate of change is found which is assumed to originate from the difference of the defect structure. On the other hand, the difference between the transverse and the longitudinal cutoff energies observed experimentally, which is not seen in the simulations, indicates a more subtle polarization dependent mechanism involving either the disorder and/or the anharmonic polarization mixing whose consideration in the SCP-like approach will surely reveal an improved κ L (T ) dependence, especially in the intermediate temperature range as seen in Fig. 1(b).
Data from inelastic neutron scattering measurements at the ILL are available at [65], and LLB measurements correspond to proposal 657.  using an incident neutron wavelength of λ = 3.2 Å and a cryofurnace. Scans covered an range of 54 • , with a sample rotation of 0.5 • in between each scan. Time-of-flight data were reduced using Mantid [68] and then processed into a four-dimensional S(Q, E ) file and further analyzed with the Horace package [69] under Matlab. Considerable effort was made to understand q and E resolution (dq, dE) of the instrument in order to make realistic data integrations with Horace. The overall instrumental resolution is governed by the incoming neutron beam energy and monochromatization, the beam divergence in horizontal and vertical directions, the sample mosaic, the receiving slit sizes in front of the detector, and the final neutron energy [70]. This then must be compared to the minimum step size that can be achieved with the instrument, as explained hereafter.
By taking into account the divergence due to the size of the detector tubes of IN5@ILL and the beam divergence while in the λ = 3.2 Å condition, we have defined the resolution for the [001] and [110] directions while near the (006) Bragg peak. More specifically, detector tubes on IN5 have a diameter of 2.54 cm, meaning that with a distance of 4 m between the sample and detectors, there is a beginning divergence of 0.36 • due to instrumental conditions. To this we add the corresponding horizontal and vertical beam divergences, 0.64 • and 0.96 • , respectively, due to the incident neutron wavelength and neutron guide horizontal and vertical super-mirror indices on IN5 [70]. Using these starting points, we then calculated the local resolution limits around our Bragg peak of interest, the (006). These are summarized in Table II. Next, as a first approximation, we consider that q and E are decoupled for a time-of-flight spectrometer, unlike on a triple-axis spectrometer. This means that the effective phonon energy resolution depends on the effect of energy broadening due to the reciprocal space, E sound vel , and on the instrumental resolution for a given energy transfer, E instr . The latter has been calculated using the incident neutron wavelength, speed of the IN5@ILL choppers (12 000 rpm), and additional spectrometer parameters [70]. For our given experimental conditions, Fig. 8 depicts the energy resolution as it changes with the energy of the scattering event.
The former, on the other hand, refers to E sound vel = v s q, where v s is the sound velocity of the particular dispersion The minimum step size that can be chosen when preparing the 4D data matrix in Horace is therefore given by the step sizes due to detector pixel sizes, step rotation, and step size in energy. These steps are generally smaller than the effective instrumental resolution. It should be noted that dq T , or the transversal dq step size, is strongly Q dependent, since it is equal to (rotation) × Q, where (rotation) is the rotation step size expressed in radians. As stated above, there is a 0.5 • step size for these experimental data. The final Horace integration limits are derived from Table II, in which the overall dq resolution is: ( q step div. ) 2 + ( q beam div. ) 2 . Figure 12 contains two-dimensional plots showing the energy spectrum in a given direction within the momentum plane were cut and then folded along symmetries using the appropriate Horace functions. For the one-dimensional scans around the (222) Bragg peak in Fig. 3(b), the same procedure as described above has been completed to reflect the resolution limits near this new position in reciprocal space. As there is less intensity in this region, however, integrations slightly larger than those of the strict minimums were taken: dq = 0.08 r.l.u. along the propagation axis [ζ ζ 0], 0.08 r.l.u. along [00ξ ], and 0.06 r.l.u. along the out-of-plane axis [ζ ζ 0]. The 1D scans used to create Fig. 3 In addition, high-resolution time-of-flight measurements were taken with λ = 4.8 Å in order to obtain the neutronweighted generalized vibrational density of states (GVDOS) plot in Fig. 3(c), calculated by the MUPHOCOR (MUlti-PHOnon CORrection) routine [71] in the LAMP program [72] for single crystal data. Scans covering an range of 46 • at 150 K and 34 • at 530 K, both with sample-rotation step sizes of 2 • , were used to make the calculation, and the atomic mass and expected neutron scattering cross section of Ba 7.81 Ge 40.67 Au 5.33 were used as initial starting parameters. Preliminary experiments, performed on the triple-axis spectrometer 1T@LLB with a fixed k f = 2.662 Å −1 , were extremely important in this work. Indeed, the polarization dependence of the low-lying optical phonon band at E Ba was formerly observed in these experiments.

Neutron Larmor diffraction
Neutron Larmor diffraction takes advantage of the neutron resonance spin echo technique for triple-axis spectrometers in order to measure the change in lattice spacing with temperature with extreme sensitivity (1.5 × 10 −6 ) [73][74][75][76]. Such a measurement was made on the (006) Bragg peak position between 3 and 300 K, and is shown in Fig. 9(a) (also see Fig. 2). Measurements were taken in increments of 5 K with a fixed k i of 2.13 Å −1 on the thermal-neutron triple-axis spectrometer TRISP@FRM-II.  [50]. Its numerical derivative is used to find α L = α V /3. The measurement of C P (black circles) from Ikeda et al. [19] was used to calculate C V (solid blue line).
Using these elastic constants, we can calculate the material bulk modulus at 0 K, B, through Eq. (B2), which gives us B = 65.60 GPa. Again, the same value is found for Ba 8 Ge 30 Ga 16 with almost no temperature dependence [19,80]: 12 3 . (B2)

APPENDIX C: CALCULATION OF EXPERIMENTAL GRÜNEISEN PARAMETER (CONTINUED)
This section expands on the experimental Grüneisen parameter results discussed in the main text, in particular on Eq. (4). Figures 9(a) and 9(b) have already been presented in Fig. 2. From Ikeda et al. [19] we have experimental C P for Ba 7.81 Ge 40.67 Au 5. 33 . To calculate C V , we use the relation [82], and this is plotted in Fig. 9(c). Finally, using the given equation for calculating the average Grüneisen parameter, γ (T ) is calculated and plotted in Fig. 9(d). This method of experimental deduction gives γ 300K = 1.38, which is temperature independent within the temperature range of study.

APPENDIX D: COMPARING TO SCP METHOD CALCULATIONS
There are two different calculations of the E Ba optical band in Ba 8 Ge 30 Ga 16 that have been summarized by Tadano and Tsuneyuki in the Supplementary Material of Ref. [43]: (1) SCP: SCP method calculations in which the real part of the quartic anharmonic term, or the Loop free energy diagram, determines phonon frequency renormalizations in the system, and (2) SCP+Bubble: the same calculation, but made by also including the real part of the cubic anharmonic term, the Bubble free energy diagram. (The second calculation will be referred to as SCPB in this discussion.) We therefore reflect that the SCPB calculation directly provides us with ( ∂E ∂T ) V of E Ba and that the self-energy terms in the SCPB calculation self-consistently include the renormalized phonon energies including the effect of the quartic anharmonicity. (This is an important distinction for Appendix E, as this means that we cannot extract the harmonic Bubble term but only the self-consistent Bubble term.) In order to move between isobaric and isochoric representations of the SCPB calculation discussed in the main text, we have calculated the SCPB-based thermal expansion: E SCPB (T )α V γ SCPB (T ). In this way, the determination of the thermal expansion is much more exact than by using the quasiharmonic approximation for which only the harmonic phonon energies are considered. As will be shown in the following paragraphs, γ SCPB (T ) and E SCPB (T ) come from the information given in the Supplementary Material of Ref. [43], and as a first approximation we use the experimental α V for Ba 8 Ge 30 Ga 16 from Ikeda et al. [19].
The isochoric rate of change in energy of E Ba with temperature for the SCP calculation is 17.9 × 10 −4 meV/K, and 14.8 × 10 −4 meV/K for the SCPB calculation. Next, with the use of the mode-specific definition of the Grüneisen parameter, the SCPB-based Grüneisen parameter for the E Ba mode can be calculated using Fig. 10(b). The mode-specific Grüneisen parameter for the SCP calculation [ Fig. 10(a)] is 4.71, 2.57, 2.06, and 1.76 for temperatures of 0, 300, 600, and 900 K, respectively, and in similar fashion for the SCPB calculation [ Fig. 10(b)]: 4.81, 2.98, 2.34, and 1.94. Although these Grüneisen parameters appear to be slightly larger than the average clathrate experimental values discussed in the main text, we note that these parameters reflect the anharmonicity found specifically in the E Ba optical band. With these mode-specific Grüneisen parameter results, the SCPB thermal expansion shown in Fig. 7(f) was calculated and, by extension, the isobaric form of the SCPB calculation was found to have a rate of change of 11.2 × 10 −4 meV/K.
We also note that there are experimental rates of change for Ba 8   a discrepancy between Raman measurements of Ba 8 Ge 30 Ga 16 measured by this group and by those of Christensen et al. [59]. We emphasize, however, that the isobaric 17.1 × 10 −4 meV/K should not be compared to the isochoric quartic-only calculation.

APPENDIX E: FIRST APPROXIMATION DECOUPLING OF CUBIC AND QUARTIC ANHARMONICITY
Even though we cannot definitively isolate the cubic and quartic anharmonicity terms from an inelastic neutron scattering ( ∂E ∂T ) P measurement, we will attempt to qualitatively interpret their weighted importance on E Ba .
The generalized vibrational density of states (GVDOS), as seen in Fig. 3(c), represents the isobaric temperature dependence of Ba 7.81 Ge 40.67 Au 5.33 . However, if the energy axis of the data at 530 K is scaled by 3% of the original values, as seen in Fig. 11, then globally, all peaks except for E Ba now align. This is already a powerful conclusion about the uniqueness of E Ba , pointing to the strong anharmonicity that governs its behavior as opposed to all other higher energy peaks, for which one scaling factor can explain the complete temperature dependencies.
This 3% scaling factor can be understood to be made up of the thermal expansion and cubic components, since FIG. 11. Generalized vibrational density of states obtained from the data collection recorded on IN5@ILL at 150 and 530 K. The energy axis of the data at 530 K has been scaled by 3%, leading to an alignment of all but the lowest peaks.
we will assume that the quartic term is localized onto only E Ba as a first approximation. This simplification is supported by Tadano and Tsuneyuki [43], who find that the quartic anharmonic phonon energy renormalizations of the SCP (not SCPB; see Appendix D for more details) calculation are focused onto modes below 9.92 meV (80 cm −1 ) for Ba 8 Ge 30 Ga 16 . We therefore extrapolate to say that modes higher than E Ba are controlled by only thermal expansion and the cubic term, while E Ba has contributions from thermal expansion, the cubic term, and the quartic term.
Let us first look at GVDOS peaks between 10 and 35 meV, for which we must consider only thermal expansion and cubic components. To reiterate, by studying the 3% scaling factor, we are studying the isobaric rate of change for peaks between 10 and 35 meV. We recall that, using Eq. (5), the peak-specific thermal expansion for each of the GVDOS peaks between 10 and 35 meV in Fig. 11 can be calculated, and the result is that thermal expansion consistently accounts for 45%-50% of the total isobaric rate of change. Given our assumption that these modes are controlled only by the thermal expansion and cubic terms, this leaves a remaining 50%-55% of the 3% scaling which must be understood as the rate of change of the cubic anharmonicity. In this manner, the 3% scaling factor of peaks between 10 and 35 meV is fully accounted for.
Expanding on this 3% scaling concept to E Ba , we recall the values given in Table I: E Ba has an isobaric rate of change of 5.4 × 10 −4 meV/K and a thermal expansion rate of change of −2.7 × 10 −4 meV/K, giving the isochoric rate of change of 8.1 × 10 −4 meV/K. However, this time we expect for the isobaric rate of change to reflect thermal expansion, cubic, and quartic contributions. As stated above, the rate of thermal expansion for E Ba is already known, leaving the cubic and quartic contributions, for which we also know the total isochoric (cubic and quartic) contribution. However, after the 3% scaling shown in Fig. 11, the rate of change for E Ba becomes 9.4 × 10 −4 meV/K. Therefore, the difference between the 3% scaling factor and isochoric rate of change, which is −1.3 × 10 −4 meV/K, must be due to cubic anharmonicity, meaning that 9.4 × 10 −4 meV/K is the quartic contribution to E Ba .
While we cannot directly compare the Bubble contribution of the "SCP+Bubble" (SCPB) calculation (see Appendix D) to the cubic contribution using this GVDOS rescaling method, we note that −1.3 × 10 −4 meV/K has the correct sign for the cubic component [43].  Fig. 12, and they were used to construct Fig. 3(a). Figures 13(a)