Uncovering design principles for amorphous-like heat conduction using two-channel lattice dynamics

The physics of heat conduction puts practical limits on many technological ﬁelds such as energy production, storage, and conversion. It is now widely appreciated that the phonon-gas model does not describe the full vibrational spectrum in amorphous materials, since this picture likely breaks down at higher frequencies. A two-channel heat conduction model, which uses harmonic vibrational states and lattice dynamics as a basis, has recently been shown to capture both crystal-like (phonon-gas channel) and amorphous-like (diffuson channel) heat conduction. While materials design principles for the phonon-gas channel are well established, similar understanding and control of the diffuson channel is lacking. In this work, in order to uncover design principles for the diffuson channel, we study structurally-complex crystalline Yb 14 (Mn,Mg)Sb 11 , a champion thermoelectric material above 800 K, experimentally using inelastic neutron scattering and computationally using the two-channel lattice dynamical approach. Our results show that the diffuson channel indeed dominates in Yb 14 MnSb 11 above 300 K. More importantly, we demonstrate a method for the rational design of amorphous-like heat conduction by considering the energetic proximity phonon modes and modifying them through chemical means. We show that increasing (decreasing) the mass on the Sb-site decreases (increases) the energy of these modes such that there is greater (smaller) overlap with Yb-dominated modes resulting in a higher (lower) thermal conductivity. This design strategy is exactly opposite of what is expected when the phonon-gas channel and/or common analytical models for the diffuson channel are considered, since in both cases an increase in atomic mass commonly leads to a decrease in thermal conductivity. This work demonstrates how two-channel lattice dynamics can not only quantitatively predict the relative importance of the phonon-gas and diffuson channels, but also lead to rational design strategies in materials where the diffuson channel is important.


