Theory of Ultrasonic Dispersion in Local Phonon Systems Coupled with Conduction Electrons

The physical origin of frequency dependence in elastic constants, which are often found in an ultrasound propagation in filled skutterudites and clathrate compounds, is investigated theoretically. This dependence arises from a coupling between the acoustic phonon and some optical phonons, which strongly interact with electrons. Using a self-consistent ladder approximation together with a pseudofermion mapping of the phonon to the single site Holstein Anderson model, a soft mode of the optical phonon at zero frequency is shown to emerge. The temperature dependence of the spectral weight of this soft mode shows an activation-type behavior, which is characterized by the optical phonon frequency. These features can generate the frequency dependence and the shoulder in the elastic constants observed in some filled skutterudites and clathrate compounds.


Introduction
Recently, the frequency dependence of ultrasonic velocity in some filled skutterudites and clathrate compounds has attracted a great deal of attention. [1][2][3][4][5] In these systems, it is believed that a rare-earth ion in a rigid cage moves rather independently in the anharmonic potential of off-center minima. In PrOs 4 Sb 12 , for example, it is expected that these ionic degrees of freedom are important to heavy fermion superconductivity and its large effective quasiparticle mass. 6 The data of elastic constants (EC) in these compounds exhibit frequency dependence at temperatures ranging from 10 to 40 K. This temperature range is on the order of the Einsteinlike modes of the rare-earth oscillation in these compounds. [7][8][9][10][11] For PrOs 4 Sb 12 , large softening of acoustic phonons is also reported. 12 Recently, Iwasa et al., observed the development of a quasi-elastic peak at low temperatures in Pr-based filled skutterudites. 13 Usually, the frequency dependence of EC manifests as a shoulder in its temperature dependence. This is well fitted by a Debye-type formula: where C(z) is the frequency-dependent EC and C ∞ ≡ C(∞) and C 0 ≡ C(0) are constants in frequencies.
In the present paper, we use z as the frequency of dynamic quantities. The lifetime τ in eq. (1) is given as τ = τ 0 exp(E ∆ /T ), where E ∆ is a characteristic energy scale and T is the temperature. The exponential T -dependence gives the shoulder at zτ ∼ 1, at approximately T ∼ 40 K for E ∆ ∼ 200 K, even though the frequency of the ultrasound is of the order of mK. This is called "ultrasonic dispersion" (UD). This form of the relaxation time appears to be related to a kind of thermal activation-type process, so that in ultrasonic experiments it has been interpreted that there is an off-center potential and the origin of the relaxation is due to the thermal hopping between the off-center sites. In this respect, we discussed a possible scenario for the realization of heavy fermions in SmOs 4 Sb 12 14 in the strong coupling limit of electron-phonon couplings, 15,16 which is a natural extension of a two-level system. [17][18][19][20][21] Since the UD is observed only in a specific mode, it is claimed that the off-center degrees of freedom have a degenerate and anisotropic ground state. On the other hand, the results of a neutron scattering experiment in PrOs 4 Sb 12 suggests that there are no anisotropic charge distributions of the Pr-nuclei at low temperatures in particular. 22 In order to discuss the frequency dependence in EC, it is convenient to treat the two frequency regions separately. One typical region is zτ ≪ 1, where the time scale of the ultrasound is much longer than that of the relaxation mode. At this stage, the origin of the relaxation mode is unknown. The relaxation mode is scattered in a much shorter time than the oscillation period of the ultrasound. Thus, the strain caused by the ultrasound can be regarded as static. The ultrasound attenuation due to the anharmonicity of the lattice in this region was discussed by Akhieser, 23 and later by Woodruff and Ehrenreich, 24 using the semi-classical approach (Boltzmann equation). For the case of zτ ≫ 1, the Landau and Rumer theory, which is "golden rule" treatment, can be applied. 25 Although the former description appears to be good even in zτ > 1 qualitatively, there is no theoretical reasoning to use the former in this region.
Using the Woodruff and Ehrenreich theory, the frequency dependence in the EC was discussed in KTaO 3 by Barrett 26 from both theoretical and experimental points of view. In KTaO 3 , the activation factor exp(ω 0 /T ) also appears in the ultrasound velocity, where ω 0 roughly corresponds to the frequency of the optical mode. This factor arises from the specific heat of the free opti-cal mode through the time-dependent Boltzmann equation. In comparison, the theory for the frequency dependence in filled skutterudite and clathrate compounds is more difficult, since we need to bridge the above two regions continuously and to take into account the electronphonon couplings.
In the present paper, we discuss a simple interpretation of this problem. In §2, we review the Green function formulation in electron-phonon systems. We then discuss two types of couplings between the optical and acoustic phonons, and give the phenomenological interpretation of the Debye formula used in the ultrasonic experiments. In §3, we investigate the Holstein-Anderson model at finite temperatures using the self-consistent ladder approximation. We then show the results of the renormalized sound velocity, giving the dispersion in EC, together with the nature of soft optical phonon modes. In §4, we discuss the applicability of the present theory to the real materials and future problems. Finally, in §5, we summarize the results.

