S1 A. Details of computational procedure

# Quartic Anharmonicity of Rattlers and its Effect on Lattice Thermal Conductivity of Clathrates from First Principles

## Abstract

We investigate the role of the quartic anharmonicity in lattice dynamics and thermal transport of type-I clathrate BaGaGe based on ab initio self-consistent phonon calculations. We show that the strong quartic anharmonicity of rattling guest atoms causes the hardening of vibrational frequencies of low-lying optical modes and thereby affects calculated lattice thermal conductivities significantly, resulting in an improved agreement with experimental results including the deviation from at high temperature. Moreover, our static simulations with various different cell volumes shows a transition from crystal-like to glasslike around 20 K. Our analyses suggest that the resonance dip of observed in clathrates with large guest-free-spaces is attributed mainly to the strong Umklapp scattering of acoustic modes along with the presence of higher-frequency dispersive optical modes.

###### pacs:
63.20.dk, 63.20.kg, 82.75.-z

Intermetallic clathrates are a class of materials possessing a cage-like structure incorporating a guest atom, which is loosely bound inside an oversized cage and makes a large amplitude thermal vibration called “rattling”. Experimental and theoretical studies have evidenced that rattling guest atoms cause characteristic lattice dynamics of clathrates including low-frequency vibrational modes showing significant temperature () dependence Takasu et al. (2006); Christensen et al. (2008); Takabatake et al. (2014); Wu et al. (2016) and very low lattice thermal conductivity (LTC, Dong et al. (2001); Takabatake et al. (2014); Tadano et al. (2015). LTCs of clathrates are not only very low but also show exceptional and diverse -dependence. For example, of electron-doped type-I clathrate BaGaGe (BGG) shows a peak around 20 K followed by a decreasing region at higher temperature Sales et al. (2001); Avila et al. (2006), just like many crystalline solids. However, in the temperature region above 100 K, the -dependence is much milder than of typical crystalline materials May et al. (2009). More exceptional -dependence has been reported for type-I clathrates XGaGe (X=Sr, Eu) Nolas et al. (1998); Sales et al. (2001) and BaGaSn Avila et al. (2008). In these materials, the LTCs behave like a typical glass, showing a plateau region or a “resonance dip” Nolas et al. (1998) near 20 K, an increasing trend in 20–100 K, and a nearly -independent region above 100 K. These unconventional and diverse thermal transport in intermetallic clathrates has been attributed to the difference in the guest-free-space of the rattling atoms Sales et al. (2001); Suekuni et al. (2007, 2010), the host-guest coupling strength Bridges and Downward (2004); Xi et al. (2017), and the magnitude of static/dynamical disorders Bridges and Downward (2004); Christensen et al. (2016). Despite these constant efforts, the origin of the unusual LTCs of clathrates still remains unclear.

Recently, ab initio calculation of LTC based on the Peierls-Boltzmann theory Peierls (1955) has established itself as a convenient way to predict/analyze thermal transport phenomena in solids Lindsay (2016). Although the validity of the Boltzmann theory is limited to the cases where the phonon quasiparticle picture is well established Allen and Feldman (1993), it has reproduced experimental LTC of BGG in a relatively low temperature region Tadano et al. (2015). However, the ab initio Boltzmann approach considerably underestimated above K. A similar underestimation has also been reported in a more recent study Lory et al. (2017). These results clearly indicate the necessity of an improved theoretical approach. One of the most problematic approximations made in the conventional Boltzmann approach is the omission of the quartic anharmonicity. Indeed, the atomic displacement factor of guest atoms is so significant that the quartic anharmonicity cannot be neglected anymore in clathrate. Therefore, including the effect of the quartic anharmonicity in theory is at least essential to develop robust understandings of thermal transport in clathrates.

Here, we report volume- and temperature-dependence of of type-I clathrate BGG obtained from first-principles calculation, where the temperature renormalization of the vibrational frequency by the quartic anharmonicity is considered using the self-consistent phonon (SCP) theory Souvatzis et al. (2008); Errea et al. (2014); Tadano and Tsuneyuki (2015); van Roekeghem et al. (2016). We show below that the strong quartic anharmonicity of rattling motions makes the deviation from at high temperature. Our computational analyses also indicate that the resonance dip observed in clathrates with large cage-sizes can be attributed to the strong Umklapp scattering of acoustic modes along with the presense of high-frequency dispersive optical modes and the phonon-boundary scatterings.

We start by giving a brief overview of the SCP method. Following our recent implementation Tadano and Tsuneyuki (2015), we calculate the -dependent anharmonic frequencies and polarization vectors by diagonalizing the matrix defined as

 Extra close brace or missing open brace (1)

Here, is the harmonic phonon frequency with the crystal momentum and the branch index , and is the shorthand notation for satisfying and . is the reciprocal representation of the fourth-order interatomic force constants (IFCs). is the mean square displacement of the normal coordinate , where is the Bose-Einstein distribution function, with the Boltzmann constant , and is the reduced Planck constant, respectively. Since the right-hand-side of the equation depends on the solution , Eq. (1) needs to be solved self-consistently. More details about the present SCP formalism including the expression of the coupling constant and the efficient numerical implementation are described elsewhere Tadano and Tsuneyuki (2015, ).

All of the DFT calculations were conducted using the vasp code Kresse and Furthmüller (1996), with the GGA-PBE functional Perdew et al. (1996) and the projector augmented wave (PAW) method Blöchl (1994); Kresse and Joubert (1999), and the phonon calculations were performed using the alamode package Tadano (); Tadano et al. (2014). An ordered unit cell containing 54 atoms (space group ) was employed for modeling lattice anharmonicity and phonon thermal transport. To investigate the effect of the guest free space on vibrational and transport properties, the lattice constant of BGG was compressed/expanded from the optimized value (10.954 Å, hereafter called “opt.”) by % to % in steps of 2%. Harmonic and anharmonic IFCs necessary for the present SCP and LTC calculations were estimated using the compressive sensing lattice dynamics method Zhou et al. (2014). More detailed computational procedures are provided as the Supplemental Information (SI) sup ().

Calculated harmonic phonon dispersion curves of BGG in the low-frequency region ( cm) are shown in Fig. 1(a) by dotted lines. With increasing the lattice constant, the frequencies of the low-lying optical modes associated with rattlers decrease. In the largest-volume case ( Å), twelve phonon modes become unstable , all of which can be well characterized as collective motions of Ba(2) atoms inside the tetrakaidecahedral cages (see Fig. 1(c)). This indicates the transition of the potential energy surface (PES) minima of the Ba(2) atoms from the on-center site to off-center sites. To see the volume dependence of the lattice anharmonicity more directly, we also calculated the PES of the Raman active guest mode by displacing atoms by with being the amplitude of the normal coordinate of the mode. The actual displacement pattern is schematically shown in Fig. 1(c) by arrows. In Fig. 1(b), we compare the PESs calculated by DFT (open circles) and those calculated from the IFCs (lines). The left panel of Fig. 1(b) shows that the harmonic approximation fails to capture the actual shape of the PES, indicating the significance of the anharmonicity. Indeed, if we include the contribution from the fourth-order IFCs, we obtain overall good agreements with the DFT results (middle panel). Moreover, the sixth-order IFCs further improve the accuracy of the Taylor expansion potential as shown in the right panel of the figure. Since the correction by the sextic terms is minor, however, we only considered the dominant quartic terms in the SCP calculations.

Next, we calculated finite-temperature phonon dispersion curves by solving the SCP equation [Eq. (1)] at various temperatures, where all of the quartic coupling coefficients between zone center phonon modes were considered. The SCP dispersion curves at 300 K are shown in Fig.1(a) by solid lines. The quartic anharmonicity generally increases the frequencies of the low-lying rattling modes, which can be attributed to the dominant and positive contribution from the diagonal term of the quartic coefficient; for the low-lying optical modes . With increasing the cage size, the quartic component of the PES becomes more important as shown in Fig. 1(b), leading to the greater frequency shifts. For the higher-frequency modes above 80 cm, the anharmonic renormalization was turned out to be relatively small (see the SI sup ()). We note here that the present result is the first realization of the ab initio SCP calculation for a complex host-guest structure.

According to the Raman study of Takasu et al. Takasu et al. (2006), the frequency of the guest mode in BGG increases from 31 cm at 2 K to 34 cm at room temperature. The SCP theory for the “opt.” case gives 29.8 cm at 0 K and 35.5 cm at 300 K, which agree well with the experimental values especially given that the present SCP theory neglects the intrinsic frequency shift by the cubic anharmonicity and the quasiharmonic effect. To investigate the significance of the intrinsic frequency shift by the cubic terms, we also calculated the lowest-order correction by the cubic anharmonicity from the real part of the bubble self-energy as with being defined as

 Σ(B)q(ω) =ℏ2N∑∑′q1,q2,s=±1|Φ(−q;q1;q2)|28ΩqΩq1Ωq2 ×[n1+n2+1sωc+Ωq1+Ωq2−n1−n2sωc+Ωq1−Ωq2]. (2)

Here, is the number of points in the first Brillouin zone, , with being a positive infinitesimal, and the summation over are restricted to the pairs satisfying the momentum conservation; . For the “opt.” case, the value of the mode calculated with the 999 point mesh was cm at 0 K and cm at 300K. Hence, the experimental -dependence of the Raman shift can be better explained with the effect of the intrinsic frequency shift by the cubic anharmonicity. We have also calculated the values for the other systems with different cage sizes and found that the frequency shift by the bubble diagram is negative ( and less significant than the hardening by the quartic anharmonicity sup ().

To elucidate the intrinsic effects of the frequency renormalization on LTC quantitatively, we have calculated based on the Boltzmann transport equation (BTE) within the relaxation-time approximation, where

 κμνL(T)=1NV∑qcq(T)vμq(T)vνq(T)τq(T). (3)

Here, is the unit cell volume, is the mode specific heat, is the group velocity, and is the lifetime of phonon . Unlike the conventional BTE approach, where harmonic frequencies and eigenvectors are used as the ground state for calculating phonon group velocities and lifetimes, we employ the SCP frequencies and eigenvectors in the present SCP+BTE method [Eq. (3)]. Therefore, the group velocity also shows intrinsic -dependence. The phonon lifetime is estimated using the Matthiessen’s rule . The anharmonic scattering rate is calculated from the imaginary part of the bubble self-energy [Eq. (2)] as , and the phonon-isotope scattering rate is evaluated perturbatively Tamura (1983). For the phonon-boundary scattering rate, we employ with the grain size of m that reproduces the experimental crystalline peak of LTC Sales et al. (2001).

Figure 2 shows -dependence of LTC calculated by the BTE and the SCP+BTE methods with 999 points. The LTC values calculated by the SCP+BTE method are generally higher than those obtained by the conventional BTE method. This tendency becomes more pronounced with increasing and . The predicted value by the SCP+BTE method is 0.97 W/mK at 300 K, which agrees well with the experimental values of 1.31 W/mK (Sales et al., Ref. Sales et al., 2001) and 1.06 W/mK (May et al., Ref. May et al., 2009). On the other hand, the conventional method gives 0.58 W/mK, which is 40% smaller than the SCP+BTE value. Moreover, the deviation from in a high temperature range can be well reproduced by the SCP+BTE approach. These results clearly reveal the essential role of the frequency renormalization on the thermal transport properties of BGG. To elucidate the origin of the increase in due to the hardening of the guest modes, we compare the LTC spectrum for the “opt.” case at 300 K. The results are shown in Fig. 3(a). With the hardening of the low-lying guest modes, increases significantly in the low-frequency region ( cm). After careful investigation, we found that this increase can be ascribed to the decrease in the phonon linewidth , whose magnitude is roughly proportional to the available scattering phase space (SPS) and the strength of the cubic coupling . In the frequency range of 45–65 cm, the SPS decreases by % due to the hardening of the low-frequency guest modes sup (). This reduction results in the enhancement of in the same frequency range. In the low-frequency region below 45 cm, however, the change of the SPS is too small to explain the corresponding enhancement of . Indeed, we found that the coupling coefficient is also suppressed strongly due to the reduction of the optical-acoustic phonon hybridization caused by the hardenings of the optical modes (see the SI sup ()). Importance of such optical-acoustic hybridization was also reported in Refs. Tadano et al., 2015; Li et al., 2016.

Next, we discuss the volume dependence of the LTC calculated by SCP+BTE. As can be inferred from Fig. 2, the LTC values decrease with increasing the unit cell volume. It is interesting to observe that the crystalline peak of near 20 K evolved into a resonance dip when the the lattice constant was expanded up to %. In the “%” system, takes the minimum value at 15 K and increases up to 50 K, which is qualitatively different from the -dependence of the other systems. These computational results agree qualitatively with the experimental results on SrGaSiGe, where the cage size was controlled by changing the value Suekuni et al. (2007). To understand the microscopic origin of the increase of from 15 K to 50 K, we compare the calculated LTC spectra of the “%” and “%” systems in Figs. 3(b) and (c), respectively. In the both systems, the dispersive phonon modes below 50 cm contribute more than 90% to the total value at 20 K. When higher frequency phonon modes are thermally excited at 50 K, the contribution from the dispersive optical modes around 100 cm becomes significant. In the “%” system, the thermal enhancement of near 100 cm was smaller than the concurrent reduction in cm, which resulted in . When the lattice constant is changed from “%” to “%”, the frequencies of the low-lying guest modes decrease as discussed earlier. This softening causes suppression of , which is particularly significant in cm while the phonon lifetimes as well as in the higher frequency region (100 cm) are less affected sup (). Such a joint effect of strong acoustic-optical scatterings in the low-frequency range and the presence of higher-frequency dispersive optical modes resulted in for the “%” system. It is important to mention that the predicted values in the low- region depends on the employed grain size . We found that the resonance dip of the “%” system persisted even when we omitted the boundary scattering term. Moreover, the resonance dip changed into a plateau when we employed a smaller value, which is consistent with the experimental fact that the presence of the resonance dip is sensitive to details of sample preparation methods Nolas et al. (1998); Sales et al. (2001); Christensen et al. (2016).

Finally, we focus on the -dependence of above 50 K in Fig. 2. The -dependence obtained by the conventional BTE method follows for all of the studied systems, which do not agree with the experimental fact that shows weaker -dependence May et al. (2009); Suekuni et al. (2007). The SCP+BTE method considerably improved the agreement with the experimental results and produced milder temperature dependence. However, it was still inadequate to reproduce the increasing trend of observed in some clathrates having large guest-free-spaces. Since the SCP+BTE method is still based on the phonon-gas model, it does not include the nondiagonal Peierls contribution Auerbach and Allen (1984); Allen and Feldman (1993) and the anharmonic contribution Sun and Allen (2010) to the heat flux operator, which become generally more important at high temperature. Therefore, our results also indicate that these contributions may not be negligible in some clathrates at high temperature.

To summarize, we performed fully ab initio calculations of of type-I clathrate BGG with various different cage sizes, where both the three-phonon scattering and the frequency renormalization by the quartic anharmonicity were taken into account. We showed that the hardening of vibrational frequencies of rattling modes caused by the quartic anharmonicity significantly affects the calculated values, leading to an improved agreement with experimental results including the -dependence weaker than . In addition, we found that the evolution from crystal-like to glasslike near 20 K can be realized by our static calculations without disorders, which can be attributed to the presence of low-frequency guest modes that strongly couple with acoustic modes along with higher-frequency dispersive optical modes and the phonon-boundary scatterings. While our simulation does not exclude the possibility of disorders to further reinforce the glasslike behavior observed in real clathrate samples, it provides a new miscroscopic insight into the exceptional thermal transport of clathrates.

We thank Koichiro Suekuni and Stéphane Pailhès for fruitful discussions. This study is partly supported by JSPS KAKENHI Grant Number 16K17724, “Materials research by Information Integration” Initiative (MI2I) project of the Support Program for Starting Up Innovation Hub from Japan Science and Technology Agency (JST), and MEXT Element Strategy Initiative to Form Core Research Center in Japan. The computation in this work has been done using the facilities of the Supercomputer Center, Institute for Solid State Physics, The University of Tokyo.

Supplemental Materials

Quartic Anharmonicity of Rattlers and its Effect on Lattice Thermal Conductivity of Clathrates from First Principles

## Appendix S1 A. Details of computational procedure

In this study, we performed DFT calculations using the VASP code Kresse and Furthmüller (1996) with the projector augmented wave (PAW) method Blöchl (1994); Kresse and Joubert (1999). We employed the cutoff energy of 320 eV and 444 Monkhorst-Pack grid Monkhorst and Pack (1976) for the Brillouin zone integration. For the exchange-correlation functional, we employed the Perdew-Burke-Ernzerhof (PBE) functional Perdew et al. (1996). Before performing phonon calculations, we fully relaxed the lattice constant and the internal coordinates, where we obtained  Å. Then, we changed the lattice constant from the optimized value by %, %, %, and %, and relaxed the internal coordinates for them. The actual values of the lattice constants and the resulting pressures are shown in table S1.

After the relaxation, we calculated harmonic interatomic force constants (IFCs) using the finite displacement approach Parlinski et al. (1997). Here, we considered all harmonic IFCs inside the unit cell and employed the displacement length of Å. The anharmonic IFCs were calculated using the compressive sensing lattice dynamics method Zhou et al. (2014). To this end, we employed the following procedure: First, we performed ab initio molecular dynamics (AIMD) for the “%” system at 300 K and generated physically relevant atomic configurations. The AIMD simulation was performed with 222 Monkhorst-Pack grid and a relatively large convergence criterion for the SCF loop to accelerate the structural sampling. Second, we uniformly sampled 140 snapshots from the AIMD trajectory. In each sampled snapshot, we further displaced all atoms by 0.1 Å in random directions to decrease cross-correlations between the sampled configurations. Third, we performed static DFT calculations for the 140 snapshots and calculated Hellmann-Feynman forces accurately. These displacement and force data sets form a training data for estimating anharmonic IFCs of the “%” system. For the other systems with different unit cell volumes, the atomic snapshots were created simply by changing the lattice constants while leaving the internal coordinates of the original structures unaltered. In total, DFT calculations were conducted for calculating anharmonic IFCs of BaGaGe. Finally, we estimated anharmonic IFCs by using the least absolute shrinkage and selection operator (LASSO) Tibshirani (1996). The hyperparameter for the LASSO was selected from cross-validation scores Hastie et al. (2009), and the optimization problem was solved by using the coordinate descent method Hastie et al. (2015).

In the present study, we considered anharmonic IFCs up to the sixth order. To reduce the number of independent IFCs, we introduced cutoff radii for anharmonic terms so that only the host-guest (host-host) terms inside the first (second) nearest neighbor shell were included. In addition, we omitted less important four-body terms for the quartic IFCs and multi-body ( body) terms for the quintic and sextic IFCs. After removing linearly-dependent IFCs by using available symmetry operations and the constraints of the translational invariance, the numbers of irreducible IFCs turned out to be 994, 2819, 252, and 345 for cubic, quartic, quintic, and sextic terms, respectively. We estimated these 4410 terms by solving the -regularized linear regression problems. The calculated cubic and quartic IFCs were employed for the subsequent phonon calculations.

All of the phonon calculations and IFC estimations were performed by using the ALAMODE package Tadano et al. (2014); ala ().

## Appendix S2 B. Anharmonic phonon dispersion

The full anharmonic phonon dispersion curves calculated by the self-consistent phonon (SCP) theory are shown in Figs. S1S5. In these figures, we also show harmonic phonon dispersion curves for comparison.

## Appendix S3 C. Frequency shift by cubic anharmonicity

In table S2, the frequencies of the guest mode calculated with different conditions and approximations are summarized. The “SCP+Bubble” corresponds to the value of where . For the “%” system, we could not reach a convergence of the SCP equation at 0 K within a reasonable computational time. This is reasonable because the present SCP method does not allow the spontaneous symmetry breaking to occur and the thermally averaged equilibrium positions of the rattlers are assumed to be the center of the cage irrespective of temperature. For off-center systems, such an assumption is invalid in a low-temperature region, but should be more or less reasonable at high temperatures.

## Appendix S4 D. Phonon lifetime

The calculated phonon lifetimes of BaGaGe are shown in Figs. S6S10 as a function of phonon frequencies. In each figure, we compare two computational results; one is based on the harmonic approximation and the other is based on the SCP theory. For the “%” system, we only show the result based on the SCP theory. To see the validity of the quasiparticle picture, we also show by dashed lines. When is not satisfied, the quasiparticle picture becomes questionable.

## Appendix S5 E. Change of scattering phase space and cubic coupling by frequency renormalization

We show the change of the three phonon scattering phase space (SPS) and the strength of the cubic coupling coefficients caused by the hardening of low-lying guest modes. Here, all of the results were calculated with the optimized lattice constant and at 300 K.

The three-phonon SPS of phonon is defined as

 Unknown environment 'array% (S1)

When the SPS was calculated based on the SCP solution, the harmonic phonon frequencies were replaced by the temperature dependent anharmonic phonon frequencies, i.e., . The calculated results are shown in Fig. S11.

To highlight the importance of the change of anharmonic coupling coefficients, we also compare the components of calculated with and without the temperature renormalization of phonon frequencies. The computational results are shown in Fig. S12. Here, we selected transverse acoustic (TA) and longitudinal acoustic (LA) modes at . The frequencies of these two acoustic modes are 13.8 and 22.1 cm within the harmonic approximation and 14.1 and 22.7 cm within the self-consistent phonon theory. Also, the components of are shown only for the pairs satisfying the momentum- and energy-conservation laws.

## References

### References

1. Y. Takasu, T. Hasegawa, N. Ogita, M. Udagawa, M. A. Avila, K. Suekuni, I. Ishii, T. Suzuki,  and T. Takabatake, Phys. Rev. B 74, 174303 (2006).
2. M. Christensen, A. B. Abrahamsen, N. B. Christensen, F. Juranyi, N. H. Andersen, K. Lefmann, J. Andreasson, C. R. H. Bahl,  and B. B. Iversen, Nature Mater. 7, 811 (2008).
3. T. Takabatake, K. Suekuni, T. Nakayama,  and E. Kaneshita, Rev. Mod. Phys. 86, 669 (2014).
4. J. Wu, K. Akagi, J. Xu, H. Shimotani, K. K. Huynh,  and K. Tanigaki, Phys. Rev. B 93, 094303 (2016).
5. J. Dong, O. F. Sankey,  and C. W. Myles, Phys. Rev. Lett. 86, 2361 (2001).
6. T. Tadano, Y. Gohda,  and S. Tsuneyuki, Phys. Rev. Lett. 114, 095501 (2015).
7. B. C. Sales, B. C. Chakoumakos, R. Jin, J. R. Thompson,  and D. Mandrus, Phys. Rev. B 63, 245113 (2001).
8. M. A. Avila, K. Suekuni, K. Umeo, H. Fukuoka, S. Yamanaka,  and T. Takabatake, Phys. Rev. B 74, 125109 (2006).
9. A. May, E. Toberer, A. Saramat,  and G. Snyder, Phys. Rev. B 80, 1 (2009).
10. G. S. Nolas, J. L. Cohn, G. A. Slack,  and S. B. Schujman, Appl. Phys. Lett. 73, 178 (1998).
11. M. A. Avila, K. Suekuni, K. Umeo, H. Fukuoka, S. Yamanaka,  and T. Takabatake, Appl. Phys. Lett. 92, 1901 (2008).
12. K. Suekuni, M. A. Avila, K. Umeo,  and T. Takabatake, Phys. Rev. B 75, 195210 (2007).
13. K. Suekuni, Y. Takasu, T. Hasegawa, N. Ogita, M. Udagawa, M. A. Avila,  and T. Takabatake, Phys. Rev. B 81, 205207 (2010).
14. F. Bridges and L. Downward, Phys. Rev. B 70, 140201 (2004).
15. Q. Xi, Z. Zhang, J. Chen, J. Zhou, T. Nakayama,  and B. Li, Phys. Rev. B 96, 064306 (2017).
16. S. Christensen, M. S. Schmøkel, K. A. Borup, G. K. H. Madsen, G. J. McIntyre, S. C. Capelli, M. Christensen,  and B. B. Iversen, J. Appl. Phys. 119, 185102 (2016).
17. R. E. Peierls, Quantum Theory of Solids (Oxford University Press, 1955).
18. L. Lindsay, Nanoscale and Microscale Thermophysical Engineering 20, 67 (2016).
19. P. B. Allen and J. L. Feldman, Phys. Rev. B 48, 12581 (1993).
20. P.-F. Lory, S. Pailhès, V. M. Giordano, H. Euchner, H. D. Nguyen, R. Ramlau, H. Borrmann, M. Schmidt, M. Baitinger, M. Ikeda, P. Tomeš, M. Mihalkovič, C. Allio, M. R. Johnson, H. Schober, Y. Sidis, F. Bourdarot, L.-P. Regnault, J. Ollivier, S. Paschen, Y. Grin,  and M. Boissieu, Nature Commun. 8, 1 (2017).
21. P. Souvatzis, O. Eriksson, M. I. Katsnelson,  and S. P. Rudin, Phys. Rev. Lett. 100, 95901 (2008).
22. I. Errea, M. Calandra,  and F. Mauri, Phys. Rev. B 89, 064302 (2014).
23. T. Tadano and S. Tsuneyuki, Phys. Rev. B 92, 054301 (2015).
24. A. van Roekeghem, J. Carrete,  and N. Mingo, Phys. Rev. B 94, 020303 (2016).
25. K. Momma and F. Izumi, J. Appl. Crystallogr. 44, 1272 (2011).
26. T. Tadano and S. Tsuneyuki, arXiv:1706.04744, to appear in J. Phys. Soc. Jpn. .
27. G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
28. J. P. Perdew, K. Burke,  and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
29. P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
30. G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
32. T. Tadano, Y. Gohda,  and S. Tsuneyuki, J. Phys: Condens. Matter 26, 225402 (2014).
33. F. Zhou, W. Nielson, Y. Xia,  and V. Ozoliņš, Phys. Rev. Lett. 113, 185501 (2014).
34. See Supplemental Material at [URL] for details of computational conditions, anharmonic phonon dispersion relations, phonon lifetimes, and a comparison of the scattering phase space and the cubic coupling coefficients between BTE and SCP+BTE.
35. S.-I. Tamura, Phys. Rev. B 27, 858 (1983).
36. W. Li, J. Carrete, G. K. H. Madsen,  and N. Mingo, Phys. Rev. B 93, 205203 (2016).
37. A. Auerbach and P. B. Allen, Phys. Rev. B 29, 2884 (1984).
38. T. Sun and P. B. Allen, Phys. Rev. B 82, 224305 (2010).
39. G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
40. P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
41. G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
42. H. Monkhorst and J. Pack, Phys. Rev. B 13, 5188 (1976).
43. J. Perdew, K. Burke,  and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
44. K. Parlinski, Z. Q. Li,  and Y. Kawazoe, Phys. Rev. Lett. 78, 4063 (1997).
45. F. Zhou, W. Nielson, Y. Xia,  and V. Ozoliņš, Phys. Rev. Lett. 113, 185501 (2014).
46. R. Tibshirani, J. R. Stat. Soc. Series B Stat. Methodol. 58, 267 (1996).
47. T. Hastie, R. Tibshirani,  and J. Friedman, “The Elements of Statistical Learning,” Springer New York, New York, NY (2009).
48. T. Hastie, R. Tibshirani,  and M. Wainwright, Statistical Learning with Sparsity: The Lasso and Generalizations (Chapman & Hall/CRC, 2015).
49. T. Tadano, Y. Gohda,  and S. Tsuneyuki, J. Phys: Condens. Matter 26, 225402 (2014).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters