between predictions from molecular dynamics simulations (a method valid at temperatures above the Debye temperature) ... A key feature of ionic materi...

1 downloads 0 Views 575KB Size

Journal

Phonon-Mediated Thermal Conductivity in Ionic Solids by Lattice Dynamics-Based Methods Aleksandr Chernatynskiy,‡ Joseph E. Turney,§ Alan J. H. McGaughey,§ Cristina H. Amon,¶ and Simon R. Phillpot‡,† ‡

Department of Materials Science and Engineering, University of Florida, Gainesville, Florida, 32611

§

Department of Mechanical Engineering, Carnegie Mellon University, Pittsburgh, Pennsylvania, 15213

¶

Department of Mechanical and Industrial Engineering, University of Toronto, Toronto, Ontario, M5S 3G8 be carefully taken into account. In particular, modiﬁcations to existing techniques for computing thermal conductivity may be required as the standard expression for the energy current, as used in Green–Kubo technique, is derived under the assumption that the interatomic interactions are shortranged.12 Electrostatic interactions lead to complexities because the presence of internal and/or macroscopic electric ﬁelds leads to the storage of electrostatic energy. The spatial dependence of this electrostatic energy leads to an eﬀective energy current that must be considered in analyzing thermal transport. Lattice dynamics methods have long been used to explore phonon properties of materials. The strengths and weakness of LD methods mirror those of MD. LD methods can account for the quantum nature of the vibrations in a natural and rigorous manner. Due to the nature of the harmonic approximation used in lattice dynamics, it is a low-temperature theory, but it allows the properties to be determined with some accuracy over a considerable range of temperatures, typically up to about one-third of the melting point.13 At higher temperatures, LD methods fail, ﬁrst in a quantitative manner and then in a qualitative manner, as they are based on the phonon properties at zero temperature. Fortunately, many ionic materials are characterized by high melting temperatures (for example, the melting temperature of MgO is 3000 K), and therefore “low temperature” for these materials can mean temperatures as high as 1000 K, which covers a wide range of their applications. The use of LD methods for determining thermal conductivity of non-ionic systems in which phonon processes dominate is well-documented. Examples include the classic test case of solid argon,9,11 diamond,14 and silicon and germanium.7,15 The thermal conductivity of carbon nanotubes and more recently that of graphene have also been explored.16 The possibility of signiﬁcantly lower computational loads in the case of small primitive cells makes LD methods very appealing for the determination of the thermal conductivity in ionic systems. In contrast, application of LD to the thermal-transport properties of ionic materials has been limited.17–19 In one study in alkali halides,17 a nearest neighbor approximation was used for the interatomic interactions, thereby completely neglecting the long-range nature of the electrostatic interaction. Similarly, Tang and Dong18 used a supercell approach for the computation of the force constants, but needed to correct their optical phonons frequencies with a spring model to give the physically observed splitting between the longitudinal optical (LO) and transverse optical (TO) phonons at the Γ-point. The treatment of electrostatic interactions in MD simulations of ionic systems has a rich history. A number of

Phonon properties predicted from lattice dynamics calculations and the Boltzmann Transport Equation (BTE) are used to elucidate the thermal-transport properties of ionic materials. It is found that a rigorous treatment of the Coulombic interactions within the harmonic analysis is needed for the analysis of the phonon structure of the solid, while a short-range approximation is suﬃcient for the third-order force constants. The eﬀects on the thermal conductivity of the relaxation time approximation, the classical approximation to the phonon statistics, the direct summation method for the electrostatic interactions, and the quasi-harmonic approximation to lattice dynamics are quantiﬁed. Quantitative agreement is found between predictions from molecular dynamics simulations (a method valid at temperatures above the Debye temperature) and the BTE result within quasi-harmonic approximation over a wide temperature range.

I.

T

Introduction

thermal-transport properties of ionically bonded materials have attracted signiﬁcant interest, fueled by applications ranging from thermal barrier coatings,1 to uranium dioxide as a nuclear fuel,2 to novel thermoelectric applications.3 In the geophysics community, the thermaltransport properties of ionic crystals under extreme conditions are of fundamental importance for theories of planetary interiors.4 There are a number of atomistic techniques available for the simulation of phonon-mediated thermal transport in bulk materials. These include equilibrium molecular dynamics (MD) simulations using the Green–Kubo method, various ﬂavors of non-equilibrium MD, (see [5,6] for a review and comparison between MD techniques) and lattice dynamics (LD) methods.7–11 In general, MD techniques suﬀer from system size eﬀects and convergence issues5,6 that require relatively large computational loads to resolve. In addition, MD techniques are classical and therefore strictly applicable only in the high-temperature limit (i.e., above the Debye temperature). A key feature of ionic materials is the presence of the long-range electrostatic interactions whose inﬂuence needs to HE

D. Smith—contributing editor

Manuscript No. 29595. Received April 13, 2011; approved June 17, 2011. The work of AC was supported by DARPA and the NSF Materials World Network Project under Grant No. DMR-0710523. † Author to whom correspondence should be addressed. e-mail: [email protected]ﬂ.edu

3523

3524

Journal of the American Ceramic Society—Chernatynskiy et al.

algorithms have been developed to perform the Coulomb sums20 (and references therein), including the Ewald summation, fast multipole algorithms, and particle–particle-mesh methods. More recently, a direct summation method has been developed21; it has the merits of computational simplicity and speed. The main idea of this method is that electrostatic interactions, while mathematically conditionally convergent and long-ranged, are physically rather short-ranged in solid-state settings. Earlier computational diﬃculties (strong oscillations of the total energy as a function of cutoﬀ radius, for example) were attributed to charge oscillations within the cutoﬀ sphere. By placing compensating charges on the cutoﬀ sphere, a method for the summation of the electrostatic interactions was derived, which is shortranged and reproduces the energy and structure of ionic compounds with good accuracy. The direct summation has been used to study such diverse phenomena as ferroelectricity in PbTiO3,22 defects and dislocation dynamics in MgO,23 and the liquid–vapor water interface.24 It has also been used for MD simulations of thermal conductivity in various systems.25–29 In this paper, we address the applicability of the lattice dynamics method for computing the bulk thermal conductivity on the technologically important ionic systems of MgO, SrTiO3, and UO2 by comparing the results of MD simulations and LD calculations. The range of applicability of the direct summation method for thermal conductivity calculations is also determined. In addition, we include a shell model into the lattice dynamics calculations of thermal conductivity on the example of the KCl.