Introduction
The flow of heat through solids is a topic of wide-spread technological importance. In simple crystals the thermal conductivity (κ) can be well understood within the phonon-gas model (PGM) 1 . However, for amorphous materials [2][3][4] , defective crystals 5 , anharmonic crystals 6 , hybrid perovskites 7 , and -as we show here -Zintl-type phases with complex crystal structures, the PGM is incomplete. Two-channel lattice dynamics simulations, which has been developed by Simoncelli et al. 6 and Isaeva et al. 4 and is applied in this work, has emerged as an effective way to quantitatively capture both crystal-like and amorphous-like thermal transport in solids. 8,9 The thermal conductivity of all solids spans roughly five orders of magnitude (Figure 1a). The PGM tends to dominate in crystalline-like heat conductors and, in general, can explain the top two to three orders of magnitude. However, there is far less understanding of the detailed physics governing heat conduction in amorphous-like heat conductors which tend to occupy the bottom two to three orders of magnitude. To build intuition regarding amorphous-like heat conduction we study the vibrational and thermal properties of crystalline Yb 14 MSb 11 (M = Mn or Mg) which has a very complex crystal structure containing 104 atoms in its primitive unit cell ( Figure 2a). It is a champion p-type material system for high temperature thermoelectric energy conversion and is well studied experimentally and computationally [10][11][12][13][14] . Yb 14 MSb 11 has a κ below that of amorphous fused silica even though it is indeed a crystal. Additionally, its temperature dependence of κ above room temperature is relatively flat, similar to that of amorphous materials, rather than κ vib ∝ T −1 which is characteristic of the PGM (crystalline-like). From this case study emerges a simple physical picture for heat conduction which enables the design of materials that exhibit crystalline-like conduction, amorphous-like conduction, as well as materials transitioning between the two extremes. orders of magnitude (above κ vib ≈ 50 W/mK) are well described by the phonon-gas model, exhibiting crystalline-like heat conduction. The phonon-gas model alone is insufficient for materials in the bottom 2 to 3 orders of magnitude. b) A one dimensional illustration of band folding when the unit cell is increased in size by, for example, modulating the interatomic force constants. The overlap of broadened modes enables heat conduction through the diffuson channel and is highlighted in red. c) The thermal conductivity viewed as a matrix with rows and columns designating the phonon branch index (s and s ). The phonon-gas channel is simply the sum of the diagonal components of this matrix. The diffuson-channel is the sum of off-diagonal components. Heat conducting though these two channels have very different physical interpretations which, as described in the text, the trucks illustrations help conceptualize.
We analyze the vibrational and thermal properties of Yb 14 MSb 11 within a two-channel lattice dynamical (LD) approach κ vib = ∑ qss C ss (q)v ss (q)v s s (q)τ ss (q). (1) Here, we follow the work of Hardy 15 , Allen et al. 16 , Simoncelli et al. 6 , Isavea et al. 4 , and Luo et al. 8 and we consider not only the phonon-gas channel (s = s , where s is the phonon branch index), but also heat conduction through the coupling of eigenstates with the same phonon wavevector q (i.e. between qs and qs ) which we call the diffuson channel.
In principle coupling between eigenstates with different q-vector should be possible (i.e. between qs and q s ), but has not been considered to date and is also neglected here. Any mechanism that limits phonon lifetime or spatial extent will lead to a q-broadening, which is inseparable from the more commonly considered energy broadening. It is unclear at the moment how important coupling between q-vectors will be. We can say that approaching the amorphous limit, we expect qs to qs coupling to be sufficient for capturing the diffuson channel since the Brillouin zone effectively collapses to a single q-point. As for Yb 14 (Mn, Mg)Sb 11 , we expect q to q coupling to be relatively less important than in simpler crystals (larger Brillouin zones) with similarly low thermal conductivity that have been considered using two-channel lattice dynamics, such as CsPbBr 3 6 .
Each term in the thermal conductivity (Eq. 1) is represented by its appropriate matrix, and they are defined in Computational methods section. C ss (q) is the heat capacity matrix; v ss (q) is the velocity matrix; and τ ss (q) is the lifetime matrix. The diagonal and off-diagonal components of these matrices have distinctly different physical interpretations. A schematic of this formalism is shown in Figure 1c.
The diagonal components of this matrix (s = s ) designate the commonly used phonon-gas model where each of the matrices reduce to the common phonon property, i.e. the heat capacity C s (q) = C ss (q) (dropping the redundant branch index subscript), the phonon group velocity v s (q) = v ss (q), and the phonon relaxation time τ s (q) = τ ss (q). The off-diagonal components (s = s ) we call the diffusonchannel since it is mathematically the same mechanism through which diffusons were defined 2, 17 , and represents a coupling between vibrational states which results in diffusive heat conduction Therefore, κ ph and κ diff simply sum together to give the total thermal conductivity due to atomic vibrations (i.e. the 2/13 lattice thermal conductivity) There have been several reports of two channel models for κ vib where one of the additive terms is the phonon-gas model and the other is an analytical model for the amorphous/minimum thermal conductivity [18][19][20] . For the amorphous thermal conductivity specifically, there have been a number of analytical models proposed [21][22][23][24] . It seems as though these analytical (or semi-analytical) treatments have been attempting to capture the effect of the diffuson channel defined in Eq. 3.
The diffuson channel, κ diff represents the diffusion of thermal energy through quasi-degenerate normal modes 4 , or normal mode mixing (also referred to as coherences 6 ). Normal mode mixing and κ diff can be understood through the concept of band folding which we schematically show for one dimension in Figure 1b. As the number of atoms per unit cell is increased, the bands are folded resulting in a smaller maximum q-vector, and more branches occupying the same qspace. Each phonon branch has a linewidth, which is inversely proportional to its lifetime, and may arise from the materials natural anharmonicity and/or defects. When the phonon branches are packed close enough together in energy and/or their linewidths are large enough such that the modes at the same q-point overlap in energy (highlighted in red in Figure  1b), normal mode mixing can occur. This is mathematically reflected by the Lorenzian form of τ ss (q) (Eq. 14). Given that normal mode mixing 'can' occur (τ ss (q) > 0), the velocity matrix determines whether or not it 'will' occur (v ss (q) > 0). For off-diagonal components, v ss (q) captures how strongly phonon mode qs and qs will be coupled, based on the mode shapes (i.e. eigenvectors, polarizations). We note that just as harmonic phonon eigenstates are used as a basis to compute the thermal conductivity in standard lattice dynamics and κ ph , they are also used as a basis to compute κ diff . In both cases, the real system can not be described within the harmonic approximation, and anharmonic contributions are required to obtain physical results (i.e. a finite κ ph and κ diff > 0). Therefore, using precise language, κ ph and κ diff are defined as the thermal conductivity of harmonic phonon eigenstates which are conducting heat through the phonongas and diffuson channels, respectively. Modes which have a significant κ diff contribution are expected to quickly loose there phonon-like, propagating plane-wave character if the true potential energy landscape was used to obtain vibrational properties.
A useful metaphor for conceptualizing the phonon-gas channel is that of trucks on a highway (Figure 1c). The truck represents the phonon mode; its cargo represents the heat capacity; its speed represents the phonon group velocity; and collisions between trucks represent phonon scattering events. If one desires an equally simplistic view of the diffuson channel, this metaphor can be extended. When the diffuson channel is important there are many trucks on the road, and presumably the majority of the trucks have a low velocity and collisions may be frequent. However, within the diffuson channel, packages can shuffle between trucks. The proximity of two trucks may represent the lifetime matrix, τ ss (q), and whether or not the cargo will fit into the receiving truck (i.e. shape of the box or the mode shape) represents the velocity matrix v ss (q). Just as in random walk theory, if there is a gradient in cargo density (i.e. a temperature gradient) there will be a net diffusion of cargo. Both mechanisms are always present and contribute additively to cargo conduction. However, one may dominate over the other depending on the state of traffic.
The two channel nature of vibrational thermal transport stems back to Hardy's early work (see Eq. 2.14 and Eqs. 3.29 to 3.31 of Ref. 15 ), and recently large strides have been made applying it in a lattice dynamical context by Simoncelli et al. 6 and Isaeva et al. 4 . Hardy's heat current operator is commonly used along with the Green-Kubo method to extract κ vib from molecular dynamics (MD) simulations. MD is an invaluable tool in pioneering heat conduction beyond the phonon-gas model 3,5,25,26 . Since MD uses Hardy's heat current operator directly, and does not explicitly assume the phonon-gas model for heat conduction, it naturally captures other mechanisms of vibrational heat conduction. However, the atomistic nature of MD approaches comes with limitations of interpretation. In an attempt to interpret MD-based computational experiments of amorphous materials, a 'descriptive' taxonomy has been established to classify vibrational modes as propagons, diffusons, or locons. Quantitative metrics have even been defined to characterize the nature of heat conduction 3 . However, an intuitive physical picture is critical for the 'predictive' design of defective and complex crystals, as well as amorphous materials 26 . Standard lattice dynamics (LD) (as computed by common computational suites such as phono3py 27 ), by design, only captures the phonon-gas channel. However, as we outline in Figure 1 the two-channel LD approach 4, 6 captures vibrational heat conduction beyond the phonon-gas model. One benefit of the two-channel LD framework is that it can provide an intuitive physical picture which can help facilitate predictive materials design. MD approaches tend to behave as

