Hopping Processes Explain T-linear Rise of Thermal Conductivity in Thermoelectric Clathrates above the Plateau

TL;DR

This paper explains the T-linear increase in thermal conductivity of off-center type-I clathrates above the plateau through a hopping mechanism of localized modes, aligning well with experimental observations.

Contribution

It introduces a hopping mechanism for localized modes that accounts for the T-linear thermal conductivity rise in off-center clathrates, supported by theoretical calculations.

Findings

Hopping of localized modes explains T-linear thermal conductivity rise.

Calculated results match experimental data for off-center clathrates.

Both magnitude and temperature dependence are accurately reproduced.

Abstract

Type-I clathrate compounds with off-center guest ions realize the phonon-glass electron-crystal concept by exhibiting almost identical lattice thermal conductivities $\kappa_{\rm L}$ to those observed in network-forming glasses. This is in contrast with type-I clathrates with on-center guest ions showing $\kappa_{\rm L}$ of conventional crystallines. Glasslike $\kappa_{\rm L}$ stems from the peculiar THz frequency dynamics in off-center type-I clathrates where there exist three kinds of modes classified into extended(EX), weakly(WL) and strongly localized(SL) modes as demonstrated by Liu $et.\,al.$, Phys. Rev. B $\textbf{93}$, 214305(2016). Our calculated results based on the hopping mechanism of SL modes via anharmonic interactions show fairly good agreement with observed ${\it T}$-linear rise of $\kappa_{\rm L}$ above the plateau. We emphasize that both the magnitude and the temperature dependence are in accord with the experimental data of off-center type-I clathrates.

Equations

$$P({\epsilon }_q)=\frac{{\left({\sum }_{\ell =1}^N{\left|{\phi }_{\ell }({\epsilon }_q)\right|}^2\right)}^2}{N{\sum }_{\ell =1}^N{\left|{\phi }_{\ell }({\epsilon }_q)\right|}^4},$$

$$H_0=K_{\text{C}}+K_{\text{G}}+V_{\text{CC}}+V_{\text{CG}}.$$

$$K=\frac{1}{2}\sum _{\ell }\left(M{\dot{\text{r}}}_{\ell }(t{)}^2+m{\dot{\text{u}}}_{\ell }(t{)}^2\right),$$

$$V_{CC}=\sum _{{\ell }^{\prime }>\ell ,\mu }\frac{f_{\ell ,{\ell }^{\prime },\mu }}{2}({\text{r}}_{\ell ,\mu }(t)-{\text{r}}_{{\ell }^{\prime },\mu }(t){)}^2,$$

$$V_{CG}=\sum _{\ell ,m=\mathrm{in},\mathrm{out}}\frac{{\xi }_m}{2}{\text{w}}_{\ell ,m}^2(t),$$

$$V_{CG}$$

$${\xi }_{r,(\varphi )}$$

$$\begin{array}{r}{V}_{CG}^{\prime }=\frac{1}{2}\sum _{\ell ,\alpha \ne \beta }(D_r{e}_{\alpha \alpha }+S_r{e}_{\alpha \beta })({\hat{\mathbf{U}}}_{\ell }\cdot {\mathbf{w}}_{\parallel ,\ell }{)}^2\\ +\frac{1}{2}\sum _{\ell ,\alpha \ne \beta }(D_{\varphi }{e}_{\alpha \alpha }+S_{\varphi }{e}_{\alpha \beta })({\hat{\mathbf{U}}}_{\ell }\times {\mathbf{w}}_{\parallel ,\ell }{)}^2.\end{array}$$

$${\text{r}}_{\ell }(t)=\sum _{{\text{k}}_{\mu }}\sqrt{\frac{\hslash }{2\rho {\omega }_{{\text{k}}_{\mu }}}}{\hat{\mathbf{e}}}_{{\text{k}}_{\mu }}\left({\varphi }_{{\text{k}}_{\mu }}({\text{R}}_{\ell })b{}_{{\text{k}}_{\mu }}^{†}(t)+h.c.\right).$$

$${\varphi }_{{\text{k}}_{\mu }}({\text{R}}_{\ell })=\sqrt{\frac{1}{V}}{e}^{i{\text{k}}_{\mu }\cdot {\text{R}}_{\ell }}.$$

$$\int \mid {\varphi }_{{\text{k}}_{\mu }}({\text{R}}_{\ell }){\mid }^{2}d{\text{R}}_{\ell }=1.$$

$${\text{w}}_{\ell }(t)=\sum _{\lambda }\sqrt{\frac{\hslash }{2m^{\prime }{\omega }_{\lambda }}}{\hat{\mathbf{e}}}_{\lambda }\left({\psi }_{\lambda }({\text{R}}_{\ell })c_{\lambda }^{†}(t)+h.c.\right).$$

$${\psi }_{\lambda }({\text{R}}_{\ell })=A\cos \left[{\text{k}}_{\lambda }\cdot ({\text{R}}_{\ell }-{\text{R}}_{\lambda })\right]e^{-|{\text{R}}_{\ell }-{\text{R}}_{\lambda }|/L_{\lambda }}.$$

$$\sum _{\ell }\mid {\psi }_{\lambda }({\text{R}}_{\ell }){\mid }^2=\frac{1}{\Omega }\int d{\text{R}}_{\ell }\mid {\psi }_{\lambda }({\text{R}}_{\ell }){\mid }^2=1,$$

$$A\cong \sqrt{\frac{2\Omega }{\pi L_{\lambda }^3}}.$$

$${\psi }_{\lambda }({\text{R}}_{\ell })=\sqrt{\frac{2\Omega }{\pi L_{\lambda }^3}}\cos \left[{\text{k}}_{\lambda }\cdot ({\text{R}}_{\ell }-{\text{R}}_{\lambda })\right]e^{-|{\text{R}}_{\ell }-{\text{R}}_{\lambda }|/L_{\lambda }}.$$

$$\begin{array}{rl}{H}_{CG}^{\prime }& =\sum _{{\text{k}}_{\mu },{\lambda }^{\prime },{\lambda }^{\prime \prime }}(A_{{\text{k}}_{\mu },{\lambda }^{\prime },{\lambda }^{\prime \prime }}{b}_{{\mathbf{k}}_{\mu }}{c}_{{\lambda }^{\prime }}c{}_{{\lambda }^{\prime \prime }}^{†}+h.c.),\\ & +\sum _{{\text{k}}_{\mu },{\lambda }^{\prime \prime \prime },{\lambda }^{\prime \prime \prime \prime }}(B_{{\mathbf{k}}_{\mu },{\lambda }^{\prime \prime \prime },{\lambda }^{\prime \prime \prime \prime }}{b}_{{\mathbf{k}}_{\mu }}{c}_{{\lambda }^{\prime \prime \prime }}c{}_{{\lambda }^{\prime \prime \prime \prime }}^{†}+h.c.)\\ & +\sum _{{\text{k}}_{\mu },{\lambda }^{\prime },{\lambda }^{\prime \prime \prime }}(C_{{\mathbf{k}}_{\mu },{\lambda }^{\prime },{\lambda }^{\prime \prime \prime }}{b}_{{\mathbf{k}}_{\mu }}{c}_{{\lambda }^{\prime }}c{}_{{\lambda }^{\prime \prime \prime }}^{†}+h.c.),\end{array}$$

$$\begin{array}{rl}& A_{{\text{k}}_{\mu },{\lambda }^{\prime },{\lambda }^{\prime \prime }}=-\frac{1}{4}\sum _{l}i\sqrt{\frac{\hslash }{2\rho {\omega }_{{\text{k}}_{\mu }}}}\sqrt{\frac{\hslash }{2m^{\prime }{\omega }_{{\lambda }^{\prime }}}}\sqrt{\frac{\hslash }{2m^{\prime }{\omega }_{{\lambda }^{\prime \prime }}}}\\ & \times {\varphi }_{{\text{k}}_{\mu }}{\psi }_{{\lambda }^{\prime }}{\psi }_{{\lambda }^{\prime \prime }}\left[(D_r+D_{\varphi })k_{\parallel }{\delta }_{\mu ,\parallel }+(S_r+S_{\varphi })k_{\perp }{\delta }_{\mu ,\perp }\right],\end{array}$$

