Vibrational hierarchy leads to dual-phonon transport in low thermal conductivity crystals

Many low-thermal-conductivity (κL) crystals show intriguing temperature (T) dependence of κL: κL ∝ T−1 (crystal-like) at intermediate temperatures whereas weak T-dependence (glass-like) at high temperatures. It has been in debate whether thermal transport can still be described by phonons at the Ioffe-Regel limit. In this work, we propose that most phonons are still well defined for thermal transport, whereas they carry heat via dual channels: normal phonons described by the Boltzmann transport equation theory, and diffuson-like phonons described by the diffusion theory. Three physics-based criteria are incorporated into first-principles calculations to judge mode-by-mode between the two phonon channels. Case studies on La2Zr2O7 and Tl3VSe4 show that normal phonons dominate low temperatures while diffuson-like phonons dominate high temperatures. Our present dual-phonon theory enlightens the physics of hierarchical phonon transport as approaching the Ioffe-Regel limit and provides a numerical method that should be practically applicable to many materials with vibrational hierarchy. Predicting thermal transport in low-thermal-conductivity (κL) materials is challenging. Here, the authors propose a dual-phonon theory, where normal phonons are treated using the Boltzmann thermal equation and diffuson-like phonons are treated within diffusion theory, yielding robust predictions of κL.

L ow-thermal conductivity (κ L ) crystals are of great interest in a variety of applications including thermal barrier coatings (TBC), thermoelectrics and nuclear reactors. They often show intriguing thermal-transport properties: κ L decreases inversely with temperature (T) at intermediate temperatures as expected for crystals; but shows weak or even no distinct dependence on T at high temperatures, which is an anomalous, glass-like behavior [1][2][3][4][5][6][7] . While the former can be explained within the scheme of standard phonon Boltzmann transport equation (BTE) by primarily considering three-phonon scattering 8 , the latter is still an open question. Advance in developing unified theories and numerical methods is just looming over the horizon to improve understanding the disparate κ L~T relationship in these crystals.
Recent studies have identified that in low-κ L crystals, some phonon modes have mean free path (l) shorter than the Ioffe-Regel limit 9 , casting questions whether these modes can still be defined as phonons and how they contribute to thermal transport. Theoretical models of various sophistication have been developed to resolve this challenge [10][11][12][13][14][15][16] . Agne et al. 10 proposed that heat transfer in low-κ L complex crystals could be reasonably described by assuming a media of diffusons according to the random-walk-based diffusion theory, as opposed to phonons, and a model of diffuson-mediated κ L was proposed to better explain the experimental results particularly at high temperatures. The diffuson model by itself, however, does not fit to pure crystals where phonons are still well-defined as ensured by the periodicity of the lattice. On the other hand, other models have considered the hierarchy of vibrational modes based on the Ioffe-Regel limit. Chen et al. 11 proposed that for weakly disordered crystals with complex unit cell (e.g., higher manganese silicides), κ L could be explained by a hybrid phonon and diffuson model. The model employed a few fitting parameters and used the inelastic neutron scattering spectra to obtain an Ioffe-Regel crossover-frequency of 20 meV, below which the vibrational modes were treated as phonons and above which were diffusons. Mukhopadhyay et al. 12 used the scale of interatomic spacing as Ioffe-Regel criterion and proposed that for Tl 3 VSe 4 crystals, the phonon modes with mean free path l smaller than the Ioffe-Regel limit no longer behave as phonons, but should be replaced by hopping modes whose frequencies or eigenvectors cannot be meaningfully defined. They hence proposed a two-channel model that combines the phonon channel treated with BTE and the hopping channel calculated using Einstein's model or Cahill's model, yielding κ L and its temperature dependence in better agreement with experimental data. However, some open questions still remain. Particularly, how do the well-defined phonons interact with the hopping modes? How to subtract the well-defined modes from the hopping channel in κ L calculation? On the other hand, we note that infrared and Raman spectroscopies of some semiconductors have shown that, the zone-center optical phonons, many of which would have extremely short or nearly zero l, still have welldefined frequencies and linewidths (scattering rates) that can be accurately predicted by first-principles calculations 17,18 . This may suggest that the phonon concept for these modes is still valid, while the failure is with BTE which does not recognize the physical lower-limit of l. We can also note that this is not the first time BTE fails for phonons; in fact, BTE was known to fail for coherent phonons in superlattices or phononic crystals; therein, it does not capture the wave effects 19,20 .
Most recently, attempt has also been made to unify the thermal-transport theory in crystals and glasses. Simoncelli et al. 13 have transmuted the BTE formulism into a κ L equation written in terms of the phonon velocity operator with diagonal and off-diagonal elements describing the particle-like propagation of phonons and the wave-like tunneling of coherence, respectively. Applying this model in perovskite CsPbBr 3 arrives at reasonable simulation of its glass-like κ L . Meanwhile, Isaeva et al. 14 developed a quasi-harmonic Green-Kubo method, as if to generalize the Allen-Feldman model 21-24 for amorphous systems into relaxation-time-approximation-based BTE for crystals by expressing energy transport in a quantum mechanical fashion, and gives reasonable κ L predictions for amorphous and crystalline Si. Nevertheless, further physical insights are expected to improve understanding the nature of hierarchical vibrations in the context of physically based theories.
Inspired by the idea of vibrational hierarchy from previous models, while attempting to resolve the open questions, in this work we sparkle a different concept that the vibrational modes with very short l could still be treated as phonons with welldefined frequencies, eigenvectors, and scattering rates, but their heat conduction should be described by the diffusion theory instead of BTE. We hence propose a dual-phonon theory, by treating the short-l phonons with the diffusion theory, and other normal phonons with BTE. Our theory does not introduce a different type of heat-carrying modes other than phonons, and eliminates the theoretical challenge of how they would interact with phonons and alter their scattering events. Also, for a sophisticated model, we introduce three different criteria, based on phonon mean free path, wavelength, and thermal diffusivity, to judge mode-by-mode between the normal phonon channel and diffuson-like phonon channel, and hence avoid the doublecounting issue. The three criteria all yield consistent results that agree quantitatively with experiments, demonstrating the robustness and predictive capability of our theory. Our theory is demonstrated on La 2 Zr 2 O 7 , a thermal-barrier-coating (TBC) candidate material for gas turbine technologies 25 , and Tl 3 VSe 4 , a potential thermoelectric material 12 . More background of La 2 Zr 2 O 7 could be found in Supplementary Note 1. Our approach is able to explain the thermal conductivity and the intriguing κ LT dependence over the entire temperature range, and is expected to help understand the thermal transport of such low-κ L crystals. Also, our model may provide physical insights toward unifying the theories of thermal conductivity in crystals and amorphous materials.