3/13
computational experiments and thus provide a more descriptive understanding. Both approaches are important.
It might be helpful to consider early simulations on diffusons in amorphous Si 2, 16, 28 as an extreme case of band folding illustrated in Figure 1b, where a 1000+ atom simulation cells are used with periodic boundary conditions. This results in extreme band folding and a very small maximum q-vector (measured from the Γ-point) meaning the q-vector, to some extent, looses its meaning. Therefore, theoretical treatments of amorphous materials tend to collapse down to spectral treatments where properties are expressed only as a function of energy, rather than as wavevector and energy. As we show here, Yb 14 MSb 11 seems to be somewhere in between these extremes.

Band structure and density of states
The phonon band structures of Yb 14 MgSb 11 and Yb 14 MnSb 11 , obtained from density functional theory (DFT) based lattice dynamics are shown in Figure 3. The many crossing points and close proximity of bands suggests that there should indeed be a significant diffuson channel contribution to κ vib . Here, we examine the computed vibration density of states and compare this to experimental inelastic neutron scattering data.
Dense polycrystalline pellets of Yb 14 (Mg, Mn)Sb 11 were synthesized (see the Synthesis section and Figure S2) and the vibrational properties were measured using time-of-flight neutron spectroscopy (details in the Inelastic neutron scattering section). The neutron scattering intensity measured versus energy and momentum transfer are shown in the left panels of Figures 3a and c. The right panels show the neutron scattering strength simulated using the DFT-based harmonic lattice dynamics, displayed in Figures 3b and d. The harmonic band structure is plotted along the Γ-N and Γ-X directions (labeled in Figure 2b). The atom projected and total vibrational density of states (DOS) are shown next to the band structures. These have been convoluted with a 1.5 meV Gaussian which is consistent with the experimental instrument resolution.
Starting at the top of Figure 3b, Yb 14 MgSb 11 is expected to have very flat high energy optical modes at around 30 meV, which are dominated by Mg character with a small amount of Sb character. Though relatively faint, this mode is clearly observed in the neutron scattering data and simulation in Figure  3a. This can be attributed to the Mg atoms oscillating at high frequency within their Sb 4 tetrahedron cage while interacting very little with the rest of the atoms (see Figure S4a). The same type of behavior is expected in Yb 14 MnSb 11 ( Figure  S4d The number of optical branches can be expected from the concept of zone folding shown in Figure 1b. With N=104 atoms per primitive unit cell, Yb 14 MSb 11 has 3N-3=309 optical modes. What may not have been obvious without the DFT-based LD simulations is their close proximity in energy, which as we outline in Figure 1c and the related text, is important for diffuson channel conduction. The LD simulations also indicate that the mode overlap observed in Figures 3a and c is not expected to disappear if large enough single crystals are able to be grown and measured. Additionally, we conclude that mode overlap observed in Figures 3a and c is not a consequence of the 1.5 meV resolution as demonstrated by the data obtained in a more narrow energy range (25 meV) with better resolution, see Figure S1.
These results are now compared to inelastic neutron scattering obtained density of states. However, before a direct comparison can be made the atom projected density of states (g i (E)) must be weighted by their neutron scattering strength, resulting in the neutron weighted density of states N i is the number of atoms of species i in the formula unit, σ i is its neutron scattering cross-section, and M i is its atomic mass. The results of this comparison are shown in the rightmost panels of Figure 3c and d. Note, the maximum energy is scaled by 96% for comparison of computational and experimental neutron-weighted density of states. The shift in energy is expected due to the use of the PBE functional which is known to slightly over-estimate lattice parameters and slightly under-estimate elastic moduli and vibrational mode energy. This level of agreement is consistent with current state of the art methods 29 . We also note that temperature effects were not considered when computing inter-atomic force constants and that the experimental data shown is at 300 K. The standard error for the experimental data is shown by the error bars, which is equal to or smaller than the data point size. The data agree well with the simulations, especially below 20 meV. In particular we note the agreement in the relative magnitude and shape of the peaks at around 8 meV and modes higher in energy. Whereas, the Mg peak seems to be stronger than predicted, the relative position in energy seems to be well represented. This level of agreement lends confidence to the band structure results shown in Figure 3, which suggest that Yb 14 MSb 11 will exhibit significant heat conduction through the diffuson conduction channel.

Thermal conductivity
Here, we apply the two-channel lattice dynamics approach given in Eq. 1 to more quantitatively estimate the contributions of the phonon-gas and diffuson channels in Yb 14 (Mn, Mg)Sb 11 . Additionally, we demonstrate how modifying the energy of phonon modes by varying atomic mass can All DOS calculations (atom and neutron weighted) shown are convoluted with a 1.5 meV full width half max Gaussian function which is consistent with the instrument resolution seen in the inelastic neutron experiments. Note, the maximum energy scales for the computational and experimental neutron-weighted density of state data are scale slightly (by 96%) for comparison.
systematically control κ diff . Our thermal conductivity measurements and simulations for Yb 14 (Mn, Mg)Sb 11 , displayed in Figure 4, indicate that the diffuson channel dominates above approximately 300 K. Yb 14 MnSb 11 exhibits exceptional thermoelectric properties above 800 K 10,11,13 . Details regarding the measurement of thermal conductivity, calculation of the vibrational thermal conductivity (κ vib ), and a discussion regarding experimental error are given in the 'Transport measurements' section.
The mathematical and computational details for our implementation of two-channel lattice dynamics are described in the Computational methods section (Eqs. 13 to 15). To summarize, we first compute harmonic inter-atomic force constants (IFCs) using density functional theory. Next, we take these IFCs and compute the dynamical matrix (D(q)) and its gradient (∇ i q D(q)) on a uniform q-mesh to sample phonon properties throughout the Brillouin zone. The eigenvalues and eigenvectors obtained by diagonalizing D(q), along with ∇ i q D(q) are then used to compute the required matrices: the heat capacity matrix C ss (q) (Eq. 13), the velocity matrix v ss (q) (Eq. 9), and the lifetime matrix τ ss (q) (Eq. 14). Finally these matrices are combined to compute the thermal conductivity matrix at each q-point κ ss (q) = C ss (q)v ss (q)v s s (q)τ ss (q).
In order to examine the results obtained via two-channel lattice dynamics, it can be helpful to sum over q revealing which phonon branches mix and how much they contribute to κ vib , κ ss = ∑ q κ ss (q).  1). a) The temperature dependence of the phonon-gas channel κ ph , diffuson channel κ diff and the vibrational (lattice) thermal conductivity κ vib = κ ph + κ diff for Yb 14 MnSb 11 . The analytical scattering rate coefficients are fixed using data below 50 K from Ribeiro et al. 30 , and the simulation is then compared to data above room temperature measured in this work and provided by Toberer et al. 11 . b) Vibrational thermal conductivity of Yb 14 MgSb 11 . In both (a) and (b), solid lines indicate τ −1 pp = AT and dashed lines indicate τ −1 pp = A ω 2 T . c) The thermal conductivity matrix κ ss (Eq. 7) for Yb 14 MnSb 11 at 1000 K (the same for Yb 14 MgSb 11 is given in Figure S10). Histograms of the for the diffuson and phonon-gas channels are shown below. The significant population of the off-diagonal (s = s ) elements throughout the density of states results in a κ diff -dominated κ vib . κ ph contributions are isolated to the low index branches, which are acoustic in character.
Finally, the vibrational thermal conductivity is obtained by summing over s and s , κ vib = ∑ ss κ ss . This is visually represented by summing each element of the matrix shown in Figure 4c.
We provide results for Si in Figure S5, where the lifetime for every mode qs is computed using third order IFCs and phono3py 27 . Our results agree with those of Simoncelli et al. 6 and we too conclude that thermal conductivity in Si is dominated by the phonon-gas channel, as expected. Additionally, we show that this dominance extends to high temperatures.
Due to its large unit cell size, third order IFCs for Yb 14 MnSb 11 could not be reasonably obtained. Therefore, two analytical forms were applied for the phonon-phonon scattering rate (inverse of phonon lifetime) Γ pp = τ −1 pp = AT and Γ pp = τ −1 pp = A ω 2 T . The coefficients A and A are fixed using thermal conductivity below 50K provided by Ribeiro et al. 30 (inset of Figure 4a). The ω-independent form ap-proximates phase space and anharmonicity contributions to the scattering rate as a constant (A) while maintaining the τ ∝ 1/T relationship which is expected for phonon-phonon scattering. The ω 2 form is more commonly used, since this is the frequency dependence that is commonly expected for low frequency acoustic phonons 31 . However, the many optical branches may invalidate the ω 2 dependence over most of the frequency range from phase space considerations. Regardless, the conclusions presented here were insensitive to the ω-dependence of τ pp , likely due to the close proximity of mode.
The treatment of higher order IFCs (3rd and 4th order) can influence the energy of phonon modes through renormalization and determines their linewidth, both of which may influence the conclusions regarding the magnitude of κ ph and κ diff 9 . The moderate Grüneisen parameter obtained here for Yb 14 MnSb 11 using the harmonic finite displacement method 6/13 ( Figure S8), and the insensitivity to the ω-dependence of τ pp suggest that these effect may be small in this case.
The thermal conductivity of Yb 14 MnSb 11 below 50 K 30 was used to fix the analytical form of the phonon scattering rate which is defined as the Matthiessen's rule combination of the phonon-phonon (Γ pp ) and crystal boundary scattering (Γ b,s (q) = |v s (q)|/d): Γ = Γ pp + Γ b . The simulation results shown in Figure 4a used A = 4 × 10 8 (sK) −1 and d = 0.8 µm when Γ pp , and A = 1.5 × 10 −17 s/K and d = 0.8 µm when Γ pp = A ω 2 T . An 8 x 8 x 8 q-mesh for Brillouin zone integration. We note that Yb 14 MnSb 11 is intrinsically p-type with a carrier concentrations in the range of 10 20 to 10 21 cm −3 (our sample was measured to have 2.9 ± 0.5 × 10 20 cm −3 ). While we do not explicitly consider electron-phonon scattering effects 32 , since we calibrate Γ to experimental data, they may influence the extracted coefficients, A and d. If the real linewidths are larger than we approximate, due to electron-phonon scattering effects, this is expected to suppress κ ph and increase κ diff .
After the analytical form of the phonon scattering rate was set, our two-channel lattice dynamics simulation captures the thermal conductivity above 300 K without any additional adjustable parameters, and indicates that κ diff is the dominant channel above room temperature. Since Isaeva et al. demonstrated that the thermal conductivity computed via twochannel lattice dynamics (i.e. quasi-harmonic Green-Kubo, QHGK) is influenced by the magnitude of the line width 4 , it was important to set phonon line widths using low temperature experimental data. The diffuson channel is expected to become important in systems with many atoms per cell, making methods for approximating the line width without the use of 3rd order IFCs increasingly important.
Histograms visualizing κ diff and κ ph are shown for 1000 K in Figure 4c. For the κ ph histogram, the large contribution from the three acoustic branches can be easily observed. The κ diff histogram reveals significant contribution across all branches, with the most heat carried by branches with s < 200 which correspond to the densely packed branches below 17.5 meV (see Figure 3d).
We highlight that the importance of κ diff in Yb 14 (Mn, Mg)Sb 11 comes mostly from crystal structure complexity and the resulting close energetic proximity of bands, rather than strong anharmonicity and large phonon line widths as is the case in CsPbBr 3 (Simoncelli et al. 6 ). This conclusion is supported by Grüneisen parameter computations (Eq. 8 in Computational Methods), where we obtained γ = 1.46 for Yb 14 MnSb 11 compared to γ = 2.6 for CsPbBr 3 . This indicates that while Yb 14 MnSb 11 is less anharmonic it is still predicted to have a large κ diff contribution. The decrease in overlap between Yb dominated modes and Sb-site dominated modes when the Sb mass is changed to that of As results in a 15 % decrease in κ diff at 850 K. The increase in overlap when increasing Sb-site mass to that of Bi results in a 15 % increase in κ diff . The trend in κ vib is predicted to follow that of κ diff .