Coupling between acoustic and optical phonons
In this section, we discuss how eq. (1) is derived in interacting electron-phonon systems. Note that the displacements due to the acoustic modes near q = 0 include contributions from optical modes at q = 0. This essential point is described in §2.1. In §2.2, we discuss the case in which the anharmonicity of the lattice is important.

Harmonic coupling
The classical Hamiltonian for the phonon system in the harmonic approximation is written using a force constant matrix K as where P µ i (X µ i ) is the momentum (displacement) variable of the ion in the µ direction located at the n-th site in the i-th unit cell. M n is the mass of the ion at the n-th site in a unit cell. Introducing variables p µ the Hamiltonian (2) is reduced to where we have introduced the Fourier components of each variable in the second line. After diagonalizing the k-matrix, we obtain all of the eigenvalues in principle as in textbooks. 28 In the usual step of moving to the corresponding quantum Hamiltonian, we use the diagonalized basis. Although we can diagonalize the k-matrix and obtain the energy eigenvalues of every phonon branch, we here introduce even and odd parts, k e and k o , with respect to q, and concentrate on the low q limit of the phonon system. In order to determine the phonon velocity, we need only retain up to the second-order contribution in q. Thus, we write k e = k (0) + k (2) and k o = k (1) . Dividing k into k e and k o is meaningful only in the region near q = 0. Near q = 0 and z = 0 (z being the frequency of phonons), k o = O(|q|) acts as the hybridization between some optical and acoustic modes at q ≃ 0. Since we are interested in the sound velocity of the system under which optical phonons strongly interact with conduction electrons (this comes from the fact that UD has an activation type relaxation), dividing k into the even and odd parts is meaningful. For large |q| > ∼ a −1 (a being the lattice constant), this characterization loses its meaning, and k o should be regarded simply as the hybridization.
In order to discuss the quantum mechanical Hamiltonian, we first diagonalize k e . We define the creation and annihilation operators of the mode λ with the energy eigenvalueω qλ (related to the eigenvalue of k e , not k) as follows: Here,x qλ is the displacement operator that diagonalizes k e . Using this set of displacement variables, Hamiltonian (3) can be written as Next, we introduce the retarded Green functionD qλλ ′ (t) in the usual manner: The equation of motion forD qλλ ′ (t) has the following form: It is useful to employ D qλλ ′ (t), which is defined as The Fourier components of D qλλ ′ (t) obey the following coupled equations: Here, we have written the mixing term 1 The energy eigenvalues are determined by the poles of D qλλ (z). It is evident that these poles coincide with those in usual basis, which diagonalize k e + k o .
If we take the electron-phonon coupling into account, we must introduce the polarization functions Π qλλ ′′ (z) into the left-hand side of eq. (10). Defining (D −1 0q (z)) λλ ′′ as the matrix in the large brackets of the left-hand side in eq. (10) divided by 2ω qλ 2ω qλ ′′ , the equation of motion for D qλ ′′ λ ′ (z) in the interacting electron-phonon system can be written as In order to clarify the importance of κ λλ ′′ (q), it is instructive to consider the simplest case in which there exist one acoustic (its energyω q1 ) and one optical mode (ω q2 ) in a system we consider. Furthermore, the electronphonon couplings are nonzero only for the optical mode and κ λλ (q) = 0 by definition. Here, we use the words "acoustic" and "optical" with respect to those in the basisx qλ . In the diagonalized basis, there is no mixing between the acoustic and optical phonons in the harmonic Hamiltonian (2). Then, eq. (11) is reduced to We find the form of the Green function for the acoustic mode as . (13) Noted that in eq. (13), κ 12 (q) ≃ κ 0 (q)|q|,ω q2 ≃Ω = const. andω q1 ≃ũ 0 (q)|q| for small q, whereq ≡ q/|q|. Although the expressionω q1 ≃ũ 0 (q)|q| is not satisfied in general, it is a reasonable assumption for filled skutterudite and clathrate compounds. Diagrammatically, eq. (13) is represented by Fig. 1(a). Thus, the frequencydependent phonon velocity u(z,q) is given formally by In the limit of ultrasound, we can ignore z ≪Ω and set Πq 22 = Π ′q − iΓqz in eq. (14), leading to where τq ≡ Γq/(Ω + 2Π ′q ). Equation (15) corresponds to the phenomenological expression eq. (1) used in the analysis of the ultrasound experiments. In phenomenological theory, the temperature dependence of τq is assumed to be proportional to exp(E ∆ /T ). The activation scale E ∆ is expected to be of the order ofΩ. The fact that the ultrasound experimental results are well explained by eq. (1) indicates the existence of a soft mode of the phonon (in this simple example, the soft mode, which appears around zero frequency, is related to the manybody effects between the optical phonon and conduction electrons). In the above simple model, we have ignored Π qnm with (n, m) = (2, 2). This is because we consider that most of the important properties come from the self-energy of the optical phonon. If these come from the degrees of freedom in the lattice vibrations, it is clear that E ∆ cannot appear without the finite energy excitation of the optical phonon.

