# Prediction of Spectral Phonon Mean Free Path and Thermal Conductivity with Applications to Thermoelectrics and Thermal Management: A Review

Authors: Tianli Feng; Xiulin Ruan
Journal: Journal of Nanomaterials
 (2014)

---

## Abstract

We give a review of the theoretical approaches for predicting spectral phonon mean free path and
thermal conductivity of solids. The methods can be summarized into two categories: anharmonic lattice dynamics calculation and molecular dynamics simulation. In the anharmonic lattice dynamics calculation, the anharmonic force constants are used first to calculate the phonon scattering rates, and then the Boltzmann transport equations are solved using either standard single mode relaxation time approximation or the Iterative Scheme method for the thermal conductivity. The MD method involves the time domain or frequency domain normal mode analysis. We present the theoretical frameworks of the methods for the prediction of phonon dispersion, spectral phonon relaxation time, and thermal conductivity of pure bulk materials, layer and tube structures, nanowires, defective materials, and superlattices. Several examples of their applications in thermal management and thermoelectric materials are given. The strength and limitations of these methods are compared in several different aspects. For more efficient and accurate predictions, the improvements of those methods are still needed.


---

## Body

## 1. Introduction

In recent years, increasing attention has been focused on seeking novel structures and materials with desired thermal properties, especially thermal conductivity. High thermal conductivity can help remove heat rapidly and reduce device temperatures so as to improve performance of nanoelectronics and optoelectronics, while low thermal conductivity is desired in thermoelectrics for improving the figures of meritZT [1] of the material: ZT=S2σT/k, where S, σ, T, and k are Seebeck coefficient, electronic conductivity, temperature, and thermal conductivity, respectively. The thermal conductivity k is a summation of the lattice contribution kl and electron contribution ke. Since, in most thermoelectric materials, the phonon mean free path is much longer than that of electrons, one major strategy to enhance ZT is to reduce kl without much affecting ke. This is made possible by the rapid development of nanofabrication techniques.Gaining a deeper physical insight into the spectral phonon properties, for example, the spectral phonon relaxation time and mean free path, is necessary to correctly explain experimental results and accurately predict and guide the further designs and applications. Analytical models have been used by Balandin and Wang to estimate frequency-dependent phonon group velocity and various phonon scattering rates including phonon-phonon, phonon-impurity, and phonon-boundary scattering processes. They used this approach to observe the strong modification of acoustic phonon group velocity and enhanced phonon scattering rate due to boundary scattering in semiconductor quantum wells, so as to successfully explain their significantly reduced lattice thermal conductivity [2]. This effect of phonon confinement was then extended to nanowires and quantum dot superlattices [3–5]. Analytical models of spectral phonon properties are advantageous in their clear physical insights, but they usually contain empirical fitting parameters, and this limitation has motivated the development of numerical methods based on first principles and molecular dynamics that can predict these spectral properties from their atomic structure, without fitting parameters and with greater accuracy. This review will be focused on these predictive simulation methods.The methods of predicting spectral phonon relaxation times and mean free paths become increasingly important for predicting the thermal properties of numerous novel materials. For instance, superlattice structure is found to be an effective way to suppress the thermal conductivity because of the interface mass mismatch scattering [6–8], but the phenomenon in which the short-period superlattice can have even higher thermal conductivity still needs the deep insight of phonon relaxation time. Doping and alloying are widely used to explore novel high performance materials [9], and the natural materials are rarely pure; thus the study of impurity scattering contributing to phonon relaxation time is important. Layer and tube structured materials, for example, graphene [10–13] and carbon nanotube (CNT) [14–17], are proved to have unusual phonon transport features such as high thermal conductivity [18–20], which still need to be further understood. Comprehensive reviews of their thermal transport can be found in [21–23]. Nanowires are most commonly studied and used in both theoretical and experimental research [24–32], and the accurate prediction of thermal conductivity needs the knowledge of spectral phonon scattering by boundaries.Many methods have been proposed and applied to predict spectral phonon relaxation time in the last half century. At the earliest, Klemens and other researchers obtained the frequency-dependent phonon relaxation time mostly by long-wave approximation (LWA) and Deybe model: Klemens obtained the phonon relaxation times by Umklapp (U) three-phonon scattering [33, 34] and defect scattering [35], Herring studied normal (N) three-phonon scattering [36], Holland extended the results to dispersive transverse mode range [37], and Casimir studied boundary scattering [38–40]. A more accurate method, the third-order anharmonic lattice dynamics (ALD) calculation which can predict the intrinsic spectral phonon relaxation times without LWA, was presented by Maradudin and the coworkers [41, 42]. ALD methods were then applied to silicon and germanium byab initio approach first by Debernardi et al. [43] and Deinzer et al. [44]. Beyond the standard ALD calculation, Omini and Sparavigna [45, 46] proposed an Iterative Scheme which gives exact solutions to the linearized Boltzmann transport equation (BTE). The Iterative Scheme has been successfully applied to many structures in the recent ten years by Broido, Lindsay, Ward, and so forth [47–65]. Other than the lattice dynamics calculation, a time domain normal mode analysis (NMA) method based on molecular dynamics (MD) simulation was proposed by Ladd et al. [66] and extended by McGaughey and Kaviany [67]. Another version of normal mode analysis is implemented in frequency domain, so called spectral energy density (SED) analysis. The normal mode analysis was early implemented by Wang et al. [68–75] to obtain the relaxation times of a few phonon modes and then extended by de Koker [76] and Thomas et al. [77, 78] to calculate lattice thermal conductivity.In this work, we present a review of the methods of predicting spectral phonon properties, discuss the applications to each method, and compare them in different aspects. Section2 gives an overview of thermal conductivity and the frequency-dependent relaxation time predicted from early long-wave approximation (LWA) and Debye model. Section 3 presents the ALD calculation which is divided into three subsections: Section 3.1 covers the standard single mode relaxation time approximation (SMRTA), Section 3.2 gives the Iterative Scheme, and Section 3.3 reviews examples of the applications to pure bulk, layer and tube structures, nanowires, defective materials, and superlattice. In Section 4, we introduce the time domain NMA and frequency NMA methods based on MD simulations and their applications. The summary is presented in Section 5. The appendix provides some derivations of ALD methods.


## 2. Theory Overview