II.

~00

00

K~k ;n~0

k;nk ;n0

~00

00~0

K~kk;n;n

k ;n0

8 1 < X 0 00 rT ¼ ~ v~k;n U~k;n þ U~k ;n0 U~k ;n00 kB T : 0 0 00 00 @T ~ ~ ;n k ;n ;k ) ~0 0 ~00 00 00 1 ~ k ;n00 k ;n ;k ;n þ ... U~k;n U~k0 ;n0 U~k00 ;n00 K~k;n K~ ~0 0 þ k;n;k ;n 2 @f ~0k;n

(1) In Eq. 1, n labels the phonon branch index, ~ k is the reciprocal lattice vector of the phonon states; ~ v~k;n is the group velocity, and Λ is the equilibrium transition rate for a three-phonon process. The right-hand side represents the change of the probability distribution function of the phonon state as a result of phonon–phonon interactions. The transition probabilities Λ are computed from the Fermi golden rule by treating anharmonicity as a perturbation to the harmonic problem. In the current investigation, we consider cubic anharmonicity; therefore, only two types of processes are described: the splitting of a phonon into two phonons and the combination of one phonon with another phonon to produce a third phonon. Processes involving four or more phonons are not included in this analysis. The three-phonon processes result in the following expressions for the transitions probability

2p 2 C 0 00 ~~0~00 f ~0k;n f ~00 0 f 0 ~00 00 þ 1 k ;n k ;n h nn n ;kk k d hx~k;n þ hx~k0 ;n0 hxk~00 ;n00

2p 2 Cnn0 n00 ;k~~k0~k00 f ~0k;n f 0 ~0 0 þ 1 f 0 ~00 00 þ 1 k k ;n ;n h d hx~k;n hxk~0 ;n0 hxk~00 ;n00 ðh 2Þ3=2 1=2 1=2 1=2 ¼i x~k;n x~0 0 x~00 00 k ;n k ;n V1=2 X BIJK;k~~k0~k00 ~ e ~;n~ eJ;k~0 ;n0~ eK;k~00 ;n00 1=2 1=2 1=2 I;k IJK mI mJ mK

¼

Cnn0 n00 ;k~~k0~k00

(2)

where B is the matrix of Fourier-transformed third-order force constants, that is the third derivatives of the total potential energy with respect to atomic positions, I, J, K enumerate atoms in the primitive cell; mI is the mass of atom of species I, while ~ eI;k~;n is the eigenvector of the phonon with frequency x~k;n . The d-functions enforce conservation of energy. Solutions methods of the BTE include the relaxation time approximation11 and a full iterative solution.7–9 Once the solution of the BTE is obtained, Fourier’s law is used to compute the thermal conductivity through the expression for ~~ the heat current in terms of the distribution function (F k;n is a proxy for the distribution function perturbation: U~k;n ¼ ~ F~k;n rT)

Thermal Conductivity from Lattice Dynamics

Our implementation of the solution of the BTE is described elsewhere.8 Here, we only provide a brief summary of the pertinent aspects. The BTE for phonons is solved under the assumption of small deviations from thermodynamic equilibrium; that is, Fourier’s law is applicable. This assumption allows the linearization of the BTE with regard to small deviations U~k;n of the probability distribution function f~k;n from the equilibrium Bose–Einstein distribution f ~0k;n

¼

Vol. 94, No. 10

k¼

Z X n

hx~k;n

f~0k;n 1 þ f~0k;n kB T

i ~ ~ v~k;n ~ F~k;n dk

(3)

As we observe from this description, the thermal conductivity depends on the second derivatives of the total potential energy with respect to atomic positions by deﬁning the phonons properties in the harmonic approximation, and on the third derivatives through phonon–phonon interactions. The total potential energy of an ionic solid can be computed from a standard classical model of the pair interactions: V¼

VCoul ij

þ

VBuck ij

1 qi qj rij ¼ þ Aij exp 4pe0 rij qij

!

Cij Dij 8 r6ij rij (4)

The ﬁrst term in Eq. 4 is the long-ranged Coulombic interactions. The second, third, and fourth terms are an empirical short-range repulsive term and the van der Waals terms; together, these three terms are known as the Buckingham potential. Note that the short-range part of the interaction does not need to come from an empirical model at all, and could be computed from the ﬁrst principles, making it possible to compute the thermal conductivity ab initio.14,18,19,30 The long-range nature of the electrostatic interactions, however, still needs to be taken into account. As this is the eﬀect of interest in the present study, we employ Eq. 4 for our investigation. Table I gives the values of potential parameters used. For the description of the core-shell model and the related eﬀects, see Section V.

III.

Thermal Conductivity of Ionic Materials by the BTE

We apply the BTE for computing the thermal conductivity of four model materials: rocksalt-structured MgO and KCl (the latter is an example of using a core-shell model and will

October 2011

3525

Elucidation of Thermal-Transport Properties of Ionic Materials Table I.

Potential Parameters for MgO, SrTiO3, and UO2 Used in this Work

MgO31 qMg = 1.7, qO = 1.7 SrTiO332 qSr = 2, qTi = 4, qO = 2 UO233 qU = 4, qO = 2 KCl34 qcK ¼ 6:9659, qsK ¼ 5:9659, j = 484.83 eV/A˚2qccl ¼ 1:6532, qscl ¼ 2:6532, j = 31.334 eV/A˚2

Mg–O O–O Sr–O Ti–O O–O U–O O–O K–K Cl–Cl K–Cl

be discussed in detail in Section V), ﬂuorite-structured UO2, and perovskite-structured SrTiO3. For MgO31 and UO22, the thermal conductivity was computed previously by non-equilibrium MD methods, while the thermal conductivity of SrTiO3 has not been previously investigated by simulation. All MD calculations (except those for KCl) take into account the size eﬀect associated with the ﬁnite simulation cell by extrapolating the thermal conductivity value to the inﬁnite system size.5 For computational simplicity, our current computational implementation of the solution of BTE is limited to tetragonal cells. We therefore use a four-atom non-primitive unit cell for MgO and a six-atom unit cell for UO2. The ﬁve-atom primitive cell for SrTiO3 has the cubic perovskite structure; we do not take into account the antiferrodistortive structure of SrTiO3 at low temperatures. A suﬃciently dense mesh in the reciprocal space is used for each material: 13 9 13 9 13 for MgO and UO2, 9 9 9 9 9 for SrTiO3. Using denser meshes produces thermal conductivity values that diﬀer by 1% or less from the values reported here. As the previous MD calculations of the thermal conductivity used the direct summation to account for the electrostatic interactions, the same method is used for this comparison. The direct summation method is characterized by two parameters: a damping factor g and a cutoﬀ radius Rc.21 A large damping factor corresponds to an aggressive approximation to the true Coulombic interaction and corresponds to a short cut-oﬀ radius, and vice versa. As g goes to zero, the full Coulombic potential is recovered. Parameters for the direct summation simulations are those used in MD simulations: g = 1.5a1and Rc = 2a, where a is a lattice constant of a cubic cell for a given material. All the calculations of the thermal conductivity are performed in the quasi-harmonic approximation, i.e., phonons frequencies are treated as being only volume-dependent. In a fully self-consistent QHA, the volumes are determined from minimization of the Gibbs free energy as determined from the lattice dynamics calculations; here, to allow direct comparison with the MD simulation results, we use the temperature-dependent zero-pressure volumes determined by MD for these potentials. As MD simulations sample phonons occupations in accord with classical Maxwell–Boltzmann distribution, we used the same distribution in our lattice dynamics calculations rather than the quantum Bose–Einstein distribution, unless otherwise noted, thereby allowing direct comparison between MD and LD methods. First, we compare the results of applying the full iterative solution of the BTE to estimate the thermal conductivity with the MD results, see Fig. 1. The MD and LD results match very well for both SrTiO3 and UO2 even at very high temperatures (~2000 K for UO2). Such good agreement at high temperature is rather unexpected, and indicates that the cubic term in the energy expansion is suﬃcient to describe the thermal conductivity. A similar conclusion was reached previously in the studies of solid argon.8,9 In MgO, however, the LD calculations yield a thermal conductivity that is larger than the MD results.

Aij (eV)

ρij (A˚)

Cij (eV·A˚6)

Dij (eV·A˚8)

929.69 35 686.18 1400.0 754.2 22 764.3 1761.775 9547.96 6.172 9 109 3361.0 4852.0

0.2991 0.2010 0.3500 0.3879 0.1490 0.3564 0.2192 0.1085 0.3564 0.3080

0.0 32.0 0.0 8.986 43.0 0.0 32.0 36.23 193.9 80.65

0.0 0.0 0.0 0.0 0.0 0.0 0.0 28.74 269.1 87.95

Fig. 1. Comparison between LD and MD methods for thermal conductivity using Maxwell–Boltzmann statistics. (Typical statistical error in MD methods is about 10%).

Fig. 2. Thermal conductivity of MgO, SrTiO3, and UO2 as computed by LD in comparison with the experimental data.

We compare our calculations with the experimental data for these materials in Fig. 2. For this comparison, we employed the Ewald summation for the electrostatics and the Bose–Einstein distribution for the phonon occupation numbers. Experimental data for MgO are obtained from [35], for SrTiO3 from [36], and for UO2 from [37]. The LD values for the MgO thermal conductivity are very close to the experimental values, a measure of the high ﬁdelity of this potential for this property. By contrast, the LD values are signiﬁcantly greater than experimental for the SrTiO3 (by factor of 2) and UO2 (by factor of 2.5). Such disagreement is common for

3526

Journal of the American Ceramic Society—Chernatynskiy et al.

Vol. 94, No. 10

empirically ﬁtted potentials, for which anharmonic eﬀects have not typically been part of the training set for the determination of the potential parameters. Similar eﬀects are also manifested in the thermal expansion, which is typically underestimated by empirical potentials, again due to an underestimate of the anharmonicity. It should be noted that our calculations consider ideal single crystal, while experimental data are obtained on crystals with some concentration of defects; such defects decrease thermal conductivity. In particular, UO2 experimental data were obtained for 95% dense material, while SrTiO3 data were from Nb-doped samples. Why should the LD method reproduce the MD results so well for two materials, but not for the third? We speculate that for MgO, the eﬀects of higher order terms in the energy expansion (i.e., processes involving more than three processes) are also important. Our lattice dynamics calculations lack these higher order terms, resulting in the overestimation of thermal conductivity. Next, we assess the applicability of the widely used relaxation time approximation for calculating thermal conductivity in these systems. Figure 3 presents the thermal conductivity as a function of the number of iteration steps in the iterative procedure (see [8] for more details). The ﬁrst step is formally identical to the relaxation time approximation (RTA). The RTA provides computational time savings: if reaching full solution takes N steps, a full solution takes N times more computer time than the relaxation approximation. As the thermal conductivities of the three materials are quite diﬀerent, we plot them normalized to their respective converged thermal conductivity value. From Fig. 3, we can observe that the RTA for UO2 yields 97% of the full solution; in SrTiO3, the RTA gives about 99%, while for MgO, the RTA gives only 93% of the converged value. Underprediction of the thermal conductivity by the RTA by a similar amount was previously found for argon9 and silicon.7 In diamond, however,14 the RTA underprediction is from 30% to 80% in the temperature range from 1000 to 200 K. These results can be understood in the following way. The more harmonic a solid is (i.e., the higher the thermal conductivity), the less the various phonon modes interact with one another, and the longer the relaxation times; that is, the phonon modes remain out of equilibrium longer. The premise of the RTA approximation, on the other hand, is that each individual mode relaxes to equilibrium while all the other modes are in equilibrium. Therefore, if the relaxation times are long, such an approximation is not justiﬁed and proper accounting of the interactions between the phonon modes that are out of