Results
Normal phonons and diffuson-like phonons. We first calculate the phonon properties of La 2 Zr 2 O 7 using the standard anharmonic lattice dynamics based on density functional theory. Immediately we find that a large percentage of vibrational modes have very small l even at room temperature, and more so at high temperatures, as shown in Fig. 1a. Similar features have been identified for Tl 3 VSe 4 (see Supplementary Fig. 5). Apparently, these modes cannot be treated as normal phonons, which by definition, are expected to propagate far enough to sample the periodicity of the transport media, i.e., comparable to the scale of phonon wavelength or several lattice spacing [21][22][23] . In developing our theory, these modes are treated as diffuson-like phonons.
In order to sophisticatedly judge whether a phonon mode should be treated as normal phonon or diffuson-like phonon, we propose three criteria, as conceptually shown in Fig. 2. The first two criteria are derived from the Ioffe-Regel limit of phonons 9 , arguing that the value of l for normal phonons should not be smaller than their wavelength (λ) (I. l-λ criterion), or the minimum interatomic spacing (a min ) of the lattice (II. l-a min criterion), respectively. The third criterion argues that a vibrational mode should be characterized as diffuson-like if its phonon thermal diffusivity (D Phon ) becomes smaller than its diffuson thermal diffusivity (D Diff ) (III. D Phon -D Diff criterion).
Here we note that, another reasonable criterion to define a normal phonon mode is based on its long-enough relaxation time (τ) as τ >> 1/ν (ν is the vibrational frequency). For acoustic phonons approaching the long-wavelength-limit, this criterion is the same with the l-λ criterion, whereas for some phonons whose phase velocity (v p ) is considerably higher than group velocity (v g ), it is less restrictive 23 . Therefore, the τ-criterion is not considered in developing our theory.
For modeling κ L in materials with such strong vibrational hierarchy, we assume a thermal-transport media with both normal phonons and diffuson-like phonons. All these phonon modes are able to transfer heat (i.e., non-localized), but following BTE theory for the former and diffusion theory for the latter, respectively. Our dual-phonon theory combines a normal phonon channel and a diffuson-like phonon channel: Heat conduction from both channels are derived by following the physical picture that heat is carried by atomic vibrations of a solid. There are three components to be considered: (1) the number of vibrational modes that are available to carry heat, summed up to be the total degrees of vibrational freedom; (2) the amount of heat that can be carried by each vibrational mode; (3) the propagation behavior of vibrations through the dual-phonon media. We write the contribution from each channel as: Here, i and j denote the indices of vibrational modes, sampled over the Brillouin zone (BZ). N Phon and N Diff are the numbers of normal phonons and diffuson-like phonons, respectively. One of the merits of our theory is that, we impose the conservation of vibrational degrees of freedom by requiring N Phon + N Diff = 3N, where N is the total number of atoms. As such, once a mode is characterized as a normal phonon, it cannot be treated as a diffuson-like phonon at the same time, and vice versa. C s is the permode specific heat following the Bose-Einstein statistics of phonons. D Phon and D Diff are defined as the per-mode thermal diffusivities for normal phonons and diffuson-like phonons, respectively. In our present model, D Phon is calculated using the phonon BTE theory; whereas D Diff could be calculated based on the random-walk theory 10 (denoted as D RW Diff ) or the Allen-Feldman theory 23 (denoted as D AF Diff ). Details are presented in the Methods section. In the present work, we only performed limited test for our dual-phonon theory coupled with D AF Diff (using the III. D Phon -D Diff criterion), mainly due to its high computational expense.
Hierarchy of lattice vibration in La 2 Zr 2 O 7 . The hierarchy of lattice vibrations for La 2 Zr 2 O 7 is shown in Fig. 1. We see that the vibrational modes that are judged as diffuson-like phonons account for 75.35%, 23.82%, 65.42%, and 69.73% of the total number of modes at T = 300 K, respectively, according to our proposed judging criteria I, II, and III coupled with the randomwalk theory, and III coupled with the Allen-Feldman theory. The fractions rise to 98.00%, 80.23%, 91.25%, and 93.55% at T = 1500 K. As illustrated in Fig. 3, this behavior originates either from small v g , especially for high-frequency vibrations, or from small τ, indicative of intense scattering among those modes. The increased population of diffuson-like phonons at higher temperatures is presumably due to increased phonon scattering.
To better visualize the degree of the vibrational hierarchy, the vibrational density of states (DOS) for normal phonons vs. diffuson-like phonons are illustrated in Fig. 1b, d, f, h. There is a crossover from normal phonon-dominated states at the lowfrequency range to diffuson-like phonon-dominated states at highν range, and the crossover shifts to lower-ν at higher temperatures. Take the case of II. l-a min Criterion as an example. At T = 300 K, normal phonons dominate below ν~13 THz; whereas for T = 1500 K, diffuson-like modes exhibit the first peak at ν~2 THz, and quickly becomes dominant at above ν~5 THz. This could be understood from increased scattering at higher temperatures, resulting in suppressed τ and l values for all vibrational modes. Furthermore, detailed comparisons among the three criteria show that, the number of vibrational modes assigned as diffuson-like phonons based on the random-walk theory increases in the order of II. l-a min Criterion, III. D Phon -D Diff Criterion, and I. l-λ Criterion, and the major difference comes from assignment of the vibrational modes in the frequency range 2 < ν < 13 THz. On the one hand, the fact that a larger number of vibrations are assigned as diffuson-like phonons by Criterion III than Criterion II could be understood from the interplay of low ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16371-w v g · l and high ν for certain modes. On the other hand, the difference between Criterion I and Criterion II is expected to be more pronounced for zone-center long-wavelength vibrations, as certain modes have a min < l < λ hence they are assigned as normal phonons by the l-a min Criterion while as diffuson-like phonons by the l-λ Criterion. Moreover, as shown in Fig. 1f, the ratio N Diff /N Phon appears to increase with the vibrational frequency under Criterion III coupled with the random-walk theory, but the trend is not so clear for the results derived from the Allen-Feldman theory (Fig. 1h). Such difference could be understood from the different frequency dependence of D RW Diff vs. D AF Diff . The random-walk theory assumes that heat is  Lattice thermal conductivity. The κ L values calculated from an iterative solution of BTE coupled with anharmonic lattice dynamics are shown in Fig. 4a. Both phonon-phonon scattering and isotope scattering are considered in our calculations, but the latter has little impact on κ L values and their temperature dependence. At room temperature, the calculated κ L shows reasonable agreement with experimental data [26][27][28] . At higher temperatures, it yields κ L~T −1 dependence as expected from the dominance of phonon-phonon Umklapp scattering, whereas experimental data goes nearly temperature-independent. The calculated κ L = 0.55 W·m −1 ·K −1 at T = 1500 K is only one-third of the measured value κ L = 1.5 W·m −1 ·K −1 , and is even lower than the reported high-temperature limit (κ min = 1.2 W·m −1 ·K −1 ) for La 2 Zr 2 O 7 crystals 29,30 . Moreover, we notice that including the state-of-the-art four-phonon scattering 31-33 is expected to give further reduction of κ L at high temperature and results in stronger-than-T −1 temperature dependence, and thus could not bridge the gap between experimental and BTE-derived results in our case. Guided by the Ioffe-Regel limit of phonons 9 , we recalculate κ L within the framework of phonon BTE theory, by assuming a lower-bound of l for all modes, i.e., the l values are manually set to a min if l < a min . This trial yields κ L = 0.61 W·m −1 ·K −1 at T = 1500 K, showing inadequate correction to the previous results. Such significant gap between the calculated temperature-dependent κ L and experimental values indicate that the conventional BTE theory becomes invalid for heat conduction in La 2 Zr 2 O 7 , and a different physical picture is needed to describe the transport behavior of the interesting small-l or small-D Phon phonon modes.
Next, we calculate κ L of La 2 Zr 2 O 7 using the two-channel model proposed by Mukhopadhyay et al. 12 , by combining a phonon conduction channel (κ phonon ) and a hopping channel (κ hop ) from modes of l < a min . In this model, κ phonon is calculated from phonon BTE theory, by excluding the contribution from small-l vibrations; whereas κ hop is calculated either via Cahill's formula 34,35 or Einstein's formula from which the Cahill's model is derived. Cahill's formula requires v g of acoustic phonons as input parameters, which could be extracted from phonon dispersions of La 2 Zr 2 O 7 (shown in Supplementary Fig. 2): 3837 and 3969 m·s −1 for transverse acoustic phonons (TA1 and TA2), and 6872 m·s −1 for longitudinal acoustic phonons (LA). Einstein's formula requires defining the Einstein temperature (θ E ) for the oscillators, which could be estimated from the calculated lowtemperature specific heat of La 2 Zr 2 O 7 (shown in Supplementary  Fig. 3): θ E = 180 K, in moderate agreement with reported values 36 . Using these parameters as input, the total κ L (κ phonon + κ hop ) are calculated to be 1.90 and 0.92 W·m −1 ·K −1 at T = 1500 K by employing Cahill's formula and Einstein's formula, respectively. The corresponding temperature dependencies are κ L~T −0.31 and κ L~T −0.65 above T = 1000 K. Now the experimental values of κ L fall within the range provided by the twochannel model, and the temperature dependence shows improved agreement. As shown in Fig. 4a, the range in between might be originated from uncertainties in defining the θ E values, which has also been pointed out in ref. 12 . Besides, we note that in order to preserve the total number of vibrational modes, the normal phonon modes need to be subtracted from the hopping channel, which is difficult to do when incorporating the Cahill's formula in the two-channel model.
The results of our dual-phonon theory for La 2 Zr 2 O 7 are shown in Fig. 4a. Our model coupled with the random-walk theory gives weakened κ L~T dependence as T increases: κ L~T relationship calculated using I. l − λ criterion (II. l-a min criterion; III. D Phon -D Diff criterion) weakens from κ L~T −0.69 (κ L~T −0.72 ; κ L~T −0.67 ) for 300 K < T < 500 K, to κ L~T −0.47 (κ L~T −0.40 ; κ L~T −0.47 ) for 500 K < T < 1000 K, and to κ L~T -0.29 (κ L~T −0.24 ; κ L~T −0.30 ) for T > 1000 K, reasonably reproducing the flattening-out behaviors of the experimental κ L~T data. Partial contributions from normal phonons (κ Phon L ) and diffuson-like phonons (κ Diff L ) to the total κ L are plotted in Fig. 4b. It is shown that, κ Diff L starts to dominate over κ Phon L at around T = 700, 1000, and 800 K based on Criterion I, II, and III, respectively. Clearly, the role of diffuson-like phonons for the heat conduction of La 2 Zr 2 O 7 is more significant at high-temperature ranges. Besides, the calculated κ Phon L decreases in the order of criteria II, III, and then I, while the trend is reversed for κ Diff L . This result is consistent with the percentage of diffuson-like phonons assigned out of all the vibrational modes (inset Fig. 4b). Evidently, despite using three different criteria to judge normal phonons vs. diffuson-like phonons, the calculated κ L values as 5 I. l -criterion, RW Expt. 26 Expt. 27 Expt. 28 Iterative ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16371-w well as the κ L~T dependence show quite similar features, demonstrating the robustness of our approach. Moreover, our dual-phonon theory is also validated on Tl 3 VSe 4 , a potential thermoelectric material with ultra-low κ L that was used as the model materials in ref. 12  Also shown in Fig. 4 are the results from our dual-phonon theory coupled with the Allen-Feldman model. The percentage of diffuson-like phonons assigned out of all the vibrational modes using the III. D Phon -D Diff criterion coupled with D AF Diff is generally consistent with the results derived with D RW Diff , and the calculated values of κ Phon L from the two routes show minor differences. Moreover, the κ L~T dependence calculated with D AF Diff gives κ LT −0.61 for 300 K < T < 500 K, κ L~T −0.42 for 500 K < T < 1000 K, and κ L~T −0.27 for T > 1000 K, consistent with that derived using D RW Diff . These results demonstrate the robustness of the dualphonon theory. Evidently, imposing the dual-phonon theory coupled with Allen-Feldman model tends to yield higher κ Diff L , and thus slightly higher values of total κ L , as compared with those derived with the random-walk picture. This is actually inherited from the different frequency-dependent formulas of D RW Diff vs. D AF Diff as we have discussed above. Interestingly, at least for La 2 Zr 2 O 7 , the random-walk-based thermal diffusivities in our dual-phonon theory perform better than the Allen-Feldman theory. Meanwhile, this route requires far less computational expense, which is potentially favorable for high-throughput κ L predictions. Nevertheless, more works are required in the future to test this issue using a larger material pool.