$$\begin{array}{rl}& C_{{\text{k}}_{\mu },{\lambda }^{\prime },{\lambda }^{\prime \prime \prime }}=-\frac{1}{\pi }\sum _{l}i\sqrt{\frac{\hslash }{2\rho {\omega }_{{\text{k}}_{\mu }}}}\sqrt{\frac{\hslash }{2m^{\prime }{\omega }_{{\lambda }^{\prime }}}}\sqrt{\frac{\hslash }{2m^{\prime }{\omega }_{{\lambda }^{\prime \prime \prime }}}}\\ & \times {\varphi }_{{\text{k}}_{\mu }}{\psi }_{{\lambda }^{\prime }}{\psi }_{{\lambda }^{\prime \prime \prime }}\left[(D_r-D_{\varphi })k_{\parallel }{\delta }_{\mu ,\parallel }+(S_r-S_{\varphi })k_{\perp }{\delta }_{\mu ,\perp }\right].\end{array}$$

$$\begin{array}{rl}{A}_{{\mathbf{k}}_{\mu },{\lambda }^{\prime },{\lambda }^{\prime \prime }}^2=& \frac{C{I}_1^2}{V{L}_{{\lambda }^{\prime }}^3{L}_{{\lambda }^{\prime \prime }}^3}\frac{1}{{\omega }_{{\text{k}}_{\mu }}{\omega }_{{\lambda }^{\prime }}{\omega }_{{\lambda }^{\prime \prime }}}\\ & {\left[(D_r+D_{\varphi })k_{\parallel }{\delta }_{\mu ,\parallel }+(S_r+S_{\varphi })k_{\perp }{\delta }_{\mu ,\perp }\right]}^2,\end{array}$$

$$C=\frac{{\hslash }^3{\Omega }^2}{2^5\rho m^{\prime 2}}.$$

$$\frac{d{n}_{{\lambda }^{\prime }}}{dt}$$

$$\frac{1}{{\tau }_{{\lambda }^{\prime }}^{\mathrm{same}}}$$

$$I_1$$

$$I_1=\sum _{\ell }f({\text{R}}_{\ell })f({\text{R}}_{\ell }-\Delta {\text{R}}_{{\lambda }^{\prime \prime }})e^{-i{\text{k}}_{\mu }\cdot {\text{R}}_{\ell }},$$

$$f({\text{X}}_{\ell })=\cos ({\text{k}}_{{\lambda }^{\prime }}\cdot {\text{X}}_{\ell })e^{-|{\text{X}}_l|/L_{{\lambda }^{\prime }}}.$$

$$\begin{array}{rl}{I}_1& \cong \frac{1}{\Omega }{\int }_{\mid {\text{X}}_{\ell }\mid <L}d{\text{R}}_{\ell }f({\text{R}}_{\ell })f({\text{R}}_{\ell }-\Delta {\text{R}}_{{\lambda }^{\prime \prime }})e^{-i{\text{k}}_{\mu }\cdot {\text{R}}_{\ell }}\\ & \cong \frac{1}{\Omega }{\int }_{\mid {\text{X}}_{\ell }\mid <L}d{\text{R}}_{\ell }{e}^{\frac{-|{\text{R}}_{\ell }|-|{\text{R}}_{\ell }-\Delta {\text{R}}_{{\lambda }^{\prime \prime }}|}{L}}{e}^{-i{\text{k}}_{\mu }\cdot {\text{R}}_{\ell }}\\ & \left[\frac{1}{2}\cos \left(2{\text{k}}_{{\lambda }^{\prime }}\cdot {\text{R}}_{\ell }-{\text{k}}_{{\lambda }^{\prime }}\cdot \Delta {\text{R}}_{{\lambda }^{\prime \prime }}\right)+\frac{1}{2}\cos \left({\text{k}}_{{\lambda }^{\prime }}\cdot \Delta {\text{R}}_{{\lambda }^{\prime \prime }}\right)\right]\\ & \cong |\Delta {\text{R}}_{{\lambda }^{\prime \prime }}|\pi L^2\frac{1}{2\Omega }{e}^{-|\Delta {\text{R}}_{{\lambda }^{\prime \prime }}|/L}.\end{array}$$

$$I_1^2\simeq {\left(\frac{\pi \Delta R_{{\lambda }^{\prime \prime }}{L}^2}{2\Omega }\right)}^2{e}^{-2\Delta R_{{\lambda }^{\prime \prime }}/L},$$

$$\begin{array}{rl}\frac{1}{{\tau }_{{\lambda }^{\prime }}^{\mathrm{same}}}& \cong \frac{2\pi k_{B}T(D_r+D_{\varphi }{)}^{2}C{I}_1^2}{{\hslash }^{3}V{L}^6{v}_{\parallel }^2}\\ & \times \sum _{{\text{k}}_{\parallel },{\lambda }^{\prime \prime }}\left[\frac{\delta ({\omega }_{{\lambda }^{\prime \prime }}-{\omega }_{{\lambda }^{\prime }}-{\omega }_{{\text{k}}_{\parallel }})}{{\omega }_{{\lambda }^{\prime \prime }}^2}+\frac{\delta ({\lambda }^{\prime }\rightleftarrows {\lambda }^{\prime \prime })}{{\omega }_{{\lambda }^{\prime \prime }}^2}\right]\\ & +\left[2\times \left(D⟶S,\parallel ⟶\perp \right)\mathrm{in}\mathrm{the}\mathrm{above}\right].\end{array}$$

$$\begin{array}{r}{\kappa }_{\mathrm{hop}}(T)=\frac{1}{3V}\sum _{{\lambda }^{\prime }}{C}_{{\lambda }^{\prime }}(T)\frac{\Delta R_{{\lambda }^{\prime \prime }}^2}{{\tau }_{{\lambda }^{\prime }}},\end{array}$$

Full text

I INTRODUCTION

Lattice thermal conductivity constitutes a key element to improve the efficiency of the thermal-to-electrical conversion in thermoelectric (TE) devices as understood from the material’s figure of merit describing the efficiency $Z=S^{2}\sigma/\kappa_{\rm tot}$ [K*-1*]. The numerator contains the Seebeck coefficient $S(T)\,$[V/K] and the electrical conductivity $\sigma(T)\,$[1/($\Omega$m)], while the denominator $\kappa_{\rm tot}(T)$ [W/(mK)] consists of the sum of electrical $\kappa_{\rm el}$ and lattice $\kappa_{\rm L}$ thermal conductivity. Hence, the high performance of thermoelectricity can be achieved for materials with the lowest possible thermal conductivity $\kappa_{\rm tot}$, the highest possible electrical conductivity $\sigma$ and the highest possible Seebeck coefficient $S$. Provided that the Wiedemann-Franz law $\kappa_{\rm el}(T)\propto\sigma(T)$ holds for, $\kappa_{\rm L}$ becomes a crucial parameter to improve the performance of TE conversion. In this framework, Slack Slack1995 has proposed the concept of “phonon-glass electron-crystal”. This has been one of guiding principles for exploring high-performance TE materials Takabatake2014 ; Beekman2015 .