equilibrium is important. Interactions between phonon modes at low temperatures are weaker; hence, the performance of the RTA becomes progressively worse as the temperature decreases. Overall, UO2 and SrTiO3 are represented by the RTA very well, while for MgO, the full iterative solution is required. We recall that, in order for a one-to-one comparison against MD simulations, all the results discussed so far were obtained by using classical statistics for the phonons occupation numbers. However, for a correct accounting of quantum eﬀects, Bose–Einstein statistics should be used. In the limit of high temperature (roughly, above the Debye temperature), both statistics must produce the same results. As temperature is reduced below the Debye temperature, quantum eﬀects start to play a role, and the Bose–Einstein statistics lead to two competing eﬀects on the thermal conductivity. One eﬀect stems from the fact that the transition rates for phonon processes are proportional to the occupation numbers of the phonons involved. Details were given previously by Turney et al.38 for the example of silicon. The net result is that low-frequency phonons have smaller scattering rates, while high-frequency phonons have larger scattering rates, in quantum statistics when compared with classical. If lowfrequency phonons dominate the thermal transport (as they do in Si), then quantum statistics produce a higher thermal conductivity than classical statistics around the Debye temperature. The second eﬀect stems from the fact that the speciﬁc heat of the individual phonon modes is overestimated by the classical statistics. This result tends to drive a quantum treatment to produce lower thermal conductivity. At temperatures just below the Debye temperature, the quantum statistics have the most eﬀect in reducing the number of high-frequency phonons. If high-frequency phonons make a noticeable contribution to the thermal conductivity, then they can compensate for, or even dominate, the eﬀect of the longer lifetimes of the low-frequency phonons. This is indeed the case for the ionic solids, as illustrated in Fig. 4(a). In this ﬁgure, we plot the ratio of the thermal conductivities computed using quantum and classical statistics as a function temperature, normalized to the Debye temperature for a given material. The results for MgO, SrTiO3, and UO2 are presented, and

Fig. 3. Normalized thermal conductivity as a function of the number of iteration step in the iterative solution (SrTiO3 is at 250 K, MgO, and UO2 are at 300 K).

Fig. 4. (a) Comparison of quantum versus classical statistics for MgO, UO2, SrTiO3, and Tersoﬀ-Si. (b) Peak in thermal conductivity and ~T3 decay as temperature goes to 0 K for SrTiO3.

October 2011

Elucidation of Thermal-Transport Properties of Ionic Materials

compared with Si described by the Tersoﬀ potential. The behavior of the silicon thermal conductivity is in agreement with the previous results (while Turney et al. in38 used Stillinger–Weber potential for silicon, the general conclusions are the same): above the Debye temperature, classical and quantum statistics agree with each other, while below the Debye temperature, the classical statistics underestimate the thermal conductivity. In contrast to this, the quantum and classical calculations for the thermal conductivity in MgO and SrTiO3 agree very well for temperatures far below the Debye temperature—all the way down to ~TD/3 (TD = 941 K for MgO, and TD = 513 K for SrTiO339). This somewhat unexpected agreement is due to the cancelation of the two eﬀects described above. It also explains why the thermal conductivity of the MgO as predicted by the MD31 is in good agreement with the experiment even at temperatures signiﬁcantly below Debye temperature. The situation for UO2 is more extreme. As one can observe from Fig. 4(a), the ratio of the thermal conductivities as computed with Bose–Einstein statistics and Maxwell– Boltzmann statistics ﬁrst decreases with decreasing temperature, passing through a minimum at about half of the Debye temperature, and then increases. At the minimum, the value of the thermal conductivity produced by the quantum statistics is 25% lower than the classical result. This indicates that the heat capacity reduction of the phonon modes in the quantum statistics is the dominant eﬀect around the Debye temperature. The above analysis describes the situation just below the Debye temperature. In the temperature range of ð1=30Þ.ðT=TD Þ.ð1=10Þ, the quantum statistics occupation numbers contain an exponential factor that eventually leads to an exponential dependence in thermal conductivity with temperature40 k AT s exp

b T

3527

potential,21 derivatives of the energy do possess distinctively long-range eﬀects. Notably, they cause the splitting of the longitudinal and transverse optical phonon modes. The direct summation method does not take the long-range nature of the derivatives into account, whereas the reciprocal space term in the Ewald summation does. As a result, the phonon band structure determined from the direct summation does not reproduce this splitting correctly. The phonon dispersion relationship for SrTiO3 is illustrated on the Fig. 5(a). It shows how the highest optical phonon branch changes as the damping parameter in the direct summation changes. This optical mode produces a macroscopic polarization, and therefore is susceptible to the long-range eﬀects of the electrostatic interaction. One can observe that its frequency increases with decreasing damping factor. However, the direct sum never reproduces the LO/TO splitting at the Γ point, leading to a mode with very large group velocity. Figure 5(b) shows the phonon density of states. Upon decreasing the damping factor, the DOS cut-oﬀ approaches the exact result, with only minor diﬀerences in the overall density of states. The diﬀerences in the phonon band structure are the largest for the high frequency phonons that do not transport much heat themselves, but provide states for the acoustic phonons to scatter into, as well as acting as a source of acoustic phonons. We can assess how the parameters of the direct summation inﬂuence thermal conductivity for our model systems. Here, we choose a very large cut-oﬀ radius for the interactions, 24 A˚, while we vary only the damping factor.

(5)