Spectral phonon mean free path (MFP), determined by phonon scattering rate, dominates the behavior of thermal properties, especially the thermal conductivityk. Based on BTE under the relaxation time approximation (RTA), thermal conductivity is determined by the spectral phonon relaxation time τλ, phonon group velocity vλ=∇kω, and phonon specific heat cλ [34]:
(1)kz=1V∑λ‍(vλ·z^)2cλτλ,
where z^ denotes the transport direction, λ is the shorthand of phonon mode (k,ν) with k representing the phonon wave vector and ν labeling phonon dispersion branch, V is the volume of the domain, and the summation is done over the resolvable phonon modes in the domain. The specific heat per mode is cλ=ℏωλ∂nλ0/∂T=kBx2ex/(ex-1)2, where nλ0 is phonon occupation number of the Bose-Einstein distribution nλ0=(ex-1)-1 and x is the shorthand of ℏω/kBT. Equation (1) can also be expressed in terms of phonon mean free path Λ→λ=vλτλ. The continuous form of (1) is, with the help of ∑k‍=(V/(2π)3)∫dk,
(2)kz=1(2π)3∑ν‍∫‍(vλ·z^)2cλτλdk.If isotropic heat transport is assumed, the integration of|vλ,z|2 in (2) gives vλ2/3, and ∫dk gives ∫4πk2dk, we get the commonly used formula
(3)kz=4π31(2π)3∑ν‍∫‍cλvλΛλk2dk.The early theoretical predictions of phonon relaxation times for different scattering processes are briefly summarized in Table1 [37, 79]. T is temperature; subscripts N, U, T, and L indicate the Umklapp scattering, normal scattering, transverse wave, and longitudinal wave, respectively; A, B’s, and C’s are constants; θ is Debye temperature; α is numerical constant in [34]; Low T means T≪θ, and high T means T≫θ. ω1 is the transverse mode frequency at which the group velocity starts to decrease, and ω2 is the maximum transverse frequency. The intrinsic three-phonon scattering rates are derived mostly in LWA or linear dispersion approximation.Table 1
Analytical models of inverse relaxation time for different scattering processes.
Scattering process
Inverse relaxation time
Intrinsic Three-Phonon
N process
Herringa
τ
L
,
N
-
1
=
B
L
ω
2
T
3, low T
τ
T
,
N
-
1
=
B
T
ω
T
4, low T
τ
L
,
N
-
1
=
B
L
′
ω
2
T, high T
τ
T
,
N
-
1
=
B
T
′
ω
T, high T
Callawayb
τ
N
-
1
=
B
N
ω
2
T
3
U process
Klemensc
τ
U
-
1
=
B
U
ω
2
T
3
exp
⁡
(
-
θ
α
T
), low T
Klemensd
τ
U
-
1
=
B
U
ω
T
3
exp
⁡
(
-
θ
α
T
), low T
τ
U
-
1
=
B
′
ω
2
T, high T
Callawayb
τ
U
-
1
=
B
U
ω
2
T
3
Hollande
τ
T
,
U
-
1
=
{
B
T
U
ω
2
sinh
⁡
(
x
)
,
ω
1
≤
ω
≤
ω
2
0
,
ω
<
ω
1
Asen-Palmer et al.f
τ
T
,
U
-
1
=
B
T
U
′
ω
2
T
exp
⁡
(
C
T
T
)
τ
L
,
U
-
1
=
B
L
U
′
ω
2
T
exp
⁡
(
C
L
T
)
Boundary
g
τ
b
-
1
=
v
λ
l
F  l: diameter, l=2l1l2/π;      F: surface roughmess
Impurity
h
τ
im
-
1
=
π
2
g
ω
2
D
(
ω
)
≈
V
0
g
4
π
〈
v
3
〉
ω
4
~
A
ω
4
References:a[36], b[80], c[33], d[34], e[37], f[81], g[38–40], and h[35].Boundary scatteringτb-1 exists anywhere, since every sample has a finite size. F captures the boundary scattering characteristic of the sample with F=1 representing completely diffusive and F→∞ meaning specular. l=2Sc/π is a measure of the size perpendicular to the transport direction, with Sc being the area of cross section. vλ is often replaced by the average phonon speed of the three acoustic branches vave for simplicity [37]:
(4)vave=[13∑ν3‍1vν]-1=[13(2vT+1vL)]-1.The last equation in Table1 takes into account the impurity scattering rate, where
(5)g=∑i‍fi(1-mim¯)2
is a measure of mass disorder, D(ω) is phonon density of states normalized to unity, fi is the concentration of the impurity species i, and mi and m¯ are the mass of i and average mass for the given composition, respectively. The exact expression for 〈v3〉 is found in [82], while, in long wave approximation, 〈v3〉 approximates the cube of acoustic phonon speed of the material. This equation was derived by Klemens for isotope scattering with only mass disorder. For crystal defects other than isotope doping, such as vacancy, interstitial, and antisite defects, the impurity scattering comes from not only the mass disorder but also the interatomic force change and link break. Klemens took into account such effect by adding a modification to g:
(6)g~=g+2∑i‍fi(Δϕiϕ-6.4γΔrir)2,
where ϕi/ϕ, and Δri/r describe the average relative variations of the local force constants and atomic displacements [35, 83–85], respectively. Some consider the dislocations by adding a scattering term τD-1~ω to the total phonon scattering rate [84], predicted from single dislocation assumption by [34, 35, 86, 87]. Although τim-1=Aω4 is derived for low frequency phonons, many works use it to predict thermal conductivity or explain data from experiments for alloys and crystals with impurities [24–26, 37, 81, 84, 85, 88–91]. In Section 3.3.4, we will give more precise expressions for isotope scattering.For the system that contains several scattering mechanisms, the Matthiessen rule is often used to evaluate the total scattering rate,(7)τ-1=∑i‍τi-1.
In most cases the Matthiessen rule gives reasonable results, although it is found to be not accurate in some cases recently [58, 92, 93].These frequency dependent relaxation time expressions in Table1 have been used in many works for thermal conductivity prediction and analysis, and the choice of those expressions looks quite arbitrary. For instance, in the choice of intrinsic phonon relaxation time in the thermal conductivity analysis of silicon, Glassbrenner and Slack [94] used τ-1~ω2T, while Asen-Palmer et al. [81] and Mingo et al. [24, 25] used τ-1~ω2Texp⁡(C/T) for all phonon modes; Martin et al. [26] used τL-1~ω2T3 for longitudinal mode, while Holland added τT-1~ω2/sinh⁡(x) to dispersive transverse range. The thermal conductivity results predicted by these expressions can be reasonable due to the adjustable fitting parameters. Therefore, it becomes important to accurately predict spectral phonon relaxation time without any fitting parameter, which allows us to understand thermal transport and examine (a) the validity of low-frequency approximation or the Debye model, (b) the importance of optical branch to thermal transport, (c) the contributions of phonons with different mean free path or different wavelength to thermal conductivity, (d) the relative importance of different scattering mechanisms in a given material, and so forth.


## 3. Anharmonic Lattice Dynamics Methods

In perturbation theory, the steady-state phonon BTE [34, 79, 95] describes the balance of phonon population between diffusive drift and scattering as
(8)vλ·∇nλ=∂nλ∂t|s,
where nλ=nλ0+nλ′ is the total phonon occupation number with nλ′ representing the deviation from the equilibrium phonon distribution nλ0. With ∇nλ=(∂nλ/∂T)∇T and assuming that nλ′ is independent of temperature: (∂nλ/∂T)≃(∂nλ0/∂T), we have
(9)vλ·∇T∂nλ0∂T=∂nλ′∂t|s.The RTA assumes that deviation of single phonon mode population decays exponentially with time:(10)nλ′~exp⁡(-tτλ),
where τλ is the relaxation time. Therefore, the collision term in BTE (9) becomes
(11)∂nλ′∂t|s≃-nλ′τλ.
Generally, the value of τλ is considered as the average time between collisions of the phonon mode λ with other modes, whereby τλ=1/Γλ, where Γλ denotes the scatting rate.Considering only three-phonon scattering, (9) becomes [95]
(12)vλ·∇T∂nλ0∂T=-∑λ′λ′′‍{[nλ(1+nλ′)(1+nλ′′)-(1+nλ)nλ′nλ′′]𝔏-12[nλ(1+nλ′)(1+nλ′′)-(1+nλ)nλ′nλ′′]𝔏-[nλnλ′(1+nλ′′)-(1+nλ)(1+nλ′)nλ′′]𝔏+hh+12[nλ(1+nλ′)(1+nλ′′)-(1+nλ)nλ′nλ′′]𝔏-},
where the summation is done over all the phonon modes λ′ and λ′′ that obey the energy conservation ωλ±ωλ′=ωλ′′ and quasimomentum conservation k±k′=k′′+G with G=0 for N processes and G≠0 for U processes, where G is a reciprocal-lattice vector. 𝔏± is the probability of λ±λ′→λ′′ scattering occurrence, determined via Fermi’s golden rule
(13)𝔏±=ℏπ4N0|V±(3)|2δ(ωλ±ωλ′-ωλ′′)ωλωλ′ωλ′′,(14)V±(3)=∑bl′b′l′′b′′‍∑αβγ‍Φ0b,l′b′,l′′b′′αβγeαbλeβb′±λ′eγb′′-λ′′m¯bm¯b′m¯b′′e±ik′·rl′e-ik′′·rl′′,
where b’s and l’s are the indexes of basis atoms and unit cells, respectively, α, β, and γ represent coordinate directions, mb is the mass of basis atom b, considering that some doping material m¯b is the average mass in the bth basis sites, eb,αλ is the α component of the bth part of the mode λ=(k,ν)’s eigenvector, and Φ is the third-order interatomic force constant (IFC). The factor “1/2” in (12) accounts for the double counting in the summation of λ′ and λ′′ for the “−” process. In (14), the factor eik·rl is often omitted, since it is a constant in the summation and thus contributes nothing to |V±(3)|2.


### 3.1. Standard Single Mode Relaxation Time Approximation

The Standard SMRTA assumes that the system is in its complete thermal equilibrium, except that one phonon modeλ has its occupation number nλ=nλ0+nλ′ differing a small amount from its equilibrium value nλ0. Therefore, on the right hand side of (12), replacing nλ by nλ0+nλ′, whilst nλ′ and nλ′′ by nλ′0 and nλ′′0, respectively, one can obtain the phonon relaxation time τλ0 of mode λ (for the derivation, see Appendix A.1):
(15)1τλ0=∑λ′λ′′+‍Γλλ′λ′′++12∑λ′λ′′-‍Γλλ′λ′′-+∑λ′‍Γλλ′ext,
where the first two terms on the right hand side are intrinsic three-phonon scattering rates (ωλ±ωλ′=ωλ′′):
(16)Γλλ′λ′′±=ℏπ4N0{nλ′0-nλ′′0nλ′0+nλ′′0+1}×|V±(3)|2δ(ωλ±ωλ′-ωλ′′)ωλωλ′ωλ′′.
The last term Γλλ′ext represents the extrinsic scattering such as boundary scattering and impurity scattering.


### 3.2.  Iterative Scheme: Exact Solution to Linearized BTE

Different from the Standard SMRTA, the other method to solve the phonon BTE allows all the modes to be in their thermal nonequilibrium states at the same time. By replacing the occupation numbersnλ, nλ′, and nλ′′ by nλ0+nλ′, nλ′0+nλ′′, and nλ′′0+nλ′′′, respectively, on the right hand side of (12), the relaxation time τλ of mode λ is obtained (for the derivation see Appendix A.2)
(17)τλ=τλ0(1+Δλ),(18)Δλ=∑λ′λ′′(+)‍Γλλ′λ′′+(ξλλ′′τλ′′-ξλλ′τλ′)+∑λ′λ′′(-)‍12Γλλ′λ′′-(ξλλ′′τλ′′+ξλλ′τλ′)+∑λ′‍Γλλ′extξλλ′τλ′,
where ξλλ′=vλ′,zωλ′/vλ,zωλ, and vz is phonon group velocity component along the transport direction.Equation (17) is solved iteratively because both the left and the right hand sides contain the unknown variable τλ, and thus the method is called Iterative Scheme. This scheme is also based on RTA; thus (10) and (11) are still valid (one can reach this by substituting (A.3), (A.4), (A.11), (A.12), (A.13), and (A.14) into (9)). The last summation in (18) is done over λ′ with λ′≠λ.


### 3.3. Discussions and Applications

ALD methods can be divided into classical method andab initiomethod, differing in how to calculate the harmonic and anharmonic IFCs, which are the only inputs to these methods. The classical approach relies on empirical interatomic potential whose nth order derivatives are taken as the nth order IFCs:
(19)Φl1b1,…,lnbnα1⋯αn=∂nΦ∂u1(l1b1),…,∂un(lnbn).
In contrast, theab initio approach is a first principle calculation in the framework of density functional perturbation theory (DFPT) [43, 96, 97] using norm-conserving pseudopotentials in the local density approximation (LDA) without introducing any adjustable parameters. The formulism of the IFCs using first principle method can be found in [44] and realized by, for example, the QUANTUM ESPRESSO package [98]. Compared to the classical method, this method can deal with new materials whose empirical interatomic potentials are unknown. Further, this method can be more accurate since the empirical interatomic potentials cannot always represent the exact nature of interatomic force.In (16), the delta function δ(ω±ω′-ω′′) is typically approximated by δ(x)=lim⁡ε→0+(1/π)(ε/(x2+ε2)). To accurately evaluate (16), the choice of ε value is critical: it must be small but larger than the smallest increment in discrete x, which results from the use of finite grid of k points in Brillouin zone. The general practice is as follows: pick the densest grid possible and start with a sufficiently small guess, and increase it gradually until the final results reach convergence.To calculate the relaxation time, one can use Standard SMRTA scheme [43, 44, 99–111] or Iterative Scheme [47–65, 112], and, in each of them, one can choose empirical interatomic potential approach [47–49, 52–56, 58–60, 99–102] orab initio-derived IFC IFC [43, 44, 50, 51, 57, 61–65, 103–112]. The methods can be used on pure bulk, nanowires, doped bulk, doped nanowires, alloys, and so forth.One way to predict thermal conductivityk without working out all the phonon modes relaxation times is the Monte Carlo integration technique [101, 113]. The protocol of this technique is as follows: (1) randomly sample some phonon modes λ, (2) for each of these modes, randomly choose two other modes λ′ and λ′′ that interact with λ to calculate the relaxation time, and (3) select as many points as necessary to ensure that the statistical error is small enough in both cases. Monte Carlo technique only works for the Standard SMRTA scheme, since the Iterative Scheme requires the relaxation times of all the phonon modes to do iteration. Monte Carlo technique reduces the computational cost but lowers the accuracy.In addition to intrinsic phonon scatteringΓ±, extrinsic scattering 1/τλext plays an important role in nanostructures,
(20)1τλext=∑λ′‍Γλλ′ext,
such as boundary scattering 1/τλbs and impurity scattering 1/τλimp.


#### 3.3.1. Intrinsic Phonon Scattering: Bulk Materials

Without any fitting parameters, Standard SMRTA withab initio approach can accurately predict spectral phonon relaxation times and thermal conductivities. Ward and Broido [110] checked the validity of some old approximations introduced in Section 2: (1) long-wave approximation for three-phonon scattering and (2) ignoring optical phonons, using silicon and germanium as examples. First, the values of matrix element |V±(3)|2, which govern the scattering strength Γ, fromab initio calculation for acoustic phonons are compared to those given by LWA. The percentage error of |V±(3)|2 is shown in Figure 1. We note that the LWA only works for the very low frequency ωeff<0.8 THz, while, for most part 0.8<ωeff<12 THz, the LWA gives large discrepancy, where ωeff≡(ωλωλ′ωλ′′)1/3 is the geometric average of the three-phonon frequencies. Second, the relaxation times of optical modes are found to only contribute less than 10% to the total thermal conductivity of silicon. However, ignoring optical modes is erroneous since the optical phonons are essential to provide channels for acoustic phonon scattering. The explicit calculation of millions of three-phonon scattering shows that optical phonons are involved in 50–60% of the total acoustic phonon-scattering processes in Si and Ge. Last, beyond the ω,T dependencies of τ listed in Table 1, which rely on many approximations, the ALD calculation can give more precise τω,T dependence. This is illustrated in the inset of Figure 1, the relaxation times of the LA phonons in Si for T=300K. By decomposing the total scattering into U process and N process, we find the U process has a stronger frequency dependence τ(U)(ω)~ω-4 than N process τ(N)(ω)~ω-2. The results also show that normal scattering governs the total relaxation time τ(0)(ω) at low frequency, while Umklapp scattering dominates at high frequency. Such ~ω-4 relation is not expected in the analytical models in Table 1.Figure 1
Percent error (color online) in|V±(3)|2 from the LWA compared to first principle for silicon at 300 K. Insert shows the normal (blue dashed curve), Umklapp (green dotted curve), and total (red solid curve) relaxation times for the LA phonons calculated from Standard SMRTA byab initio approach. Adapted with permission from [110]. Copyrighted by the American Physical Society.One flaw of the Standard SMRTA is that it does not grasp the interplay between theN process and U process. The right hand side of (15) can be decomposed as Γλ(N)+Γλ(U) (only consider intrinsic phonon scattering), according to whether they are N or U scattering events. The Standard SMRTA scheme treats the N process and U process as two independent scattering events and use Matthiessen’s rule to account for the total relaxation time 1/τλ0=1/τλ(N)+1/τλ(U), where τλ(N,U) is defined as τλ(N,U)≡1/Γλ(N,U). However, it is well know that N process does not contribute to thermal resistance directly. Instead, it affects the U process (low-frequency N scattering produces high-frequency phonons which boosts U process), and then the U process produces thermal resistance. This error can be remedied in the Iterative Scheme by doing the iteration in (17). Therefore, the Standard SMRTA scheme only works for the system, where U process dominates so that the N scattering makes little difference to U process as well as to thermal resistance [51, 110].For Si and Ge at room temperature where theU process is strong, the thermal conductivity predicted by Standard SMRTA scheme is only 5–10% smaller than that by Iterative Scheme [110], the latter shows excellent agreement with experiment (see Figure 1 of [61]).In contrast, theU scattering in diamond is much weaker [114–116] due to the much smaller phase space [117]. As a result, the thermal conductivity given by these two methods can differ by 50% at room temperature [110]. As shown in Figure 2, this discrepancy increases with decreasing temperature since the Umklapp scattering is weakened when temperature decreases. The thermal conductivity of diamond predicted by Iterative Scheme withab initio approach agrees excellently with experiment as shown in [110]. It is also noted that the Standard SMRTA scheme always underpredicts the thermal conductivity because it treats N process as an independent channel for thermal resistance. On the other hand, if the relaxation time for U process only is used in the calculation, the thermal conductivity is always overpredicted. This again confirms that the N process has an indirect and partial contribution to the thermal resistance.Figure 2
The calculated intrinsic lattice thermal conductivity of diamond for the Standard SMRTA (dashed line) and the Iterative Scheme (solid line), both byab initio approach. Dotted line shows percent error of the Standard SMRTA result compared to the Iterative Scheme solution. Reprinted with permission from [51]. Copyrighted by the American Physical Society.One important application of ALD calculation is to predict and understand the thermal conductivity of thermoelectric materials and help to design higher thermoelectric performance structures. Based on first principle calculation, Shiga et al. [104] obtain the frequency-dependent relaxation times of pristine PbTe bulk at 300 K as shown in Figure 3. At low-frequency region, TA phonons have longer relaxation times than LA phonons with τ’s exhibiting ~ω-2 dependence. Separating the scattering rates into those of normal and Umklapp processes, they find the relations τNormal~ω-2 and τUmklapp~ω-3, which again indicate that the normal process dominates low-frequency region while the Umklapp dominates high-frequency part. By further studying the participation of each phonon mode to the total scattering rates, they find that the low thermal conductivity of PbTe is attributed to the strong scattering of LA phonons by TO phonons and the small group velocity of TA phonons. Figure 4 compares phonon relaxation times of PbTe and PbSe [106]. Although the anharmonicity of PbSe is normally expected to be larger due to the larger average Gruneisen parameter reported from experiments [121], in this work, it is found that, for TA mode, the relaxation times of PbSe are substantially longer than those of PbTe. Surprisingly, the optical phonons are found to contribute as much as 25% for PbSe and 22% for PbTe to the total thermal conductivity at the temperature range 300–700 K. Motivated by the question that phonons with what kind of MFP contribute the most to the total thermal conductivity, the cumulative k’s as functions of phonon MFP are calculated by ALD method with first principle approach as shown in Figure 5. Silicon is found to have phonon MFPs which span 6 orders of magnitude (0–106 nm), while the thermal transport in diamond is dominated by the phonon with narrow range of MFP (0.4–2 μm). It is found that the phonons with MFP below 4 μm for silicon, 1.6 μm for GaAs, 120 nm for ZrCoSb, 20nm for PbSe, and 10 nm for PbTe contribute 80% of total thermal conductivity. GaAs/AlAs superlattice is found to have similar phonon MFP with bulk GaAs. The curves of the alloy Mg2Si0.6Sn0.4 and its pure phases Mg2Si and Mg2Sn cross at the intermediate MFPs. These results provide great guidance for experimental works. For example, the PbTe-PbSe alloys with size of nanoparticle below 10 nm are synthesized and found to lead to as much as 60% reduction to the thermal conductivity which provides large space for improving ZT [122].(a) Spectral phonon relaxation times of pristine PbTe bulk at 300 K by Standard SMRTA scheme with first principle IFCs. Relaxation times of (b) normal and (c) Umklapp processes, respectively. The solid lines plot (a)τ=6×1012ω-2, (b) τ=8×1012ω-2, and (c) τ=8×1024ω-3. Reprinted with permission from [104]. Copyrighted by the American Physical Society.
(a)
(b)
(c)Spectral phonon relaxation times of PbSe bulk (squares) and PbTe bulk (crosses) at 300 K by Standard SMRTA scheme with IFCs from first principle calculation: (a) TA, (b) LA, (c) TO, and (d) LO. Reprinted with permission from [106]. Copyrighted by the American Physical Society.
(a)
(b)
(c)
(d)Figure 5
The accumulated thermal conductivity of different bulk materials as a function of phonon mean free path at room temperature calculated from ALD by first principle approach. (Diamond [64]; Si (iterative) [64]; Si (standard) [107]; GaAs [118]; FaAs/AlAs [119]; PbTe [106]; PbSe [106]; PbTe0.5Se0.5 [106]; Mg2Si [63]; Mg2Sn [63]; Mg2Si0.6Sn0.4 [63]; ZrCoSb [120]).


#### 3.3.2. Single and Few-Layer 2D Materials, Nanoribbons, and Nanotubes

For single- and multilayer 2D materials, the boundary scattering from the sides perpendicular to the transport direction is much weaker than for 3D systems [123], making the boundary scattering expression in Table 1 unsuitable. Instead, when studying single-/multilayer graphene (SLG/MLG) and graphite [54, 55], single-wall carbon nanotubes (SWCNTs) [52, 53], single-/multilayer boron nitride (SLBN/MLBN), and boron nitride nanotubes (BNNTs) [56, 58], Lindsay and Broido only consider the boundary scattering from the two ends in the transport direction and show that
(21)1τbs=2|vλ·z^|L
works well in accounting for the boundary scattering, with L being the length between boundaries in the transport direction z^. Such formula has been shown to give correct thermal conductivity values of nanotubes [124] and nanoribbons [125] in the ballistic limit (L→0) and diffusive limit (L→∞).Vibrations in 2D lattices are characterized by two types of phonons: those vibrating in the plane of layer (TA and LA) and those vibrating out of plane, so called flexural phonons (ZA and ZO). Lindsay et al. [54] find the selection rule for all orders in anharmonic phonon-phonon scattering in the 2D crystals:only even numbers (including zero) of flexural phonons can be involved, arising from the reflection symmetry perpendicular to the plane of layer. This selection rule has forbidden about 60% of both N and U three-phonon scattering phase space of ZA phonons for single layer graphene. They show that such suppressed scattering yields long relaxation time and mean free path for ZA phonons, leading to ZA phonons contributing most of the thermal conductivity of SLG, about 70% at room temperature (another cause being the large density of states and occupation number of ZA modes). However, this conclusion is still under debate since this approach does not include the fourth- and higher-order phonon scattering rates, which are not necessarily low since the reflection symmetry allows more 4-phonon processes than 3-phonon processes. Actually, the method of spectral energy analysis based on MD (discussed in Section 4) indicates that only 25%–30% of the total k is contributed by ZA mode at room temperature [126–128]. It should be noted that MD has its own drawback of not reproducing the Bose-Einstein distribution for graphene phonons at room temperature. Hence, the discrepancies between the two methods still need further study.The selection rule mentioned above does not hold for multilayer graphene, twisted graphene, graphite (because of the interlayer coupling), CNT (due to the curvature), graphene nanoribbon (GNR) (due to boundary scattering), substrate-supported graphene (due to scattering with the substrate), and defective graphene (due to defective scattering). Therefore, the thermal conductivity of these structures is typically lower than that of single layer graphene, and the contribution of each phonon mode changes [54, 129–132]. In Figure 6, single-layer graphene, GNR, and SWCNT are compared, where graphene has an infinite width and finite length L, SWCNT has a finite diameter d and length L, and GNR has a finite width W=πd with artificial periodic boundary condition applied. As expected, kRTA underpredicted thermal conductivity. SWCNT is found to have a lower thermal conductivity than graphene with a minimum value of 77% of kgraphene at a critical diameter d≈1.5nm. From this critical diameter, kSWCNT increases with increasing diameter and reaches 90% of kgraphene at d≈4 nm. On the other hand, if d goes small enough, phonon-phonon scattering decreases and the thermal conductivity increases. At this short limit of d, the system becomes more like a 1D chain which generally has much larger thermal conductivity than 2D and 3D systems. The increasing trend of kGNR with decreasing d comes from the reason that the decrease of the width πd pushes the optical modes to higher frequencies and thus the U scattering by optical phonons becomes weaker.Figure 6
Thermal conductivityk versus diameter d (color online) for single wall carbon nanotubes. Solid red circles, blue squares, and green triangles represent zigzag, armchair, and chiral predicted by Iterative Scheme; open red circles, blue squares, and green triangles are those from Standard SMRTA. Black line shows kgraphene while the black open squares give kGNR. For all cases, length L=3μm and T=300K. Reprinted with permission from [53]. Copyrighted by the American Physical Society.For 2D materials and nanotube structures, theN scattering is usually strong. For example, for CNT, Lindsay et al. [52] find that all the three-acoustic-phonon scatterings are N processes, so that the U process is, respectively, weak because it must involve optical phonons, which are less likely to be thermally excited. Thus, the Iterative Scheme can be used to layer and tube structure, rather than Standard SMRTA whose results are less accurate. Figure 7 shows the ratio between thermal conductivities kL (from Iterative Scheme) and kRTA (from Standard SMRTA), where the discrepancy is typically larger than 100%. The ZA mode shows the largest divergence which can reach 8-fold at length of 10 μm, because the flexural phonons have lower frequencies than other modes and thus stronger N process than U process.Figure 7
The ratio (color online) between thermal conductivity predicted from Iterative Scheme (kL) and Standard SMRTA (kRTA) of graphene as a function of length L for temperature 300 K, kZA (solid red), kTA (blue dashed), and kLA (green dotted) being the contributions from different branches. Reprinted with permission from [54]. Copyrighted by the American Physical Society.


#### 3.3.3. Boundary Scattering: Nanowires

For nanowires, the Casimir model (Table1) has been applied to predict thermal conductivity in many works [24–27, 27–32, 133] recently. Generally, k decreases with decreasing nanowire diameter; however, at some point, as the diameter continues to decrease, k will increase due to the 3D-1D transition. The main problems are that the results strongly rely on fitting parameters and that the use of Matthiessen approximation is still questioned. Instead, Ziman [95] presents an approach of solving space-dependent BTE (Peierls-BTE [134]). The final result of this Peierls-BTE approach gives, according to the simplification by Li et al. [63, 64], the position-dependent spectral phonon relaxation time
(22)τr,λ=τλ0(1+Δ¯λ){1-e-|(r-rb)/τλ0vλ|Gr,λ},
where Δ¯λ is the average value of Δλ over the cross section, rb is the point on the surface with r-rb being the same direction with group velocity vector vλ, and Gr,λ describes the boundary condition with Gr,λ=1 for completely diffusive and Gr,λ=0 for mirror like. The average of τr,λ over cross section is
(23)τ¯λ=τλ0(1+Δ¯λ)[1Sc∫Sc{1-e-|(r-rb)/τλ0vλ|Gr,λ}ds].
So far, the calculation for nanowire still needs an adjustable parameter to account for the boundary scattering.


#### 3.3.4. Impurity-Isotope Scattering: Doping and Alloys

From second-order perturbation theory [34, 95, 135], assuming that the isotopes are distributed randomly, the single scattering rate by the isotopes [82, 136] is given by
(24)Γλλ′iso=π2Ncωλωλ′∑bn‍gb|eλb·eλ′b*|2δ(ωλ-ωλ′),
where Nc and n stand for the number of unit cells and the number of atoms per unit cell, respectively, eλb is the eigenvector of λ mode in the basis atom b part, * denotes complex conjugate, and
(25)gb=∑i‍fib(1-mibm¯b)2
characterizes the magnitude of mass disorder, where i indicates isotope types, fib is the fraction of isotope i in lattice sites of basis atom b, mib is the mass of isotope i, and m¯b is the average atom mass of basis b sites. Sum over all the modes λ′ with λ′≠λ, the total scattering rate, or the inverse relaxation time of mode λ is
(26)1τλiso=∑λ′‍Γλλ′iso=π2Ncωλ2∑bn‍∑λ′‍gb|eλb·eλ′b*|2δ(ωλ-ωλ′).
For cubic symmetry system [82] such as Si, Ge, and Ar, the summation of eigenvectors in (26) can be reduced to
(27)1τλiso=π6nNcgωλ2∑λ′‍δ(ωλ-ωλ′)=πV06gω2D¯(ω)(28)=π2gω2D(ω),
where g is given by (5), V0=V/(nNc) is the volume per atom, D¯(ω)=∑λ′‍δ(ωλ-ωλ′)/V is density of states per unit volume, and D(ω) is density of states normalized to unity, noticing that the total amount of states is 3nNc and the total volume is V=V0nNc. Rewriting the D(ω) [82], one can obtain the formula in Table 1. From (27), the relaxation time of mode λ only depends on the frequency rather than the phonon branches.Equation (24) or (26) combined with the Standard SMRTA or Iterative Scheme have been used to predict the spectral phonon relaxation times of doped materials and even alloys. For isotope-doped system, this method can be directly applied, such as silicon and germanium [47, 60, 91], hexagonal boron nitride layers and nanotube [56, 58], GaN [57], and SiC [48]. For alloy, the disordered crystal is treated as an ordered one of the average atomic mass, lattice parameter, and force constants. This approach, the so-called virtual crystal approach, first introduced by Abeles [137], has been applied to Si-Ge alloys [65, 109] and PbTe(1−x)Sex alloys [106] byab initio IFCs, (Bi(1−x)Sbx)2Te3 alloys [138] with classical potential, and Ni0.55Pd0.45 alloys [139] for comparison with experiment. It turns out that the second-order perturbation ((24) or (26)) can give good prediction even for large mass disorder.


#### 3.3.5. Superlattices

Superlattices (SLs), composed of periodically arranged layers of two or more materials, have been extensively investigated in the aspect of thermal transport. Because of the heat transport suppression by interfaces and mass mismatch, superlattice has been designed to have a lower thermal conductivity than pure bulk. SLs are classified into two categories: diffuse and specular interfaces. The phonons in the first case are diffusively scattered by interfaces, while the phonons in the latter one propagate through the whole structure as if in one material, so-called coherent phonon transport [119]. Although proposed in theoretical studies long ago, the coherent phonon transport in SLs was not observed in experiments until Luckyanova et al. studied finite-thickness GaAs/AlAs SLs by time-domain thermoreflectance measurements [119]. It is found that the cross-plane k increases linearly with the number of periods when keeping the periods constants. Such phenomenon suggests that the phonon MFPs are equal to the sample thickness and the phonons do not “see” interfaces. Luckyanova et al. performed a first principle calculation (Standard SMRTA scheme) of GaAs/AlAs SLs to support their experimental results. They found that the anharmonic scattering rates and interface scattering rates, within the low-frequency region, had the frequency dependence as ~ω-2 and ~ω-4, respectively. The high-frequency phonons are scattered by interfaces, while the low-frequency phonons have long MFPs and thus can propagate though the entire SLs. Another evidence that the phonons in SL do not “see” interfaces is the fact that the accumulated k of GaAs/ALAs SL is similar to bulk GaAs as shown in Figure 5. All the calculated results support the experimental finding of coherent phonon propagations. In the following discussion of SLs, we only consider coherent phonon transport.Generally,k increases with increasing period length PL (at the limit PL→∞k increases to that of the pure bulk material). However, it is found that, for extremely short period length, k even increases with decreasing PL. This leads to a phenomenon that k as a function of PL reaches its minimum at a critical PL, and calculation of such value of PL is crucially important for designing low thermal conductivity materials. For instance, Yang et al. found that the isotope silicon superlattice isoSi/28Si nanowire had its lowest k at PL≈1nm [6]; Hu and Poulikakos noticed that the Si/Ge superlattice nanowire with 3.07nm diameter had its lowest k at PL≈4nm [7]; k of GaAs/AlAs superlattice [8] as a function of periodic length also obeys this principle. The exact phonon relaxation time explanation for such phenomenon is not available until the ALD method is explored [50, 59, 101, 103, 108].Garg et al. [108] studied short-period (0.3 nm) Si/Ge[001]1+1 superlattice using Standard SMRTA withab initio IFCs. They find that the thermal conductivity and phonon relaxation time of such superlattice are even greater than those of the two composition materials: pristine Si and Ge bulks. To understand this unusual behavior, the inverse relaxation time of TA mode is calculated and shown in Figure 8. Also plotted are the detailed three-phonon scattering rates for (a) TA + A→A, (b) TA + A→O, and (c) TA + O→O, where A and O stand for acoustic and optical, respectively. The “average material” is an imaginary material with averaged mass and potential of Si and Ge Bulk, for comparison with SiGe[001]1+1 superlattice. We note that only the (a) component can provide scattering for TA phonons, and that both (b) and (c) which affiliate with optical modes are almost completely absent. This indicates that the gap between optical and acoustic modes becomes so larger that the acoustic phonon can hardly be scattered by optical phonons. Such reduced scattering makes the relaxation times and thermal conductivity much larger than the two composition bulk materials.Figure 8
Scattering rate (color online) of TA mode at 300 K due to (a) absorption of an acoustic phonon to yield another acoustic phonon, (b) absorption of an acoustic phonon to yield an optical phonon, (c) absorption of an optical phonon to yield another optical phonon, and (d) total scattering rate, along (ξ,ξ,ξ) in 28Si, “average material,” and SiGe[001]1+1 superlattice. ξ is reduced wave vector, with k=(ξ,ξ,ξ)·2π/a. Reprinted with permission from [108]. Copyright (2011) American Chemical Society.More generally, Broido and Reinecke [59] and Ward and Broido [50] studied M1/M2[n+n] superlattice (the diamond structure with periodical n layers of mass M1 atoms and mass M2 atoms in [0,0,1] direction) using Iterative Scheme. In these two works, the IFCs are determined using Keating model [144, 145] and adiabatic bond charge (ABC) model [146, 147], respectively. In such M1/M2[n+n] superlattice, the high thermal conductivity is also found for n=1 (Figure 9). When n and mass ratio M1/M2 are increasing, k is determined by the competition between the decease of phonon group velocity and the increase of phonon relaxation time. It turns out that, from about M1/M2=2.3, the latter competitor dominates and thus k increase with increasing mass ratio. In Figure 9, as expected the k’s from Iterative Scheme are generally larger than those from Standard SMRTA method. This difference increases with increasing mass ratio because the occurrence of N process increases when the acoustic-optical gap gets larger.Figure 9
Thermal conductivities calculated by Standard SMRTA method (dashed lines) and Iterative Scheme (solid lines) of1+1, 2+2, and 4+4 Si/Ge-based SLs as a function of M1/M2, the mass ratio of constituent atoms. The thin vertical line shows the Ge/Si mass ratio. Reprinted with permission from [50]. Copyrighted by the American Physical Society.


## 3.1. Standard Single Mode Relaxation Time Approximation

The Standard SMRTA assumes that the system is in its complete thermal equilibrium, except that one phonon modeλ has its occupation number nλ=nλ0+nλ′ differing a small amount from its equilibrium value nλ0. Therefore, on the right hand side of (12), replacing nλ by nλ0+nλ′, whilst nλ′ and nλ′′ by nλ′0 and nλ′′0, respectively, one can obtain the phonon relaxation time τλ0 of mode λ (for the derivation, see Appendix A.1):
(15)1τλ0=∑λ′λ′′+‍Γλλ′λ′′++12∑λ′λ′′-‍Γλλ′λ′′-+∑λ′‍Γλλ′ext,
where the first two terms on the right hand side are intrinsic three-phonon scattering rates (ωλ±ωλ′=ωλ′′):
(16)Γλλ′λ′′±=ℏπ4N0{nλ′0-nλ′′0nλ′0+nλ′′0+1}×|V±(3)|2δ(ωλ±ωλ′-ωλ′′)ωλωλ′ωλ′′.
The last term Γλλ′ext represents the extrinsic scattering such as boundary scattering and impurity scattering.


## 3.2.  Iterative Scheme: Exact Solution to Linearized BTE

Different from the Standard SMRTA, the other method to solve the phonon BTE allows all the modes to be in their thermal nonequilibrium states at the same time. By replacing the occupation numbersnλ, nλ′, and nλ′′ by nλ0+nλ′, nλ′0+nλ′′, and nλ′′0+nλ′′′, respectively, on the right hand side of (12), the relaxation time τλ of mode λ is obtained (for the derivation see Appendix A.2)
(17)τλ=τλ0(1+Δλ),(18)Δλ=∑λ′λ′′(+)‍Γλλ′λ′′+(ξλλ′′τλ′′-ξλλ′τλ′)+∑λ′λ′′(-)‍12Γλλ′λ′′-(ξλλ′′τλ′′+ξλλ′τλ′)+∑λ′‍Γλλ′extξλλ′τλ′,
where ξλλ′=vλ′,zωλ′/vλ,zωλ, and vz is phonon group velocity component along the transport direction.Equation (17) is solved iteratively because both the left and the right hand sides contain the unknown variable τλ, and thus the method is called Iterative Scheme. This scheme is also based on RTA; thus (10) and (11) are still valid (one can reach this by substituting (A.3), (A.4), (A.11), (A.12), (A.13), and (A.14) into (9)). The last summation in (18) is done over λ′ with λ′≠λ.


## 3.3. Discussions and Applications

ALD methods can be divided into classical method andab initiomethod, differing in how to calculate the harmonic and anharmonic IFCs, which are the only inputs to these methods. The classical approach relies on empirical interatomic potential whose nth order derivatives are taken as the nth order IFCs:
(19)Φl1b1,…,lnbnα1⋯αn=∂nΦ∂u1(l1b1),…,∂un(lnbn).
In contrast, theab initio approach is a first principle calculation in the framework of density functional perturbation theory (DFPT) [43, 96, 97] using norm-conserving pseudopotentials in the local density approximation (LDA) without introducing any adjustable parameters. The formulism of the IFCs using first principle method can be found in [44] and realized by, for example, the QUANTUM ESPRESSO package [98]. Compared to the classical method, this method can deal with new materials whose empirical interatomic potentials are unknown. Further, this method can be more accurate since the empirical interatomic potentials cannot always represent the exact nature of interatomic force.In (16), the delta function δ(ω±ω′-ω′′) is typically approximated by δ(x)=lim⁡ε→0+(1/π)(ε/(x2+ε2)). To accurately evaluate (16), the choice of ε value is critical: it must be small but larger than the smallest increment in discrete x, which results from the use of finite grid of k points in Brillouin zone. The general practice is as follows: pick the densest grid possible and start with a sufficiently small guess, and increase it gradually until the final results reach convergence.To calculate the relaxation time, one can use Standard SMRTA scheme [43, 44, 99–111] or Iterative Scheme [47–65, 112], and, in each of them, one can choose empirical interatomic potential approach [47–49, 52–56, 58–60, 99–102] orab initio-derived IFC IFC [43, 44, 50, 51, 57, 61–65, 103–112]. The methods can be used on pure bulk, nanowires, doped bulk, doped nanowires, alloys, and so forth.One way to predict thermal conductivityk without working out all the phonon modes relaxation times is the Monte Carlo integration technique [101, 113]. The protocol of this technique is as follows: (1) randomly sample some phonon modes λ, (2) for each of these modes, randomly choose two other modes λ′ and λ′′ that interact with λ to calculate the relaxation time, and (3) select as many points as necessary to ensure that the statistical error is small enough in both cases. Monte Carlo technique only works for the Standard SMRTA scheme, since the Iterative Scheme requires the relaxation times of all the phonon modes to do iteration. Monte Carlo technique reduces the computational cost but lowers the accuracy.In addition to intrinsic phonon scatteringΓ±, extrinsic scattering 1/τλext plays an important role in nanostructures,
(20)1τλext=∑λ′‍Γλλ′ext,
such as boundary scattering 1/τλbs and impurity scattering 1/τλimp.


### 3.3.1. Intrinsic Phonon Scattering: Bulk Materials

Without any fitting parameters, Standard SMRTA withab initio approach can accurately predict spectral phonon relaxation times and thermal conductivities. Ward and Broido [110] checked the validity of some old approximations introduced in Section 2: (1) long-wave approximation for three-phonon scattering and (2) ignoring optical phonons, using silicon and germanium as examples. First, the values of matrix element |V±(3)|2, which govern the scattering strength Γ, fromab initio calculation for acoustic phonons are compared to those given by LWA. The percentage error of |V±(3)|2 is shown in Figure 1. We note that the LWA only works for the very low frequency ωeff<0.8 THz, while, for most part 0.8<ωeff<12 THz, the LWA gives large discrepancy, where ωeff≡(ωλωλ′ωλ′′)1/3 is the geometric average of the three-phonon frequencies. Second, the relaxation times of optical modes are found to only contribute less than 10% to the total thermal conductivity of silicon. However, ignoring optical modes is erroneous since the optical phonons are essential to provide channels for acoustic phonon scattering. The explicit calculation of millions of three-phonon scattering shows that optical phonons are involved in 50–60% of the total acoustic phonon-scattering processes in Si and Ge. Last, beyond the ω,T dependencies of τ listed in Table 1, which rely on many approximations, the ALD calculation can give more precise τω,T dependence. This is illustrated in the inset of Figure 1, the relaxation times of the LA phonons in Si for T=300K. By decomposing the total scattering into U process and N process, we find the U process has a stronger frequency dependence τ(U)(ω)~ω-4 than N process τ(N)(ω)~ω-2. The results also show that normal scattering governs the total relaxation time τ(0)(ω) at low frequency, while Umklapp scattering dominates at high frequency. Such ~ω-4 relation is not expected in the analytical models in Table 1.Figure 1
Percent error (color online) in|V±(3)|2 from the LWA compared to first principle for silicon at 300 K. Insert shows the normal (blue dashed curve), Umklapp (green dotted curve), and total (red solid curve) relaxation times for the LA phonons calculated from Standard SMRTA byab initio approach. Adapted with permission from [110]. Copyrighted by the American Physical Society.One flaw of the Standard SMRTA is that it does not grasp the interplay between theN process and U process. The right hand side of (15) can be decomposed as Γλ(N)+Γλ(U) (only consider intrinsic phonon scattering), according to whether they are N or U scattering events. The Standard SMRTA scheme treats the N process and U process as two independent scattering events and use Matthiessen’s rule to account for the total relaxation time 1/τλ0=1/τλ(N)+1/τλ(U), where τλ(N,U) is defined as τλ(N,U)≡1/Γλ(N,U). However, it is well know that N process does not contribute to thermal resistance directly. Instead, it affects the U process (low-frequency N scattering produces high-frequency phonons which boosts U process), and then the U process produces thermal resistance. This error can be remedied in the Iterative Scheme by doing the iteration in (17). Therefore, the Standard SMRTA scheme only works for the system, where U process dominates so that the N scattering makes little difference to U process as well as to thermal resistance [51, 110].For Si and Ge at room temperature where theU process is strong, the thermal conductivity predicted by Standard SMRTA scheme is only 5–10% smaller than that by Iterative Scheme [110], the latter shows excellent agreement with experiment (see Figure 1 of [61]).In contrast, theU scattering in diamond is much weaker [114–116] due to the much smaller phase space [117]. As a result, the thermal conductivity given by these two methods can differ by 50% at room temperature [110]. As shown in Figure 2, this discrepancy increases with decreasing temperature since the Umklapp scattering is weakened when temperature decreases. The thermal conductivity of diamond predicted by Iterative Scheme withab initio approach agrees excellently with experiment as shown in [110]. It is also noted that the Standard SMRTA scheme always underpredicts the thermal conductivity because it treats N process as an independent channel for thermal resistance. On the other hand, if the relaxation time for U process only is used in the calculation, the thermal conductivity is always overpredicted. This again confirms that the N process has an indirect and partial contribution to the thermal resistance.Figure 2
The calculated intrinsic lattice thermal conductivity of diamond for the Standard SMRTA (dashed line) and the Iterative Scheme (solid line), both byab initio approach. Dotted line shows percent error of the Standard SMRTA result compared to the Iterative Scheme solution. Reprinted with permission from [51]. Copyrighted by the American Physical Society.One important application of ALD calculation is to predict and understand the thermal conductivity of thermoelectric materials and help to design higher thermoelectric performance structures. Based on first principle calculation, Shiga et al. [104] obtain the frequency-dependent relaxation times of pristine PbTe bulk at 300 K as shown in Figure 3. At low-frequency region, TA phonons have longer relaxation times than LA phonons with τ’s exhibiting ~ω-2 dependence. Separating the scattering rates into those of normal and Umklapp processes, they find the relations τNormal~ω-2 and τUmklapp~ω-3, which again indicate that the normal process dominates low-frequency region while the Umklapp dominates high-frequency part. By further studying the participation of each phonon mode to the total scattering rates, they find that the low thermal conductivity of PbTe is attributed to the strong scattering of LA phonons by TO phonons and the small group velocity of TA phonons. Figure 4 compares phonon relaxation times of PbTe and PbSe [106]. Although the anharmonicity of PbSe is normally expected to be larger due to the larger average Gruneisen parameter reported from experiments [121], in this work, it is found that, for TA mode, the relaxation times of PbSe are substantially longer than those of PbTe. Surprisingly, the optical phonons are found to contribute as much as 25% for PbSe and 22% for PbTe to the total thermal conductivity at the temperature range 300–700 K. Motivated by the question that phonons with what kind of MFP contribute the most to the total thermal conductivity, the cumulative k’s as functions of phonon MFP are calculated by ALD method with first principle approach as shown in Figure 5. Silicon is found to have phonon MFPs which span 6 orders of magnitude (0–106 nm), while the thermal transport in diamond is dominated by the phonon with narrow range of MFP (0.4–2 μm). It is found that the phonons with MFP below 4 μm for silicon, 1.6 μm for GaAs, 120 nm for ZrCoSb, 20nm for PbSe, and 10 nm for PbTe contribute 80% of total thermal conductivity. GaAs/AlAs superlattice is found to have similar phonon MFP with bulk GaAs. The curves of the alloy Mg2Si0.6Sn0.4 and its pure phases Mg2Si and Mg2Sn cross at the intermediate MFPs. These results provide great guidance for experimental works. For example, the PbTe-PbSe alloys with size of nanoparticle below 10 nm are synthesized and found to lead to as much as 60% reduction to the thermal conductivity which provides large space for improving ZT [122].(a) Spectral phonon relaxation times of pristine PbTe bulk at 300 K by Standard SMRTA scheme with first principle IFCs. Relaxation times of (b) normal and (c) Umklapp processes, respectively. The solid lines plot (a)τ=6×1012ω-2, (b) τ=8×1012ω-2, and (c) τ=8×1024ω-3. Reprinted with permission from [104]. Copyrighted by the American Physical Society.
(a)
(b)
(c)Spectral phonon relaxation times of PbSe bulk (squares) and PbTe bulk (crosses) at 300 K by Standard SMRTA scheme with IFCs from first principle calculation: (a) TA, (b) LA, (c) TO, and (d) LO. Reprinted with permission from [106]. Copyrighted by the American Physical Society.
(a)
(b)
(c)
(d)Figure 5
The accumulated thermal conductivity of different bulk materials as a function of phonon mean free path at room temperature calculated from ALD by first principle approach. (Diamond [64]; Si (iterative) [64]; Si (standard) [107]; GaAs [118]; FaAs/AlAs [119]; PbTe [106]; PbSe [106]; PbTe0.5Se0.5 [106]; Mg2Si [63]; Mg2Sn [63]; Mg2Si0.6Sn0.4 [63]; ZrCoSb [120]).


### 3.3.2. Single and Few-Layer 2D Materials, Nanoribbons, and Nanotubes

For single- and multilayer 2D materials, the boundary scattering from the sides perpendicular to the transport direction is much weaker than for 3D systems [123], making the boundary scattering expression in Table 1 unsuitable. Instead, when studying single-/multilayer graphene (SLG/MLG) and graphite [54, 55], single-wall carbon nanotubes (SWCNTs) [52, 53], single-/multilayer boron nitride (SLBN/MLBN), and boron nitride nanotubes (BNNTs) [56, 58], Lindsay and Broido only consider the boundary scattering from the two ends in the transport direction and show that
(21)1τbs=2|vλ·z^|L
works well in accounting for the boundary scattering, with L being the length between boundaries in the transport direction z^. Such formula has been shown to give correct thermal conductivity values of nanotubes [124] and nanoribbons [125] in the ballistic limit (L→0) and diffusive limit (L→∞).Vibrations in 2D lattices are characterized by two types of phonons: those vibrating in the plane of layer (TA and LA) and those vibrating out of plane, so called flexural phonons (ZA and ZO). Lindsay et al. [54] find the selection rule for all orders in anharmonic phonon-phonon scattering in the 2D crystals:only even numbers (including zero) of flexural phonons can be involved, arising from the reflection symmetry perpendicular to the plane of layer. This selection rule has forbidden about 60% of both N and U three-phonon scattering phase space of ZA phonons for single layer graphene. They show that such suppressed scattering yields long relaxation time and mean free path for ZA phonons, leading to ZA phonons contributing most of the thermal conductivity of SLG, about 70% at room temperature (another cause being the large density of states and occupation number of ZA modes). However, this conclusion is still under debate since this approach does not include the fourth- and higher-order phonon scattering rates, which are not necessarily low since the reflection symmetry allows more 4-phonon processes than 3-phonon processes. Actually, the method of spectral energy analysis based on MD (discussed in Section 4) indicates that only 25%–30% of the total k is contributed by ZA mode at room temperature [126–128]. It should be noted that MD has its own drawback of not reproducing the Bose-Einstein distribution for graphene phonons at room temperature. Hence, the discrepancies between the two methods still need further study.The selection rule mentioned above does not hold for multilayer graphene, twisted graphene, graphite (because of the interlayer coupling), CNT (due to the curvature), graphene nanoribbon (GNR) (due to boundary scattering), substrate-supported graphene (due to scattering with the substrate), and defective graphene (due to defective scattering). Therefore, the thermal conductivity of these structures is typically lower than that of single layer graphene, and the contribution of each phonon mode changes [54, 129–132]. In Figure 6, single-layer graphene, GNR, and SWCNT are compared, where graphene has an infinite width and finite length L, SWCNT has a finite diameter d and length L, and GNR has a finite width W=πd with artificial periodic boundary condition applied. As expected, kRTA underpredicted thermal conductivity. SWCNT is found to have a lower thermal conductivity than graphene with a minimum value of 77% of kgraphene at a critical diameter d≈1.5nm. From this critical diameter, kSWCNT increases with increasing diameter and reaches 90% of kgraphene at d≈4 nm. On the other hand, if d goes small enough, phonon-phonon scattering decreases and the thermal conductivity increases. At this short limit of d, the system becomes more like a 1D chain which generally has much larger thermal conductivity than 2D and 3D systems. The increasing trend of kGNR with decreasing d comes from the reason that the decrease of the width πd pushes the optical modes to higher frequencies and thus the U scattering by optical phonons becomes weaker.Figure 6
Thermal conductivityk versus diameter d (color online) for single wall carbon nanotubes. Solid red circles, blue squares, and green triangles represent zigzag, armchair, and chiral predicted by Iterative Scheme; open red circles, blue squares, and green triangles are those from Standard SMRTA. Black line shows kgraphene while the black open squares give kGNR. For all cases, length L=3μm and T=300K. Reprinted with permission from [53]. Copyrighted by the American Physical Society.For 2D materials and nanotube structures, theN scattering is usually strong. For example, for CNT, Lindsay et al. [52] find that all the three-acoustic-phonon scatterings are N processes, so that the U process is, respectively, weak because it must involve optical phonons, which are less likely to be thermally excited. Thus, the Iterative Scheme can be used to layer and tube structure, rather than Standard SMRTA whose results are less accurate. Figure 7 shows the ratio between thermal conductivities kL (from Iterative Scheme) and kRTA (from Standard SMRTA), where the discrepancy is typically larger than 100%. The ZA mode shows the largest divergence which can reach 8-fold at length of 10 μm, because the flexural phonons have lower frequencies than other modes and thus stronger N process than U process.Figure 7
The ratio (color online) between thermal conductivity predicted from Iterative Scheme (kL) and Standard SMRTA (kRTA) of graphene as a function of length L for temperature 300 K, kZA (solid red), kTA (blue dashed), and kLA (green dotted) being the contributions from different branches. Reprinted with permission from [54]. Copyrighted by the American Physical Society.


### 3.3.3. Boundary Scattering: Nanowires

For nanowires, the Casimir model (Table1) has been applied to predict thermal conductivity in many works [24–27, 27–32, 133] recently. Generally, k decreases with decreasing nanowire diameter; however, at some point, as the diameter continues to decrease, k will increase due to the 3D-1D transition. The main problems are that the results strongly rely on fitting parameters and that the use of Matthiessen approximation is still questioned. Instead, Ziman [95] presents an approach of solving space-dependent BTE (Peierls-BTE [134]). The final result of this Peierls-BTE approach gives, according to the simplification by Li et al. [63, 64], the position-dependent spectral phonon relaxation time
(22)τr,λ=τλ0(1+Δ¯λ){1-e-|(r-rb)/τλ0vλ|Gr,λ},
where Δ¯λ is the average value of Δλ over the cross section, rb is the point on the surface with r-rb being the same direction with group velocity vector vλ, and Gr,λ describes the boundary condition with Gr,λ=1 for completely diffusive and Gr,λ=0 for mirror like. The average of τr,λ over cross section is
(23)τ¯λ=τλ0(1+Δ¯λ)[1Sc∫Sc{1-e-|(r-rb)/τλ0vλ|Gr,λ}ds].
So far, the calculation for nanowire still needs an adjustable parameter to account for the boundary scattering.


### 3.3.4. Impurity-Isotope Scattering: Doping and Alloys

From second-order perturbation theory [34, 95, 135], assuming that the isotopes are distributed randomly, the single scattering rate by the isotopes [82, 136] is given by
(24)Γλλ′iso=π2Ncωλωλ′∑bn‍gb|eλb·eλ′b*|2δ(ωλ-ωλ′),
where Nc and n stand for the number of unit cells and the number of atoms per unit cell, respectively, eλb is the eigenvector of λ mode in the basis atom b part, * denotes complex conjugate, and
(25)gb=∑i‍fib(1-mibm¯b)2
characterizes the magnitude of mass disorder, where i indicates isotope types, fib is the fraction of isotope i in lattice sites of basis atom b, mib is the mass of isotope i, and m¯b is the average atom mass of basis b sites. Sum over all the modes λ′ with λ′≠λ, the total scattering rate, or the inverse relaxation time of mode λ is
(26)1τλiso=∑λ′‍Γλλ′iso=π2Ncωλ2∑bn‍∑λ′‍gb|eλb·eλ′b*|2δ(ωλ-ωλ′).
For cubic symmetry system [82] such as Si, Ge, and Ar, the summation of eigenvectors in (26) can be reduced to
(27)1τλiso=π6nNcgωλ2∑λ′‍δ(ωλ-ωλ′)=πV06gω2D¯(ω)(28)=π2gω2D(ω),
where g is given by (5), V0=V/(nNc) is the volume per atom, D¯(ω)=∑λ′‍δ(ωλ-ωλ′)/V is density of states per unit volume, and D(ω) is density of states normalized to unity, noticing that the total amount of states is 3nNc and the total volume is V=V0nNc. Rewriting the D(ω) [82], one can obtain the formula in Table 1. From (27), the relaxation time of mode λ only depends on the frequency rather than the phonon branches.Equation (24) or (26) combined with the Standard SMRTA or Iterative Scheme have been used to predict the spectral phonon relaxation times of doped materials and even alloys. For isotope-doped system, this method can be directly applied, such as silicon and germanium [47, 60, 91], hexagonal boron nitride layers and nanotube [56, 58], GaN [57], and SiC [48]. For alloy, the disordered crystal is treated as an ordered one of the average atomic mass, lattice parameter, and force constants. This approach, the so-called virtual crystal approach, first introduced by Abeles [137], has been applied to Si-Ge alloys [65, 109] and PbTe(1−x)Sex alloys [106] byab initio IFCs, (Bi(1−x)Sbx)2Te3 alloys [138] with classical potential, and Ni0.55Pd0.45 alloys [139] for comparison with experiment. It turns out that the second-order perturbation ((24) or (26)) can give good prediction even for large mass disorder.


### 3.3.5. Superlattices

Superlattices (SLs), composed of periodically arranged layers of two or more materials, have been extensively investigated in the aspect of thermal transport. Because of the heat transport suppression by interfaces and mass mismatch, superlattice has been designed to have a lower thermal conductivity than pure bulk. SLs are classified into two categories: diffuse and specular interfaces. The phonons in the first case are diffusively scattered by interfaces, while the phonons in the latter one propagate through the whole structure as if in one material, so-called coherent phonon transport [119]. Although proposed in theoretical studies long ago, the coherent phonon transport in SLs was not observed in experiments until Luckyanova et al. studied finite-thickness GaAs/AlAs SLs by time-domain thermoreflectance measurements [119]. It is found that the cross-plane k increases linearly with the number of periods when keeping the periods constants. Such phenomenon suggests that the phonon MFPs are equal to the sample thickness and the phonons do not “see” interfaces. Luckyanova et al. performed a first principle calculation (Standard SMRTA scheme) of GaAs/AlAs SLs to support their experimental results. They found that the anharmonic scattering rates and interface scattering rates, within the low-frequency region, had the frequency dependence as ~ω-2 and ~ω-4, respectively. The high-frequency phonons are scattered by interfaces, while the low-frequency phonons have long MFPs and thus can propagate though the entire SLs. Another evidence that the phonons in SL do not “see” interfaces is the fact that the accumulated k of GaAs/ALAs SL is similar to bulk GaAs as shown in Figure 5. All the calculated results support the experimental finding of coherent phonon propagations. In the following discussion of SLs, we only consider coherent phonon transport.Generally,k increases with increasing period length PL (at the limit PL→∞k increases to that of the pure bulk material). However, it is found that, for extremely short period length, k even increases with decreasing PL. This leads to a phenomenon that k as a function of PL reaches its minimum at a critical PL, and calculation of such value of PL is crucially important for designing low thermal conductivity materials. For instance, Yang et al. found that the isotope silicon superlattice isoSi/28Si nanowire had its lowest k at PL≈1nm [6]; Hu and Poulikakos noticed that the Si/Ge superlattice nanowire with 3.07nm diameter had its lowest k at PL≈4nm [7]; k of GaAs/AlAs superlattice [8] as a function of periodic length also obeys this principle. The exact phonon relaxation time explanation for such phenomenon is not available until the ALD method is explored [50, 59, 101, 103, 108].Garg et al. [108] studied short-period (0.3 nm) Si/Ge[001]1+1 superlattice using Standard SMRTA withab initio IFCs. They find that the thermal conductivity and phonon relaxation time of such superlattice are even greater than those of the two composition materials: pristine Si and Ge bulks. To understand this unusual behavior, the inverse relaxation time of TA mode is calculated and shown in Figure 8. Also plotted are the detailed three-phonon scattering rates for (a) TA + A→A, (b) TA + A→O, and (c) TA + O→O, where A and O stand for acoustic and optical, respectively. The “average material” is an imaginary material with averaged mass and potential of Si and Ge Bulk, for comparison with SiGe[001]1+1 superlattice. We note that only the (a) component can provide scattering for TA phonons, and that both (b) and (c) which affiliate with optical modes are almost completely absent. This indicates that the gap between optical and acoustic modes becomes so larger that the acoustic phonon can hardly be scattered by optical phonons. Such reduced scattering makes the relaxation times and thermal conductivity much larger than the two composition bulk materials.Figure 8
Scattering rate (color online) of TA mode at 300 K due to (a) absorption of an acoustic phonon to yield another acoustic phonon, (b) absorption of an acoustic phonon to yield an optical phonon, (c) absorption of an optical phonon to yield another optical phonon, and (d) total scattering rate, along (ξ,ξ,ξ) in 28Si, “average material,” and SiGe[001]1+1 superlattice. ξ is reduced wave vector, with k=(ξ,ξ,ξ)·2π/a. Reprinted with permission from [108]. Copyright (2011) American Chemical Society.More generally, Broido and Reinecke [59] and Ward and Broido [50] studied M1/M2[n+n] superlattice (the diamond structure with periodical n layers of mass M1 atoms and mass M2 atoms in [0,0,1] direction) using Iterative Scheme. In these two works, the IFCs are determined using Keating model [144, 145] and adiabatic bond charge (ABC) model [146, 147], respectively. In such M1/M2[n+n] superlattice, the high thermal conductivity is also found for n=1 (Figure 9). When n and mass ratio M1/M2 are increasing, k is determined by the competition between the decease of phonon group velocity and the increase of phonon relaxation time. It turns out that, from about M1/M2=2.3, the latter competitor dominates and thus k increase with increasing mass ratio. In Figure 9, as expected the k’s from Iterative Scheme are generally larger than those from Standard SMRTA method. This difference increases with increasing mass ratio because the occurrence of N process increases when the acoustic-optical gap gets larger.Figure 9
Thermal conductivities calculated by Standard SMRTA method (dashed lines) and Iterative Scheme (solid lines) of1+1, 2+2, and 4+4 Si/Ge-based SLs as a function of M1/M2, the mass ratio of constituent atoms. The thin vertical line shows the Ge/Si mass ratio. Reprinted with permission from [50]. Copyrighted by the American Physical Society.


## 3.3.1. Intrinsic Phonon Scattering: Bulk Materials

Without any fitting parameters, Standard SMRTA withab initio approach can accurately predict spectral phonon relaxation times and thermal conductivities. Ward and Broido [110] checked the validity of some old approximations introduced in Section 2: (1) long-wave approximation for three-phonon scattering and (2) ignoring optical phonons, using silicon and germanium as examples. First, the values of matrix element |V±(3)|2, which govern the scattering strength Γ, fromab initio calculation for acoustic phonons are compared to those given by LWA. The percentage error of |V±(3)|2 is shown in Figure 1. We note that the LWA only works for the very low frequency ωeff<0.8 THz, while, for most part 0.8<ωeff<12 THz, the LWA gives large discrepancy, where ωeff≡(ωλωλ′ωλ′′)1/3 is the geometric average of the three-phonon frequencies. Second, the relaxation times of optical modes are found to only contribute less than 10% to the total thermal conductivity of silicon. However, ignoring optical modes is erroneous since the optical phonons are essential to provide channels for acoustic phonon scattering. The explicit calculation of millions of three-phonon scattering shows that optical phonons are involved in 50–60% of the total acoustic phonon-scattering processes in Si and Ge. Last, beyond the ω,T dependencies of τ listed in Table 1, which rely on many approximations, the ALD calculation can give more precise τω,T dependence. This is illustrated in the inset of Figure 1, the relaxation times of the LA phonons in Si for T=300K. By decomposing the total scattering into U process and N process, we find the U process has a stronger frequency dependence τ(U)(ω)~ω-4 than N process τ(N)(ω)~ω-2. The results also show that normal scattering governs the total relaxation time τ(0)(ω) at low frequency, while Umklapp scattering dominates at high frequency. Such ~ω-4 relation is not expected in the analytical models in Table 1.Figure 1
Percent error (color online) in|V±(3)|2 from the LWA compared to first principle for silicon at 300 K. Insert shows the normal (blue dashed curve), Umklapp (green dotted curve), and total (red solid curve) relaxation times for the LA phonons calculated from Standard SMRTA byab initio approach. Adapted with permission from [110]. Copyrighted by the American Physical Society.One flaw of the Standard SMRTA is that it does not grasp the interplay between theN process and U process. The right hand side of (15) can be decomposed as Γλ(N)+Γλ(U) (only consider intrinsic phonon scattering), according to whether they are N or U scattering events. The Standard SMRTA scheme treats the N process and U process as two independent scattering events and use Matthiessen’s rule to account for the total relaxation time 1/τλ0=1/τλ(N)+1/τλ(U), where τλ(N,U) is defined as τλ(N,U)≡1/Γλ(N,U). However, it is well know that N process does not contribute to thermal resistance directly. Instead, it affects the U process (low-frequency N scattering produces high-frequency phonons which boosts U process), and then the U process produces thermal resistance. This error can be remedied in the Iterative Scheme by doing the iteration in (17). Therefore, the Standard SMRTA scheme only works for the system, where U process dominates so that the N scattering makes little difference to U process as well as to thermal resistance [51, 110].For Si and Ge at room temperature where theU process is strong, the thermal conductivity predicted by Standard SMRTA scheme is only 5–10% smaller than that by Iterative Scheme [110], the latter shows excellent agreement with experiment (see Figure 1 of [61]).In contrast, theU scattering in diamond is much weaker [114–116] due to the much smaller phase space [117]. As a result, the thermal conductivity given by these two methods can differ by 50% at room temperature [110]. As shown in Figure 2, this discrepancy increases with decreasing temperature since the Umklapp scattering is weakened when temperature decreases. The thermal conductivity of diamond predicted by Iterative Scheme withab initio approach agrees excellently with experiment as shown in [110]. It is also noted that the Standard SMRTA scheme always underpredicts the thermal conductivity because it treats N process as an independent channel for thermal resistance. On the other hand, if the relaxation time for U process only is used in the calculation, the thermal conductivity is always overpredicted. This again confirms that the N process has an indirect and partial contribution to the thermal resistance.Figure 2
The calculated intrinsic lattice thermal conductivity of diamond for the Standard SMRTA (dashed line) and the Iterative Scheme (solid line), both byab initio approach. Dotted line shows percent error of the Standard SMRTA result compared to the Iterative Scheme solution. Reprinted with permission from [51]. Copyrighted by the American Physical Society.One important application of ALD calculation is to predict and understand the thermal conductivity of thermoelectric materials and help to design higher thermoelectric performance structures. Based on first principle calculation, Shiga et al. [104] obtain the frequency-dependent relaxation times of pristine PbTe bulk at 300 K as shown in Figure 3. At low-frequency region, TA phonons have longer relaxation times than LA phonons with τ’s exhibiting ~ω-2 dependence. Separating the scattering rates into those of normal and Umklapp processes, they find the relations τNormal~ω-2 and τUmklapp~ω-3, which again indicate that the normal process dominates low-frequency region while the Umklapp dominates high-frequency part. By further studying the participation of each phonon mode to the total scattering rates, they find that the low thermal conductivity of PbTe is attributed to the strong scattering of LA phonons by TO phonons and the small group velocity of TA phonons. Figure 4 compares phonon relaxation times of PbTe and PbSe [106]. Although the anharmonicity of PbSe is normally expected to be larger due to the larger average Gruneisen parameter reported from experiments [121], in this work, it is found that, for TA mode, the relaxation times of PbSe are substantially longer than those of PbTe. Surprisingly, the optical phonons are found to contribute as much as 25% for PbSe and 22% for PbTe to the total thermal conductivity at the temperature range 300–700 K. Motivated by the question that phonons with what kind of MFP contribute the most to the total thermal conductivity, the cumulative k’s as functions of phonon MFP are calculated by ALD method with first principle approach as shown in Figure 5. Silicon is found to have phonon MFPs which span 6 orders of magnitude (0–106 nm), while the thermal transport in diamond is dominated by the phonon with narrow range of MFP (0.4–2 μm). It is found that the phonons with MFP below 4 μm for silicon, 1.6 μm for GaAs, 120 nm for ZrCoSb, 20nm for PbSe, and 10 nm for PbTe contribute 80% of total thermal conductivity. GaAs/AlAs superlattice is found to have similar phonon MFP with bulk GaAs. The curves of the alloy Mg2Si0.6Sn0.4 and its pure phases Mg2Si and Mg2Sn cross at the intermediate MFPs. These results provide great guidance for experimental works. For example, the PbTe-PbSe alloys with size of nanoparticle below 10 nm are synthesized and found to lead to as much as 60% reduction to the thermal conductivity which provides large space for improving ZT [122].(a) Spectral phonon relaxation times of pristine PbTe bulk at 300 K by Standard SMRTA scheme with first principle IFCs. Relaxation times of (b) normal and (c) Umklapp processes, respectively. The solid lines plot (a)τ=6×1012ω-2, (b) τ=8×1012ω-2, and (c) τ=8×1024ω-3. Reprinted with permission from [104]. Copyrighted by the American Physical Society.
(a)
(b)
(c)Spectral phonon relaxation times of PbSe bulk (squares) and PbTe bulk (crosses) at 300 K by Standard SMRTA scheme with IFCs from first principle calculation: (a) TA, (b) LA, (c) TO, and (d) LO. Reprinted with permission from [106]. Copyrighted by the American Physical Society.
(a)
(b)
(c)
(d)Figure 5
The accumulated thermal conductivity of different bulk materials as a function of phonon mean free path at room temperature calculated from ALD by first principle approach. (Diamond [64]; Si (iterative) [64]; Si (standard) [107]; GaAs [118]; FaAs/AlAs [119]; PbTe [106]; PbSe [106]; PbTe0.5Se0.5 [106]; Mg2Si [63]; Mg2Sn [63]; Mg2Si0.6Sn0.4 [63]; ZrCoSb [120]).


## 3.3.2. Single and Few-Layer 2D Materials, Nanoribbons, and Nanotubes

For single- and multilayer 2D materials, the boundary scattering from the sides perpendicular to the transport direction is much weaker than for 3D systems [123], making the boundary scattering expression in Table 1 unsuitable. Instead, when studying single-/multilayer graphene (SLG/MLG) and graphite [54, 55], single-wall carbon nanotubes (SWCNTs) [52, 53], single-/multilayer boron nitride (SLBN/MLBN), and boron nitride nanotubes (BNNTs) [56, 58], Lindsay and Broido only consider the boundary scattering from the two ends in the transport direction and show that
(21)1τbs=2|vλ·z^|L
works well in accounting for the boundary scattering, with L being the length between boundaries in the transport direction z^. Such formula has been shown to give correct thermal conductivity values of nanotubes [124] and nanoribbons [125] in the ballistic limit (L→0) and diffusive limit (L→∞).Vibrations in 2D lattices are characterized by two types of phonons: those vibrating in the plane of layer (TA and LA) and those vibrating out of plane, so called flexural phonons (ZA and ZO). Lindsay et al. [54] find the selection rule for all orders in anharmonic phonon-phonon scattering in the 2D crystals:only even numbers (including zero) of flexural phonons can be involved, arising from the reflection symmetry perpendicular to the plane of layer. This selection rule has forbidden about 60% of both N and U three-phonon scattering phase space of ZA phonons for single layer graphene. They show that such suppressed scattering yields long relaxation time and mean free path for ZA phonons, leading to ZA phonons contributing most of the thermal conductivity of SLG, about 70% at room temperature (another cause being the large density of states and occupation number of ZA modes). However, this conclusion is still under debate since this approach does not include the fourth- and higher-order phonon scattering rates, which are not necessarily low since the reflection symmetry allows more 4-phonon processes than 3-phonon processes. Actually, the method of spectral energy analysis based on MD (discussed in Section 4) indicates that only 25%–30% of the total k is contributed by ZA mode at room temperature [126–128]. It should be noted that MD has its own drawback of not reproducing the Bose-Einstein distribution for graphene phonons at room temperature. Hence, the discrepancies between the two methods still need further study.The selection rule mentioned above does not hold for multilayer graphene, twisted graphene, graphite (because of the interlayer coupling), CNT (due to the curvature), graphene nanoribbon (GNR) (due to boundary scattering), substrate-supported graphene (due to scattering with the substrate), and defective graphene (due to defective scattering). Therefore, the thermal conductivity of these structures is typically lower than that of single layer graphene, and the contribution of each phonon mode changes [54, 129–132]. In Figure 6, single-layer graphene, GNR, and SWCNT are compared, where graphene has an infinite width and finite length L, SWCNT has a finite diameter d and length L, and GNR has a finite width W=πd with artificial periodic boundary condition applied. As expected, kRTA underpredicted thermal conductivity. SWCNT is found to have a lower thermal conductivity than graphene with a minimum value of 77% of kgraphene at a critical diameter d≈1.5nm. From this critical diameter, kSWCNT increases with increasing diameter and reaches 90% of kgraphene at d≈4 nm. On the other hand, if d goes small enough, phonon-phonon scattering decreases and the thermal conductivity increases. At this short limit of d, the system becomes more like a 1D chain which generally has much larger thermal conductivity than 2D and 3D systems. The increasing trend of kGNR with decreasing d comes from the reason that the decrease of the width πd pushes the optical modes to higher frequencies and thus the U scattering by optical phonons becomes weaker.Figure 6
Thermal conductivityk versus diameter d (color online) for single wall carbon nanotubes. Solid red circles, blue squares, and green triangles represent zigzag, armchair, and chiral predicted by Iterative Scheme; open red circles, blue squares, and green triangles are those from Standard SMRTA. Black line shows kgraphene while the black open squares give kGNR. For all cases, length L=3μm and T=300K. Reprinted with permission from [53]. Copyrighted by the American Physical Society.For 2D materials and nanotube structures, theN scattering is usually strong. For example, for CNT, Lindsay et al. [52] find that all the three-acoustic-phonon scatterings are N processes, so that the U process is, respectively, weak because it must involve optical phonons, which are less likely to be thermally excited. Thus, the Iterative Scheme can be used to layer and tube structure, rather than Standard SMRTA whose results are less accurate. Figure 7 shows the ratio between thermal conductivities kL (from Iterative Scheme) and kRTA (from Standard SMRTA), where the discrepancy is typically larger than 100%. The ZA mode shows the largest divergence which can reach 8-fold at length of 10 μm, because the flexural phonons have lower frequencies than other modes and thus stronger N process than U process.Figure 7
The ratio (color online) between thermal conductivity predicted from Iterative Scheme (kL) and Standard SMRTA (kRTA) of graphene as a function of length L for temperature 300 K, kZA (solid red), kTA (blue dashed), and kLA (green dotted) being the contributions from different branches. Reprinted with permission from [54]. Copyrighted by the American Physical Society.


## 3.3.3. Boundary Scattering: Nanowires

For nanowires, the Casimir model (Table1) has been applied to predict thermal conductivity in many works [24–27, 27–32, 133] recently. Generally, k decreases with decreasing nanowire diameter; however, at some point, as the diameter continues to decrease, k will increase due to the 3D-1D transition. The main problems are that the results strongly rely on fitting parameters and that the use of Matthiessen approximation is still questioned. Instead, Ziman [95] presents an approach of solving space-dependent BTE (Peierls-BTE [134]). The final result of this Peierls-BTE approach gives, according to the simplification by Li et al. [63, 64], the position-dependent spectral phonon relaxation time
(22)τr,λ=τλ0(1+Δ¯λ){1-e-|(r-rb)/τλ0vλ|Gr,λ},
where Δ¯λ is the average value of Δλ over the cross section, rb is the point on the surface with r-rb being the same direction with group velocity vector vλ, and Gr,λ describes the boundary condition with Gr,λ=1 for completely diffusive and Gr,λ=0 for mirror like. The average of τr,λ over cross section is
(23)τ¯λ=τλ0(1+Δ¯λ)[1Sc∫Sc{1-e-|(r-rb)/τλ0vλ|Gr,λ}ds].
So far, the calculation for nanowire still needs an adjustable parameter to account for the boundary scattering.


## 3.3.4. Impurity-Isotope Scattering: Doping and Alloys

From second-order perturbation theory [34, 95, 135], assuming that the isotopes are distributed randomly, the single scattering rate by the isotopes [82, 136] is given by
(24)Γλλ′iso=π2Ncωλωλ′∑bn‍gb|eλb·eλ′b*|2δ(ωλ-ωλ′),
where Nc and n stand for the number of unit cells and the number of atoms per unit cell, respectively, eλb is the eigenvector of λ mode in the basis atom b part, * denotes complex conjugate, and
(25)gb=∑i‍fib(1-mibm¯b)2
characterizes the magnitude of mass disorder, where i indicates isotope types, fib is the fraction of isotope i in lattice sites of basis atom b, mib is the mass of isotope i, and m¯b is the average atom mass of basis b sites. Sum over all the modes λ′ with λ′≠λ, the total scattering rate, or the inverse relaxation time of mode λ is
(26)1τλiso=∑λ′‍Γλλ′iso=π2Ncωλ2∑bn‍∑λ′‍gb|eλb·eλ′b*|2δ(ωλ-ωλ′).
For cubic symmetry system [82] such as Si, Ge, and Ar, the summation of eigenvectors in (26) can be reduced to
(27)1τλiso=π6nNcgωλ2∑λ′‍δ(ωλ-ωλ′)=πV06gω2D¯(ω)(28)=π2gω2D(ω),
where g is given by (5), V0=V/(nNc) is the volume per atom, D¯(ω)=∑λ′‍δ(ωλ-ωλ′)/V is density of states per unit volume, and D(ω) is density of states normalized to unity, noticing that the total amount of states is 3nNc and the total volume is V=V0nNc. Rewriting the D(ω) [82], one can obtain the formula in Table 1. From (27), the relaxation time of mode λ only depends on the frequency rather than the phonon branches.Equation (24) or (26) combined with the Standard SMRTA or Iterative Scheme have been used to predict the spectral phonon relaxation times of doped materials and even alloys. For isotope-doped system, this method can be directly applied, such as silicon and germanium [47, 60, 91], hexagonal boron nitride layers and nanotube [56, 58], GaN [57], and SiC [48]. For alloy, the disordered crystal is treated as an ordered one of the average atomic mass, lattice parameter, and force constants. This approach, the so-called virtual crystal approach, first introduced by Abeles [137], has been applied to Si-Ge alloys [65, 109] and PbTe(1−x)Sex alloys [106] byab initio IFCs, (Bi(1−x)Sbx)2Te3 alloys [138] with classical potential, and Ni0.55Pd0.45 alloys [139] for comparison with experiment. It turns out that the second-order perturbation ((24) or (26)) can give good prediction even for large mass disorder.


## 3.3.5. Superlattices

Superlattices (SLs), composed of periodically arranged layers of two or more materials, have been extensively investigated in the aspect of thermal transport. Because of the heat transport suppression by interfaces and mass mismatch, superlattice has been designed to have a lower thermal conductivity than pure bulk. SLs are classified into two categories: diffuse and specular interfaces. The phonons in the first case are diffusively scattered by interfaces, while the phonons in the latter one propagate through the whole structure as if in one material, so-called coherent phonon transport [119]. Although proposed in theoretical studies long ago, the coherent phonon transport in SLs was not observed in experiments until Luckyanova et al. studied finite-thickness GaAs/AlAs SLs by time-domain thermoreflectance measurements [119]. It is found that the cross-plane k increases linearly with the number of periods when keeping the periods constants. Such phenomenon suggests that the phonon MFPs are equal to the sample thickness and the phonons do not “see” interfaces. Luckyanova et al. performed a first principle calculation (Standard SMRTA scheme) of GaAs/AlAs SLs to support their experimental results. They found that the anharmonic scattering rates and interface scattering rates, within the low-frequency region, had the frequency dependence as ~ω-2 and ~ω-4, respectively. The high-frequency phonons are scattered by interfaces, while the low-frequency phonons have long MFPs and thus can propagate though the entire SLs. Another evidence that the phonons in SL do not “see” interfaces is the fact that the accumulated k of GaAs/ALAs SL is similar to bulk GaAs as shown in Figure 5. All the calculated results support the experimental finding of coherent phonon propagations. In the following discussion of SLs, we only consider coherent phonon transport.Generally,k increases with increasing period length PL (at the limit PL→∞k increases to that of the pure bulk material). However, it is found that, for extremely short period length, k even increases with decreasing PL. This leads to a phenomenon that k as a function of PL reaches its minimum at a critical PL, and calculation of such value of PL is crucially important for designing low thermal conductivity materials. For instance, Yang et al. found that the isotope silicon superlattice isoSi/28Si nanowire had its lowest k at PL≈1nm [6]; Hu and Poulikakos noticed that the Si/Ge superlattice nanowire with 3.07nm diameter had its lowest k at PL≈4nm [7]; k of GaAs/AlAs superlattice [8] as a function of periodic length also obeys this principle. The exact phonon relaxation time explanation for such phenomenon is not available until the ALD method is explored [50, 59, 101, 103, 108].Garg et al. [108] studied short-period (0.3 nm) Si/Ge[001]1+1 superlattice using Standard SMRTA withab initio IFCs. They find that the thermal conductivity and phonon relaxation time of such superlattice are even greater than those of the two composition materials: pristine Si and Ge bulks. To understand this unusual behavior, the inverse relaxation time of TA mode is calculated and shown in Figure 8. Also plotted are the detailed three-phonon scattering rates for (a) TA + A→A, (b) TA + A→O, and (c) TA + O→O, where A and O stand for acoustic and optical, respectively. The “average material” is an imaginary material with averaged mass and potential of Si and Ge Bulk, for comparison with SiGe[001]1+1 superlattice. We note that only the (a) component can provide scattering for TA phonons, and that both (b) and (c) which affiliate with optical modes are almost completely absent. This indicates that the gap between optical and acoustic modes becomes so larger that the acoustic phonon can hardly be scattered by optical phonons. Such reduced scattering makes the relaxation times and thermal conductivity much larger than the two composition bulk materials.Figure 8
Scattering rate (color online) of TA mode at 300 K due to (a) absorption of an acoustic phonon to yield another acoustic phonon, (b) absorption of an acoustic phonon to yield an optical phonon, (c) absorption of an optical phonon to yield another optical phonon, and (d) total scattering rate, along (ξ,ξ,ξ) in 28Si, “average material,” and SiGe[001]1+1 superlattice. ξ is reduced wave vector, with k=(ξ,ξ,ξ)·2π/a. Reprinted with permission from [108]. Copyright (2011) American Chemical Society.More generally, Broido and Reinecke [59] and Ward and Broido [50] studied M1/M2[n+n] superlattice (the diamond structure with periodical n layers of mass M1 atoms and mass M2 atoms in [0,0,1] direction) using Iterative Scheme. In these two works, the IFCs are determined using Keating model [144, 145] and adiabatic bond charge (ABC) model [146, 147], respectively. In such M1/M2[n+n] superlattice, the high thermal conductivity is also found for n=1 (Figure 9). When n and mass ratio M1/M2 are increasing, k is determined by the competition between the decease of phonon group velocity and the increase of phonon relaxation time. It turns out that, from about M1/M2=2.3, the latter competitor dominates and thus k increase with increasing mass ratio. In Figure 9, as expected the k’s from Iterative Scheme are generally larger than those from Standard SMRTA method. This difference increases with increasing mass ratio because the occurrence of N process increases when the acoustic-optical gap gets larger.Figure 9
Thermal conductivities calculated by Standard SMRTA method (dashed lines) and Iterative Scheme (solid lines) of1+1, 2+2, and 4+4 Si/Ge-based SLs as a function of M1/M2, the mass ratio of constituent atoms. The thin vertical line shows the Ge/Si mass ratio. Reprinted with permission from [50]. Copyrighted by the American Physical Society.


## 4. MD Simulation

### 4.1. Time Domain Normal Mode Analysis

The time domain normal mode analysis based on MD simulation was first proposed by Ladd et al. [66] and then modified by McGaughey and Kaviany [67]. From (10), a result of SMRTA, the relaxation time τλ can be obtained by
(29)τλ=∫0∞‍〈nλ′(t)nλ′(0)〉dt〈nλ′2(0)〉.
According to the analysis by Ladd et al. [66], the fluctuation nλ′ in (29) can be replaced by the total phonon occupation number nλ, which does not influence the calculation of thermal conductivity when considering that the ensemble-average heat current is zero. From lattice dynamics [34, 148], the occupation number nλ is proportional to the energy of single phonon mode λ described by the normal mode amplitude:
(30)Eλ=12(q˙λq˙λ*+ωλ2qλqλ*),
where qλ is the normal mode coordinate. Thus, (29) is transformed to
(31)τλ=∫0∞‍〈Eλ(t)Eλ(0)〉dt〈Eλ2(0)〉,
which is exactly what McGaughey and Kaviany [67] got, with the equivalent form 〈Eλ(t)Eλ(0)〉/〈Eλ2(0)〉=exp⁡(-t/τλ). Originally, Ladd et al. [66] only considered the potential energy and assumed Eλ~qλqλ*, which does not influence the result since the normal mode has the form [66, 149]
(32)qλ(t)=qλ,0exp⁡[i(ωλA+iηλ)t],
where qλ,0 is the vibration amplitude, a constant for a given mode λ, ωλA=ωλ+Δλ is the anharmonic frequency, and ηλ is linewidth. With the help of this equation, both Eλ~qλqλ* and Eλ~q˙λq˙λ* give the equivalent value of τλ=1/2ηλ.The calculation of normal mode coordinateqλ(t) is required to evaluate Eλ in (30) and further predict τλ in (31). From lattice dynamics [148],
(33)qλ(t)=∑α3‍∑bn‍∑lNc‍mbNcuαl,b(t)eb,αλ*exp⁡[ik·r0l]=∑α3‍∑bn‍eb,αλ*qαb(k,t),
where α indicates x, y, z directions, uαl,b(t) is α component of the displacement of the bth atom in lth unit cell from its equilibrium position, r0l is the equilibrium position of unit cell l, the star denotes complex conjugate, and qαb(k,t) denotes the contribution of the bth basis atom in α direction to the total normal mode with
(34)qαb(k,t)=∑lNc‍mbNcuαl,b(t)exp⁡[ik·r0l].
In (33), the time history of the atomic position displacement u(t) is extracted from MD simulation, and the eigenvector e is obtained from LD calculations.


### 4.2. Frequency Domain Normal Mode Analysis

Here, the frequency domain normal mode analysis is demonstrated by a simplified version; for detailed derivation, see [76, 77, 150]. Staring from (32), we have the spectral energy density (SED):
(35)Φλ(ω)=|q˙λ(ω)|2=|∫0+∞‍q˙λ(t)e-iωtdt|2=Cλ(ω-ωλA)2+ηλ2,
where Cλ=(ωλA2+Γλ2)qλ,02 is a constant related to λ. Physically, Φλ(ω) is the kinetic energy of single-phonon mode λ in the frequency domain, in contrast to (30) which is the energy in time domain. Equation (35) is actually a Lorentzian function with peak position ωk,νA and full width at half maximum 2ηλ. By fitting this SED function as Lorentzian form, the relaxation time τλ=1/2ηλ can be obtained.In some works, the total SED function for a given wave vectorΦ(k,ω)=∑ν‍Φλ(ω), which is the summation of the SEDs of phonons with the same k but from different phonon branches, is evaluated instead of that of each mode. Thomas et al. [77, 150] and Feng et al. [151] pointed out that the eigenvectors are unnecessary due to the orthogonality; thus,
(36)Φ(k,ω)=∑ν3n‍Φλ(ω)=14πτ0∑α,b3,n‍mbNc|∑lN‍∫0τ0‍u˙αl,b(t)exp⁡[ik·r0l-iωt]|2=∑α3‍∑bn‍|q˙αb(k,ω)|2,
where q˙αb(k,ω) is time derivative and Fourier Transform of (34). From (33) to (36), the eigenvector has been abandoned, and the mathematical proof of this is presented in [151].According to Ong et al. [152], the expression (36) is equivalent to the SED functions in [68–76]. In some of the works, the mass mb and unit cell number Nc in (34) are discarded for single-mass system, since the constants do not influence the fitting results of the Lorentzian function in (35), so that only atomic velocities are needed.


### 4.3. Discussion and Applications

Figures10 and 11 show two examples of time-domain NMA and frequency-domain NMA methods. Figure 10 presents the autocorrelation functions of total energy and potential energy of normal mode as functions of time of the TA mode of argon at 50 K. The oscillation of the potential energy indicates that the phonon frequency and the decay rate of total energy gives the relaxation time. Figure 11 shows the SED functions of empty CNT and water-filled CNT. From the fitting of these peaks as Lorentzian functions, the phonon frequencies and relaxation times are obtained. The linewidth broadening caused by the water filled is clear from Figure 11(b).Figure 10
Autocorrelation functions of potential and total energies of time-dependent normal modes of TA mode atk*=0.5 in [1,0,0] direction of argon at 50 K. Reprinted with permission from [67]. Copyrighted by the American Physical Society.SED functions (color online) of (a) empty CNT and (b) water-filled CNT alongkz=0 for frequencies below 5 THz at T = 298 K. Reprinted with permission from [77]. Copyrighted by the American Physical Society.
(a)
(b)The relaxation times predicted from MD simulation includes the effects of three-, four-, and higher-order phonon scattering processes; in contrast, ALD calculation only considers the lowest one. Thus, the ALD calculation may lose its accuracy when temperature increases, since the higher-order anharmonicity of lattice becomes greater for higher temperature due to thermal expansion. For instance, Turney et al. [102] compared the relaxation times of argon bulk predicted from the Standard SMRTA ALD calculation and the time-domain NMA at different temperatures. Figure 12 shows the inverse relaxation times of LA and TA phonon modes for argon at 20 K and 50 K. We note that, at 20 K, these two methods give reasonable agreement, whilst, at 50 K, the ALD calculation underpredicts the scattering rate by as much as 2 or more times.The inverse relaxation times (color online) of LA and TA mode of argon as functions of wave vector in [109] direction, predicted by Standard SMRTA method at 20 K (solid curve) and 50 K (dashed curve) and time domain NMA at 20 K (open triangle) and 50 K (solid triangle). Adapted with permission from [102]. Copyrighted by the American Physical Society.
(a)
(b)Compared to ALD calculation, MD simulation is a better tool for predicting the phonon properties of complex systems, such as the CNT filled with water and the graphene supported by substrate. So far, it is hard for ALD method to handle the extrinsic phonon scattering processes other than the Umklapp scattering without fitting parameters. However, in the MD simulation, the surrounding influence is reflected by the atomic vibrating trajectory of the studied system. Qiu and Ruan [127, 128, 140] studied the phonon transport in suspended and silicon dioxide supported SLG by frequency domain NMA with the results shown in Figures 13 and 14. We note that the flexural phonon modes (ZA and ZO) have much longer relaxation times than the other modes for suspended SLG, which qualitatively agree with the ALD calculation results discussed in Section 3.3.2. The MD result indicates that ZA mode contributes about 29% to the total k for suspended SLG, while TA and LA modes contribute 33% and 26%, respectively. Chen and Kumar [126] performed the same NMA method and obtained the similar results that ZA, TA, and LA modes contribute 23%, 21%, and 41%, respectively. The relaxation times of supported SLG are found generally shorter, by about 10 ps, than suspended SLG. This indicates that the SiO2 substrate provides strong phonon scattering by the interface and breaks down the reflection symmetry in suspended SLG. As a result, the percentage thermal conductivity contribution from ZA mode decreases about 10%, while those of TA and LA modes increase about 3% and 8%, respectively.Figure 13
The relaxation time (color online) of suspended (“s”) andSiO2 supported (“p”) single-layer graphene as a function of phonon frequency for different phonon branches at temperature 300 K. Reprinted with permission from [140]. Copyright 2012, AIP Publishing LLC.(a) Contributions (color online) and (b) the corresponding percentages of thermal conductivity from ZA, TA, LA, and TO modes in suspended (“s”) and supported (“p”) SLG at different temperatures. Reprinted with permission from [140].
(a)
(b)Due to the low computational complexity, the NMA methods have been applied to many cases. Time-domain NMA was used for Ar [66, 67, 102], Si [143, 153], Ge [154], and polyethylene [155, 156]; in the meanwhile, frequency domain NMA has been applied to Ar [150], Ge [151], MgO [76], CNT [77, 78, 157], supported CNT [152], suspended and supported graphene [140], and thermoelectric materials such as PbTe [141] and Bi2Te3 [142]. So far, only few works applied NMA to defective bulk, nanowires [153], and nanoribbons. As a representative application of frequency domain NMA to bulk material, the spectral phonon relaxation times of pristine PbTe bulk at different temperatures in different directions are given in Figure 15. The results reveal typical features of phonon relaxation time in bulk materials: (a) acoustic phonons generally have much higher relaxation times than optical phonons, (b) for acoustic modes, the relaxation times always decrease with increasing frequency except for the high-frequency ranges which often show opposite trend, such phenomenon is also found in other materials such as argon [67, 102], silicon [143], and germanium [151], (c) the value of α in frequency dependence relation τ~ω-α of the acoustic phonon often deviates from 2 and ranges from 0.5 to 4, (d) τ of optical mode has weak frequency dependence, and (e) increasing temperature typically shortens the phonon relaxation time and mean free path. The order of relaxation time amplitudes of PbTe bulk at 300 K obtained by frequency-domain NMA agrees well with those obtained from ALD calculation in previous sections [104, 106]. It is important to note that NMA methods do not distinguish between N scattering and U scattering but give a total scattering rate just as the method of ALD calculation based on Standard SMRTA. Figure 16(a) gives the contribution of each phonon mode to total thermal conductivity of PbTe bulk at 300 K and 600 K. The results show that optical modes only contribute 5% to the total k, different with 20% given by first principle ALD calculation. The discrepancy may come from the ignorance of higher-order phonon scattering in ALD calculation. Figure 16(b) gives the accumulated k’s as functions phonon MFP for Silicon at 300 K and PbTe at 300 K and 600 K. 80% of the total k of PbTe is contributed by the phonons with MFP below 50 nm, different from the value of 10 nm in ALD calculation [106]. This suggests that the relaxation times of low-frequency phonons predicted from ALD are longer than those from NMA, since both ALD and NMA results give reasonable total thermal conductivity. The MFPs of phonons of PbTe decrease roughly by a factor of 2 when temperature increases from 300 K to 600 K. It is found that the phonons with MFP below 10 nm contribute about 32% of k at 300 K while about 65% of k at 600 K.Spectral phonon relaxation time bulk PbTe (color online) at 300 K and 600 K. Reprinted with permission from [141]. Copyright (2011), with permission from Elsevier.
(a)
(b)
(c)
(d)Contribution (color online) of each phonon mode to total thermal conductivity of PbTe bulk at 300 K and 600 K. Reprinted with permission from [141]. Copyright (2011), with permission from Elsevier.
(a)
(b)The phonon properties of Bi2Te3 are studied by time-domain NMA [142]. The relaxation times and power law fitting of the low-frequency range are presented in Figure 17. The phonons with wavelength of 125 nm have relaxation time 16.9 ns, which indicate that those phonons do not experience obvious scattering when traveling for about 400 ps in Bi2Te3, consistent with experimental measurements [158]. The normalized accumulated thermal conductivity of Bi2Te3 as a function of phonon MFP is plotted in Figure 18. It is found that 90% and 50% of total thermal conductivity are contributed by the phonons with MFPs shorter than 10 nm and 3 nm, respectively. Also shown in Figure 18 is comparison between the results from ALD calculation and MD simulation. The two curves for Si agree well with each other, while a discrepancy is found for bulk PbTe. This discrepancy may come from the inaccuracy of the interatomic potential used in performing MD simulation. These results are useful for the nanodesign of Bi2Te3/PbTe/Si based thermoelectric materials in the future.(a) Phonon relaxation times (color online) of of Bi2Te3 along the Γ-Z direction computed using time-domain NMA. L denotes the number of cells along axis at 300 K. (b) Phonon relaxation times of low frequency acoustic phonons along Γ-Z and the power law fitting. Reprinted with permission from [142]. Copyright 2013 by ASME.
(a)
(b)Figure 18
The normalized accumulated thermal conductivity (color online) of several bulk materials at room temperature as a function of phonon MFP. (Si (MD) [143]; Si (ALD) [107]; PbTe (MD) [141]; PbTe (MD 600 K) [141]; PbTe (ALD) [106]; Bi2Te3 (MD) [142]).


## 4.1. Time Domain Normal Mode Analysis

The time domain normal mode analysis based on MD simulation was first proposed by Ladd et al. [66] and then modified by McGaughey and Kaviany [67]. From (10), a result of SMRTA, the relaxation time τλ can be obtained by
(29)τλ=∫0∞‍〈nλ′(t)nλ′(0)〉dt〈nλ′2(0)〉.
According to the analysis by Ladd et al. [66], the fluctuation nλ′ in (29) can be replaced by the total phonon occupation number nλ, which does not influence the calculation of thermal conductivity when considering that the ensemble-average heat current is zero. From lattice dynamics [34, 148], the occupation number nλ is proportional to the energy of single phonon mode λ described by the normal mode amplitude:
(30)Eλ=12(q˙λq˙λ*+ωλ2qλqλ*),
where qλ is the normal mode coordinate. Thus, (29) is transformed to
(31)τλ=∫0∞‍〈Eλ(t)Eλ(0)〉dt〈Eλ2(0)〉,
which is exactly what McGaughey and Kaviany [67] got, with the equivalent form 〈Eλ(t)Eλ(0)〉/〈Eλ2(0)〉=exp⁡(-t/τλ). Originally, Ladd et al. [66] only considered the potential energy and assumed Eλ~qλqλ*, which does not influence the result since the normal mode has the form [66, 149]
(32)qλ(t)=qλ,0exp⁡[i(ωλA+iηλ)t],
where qλ,0 is the vibration amplitude, a constant for a given mode λ, ωλA=ωλ+Δλ is the anharmonic frequency, and ηλ is linewidth. With the help of this equation, both Eλ~qλqλ* and Eλ~q˙λq˙λ* give the equivalent value of τλ=1/2ηλ.The calculation of normal mode coordinateqλ(t) is required to evaluate Eλ in (30) and further predict τλ in (31). From lattice dynamics [148],
(33)qλ(t)=∑α3‍∑bn‍∑lNc‍mbNcuαl,b(t)eb,αλ*exp⁡[ik·r0l]=∑α3‍∑bn‍eb,αλ*qαb(k,t),
where α indicates x, y, z directions, uαl,b(t) is α component of the displacement of the bth atom in lth unit cell from its equilibrium position, r0l is the equilibrium position of unit cell l, the star denotes complex conjugate, and qαb(k,t) denotes the contribution of the bth basis atom in α direction to the total normal mode with
(34)qαb(k,t)=∑lNc‍mbNcuαl,b(t)exp⁡[ik·r0l].
In (33), the time history of the atomic position displacement u(t) is extracted from MD simulation, and the eigenvector e is obtained from LD calculations.


## 4.2. Frequency Domain Normal Mode Analysis

Here, the frequency domain normal mode analysis is demonstrated by a simplified version; for detailed derivation, see [76, 77, 150]. Staring from (32), we have the spectral energy density (SED):
(35)Φλ(ω)=|q˙λ(ω)|2=|∫0+∞‍q˙λ(t)e-iωtdt|2=Cλ(ω-ωλA)2+ηλ2,
where Cλ=(ωλA2+Γλ2)qλ,02 is a constant related to λ. Physically, Φλ(ω) is the kinetic energy of single-phonon mode λ in the frequency domain, in contrast to (30) which is the energy in time domain. Equation (35) is actually a Lorentzian function with peak position ωk,νA and full width at half maximum 2ηλ. By fitting this SED function as Lorentzian form, the relaxation time τλ=1/2ηλ can be obtained.In some works, the total SED function for a given wave vectorΦ(k,ω)=∑ν‍Φλ(ω), which is the summation of the SEDs of phonons with the same k but from different phonon branches, is evaluated instead of that of each mode. Thomas et al. [77, 150] and Feng et al. [151] pointed out that the eigenvectors are unnecessary due to the orthogonality; thus,
(36)Φ(k,ω)=∑ν3n‍Φλ(ω)=14πτ0∑α,b3,n‍mbNc|∑lN‍∫0τ0‍u˙αl,b(t)exp⁡[ik·r0l-iωt]|2=∑α3‍∑bn‍|q˙αb(k,ω)|2,
where q˙αb(k,ω) is time derivative and Fourier Transform of (34). From (33) to (36), the eigenvector has been abandoned, and the mathematical proof of this is presented in [151].According to Ong et al. [152], the expression (36) is equivalent to the SED functions in [68–76]. In some of the works, the mass mb and unit cell number Nc in (34) are discarded for single-mass system, since the constants do not influence the fitting results of the Lorentzian function in (35), so that only atomic velocities are needed.


## 4.3. Discussion and Applications

Figures10 and 11 show two examples of time-domain NMA and frequency-domain NMA methods. Figure 10 presents the autocorrelation functions of total energy and potential energy of normal mode as functions of time of the TA mode of argon at 50 K. The oscillation of the potential energy indicates that the phonon frequency and the decay rate of total energy gives the relaxation time. Figure 11 shows the SED functions of empty CNT and water-filled CNT. From the fitting of these peaks as Lorentzian functions, the phonon frequencies and relaxation times are obtained. The linewidth broadening caused by the water filled is clear from Figure 11(b).Figure 10
Autocorrelation functions of potential and total energies of time-dependent normal modes of TA mode atk*=0.5 in [1,0,0] direction of argon at 50 K. Reprinted with permission from [67]. Copyrighted by the American Physical Society.SED functions (color online) of (a) empty CNT and (b) water-filled CNT alongkz=0 for frequencies below 5 THz at T = 298 K. Reprinted with permission from [77]. Copyrighted by the American Physical Society.
(a)
(b)The relaxation times predicted from MD simulation includes the effects of three-, four-, and higher-order phonon scattering processes; in contrast, ALD calculation only considers the lowest one. Thus, the ALD calculation may lose its accuracy when temperature increases, since the higher-order anharmonicity of lattice becomes greater for higher temperature due to thermal expansion. For instance, Turney et al. [102] compared the relaxation times of argon bulk predicted from the Standard SMRTA ALD calculation and the time-domain NMA at different temperatures. Figure 12 shows the inverse relaxation times of LA and TA phonon modes for argon at 20 K and 50 K. We note that, at 20 K, these two methods give reasonable agreement, whilst, at 50 K, the ALD calculation underpredicts the scattering rate by as much as 2 or more times.The inverse relaxation times (color online) of LA and TA mode of argon as functions of wave vector in [109] direction, predicted by Standard SMRTA method at 20 K (solid curve) and 50 K (dashed curve) and time domain NMA at 20 K (open triangle) and 50 K (solid triangle). Adapted with permission from [102]. Copyrighted by the American Physical Society.
(a)
(b)Compared to ALD calculation, MD simulation is a better tool for predicting the phonon properties of complex systems, such as the CNT filled with water and the graphene supported by substrate. So far, it is hard for ALD method to handle the extrinsic phonon scattering processes other than the Umklapp scattering without fitting parameters. However, in the MD simulation, the surrounding influence is reflected by the atomic vibrating trajectory of the studied system. Qiu and Ruan [127, 128, 140] studied the phonon transport in suspended and silicon dioxide supported SLG by frequency domain NMA with the results shown in Figures 13 and 14. We note that the flexural phonon modes (ZA and ZO) have much longer relaxation times than the other modes for suspended SLG, which qualitatively agree with the ALD calculation results discussed in Section 3.3.2. The MD result indicates that ZA mode contributes about 29% to the total k for suspended SLG, while TA and LA modes contribute 33% and 26%, respectively. Chen and Kumar [126] performed the same NMA method and obtained the similar results that ZA, TA, and LA modes contribute 23%, 21%, and 41%, respectively. The relaxation times of supported SLG are found generally shorter, by about 10 ps, than suspended SLG. This indicates that the SiO2 substrate provides strong phonon scattering by the interface and breaks down the reflection symmetry in suspended SLG. As a result, the percentage thermal conductivity contribution from ZA mode decreases about 10%, while those of TA and LA modes increase about 3% and 8%, respectively.Figure 13
The relaxation time (color online) of suspended (“s”) andSiO2 supported (“p”) single-layer graphene as a function of phonon frequency for different phonon branches at temperature 300 K. Reprinted with permission from [140]. Copyright 2012, AIP Publishing LLC.(a) Contributions (color online) and (b) the corresponding percentages of thermal conductivity from ZA, TA, LA, and TO modes in suspended (“s”) and supported (“p”) SLG at different temperatures. Reprinted with permission from [140].
(a)
(b)Due to the low computational complexity, the NMA methods have been applied to many cases. Time-domain NMA was used for Ar [66, 67, 102], Si [143, 153], Ge [154], and polyethylene [155, 156]; in the meanwhile, frequency domain NMA has been applied to Ar [150], Ge [151], MgO [76], CNT [77, 78, 157], supported CNT [152], suspended and supported graphene [140], and thermoelectric materials such as PbTe [141] and Bi2Te3 [142]. So far, only few works applied NMA to defective bulk, nanowires [153], and nanoribbons. As a representative application of frequency domain NMA to bulk material, the spectral phonon relaxation times of pristine PbTe bulk at different temperatures in different directions are given in Figure 15. The results reveal typical features of phonon relaxation time in bulk materials: (a) acoustic phonons generally have much higher relaxation times than optical phonons, (b) for acoustic modes, the relaxation times always decrease with increasing frequency except for the high-frequency ranges which often show opposite trend, such phenomenon is also found in other materials such as argon [67, 102], silicon [143], and germanium [151], (c) the value of α in frequency dependence relation τ~ω-α of the acoustic phonon often deviates from 2 and ranges from 0.5 to 4, (d) τ of optical mode has weak frequency dependence, and (e) increasing temperature typically shortens the phonon relaxation time and mean free path. The order of relaxation time amplitudes of PbTe bulk at 300 K obtained by frequency-domain NMA agrees well with those obtained from ALD calculation in previous sections [104, 106]. It is important to note that NMA methods do not distinguish between N scattering and U scattering but give a total scattering rate just as the method of ALD calculation based on Standard SMRTA. Figure 16(a) gives the contribution of each phonon mode to total thermal conductivity of PbTe bulk at 300 K and 600 K. The results show that optical modes only contribute 5% to the total k, different with 20% given by first principle ALD calculation. The discrepancy may come from the ignorance of higher-order phonon scattering in ALD calculation. Figure 16(b) gives the accumulated k’s as functions phonon MFP for Silicon at 300 K and PbTe at 300 K and 600 K. 80% of the total k of PbTe is contributed by the phonons with MFP below 50 nm, different from the value of 10 nm in ALD calculation [106]. This suggests that the relaxation times of low-frequency phonons predicted from ALD are longer than those from NMA, since both ALD and NMA results give reasonable total thermal conductivity. The MFPs of phonons of PbTe decrease roughly by a factor of 2 when temperature increases from 300 K to 600 K. It is found that the phonons with MFP below 10 nm contribute about 32% of k at 300 K while about 65% of k at 600 K.Spectral phonon relaxation time bulk PbTe (color online) at 300 K and 600 K. Reprinted with permission from [141]. Copyright (2011), with permission from Elsevier.
(a)
(b)
(c)
(d)Contribution (color online) of each phonon mode to total thermal conductivity of PbTe bulk at 300 K and 600 K. Reprinted with permission from [141]. Copyright (2011), with permission from Elsevier.
(a)
(b)The phonon properties of Bi2Te3 are studied by time-domain NMA [142]. The relaxation times and power law fitting of the low-frequency range are presented in Figure 17. The phonons with wavelength of 125 nm have relaxation time 16.9 ns, which indicate that those phonons do not experience obvious scattering when traveling for about 400 ps in Bi2Te3, consistent with experimental measurements [158]. The normalized accumulated thermal conductivity of Bi2Te3 as a function of phonon MFP is plotted in Figure 18. It is found that 90% and 50% of total thermal conductivity are contributed by the phonons with MFPs shorter than 10 nm and 3 nm, respectively. Also shown in Figure 18 is comparison between the results from ALD calculation and MD simulation. The two curves for Si agree well with each other, while a discrepancy is found for bulk PbTe. This discrepancy may come from the inaccuracy of the interatomic potential used in performing MD simulation. These results are useful for the nanodesign of Bi2Te3/PbTe/Si based thermoelectric materials in the future.(a) Phonon relaxation times (color online) of of Bi2Te3 along the Γ-Z direction computed using time-domain NMA. L denotes the number of cells along axis at 300 K. (b) Phonon relaxation times of low frequency acoustic phonons along Γ-Z and the power law fitting. Reprinted with permission from [142]. Copyright 2013 by ASME.
(a)
(b)Figure 18
The normalized accumulated thermal conductivity (color online) of several bulk materials at room temperature as a function of phonon MFP. (Si (MD) [143]; Si (ALD) [107]; PbTe (MD) [141]; PbTe (MD 600 K) [141]; PbTe (ALD) [106]; Bi2Te3 (MD) [142]).


## 5. Summary

The three methods, anharmonic lattice dynamics based on Standard SMRTA, iterative anharmonic lattice dynamics, and normal mode analysis, can all predict thermal conductivity by calculating the velocities, relaxation times, and specific heats of all phonon modes. The applications are listed in Table2, and the features of these methods are compared and listed in Table 3.Table 2
Applications of the numerical methods in predicting spectral phonon properties and thermal conductivity.
Materials
Methods*
FP†
Reference
Year
Ar
1 and 3
[66]
1986
Ar
2
[46]
1995
Ar and Kr
2
[45]
1996
Ar
3
[67]
2004
Ar
1 and 3
[102]
2009
Ar, Si thin films
1
[92]
2010
C, Si, and Ge
1
√
[43]
1995
C, isotope-doped C, Si, and Ge
1 and 2
√
[51]
2009
C (Pure and natural) (extreme pressure)
1 and 2
[62]
2012
C nanowire
1 and 2
√
[64]
2012
Si (isotope-doped)
1
[136]
1999
Si (isotope-doped)
1
[91]
2001
SiC
2
[48]
2002
Si and Ge
1
√
[44]
2003
Si
1 and 2
[60]
2005
Si and Ge
2
√
[61]
2007
Si
3
[143]
2008
Si (isotope-doped)
2
[47]
2009
Si
1
√
[107]
2011
Si
1
[99]
2012
Si Nanowire
3
[153]
2009
Si, Ge
1
√
[110]
2010
Si/Ge,M1/M2 SLs
1 and 2
[59]
2004
Si/Ge SLs
1
√
[108]
2011
Si/Ge SLs
1
[101]
2013
Si/Ge SLs
1
√
[103]
2013
SiGe alloys with embedded nanoparticles
2
√
[65]
2011
SiGe alloys
1
√
[109]
2011
Si, Ge, and Si0.5Ge0.5
1
[100]
2012
Si/Ge, GaAs/AlAs, andM1/M2 SLs
1 and 2
√
[50]
2008
Ge
3
[154]
2010
Ge
4
[151]
2013
Semiconductors (Groups IV, III–V, and II–VI)
1
[117]
2008
Graphene
1 and 2
[54]
2010
Graphene and graphite
1 and 2
[55]
2011
Graphene (supported and suspended)
4
[140]
2012
Graphene (free-standing and strained)
1
√
[105]
2012
Graphene and graphite
1
√
[111]
2013
CNT to graphene (diameter dependence)
1 and 2
[53]
2010
CNT
4
[157]
2006
CNT
1 and 2
[52]
2009
CNT (empty and water-filled)
4
[77]
2010
CNT (on amorphous silica)
4
[152]
2011
BN (pristine and isotope-doped)
1 and 2
[56]
2011
BN (multilayer and nanotubes) (pristine and  isotope-doped)
2
[58]
2012
Mg2SixSn1-x alloys (bulk and nanowire)
2
[63]
2012
Compound semiconductors (Si, Ge, GaAs, Al-V, Ga-V, In-V, SiC, AlN, etc.)
1 and 2
√
[112]
2013
Ionic solids (MgO,UO2, and SrTiO3)
2
[49]
2011
GaN (GaAs, GaSb, and GaP) (pristine and isotope-doped)
1 and 2
√
[57]
2012
PbTe
1
√
[104]
2012
PbTe
4
[141]
2011
PbTe, PbSe, andPbTe1-xSex
1
√
[106]
2012
Bi
2
Te
3
4
[142]
2013
Heusler
3
[120]
2011
MgO
4
√
[76]
2009
Polyethylene
3
[155]
2009
Polyethylene
3
[156]
2009
GaAs
1
√
[118]
2013
GaAs/AlAs SLs
1
√
[119]
2012
*Methods: 1: Standard SMRTA scheme; 2: Iterative Scheme; 3: time-domain NMA; 4: frequency-domain NMA.
†FP: first principle.Table 3
Comparison of different methods for predicting spectral phonon relaxation time and thermal conductivity.
Methods
Analytical model
ALD calculation
MD simulation
Standard SMRTA
Iterative Scheme
Time NMA
Frequency NMA
Equations
Table1
Equations (14), (15), and (16)
Equations (14), (15), (16), (17), and (18)
Equations (30), (31), and (33)
Equations (33) and (35)
Equations (33) and (36)
Characteristics
Lots of approximations, need fitting parameters
N and U processes are independent thermal resistant sources
N process does not provide thermal resistant itself
Eigenvectors needed
Eigenvectors not needed
Need 2nd- and 3rd-order IFCs
Need interatomic potential (orab initioMD)
Suitable for
Long wavelength, Debye model
Low temperature higher-order anharmonicity not large
Higher than Debye temperature
Temperature not too low, quantum effect negligible
Accuracy
Low
Medium
Higher
Higher
Computational complexity
High
Higher
Low
Applications so far
Some thermal conductivity analysis and prediction
Pure and isotope-doped bulk, alloy superlattice, nanostructures
Pure lattice Materials with surrounding influences
Further research
Temperature dependent IFCs, 4th- and higher-order phonon scattering
Accurate interatomic potential, large domain first    principle MD, defects, boundariesAll the three methods are based on phonon Boltzmann Transport Equation and relaxation time approximation. To obtain the spectral phonon relaxation time, the first two methods calculate three-phonon scattering rates from anharmonic interatomic force constants, while the last method calculate the linewidth of spectral energy in frequency domain or the decay rate of spectral energy in time domain from molecular dynamics. Since the first two methods ignore the 4th- and higher-order phonon scattering processes, they are only valid at low temperature. The first two methods differ with each other at solving the phonon BTE: the first method assumes single mode RTA, while the second one solves the linearized BTE iteratively instead. As a result, the first method treatsN scattering and U scattering as two independent processes that provide thermal resistance individually. However it is well known that the N scattering only contribute to thermal resistance by influencing the U scattering rate. The Iterative ALD remedies this error by recording all the phonon scattering processes step by step and evaluates the U scattering rates in the end. Compared to Green-Kubo MD (GK-MD) and Nonequilibrium MD (NEMD), these three methods give deeper insight into the thermal conductivity: the spectral phonon velocity, relaxation time, and mean free path, and the contribution of each phonon mode to thermal conductivity, which can guide the nanodesign. For accuracy and capability, theab initio ALD calculations are better than GK-MD and NEMD, since calculatingab initio 3rd-order IFCs is much easier than implementingab initio MD. The limitations of the normal mode analysis are as follows: (1) it cannot distinguish U and N processes and (2) it is of classical nature so it cannot accurately capture the quantum distribution function (Bose-Einstein distribution) for high Debye-temperature materials at relatively low temperatures (such as graphene and CNT at room temperature). The disadvantage of these three methods is the much computational cost. Compared to analytical models, these methods do not rely on adjustable fitting parameters and thus give more reliable and accurate predictions.These numerical methods have been applied to numerous materials and structures and revealed lots of physical nature that has never been reached before. The acoustic phonons are verified to have the ~ω-α frequency dependence which agrees with earlier analytical models, while the facts that the value of α varies from 0 to 4 at low frequency and that the frequency dependence becomes weak and abnormal at high frequency were not observed clearly before. The optical modes are found to carry very little heat but contribute much to the scattering of acoustic phonons and thus are essential to thermal transport. In layer-/tube-structured materials, the strict selection rule of phonon scattering because reflection symmetry severely blocks the scattering of flexural acoustic phonons and thus causes extremely high relaxation time and then high thermal conductivity. In short-period superlattice, the large gaps between acoustic and optical phonon branches make the scattering rarely happen and thus lead to high thermal conductivity, even higher than its corresponding pure materials. These methods are also applied to defected and alloy materials using virtual crystal approach. Despite these applications, further work is still needed to predict spectral phonon properties more accurately and efficiently, such as considering the temperature-dependent IFCs and higher-order anharmonicities in ALD calculations, implementing large domainab initio molecular dynamics for normal mode analysis.



---

*Source: 206370-2014-03-31.xml*