Rational design of diffuson channel transport
The experimental and computational results presented above reveal that the energetic proximity of modes is an important metric governing the diffuson channel. Therefore, isolating a vibrational mode from others energetically seems to be a valid avenue for decreasing κ diff . To test this hypothesis we changed the mass of Sb to that of As, keeping the structure and IFCs unchanged. The results shown in Figure 5 indicate that the Sb-site dominated modes were shifted up in energy away from the Yb dominated modes (band structures shown in Figure S6) and a 15% drop in κ diff was indeed observed. Conversely, when the mass on the Sb-site was increased to that of Bi, the Sb-site dominated modes shifted down in energy increasing their overlap with Yb dominated modes and a corresponding 15% increase in κ diff was seen. Overall, when increasing the mass on the Sb-site from As to Bi the simula-tions showed a 35% increase of κ diff and a 27% reduction in κ ph . Nonetheless, an 18% increase in κ = κ diff + κ ph was observed, which is consistent with our hypothesis that modifying the energetic proximity of vibrational modes is a valid design strategy for amorphous-like thermal conductivity. The lattice thermal conductivity of Yb 14 MnBi 11 33 has been measured to be indeed larger than that of Yb 14 MnSb 11 11 , i.e. 0.35 W/mK versus 0.27 W/mK at 850 K, respectively.
We highlight that these findings are opposite of what is expected when only the phonon-gas channel is considered, where increasing atomic mass results in a lower frequencies, group velocities and lower thermal conductivity. Additionally, this trend is exactly opposite of what is expected for analytical (spectral) models of amorphous thermal conductivity, since these models have a direct correlation between κ vib and average vibrational frequency and/or speed of sound [21][22][23] . However, it is quite intuitive when the diffuson channel is considered in the mode-specific, two-channel lattice dynamical context ( Figure 1).
It is often desired to minimizing the total vibrational thermal conductivity κ vib . The approach to decreasing κ vib depends on whether κ ph or κ diff dominates. If the κ ph dominates, traditional defect scattering approaches can be used. If κ diff dominates one could focus on increasing the energetic spacing between the phonon modes. Unfortunately, mechanisms that decrease κ ph tend to increase κ diff . However, the two-channel lattice dynamics framework indicates that if the average energetic spacing between optical modes and their lifetime can be can be increased, while the group velocity and lifetime of acoustic modes are decreased, it may be possible to decrease both simultaneously.