In classical statistics, this exponential factor is absent, and therefore the thermal conductivity increases only as a power law of temperature; hence, in this region, classical statistics signiﬁcantly underpredicts the thermal conductivity. Finally, at even lower temperatures, ðT=TD Þ.ð1=30Þ, boundary scattering limits the free path to the sample size. The only temperature dependence remains in the phonons heat capacitance and therefore, as calculated with quantum statistics, the thermal conductivity exhibits a peak and a ~T3 decay as the temperature approaches absolute zero. As the classical heat capacitance is temperature-independent, the classical thermal conductivity approaches a constant value at absolute zero. So yet again, the classical result overestimates the thermal conductivity. Boundary scattering might be considered within the RTA by introducing a free path of a ﬁxed length (simulating the sample size) and combining it with the individual phonon–phonon scattering length by the Matthiessen rule. Figure 4(b) demonstrates the ~T3 behavior and a peak in thermal conductivity in quantum statistics for the example of the SrTiO3. We also note that one needs to consider very dense k-meshes to correctly reproduce the thermal conductivity at this low temperature as the major contribution to thermal conductivity comes from phonons with the very small wave vectors (i.e., near the zone-center).

IV.

Electrostatic Interactions: Direct Summation versus Ewald Summation

In this section, we address the applicability of the direct summation method for the treatment of the electrostatic interactions in calculations of thermal conductivity. While the energy of an ionic system appears to be approximated very well using a short-ranged approximation for the electrostatic

Fig. 5. (a) Dispersion relation for the highest optical mode in SrTiO3 for various damping factors in the direct summation. (b) Phonon density of states. (c) Diﬀerence in PDOS for diﬀerent direct sum parameters.

3528

Journal of the American Ceramic Society—Chernatynskiy et al.

As discussed above, thermal conductivity calculations require second and third derivatives of the total energy with respect to atomic positions. For the Ewald summation, second derivatives are directly computed in reciprocal space to take into account the well-known non-analyticity of the dynamical matrix at the Γ point.41 Third derivatives are computed in real space within a deﬁned cut-oﬀ radius, and then Fourier-transformed. The eﬀect of the third derivative cutoﬀ is investigated in the following section. All the other parameters for thermal conductivity calculations are the same as in the previous one. The results of these calculations for UO2, SrTiO3, and MgO are presented in Fig. 6. Large values of the damping factor g, i.e., small values of g1, which correspond to a shorter range in the direct summation and, thus, a more aggressive approximation of electrostatic interactions, lead to an underestimation of the thermal conductivity. Upon decreasing the damping factor, the thermal conductivity increases monotonically at ﬁrst, then oscillates around a saturated value. For the three materials considered here, these asymptotic values predict thermal conductivities that are 10%–20% higher than those determined with the exact Ewald summation result. To understand these results, we draw attention to the density of states and phonon band structure at large damping factors (Fig. 5). There are a number of noteworthy diﬀerences. First, the DOS demonstrates that the direct summation in general reproduces the low-energy part of the spectrum (0–10 THz) reasonably well, but that the cut-oﬀ frequency is signiﬁcantly lower than that of the spectrum computed with Ewald summation (20 THz vs. 30 THz). This smaller cutoﬀ completely removes the relatively small contribution to the thermal conductivity from the high-frequency modes. As the

Fig. 6. Thermal conductivity at 500 K for SrTiO3 (a), and 300 K for MgO (b) and UO2 (c) as computed by the Ewald summation (solid line) and direct summation (solid circles) as a function of the direct summation damping factor.

Vol. 94, No. 10

number of phonon modes is determined by the number of atoms in the system, this reduction in the number of highfrequency phonons means that there is a correspondingly higher phonon density of state in the 10–20 THz frequency range. This results in more scattering channels available for phonons and therefore the thermal conductivity is further reduced. These two eﬀects explain the underestimation of the thermal conductivity by the direct sum with large damping factor. Second, the phonon band structure presented in Fig. 5(b) and (c) helps us to understand the behavior of thermal conductivity at small damping factors. The strong dispersion in the optical modes means that these phonons have very large group velocities, and thus can contribute signiﬁcantly to the total thermal conductivity. Indeed, there are peaks in the frequency-dependent contribution to the thermal conductivity at 24 and 26 THz (Fig. 7). Upon further decrease in the damping factor, the group velocity of the LO mode increases near the Γ-point further, tending toward a further increase in the thermal conductivity. However, it becomes ﬂatter and ﬂatter away from the Γ-point, where it resembles more and more the true state, tending to decrease the thermal conductivity. The competition between these two eﬀects produces the oscillations in thermal conductivity as a function of damping factor. We conclude that the direct summation produces thermal conductivities within ~40% of the Ewald summation results, with the direction of the error depending on the damping parameter: a too aggressive approximation (too large damping factors) will underestimate thermal conductivity, while a too small damping factor will overestimate it. The ﬁrst case is more relevant for the actual molecular dynamics simulations, as a large damping factor means a shorter range for the direct summation and therefore lower computational load. Finally, we address the question of the range of the third derivative. As mentioned previously, third derivatives are computed in real space within a certain cut-oﬀ radius and then Fourier-transformed. Fortunately, they decay with distance as 1/r4; therefore, their inﬂuence should be shortranged. This is indeed the case, as demonstrated in Fig. 8, which plots the thermal conductivity as a function of the cutoﬀ radius for the third derivatives calculations. One can observe clearly that third derivatives are important out to about the third-nearest neighbors, at which point, the thermal conductivity is saturated to within 2%. It is interesting to note that if only nearest neighbors are taken into account, thermal conductivity is very large (which is consistent with the limited phonon–phonon scattering produced by only nearest neighbors third derivatives). However, with the increased range, thermal conductivity approaches limiting value from below: larger interaction range produces larger thermal conductivity.

Fig. 7. Contribution to thermal conductivity from phonons of given frequency as computed by direct summation and the Ewald method in SrTiO3.

October 2011

Elucidation of Thermal-Transport Properties of Ionic Materials

Fig. 8. Normalized thermal conductivity in UO2, MgO, and SrTiO3 as a function of cut-oﬀ radius (in units of lattice constant for a given material) for the third order derivatives.