Type-I clathrates with “off-center” guest ions, such as R8Ga16Ge30(R=Ba, Sr, Eu) Nolas:1998aa ; Cohn:1999a ; Christensen:2016a ; Paschen:2001a ; Sales:2001a ; Avila:2006a , Ba8Ga16Sn30 Avila:2008aa ; Suekuni:2008a , Sr8Ga16Si30-xGex Suekuni:2007a , are particularly interesting in this respect since these systems exhibit almost identical lattice thermal conductivities to those of structural glasses, which consist of four specific regions characterized by: (i) T∼2-dependence below a few Kelvin, (ii) the plateau region between a few K and a few 10K, and (iii) the subsequent rise proportional to T, and (iv) its saturation above *T$\sim$*100K. These characteristics of $\kappa_{\rm L}$ exhibit a remarkable uniformity which appears to be insensitive to chemical compositions, suggesting the existence of a unified mechanism Nakayama2002 . However, this issue remains as an open and challenging problem of long-standing due to the difficulty to identify relevant entities or elements at atomistic level caused by their complex microscopic structures. Surprisingly enough, though “off-center” clathrates are crystalline with regularly network structure, the temperature dependence as well as the magnitudes of their thermal conductivities are almost identical to those of structural glasses over the full temperature range. In contrast, type-I clathrates with “on-center” guest ions show conventional crystalline $\kappa_{\rm L}$ Takabatake2014 .

This paper is organized as follows. Section II surveys the characteristics of vibrational modes according to the results of the spectral density of states, eigenvalues and their eigenvectors Liu2016 . We claim in this Section that the onset of the plateau is due to the delocalization-localization (weak localization) transition of acoustic modes. In addition, we point out that the temperature region showing the T-linear rise subsequent to the plateau is associated with the energy range where SL modes are fully excited. Section III describes the construction of anharmonic interaction Hamiltonian between SL and EX modes. The second quantized form of anharmonic Hamiltonian is given in Section IV. Section V develops a theory on the mechanism governing the T-linear rise of $\kappa_{\rm L}(T)$ above a few 10 K. Excited modes in this temperature region are mostly strongly-localized (SL) modes satisfying the Ioffe-Regel condition as evident from the mode pattern obtained by large-scale numerical simulations Liu2016 . These are hybridized modes between acoustic phonons associated with network cages and local vibrations of guest ions in cages. Based on these numerical evidences, we explain in quantitative manner $\kappa_{\rm L}(T)$ proportional to T observed above the plateau, by introducing the quantum mechanical process of hopping of SL modes due to anharmonic interactions, first proposed for fracton excitations Alexander1987 . Summary and conclusions are given in Sec. VI.

II CHARACTERISTICS OF EXCITED PHONONS AT THz FREQUENCY REGION

Type-I clathrates form a primitive cubic structure ($Pm\bar{3}n$) consisting of 6 tetrakaidecahedron (14-hedrons) and 2 dodecahedron (12-hedrons) per unit cell, in which the group-I or -II elements in the periodic table are encaged in the polyhedrons as guest ions. See Fig. 1. The THz frequency phonon dynamics of off-center type-I clathrates has been investigated in terms of large-scale numerical simulations. They have illustrated type-I Ba8Ga16Sn30  (BGS) exhibiting glasslike $\kappa_{\rm L}(T)$ as a prototype material with off-center guest ions, in which the guest ion Ba(2) in tetrakaidecahedron cage has the mass $m$ and the molecular unit composed of one tetrakaidecahedron and 1/3 dodecahedron does the total mass $M$ excluding the off-center guest ion. The coarse-grained picture, an operation of reducing the degrees of freedom of the original system, is valid for our purpose from the following reasons. First of all, EX acoustic modes at THz frequencies play a dominant role in heat transport since optical modes concerning to the vibrations of cages themselves do not contribute to thermal conductivity. Second, the wave-length $\lambda$ of phonons in the frequency regime $\nu$ $\leq$ 2.5 THz ($E$ $\leq$ 10 meV) becomes $\lambda$ $\geq$ 1.6 nm, which is larger than the size of a unit cell of $a_{0}$ $\simeq$1 nm in type-I clathrates, as estimated from the relation $\lambda$=$v/\nu$ using the sound velocity $v\approx 4\times 10^{3}$[m/sec]. These validate the coarse-grained Hamiltonian for describing THz frequency dynamics rather than treating all microscopic constituents as equally relevant degrees of freedom.

Extremely large system-sizes are required in computer simulations on disorder systems in order to distinguish localized modes from extended modes. However, the present status of first-principles calculations (FPC) are limited to insufficient system-sizes for properly incorporating the disorder attributing to off-centeredness of guest atoms in off-center type-I clathrates consisting of a unit cell with ‘54’ atoms. Thus, it is difficult not only to include realistic disorder reproducing glasslike thermal conductivities, but also to exclude finite size effect for propagating acoustic phonons. Liu $et.~{}al.$ Liu2016 have performed calculations for 3D systems of (20$\times$20$\times$20)$\sim$(100$\times$100$\times$100) molecular units, for which they have employed a powerful numerical method called the forced oscillator method.Williams1985 ; Nakayama2001 They have also studied the localization nature of exited modes by taking the participation ratio (PR) as a criterion. The PR of a relevant mode {$\varphi_{\ell}(\varepsilon_{q});\ell=1,2,...N\}$ belonging to the eigenenergy $\varepsilon_{q}$ is defined by

[TABLE]

where $\ell$ denotes the $\ell$-th molecular unit depicted in Fig. 1 (d) and $N$ is the total mode number. For EX modes in a finite system, $P(\varepsilon_{q})$ take values close to $\approx$0.6 when $\varepsilon_{q}\neq 0$, and $P(\varepsilon_{q})$ becomes $\approx 1/N$ for SL modes Quasicrystals . Figure 2 (a) is the calculated phonon density of states (DOS), and (b)  the results of $P(\varepsilon_{q})$ for the size of 20$\times$20$\times$20 lattice of off-center type-I BGS. It is remarkable that $P(\varepsilon_{q})$ ranges from a value of SL modes $P(\varepsilon_{q})\approx 0$ to EX modes of $P(\varepsilon_{q})\approx 0.6$. We should emphasize that there appear three kinds of modes in the THz frequency region and below classified into EX, WL and SL modes. SL modes with PR values much smaller than unity are realized in the energy range from 2 to 3 meV as found from calculated mode patterns. Figure  3 depicts the mode patterns of SL mode at $\varepsilon_{q}$=2.6 meV.

The calculations of the PR for excited modes depicted in Fig. 2 have demonstrated that there exists the delocalization-localization transition at a “finite” frequency $\omega_{c}$ distinguishing EX and WL modes with the nature of acoustic modes vibrating “in-phase” between guest ions and cages. Furthermore, it has been found Liu2016 that WL modes convert to SL modes at higher frquencies with the nature of optical modes vibrating “out-of-phase” between guest ions and cages. In this aspect, we note that Nakayama Nakayama1998 had demonstrated the clear existence of the transition from WL to SL modes for the quasi-one-dimensional (1D) coarse-grained model consisting of host network and guest atoms connected by random springs. It was found Nakayama1998 that WL modes vibrate in-phase between network atoms and guest atoms, while SL modes manifest optical modes vibrating out-of-phase. However, there is no EX modes due to “quasi-1D” model. This manifests the Anderson weak localization criteria where the critical frequency $\omega_{c}$ takes a finite value in three dimensional (3D) systems, while it vanishes for 1D and 2D systems suggesting no EX modes in 1D and 2D disordered systems. The quasi-1D model Nakayama1998 should be thought as the simplest theoretical model for cage-guest systems with broad implication for the dynamics of cage-guest systems.

The observed delocalization-localization transition at $\varepsilon_{q}\thickapprox$ 1.3 meV accords with the observed onset temperature of the plateau of $\kappa_{\rm L}$ in BGS at $T_{\rm P}\thickapprox$1.3 meV/3.84$k_{\rm{B}}\thickapprox$ 3.9 K as estimated from the Wien’s displacement law for lattice thermal conductivities. Thus, the onset of the plateau is apparently due to the weak localization of acoustic modes. The plateau region should be interpreted as the contribution of EX phonons “saturates” at $T_{\rm P}$ for off-center type-I BGS. We note here that the random orientation of guest ions in cages plays a crucial role to the localization.