Discussion
By examining the vibrational properties of the Yb 14 MSb 11 system in detail via computational and experimental methods, we show that for complex crystals one should not consider the phonon-gas model as incorrect, but rather as incomplete. The physical picture that emerges suggests that this general conclusion extends to defective and anharmonic crystals, as well as amorphous materials. To fully explain the vibrational thermal conductivity within a lattice dynamical framework, two conduction channels must be considered. The first is the standard phonon-gas channel (κ ph ), and the second we call the diffuson channel (κ diff ). By applying the two-channel lattice dynamical approach to Yb 14 (Mn, Mg)Sb 11 we find that the κ diff dominates above 300 K. Therefore, this material can be thought of as a crystal which conducts heat like an amorphous material above room temperature. Using Yb 14 (Mn, Mg)Sb 11 as a model system, two-channel lattice dynamics is used to demonstrated materials design strategies for controlling κ diff . We demonstrate how the κ diff can be decreased by up to 35 % by increasing the energetic spacing between vibrational modes which was achieved computationally by reducing the atomic mass on the Sb-site from Bi, to Sb, to As.

Synthesis
Cylindrical pellets of polycrystalline Yb 14 MgSb 11 and Yb 14 MnSb 11 with final dimensions of approximately 12 mm in diameter and 15 mm in height were synthesized as follows. Raw elements were weighed in stoichiometric amounts in an Ar filled glove box and sealed in a ball mill vial. High purity Mn (99.99 at.% plates Alfa Aesar) Granules Alfa Aesar), Mg( at. 99.99% filings Alfa Aesar), and Sb (at. 99.999% shot 5N plus) were used as received from suppliers. Yb(99.9 at.% ingot Alfa Aesar) was arc melted five times to further purify it of oxide prior to synthesis. These raw elements were ball milled together for 5 hours, re-mixing the constituents in the glove box every hour. The resulting powder was pressed at 900 C for 20 minutes under 45 MPa of pressure in a uniaxial hot press. The neutron diffraction histogram was refined to examine the phase purity of the samples. The diffraction signal (data points in Figure S2) was extracted from the elastic line (energy transfer between -1 and 1 meV) of time-of-flight neutron scattering intensity, with an incident neutron energy of 25 meV ( Figure S1). This diffraction signal contained contributions from the Al sample can, while in the inelastic neutron scattering intensities shown in Figures 3 and S1, the Al sample can contribution has been removed. The diffraction signal was refined in GSAS-II 34 using Yb 14 MSb 11 (I4 1 /acd no. 142 with c/a>1) and Al (Fm3m no. 225) phases. The diffraction from the empty Al can was used to set initial values for the instrument line broadening parameters, which were then refined along with the lattice parameters, atomic displacement parameters, and phase fractions of both phases. Therefore, these refinements are only used for phase analysis and no physical interpretation is given to the peak broadening. These refinements and phase analysis indicate phase purity of the samples in the sense that the vibrational density of states obtained is dominated by the Yb 14 MSb 11 phase.

Transport measurements
Thermal conductivity was calculated from the relation κ = DdC p , where D is the thermal diffusivity measured with a Netzsch LFA 457 laser flash apparatus, d is the geometrical density of the material, and C p is the heat capacity at constant pressure. d = 8.23g/cm 3 for Yb 14 MnSb 11 and d = 8.13g/cm 3 for Yb 14 MnSb 11 , which is 98% of the theoretical density for both, based on diffraction data. C p of the compounds were calculated via the Dulong Petit law, where C p = 3RN/M, where N is the number of atoms in a formula unit and M is the molar mass of the formula unit, and R is the molar gas constant. The measured κ tot can be obtained by adding κ e = Lσ T to the data in Figures 4a and b, where the polynomial expression for the measured σ is given in Figure S9.
The electrical conductivity (σ ) and Hall coefficient measurements were determined using the 4-point probe Van der Pauw technique with a 0.8 T magnetic field under high vacuum. 35 The electrical conductivity measurements are shown in Figure S9, along with polynomial fits which are used to compute the electronic contribution to the thermal conductivity, κ e = Lσ T , with L = 2.44 × 10 −8 WΩK −2 . This use of L is consistent with Toberer et al. 11 , but may be an over estimation in the Yb 14 (Mn, Mg)Sb 11 . The error bars on our data in Figure 4a and b embody laser flash measurement error, along with and estimation of the potential errors in subtracting κ e .