V.

Eﬀects of the Core-Shell Model

Potentials of the type described by Eq. 4 are rigid-ion potentials, reﬂecting the fact that atomic polarizability is not taken into account. However, optical phonon modes produce macroscopic electric ﬁelds; hence, contribution from the atomic polarizabilities can be important. Their neglect results in overestimated frequencies of the optical phonon modes and this is one of the reasons why rigid-ion potentials typically overestimate the thermal conductivity of the ionic solids. A simple approach to include polarizability in calculations, the core-shell model, was proposed by Dick and Overhauser.42 In this approach, each atom is modeled as a core that represents the nucleus and core electrons, and a massless shell that represent the valence electrons. The core and shell of a single atom do not interact coulombically; rather, they interact through a harmonic potential characterized by the spring constant j. The short-range interaction is considered to act between shells. This model signiﬁcantly improves the dielectric constants and properties of the optical phonon modes. From the perspective of MD simulations, however, the shell model is computationally more expensive than the rigid-ion model: shells need to be relaxed to their equilibrium position at every time step (this is an analogy of Born–Oppenheimer approximation in ﬁrst principles MD). An alternative approach ascribes a mass on the shells and runs a straightforward MD simulation, but in that case, the shell mass has to be small, requiring a smaller time step.43 As a result, there are very few simulations of the thermal conductivity in the literature employing core-shell potentials.34 Sometimes, the core-shell part of the potential is stripped oﬀ, and the potential is employed as a rigid-ion potential.2 From the perspective of lattice dynamics, the shell aﬀects the dynamical matrix of the cores, which are the only true degrees of freedom in the system. This contribution can be computed exactly.44 Third-order force constants are aﬀected in a non-trivial way: a full analytical treatment of the third order constants is algebraically cumbersome and one generally relies on suitable approximations.44 A complete description of the procedure is beyond the scope of this paper; hence, we merely outline the necessary steps. One starts by deﬁning the Hamiltonian of the system in terms of the deviations from the equilibrium positions of both cores and shells up to some order. It is convenient to use core coordinates and the diﬀerences between shell and core coordinates as independent variables. Equations of motion are then written for the cores and the shells. Equations for the cores are the standard Newtonian diﬀerential equation of motion. Equations for the shells, on the other hand, due to their zero mass, are algebraic. One uses this second set of equations to eliminate the shell coordinates from the ﬁrst, renormalizing

3529

the force constants in the process. If one keep the terms in the Hamiltonian up to second order (harmonic approximation), then the shell equations are linear in terms of coordinates, and solution to the set of linear equations can be easily found. However, if cubic terms are also retained, then the shell equations are quadratic in terms of coordinates, and exact solution is diﬃcult. At this point, we recall that the core-shell displacement is generally small. One therefore simply neglects all terms higher than the ﬁrst order of this diﬀerence in both sets of equations. Shell equations are linear again, and shell coordinates can be eliminated. We implemented the correction to the third order force constants resulting from the core-shell model, as described in,44 and tested it on the only available case study in the literature of thermal conductivity in KCl.34 Details of the interatomic potentials for KCl are presented in Table I. Figure 9 shows the comparison between MD data obtained in34 and our LD calculations. The MD results were obtained by the Green–Kubo technique using the Ewald summation for the electrostatic interactions; hence, we also use the Ewald summation here. As one can observe, the thermal conductivity for the rigid-ion potential as calculated by BTE is signiﬁcantly higher the value produced by the MD method, especially at lower temperature. We believe that there are two reasons for that. First, the size of the system considered in34 was limited to the 64 atoms, which might not be suﬃciently large, especially at the low temperature, where thermal conductivity is larger, and therefore the phonon mean three path is also larger. Another investigation of the thermal conductivity in KCl45 by Green–Kubo MD compared two diﬀerent system sizes: 64 and 216 atoms. At 1000 K, the thermal conductivity was quite insensitive to the system size (0.88 and 0.78 W/mK), while at 500 K, there was a signiﬁcant increase (2.78 vs. 3.52 W/mK). Secondly, the crystal structure of KCl is rocksalt, the same as MgO, which, we found, also showed overestimation of the thermal conductivity by BTE. It therefore might be a structural feature of the rocksalt compounds that make higher order anharmonic contributions important. Interestingly, agreement between MD and LD is signiﬁcantly better for the core-shell potential. This might be a result of the fact that core-shell potential introduces many-body interactions, increases the anharmonicity in the system; consequently, the third-order force constants might become a dominant contribution to the thermal conductivity over higher order terms.

VI.

Conclusions

In this paper, we applied the BTE to compute the thermal conductivity of model ionic crystals. The RTA was found to

Fig. 9. Comparison of the core-shell model and rigid ion results in KCl. MD data are from reference [34].

3530

Journal of the American Ceramic Society—Chernatynskiy et al.

be a very good approximation in ionic solids, only slightly underpredicting the thermal conductivity compared with the full solution of the BTE. We also found that the thermal conductivity computed within quasiharmonic approximations is in exceptionally good agreement with MD simulations over a wide temperature range, even signiﬁcantly above the Debye temperature. At temperatures close to the Debye temperature or below it, ionic compounds were found to exhibit more complex behavior than covalent compounds in relationship with the importance of quantum eﬀects. In covalent compounds, lowfrequency phonon modes typically dominate thermal transport and classical statistics underestimate thermal conductivity when compared with quantum statistics. By contrast, in ionic compounds, optical modes make a larger contribution to the thermal conductivity resulting in either a fortuitous agreements between classical and quantum treatment or even an overestimate by the classical statistics compared with quantum. At temperatures signiﬁcantly lower than the Debye temperature in which quantum eﬀects dominate, the situation is similar to that of covalent compounds: namely, the classical description produces inﬁnite thermal conductivity at absolute zero, due to the constant heat capacitance, while the quantum description produces the correct limit of zero thermal conductivity. However, at the temperatures where the quantum description produces a peak in the thermal conductivity, its value is signiﬁcantly higher in the quantum description than in classical description. The direct summation method for the electrostatic interactions was found to underpredict or overpredict thermal conductivity, when compared with exact Ewald summation depending on the speciﬁc parameters used. Deﬁciencies in the direct summation method in reproducing the phonon structure correctly were found to be responsible for the discrepancy. Finally, third order force constants were found to be quite short-ranged, in spite of the long-range nature of the Coulombic interactions. The core-shell model was also successfully implemented in the framework of calculating thermal conductivity with the BTE. In summary, approximate methods for Columbic interactions and a classical treatment of the thermodynamic phonons properties result in erroneous values of the thermal conductivity with the direction of error that is diﬃcult to determine a priori for ionic compounds. These approximations therefore should be taken into account rigorously in calculations. The RTA always only slightly underpredicts the thermal conductivity, and therefore is often suﬃciently accurate. Rigorous lattice dynamics (if we use only the zero-temperature volume) always overpredicts the thermal conductivity; however, the quasiharmonic approximation produces exceptionally good results over a wide temperature range.