With increasing temperature further above a few 10K, $\kappa_{\rm L}$ show a linear rise on temperature Takabatake2014 . This type of anomalous thermal conductivities characterized by the plateau and the subsequent T-linear rise of thermal conductivities have been clearly observed for off-center type-I clathrates Nolas:1998aa ; Cohn:1999a ; Sales:2001a ; Paschen:2001a ; Suekuni:2007a ; Avila:2008aa ; Suekuni:2008a . SL modes are fully excited above the temperature *T$\simeq$*10K$\approx$3 meV/3.84$k_{\mathrm{B}}$ from the Wien’s displacement law. This indicates that T-linear rise subsequent to the plateau attributes to the excitations of SL modes. In the following Sections, we present the theoretical interpretation on the underlying mechanism of the linear rise on temperature above the plateau region for $\kappa_{\rm L}$.

III Coarse-grained Hamiltonian for Type-I Off-center Clathrates

III.1 Harmonic Hamiltonian

The Hamiltonian for off-center type-I clathrates under a coarse-grained picture consists of the kinetic energy of networked cages $K_{C}$ and off-center guest ions in cages $K_{G}$ in addition to the potential energy of the cage-cage interaction $V_{CC}$ and the cage-guest interaction $V_{CG}$. This is expressed by

[TABLE]

The explicit form of the total kinetic energy is given by the sum of $K_{C}$ and $K_{G}$ such as

[TABLE]

where $m$ and $M$ are masses of the guest ion in tetrakaidecahedron cage and the remained molecular unit, respectively. The vectors $\textbf{{r}}_{\ell}(t)$ and $\textbf{{u}}_{\ell}(t)$ represent small displacements of cage and guest ion from their equilibrium positions, $\textbf{{R}}_{\ell}$ and $\textbf{{R}}_{\ell}+\textbf{{U}}_{\ell}$, at the site $\ell$ as depicted in Fig.  4. Note here that guest ions take random orientation $\textbf{{U}}_{\ell}(\phi_{\ell})$ in tetrakaidecahedron cages.

The molecular unit composed of tetrakaidecahedron and dodecahedron is elastically connected with neighboring ones by the force constants $f_{\|},f_{\perp}$. These are related to the sound velocities of longitudinal ($\mu=\|$) and transverse ($\mu=\perp$) acoustic modes via the relation $v_{\mu}=a[f_{\mu}/(m+M)]^{1/2}$ with $a=a_{0}/2$ where $a_{0}$ is the lattice spacing of primitive cubic structure ($Pm\bar{3}n$) of type-I clathrates. Thus, we can estimate the force constants $f_{\|},f_{\perp}$ from the observed data of sound velocities. Note here that 6 molecular units are included in unit cell in type-I clathrates. In terms of these quantities, the potential energy of network cages becomes

[TABLE]

where $\mu=\|,\perp,\perp^{\prime}$. Hereafter, we keep up to the nearest neighbor coupling ($\ell^{\prime}=\ell+1$) between molecular units, which are denoted by $f_{\|}$, $f_{\perp}$ and $f_{\perp^{\prime}}$. The effect of randomly orientated guest ions are included in the following cage-guest interaction Hamiltonian.

The Hamiltonian should satisfy the symmetry of infinitesimal translation-invariance as a whole, $i.e.$, $\textbf{{r}}_{\ell}=\textbf{{u}}_{\ell}=\delta\textbf{{a}}$, which guarantees acoustic phonons as the Nambu-Goldstone boson with the eigenfrequency $\omega_{\textbf{{k}}}\rightarrow 0$ for $\textbf{{k}}\rightarrow 0$. This symmetry principle also holds for the potential of cage-guest interaction. Hence, the potential function for the cage-guest interaction $V_{CG}$ should be given by relative coordinates between the cage and the guest ion of $\textbf{{w}}_{\ell}(t)=\textbf{{u}}_{\ell}(t)-\textbf{{r}}_{\ell}(t)$, which is expressed by

[TABLE]

where $\xi_{m}$ represents the force constants between cage and guest ion depending on in-plane (parallel) or out-of-plane motion (perpendicular) to the hexagonal face in the tetrakaidecahedron cage. The guest ions execute in-plane vibration parallel to $x-y$ plane in addition to out-of-plane motions Avila:2008aa because of the anisotropic shape of tetrakaidecahedron cages. This is because off-center guest ions are involved in tetrakaidecahedron cages whose shape distinguishes the vibrations of off-center guest ion(2) in the plane parallel and perpendicular to the hexagonal face of the cage. Mori et al. Mori2011 observed by means of THz time-domain spectroscopy that the lowest-lying peak of off-center BGS at 0.71 THz splits into double peaks, $\omega_{0}^{\phi}/2\pi$=0.5THz and $\omega_{0}^{r}/2\pi$=0.72THz for off-center type-I BGS below T$\simeq$100 K. These spectra should be assigned to the libration and stretching modes of Ba(2) associated with $\xi_{\phi}$ and $\xi_{r}$. The peak around 1.35 THz is assigned as the out-of-plane motion of Ba(2) to the hexagonal faces of tetrakaidecahedron, which should be concerned with $\xi_{\theta}$. The Raman spectra of off-center $\rm Sr_{8}Ga_{16}Ge_{30}$ (SGG) have observed A1g stretching mode as 48 cm*-1*, and for off-center $\rm Eu_{8}Ga_{16}Ge_{30}$ (EGG) as 36 cm*-1* at 2 K Takasu2006 . Using these data, we can estimate the force constants via the relation $\xi_{r,(\phi,\theta)}=m^{\prime}\omega_{r,(\phi,\theta)}^{2}$, where $m^{\prime}$ is the reduced mass defined by $1/m^{\prime}=1/M+1/m$.

By taking account of this aspect, the quasi-harmonic Hamiltonian valid at *T$\lesssim$*100 K, attributing to coupled vibrations between cages and guest atoms, can be expressed in the vector form as

[TABLE]

where $\hat{\bm{U}}_{\ell}=(\hat{U}_{\ell}^{x},\hat{U}_{\ell}^{y})$ is the unit vector for the vector $\bm{U}_{\ell}$. $\{\phi_{\ell}\}$ and $\{\theta_{\ell}\}$ represent the azimuthal and the polar angle in spherical coordinates. The effect of “random” orientation of guest ions $\{\phi_{\ell}\}$ induced by off-centeredness are involved in $\left\{\textbf{{U}}_{\ell}\right\}$. The relation between off-centeredness and disorder in Eq. (6) is described in details in Supplemental Material (SM).

III.2 Anharmonic coupling between acoustic phonons and SL modes

When acoustic modes (LA and TA) are propagating along networked cages, the cages are distorted and these change the states of guest ions, which are realized via the change of the force constants $\xi_{r}$ and $\xi_{\phi}$ in Eq. (6). The in-plane (stretching and libration) modes are sensitive to temperature/pressure compared with out-of-plane modes as shown in the optic spectroscopy data below *T$\simeq$*100 K .Mori2011 ; Takasu2006 Thus, the anharmonic effect between acoustic modes and in-plane modes in the first and the second terms in Eq. (6) becomes relevant in comparison with the third term. The expansions of $\xi_{r}$ and $\xi_{\phi}$ with respect to the strain tensor $e_{\alpha\beta}$ for $\alpha,\beta=x,y,z$ provide

[TABLE]

Here the coefficients are defined by $D_{r,(\phi)}=\partial\xi_{r,(\phi)}/\partial e_{\alpha\alpha}$, $S_{r,(\phi)}=\partial\xi_{r,(\phi)}/\partial e_{\alpha\beta(\alpha\neq\beta)}$ where $e_{\alpha\beta}=1/2\left(\partial u_{\alpha}/\partial x_{\beta}+\partial u_{\beta}/\partial x_{\alpha}\right)$ is the component of strain tensor. It should be noted that $e_{\alpha\alpha}$ expresses the compression or expansion, and $e_{\alpha\beta(\alpha\neq\beta)}$ does the shear destorsion. The expansion in Eq. (7) leads to the following anharmonic interaction expressed in the vector form as