Inelastic neutron scattering
Inelastic neutron scattering experiments were conducted on the ARCS spectrometer at the Spallation Neutron Source (SNS), Oak Ridge National Laboratory, with incident neutron energies of 25 meV and 45 meV, and a sample temperature of 300 K. The samples were dense polycrystalline cylindrical pellets approximately 12 mm in diameter and 15 mm in height. All data are first represented as a function of neutron energy transfer, E, and momentum transfer, |Q| = 4π sin(θ )/λ , where θ is the scattering angle and λ is the neutron wavelength. The energy resolution is approximately 1.5 meV. An example of the scattering function (dynamic structure factor) S(|Q|, E) data obtained via the program Mslice in the Data Analysis and Visualization Environment (DAVE) 36 is shown in Figure 3. Due primarily to the large neutron absorption cross-section of Yb, neutron absorption corrections were included in the analysis as implemented by Mslice. This data was used to compute the vibrational density of state by first integrating S(|Q|, E) from |Q| = 4 to 8 Å −1 , obtaining S(E).
Data at the highest scattering angles, from 120 to 136 o were masked and data above 36 meV were removed to avoid any instrument related artifacts in S(|Q|, E). The background was subtracted from S(E) and the data was converted from arbitrary units to counts by recognizing that the intensity and error bars obtained from Mslice are An and A √ n, where n is counts and A is an arbitrary scaling factor. After these corrections, S(E) was converted to the vibrational density of state g(E) by removing the elastic scattering peak and multi-phonon contributions in the program isdos10 37 , similar to the approach used in Ref. 12 .