Acknowledgments This work was co-authored by subcontractors (AC, SRP) of the U.S. Government under DOE Contract No. DE-AC07-05ID14517, under the Energy Frontier Research Center (Oﬃce of Science, Oﬃce of Basic Energy Science, FWP 1356). Accordingly, the U.S. Government retains and the publisher (by accepting the article for publication) acknowledges that the U.S. Government retains a nonexclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for U.S. Government purposes.

References 1

D. R. Clarke and C. G. Levi, “Materials Design for the Next Generation Thermal Barrier Coatings,” Annu. Rev. Mater. Sci., 33, 383–417 (2003). 2 T. Watanabe, S. G. Srivilliputhur, P. K. Schelling, J. S. Tulenko, S. B. Sinnott, and S. R. Phillpot, “Thermal Transport in oﬀ-Stoichiometric Uranium Dioxide by Atomic Level Simulation,” J. Am. Ceram. Soc., 92 [4] 850–6 (2009). 3 K. Koumoto, Y. Wang, R. Zhang, A. Kosuga, and R. Funahashi, “Oxide Thermoelectric Materials: A Nanostructuring Approach,” Annu. Rev. Mater. Res., 40 [1] 363–94 (2010). 4 S. Stackhouse and L. Stixrude, “Theoretical Methods for Calculating the Lattice Thermal Conductivity of Minerals,” Rev. Mineral. Geochem., 71 [1] 253–69 (2010).

Vol. 94, No. 10