[TABLE]

Here we note that Eq. (8) satisfies the condition of infinitesimal translational invariance as a whole; $V^{\prime}_{CG}\rightarrow 0$ under the long wavelength limit $k_{\mu}\rightarrow 0$. We emphasize again that Eq. (8) is valid at temperatures *T$\lesssim$*100 K where the guest atoms execute coupled vibrations with cages.Mori2011 ; Takasu2006 While, at *T$\gtrsim$*100 K, $\kappa_{L}$(T) saturates without exhibiting the appreciable T-dependence, where guest atoms behave like rattlers in cages termed by the ”rattling” motion, where the concept of vibrational modes is invalid.Mori2011 ; Takasu2006

IV The 2nd quantized form of interaction Hamiltonian

IV.1 Acoustic phonons causing from networked cages

Provided that EX acoustic phonons with wavelengths $\lambda$ much larger than the lattice spacing $a_{0}$ propagate through networked cages, the molecular units and guest ions vibrate “in phase”. The displacement at the site $\ell$ is expressed by the sum of plane waves as given by

[TABLE]

Here the symbols $b{{}_{\textbf{{k}}_{\mu}}^{\dagger}}\,(b_{\textbf{{k}}_{\mu}})$ express the creation (annihilation) operator for acoustic phonon of the mode $(\textbf{{k}}_{\mu})$ with $\mu=\parallel,\perp$, which represent longitudinal and transverse modes, respectively. The vector $\bm{R}_{\ell}$ expresses the equilibrium position of the $\ell$th molecular unit as depicted in Fig. 4, and $h.\,c.$ indicates the Hermitian conjugate. The mass density is defined as $\rho=6(m+M)/a_{0}^{3}$ with the size of unit cell of $a_{0}$ since 6 molecuar units are involved in unit cell of type-I clathrates. See Sec. I in Supplemental Material (SM) about the definitions employed in this paper.

The function $\phi_{\textbf{{k}}_{\mu}}(\textbf{{R}}_{\ell})$ in Eq. (9) takes the form of

[TABLE]

The normalization condition for $\phi_{\textbf{{k}}_{\mu}}(\textbf{{R}}_{\ell})$ is given by

[TABLE]

IV.2 Strongly localized modes due to guest ions

Figure 3 provides the mode belonging to the eigenenergy $\varepsilon_{q}$=2.6 meV obtained for the system size 20$\times$20$\times$20. This mode pattern indicates that the localization length $L_{\lambda}$ is comparable with the wavelength $2\pi/k_{\lambda}$, $i.e.$, localized within several molecular units, manifesting the Ioffe-Regel condition of the strong localization. On the basis of these numerical findings, we can express the form of SL modes in terms of the relative coordinate $\textbf{{w}}_{\ell}(t)=\textbf{{u}}_{l}(t)-\textbf{{r}}_{l}(t)$ as

[TABLE]

Here the mass $m^{\prime}$ is the reduced mass defined by $1/m^{\prime}=1/M+1/m$, where $M$ is the mass of the molecular unit given in Fig. 1, much larger than the mass of guest ion $m$, for example, $M=6.01m$ for off-center type-I BGS. The symbol $c{{}_{\lambda}^{\dagger}}\,(c_{\lambda})$ represents the creation (annihilation) operator for the localized mode $\lambda$. We put forward the Ansatz for the amplitude $\psi_{\lambda}(\textbf{{R}}_{\ell})$ of the form

[TABLE]

where $\textbf{{R}}_{\lambda}$ represents the center of SL mode $\lambda$. This wave function has vanishing group-velocities $v_{g}$ characterizing localized modes.

The prefactor $A$ in Eq. (13) can be determined from the normalization condition of

[TABLE]

where $\Omega=V/N$ is the volume of the molecular unit depicted in Fig. 1(d). This yields, by combing with the Ioffe-Regel condition,

[TABLE]

The above has been obtained by using the formula $\mathop{\mathrm{cos}}\nolimits^{2}(\textbf{{k}}\cdot\textbf{{R}})=(\mathop{\mathrm{cos}}\nolimits(2\textbf{{k}}\cdot\textbf{{R}})+1)/2$. According to the Ioffe-Regel condition $k\approx 2\pi/L_{\lambda}$, the 1st term in the integral becomes negligible compared with the 2nd term since the 1st term yields rapidly oscillating function in the integrand. This leads to Eq. (15). Thus, the normalized wave function of the SL mode $\lambda$ becomes

[TABLE]

IV.3 Anharmonic Hamiltonian between SL and EX modes

We consider here the effect of incoming EX acoustic phonons with the polarization vector $\hat{\textbf{{e}}}_{\textbf{{k}}_{\mu}}$ to SL modes with the polarization vectors $\hat{\textbf{{e}}}_{\lambda^{\prime}}$ and $\hat{\textbf{{e}}}_{\lambda^{\prime\prime}}$. These are included in Eq. (8) as the scalar product $(\hat{\textbf{{e}}}_{\lambda^{\prime}}\cdot\hat{\textbf{{U}}}_{\ell})(\hat{\textbf{{e}}}_{\lambda^{\prime\prime}}\cdot\hat{\textbf{{U}}}_{\ell})$ and the product $(\hat{\textbf{{e}}}_{\lambda^{\prime}}\times\hat{\textbf{{U}}}_{\ell})\cdot(\hat{\textbf{{e}}}_{\lambda^{\prime\prime}}\times\hat{\textbf{{U}}}_{\ell})$. At first, we fix the direction of the wave vector of incoming EX phonons $\textbf{{k}}_{\mu}$ and later we include the contributions from 3 components of the wave vector $\textbf{{k}}_{\mu}$. We should note that the deformation (normal or shear strain) of cages causing from incoming acoustic phonons responses to every directions of the polarization vector of SL modes, which provides both the interaction between the same polarization and different polarizations of SL modes as shown below.

The second quantized anharmonic Hamiltonian is obtained by substituting Eqs. (9) and (12) into Eq. (8) by using the relations given in Sec. II in SM. The product of the field operators $b_{\textbf{{k}}_{\mu}}c_{\lambda^{\prime}}c_{\lambda^{\prime\prime}}$ consists of eight terms. The two involve the combinations $b_{\textbf{{k}}_{\mu}}^{\dagger}c_{\lambda^{\prime}}^{\dagger}c_{\lambda^{\prime\prime}}^{\dagger}$ and $b_{\textbf{{k}}_{\mu}}c_{\lambda^{\prime}}c_{\lambda^{\prime\prime}}$ are irrelevant to the hooping processes because of not conserving the total energy. Furthermore the other two terms $b_{\textbf{{k}}_{\mu}}^{\dagger}c_{\lambda^{\prime}}c_{\lambda^{\prime\prime}}$ and $b_{\textbf{{k}}_{\mu}}c_{\lambda^{\prime}}^{\dagger}c_{\lambda^{\prime\prime}}^{\dagger}$ do not contribute to the scattering processes since the energies of EX modes are smaller than those of SL modes. Hence, the relevant second quantized anharmonic Hamiltonian for the process on EX + SL $\rightarrow$ SL is given by

[TABLE]

where $A_{\textbf{{k}}_{\mu},\lambda^{\prime},\lambda^{\prime\prime}}$ is associated with the interaction between the modes with $x$-polarization, $B_{\textbf{{k}}_{\mu},\lambda^{\prime\prime\prime},\lambda^{\prime\prime\prime\prime}}$ corresponds to the interaction between $y$-polarization, and $C_{\textbf{{k}}_{\mu},\lambda^{\prime},\lambda^{\prime\prime\prime}}$ does the interaction between two different polarizations. See Fig. 5.