Computational methods
Interatomic force constants. The computation of harmonic interatomic force constants was performed using density functional theory (DFT) in the Vienna Ab initio Simulation Package (VASP) [38][39][40][41] and Phonopy 42 . First, the optimization of the electronic and ionic structures of Yb 14 MgSb 11 , Yb 14 MnSb 11 , Si, and CsPbBr 3 were performed with periodic DFT in VASP [38][39][40][41] with strict convergence criteria. For electronic and structural optimization, a criterion of ∆E < 10 −7 and ∆E < 10 −5 eV per cell was used, respectively. The projector augmented-wave method 43,44 with a plane-wave cutoff of 520 eV and the Perdew-Burke-Ernzerhof (PBE) functional were used during these calculations 45 . For Yb, we used Yb 2+ pseudo-potentials where the f electrons are considered to be in the core. Furthermore, the optimizations of the primitive unit cells were performed with a 4x4x4 k-space mesh (where k is the electron wave-vector) for Yb 14 MgSb 11 and Yb 14 MnSb 11 , a 17x17x17 k-space mesh for Si, and a 7x7x5 k-space mesh for CsPbBr 3 . Since Mn in this structure is expected to exhibit magnetic properties 46 , we utilized a ferromagnetic model for Yb 14 MnSb 11 during the optimization with magnetic moments on all Mn atoms in the primitive unit cell. We compare the results to an anti-ferromagnetic model in the Supporting Information ( Figure S7). The phonon frequencies are only influenced slightly by the magnetic model (root mean square deviation between the phonon frequencies is smaller than 0.12 meV on a 8x8x8 q-point grid).
Once the optimized ionic and electronic structures were obtained, the harmonic interatomic force constants (IFCs) were computed using the finite displacement method as implemented in Phonopy 42,47 , with a displacement of 0.01 Å and with the help of the conventional cell of Yb 14 MgSb 11 and Yb 14 MnSb 11 , a 3x3x2 supercell of the conventional cell of CsPbBr 3 , and a 4x4x4 supercell of the conventional cell of Si as the supercell. The forces for this evaluation were calculated at the Γ-point for Yb 14 MgSb 11 , Yb 14 MnSb 11 , and CsPbBr 3 , and with a 2x2x2 k-space mesh for Si. Again, a ferromagnetic model was used for Yb 14 MnSb 11 . Furthermore, a higher plane-wave cutoff of 900 eV was used to calculate the second-order force constants for Si.
To obtain third-order force constants for Si, Phono3py 27 in combination with VASP was used. This was done with the default displacement of 0.03 Å and a 2x2x2 supercell of the conventional cell. The forces were evaluated with the help of a 5x5x5 k-space mesh and a plane-wave cutoff of 520 eV.
To calculate the mode specific Grüneisen parameters (γ s (q) = −d ln ω s (q)/d lnV ) for Yb 14 MnSb 11 and CsPbBr 3 , two additional structural optimizations were performed at a constant volume of 0.997 3 ×V 0 and 1.003 3 ×V 0 with V 0 as the DFT ground-state volume, allowing the lattice parameters to relax. The computational details were the same as for the structural optimization of the DFT ground-state structures. Furthermore, harmonic phonon calculations were performed with the same computational details as for the DFT ground-state structure of Yb 14 MnSb 11 and CsPbBr 3 , respectively. Based on this set of phonon calculations of the DFT ground-state structure and the compressed and expanded structures, mode Grüneisen parameters were calculated as implemented in Phonopy. The mode specific Grüneisen parameters were used to compute the thermodynamic average Grüeisen parameter, which is defined as a heat capacity weighted average of γ s (q), We obtained γ = 1.46 for Yb 14 MnSb 11 and γ = 2.6 for CsPbBr 3 at 300 K. A plot of the mode specific Grüneisen parameters is shown in the Supporting Information Figure  S8. The mode specific Grüneisen parameters of CsPbBr 3 are significantly larger than those of Yb 14 MnSb 11 . Raw 9/13 data for the phonon calculations can be accessed on https: //doi.org/10.5281/zenodo.3951927. Band structure, DOS, and INS. These harmonic IFCs were implemented in a lattice dynamical approach to obtain the phonon band structure. For the computation of the vibrational density of states (DOS), an 8x8x8 uniform q-mesh and the tetrahedron method was used to approximate the integration over the Brillouin zone. Additionally, the DOS was convoluted with a 1.5 meV full width half max Gaussian function which is consistent with the instrument resolution of the inelastic neutron scattering experiments.
Inelastic neutron scattering (INS) spectra for Yb 14 MSb 11 were simulated using OCLIMAX 24 . DFT phonon frequencies and polarization vectors sampled on a 5x5x5 gamma-centered q-point mesh were used as input. Full scattering (coherent, incoherent, elastic, inelastic) for powder samples was simulated, using the ARCS experiment parameters such as incident energy, detector coverage (scattering angles), resolution, and temperature. Multi-phonon scattering contributions up to 10-phonon events were also calculated based on incoherent approximation and included in the total spectra.
Thermal conductivity. Here we define the matrices required to compute κ vib as defined in Eq. 1.
The velocity matrix is written as This is a complex matrix which obeys the relationship, v ss (q) = v s s (q) * . Its diagonal components (s = s ) have negligible imaginary contributions. Eq. 9 is equivalent to Eq. 8 in Simoncelli et al. 6 , by using Eq. 30 and Eq. 4 of Auerbach and Allen 48 . The gradient of the dynamical matrix was computed at each q-point, where the gradient is taken analytically as implemented by phonopy. The diagonal elements of v i ss (q) are simply the group velocity The matrix element can be written explicitly as where ε s (q) is the eigenvector of normal mode at q in branch s. The star indicates a complex conjugate. D(q) is the dynamical matrix at q defined as the following ε * s (q)D(q)ε s (q) = ω s (q) 2 δ s,s where ω s (q) 2 are the eigenvalues. The heat capacity matrix, in units of J m −3 K −1 is given as C ss (q) = 2 (ω s (q) + ω s (q)) 4V Nk B T 2 [ω s (q)n o s (q)(n o s (q) + 1) + ω s (q)n o s (q)(n o s (q) + 1)] where n o s (q) is the Bose-Einstein distribution function of mode s at temperature T , ω s (q) is the phonon angular frequency of a phonon at wave vector q in branch s, V is the volume of the unit cell, N is the number of q-points considered on a discrete, uniform q-mesh.
When third order IFCs are available Γ s can be computed from first principles. However, due to its large unit cell size, third order IFCs for Yb 14 MnSb 11 could not be reasonably obtained. Therefore, the simple analytical form Γ pp = τ −1 pp = AT was used for the phonon-phonon scattering rate (inverse of phonon lifetime). The coefficient A, which is meant to approximate anharmonicity and phase space effects, we set using low temperature thermal conductivity data. A Mattheissen's rule combination of the phonon-phonon (Γ pp ) and crystal boundary scattering was used (Γ b,s (q) = |v s (q)|/d): Γ = Γ pp + Γ b , and the coefficients A = 4 × 10 8 (sK) −1 , d = 0.8 µm were obtained Eqs. 13, 9, and 14 combine to give the thermal conductivity tensor (using index notation) κ i j = ∑ qss C ss (q)v i ss (q)v j s s (q)τ ss (q).
In matrix notation κ i j is represented as κ. By omitting the tensor superscripts in Eq. 1, we imply κ = κ xx .
To test our numerical implementation, the phonon-gas channel of Si computed via our code was compared to the thermal conductivity computed directly from phono3py and the results converged to the same value.
The phonon-gas and diffuson channel thermal conductiv- The xx components are reported in the Results section. Allen for helpful discussions surrounding the physical interpretation of simulated results and terminology. RH thanks Michele Simoncelli for correspondence regarding the computational details of the unified twochannel lattice dynamical approach.