Hopping Processes Explain T-linear Rise of Thermal Conductivity in Thermoelectric Clathrates above the Plateau
Qing Xi, Zhongwei Zhang, Jie Chen, Jun Zhou, Tsuneyoshi Nakayama, Baowen Li
aa r X i v : . [ c ond - m a t . o t h e r] J u l Hopping Processes Explain T -linear Rise of Thermal Conductivity in Thermoelectric Clathratesabove the Plateau Qing Xi,
1, 2, 3
Zhongwei Zhang,
1, 2, 3
Jie Chen,
1, 2, 3
Jun Zhou, ∗ Tsuneyoshi Nakayama, † and Baowen Li ‡ Center for Phononics and Thermal Energy Science, School of PhysicsScience and Engineering, Tongji University, 200092 Shanghai, P. R. China China-EU Joint Center for Nanophononics, School of Physics Science and Engineering, Tongji University, 200092 Shanghai, P. R. China Shanghai Key Laboratory of Special Artificial Microstructure Materials and Technology,School of Physics Science and Engineering, Tongji University, 200092 Shanghai, P. R. China Hokkaido University, 060-0826 Sapporo, Japan Department of Mechanical Engineering, University of Colorado, Boulder, Colorado 80309, USA (Dated: September 27, 2018)Type-I clathrate compounds with off-center guest ions realize the phonon-glass electron-crystal concept byexhibiting almost identical lattice thermal conductivities κ L to those observed in network-forming glasses. Thisis in contrast with type-I clathrates with on-center guest ions showing κ L of conventional crystallines. Glasslike κ L stems from the peculiar THz frequency dynamics in off-center type-I clathrates where there exist three kindsof modes classified into extended(EX), weakly(WL) and strongly localized(SL) modes as demonstrated by Liu et. al. , Phys. Rev. B , 214305(2016). Our calculated results based on the hopping mechanism of SL modesvia anharmonic interactions show fairly good agreement with observed T -linear rise of κ L above the plateau.We emphasize that both the magnitude and the temperature dependence are in accord with the experimental dataof off-center type-I clathrates. PACS numbers: 63.20.Pw Localized modes 63.20.Ry Anharmonic lattice modes 63.50.+x Vibrational states in disorderedsystems
I. INTRODUCTION
Lattice thermal conductivity constitutes a key elementto improve the efficiency of the thermal-to-electrical con-version in thermoelectric (TE) devices as understood fromthe material’s figure of merit describing the efficiency Z = S σ/κ tot [K − ]. The numerator contains the See-beck coefficient S ( T ) [V/K] and the electrical conductiv-ity σ ( T ) [1/( Ω m)], while the denominator κ tot ( T ) [W/(mK)]consists of the sum of electrical κ el and lattice κ L thermal con-ductivity. Hence, the high performance of thermoelectricitycan be achieved for materials with the lowest possible thermalconductivity κ tot , the highest possible electrical conductivity σ and the highest possible Seebeck coefficient S . Providedthat the Wiedemann-Franz law κ el ( T ) ∝ σ ( T ) holds for, κ L becomes a crucial parameter to improve the performance ofTE conversion. In this framework, Slack has proposed theconcept of “phonon-glass electron-crystal”. This has been oneof guiding principles for exploring high-performance TE ma-terials .Type-I clathrates with “off-center” guest ions, suchas R Ga Ge (R=Ba, Sr, Eu) , Ba Ga Sn ,Sr Ga Si − x Ge x , are particularly interesting in thisrespect since these systems exhibit almost identical latticethermal conductivities to those of structural glasses, whichconsist of four specific regions characterized by: (i) T ∼ -dependence below a few Kelvin, (ii) the plateau regionbetween a few K and a few 10K, and (iii) the subsequentrise proportional to T , and (iv) its saturation above T ∼ κ L exhibit a remarkable uniformitywhich appears to be insensitive to chemical compositions,suggesting the existence of a unified mechanism . However,this issue remains as an open and challenging problem of long-standing due to the difficulty to identify relevant entitiesor elements at atomistic level caused by their complex micro-scopic structures. Surprisingly enough, though “off-center”clathrates are crystalline with regularly network structure,the temperature dependence as well as the magnitudes oftheir thermal conductivities are almost identical to thoseof structural glasses over the full temperature range. Incontrast, type-I clathrates with “on-center” guest ions showconventional crystalline κ L .This paper is organized as follows. Section II surveys thecharacteristics of vibrational modes according to the results ofthe spectral density of states, eigenvalues and their eigenvec-tors . We claim in this Section that the onset of the plateauis due to the delocalization-localization (weak localization)transition of acoustic modes. In addition, we point out thatthe temperature region showing the T -linear rise subsequentto the plateau is associated with the energy range where SLmodes are fully excited. Section III describes the constructionof anharmonic interaction Hamiltonian between SL and EXmodes. The second quantized form of anharmonic Hamilto-nian is given in Section IV. Section V develops a theory onthe mechanism governing the T -linear rise of κ L ( T ) above afew 10 K. Excited modes in this temperature region are mostlystrongly-localized (SL) modes satisfying the Ioffe-Regel con-dition as evident from the mode pattern obtained by large-scale numerical simulations . These are hybridized modesbetween acoustic phonons associated with network cages andlocal vibrations of guest ions in cages. Based on these nu-merical evidences, we explain in quantitative manner κ L ( T ) proportional to T observed above the plateau, by introducingthe quantum mechanical process of hopping of SL modes dueto anharmonic interactions, first proposed for fracton excita-tions . Summary and conclusions are given in Sec. VI. a a a AB C (2) (1) (a) (b) (c) (d) BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB CC A x y FIG. 1. (Color online) (a) Illustration of type-I clathrate. The four-fold inversion axes are directed along the axises x, y, z . Red and blueballs represent off-center guest ions in tetrakaidecahedron cages andcentered guest ions in dodecahedron cages, respectively. (b) Two off-center guest ions along y -axis are depicted. (c) The configuration ofeight nearest neighbor guest ions connected by equilateral triangle.The sites A, B, and C in (c) are seated on the chains parallel to x, y, and z , respectively: A=(a/4, 0, a/4), B = (0, a/4, 3a/4), C =(a/2, a/2,a/2). (d) The molecular unit composed of tetrakaidecahedron cagewith off-center guest ion(2) at 24k site and smaller dodecahedral cagewith guest ion(1) at 2a site. II. CHARACTERISTICS OF EXCITED PHONONS AT THZFREQUENCY REGION
Type-I clathrates form a primitive cubic structure (
P m ¯3 n )consisting of 6 tetrakaidecahedron (14-hedrons) and 2 dodec-ahedron (12-hedrons) per unit cell, in which the group-I or -IIelements in the periodic table are encaged in the polyhedronsas guest ions. See Fig. 1. The THz frequency phonon dy-namics of off-center type-I clathrates has been investigated interms of large-scale numerical simulations. They have illus-trated type-I Ba Ga Sn (BGS) exhibiting glasslike κ L ( T ) as a prototype material with off-center guest ions, in which theguest ion Ba(2) in tetrakaidecahedron cage has the mass m and the molecular unit composed of one tetrakaidecahedronand 1/3 dodecahedron does the total mass M excluding theoff-center guest ion. The coarse-grained picture, an operationof reducing the degrees of freedom of the original system, isvalid for our purpose from the following reasons. First of all,EX acoustic modes at THz frequencies play a dominant rolein heat transport since optical modes concerning to the vibra-tions of cages themselves do not contribute to thermal con-ductivity. Second, the wave-length λ of phonons in the fre-quency regime ν ≤ E ≤
10 meV) becomes λ ≥ a ≃ λ = v/ν us-ing the sound velocity v ≈ × [m/sec]. These validatethe coarse-grained Hamiltonian for describing THz frequencydynamics rather than treating all microscopic constituents asequally relevant degrees of freedom. Extremely large system-sizes are required in computer sim-ulations on disorder systems in order to distinguish local-ized modes from extended modes. However, the presentstatus of first-principles calculations (FPC) are limited toinsufficient system-sizes for properly incorporating the dis-order 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 dis-order reproducing glasslike thermal conductivities, but also toexclude finite size effect for propagating acoustic phonons.Liu et. al. have performed calculations for 3D systemsof (20 × × ∼ (100 × × They have also studied the lo-calization nature of exited modes by taking the participa-tion ratio (PR) as a criterion. The PR of a relevant mode { ϕ ℓ ( ε q ); ℓ = 1 , , ...N } belonging to the eigenenergy ε q isdefined by P ( ε q ) = (cid:16)P Nℓ =1 | ϕ ℓ ( ε q ) | (cid:17) N P Nℓ =1 | ϕ ℓ ( ε q ) | , (1)where ℓ denotes the ℓ -th molecular unit depicted in Fig. (d)and N is the total mode number. For EX modes in a finitesystem, P ( ε q ) take values close to ≈ ε q = 0 , and P ( ε q ) becomes ≈ /N for SL modes . Figure 2 (a) is thecalculated phonon density of states (DOS), and (b) the resultsof P ( ε q ) for the size of 20 × ×
20 lattice of off-center type-IBGS. It is remarkable that P ( ε q ) ranges from a value of SLmodes P ( ε q ) ≈ to EX modes of P ( ε q ) ≈ . . We shouldemphasize that there appear three kinds of modes in the THzfrequency region and below classified into EX, WL and SLmodes. SL modes with PR values much smaller than unity arerealized in the energy range from 2 to 3 meV as found fromcalculated mode patterns. Figure 3 depicts the mode patternsof SL mode at ε q =2.6 meV.The calculations of the PR for excited modes depicted inFig. 2 have demonstrated that there exists the delocalization-localization transition at a “finite” frequency ω c distinguish-ing EX and WL modes with the nature of acoustic modesvibrating “in-phase” between guest ions and cages. Further-more, it has been found that WL modes convert to SLmodes at higher frquencies with the nature of optical modesvibrating “out-of-phase” between guest ions and cages. Inthis aspect, we note that Nakayama had demonstrated theclear existence of the transition from WL to SL modes forthe quasi-one-dimensional (1D) coarse-grained model con-sisting of host network and guest atoms connected by randomsprings. It was found that WL modes vibrate in-phase be-tween network atoms and guest atoms, while SL modes man-ifest optical modes vibrating out-of-phase. However, thereis no EX modes due to “quasi-1D” model. This manifeststhe Anderson weak localization criteria where the critical fre-quency ω c takes a finite value in three dimensional (3D) sys-tems, while it vanishes for 1D and 2D systems suggesting noEX modes in 1D and 2D disordered systems. The quasi-1Dmodel should be thought as the simplest theoretical modelfor cage-guest systems with broad implication for the dynam- D O S ε q (meV)(a) 0 0.2 0.4 0 1 2 3 4 5 6 7 D O S ε q (meV)(a) 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3 3.5 4 P R ε q (meV)(a) (b) FIG. 2. (Color online) (a). Calculated phonon density of states(DOS) of off-center type-I BGS for the system size of 99 × × P ( ε q ) defined in Eq. (1) as a function of eigenenergy ε q in theenergy range marking by the blue shadow in (a) for the system sizeof 20 × ×
20 under periodic boundary condition. ics of cage-guest systems.The observed delocalization-localization transition at ε q ≈ κ L in BGS at T P ≈ k B ≈ T P for off-center type-I BGS. We notehere that the random orientation of guest ions in cages plays acrucial role to the localization.With increasing temperature further above a few 10K, κ L show a linear rise on temperature . This typeof anomalous thermal conductivities characterized by theplateau and the subsequent T -linear rise of thermal con-ductivities have been clearly observed for off-center type-I clathrates . SL modes are fully excited abovethe temperature T ≃ ≈ k B from the Wien’s dis-placement law. This indicates that T -linear rise subsequent tothe plateau attributes to the excitations of SL modes. In thefollowing Sections, we present the theoretical interpretationon the underlying mechanism of the linear rise on tempera-ture above the plateau region for κ L . III. COARSE-GRAINED HAMILTONIAN FOR TYPE-IOFF-CENTER CLATHRATESA. Harmonic Hamiltonian
The Hamiltonian for off-center type-I clathrates under acoarse-grained picture consists of the kinetic energy of net-worked cages K C and off-center guest ions in cages K G inaddition to the potential energy of the cage-cage interaction V CC and the cage-guest interaction V CG . This is expressed FIG. 3. (Color online) The mode pattern of SL modes belongingto the eigenenergy ε q =2.6 meV. Both the color scale and cubic sizeindicate the strength of amplitudes at each site. The mode patternis obtained from the system size 20 × ×
20 under fixed boundarycondition. by H = K C + K G + V CC + V CG . (2)The explicit form of the total kinetic energy is given by thesum of K C and K G such as K = 12 X ℓ (cid:0) M ˙ r ℓ ( t ) + m ˙ u ℓ ( t ) (cid:1) , (3)where m and M are masses of the guest ion in tetrakaidecahe-dron cage and the remained molecular unit, respectively. Thevectors r ℓ ( t ) and u ℓ ( t ) represent small displacements of cageand guest ion from their equilibrium positions, R ℓ and R ℓ + U ℓ ,at the site ℓ as depicted in Fig. 4. Note here that guest ionstake random orientation U ℓ ( φ ℓ ) in tetrakaidecahedron cages.The molecular unit composed of tetrakaidecahedron anddodecahedron is elastically connected with neighboring onesby the force constants f k , f ⊥ . These are related to the soundvelocities of longitudinal ( µ = k ) and transverse ( µ = ⊥ )acoustic modes via the relation v µ = a [ f µ / ( m + M )] / with a = a / where a is the lattice spacing of primitive cubicstructure ( P m ¯3 n ) of type-I clathrates. Thus, we can estimatethe force constants f k , f ⊥ from the observed data of soundvelocities. Note here that 6 molecular units are included inunit cell in type-I clathrates. In terms of these quantities, thepotential energy of network cages becomes V CC = X ℓ ′ >ℓ,µ f ℓ,ℓ ′ ,µ r ℓ,µ ( t ) − r ℓ ′ ,µ ( t )) , (4)where µ = k , ⊥ , ⊥ ′ . Hereafter, we keep up to the nearestneighbor coupling ( ℓ ′ = ℓ +1 ) between molecular units, whichare denoted by f k , f ⊥ and f ⊥ ′ . The effect of randomly ori-entated guest ions are included in the following cage-guestinteraction Hamiltonian.The Hamiltonian should satisfy the symmetry of infinites-imal translation-invariance as a whole, i.e. , r ℓ = u ℓ = δ a ,which guarantees acoustic phonons as the Nambu-Goldstoneboson with the eigenfrequency ω k → for k → . Thissymmetry principle also holds for the potential of cage-guestinteraction. Hence, the potential function for the cage-guestinteraction V CG should be given by relative coordinates be-tween the cage and the guest ion of w ℓ ( t ) = u ℓ ( t ) − r ℓ ( t ) ,which is expressed by V CG = X ℓ,m =in , out ξ m w ℓ,m ( t ) , (5)where ξ m represents the force constants between cage andguest ion depending on in-plane (parallel) or out-of-plane mo-tion (perpendicular) to the hexagonal face in the tetrakaideca-hedron cage. The guest ions execute in-plane vibration paral-lel to x − y plane in addition to out-of-plane motions becauseof the anisotropic shape of tetrakaidecahedron cages. Thisis because off-center guest ions are involved in tetrakaidec-ahedron cages whose shape distinguishes the vibrations ofoff-center guest ion(2) in the plane parallel and perpendicu-lar to the hexagonal face of the cage. Mori et al . observedby means of THz time-domain spectroscopy that the lowest-lying peak of off-center BGS at 0.71 THz splits into doublepeaks, ω φ / π =0.5THz and ω r / π =0.72THz for off-centertype-I BGS below T ≃
100 K. These spectra should be assignedto the libration and stretching modes of Ba(2) associated with ξ φ and ξ r . The peak around 1.35 THz is assigned as the out-of-plane motion of Ba(2) to the hexagonal faces of tetrakaidec-ahedron, which should be concerned with ξ θ . The Ramanspectra of off-center Sr Ga Ge (SGG) have observed A stretching mode as 48 cm − , and for off-center Eu Ga Ge (EGG) as 36 cm − at 2 K . Using these data, we can esti-mate the force constants via the relation ξ r, ( φ,θ ) = m ′ ω r, ( φ,θ ) ,where m ′ is the reduced mass defined by /m ′ = 1 /M +1 /m .By taking account of this aspect, the quasi-harmonicHamiltonian valid at T .
100 K, attributing to coupled vibra-tions between cages and guest atoms, can be expressed in thevector form as V CG = 12 X ℓ ξ r ( ˆ U ℓ · w k , ℓ ) + 12 X ℓ ξ φ ( ˆ U ℓ × w k , ℓ ) + 12 X ℓ ξ θ ( w ⊥ , ℓ ) , (6)where ˆ U ℓ = ( ˆ U xℓ , ˆ U yℓ ) is the unit vector for the vector U ℓ . { φ ℓ } and { θ ℓ } represent the azimuthal and the polar angle inspherical coordinates. The effect of “random” orientation ofguest ions { φ ℓ } induced by off-centeredness are involved in { U ℓ } . The relation between off-centeredness and disorder inEq. (6) is described in details in Supplemental Material (SM). B. Anharmonic coupling between acoustic phonons and SLmodes
When acoustic modes (LA and TA) are propagating alongnetworked cages, the cages are distorted and these change thestates of guest ions, which are realized via the change of the R l + r l (t)R l + U l + u l (t) R l+1 + U l+1 + u l+1 (t) RR R l+1 + r l+1 (t)
FIG. 4. (Color online) The definition of the position vectors: R ℓ + r ℓ ( t ) is the position vector of the ℓ -th molecular unit at time t , where R ℓ is the equilibrium position of the ℓ -th cage center and the vector r ℓ ( t ) represents a small displacement from R ℓ at time t . The positionvector of the guest ion(2) is defined by the vector R ℓ + U ℓ + u ℓ ( t ) ,where U ℓ is the equilibrium position of guest ion(2) from R ℓ , and u ℓ ( t ) is a small displacement from R ℓ + U ℓ . force constants ξ r and ξ φ in Eq. (6). The in-plane (stretch-ing and libration) modes are sensitive to temperature/pressurecompared with out-of-plane modes as shown in the optic spec-troscopy data below T ≃
100 K .
Thus, the anharmonic ef-fect between acoustic modes and in-plane modes in the firstand the second terms in Eq. (6) becomes relevant in compar-ison with the third term. The expansions of ξ r and ξ φ withrespect to the strain tensor e αβ for α, β = x, y, z provide ξ r, ( φ ) = ξ (0) r, ( φ ) + X α = x,y,z D r, ( φ ) e αα + X α,β = x,y,zα = β S r, ( φ ) e αβ + ..., (7)Here the coefficients are defined by D r, ( φ ) = ∂ξ r, ( φ ) /∂e αα , S r, ( φ ) = ∂ξ r, ( φ ) /∂e αβ ( α = β ) where e αβ = 1 / ∂u α /∂x β + ∂u β /∂x α ) is the componentof strain tensor. It should be noted that e αα expresses thecompression or expansion, and e αβ ( α = β ) does the sheardestorsion. The expansion in Eq. (7) leads to the followinganharmonic interaction expressed in the vector form as V ′ CG = 12 X ℓ,α = β ( D r e αα + S r e αβ )( ˆ U ℓ · w k ,ℓ ) + 12 X ℓ,α = β ( D φ e αα + S φ e αβ )( ˆ U ℓ × w k ,ℓ ) . (8)Here we note that Eq. (8) satisfies the condition of infinites-imal translational invariance as a whole; V ′ CG → underthe long wavelength limit k µ → . We emphasize againthat Eq. (8) is valid at temperatures T .
100 K where the guestatoms execute coupled vibrations with cages.
While, at T &
100 K, κ L ( T ) saturates without exhibiting the appreciable T -dependence, where guest atoms behave like rattlers in cagestermed by the ”rattling” motion, where the concept of vibra-tional modes is invalid. IV. THE 2ND QUANTIZED FORM OF INTERACTIONHAMILTONIANA. Acoustic phonons causing from networked cages
Provided that EX acoustic phonons with wavelengths λ much larger than the lattice spacing a propagate through net-worked cages, the molecular units and guest ions vibrate “inphase”. The displacement at the site ℓ is expressed by the sumof plane waves as given by r ℓ ( t ) = X k µ s ~ ρω k µ ˆ e k µ (cid:16) φ k µ ( R ℓ ) b † k µ ( t ) + h.c. (cid:17) . (9)Here the symbols b † k µ ( b k µ ) express the creation (annihilation)operator for acoustic phonon of the mode ( k µ ) with µ = k , ⊥ ,which represent longitudinal and transverse modes, respec-tively. The vector R ℓ expresses the equilibrium position ofthe ℓ th molecular unit as depicted in Fig. 4, and h. c. indicatesthe Hermitian conjugate. The mass density is defined as ρ = 6( m + M ) /a with the size of unit cell of a since 6molecuar units are involved in unit cell of type-I clathrates.See Sec. I in Supplemental Material (SM) about the defini-tions employed in this paper.The function φ k µ ( R ℓ ) in Eq. (9) takes the form of φ k µ ( R ℓ ) = r V e i k µ · R ℓ . (10)The normalization condition for φ k µ ( R ℓ ) is given by Z | φ k µ ( R ℓ ) | d R ℓ = 1 . (11) B. Strongly localized modes due to guest ions
Figure 3 provides the mode belonging to the eigenenergy ε q =2.6 meV obtained for the system size 20 × ×
20. Thismode pattern indicates that the localization length L λ is com-parable with the wavelength π/k λ , i.e. , localized within sev-eral molecular units, manifesting the Ioffe-Regel condition ofthe strong localization. On the basis of these numerical find-ings, we can express the form of SL modes in terms of therelative coordinate w ℓ ( t ) = u l ( t ) − r l ( t ) as w ℓ ( t ) = X λ r ~ m ′ ω λ ˆ e λ (cid:16) ψ λ ( R ℓ ) c † λ ( t ) + h.c. (cid:17) . (12)Here the mass m ′ is the reduced mass defined by /m ′ =1 /M + 1 /m , where M is the mass of the molecular unit givenin Fig. 1, much larger than the mass of guest ion m , for ex-ample, M = 6 . m for off-center type-I BGS. The symbol c † λ ( c λ ) represents the creation (annihilation) operator for thelocalized mode λ . We put forward the Ansatz for the ampli-tude ψ λ ( R ℓ ) of the form ψ λ ( R ℓ ) = A cos [ k λ · ( R ℓ − R λ )] e −| R ℓ − R λ | /L λ . (13) FIG. 5. (Color online) The diagrams showing the hopping processfor strongly SL modes arising from anharmonic interaction betweenSL modes and EX modes: (a) SL → EX + SL, and (b) EX + SL → SL. The solid lines denote SL mode and the wavy lines EX mode. where R λ represents the center of SL mode λ . This wavefunction has vanishing group-velocities v g characterizing lo-calized modes.The prefactor A in Eq. (13) can be determined from the nor-malization condition of X ℓ | ψ λ ( R ℓ ) | = 1Ω Z d R ℓ | ψ λ ( R ℓ ) | = 1 , (14)where Ω =
V /N is the volume of the molecular unit depictedin Fig. 1(d). This yields, by combing with the Ioffe-Regel con-dition, A ∼ = s πL λ . (15)The above has been obtained by using the formula cos ( k · R ) = (cos(2 k · R ) + 1) / . According to the Ioffe-Regel con-dition k ≈ π/L λ , the 1st term in the integral becomes negli-gible compared with the 2nd term since the 1st term yieldsrapidly oscillating function in the integrand. This leads toEq. (15). Thus, the normalized wave function of the SL mode λ becomes ψ λ ( R ℓ ) = s πL λ cos [ k λ · ( R ℓ − R λ )] e −| R ℓ − R λ | /L λ . (16) C. Anharmonic Hamiltonian between SL and EX modes
We consider here the effect of incoming EX acousticphonons with the polarization vector ˆ e k µ to SL modes withthe polarization vectors ˆ e λ ′ and ˆ e λ ′′ . These are included inEq. (8) as the scalar product (ˆ e λ ′ · ˆ U ℓ )(ˆ e λ ′′ · ˆ U ℓ ) and the prod-uct (ˆ e λ ′ × ˆ U ℓ ) · (ˆ e λ ′′ × ˆ U ℓ ) . At first, we fix the directionof the wave vector of incoming EX phonons k µ and later weinclude the contributions from 3 components of the wave vec-tor k µ . We should note that the deformation (normal or shearstrain) of cages causing from incoming acoustic phonons re-sponses to every directions of the polarization vector of SLmodes, which provides both the interaction between the samepolarization and different polarizations of SL modes as shownbelow.The second quantized anharmonic Hamiltonian is obtainedby substituting Eqs. (9) and (12) into Eq. (8) by using the rela-tions given in Sec. II in SM. The product of the field operators b k µ c λ ′ c λ ′′ consists of eight terms. The two involve the combi-nations b † k µ c † λ ′ c † λ ′′ and b k µ c λ ′ c λ ′′ are irrelevant to the hoopingprocesses because of not conserving the total energy. Fur-thermore the other two terms b † k µ c λ ′ c λ ′′ and b k µ c † λ ′ c † λ ′′ do notcontribute to the scattering processes since the energies of EXmodes are smaller than those of SL modes. Hence, the rel-evant second quantized anharmonic Hamiltonian for the pro-cess on EX + SL → SL is given by H ′ CG = X k µ ,λ ′ ,λ ′′ ( A k µ ,λ ′ ,λ ′′ b k µ c λ ′ c † λ ′′ + h.c. ) , + X k µ ,λ ′′′ ,λ ′′′′ ( B k µ ,λ ′′′ ,λ ′′′′ b k µ c λ ′′′ c † λ ′′′′ + h.c. )+ X k µ ,λ ′ ,λ ′′′ ( C k µ ,λ ′ ,λ ′′′ b k µ c λ ′ c † λ ′′′ + h.c. ) , (17)where A k µ ,λ ′ ,λ ′′ is associated with the interaction between themodes with x -polarization, B k µ ,λ ′′′ ,λ ′′′′ corresponds to the in-teraction between y -polarization, and C k µ ,λ ′ ,λ ′′′ does the in-teraction between two different polarizations. See Fig. 5.By taking the unit vectors ˆ x, ˆ y, ˆ z the same as the directionsof the polarizations ˆ e k , ˆ e ⊥ , ˆ e ⊥ ′ of EX acoustic modes, we have A k µ ,λ ′ ,λ ′′ = − X l i s ~ ρω k µ r ~ m ′ ω λ ′ r ~ m ′ ω λ ′′ × φ k µ ψ λ ′ ψ λ ′′ (cid:2) ( D r + D φ ) k k δ µ, k + ( S r + S φ ) k ⊥ δ µ, ⊥ (cid:3) , (18)and the term on B k µ ,λ ′′′ ,λ ′′′′ becomes the same as A k µ ,λ ′ ,λ ′′ by setting [ λ ′ , λ ′′ −→ λ ′′′ , λ ′′′′ ]. The last one should be C k µ ,λ ′ ,λ ′′′ = − π X l i s ~ ρω k µ r ~ m ′ ω λ ′ r ~ m ′ ω λ ′′′ × φ k µ ψ λ ′ ψ λ ′′′ (cid:2) ( D r − D φ ) k k δ µ, k + ( S r − S φ ) k ⊥ δ µ, ⊥ (cid:3) . (19)The squared quantity on Eq. (18) is given by A k µ ,λ ′ ,λ ′′ = CI V L λ ′ L λ ′′ ω k µ ω λ ′ ω λ ′′ (cid:2) ( D r + D φ ) k k δ µ, k + ( S r + S φ ) k ⊥ δ µ, ⊥ (cid:3) , (20)where the coefficient C is defined as C = ~ Ω ρm ′ . (21)The expression of B k µ ,λ ′′′ ,λ ′′′′ takes the same form as A k µ ,λ ′ ,λ ′′ , since they both correspond to the interactionbetween SL modes with the same polarization. While C k µ ,λ ′ ,λ ′′′ corresponding to interaction between differ-ent polarizations has an additional factor (4 /π ) and (cid:2) ( D r − D φ ) k k δ µ, k + ( S r − S φ ) k ⊥ δ µ, ⊥ (cid:3) . V. HOPPING PROCESSA. Relaxation time of SL modes
This subsection gives the formula for the relaxation time ofSL mode due to the scattering process EX+SL → SL (hoppingprocess) together with its reverse process shown in Fig. 5 byapplying the Fermi golden rule. To obtain the total transitionrate of the SL mode in λ ′ , we have to incorporate all of fourprocesses for each polarization as given below. These providethe decay of the Bose-Einstein distribution function n λ ′ forthe occupied state λ ′ , dn λ ′ dt = 2 π ~ X k µ ,λ ′′ | A k µ ,λ ′ ,λ ′′ | [ n λ ′′ (1 + n k µ )(1 + n λ ′ ) − n k µ n λ ′ (1 + n λ ′′ )] δ ( ω λ ′′ − ω λ ′ − ω k µ )+ | A k µ ,λ ′′ ,λ ′ | [ n k µ n λ ′′ (1 + n λ ′ ) − n λ ′ (1 + n k µ ) × (1 + n λ ′′ )] δ ( ω λ ′ − ω λ ′′ − ω k µ )+ (cid:2) A k µ ,λ ′ ,λ ′′ −→ C k µ ,λ ′ ,λ ′′′ , λ ′′ −→ λ ′′′ (cid:3) (22)We consider, at first, the decay due to the hopping processbetween the same polarization, i.e. , the contribution from thefirst two terms of Eq. (22). By separating the distribution func-tion into two parts. n = n (0) + n (1) , where n (0) is the Bose-Einstein distribution function in equilibrium state and n (1) isits deviation due to the scattering processes, and by employingthe relaxation time approximation, dn (1) λ ′ /dt = − n λ ′ (1) /τ λ ′ ,we have the inverse of relaxation time from Eq. (22) for thesame polarization process, τ same λ ′ ∼ = 2 π ~ CI V L X k µ ,λ ′′ ω k µ ω λ ′ ω λ ′′ × (cid:2) ( D r + D φ ) k k δ µ, k + ( S r + S φ ) k ⊥ δ µ, ⊥ (cid:3) × [ δ ( ω λ ′′ − ω λ ′ − ω k µ )( n (0) k µ − n (0) λ ′′ )+ δ ( ω λ ′ − ω λ ′′ − ω k µ )(1 + n (0) k µ + n (0) λ ′′ )] , (23)where the explicit form of the summation I arising from theoverlapping of wave functions ψ λ ′ and ψ λ ′′ is given by I = X ℓ e − i k µ · R l cos [ k λ ′ · ( R l − R λ ′ )] e −| R ℓ − R λ ′ | /L λ ′ × cos [ k λ ′′ · ( R l − R λ ′′ )] e −| R ℓ − R λ ′′ | /L λ ′′ . (24)The above sum I can be reduced to, by taking the origin ofthe sum as R λ ′ = 0 and the nearest neighbor position fromthe origin as R λ ′′ = ∆ R λ ′′ , I = X ℓ f ( R ℓ ) f ( R ℓ − ∆ R λ ′′ ) e − i k µ · R ℓ , (25)where the even function f ( X ℓ ) is defined as f ( X ℓ ) = cos( k λ ′ · X ℓ ) e −| X l | /L λ ′ . (26)Since the localization lengths of SL modes are the same, e.g. , L λ ′ ∼ = L λ ′′ , hereafter we denote this as L . As f ( X ℓ ) is con-cerned with SL modes, the relevant sum should be made inthe region | X ℓ |≤ L , so we can approximate the summationby I ∼ = 1Ω Z | X ℓ | 12 cos (2 k λ ′ · R ℓ − k λ ′ · ∆ R λ ′′ ) + 12 cos ( k λ ′ · ∆ R λ ′′ ) (cid:21) ∼ = | ∆ R λ ′′ | πL e −| ∆ R λ ′′ | /L . (27)where we have used the approximation cos( k λ ′ · ∆ R λ ′′ ) ≈ cos( k λ ′ nL ) ≈ from the Ioffe-Regel condition L ≈ π/k λ ′ for SL modes and e − i k µ · R ℓ ≈ due to | k µ |≪ π/L forthe wave number of EX acoustic modes. The term contain-ing cos (2 k λ ′ · R ℓ − k λ ′ · ∆ R λ ′′ ) becomes negligible since ityields rapidly oscillating function in the integrand.This gives the squared hopping integral of the form I ≃ (cid:18) π ∆ R λ ′′ L (cid:19) e − R λ ′′ /L , (28)where ∆ R λ ′′ is the hopping distance.In the temperature regime T ≃ a few 10 K, i.e. , k B T > ~ ω λ ′ , ~ ω λ ′′ > ~ ω k µ , the inverse of the relaxation time takesthe following form under the above conditions and by em-ploying the linear dispersion relation for EX phonon mode ω k µ = v µ k µ , τ same λ ′ ∼ = 2 πk B T ( D r + D φ ) CI ~ V L v k × X k k ,λ ′′ " δ ( ω λ ′′ − ω λ ′ − ω k k ) ω λ ′′ + δ ( λ ′ ⇄ λ ′′ ) ω λ ′′ + [2 × ( D −→ S, k−→⊥ ) in the above] . (29) Here the coefficient C is defined in Eq. (21). We have omittedthe temperature independent term providing only small contri-butions. B. Thermal conductivity due to the hopping of SL modes In the previous subsection, we have formulated the relax-ation rate of SL modes due to the anharmonic interaction be-tween SL modes and EX modes. This is a quantum processrealizing the decay of SL ′ mode to SL ′′ mode assisted by EXmode: SL ′ +EX → SL ′′ . Without anharmonic interaction, SLmodes cannot diffuse/contribute to thermal transport. Thismeans that the plateau region should continue over at highertemperatures after exhibiting the onset of the plateau, i.e. , thecontribution from EX modes to lattice thermal conductivityis saturated at higher temperatures. This is because the onsetof the plateau arises from the weak localization of acousticmodes as explained in Sec.II. Thus, the T -linear rise of κ L ( T ) cannot recover without anharmonic interaction between SLmodes and EX modes. In addition, we emphasize that disorder, induced by off-centeredness as shown in Supplemental Material, is essentialto generate the hopping of SL modes. This occurs only inthe case that SL ′ mode belonging to the eigenfrequency ω SL ′ can hop to a site of SL ′′ mode with a different eigenfrequency ω SL ′′ via absorption or emission of EX mode with finite fre-quency ± ( ω SL ′ − ω SL ′′ ) . This finite frequency is created bylevel repulsion between eigenfrequencies due to disorder, i.e. ,localized modes never belong to the same eigenfrequency ac-cording to the level repulsion.Let us provide the formula of κ L ( T ) due to the diffusionprocess where SL modes serve as primary heat carriers. Inthis process, the characteristic length-scale should be the hop-ping distance ∆ R λ ′′ from the site of SL ′ mode to that of SL ′′ mode, and the characteristic time-scale is the relaxation time τ λ ′ of the SL ′ mode. This leads to the following formula ofthe lattice thermal conductivity due to the hopping process,which was first proposed for fracton excitations by Alexander et. al. , κ hop ( T ) = 13 V X λ ′ C λ ′ ( T ) ∆ R λ ′′ τ λ ′ , (30)where ∆ R λ ′′ /τ λ ′ is the thermal diffusivity of SL mode λ ′ , C λ ′ ( T ) is the specific heat associated with the SL mode λ ′ .In the high temperature regime above the plateau region T & afew 10 K, the specific heat follows the Dulong-Petit relationof the form C λ ′ ( T ) = k B per one polarization of SL mode λ ′ . Note that /τ λ ′ = 1 /τ same λ ′ + 1 /τ dif λ ′ , we first calculate thehopping process between the same polarization by, κ samehop ( T ) = k B V X λ ′ ∆ R λ ′′ τ same λ ′ . (31)The substitution of Eq. (29) into Eq. (31) together withEq. (28) yields κ samehop ( T ) ∼ = k B V π k B T ( D r + D φ ) C ~ v k L Ω X k k ,λ ′ ,λ ′′ ∆ R λ ′′ ω λ ′′ × e − R λ ′′ /L (cid:2) δ ( ω λ ′′ − ω λ ′ − ω k , k ) + δ ( λ ′ ⇄ λ ′′ ) (cid:3) + [2 × ( D −→ S, k−→⊥ ) in the above] (32) Transforming the sum P k µ for EX phonon modes to theintegral V / (2 π ) R d k µ = V / (2 π v µ ) R ω k µ dω k µ , we have κ samehop ( T ) = πk B T C ~ V Ω L " ( D r + D φ ) v k + 2 ( S r + S φ ) v ⊥ × X λ ′′ ,λ ′ ∆ R λ ′′ e − R λ ′′ /L ( ω λ ′′ − ω λ ′ ) ω λ ′′ (33)The sum on λ ′ and λ ′′ above should include the density ofstates of SL modes D SL ( ω λ ′ ) and D SL ( ω λ ′′ (∆ R λ ′′ )) for thesame polarization process. The volume Ω should contain twoindependent SL modes corresponding to two independent in-plane mode, say, stretching or libration, in the band width of ∆ ω sl , which leads to D SL ( ω λ ′ )Ω∆ ω sl = 2 . (34)and D SL ( ω λ ′′ (∆ R λ ′′ ))Ω∆ ω sl = 1 . (35)where the volume Ω contains at least one possible SL mode λ ′′ with the same/different polarization as/from mode λ ′ . Sincethe term ∆ R λ ′′ e − R λ ′′ /L in Eq. (35) achieves its maximumat ∆ R λ ′′ = 2 L and it decays fast with the further increasingof ∆ R λ ′′ , the sum of λ ′′ could be estimated within the sphereregion ∆ R λ ′′ ≤ ∆ R . X λ ′′ ,λ ′ ∆ R λ ′′ e − R λ ′′ /L ( ω λ ′′ − ω λ ′ ) ω λ ′′ ∼ = π ∆ R V Ω ∆ R e − R/L × (10 − ) (36)Here the sum on SL modes are done by P λ ′′ =4 π ∆ R / R ω sl +∆ ω sl ω sl D ( ω λ ′′ (∆ R λ ′′ )) dω λ ′′ and P λ ′ = V R ω sl +∆ ω sl ω sl D ( ω λ ′ ) dω λ ′ , where the factor π ∆ R / from Eq. (35) means the total number of hopping sitesfrom λ ′ to λ ′′ for the same polarization process, and V / Ω from Eq. (34) is the total number of λ ′ contribut-ing the thermal conductivity κ hop . The numerical fac-tor − arises from the magnitude estimation of integral R ω sl +∆ ω sl ω sl dω λ ′ R ω sl +∆ ω sl ω sl dω λ ′′ ( ω λ ′′ − ω λ ′ ) ∆ ω ω λ ′′ .The formula of the thermal conductivity due to the hoppingmechanism is given by κ samehop ( T ) = π k B T ∆ R ρm ′ Ω L e − R/L × (10 − ) " ( D r + D φ ) v k + 2 ( S r + S φ ) v ⊥ (37)The same procedure for the hopping process due to anhar-monic interaction between different polarizations leads to κ difhop ( T ) = 4 k B T ∆ R ρm ′ Ω L e − R/L × (10 − ) " ( D r − D φ ) v k + 2 ( S r − S φ ) v ⊥ (38)The total thermal conductivity due to the hopping mechanismis given by the sum of these components as κ hop ( T ) = κ samehop ( T ) + κ difhop ( T ) , (39) C. Evaluation of anharmonic coupling D and S Here we estimate the anharmonic coupling constants D r ( φ ) and S r ( φ ) by illustrating type-I BGS. The coupling constants D r ( S r ) and D φ ( S φ ) are associated with the stretching and li-bration motion of guest-cage vibrations identified by the forceconstant ξ r and ξ φ in Eq. (6) by the relation ξ r ( φ ) = m ′ ω r ( φ ) ,where m ′ is the reduced mass defined by /m ′ = 1 /m +1 /M .In our coarse-grained Hamiltonian introduced in Sec. III, theguest ion Ba(2) in tetrakaidecahedron cage has the mass m and the molecular unit composed of 1 tetrakaidecahedron and1/3 dodecahedron does the total mass M excluding the off-center guest ion.We first evaluate the coupling constants D r ( φ ) from the Ra-man spectroscopy data of pressure dependence . The D r can be related to the pressure P by D r = ∂ξ r ∂e αα = 3 B ∂ξ r ∂ω r ∂ω r ∂P = 3 B (2 m ′ ω r ) ∂ω r ∂P . (40)Here B = ∆ P (∆ V/V ) is the linear thermal expansion coefficient,where the dilation is given by ∆ V /V = P α e αα for cubicstructure. The coupling constant D φ can be defined in a simi-lar manner to Eq. (40) as D φ = ∂ξ φ ∂e αα = 3 B (2 m ′ ω φ ) ∂ω φ ∂P . (41)In the pressure range from 0.8 GPa to 5.8 GPa, E g mode spansfrom 20 cm − to 27 cm − . While, for T mode, it rangesfrom 17 cm − to 27 cm − . The observed spectra of thesetwo modes are overlapped/mixed. Taking account of these as-pects, we have ∂ω r /∂P = 2 π × . × [sec − GPa − ] and ∂ω φ /∂P = 2 π × . × [sec − GPa − ] . We then obtainthe coupling constants D r = m ′ π × . × [kg · sec − ] and D φ = m ′ π × . × [kg · sec − ] using the observedbulk modulus B = 41 . . Within our knowledge, theexperiment data for estimating the coupling coefficients S r ( φ ) are not available, so we assume as S r ≈ D r and S φ ≈ D φ atthe present stage. The above coupling constants yield κ hop = 3 . × − T (Wm − K − ) , (42)where we have employed the values of parameters in Eq. (39)as the localization length L = 2 a , the hopping distance ∆ R = 3 . L , the volume of molecular unit Ω = ( a ) / ,the lattice spacing a = 11 . ˚A, the mass density ρ =6 . × kg/m , in addition to the velocities of acousticphonons v k = 3369 m/s and v ⊥ = 1936 m/s . The valueof κ hop in Eq. (42) is smaller than the observed one of κ hop =9 . × − T (Wm − K − ) for type-I BGS. This mainly arisesfrom, as will be demonstrated below by means of FPC, theunderestimated shear coupling constants S r ( φ ) obtained by as-suming the relations S r ( φ ) ≈ D r ( φ ) .Due to the lack of experiment data for the shear couplingcoefficients S r ( φ ) , we have performed FPC for type-I BGSto obtain the coupling constants from the shift of eigenfre-quencies at Γ -point of low-lying optical mode by imposingstrain to the cage structure. The normal strain is isotropicand defined as e αα = ( a − a ) /a where a and a are thelattice constant for the unstrained and strained unit cell , re-spectively. The shear strain is also isotropic and defined as e αβ = (1 − p − (2 cos θ − 1) cos θ ) / (2 cos θ − where θ is the acute angle between edges after deformation.We have performed the FPC by the VASP code with thePerdew-Burke-Ernzerhof functional and the PAW method ,plane wave cut-off energy 250 eV and the force convergenceless than − eV / ˚ A . The phonon frequencies are calculatedby PHONOPY code with the × × Monkhorst-Pack k grids and for a unit cell containing 54 atoms. The couplingconstants obtained from normal strain are D r = m ′ π × . × [kg · sec − ] , D φ = m ′ π × . × [kg · sec − ] , andfrom sheared unit cell are S r = m ′ π × . × [kg · sec − ] , S φ = m ′ π × . × [kg · sec − ] , respectively. The D r ( φ ) are smaller than those estimated from the Raman spectroscopydata of pressure dependence, though S r ( φ ) are larger than thevalues obtained from the assumption S r ( φ ) ≈ D r ( φ ) . Theabove coupling constants yield the thermal conductivity dueto the hopping of SL modes of κ hop = 4 . × − T (Wm − K − ) . (43)We remark here that our FPC provides the results for theon-center positioned Ba(2) because the optimization for off-center structure is quite time-consuming and may require totake into account the dipole-dipole interaction due to off-centeredness and temperature effect. The on-center struc-ture gives rise to the underestimated coupling constants S since on-center guest ions should more weakly response toshear distortion than the case of off-center. Then, the ac-tual S r ( φ ) should be larger than the above estimation. Un-der these situations, the calculated value in Eq. (43) providesreasonable agreement, to claim the relevance of the hoppingprocess of SL modes, with the observed κ hop = γT with γ = 9 . × − Wm − K − for type-I BGS , and γ =9 . × − Wm − K − for type-I EGG . For type-I SGG, sev-eral different values around γ ∼ . × − Wm − K − havebeen reported , indicating that the experimental data ofSGG depend on sample qualities according to synthesis meth-ods. In that respect, it has been reported that a flux-grownsample shows a glasslike plateau, while a zone-melted samplehas a crystalline peak. VI. SUMMARY AND CONCLUSIONS Off-center type-I clathrates show almost identical lat-tice thermal conductivities κ L to those of structural glasses . In addition, off-center type-I clathrates show theexcess density of states at THz frequencies manifesting theboson peak identical to those of network-forming glasses .These indicate that the symmetry broken guest ions in cagestake charge of the emergence of glasslike κ L ( T ) . In structural glasses, many key aspects of a detailed quantitative descrip-tion are still missing. This is due to the difficulty to identifyrelevant entities or elements at atomic scale caused by theircomplex microscopic structures.In Sec. II, we have pointed out that the PR shown inFig. 2 provides the evidence that EX acoustic phonons car-rying heat convert to WL modes modes at ∼ ≈ B from the Wien’s displacement law, sothat this conversion should be associated with the onset of theplateau thermal conductivities observed at several K in off-center type-I clathrates .With increasing temperature further, thermal conductivitiesabove a few 10 K show a linear rise on temperature. This typeof anomalous thermal conductivities with the plateau and thesubsequent T -linear rise have been clearly observed for off-center type-I clathrates . This is the prominent hallmarkof glasslike thermal conductivity since crystals with transla-tional invariance never show these features. Rather, latticethermal conductivities of crystallines decrease with increas-ing temperature proportional to κ ( T ) ∝ /T known as theUmklapp process .The theoretical elucidation on the linear rise on temperature“above” the plateau region has been the main subject of thepresent paper. Our calculated results given in Sec. V, based onhopping process, show fairly good agreement with observedthermal conductivities above the plateau. We particularly em-phasize that both the magnitude and the temperature depen-dence of κ ( T ) are in accord with the experimental data .At much higher temperatures, the T -linear rise in κ ( T ) doesnot continue, but κ ( T ) saturates above T ≃ 100 K . Inthis temperature regime, the treatment based on quantum me-chanical process does not hold for since the life-time of ex-cited modes becomes much smaller than the inverse of theirangular frequencies, where the guest ions become free fromthe constraint of atoms constituting cages. This subject willbe discussed in detail elsewhere .In conclusion, the phenomenon of T -linear rise of κ L ( T ) above a few 10K in off-center type-I clathrates has beenquantitatively explained by analytic theory, on the groundsthat off-center clathrates possess definite microscopic struc-ture. Our successful clarification in quantitative manner isowing to the fact that the systems are more tractable thannetwork-forming glasses with the difficulty to identify rele-vant constituents at atomistic level caused by their complexmicroscopic structures. Acknowledgments. This work is supported by the Na-tional Natural Science Foundation of China Grant No.11334007 and No. 51506153. J. Z. is supported by theprogram for Professor of Special Appointment (EasternScholar) 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, GrandNo.26400381.0 ∗ [email protected] † [email protected] ‡ [email protected] G. A. Slack, in CRC Handbook of Thermoelectrics , edited by D.M. Rowe (CRC Press, Boca Raton, FL,1995), pp.407-440. See a review, for example, and references therein, T. Takabatake,K. Suekuni, T. Nakayama, and E. Kaneshita, Rev. Mod. Phys. ,669 (2014). M. Beekman, D. T. Morelli, and G. S. Nolas, Nat. Mater. , 1182(2015). G. S. Nolas, J. L. Cohn, G. A. Slack, and S. B. Schujman, Appl.Phys. Lett. , 178 (1998). J. S. Cohn, G. S. Nolas, V. Fessatidis, T. H. Metcalf, and G. A.Slack, Phys. Rev. Lett. , 779 (1999). S. Christensen, M. S. Schmokel, K. A. Borup, G. K. H. Madsen,G. J. Mclntyre, S. C. Capelli, M. Christensen, and B. B. Iversen,J. Appl. Phys. , 185102 (2016). S. Paschen, W. Carrillo-Cabrera, A. Bentien, V. H. Tran, M.Baenitz, Y. Grin, and F. Steglich, Phys. Rev. B , 214404 (2001). B. C. Sales, B. C. Chakoumakos, R. Jin, J. R. Thompson, and D.Mandrus, Phys. Rev. B , 245113 (2001). M. A. Avila, K. Suekuni, K. Umeo, H. Fukuoka, S. Yamanaka,and T. Takabatake, Phys. Rev. B , 125109 (2006). M. A. Avila, K. Suekuni, K. Umeo, H. Fukuoka, S. Yamanaka,and T. Takabatake, Appl. Phys. Lett. , 041901 (2008). K. Suekuni, M. A. Avila, K. Umeo, H. Fukuoka, S. Yamanaka, T.Nakagawa, and T. Takabatake, Phys. Rev. B , 235119 (2008). K. Suekuni, M. A. Avila, K. Umeo, and T. Takabatake, Phys. Rev.B , 195210 (2007). See, for example, a review, T. Nakayama, Rep. Prog. Phys. ,1195 (2002). Y. Liu, Q. Xi, J. Zhou, T. Nakayama, and B.Li, Phys. Rev. B ,214305 (2016). S. Alexander, O. Entin-Wohlman, and R. Orbach, Phys. Rev. B , 2726 (1986). M. L. Williams and H. J. Maris, Phys. Rev. B , 4508(1985);K. Yakubo, T. Nakayama, and H. J. Maris, J. Phys. Soc. Jpn. ,3249 (1991). See a review, for example, T. Nakayama and K. Yakubo, Phys.Rep. , 239 (2001). J. B. Suck, M.Schreiber and P.H¨aussler, Quasicrystals: Anintroduction to structure, physical properties and applications (Springer, Berlin, 2002), pp.403. T. Nakayama, Phys. Rev. Lett. , 1244 (1998), T. Nakayama andN. Sato, J. Phys. Condens. Matter , L41 (1998). T. Mori, K. Iwamoto, S. Kushibiki, H. Honda, H. Matsumoto, N.Toyota, M. A. Avila, K. Suekuni, and T. Takabatake, Phys. Rev.Lett. , 015501 (2011). Y. Takasu, T. Hasegawa, N. Ogita, M. Udagawa, M. A. Avila, K.Suekuni, I. Ishii, T. Suzuki, and T. Takabatake, Phys. Rev. B ,174303 (2006). T. Kume, T. Sukemura, S. Nakano, S. Sasaki. K. Suekuni, and T.Takabatake, Photon Factory Activity Report 2014, , B (2015).;T. Sukemura, T. Kume, T. Matsuoka, S. Sasaki, T. Onimaru, andT. Takabatake, J. Phys.: Conf. Ser., , 182022 (2014). Isao Ishii, Yasuhiko Suetomi, Takahiro K. Fujita, KoichiroSuekuni, Tomoo Tanaka, Toshiro Takabatake, and TakashiSuzuki, Phys. Rev. B , 085101 (2012). J. Chen, J H. Walther, and P. Koumoutsakos, Nano Lett. , 819(2014). G. Kresse and J. Furthm¨uller, Phys. Rev. B , 11169 (1996). J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. ,3865 (1996). A. Togo and I. Tanaka, Scr. Mater. , 1 (2015). E. M. Lifshitz and L. P. Pitaevskii, in Physical Kinetics (Elsevier,Amsterdam, 1979), Chapter 68.29