By taking the unit vectors $\hat{x},\hat{y},\hat{z}$ the same as the directions of the polarizations $\hat{\textbf{{e}}}_{\parallel},\hat{\textbf{{e}}}_{\perp},\hat{\textbf{{e}}}_{\perp^{\prime}}$ of EX acoustic modes, we have

[TABLE]

and the term on $B_{\textbf{{k}}_{\mu},\lambda^{\prime\prime\prime},\lambda^{\prime\prime\prime\prime}}$ becomes the same as $A_{\textbf{{k}}_{\mu},\lambda^{\prime},\lambda^{\prime\prime}}$ by setting [$\lambda^{\prime},\lambda^{\prime\prime}\longrightarrow\lambda^{\prime\prime\prime},\lambda^{\prime\prime\prime\prime}$]. The last one should be

[TABLE]

The squared quantity on Eq. (18) is given by

[TABLE]

where the coefficient $C$ is defined as

[TABLE]

The expression of $B_{\textbf{{k}}_{\mu},\lambda^{\prime\prime\prime},\lambda^{\prime\prime\prime\prime}}^{2}$ takes the same form as $A^{2}_{\bm{k}_{\mu},\lambda^{\prime},\lambda^{\prime\prime}}$, since they both correspond to the interaction between SL modes with the same polarization. While $C_{\textbf{{k}}_{\mu},\lambda^{\prime},\lambda^{\prime\prime\prime}}^{2}$ corresponding to interaction between different polarizations has an additional factor $(4/\pi)^{2}$ and $\left[(D_{r}-D_{\phi})k_{\parallel}\delta_{\mu,\parallel}+(S_{r}-S_{\phi})k_{\perp}\delta_{\mu,\perp}\right]^{2}$.

V Hopping Process

V.1 Relaxation time of SL modes

This subsection gives the formula for the relaxation time of SL mode due to the scattering process EX+SL$\rightarrow$ SL (hopping process) together with its reverse process shown in Fig. 5 by applying the Fermi golden rule. To obtain the total transition rate of the SL mode in $\lambda^{\prime}$, we have to incorporate all of four processes for each polarization as given below. These provide the decay of the Bose-Einstein distribution function $n_{\lambda^{\prime}}$ for the occupied state $\lambda^{\prime}$,

[TABLE]

We consider, at first, the decay due to the hopping process between the same polarization, $i.e.$, the contribution from the first two terms of Eq. (22). By separating the distribution function into two parts. $n=n^{(0)}+n^{(1)}$, where $n^{(0)}$ is the Bose-Einstein distribution function in equilibrium state and $n^{(1)}$ is its deviation due to the scattering processes, and by employing the relaxation time approximation, $dn{{}_{\lambda^{\prime}}^{(1)}}/dt=-n{{}_{\lambda^{\prime}}}^{(1)}/\tau_{\lambda^{\prime}}$, we have the inverse of relaxation time from Eq. (22) for the same polarization process,

[TABLE]

where the explicit form of the summation $I_{1}$ arising from the overlapping of wave functions $\psi_{\lambda^{\prime}}$ and $\psi_{\lambda^{\prime\prime}}$ is given by

[TABLE]

The above sum $I_{1}$ can be reduced to, by taking the origin of the sum as $\textbf{{R}}_{\lambda^{\prime}}=0$ and the nearest neighbor position from the origin as $\textbf{{R}}_{\lambda^{\prime\prime}}=\Delta\textbf{{R}}_{\lambda^{\prime\prime}}$,

[TABLE]

where the even function $f(\textbf{{X}}_{\ell})$ is defined as

[TABLE]

Since the localization lengths of SL modes are the same, $e.g.$, $L_{\lambda^{\prime}}\cong L_{\lambda^{\prime\prime}}$, hereafter we denote this as $L$. As $f(\textbf{{X}}_{\ell})$ is concerned with SL modes, the relevant sum should be made in the region $\mid\textbf{{X}}_{\ell}\mid\leq L$, so we can approximate the summation by

[TABLE]

where we have used the approximation $\mathop{\mathrm{cos}}\nolimits(\textbf{{k}}_{\lambda^{\prime}}\cdot\Delta\textbf{{R}}_{\lambda^{\prime\prime}})\approx\mathop{\mathrm{cos}}\nolimits(k_{\lambda^{\prime}}nL)\approx 1$ from the Ioffe-Regel condition $L\approx 2\pi/k_{\lambda^{\prime}}$ for SL modes and $e^{-i\textbf{{k}}_{\mu}\cdot\textbf{{R}}_{\ell}}\approx 1$ due to $\mid\textbf{{k}}_{\mu}\mid\ll 2\pi/L$ for the wave number of EX acoustic modes. The term containing $\mathop{\mathrm{cos}}\nolimits\left(2\textbf{{k}}_{\lambda^{\prime}}\cdot\textbf{{R}}_{\ell}-\textbf{{k}}_{\lambda^{\prime}}\cdot\Delta\textbf{{R}}_{\lambda^{\prime\prime}}\right)$ becomes negligible since it yields rapidly oscillating function in the integrand.

This gives the squared hopping integral of the form

[TABLE]

where $\Delta R_{\lambda^{\prime\prime}}$ is the hopping distance.

In the temperature regime T$\simeq$ a few 10 K, $i.e.$, $k_{B}T>\hbar\omega_{\lambda^{\prime}},\hbar\omega_{\lambda^{\prime\prime}}>\hbar\omega_{\textbf{{k}}\mu}$, the inverse of the relaxation time takes the following form under the above conditions and by employing the linear dispersion relation for EX phonon mode $\omega_{\textbf{{k}}\mu}=v_{\mu}k_{\mu}$,

[TABLE]

Here the coefficient $C$ is defined in Eq. (21). We have omitted the temperature independent term providing only small contributions.

V.2 Thermal conductivity due to the hopping of SL modes

In the previous subsection, we have formulated the relaxation rate of SL modes due to the anharmonic interaction between SL modes and EX modes. This is a quantum process realizing the decay of SL*′* mode to SL*′′* mode assisted by EX mode: SL*′+EX$\rightarrow$SL′′*. Without anharmonic interaction, SL modes cannot diffuse/contribute to thermal transport. This means that the plateau region should continue over at higher temperatures after exhibiting the onset of the plateau, $i.e.$, the contribution from EX modes to lattice thermal conductivity is saturated at higher temperatures. This is because the onset of the plateau arises from the weak localization of acoustic modes as explained in Sec.II. Thus, the T-linear rise of $\kappa_{\rm L}(T)$ cannot recover without anharmonic interaction between SL modes and EX modes.

In addition, we emphasize that disorder, induced by off-centeredness as shown in Supplemental Material, is essential to generate the hopping of SL modes. This occurs only in the case that SL*′* mode belonging to the eigenfrequency $\omega_{\rm SL^{\prime}}$ can hop to a site of SL*′′* mode with a different eigenfrequency $\omega_{\rm SL^{\prime\prime}}$ via absorption or emission of EX mode with finite frequency $\pm(\omega_{\rm SL^{\prime}}-\omega_{\rm SL^{\prime\prime}})$. This finite frequency is created by level repulsion between eigenfrequencies due to disorder, $i.e.$, localized modes never belong to the same eigenfrequency according to the level repulsion.

Let us provide the formula of $\kappa_{\rm L}(T)$ due to the diffusion process where SL modes serve as primary heat carriers. In this process, the characteristic length-scale should be the hopping distance $\Delta R_{\lambda^{\prime\prime}}$ from the site of SL*′* mode to that of SL*′′* mode, and the characteristic time-scale is the relaxation time $\tau_{\lambda^{\prime}}$ of the SL*′* mode. This leads to the following formula of the lattice thermal conductivity due to the hopping process, which was first proposed for fracton excitations by Alexander $et.~{}al.$ Alexander1987 ,