Anharmonic coupling
The discussion in the previous subsection, however, is not directly applicable to the rare-earth mode in filled skutterudites because the optical phonons that hybridize with the acoustic phonons near q = 0 are those with even parity in that system. The rare-earth mode couples with the acoustic phonon only through the anharmonic coupling so long as we ignore the electron-acoustic phonon interactions, as shown in Fig. 1(b). Since the ultrasound generates strain ε Γ (q) with the wave number q and the symmetry Γ = Γ + 1 , Γ + 3 and Γ + 5 in the cubic O h point group, the rare-earth mode can couple with these strain-fields by constructing the direct product: is expressed by the acoustic phonon fields.) Thus, we have coupling due to the anharmonicity as where φ − is the phonon field with the odd parity and g − Γ is the coupling constant. Here we use the symbolic should be read as the Fourier component of (φ ix φ iy ) with the wave number −q. Corresponding to V − , the coupling between the optical and acoustic phonons in the previous subsection can be rewritten as The EC's are renormalized by the susceptibility of (φ − φ − ) Γ (denoted as χ Γ φφ (z, q)). A diagrammatic expression is shown in Fig. 1(b). Information of χ Γ φφ is required in order to obtain an explicit form of the sound velocity. This is more difficult than in the case of the previous subsection. In phenomenological treatment, we can assume χ Γ φφ ∝ 1/(1 − izτ q ), which is the same form as that in eq. (15).
We summarize this section as follows: • High-temperature region: the sound velocity is given by "original" sound velocityũ 0 (q) subtracted by the contribution from the "non-interacting" optical phonon |κ 0 (q)| 2 /(Ω 2 + 2ΩΠ ′q ) (see eq. (15)) because of the smallness of zτq. This corresponds to the sound velocity obtained by diagonalizing k e + k o with renormalizedΩ. Thus, the sound velocity is observed as a smaller value than that at lower temperatures as explained below. • Low-temperature region: the relaxation time τ q becomes very large around z = 0, and eventually zτq exceeds 1. Note that the relaxation time of the optical phonon at the original pole is much smaller than τq. As a result, the second term in eq. (15) does not affect the sound velocity. Thus, we observe only the "original" sound velocityũ 0 (q). These aspects are similar to those discussed by Yamada for the critical fluctuations of lattice systems. 29 In these studies, the system is assumed to be located near the structural phase transition and the imaginary part of the polarization was phenomenologically introduced. In this paper, however, we will calculate the explicit temperature dependence of τq, namely ∝ exp(E ∆ /T ), based on a relevant microscopic model. In the next section, we investigate a simple model with respect to the discussions in §2.1 for simplicity. However, in principle, it is straightforward to discussion the same topic based on §2.2.