Discussion
We can see that our dual-phonon theory emphasizes the physics of hierarchical phonon transport in low-κ L materials, by treating all vibrational modes within the phonon picture yet with different thermal-transport behaviors, i.e., normal phonons with the BTE theory and diffuson-like phonons with the diffusion theory. Basically, our dual-phonon theory proposes a conceptual change that the vibrational modes are still phonons upon approaching the Ioffe-Regel limit, and the phonon frequency, eigenvector, and relaxation time could be rigorously described by first principles. Specifically, the identified diffuson-like phonons with small l or small D Phon still fall in the physical picture of phonons with welldefined scattering rates that can be reliably predicted from first principles; just that their heat conduction cannot be treated using the scheme of mean free path or BTE. In this way, the scattering and thermal transport of diffuson-like phonons are decoupled. The three proposed physics-based judging criteria, without relying on any fitting parameters, enable a per-mode-based judging between normal phonon and diffuson-like phonon channels, and give consistent κ L predictions in agreement with experimental data, demonstrating the robustness of our approach. Moreover, while this paper was under review, we learned that progressive efforts have been made to unify thermal transport in crystals and amorphous materials, which involve sophisticated analytical and mathematical derivations, yet the corresponding physical insights are still under development 13,14 . In parallel with these efforts, our dual-phonon theory provides physical insights in understanding the vibrational hierarchy of crystals having low and glass-like κ L , and leads up to a κ L -prediction method with good robustness and favorable computational expense.   It would be interesting to explore future opportunities to expand and advance our proposed dual-phonon theory. First, the current study does not include the effects of four-phonon scattering or phonon renormalization, which are important effects addressed elsewhere but are beyond the scope of the present study, partly due to the very high computational cost for La 2 Zr 2 O 7 (with complex and large unit cell). Four-phonon scattering is important at high temperature for all materials and even at room temperature for strongly anharmonic materials and materials with exceptionally weak three-phonon scattering [31][32][33] . Phonon renormalization is important for strongly anharmonic materials, especially those exhibiting dynamic instability [37][38][39][40] . Indeed, consideration of four-phonon effect and the temperatureeffect might arrive at more accurate simulation of the phonon group velocities and scattering, leading to better placement of normal phonons vs. diffuson-like phonons, and more accurate κ L and its temperature dependence.
Second, the primary scope of our dual-phonon theory is to differentiate the thermal-transport behaviors of normal phonons against diffuson-like phonons in crystals based on whether or not they are ill-defined in the space scale (in terms of small l or D Phon ), given that they fall within the phonon picture in terms of scattering. Interestingly, we find that some vibrational modes in La 2 Zr 2 O 7 and Tl 3 VSe 4 exhibit 1/τ > ν; in other words, their lifetime is too short and thus should be considered ill-defined in the time scale. Similar behavior has been universally observed in crystals having certain types of perturbation or disorder 12,13,[41][42][43][44] . As shown in Fig. 5a, c, the number of modes that are ill-defined in the time scale (1/τ > ν) account for 0.21% for La 2 Zr 2 O 7 and 11.58% for Tl 3 VSe 4 at T = 300 K, which are much smaller than the percentages of the ill-defined-in-space phonons. They rise to 45.36% at T = 1500 K for La 2 Zr 2 O 7 and 59.71% at T = 500 K for Tl 3 VSe 4 , respectively, though. We then further evaluate the effects of these modes on the scattering rates of phonons that are welldefined in the time scale (1/τ < ν), by calculating the percentages of scattering rates of well-defined-in-time modes contributed by three-phonon scattering processes involving at least one illdefined-in-time mode. Figure 5b, d shows that the percentages are mostly below 20% for La 2 Zr 2 O 7 and below 60% for Tl 3 VSe 4 at T = 300 K, which rise to mostly above 70% at T = 1500 K for La 2 Zr 2 O 7 , and above 90% at T = 500 K for Tl 3 VSe 4 . Therefore, although our approach makes a step closer to a sound theory and appears to agree with experiments, whether we can treat scattering by these ill-defined-in-time phonons with standard anharmonic lattice dynamics still remains a challenging open question, especially at high temperatures. Future studies are warranted on how to better understand the interactions between the well-defined-in-time phonons and ill-defined-in-time phonons in our model and other lately developed models 13,14 .
The present results have important implications for a wide variety of low-κ L materials, including thermal barrier coatings (TBC), thermoelectrics, and nuclear materials. For these low-κ L materials, there might be large number of diffuson-like phonons characterized by small l or small D Phon , and the significant contribution from κ Diff L to total κ L manifest with increased temperature, and eventually dominate over the contribution from κ Phon L when approaching the high-temperature limit. In the context of our analysis on the TBC material La 2 Zr 2 O 7 , the diffusonlike phonons mainly stem from low v g and/or high scattering rates, which could be further linked with the complexity of crystal structure and the heterogeneity of interatomic bonds. Such structural characteristics result in folding-in of phonon dispersion, suppression of phonon frequencies and thus group velocities, and serious tangling and scatterings among low-frequency acoustic and optical phonons (see Supplementary Note 4 and Supplementary Fig. 8 for details). This mechanism is expected to be a signature for many TBC candidates, i.e., oxide ceramics with complex crystal structure and vibrational hierarchy, probably including other RE 2 Zr 2 O 7 pyrochlores 25 , silicates 45,46 , and acuminate-silicates 47 etc. In this sense, our proposed dualphonon model is expected to inspire deeper understanding and practical calculation methodology of κ L for TBC materials (see Supplementary Note 5 for details).
In summary, a dual-phonon theory is proposed for the κ L of crystals with vibrational hierarchy, by considering normal phonons described by the BTE theory and diffuson-like phonons described by the diffusion theory. Three physics-based criteria are used to judge mode-by-mode between normal phonons and diffuson-like phonons. Applying this theory on La 2 Zr 2 O 7 and Tl 3 VSe 4 shows that κ L is mainly contributed by normal phonons at low temperatures, whereas by diffuson-like phonons at high temperatures. This theory successfully predicts the flattening-out of κ L~T trend upon temperature increment, in much better agreement with experiments than the conventional BTE theory. Meanwhile, it resolves the limitations of other existing models and leads to a computational procedure showing promising applicability. The improvement of heat conduction theory in lowκ L crystals will provide important insights in the development of TBC materials, thermoelectric materials, and nuclear materials. Future studies can include a comparison of our approach with those in refs. 13,14 , a consideration of how to combine the mathematical rigor of those approaches and the physical insights from our approach and other heuristic models, and a better understanding of the interactions between the well-defined-intime phonons and ill-defined-in-time phonons.