[TABLE]

where $\Delta R_{\lambda^{\prime\prime}}^{2}/\tau_{\lambda^{\prime}}$ is the thermal diffusivity of SL mode $\lambda^{\prime}$, $C_{\lambda^{\prime}}(T)$ is the specific heat associated with the SL mode $\lambda^{\prime}$. In the high temperature regime above the plateau region *T$\gtrsim$*a few 10 K, the specific heat follows the Dulong-Petit relation of the form $C_{\lambda^{\prime}}(T)=k_{\mathrm{B}}$ per one polarization of SL mode $\lambda^{\prime}$. Note that $1/\tau_{\lambda^{\prime}}=1/\tau_{\lambda^{\prime}}^{\rm same}+1/\tau_{\lambda^{\prime}}^{\rm dif}$, we first calculate the hopping process between the same polarization by,

[TABLE]

The substitution of Eq. (29) into Eq. (31) together with Eq. (28) yields

[TABLE]

Transforming the sum $\sum_{\textbf{{k}}_{\mu}}$ for EX phonon modes to the integral $V/(2\pi)^{3}\int d\textbf{{k}}_{\mu}=V/(2\pi^{2}v_{\mu}^{3})\int\omega_{k_{\mu}}^{2}d\omega_{k_{\mu}}$, we have

[TABLE]

The sum on $\lambda^{\prime}$ and $\lambda^{\prime\prime}$ above should include the density of states of SL modes $D_{\rm SL}(\omega_{\lambda^{\prime}})$ and $D_{\rm SL}(\omega_{\lambda^{\prime\prime}}(\Delta R_{\lambda^{\prime\prime}}))$ for the same polarization process. The volume $\Omega$ should contain two independent SL modes corresponding to two independent in-plane mode, say, stretching or libration, in the band width of $\Delta\omega_{\rm sl}$, which leads to

[TABLE]

and

[TABLE]

where the volume $\Omega$ contains at least one possible SL mode $\lambda^{\prime\prime}$ with the same/different polarization as/from mode $\lambda^{\prime}$. Since the term $\Delta R_{\lambda^{\prime\prime}}^{4}e^{-2\Delta R_{\lambda^{\prime\prime}}/L}$ in Eq. (35) achieves its maximum at $\Delta R_{\lambda^{\prime\prime}}=2L$ and it decays fast with the further increasing of $\Delta R_{\lambda^{\prime\prime}}$, the sum of $\lambda^{\prime\prime}$ could be estimated within the sphere region $\Delta R_{\lambda^{\prime\prime}}\leq\Delta R$.

[TABLE]

Here the sum on SL modes are done by $\sum_{\lambda^{\prime\prime}}=4\pi\Delta R^{3}/3\int_{\omega_{\rm sl}}^{\omega_{\rm sl}+\Delta\omega_{\rm sl}}D(\omega_{\lambda^{\prime\prime}}(\Delta R_{\lambda^{\prime\prime}}))d\omega_{\lambda^{\prime\prime}}$ and $\sum_{\lambda^{\prime}}=V\int_{\omega_{\rm sl}}^{\omega_{\rm sl}+\Delta\omega_{\rm sl}}D(\omega_{\lambda^{\prime}})d\omega_{\lambda^{\prime}}$, where the factor $4\pi\Delta R^{3}/3\Omega$ from Eq. (35) means the total number of hopping sites from $\lambda^{\prime}$ to $\lambda^{\prime\prime}$ for the same polarization process, and $2V/\Omega$ from Eq. (34) is the total number of $\lambda^{\prime}$ contributing the thermal conductivity $\kappa_{\rm hop}$. The numerical factor $10^{-2}$ arises from the magnitude estimation of integral $\int_{\omega_{\rm sl}}^{\omega_{\rm sl}+\Delta\omega_{\rm sl}}d\omega_{\lambda^{\prime}}\int_{\omega_{\rm sl}}^{\omega_{\rm sl}+\Delta\omega_{\rm sl}}d\omega_{\lambda^{\prime\prime}}\frac{(\omega_{\lambda^{\prime\prime}}-\omega_{\lambda^{\prime}})^{2}}{\Delta\omega_{\rm sl}^{2}\omega_{\lambda^{\prime\prime}}^{2}}$.

The formula of the thermal conductivity due to the hopping mechanism is given by

[TABLE]

The same procedure for the hopping process due to anharmonic interaction between different polarizations leads to

[TABLE]

The total thermal conductivity due to the hopping mechanism is given by the sum of these components as

[TABLE]

V.3 Evaluation of anharmonic coupling $D$ and $S$

Here we estimate the anharmonic coupling constants $D_{r(\phi)}$ and $S_{r(\phi)}$ by illustrating type-I BGS. The coupling constants $D_{r}(S_{r})$ and $D_{\phi}(S_{\phi})$ are associated with the stretching and libration motion of guest-cage vibrations identified by the force constant $\xi_{r}$ and $\xi_{\phi}$ in Eq. (6) by the relation $\xi_{r(\phi)}=m^{\prime}\omega_{r(\phi)}^{2}$, where $m^{\prime}$ is the reduced mass defined by $1/m^{\prime}=1/m+1/M$. In our coarse-grained Hamiltonian introduced in Sec. III, the guest ion Ba(2) in tetrakaidecahedron cage has the mass $m$ and the molecular unit composed of 1 tetrakaidecahedron and 1/3 dodecahedron does the total mass $M$ excluding the off-center guest ion.

We first evaluate the coupling constants $D_{r(\phi)}$ from the Raman spectroscopy data of pressure dependence  Kume2015 . The $D_{r}$ can be related to the pressure $P$ by

[TABLE]

Here $B=\frac{\Delta P}{(\Delta V/V)}$ is the linear thermal expansion coefficient, where the dilation is given by $\Delta V/V=\sum_{\alpha}e_{\alpha\alpha}$ for cubic structure. The coupling constant $D_{\phi}$ can be defined in a similar manner to Eq. (40) as

[TABLE]

In the pressure range from 0.8 GPa to 5.8 GPa, Eg mode spans from  20 cm*-1* to  27 cm*-1*. While, for T2g mode, it ranges from  17 cm*-1* to  27 cm*-1*. The observed spectra of these two modes are overlapped/mixed. Taking account of these aspects, we have $\partial\omega^{r}/\partial P=2\pi\times 4.2\times 10^{10}\,\rm{[sec^{-1}GPa^{-1}]}$ and $\partial\omega^{\phi}/\partial P=2\pi\times 6.0\times 10^{10}\,\rm{[sec^{-1}GPa^{-1}]}$. We then obtain the coupling constants $D_{r}=m^{\prime}\pi^{2}\times 3.0\times 10^{25}\,\rm{[kg\cdot sec^{-2}]}$ and $D_{\phi}=m^{\prime}\pi^{2}\times 3.0\times 10^{25}\,\rm{[kg\cdot sec^{-2}]}$ using the observed bulk modulus $B=41.3\rm{GPa}$ Ishii2012 . Within our knowledge, the experiment data for estimating the coupling coefficients $S_{r(\phi)}$ are not available, so we assume as $S_{r}\approx D_{r}$ and $S_{\phi}\approx D_{\phi}$ at the present stage. The above coupling constants yield

[TABLE]

where we have employed the values of parameters in Eq. (39) as the localization length $L=2a_{0}$, the hopping distance $\Delta R=3.5L$, the volume of molecular unit $\Omega=(a_{0})^{3}/6$, the lattice spacing $a_{0}=11.68$Å, the mass density $\rho=6.01\times 10^{3}$ kg/m3, in addition to the velocities of acoustic phonons $v_{\parallel}=3369$ m/s and $v_{\perp}=1936$ m/s Avila:2006a . The value of $\kappa_{\rm hop}$ in Eq. (42) is smaller than the observed one of $\kappa_{\rm hop}=9.2\times 10^{-3}T(\mathrm{Wm^{-1}K^{-1}})$ for type-I BGS. This mainly arises from, as will be demonstrated below by means of FPC, the underestimated shear coupling constants $S_{r(\phi)}$ obtained by assuming the relations $S_{r(\phi)}\approx D_{r(\phi)}$.