Model and Calculation
In this section, we concentrate on a model in which there is only one electron band and one local Einstein phonon. To carry out the complete calculation, it is desirable to evaluate the properties of the model with two kinds of phonons (whose Green functions are formally given by eq. (12)) and conduction electrons. However, in order to obtain a qualitative understanding, we restrict ourselves in one local phonon coupled with conduction electrons and approximate the frequency dependence of the sound velocity in §3.3.
One of the simplest models is the Holstein model, 27 which is given as where t ij represents the hopping of conduction electrons, c † iσ is the creation operator of electrons at site i and spin σ, b † i is the phonon creation operator at site i with Einstein energy Ω E , and g is the electron-phonon coupling constant. This model has been discussed for the past five decades. [30][31][32] Recently, the phase diagram of this model at T = 0 was discussed based on the dynamical mean field theory (DMFT) 33 using the numerical renormalization group (NRG) as an impurity solver. [34][35][36][37] The phase diagram of a somewhat different model 38 was discussed by DMFT, but with exact diagonalization. For finite temperatures, data is lacking for frequencies smaller than the temperature in the NRG method. Recent developments in calculating spectral functions in NRG enable us to obtain the spectral functions roughly in the range of z > ∼ T × 1/10. 39,40 However, it is insufficient to discuss the soft mode at finite temperatures. Thus, in the following subsections, we introduce a simple self-consistent treatment of the phonon system, which can capture the essential aspect of the emergence of the soft mode at low temperatures. For the connection with the single site approximation, which is widely used for experimental analysis of the quadrupolar degrees of freedom in f-electron systems, 41 we restrict ourselves to the single-site problem of a local phonon interacting with conduction electrons. This model is called the Holstein-Anderson model, which will be explained below.
3.1 Self-consistent theory in pseudo-fermion representation of phonon First, we explain a pseudo-fermion mapping of phonon operators. The displacement operatorX = b † + b has the matrix elements in the phonon Hilbert space as follows: where n and m are eigenvalues of the number of the phonons (b † b). We introduce pseudo-fermions a † n , which create the state with phonon numbers n = 0, 1, 2, · · · , as reported by Abrikosov for the case of "spin" in the Kondo problem. 42 Using the pseudo-fermions, the displacement operator is represented bŷ In this representation, the Holstein-Anderson model with the Coulomb repulsion of the electrons U = 0 becomes where c † 0σ (a † n ) is the on-site conduction electron (pseudofermion) creation operator with spin σ, ǫ k and v characterize the conduction electron dispersion and the hybridization, respectively, and c † kσ is the creation operator of the conduction electron with the wave vector k and spin σ. We have introduced a Legendre multiplier λ in order to prohibit double occupancy of the pseudo-fermions. Next, we explain a self-consistent treatment of this impurity model. The method explained here is similar to the self-consistent ladder approximation in the Coqblin-Schrieffer model with crystalline-electric-field states of Ce impurities 45,46 and non-crossing approximation (NCA) of Anderson model. 43,44 We define the Matsubara Green functions of a n and c 0σ in the imaginary time τ as follows: where T τ is the time-order operator, and we omit the spin dependence in G(τ ) hereafter.
The non-interacting Green functions A 0 nm (iω n ) and G 0 (iω n ), are given as where where ω n is the fermionic Matsubara frequency ω n = (2n + 1)πT , In order to take into account the electronphonon coupling, we consider the self-energy diagrams Σ nm A (iω n ) and Σ G (iω) shown in Figs. 2(a) and 2(b), respectively, as: Here, χ GA is the solution for the diagrammatic equation given in Fig. 2(c), which is reduced to the summation of the infinite series of the ladder diagram χ 0 GA constituted of G and A. This is analytically given by where ν l = 2πlT is a bosonic Matsubara frequency. In matrix form, eq. (29) is expressed aŝ The explicit form of χ 0 GA is given by Note that the first-order term in g in Σ A and Σ G vanishes, together with the term 1 2 in the interaction term in eq. (21).
Using Σ nm A and Σ G together with Dyson's equation: we can obtain the self-consistent Green functions. In actual calculations, it is better to treat the retarded Green functions Σ nm A (z + i0) and Σ G (z + i0). Carrying out the analytic continuations, we can obtain the set of equations for Σ nm A (z + i0) and Σ G (z + i0), as shown in the Appendix (we hereinafter write z + i0 as simply z). The phonon Green function D(τ ), which is the most crucial quantity in the present paper, is given by where we ignore the vertex corrections in the last line. In this approximation, the phonon Green function in Matsubara frequency ν l is given by The main feature of the approximation scheme explained in this subsection is to take into account "phonon fluctuations". As such, the results obtained in this method might overestimate the fluctuations because we ignore the vertex corrections in eq. (34). Despite this, the present approximation can visualize the essential points of the low-energy phenomena in the present electronphonon system.