5 P. K. Schelling, S. R. Phillpot, and P. Keblinski, “Comparison of AtomicLevel Simulation Methods for Computing Thermal Conductivity,” Phys. Rev. B, 65 [14] 144306 (2002). 6 D. P. Sellan, E. S. Landry, J. E. Turney, A. J. H. McGaughey, and C. H. Amon, “Size Eﬀects in Molecular Dynamics Thermal Conductivity Predictions,” Phys. Rev. B, 81 [21] 214305 (2010). 7 D. A. Broido, A. Ward, and N. Mingo, “Lattice Thermal Conductivity of Silicon From Empirical Interatomic Potentials,” Phys. Rev. B, 72 [1] 014308 (2005). 8 A. Chernatynskiy and S. R. Phillpot, “Evaluation of Computational Techniques for Solving the Boltzmann Transport Equation for Lattice Thermal Conductivity Calculations,” Phys. Rev. B, 82 [13] 134301 (2010). 9 M. Omini and A. Sparavigna, “An Iterative Approach to the Phonon Boltzmann Equation in the Theory of Thermal Conductivity,” Physica B, 212 [2] 101–12 (1995). 10 S. Pettersson, “Solving the Phonon Boltzmann Equation With the Variational Method,” Phys. Rev. B, 43 [11] 9238–46 (1991). 11 J. E. Turney, E. S. Landry, A. J. H. McGaughey, and C. H. Amon, “Predicting Phonon Properties and Thermal Conductivity From Anharmonic Lattice Dynamics Calculations and Molecular Dynamics Simulations,” Phys. Rev. B, 79 [6] 064301 (2009). 12 N. Ohtori, M. Salanne, and P. A. Madden, “Calculations of the Thermal Conductivities of Ionic Materials by Simulation with Polarizable Interaction Potentials,” J. Chem. Phys., 130 [10] 104507 (2009). 13 G. Horton and E. Cowley, “Lattice Dynamics at High Temperature”; pp. 50–80 in Physics of Phonons, vol. 285 of Lecture Notes in Physics, Edited by T. Paszkiewicz, Springer, Berlin, Heidelberg, 1987, doi: 10.1007/3-54018244-6_33. 14 A. Ward, D. A. Broido, D. A. Stewart, and G. Deinzer, “Ab Initio Theory of the Lattice Thermal Conductivity in Diamond,” Phys. Rev. B, 80 [12] 125203 (2009). 15 A. Ward and D. A. Broido, “Intrinsic Phonon Relaxation Times from First-Principles Studies of the Thermal Conductivities of Si and Ge,” Phys. Rev. B, 81 [8] 085205 (2010). 16 L. Lindsay and D. A. Broido, “Optimized Tersoﬀ and Brenner Empirical Potential Parameters for Lattice Dynamics and Phonon Thermal Transport in Carbon Nanotubes and Graphene,” Phys. Rev. B, 81 [20] 205441 (2010). 17 S. Pettersson, “Calculation of the Thermal Conductivity of Alkali Halide Crystals,” J. Phys. C Solid State, 20 [8] 1047 (1987). 18 X. Tang and J. Dong, “Pressure Dependence of Harmonic and Anharmonic Lattice Dynamics in MgO: A First-Principles Calculation and Implications for Lattice Thermal Conductivity,” Phys. Earth Planet. In., 174 [1–4] 33– 8 (2009). 19 X. Tang and J. Dong, “Lattice Thermal Conductivity of MgO at Conditions of Earth’s Interior,” Proc. Natl. Acad. Sci. USA, 107 [10] 4539–43 (2010). 20 A. Y. Toukmaji and J. A. Board, “Ewald Summation Techniques in Perspective: A Survey,” Comput. Phys. Commun., 95 [2–3] 73–92 (1996). 21 D. Wolf, P. Keblinski, S. Phillpot, and J. Eggebrecht, “Exact Method for the Simulation of Coulombic Systems by Spherically Truncated, Pairwise r(-1) Summation,” J. Chem. Phys., 110 [17] 8254–82 (1999). 22 A. Angoshtari and A. Yavari, “Eﬀect of External Normal and Parallel Electric Fields on 180Aˆ° Ferroelectric Domain Walls in PbTiO 3,” J. Phys. Condens. Matter., 23 [3] 035901 (2011). 23 F. Zhang, A. M. Walker, K. Wright, and J. D. Gale, “Defects and Dislocations in MgO: Atomic Scale Models of Impurity Segregation and Fast Pipe Diﬀusion,” J. Mater. Chem., 20, 10445–51 (2010). 24 F. N. Mendoza, J. Lopez-Lemus, G. A. Chapela, and J. Alejandre, “The Wolf Method Applied to the Liquid-Vapor Interface of Water,” J. Chem. Phys., 129 [2] 024706 (2008). 25 A. J. H. McGaughey and M. Kaviany, “Thermal Conductivity Decomposition and Analysis Using Molecular Dynamics Simulations: Part II. Complex Silica Structures,” Int. J. Heat Mass Tran., 47 [8–9] 1799–816 (2004). 26 M. Fe`vre, A. Finel, and R. Caudron, “Local Order and Thermal Conductivity in Yttria-Stabilized Zirconia. I. Microstructural Investigations Using Neutron Diﬀuse Scattering and Atomic-Scale Simulations,” Phys. Rev. B, 72 [10] 104117 (2005). 27 A. J. Kulkarni and M. Zhou, “Size-Dependent Thermal Conductivity of Zinc Oxide Nanobelts,” Appl. Phys. Lett., 88 [14] 141921 (2006). 28 A. J. Kulkarni and M. Zhou, “Tunable Thermal Response of ZnO Nanowires,” Nanotechnology, 18 [43] 435706 (2007). 29 B.-L. Huang and M. Kaviany, “Ab Initio and Molecular Dynamics Predictions for Electron and Phonon Transport in Bismuth Telluride,” Phys. Rev. B, 77 [12] 125209 (2008). 30 N. de Koker, “Thermal Conductivity of MgO Periclase at High Pressure: Implications for the D” Region,” Earth Planet. Sci. Lett., 292 [3–4] 392–8 (2010). 31 P. Shukla, T. Watanabe, J. C. Nino, J. S. Tulenko, and S. R. Phillpot, “Thermal Transport Properties of MgO and Nd2Zr2O7 Pyrochlore by Molecular Dynamics Simulation,” J. Nucl. Mater., 380 [1–3] 1–7 (2008). 32 K. Udayakumar and A. Cormack, “Non-Stoichiometry in Alkaline-Earth Excess Alkaline-Earth Titanates,” J. Phys. Chem. Solids, 50 [1] 55–60 (1989). 33 M. Abramowski, R. Grimes, and S. Owens, “Morphology of UO2,” J. Nucl. Mater., 275 [1] 12–8 (1999). 34 G. Paolini, P. Lindan, and J. Harding, “The Thermal Conductivity of Defective Crystals,” J. Chem. Phys., 106 [9] 3681–7 (1997). 35 C. Ronchi, J. P. Ottaviani, C. Degueldre, and R. Calabrese, “Thermophysical Properties of Inert Matrix Fuels for Actinide Transmutation,” J. Nucl. Mater., 320 [1–2] 54–65 (2003).

October 2011

Elucidation of Thermal-Transport Properties of Ionic Materials

36 J. Liu, C. Wang, W. Su, H. Wang, J. Li, J. Zhang, and L. Mei, “Thermoelectric Properties of Sr1-XNdxTiO3 Ceramics,” J. Alloys Compd., 492 [1–2] L54–6 (2010). 37 C. Ronchi, M. Sheindlin, M. Musella, and G. J. Hyland, “Thermal Conductivity of Uranium Dioxide up to 2900 K From Simultaneous Measurement of the Heat Capacity and Thermal Diﬀusivity,” J. Appl. Phys., 85 [2] 776–89 (1999). 38 J. E. Turney, A. J. H. McGaughey, and C. H. Amon, “Assessing the Applicability of Quantum Corrections to Classical Thermal Conductivity Predictions,” Phys. Rev. B, 79 [22] 224305 (2009). 39 M. Ahrens, R. Merkle, B. Rahmati, and J. Maier, “Eﬀective Masses of Electrons in n-Type SrTiO3 Determined from Low-Temperature Speciﬁc Heat Capacities,” Physica B, 393 [1–2] 239–48 (2007).

3531

40 P. G. Klemens, “The Scattering of low-Frequency Lattice Waves by Static Imperfections,” Proc. Phys. Soc. A, 68 [12] 1113 (1955). 41 G. Venkataraman, F. A. Feldkamp, and V. C. Sahnim, Dynamics of the Perfect Crystals. The MIT Press, Cambridge, 1975. 42 B. G. Dick and A. W. Overhauser, “Theory of the Dielectric Constants of Alkali Halide Crystals,” Phys. Rev., 112 [1] 90 (1958). 43 P. J. Mitchell and D. Fincham, “Shell Model Simulations by Adiabatic Dynamics,” J. Phys.: Condens. Matter., 5 [8] 1031 (1993). 44 J. Oitmaa, “Lattice Dynamics of Simple Harmonic and Anharmonic Shell Models,” Aust. J. Phys., 20 [5] 495 (1967). 45 P. Sindzingre and M. Gillan, “A Computer-Simulation Study of Transport-Coeﬃcients in Alkali-Halides,” J. Phys.-Condens. Mat., 2 [33] 7033–42 (1990). h