Due to the lack of experiment data for the shear coupling coefficients $S_{r(\phi)}$, we have performed FPC for type-I BGS to obtain the coupling constants from the shift of eigenfrequencies at $\Gamma$-point of low-lying optical mode by imposing strain to the cage structure. The normal strain is isotropic and defined as $e_{\alpha\alpha}=(a_{0}-a)/a_{0}$ where $a_{0}$ and $a$ are the lattice constant for the unstrained and strained unit cell Chen2014 , respectively. The shear strain is also isotropic and defined as $e_{\alpha\beta}=(1-\sqrt{1-(2\mathop{\mathrm{cos}}\nolimits\theta-1)\mathop{\mathrm{cos}}\nolimits\theta})/(2\mathop{\mathrm{cos}}\nolimits\theta-1)$ where $\theta$ is the acute angle between edges after deformation.

We have performed the FPC by the VASP code Kresse1999 with the Perdew-Burke-Ernzerhof functional and the PAW method Perdew1996 , plane wave cut-off energy 250 eV and the force convergence less than $10^{-7}\rm{eV/\AA}$. The phonon frequencies are calculated by PHONOPY code Togo2015 with the $4\times 4\times 4$ Monkhorst-Pack $k$ grids and for a unit cell containing 54 atoms. The coupling constants obtained from normal strain are $D_{r}=m^{\prime}\pi^{2}\times 2.1\times 10^{25}\,\rm{[kg\cdot sec^{-2}]}$, $D_{\phi}=m^{\prime}\pi^{2}\times 1.5\times 10^{25}\,\rm{[kg\cdot sec^{-2}]}$, and from sheared unit cell are $S_{r}=m^{\prime}\pi^{2}\times 4.2\times 10^{25}\,\rm{[kg\cdot sec^{-2}]}$, $S_{\phi}=m^{\prime}\pi^{2}\times 2.9\times 10^{25}\,\rm{[kg\cdot sec^{-2}]}$, respectively. The $D_{r(\phi)}$ are smaller than those estimated from the Raman spectroscopy data of pressure dependence, though $S_{r(\phi)}$ are larger than the values obtained from the assumption $S_{r(\phi)}\approx D_{r(\phi)}$. The above coupling constants yield the thermal conductivity due to the hopping of SL modes of

[TABLE]

We remark here that our FPC provides the results for the on-center positioned Ba(2) because the optimization for off-center structure is quite time-consuming and may require to take into account the dipole-dipole interaction due to off-centeredness and temperature effect. The on-center structure gives rise to the underestimated coupling constants $S$ since on-center guest ions should more weakly response to shear distortion than the case of off-center. Then, the actual $S_{r(\phi)}$ should be larger than the above estimation. Under these situations, the calculated value in Eq. (43) provides reasonable agreement, to claim the relevance of the hopping process of SL modes, with the observed $\kappa_{\rm hop}=\gamma T$ with $\gamma=9.2\times 10^{-3}\,\mathrm{Wm^{-1}K^{-2}}$ for type-I BGS Avila:2006a ; Avila:2008aa , and $\gamma=9.0\times 10^{-3}\,\mathrm{Wm^{-1}K^{-2}}$ for type-I EGG Sales:2001a . For type-I SGG, several different values around $\gamma\sim 8.0\times 10^{-3}\,\mathrm{Wm^{-1}K^{-2}}$ have been reported Nolas:1998aa ; Cohn:1999a ; Christensen:2016a ; Suekuni:2007a , indicating that the experimental data of SGG depend on sample qualities according to synthesis methods. In that respect, it has been reported Christensen:2016a that a flux-grown sample shows a glasslike plateau, while a zone-melted sample has a crystalline peak.

VI Summary and Conclusions

Off-center type-I clathrates show almost identical lattice thermal conductivities $\kappa_{\rm L}$ to those of structural glasses  Nolas:1998aa ; Cohn:1999a ; Sales:2001a ; Paschen:2001a ; Suekuni:2007a ; Avila:2008aa ; Suekuni:2008a . In addition, off-center type-I clathrates show the excess density of states at THz frequencies manifesting the boson peak identical to those of network-forming glasses Avila:2006a ; Avila:2008aa ; Suekuni:2008a . These indicate that the symmetry broken guest ions in cages take charge of the emergence of glasslike $\kappa_{\rm L}(T)$. In structural glasses, many key aspects of a detailed quantitative description are still missing. This is due to the difficulty to identify relevant entities or elements at atomic scale caused by their complex microscopic structures.

In Sec. II, we have pointed out that the PR shown in Fig. 2 provides the evidence that EX acoustic phonons carrying heat convert to WL modes modes at $\sim$ 1.3 meV in off-center BGS. This energy corresponds to the temperature 3.9 K $\approx$1.3 meV/3.84kB from the Wien’s displacement law, so that this conversion should be associated with the onset of the plateau thermal conductivities observed at several K in off-center type-I clathrates Nolas:1998aa ; Cohn:1999a ; Sales:2001a ; Paschen:2001a ; Avila:2006a ; Suekuni:2007a ; Avila:2008aa ; Suekuni:2008a .

With increasing temperature further, thermal conductivities above a few 10 K show a linear rise on temperature. This type of anomalous thermal conductivities with the plateau and the subsequent T-linear rise have been clearly observed for off-center type-I clathrates Nolas:1998aa ; Cohn:1999a ; Sales:2001a ; Paschen:2001a ; Avila:2006a ; Suekuni:2007a ; Avila:2008aa ; Suekuni:2008a . This is the prominent hallmark of glasslike thermal conductivity since crystals with translational invariance never show these features. Rather, lattice thermal conductivities of crystallines decrease with increasing temperature proportional to $\kappa(T)\propto 1/T$ known as the Umklapp process Landau1979 .

The theoretical elucidation on the linear rise on temperature “above” the plateau region has been the main subject of the present paper. Our calculated results given in Sec. V, based on hopping process, show fairly good agreement with observed thermal conductivities above the plateau. We particularly emphasize that both the magnitude and the temperature dependence of $\kappa(T)$ are in accord with the experimental data Nolas:1998aa ; Cohn:1999a ; Sales:2001a ; Paschen:2001a ; Avila:2006a ; Suekuni:2007a ; Avila:2008aa ; Suekuni:2008a . At much higher temperatures, the T-linear rise in $\kappa(T)$ does not continue, but $\kappa(T)$ saturates above $T\simeq$ 100 K Sales:2001a ; Avila:2006a ; Suekuni:2007a ; Avila:2008aa . In this temperature regime, the treatment based on quantum mechanical process does not hold for since the life-time of excited modes becomes much smaller than the inverse of their angular frequencies, where the guest ions become free from the constraint of atoms constituting cages. This subject will be discussed in detail elsewhere Nakayama2017 .

In conclusion, the phenomenon of T-linear rise of $\kappa_{L}(T)$ above a few 10K in off-center type-I clathrates has been quantitatively explained by analytic theory, on the grounds that off-center clathrates possess definite microscopic structure. Our successful clarification in quantitative manner is owing to the fact that the systems are more tractable than network-forming glasses with the difficulty to identify relevant constituents at atomistic level caused by their complex microscopic structures.

Acknowledgments.  This work is supported by the National Natural Science Foundation of China Grant No. 11334007 and No. 51506153. J. Z. is supported by the program for Professor of Special Appointment (Eastern Scholar) at Shanghai Institutions of Higher Learning No. TP2014012. T. N. acknowledges the support from Grand-in-Aid for Scientific Research from the MEXT in Japan, Grand No.26400381.