Numerical results
In this subsection, we show the numerical results of the Holstein-Anderson model (21). For the conduction electron, we assume the Gaussian density of states ρ(z): which corresponds to an infinite dimensional hyper-cubic lattice with the nearest neighbor hopping t. We set t = 0.2 and Ω E = 0.01 throughout this subsection. The hybridization width Im∆(z) is given by using ρ(z) as −Im∆(z) = πv 2 ρ(z), where v is a constant. The real part of ∆(z) is calculated by the Kramers-Kronig relation where the integral is taken as its principle value.
In the numerical calculations below, we introduce the cutoff number for the pseudo-fermion: N cut , i.e., a 0 , a 1 , · · · , a Ncut−1 . Due to this restriction, we can calculate the dynamical quantities only in T < ∼ (N cut − 1)Ω E . However, as shown below, this restriction does not become serious in the temperature regions of interest.

-Electron self-energy -
We show the result of the imaginary part of the electron self-energy ImΣ G in Fig. 3. As the temperatures decrease, a dip structure is developed around |z| < Ω E . This simply represents the fact that an electron with an energy below Ω E cannot emit the Einstein phonon with energy Ω E . Thus, the temperature dependence of ImΣ G (z) for z ∼ 0 becomes the activation-type, ∝ exp(α/T ), where α is a constant on the order of Ω E , as shown in Fig. 4.

-Phonon Green functions -
In Figs. 5 and 6, we show the frequency and temperature dependence of the phonon Green function −ImD(z) in high- (Fig. 5) and low- (Fig. 6) frequency regions. The position of the main peak around Ω E decreases as temperatures decrease, as shown in Fig. 7. The peak width also decreases at low temperatures, making the peak sharper. This main peak corresponds to the excitations a 0 ↔ a 1 , a 1 ↔ a 2 , · · · in the pseudo-fermion picture. At low temperatures, below T ∼ 0.3Ω E , a low-lying excitation becomes prominent, as shown in Fig. 6. This low-energy soft mode arises from the development in the off-diagonal elements of pseudo-fermion Green functions A 01 (0),Ā 12 (0), · · · , andĀ nn (0) for n = 0.
Note that the slope at the origin maintains its linear frequency dependence, i.e., −ImD(z)/z = const. near z = 0, as shown in Fig. 8. Thus, this peak is not generated by fictitious increases of the dynamical susceptibility Imχ(z)/z that frequently arise in NCA at  low temperatures. 47,48 The temperature dependence of this slope is shown in Fig. 9. At intermediate tempera- , the temperature dependence of the slope is well fitted by the activation-type function, exp(α ′ /T )/T , where α ′ is a constant on the order of Ω E . At low temperatures (T /Ω E < ∼ 0.1), the data deviates from the exponential dependence and exhibits 1/T 2 behavior. These behaviors are discussed below.
In order to clarify the origin of the temperature dependence of ImD(z), it is useful to consider ImΣ nm A (0). In the lowest order ladder diagram of χ GA , ImΣ nm A is written as Here, n B is the Bose function and Imχ G is given by × ImG(y + z)ImG(y).
At low temperatures (T < ∼ 0.1Ω E in Fig. 9, referred to herein as region I), taking into account the fact that the temperature dependence in lim z→0 ImD(z)/z comes from the contribution of the low-energy peak in ImĀ nm , we can replace ImĀ n ′ m ′ in eq. (39) by ∝ δ(y). This yields ImΣ nm A (0) ∝ T . On the other hand, at higher temperatures (0.1Ω E < ∼ T < ∼ Ω E in Fig. 9, referred to herein as region II), the strength of ImĀ nm (0) is weak and the temperature dependence comes from where Ω ′ E is the renormalized Einstein energy. Although Ω ′ E is weakly temperature dependent, for simplicity, we regard Ω ′ E as a temperature independent parameter. These observations yield ImΣ nm For the phonon Green function, eq. (A·4) is reduced, at low temperature and frequency with z ≪ T to where we have used − l dx π e −x/T ImĀ ll (x) ≃ 1 for T < ∼ Ω E . To simplify eq. (41) further, we omit the indices n, m, n ′ , and m ′ . Taking into account the fact that the dominant contribution comes from y = 0 in the integrand, and −ImĀ(y) ∼ (γ/2)/(y 2 + (γ/2) 2 ) with γ = −2ImΣ A (0) (this is of course too simplified an approximation), we can estimate −ImD(z)/z as These estimations reproduce the numerical results shown in Fig. 9. The expression eq. (42) is expected to be a rele-