Methods
Calculation details of the dual-phonon theory. In this study, the per-mode specific heat is calculated by following the Bose-Einstein statistics for phonon: where k B is the Boltzmann constant; ћ is the reduced Planck constant; T is the absolute temperature; ω is the angular frequency of a vibrational mode; V is the volume of simulated unit cell. D Phon is calculated using the phonon BTE theory, by incorporating the phonon group velocity (v g ), mean free path (l), and relaxation time (τ) resulted from the scattering process. D Diff could be calculated by following the random-walk scheme proposed in Agne et al.'s model 10 , which assumes that heat is transferred by discrete jumps between independent harmonic oscillators. Its linear dependence with respect to ω could be understood from the diffusive picture, as if oscillators jumping at more steps per time-interval contribute to higher rates of energy transfer. Basically, the diffusion term D Diff has the units of D Phon (m 2 ·s −1 ), while avoiding independent definition of v g and l in the diffuson picture.
where n is the number density of atoms of the unit cell. In this study, moderesolved vibrational parameters (v g , l, λ, ω, τ, etc.) are gathered from outputs of DFT-based harmonic and anharmonic lattice dynamics calculations. Scattering information (τ and l) is gathered from iterative solution of the BTE; however, we tested that using parameters derived under the relaxation-time approximation (RTA) would have negligible influence on the present results. For comparison, D Diff could also be calculated using the Allen-Feldman formula 23 , which is grounded in the Green-Kubo theory and is at the reach of ab initio simulations.
where S jp is the heat-current operator measuring the thermal coupling between vibrational mode j and p based on their frequencies and spatial overlap of eigenvectors, and could be calculated from harmonic lattice dynamics. δ is the Dirac Delta function that could be approximated using Lorentzian broadening of width greater than the average mode frequency interval (Δ avg ).
where the wave vector K is summed over the Brillouin zone; e is the corresponding eigenvector; α and β denote the Cartesian directions. R s is the distance between each unit cell (labeled s) and the basis unit cell (labeled 0) within a periodic supercell system; and R kk′ is the distance between atom k and atom k′ within a unit cell. D k 0 k βα 0; s ð Þ is the dynamical matrix element, derived from the second-order interatomic force constant (Φ k 0 k βα 0; s ð Þ).
where m k and m k′ are the masses of the atom k and k′. The harmonic and anharmonic interatomic force constants (IFCs) are calculated via the real-space finite displacement difference method, where 2 × 2 × 2 supercells containing 88 atoms are constructed, and Monkhorst-Pack k-mesh are set as 3 × 3 × 3. The phonon frequencies and eigenvectors are obtained using the Phonopy 53 package, by diagonalizing the dynamical matrix constructed from the harmonic IFC matrices, and sampling on a 21 × 21 × 21 q-mesh. These are typical settings for the lattice dynamics calculations of rare-earth pyrochlore systems 25 . The anharmonic IFCs are obtained using the thirdorder scripts 54 . Interatomic interactions up to the twelfth nearest neighbors (12th NN) are taken into account, corresponding to a cutoff radius (r cutoff ) of 7.72 Å; whereas interactions beyond this range are taken to be zero. In fact, increasing the r cutoff from the 5th NN (corresponding to r cutoff = 5.19 Å) to the 12th NN could sufficiently converge the room-temperature κ L within 7%. It is noteworthy that La 2 Zr 2 O 7 is a complex oxide ceramic with long-range interatomic interactions 55 , which might lead to strong dependence of anharmonic IFCs and κ L on the number of neighbor shells. Qin et al. 56 reported that for such material systems, the extent to an adequate inclusion of long-range effect could be estimated by looking into how the interaction strength changes with increased distance between an atomic pair, based on analyzing the root mean square (RMS) of the elements of the harmonic IFC tensor (Frobenius norm): RMS ϕ mn À Á ¼ 1 9 where ϕ mn is the harmonic IFC between atom m and n; and ϕ αβ mn is the harmonic response of the force for atom m on the α-direction resulted from the displacement of atom n on the β-direction. Following this approach, the RMS(ϕ mn ) for all atomic pairs are analyzed as a function of the interatomic distance. It shows that for La 2 Zr 2 O 7 , setting the truncation at the 12th NN is expected to include most strong interatomic interactions up to RMS(ϕ mn ) = 2, and beyond this range it decays below RMS(ϕ mn ) < 0.15. For these reasons, we chose to include up to the 12th NN interatomic interactions in the present calculations to pursue higher precision.
The lattice thermal conductivity (κ Phon L ) is calculated by the iterative solution of the BTE as implemented in the ShengBTE 54 package, with integrations using a 16 × 16 × 16 q-mesh. The convergence of κ L with respect to the size of q-mesh is tested, and the results show that increasing the q-mesh up to 25 × 25 × 25 yields less than 2% difference to the room-temperature κ L of La 2 Zr 2 O 7 . Furthermore, non-analytical corrections are applied to the dynamical matrix to take into account long-range electrostatic interactions, based on calculations of Born effective charges (Z*) and dielectric constants (ε) via density functional perturbation theory (DFPT) 57 .
The Allen-Feldman calculations are performed with the same 2 × 2 × 2 supercell and 16 × 16 × 16 q-mesh, to ensure a consistent comparison between D Phon vs. D AF Diff on a per-mode basis. The Delta function in Eq. (8) is broadened into the Lorentzian form η=π ω j ð ÞÀω p ð Þ ð Þ 2 þη 2 . After convergence tests (details in Supplementary Note 7 and Supplementary Fig. 9), the Lorentzian broadening factor η is set to be 3.3Δ avg in our study, where Δ avg is the average mode frequency interval (Δ avg = 0.35 THz for La 2 Zr 2 O 7 ). Due to the high computational expense, we use per-mode thermal diffusivities on 1000 out of the 4096 wave vectors for our dual-phonon theory under the Criterion III coupled with Allen-Feldman formula. A convergence test shows that including up to 2000 wave vectors causes only marginal difference (<1%) on the calculated fraction of diffuson-like phonons (N Diff /N Total ) and the final κ L values.

Data availability
The source data of Figs. 1, 3-5, and Supplementary Figs. 1-9 are provided as a Source data file at https://archive.materialscloud.org/2020.0036/v1, and are further available from the corresponding author upon reasonable request.