Frequency dependence on sound velocity
In order to estimate the frequency dependence on the sound velocity correction with the use of the Holstein-Anderson model, we note that the high-energy spectra ImD(z) does not contribute to the frequency dependence in ReD(z) for z ≪ Ω E through the Kramers-Kronig transformation. The high-frequency part only gives weak monotonic temperature dependence. Thus, we define ReD low (z) by introducing the cutoff parameter Λ < Ω E as: From the discussions in §2, we can estimate the effect of the optical phonon on the sound velocity. We approximately replace the second term in eq. (14) by the phonon Green function obtained in the previous single-site problem: Although the present results cannot be directly compared to the experimental data in the case in which the parity of the optical mode is −1, it is expected that their qualitative behaviors do not differ greatly. In Fig. 10, we show the temperature dependence of ReD low (z) with different frequencies with Λ = 0.0035 = 0.35×Ω E and g = Ω E . Although the choice of Λ includes a certain arbitrariness, the essential feature of ReD low does not change if we take Λ to be smaller than ≃ 0.5Ω E . In the case of harmonic coupling, ReD low (z) is directly related to the sound velocity and the EC (we cannot discuss the anisotropy in EC's, see §4). For this case, the frequency-dependent term C(z) in the EC's is simply given by Clear hardenings, the positions of which shift to higher temperatures as the frequency increases, can be observed in Fig. 10. This feature is consistent with the phenomenological expression (1). In addition to ReD low , there are contributions from the high-energy spectra: ReD high (z) ≡ ReD(z)−ReD low (z) and the anharmonicity of acoustic phonons, both of which are essentially frequency independent. The elastic constant observed in the experiments is basically given by these two contributions unless we take into account electric contributions. The former makes the elastic constant soften at low temperatures and low frequency. On the other hand, the latter makes the elastic constant harden. We show the temperature, g (coupling constant), and frequency dependence of ReD low and ReD high in Figs. 11(a)-11(d). The magnitude of ReD low increases as g increases compared to the ReD high . Although, in Fig. 11(a), there exist shoulders as a function of temperatures in ReD low (z), as shown in the inset of Fig. 11(a), their magnitudes are quite small relative to ReD high . On the other hand, these features disappear at larger g, as shown in Fig. 11(c). The shoulder in the parameters in Fig. 11(c) can be seen at the higher frequency (not shown in Fig. 11(c)). In all of the data in Figs. 11(a)-11(c), softening behaviors are observed at ReD low (z) with small z. These features originate from the development of the soft mode being too strong in our calculation. More reliable methods are needed for the low-temperature region in order to discuss this point further.
In order to observe the temperature dependence of the shoulder in ReD(z), it is useful to differentiate ReD low (z) by T . In Fig. 12, we plot the temperatures T * , where |dReD low (z)/dT | takes a maximum value. Roughly speaking, T * corresponds to the midpoint of the shoulder. In Fig. 12, we show the results for three different parameter sets, keeping the renormalized Einstein frequency Ω ′ E the same for each set. As shown in Fig. 12, the frequency dependence of Ω ′ E /T * is well described by the condition zτ = 1 with τ = τ 0 exp(α ′′ /T ). Here, as before, α ′′ ≃ 0.7Ω ′ E is a parameter on the order of Ω ′ E . Interestingly, the estimated values of τ −1 0 are approximately α ′′ × 1/100. This is in good qualitative agreement with the experimental values τ −1 0 = 0.006α ′′ − 0.1α ′′ for filled skutterudites and clathrate compounds. We can find that τ −1 0 becomes small as g decreases. This arises from the fact that the soft mode amplitude becomes smaller and its frequency region moves to the lower side as the electron-phonon interaction g becomes small. This feature is sensitive to the activation energy α ′′ . A small change in α ′′ can affect the magnitude of τ 0 . This is the reason why we show the data with the same Ω ′ E values in Fig. 12.

Discussions
Thus far, we have restricted our study to the examination of a system with one component phonon. Here, we For this material, our results suggest that this T 2g mode strongly interacts with conduction electrons, and the observed UD in C 44 can be naturally interpreted according to the discussion in §2.1. For the optical phonon of the rare-earth dominated mode (T 1u ) in filled-skutterudites, there are three degrees of freedom, namely x, y, and z. In order to discuss the anisotropy in the UD, it is desired to take into account these aspects of the system. For example, we must estimate the wave number dependence in the electron-phonon coupling g, which has a crucial role for the anisotropy in the UD. The results obtained in the present paper can be thought of as simplified but essential for describing the frequency dependence in the EC. Experimentally, there is an inconsistent result regarding the anisotropy of the UD in the different samples of PrOs 4 Sb 12 . 51 This discrepancy should be investigated by further experimental efforts.
At this stage, the question as to whether the rare-earth T 1u mode is the relevant mode for the realization of the UD in filled skutterudites arises. It might be interesting to determine whether the temperature dependence of the E g modes in filled skutterudites is unusual, e.g., whether softening occurs. It is important to carry out high-resolution Raman scattering experiments to clarify this point. Recently, Ogita et al., observed the secondorder Raman spectra in filled skutterudites (especially Sb compounds) as an anomalous property of the rareearth mode. 11 In the case in which the relevant mode is actually the rare-earth T 1u mode in filled skutterudites, the mechanism of the anisotropy in the UD might be more complicated, as mentioned above. The T 1u optical phonon should be investigated as a relevant optical phonon for the UD as a further theoretical study. Experiments in unfilled skutterudites with electric bands similar to, e.g., LaOs 4 Sb 12 are also desired in order to identify which mode is important for the UD. Applying this theory to optical phonons with even parity and comparing the energies of the optical phonons estimated by Raman scattering and the activation energy obtained by the ultrasound experiments for each materials will provide various interesting properties.
In filled skutterudites, anomalous phonon contributions are observed in various quantities, especially in ROs 4 Sb 12 (R=Pr, La, Sm, etc.). One reason why the effect is prominent in ROs 4 Sb 12 may be that the size of the Sb 12 cage is the largest among existing filled skutterudites, including the rare-earth ion R. Another reason may be that the conduction electrons near the Fermi level in Os-compounds consist of a molecular orbital of the cage (A 1u ) and the d-electron in Os. Unlike the Rucompound, in which the levels of 4d electrons are deep compared to that of 5d or 3d electrons, so that their components near the Fermi level are negligible, the contributions from the d-electrons are much larger. It is thus expected that d-electrons play an important role in Os compounds. The latter point, however, requires more sophisticated analysis, because Iwasa et al., recently observed the quasi-elastic peak even in Ru compounds. 13 The observed width of the quasi-elastic peak is too large for UD to be realized in the MHz frequency range. In order to clarify the reason why the UD has been observed in Os-compounds but not in Ru-compounds, it may be necessary to carry out experiments in a wider frequency range. Note that if we assume the dominant electronphonon coupling between the A 1u molecular orbital and the T 1u optical phonon is local and linear in the displacements, it is impossible for the A 1u molecular orbital to interact with the T 1u optical phonon due to symmetry.
Recently, nuclear magnetic resonance (NMR) experiments on KOs 2 O 6 52 and LaOs 4 Sb 12 53 were carried out. The results indicate that the relaxation is due to the quadrupolar coupling between the nuclei and the ionic motion (through a direct process) at K and La sites. In this case, 1/(T 1 T ) is given by a phonon spectral weight: 52 1/(T 1 T ) ∝ −ImD(z)/z. The results are interpreted by a strongly damped oscillator with an activation-type relaxation time. Although the energy scale of the ultrasonic and NMR experiments is of the same order as the magnitude (∼ MHz), the results show no frequency dependence in 1/(T 1 T ). This might be due to the difference between the q = 0 and q-integrated spectrum, or the difference of the phonon mode in the ultrasonic and NMR experiments, as mentioned above (the ultrasonic experiment for KOs 2 O 6 has not yet been performed). Based on Fig.  6, it is expected that 1/(T 1 T ) has frequency dependence. In the case in which the relaxation time is indeed an activation-type, i.e., ∝ exp(E ∆ /T ), the small change in z does not indicate the same order change occurs in T , due to the exponential dependence. The precise origin of this problem is not yet understood and further experimental and theoretical studies are required. 54 For the theoretical aspects, the approximated method used in the present paper becomes worse in the lowtemperature regime, because − lim z→0 ImD(z)/z diverges as T −2 . This might cause the unphysical development of the soft mode. Intuitively, the divergence must be stopped at a certain temperature. However, it is important for this method to be able to capture the temperature dependence of the −ImD(z) for z ≪ Ω E at the intermediate temperature regions, where UD is observed.
In our treatment of pseudo-fermion representations, it is straightforward to use higher-order interactions e.g., X 2 , X 3 , · · · as interaction parts of the Hamiltonian. This is because the same argument can be applicable by simply replacing the matrix elements X with X n . The model including these higher-order electron-phonon couplings becomes important for ions with very large displacements, such as the rare-earth ion in filled skutterudites. The variations of hybridization between the felectron and conduction electron by these ionic motions should be taken into account in a realistic model for the further clarification of electron-phonon systems. For the anharmonicity of the potential for the local phonon, it is straightforward to include higher-order terms such as X 4 . Although the same line of discussion can be introduced in the present paper from an off-center potential, the presence of the off-center potential does not play a fundamental role with respect to our results.
Finally, we mention the relationship between the present theory and our previous studies. 15,16 In the present paper, we have discussed primarily the properties of a local phonon that interacts with conduction electrons. In our previous studies, however, we discussed the properties of electrons interacting with off-center configurations of an ion, which is an effective theory for the strong coupling limit of an electron-phonon system. 55 Further theoretical studies are needed in order to unify these aspects.

Summary
In conclusion, we have investigated the acoustic and optical phonon spectra in a strongly coupled electronphonon system at finite temperatures. We have applied the self-consistent ladder approximation, which was used successfully in the Coqblin-Schrieffer model with crystalline-electric-field states of Ce impurities, to the electron-phonon system, and have determined that the low-energy peak in the optical phonon spectral function develops at low-temperatures and causes the ultrasonic dispersion observed in filled-skutterudites and clathrate compounds. The temperature dependence of this peak shows activation type behaviors at the intermediate temperature regions, leading to the Debye-type formula for the ultrasonic dispersion with activation-type relaxation time. These results qualitatively explain the experimental results of the ultrasonic region observed in the filled skutterudites and related compounds. Science (JSPS). One of the authors (K.H.) is also supported by a Research Fellowship for Young Scientists from JSPS.

Appendix: Self-consistent Equations for the Self-Energies and Phonon Green Function in a Real Frequency
In the Appendix, we list the self-consistent equation for the imaginary part of the retarded self-energies and the form of the retarded phonon Green function. Self-Energies: ImΣ nm A (z) = 2g 2 dy π n F (y)ImG(−y) × n ′ m ′ X nm ′ Imχ m ′ n ′ GA (y + z)X n ′ m (A·1) and ImΣ G (z) = (1 + e −z/T )g 2 dy π e −y/T nmn ′ m ′ × ImĀ mn ′ (y + z)Imχ m ′ n GA (y)X nm X n ′ m ′ × − l dx π e −x/T ImĀ ll (x) with Imχ 0 GA nm (z) = dy π n F (y)ImG(y)ImĀ nm (z + y), (A·3) whereχ GA in eq. (A·2) is obtained using eq. (30). In eqs. (A·1), (A·2), and (A·3), n F represents the Fermi distribution function. The function with an overbar is defined asf (x) ≡ f (x + λ), taking λ → ∞ subject to the local constraint n a † n a n = 1. In eq. (A·2), we have divided the result, which is obtained by the analytic continuation of eq. (28), by n a † n a n . 44