Size Effects in Thermal Conduction by Phonons
SSize Effects in Thermal Conduction by Phonons
Philip B. Allen ∗ Physics and Astronomy Department, Stony Brook University, Stony Brook, NY 11794-3800, USA (Dated: October 11, 2018)Heat transport in nanoscale systems is both hard to measure microscopically, and hard to inter-pret. Ballistic and diffusive heat flow coexist, adding confusion. This paper looks at a very simplecase: a nanoscale crystal repeated periodically. This is a popular model for simulation of bulkheat transport using classical molecular dynamics (MD), and is related to transient thermal gratingexperiments. Nanoscale effects are seen in perhaps their simplest form. The model is solved byan extension of standard quasiparticle gas theory of bulk solids. Both structure and heat flow areconstrained by periodic boundary conditions. Diffusive transport is fully included, while ballistictransport by phonons of long mean free path is diminished in a specific way. Heat current J ( x ) andtemperature gradient ∇ T ( x (cid:48) ) have a non-local relationship, via κ ( x − x (cid:48) ), over a distance | x − x (cid:48) | determined by phonon mean free paths. In MD modeling of bulk conductivity, finite computerresources limit system size. Long mean free paths, comparable to the scale of heating and cooling,cause undesired finite-size effects that have to be removed by extrapolation. The present modelallows this extrapolation to be quantified. Calculations based on the Peierls-Boltzmann equation,using a generalized Debye model, show that extrapolation involves fractional powers of 1 /L . It isalso argued that heating and cooling should be distributed sinusoidally ( ˙ e ∝ cos(2 πx/L )) to improveconvergence of numerics. PACS numbers: 66.70.-f, 63.22.-m, 65.80.-g
I. INTRODUCTION
The linear relation between heat current (cid:126)J and tem-perature gradient (cid:126) ∇ T is J α ( (cid:126)r ) = − (cid:90) d(cid:126)r (cid:48) κ α,β ( (cid:126)r, (cid:126)r (cid:48) ) ∇ β T ( (cid:126)r (cid:48) ) . (1)This non-local expression defines the general linear re-sponse for a time-independent (steady state) heat flow.This paper concerns cases where heat flows in a single( x ) direction in response to a temperature gradient inthe same direction. Then the spatial variation of (cid:126)J , (cid:126) ∇ T ,and κ can be simplified by averaging over y and z di-rections. The conductivity κ α,β ( (cid:126)r, (cid:126)r (cid:48) ) becomes the scalar κ ( x, x (cid:48) ), J x ( x ) = − (cid:90) dx (cid:48) κ ( x, x (cid:48) ) ∇ x T ( x (cid:48) ) . (2)When the distance scale of variation of the temperaturegradient is longer than the carrier mean-free path Λ, a lo-cal approximation (the usual Fourier law J x = − κ ∇ x T )works. This is the usual situation in a macroscopic mea-surement on a spatially homogeneous sample; the tem-perature gradient typically has negligible spatial varia-tion on the scale of Λ. In this limit, κ ( x, x (cid:48) ) = κ ( x − x (cid:48) ),and κ = (cid:82) dxκ ( x ). In the homogeneous case, assuminga long sample, Fourier variables are appropriate. The lin-ear relation is then J x ( k ) = − κ ( k ) ∇ x T ( k ). I will use thesame symbol ( ∇ x T for example) to indicate both coordi-nate space and reciprocal space representations of func-tions, the coordinate x or k being explicitly shown. Fora homogeneous sample and constant temperature gradi-ent, the experiment is described by the k = 0 Fourier component, and κ = lim k → κ ( k ).Nanoscale systems have boundaries that addcomplexity . There is an argument that saysthe generalized Fourier law Eq.(1) does not apply inall cases. This paper avoids such issues, and addressesa case where the complexity is minimized, namely, ahomogeneous nanoscale system with periodic boundaryconditions. This geometry is experimentally realizedin transient thermal grating experiments . It is alsoa preferred geometry for simulation of heat transportby classical molecular dynamics (MD) using the “directmethod” . Heat is introduced locally, extractedlocally some distance away, and temperature is moni-tored in between. The aim of the modeling is usuallyto extract the conductivity of a macroscopic sample.Nanoscales automatically enter, because the computercannot process atomic information on a macroscopicscale.A very nice example was given by Zhou et al. , whocarefully analyze computational accuracy, using the semi-conductor GaN as an example. The large thermal con-ductivity (several hundred W/mK at room temperature)indicates that the dominant acoustic phonons have meanfree paths Λ Q of order hundreds of interatomic spac-ings a . Since acoustic-phonon mean-free paths increaserapidly as wave-vector (cid:126)Q →
0, significant heat is carriedby phonons with much longer mean free paths, thousandsin units of a . The simulations were done for periodiccells up to length L ≈ a . Probably another factorof 10 in length would be required to fully converge theanswers. At a series of lengths L < a , conductivity κ eff ( L ) was extracted as the ratio − J x / ∇ x T , where thetemperature gradient was determined at x = ± L/
4, halfway between the heat source and sink. Extrapolation of a r X i v : . [ c ond - m a t . m e s - h a ll ] J un κ eff ( L ) to L → ∞ (assuming κ eff ( L ) − κ bulk ∝ /L ) in-dicated that the value of κ bulk was typically twice biggerthan the maximum achieved value κ eff ( L ≈ a ). Zhou et al. discover evidence of a breakdown in this methodof extrapolation. The breakdown was analyzed by Sel-lan et al. ; the alternate analysis given here uses similarDebye-model simplifications, but different boundary as-sumptions. Both offer ways of improving extrapolationneeded in MD computation of κ .The simplicity of the periodic boundary means thatphonon gas theory can be easily adapted to this nanoscalesituation. Here I offer such an analysis. The result sup-ports the idea of Sellan et al. that κ eff ( L ) − κ bulk ∝ / √ L should give a better extrapolation. The analysisalso points to a better algorithm for adding and removingheat. II. SEGMENTED PERIODIC SLAB MODEL
Fig. 1 is a cartoon system (“simulation cell”) of N “slabs,” each of width d , corresponding to a total length L = N d . The system is repeated periodically, as some-times used in classical MD simulations. The primaryvariables are the imposed rate of heating per unit vol-ume, ˙ e ( (cid:96) ), of the (cid:96) ’th slab, and the “measured” slab tem-perature T ( (cid:96) ) (mean kinetic energy of the atoms in theslab, divided by 3 k B / J x ( (cid:96) + ) and the temperature gradient ∇ x T ( (cid:96) + ); bothare defined at the junction of slabs (cid:96) and (cid:96) + 1. The sec-ondary variables are related to the primary variables bythe two fundamental slab-ring equations, ∇ x T ( (cid:96) + 12 ) = [ T ( (cid:96) + 1) − T ( (cid:96) )] /d (3)˙ e ( (cid:96) ) = [ J x ( (cid:96) + 12 ) − J x ( (cid:96) −
12 )] /d (4)Equation 4 is energy conservation. To allow a steadystate, it is required that (cid:80) (cid:96) ˙ e ( (cid:96) ) = 0, meaning that nonet heating occurs.The relation between heat flux and temperature gra-dient, in linear approximation, is J x ( (cid:96) + 12 ) = − N − (cid:88) m =0 κ ( (cid:96) − m ) ∇ x T ( m + 12 ) (5)This is the analog of Eq.(2). The position variable x has been discretized into the slab index (cid:96) . The thermalconductivity is similarly discretized. In terms of thesediscrete and periodic variables, the form of Eq.(5) is re-quired by the linear approximation. The matrix κ ( (cid:96), m )is rigorously defined. It is independent of the heating˙ e ( (cid:96) ). If the slabs are identical, then κ ( (cid:96), m ) retains the
0 1 2 3 4 5 6 7 8
FIG. 1. Two schematics of a periodic supercell used for MDmodeling of thermal conductivity. The cell is divided into N slabs (here N = 16) in the direction x of heat flow. Itis periodically repeated in all three directions. Heat is ran-domly inserted as impulses on randomly chosen atoms in slab0. Equal random extraction of heat occurs in slab N/
2. Thusthe heat current density J is controlled. Temperature is mea-sured in slabs 1 , , . . . , N/ −
1, and also in the similar slabs − , − , . . . . The temperature gradient dT /dx is minimumin slab N/ N/
4. From these gradients, the value of κ eff ( N ) = J/ | dT /dx | min is evaluated. N -fold translational symmetry of the ring. Therefore κ ( (cid:96) + n, m + n ) = κ ( (cid:96), m ) = κ ( (cid:96) − m, ≡ κ ( (cid:96) − m ).All variables are N -fold periodic in (cid:96) , including the non-local conductivity κ ( (cid:96) ) = κ ( (cid:96) + N ).Conjugate to the N positions on the ring are N Fouriervectors q = (2 π/L ) n q , with L = N d and n q defined mod-ulo N . I will always use k for the Fourier transform ofthe continuous variable x , and q for the discrete case.This notation then clarifies that κ ( k ) refers to thermalconductivity of a large homogeneous system, and κ ( q ) tothe small but periodically repeated supercell.Just as the integers (cid:96) are defined modulo N , similarlythe integers n q have the same periodicity, n q + N = n q .The Fourier representation of the discretized variables isparticularly simple and convenient:˙ e ( (cid:96) ) = (cid:88) q e iqd(cid:96) ˙ e ( q ) T ( (cid:96) ) = (cid:88) q e iqd(cid:96) T ( q ) J x ( (cid:96) + 12 ) = (cid:88) q e iqd ( (cid:96) + ) J x ( q ) ∇ x T ( (cid:96) + 12 ) = (cid:88) q e iqd ( (cid:96) + ) ∇ x T ( q ) κ ( (cid:96) ) = (cid:88) q e iqd(cid:96) κ ( q ) , (6)where sums go over the N distinct values of q . The re-verse transforms are˙ e ( q ) = 1 N (cid:88) (cid:96) e − iqd(cid:96) ˙ e ( (cid:96) ) J x ( q ) = 1 N (cid:88) (cid:96) e − iqd ( (cid:96) + ) J x ( (cid:96) + 12 ) , (7)and similar for T , ∇ x T , and κ . Sums go over the N distinct values of (cid:96) . The reciprocal space version of Eq.(5)is J x ( q ) = − κ ( q ) ∇ x T ( q ) (8) III. PHONON GAS THEORY
This section assumes a macroscopic homogeneous(“continuum”) solid. The next section translates thisto a segmented supercell of periodic (“discrete”) slabs.The short-hand Q , in both continuum and discrete mod-els, denotes ( (cid:126)Q, s ), the (three-dimensional) wavevector (cid:126)Q and the branch index s of the phonons. Spatial vari-ation is driven by external heating and is assumed tooccur only in the x direction. Thermal properties aredescribed by N Q ( x ), the mean occupation of mode Q at(one-dimensional) position x . The heat current is J x ( x ) = 1Ω (cid:88) Q (cid:126) ω Q v Qx N Q ( x ) , (9)The Peierls-Boltzmann equation (PBE) describesthe dynamics of N Q ( x ). Anharmonic phonon events drive N Q to the local thermal equilibrium Bose-Einstein distri-bution n Q = 1 / (exp( (cid:126) ω Q /k B T ( x )) − N Q ( x, t ) change by phonon drift with velocity (cid:126)v Q = ∂ω Q /∂ (cid:126)Q . In steady state, N Q is stationary, ∂N Q ∂t = 0 = (cid:18) ∂N Q ∂t (cid:19) drift + (cid:18) ∂N Q ∂t (cid:19) collisions . (10)Writing N Q as the sum of a local equilibrium n Q ( T ( x ))(unaffected by collisions) and a deviation Φ Q ( x ),the drift term of the PBE has two contributions, − ( ∂n Q /∂T ) v Qx ∇ x T and − v Qx ∇ x Φ Q . The second ofthese vanishes when the temperature gradient is con-stant, but becomes important when there is spatial in-homogeneity. The collision term, after linearization,has the form − (cid:80) Q (cid:48) C QQ (cid:48) Φ Q (cid:48) , where C QQ (cid:48) is a com-plicated collision operator. To first approximation (the“relaxation time approximation,” RTA) this can be writ-ten as − Φ Q /τ Q , where 1 /τ Q is the “single mode relax-ation rate,” meaning the thermalization rate (or life-timebroadening) of mode Q that appears if only that onemode is out of equilibrium. The PBE then takes theform [ v Qx ∇ x + 1 /τ Q ]Φ Q ( x ) = − v Qx ∂n Q ∂T ∇ x T (11) In the continuum picture, the direct andreciprocal-space representations are related byΦ Q ( x ) = (1 / π ) (cid:82) dk exp( ikx )Φ Q ( k ) and Φ Q ( k ) = (cid:82) dx exp( − ikx )Φ Q ( x ). The PBE becomes( ikv Qx + 1 /τ Q )Φ Q ( k ) = − v Qx ∂n Q ∂T ∇ x T ( k ) . (12)The thermal conductivity κ ( k ) in continuous k spaceobeys the reciprocal-space version of Eq.(2), J x ( k ) = − κ ( k ) ∇ x T ( k ), and has the form κ ( k ) = 1Ω (cid:88) Q (cid:126) ω Q v Qx ( ∂n Q /∂T )1 /τ Q + ikv Qx (13)This can be Fourier-transformed back to x -space: κ ( x − x (cid:48) ) = 1Ω v Qx > (cid:88) Q (cid:126) ω Q v Qx ∂n Q ∂T exp (cid:20) − | x − x (cid:48) | Λ Qx (cid:21) , (14)where the mean free path is Λ Qx = v Qx τ Q . This uses theproperty that v Qx and Λ Qx change sign when (cid:126)Q changessign.The result Eq.(14) shows explicitly that heat trans-port is influenced non-locally by temperature variations.The velocity | v Qx | is bounded, but the mean free pathΛ Qx = v Qx τ Q is not. Relaxation times of long wave-length acoustic phonons diverge as | (cid:126)Q | →
0. In the over-simplified (“gray”) model where every phonon Q has thesame relaxation rate 1 /τ Q = 1 /τ , the Q -integrated κ ( x )decays rapidly. But because the spatial decay rate inEq.(14) involves 1 / Λ | cos θ Q | , no simple exponential for-mula works. The decay is more rapid than the simple ex-ponential exp( −| x | / Λ ). However, in more realistic mod-els, the divergence of τ Q and Λ Qx at small Q causes, after Q -integration, a slower than exponential decay. Still, inEq.(14), even long-wavelength phonons are diffusive be-cause the sample is “macroscopic.” IV. GAS THEORY ON THE SEGMENTEDPERIODIC SLAB
Gas theory translates to the segmented supercell in aslightly awkward way. Temperature T ( (cid:96) ) is a slab prop-erty, but current J x ( (cid:96) + ) is a junction property. Gas the-ory has temperature as a primary variable, so T ( (cid:96) ) shouldcorrespond to T ( x ) averaged over the (cid:96) ’th slab. Evidently N Q ( (cid:96) ) is a slab property, not a junction property. Butsince the current in gas theory is (cid:80) Q (cid:126) ω Q v Qx N Q / Ω, weare forced to compromise and define J x ( (cid:96) + 12 ) = 1Ω (cid:88) Q (cid:126) ω Q v Qx N Q ( (cid:96) ) + N Q ( (cid:96) + 1)2 . (15)Transforming to the discrete slab Fourier representation,this becomes J x ( q ) = 1Ω (cid:88) Q (cid:126) ω Q v Qx cos( qd/ Q ( q ) . (16)Similarly, the temperature gradient, Eq.(3), in Fouriervariables, is ∇ x T ( q ) = 2 i sin( qd/ T ( q ) /d, (17)and the fundamental connection Eq.(4) between heat in-put and current is˙ e ( q ) = 2 i sin( qd/ J x ( q ) /d. (18)The remaining task is to translate Eq.(13) to slab lan-guage. The ingredient needing translation is ∇ x whichappears twice in Eq.(11). The translation of this equa-tion isΦ Q ( (cid:96) ) = θ ( v Qx ) v Qx τ Q /d × (cid:20) ∂n Q ∂T ( T ( (cid:96) − − T ( (cid:96) )) + (Φ Q ( (cid:96) − − Φ Q ( (cid:96) )) (cid:21) + θ ( − v Qx ) v Qx τ Q /d × (cid:20) ∂n Q ∂T ( T ( (cid:96) ) − T ( (cid:96) + 1)) + (Φ Q ( (cid:96) ) − Φ Q ( (cid:96) + 1)) (cid:21) , (19)where θ ( x ) is the unit step function. The meaning is thatdrift entering slab (cid:96) from the left uses positive velocityphonons which carry information from the slab on theleft, while drift entering slab (cid:96) from the right uses nega-tive velocity phonons bringing information from the slabon the right.The discretized PBE, Eq.(19) can be solved in discretereciprocal space, givingΦ Q ( q ) = − v Qx τ Q ( ∂n Q /∂T ) ∇ x T ( q ) S ( q ) + 2 i sin( qd/ v Qx τ Q /d (20)where S ( q ) is exp( iqd/
2) if v Qx is positive andexp( − iqd/
2) if v Qx is negative. This gives the answerfor the discrete slab gas theory thermal conductivity, κ ( q ) = 1Ω (cid:88) Q (cid:126) ω Q ∂n Q ∂T v Qx τ Q cos( qd/ × (cid:20) S ( q ) + 2 i sin( qd/ v Qx τ Q /d (cid:21) (21)This agrees with the continuum formula Eq.(13) if thesmall q limit is taken, exp( iqd/ ≈ ≈ cos( qd/
2) and2 sin( qd/ ≈ qd . Since v Qx changes sign and S ( q ) be-comes S ( q ) ∗ when (cid:126)Q goes to − (cid:126)Q , the sum in Eq.(21)is real. Taking the real part of the factor in brackets, Eq.(21) becomes κ ( q ) = 1Ω (cid:88) Q (cid:126) ω Q ∂n Q ∂T v Qx τ Q cos ( qd/ F ( q, Λ Qx ) (22) F ( q, Λ) = (cid:34) ( qd/ (cid:40)(cid:18) Λ d (cid:19) + (cid:18) Λ d (cid:19) (cid:41)(cid:35) − (23)where Λ Qx = v Qx τ Q . Now consider what happensat the smallest q , namely q min = 2 π/L . The fac-tor cos ( q min d/
2) = cos ( π/N ) can be set to 1 andsin ( q min d/
2) can be set to ( q min d/ when the num-ber N of slabs is large. The correction factor in Eq.(23)becomes, in the large N case, F ( q min , Λ Qx ) ≈ (cid:34) (cid:18) π Λ Qx L (cid:19) (cid:35) − (24)The part of Eq.(23) linear in Λ /d is neglected comparedto the quadratic part since the correction is only impor-tant when Λ /d (cid:29)
1. Equations 22 and 24 show thatwhen a mean free path | Λ Qx | becomes comparable to L/ π , the contribution of that phonon to κ ( q min ) startsto be suppressed, the suppression becoming complete forphonons Q with 2 π Λ Qx (cid:29) L . The reason is, if a phonon’smean free path is as large as the period (slab ring circum-ference), that phonon is now carrying heat from hotterregions to random regions (after cycling around the ringmultiple times.) This is a different version of a well-known phenomenon in mesoscale heat transport, where κ in a sample of size L is diminished because long wave-length phonons travel ballistically a shorter distance L ,rather than diffusing the longer distance Λ that theywould exhibit in bulk . Computed behavior basedon Eqs. 22 and 23 will be shown in the next section,using a Debye model.The results developed above enable predictions of ∇ x T ( (cid:96) + ) and T ( (cid:96) ) for any given input ˙ e ( (cid:96) ). The ideais to use Eqs.(6, 8) to give ∇ x T ( (cid:96) + 12 ) = − (cid:88) q (cid:54) =0 e iqd ( (cid:96) + ) J x ( q ) /κ ( q )= − (cid:88) q (cid:54) =0 e iqd(cid:96) ˙ e ( q ) d (1 − e − iqd ) κ ( q ) , (25)where the second line follows from the first by usingEq.(18). This can be evaluated using Eq.(21) for κ ( q ).Using Eq.(17), the temperature T ( (cid:96) ) can also be evalu-ated, from T ( (cid:96) ) = T + (cid:88) q (cid:54) =0 e iqd(cid:96) ˙ e ( q ) d ( qd/ κ ( q ) , (26)where T is the average temperature. The q = 0 termof these sums is omitted because Eqs.(3,4) make it clearthat ∇ x T ( q = 0) = (cid:80) (cid:96) ∇ x T ( (cid:96) ) = 0, and similarly for˙ e ( q = 0).Two approximations have been made. One is the RTA.The other is discretization error. Temperature is a statis-tical variable, not definable except by averaging over a fi-nite volume. Therefore, discretization over a small width d should not cause noticeable error. However, a problemarises because gas theory has been forced to conform tothe slab/junction dichotomy of the discrete picture. Thiscauses the current in the (cid:96) + junction to be tied to thetemperature both of the two slabs (cid:96) + 1 and (cid:96) + 2 to theright, minus the temperature of both of the two slabs (cid:96) and (cid:96) − e ( (cid:96) ) is confined to a singlesite ( (cid:96) = 0 for example.) Then discretization introducessingular responses in the form of absurd oscillations inthe (unphysical) limit of very short mean free path. Thisproblem can be cured by distributing the heat input overthree slabs. Specifically, the cure is to use ˙ e ( (cid:96) ) = ˙ e/ (cid:96) = 0, and ˙ e/ (cid:96) = ±
1. This modificationhas more “realism” than a single-site input. Still, it issurprising that it is required in order for gas theory towork smoothly in a slab model. Within the homogeneousslab model, the discrete non-local conductivity (Eq. 5) isan exactly defined concept, as is its Fourier representa-tion κ ( q ). Equations (25,26) are exact connections withinlinear response, while Eq. 21 uses the PBE, plus a further(and not essential) simplification, the RTA. V. DEBYE MODEL
Full solution of the PBE, using accurate phonon prop-erties from density functional theory (DFT), is nowwidely available . Nevertheless, it is useful to have asimplified model as a standard to compare real calcula-tions against. The Debye model replaces the Brillouinzone by a sphere of radius Q D = (6 π n ) / , where n isthe number of atoms per unit volume. There are threebranches of phonons, approximated by ω Q = v | (cid:126)Q | withthe velocity v the same for each branch, and three or-thogonal directions (cid:126)v Q . The maximum frequency phonon( ω Q, max = ω D = vQ D ) occurs at the edge of the Debyesphere where | (cid:126)Q | = Q D . To establish a notation, thespecific heat can be written as C ( T ) = C ∞ (cid:90) ω D dωc ( ω ) (27)where c ( ω ) = D ( ω ) E ( ω ), and C ∞ = 3 nk B is the high T classical value of C ( T ). The factor D ( ω ) is the phonondensity of states, normalized to 1. In Debye approxi- mation, this is D D ( ω ) = 3 ω /ω D . The other factor is E ( ω ) = [( (cid:126) ω/ k B T ) / sinh( (cid:126) ω/ k B T )] , the Einstein for-mula for the specific heat (in units k B ) of one vibra-tional mode. This is replaced by the classical limit E = 1( C ( T ) = C ∞ ) when comparing with classical MD results.In the spirit of the Debye model, the mean free pathΛ Q = vτ Q can be modeled as Λ / Λ min ( T ) = ( Q D / | (cid:126)Q | ) p ,or equivalently, ( ω D /ω ) p . The exponent p depends ondetails of scattering. For point impurities such as iso-topic substitutions, p takes the Rayleigh value p = 4.For “Normal” (N, not Umklapp, or U) anharmonic scat-tering, Herring found p = 1 for transverse and p = 2for longitudinal acoustic modes at low ω . For anhar-monic U scattering, which is more relevant here, it isusually argued that p = 2, but both p = 3 and p = 4have some support from numerical calculations . Theminimum mean free path, Λ min is found at | (cid:126)Q | = Q D . Ithas a temperature-dependent value, scaling as 1 /T fromanharmonic scattering in the classical high T limit.First consider the continuum theory. The Debye ver-sion has a simple answer in the unphysical case of p = 0,where all phonons have the same mean free path, Λ min .From Eq.(13), the answer is κ D,p =0 ( k ) = κ ( T ) g ( k Λ min ) , (28) κ D,p =0 ( k → ≡ κ ( T ) = 13 C ( T ) v Λ min ( T ) . (29)The function g ( k Λ min ) = g ( u ) is defined as g ( u ) = 32 (cid:90) − d cos θ cos θ iu cos θ = 3( u − tan − u ) u . (30)The variable of integration, cos θ , is the cosine of theangle between the (3d) phonon wavevector (cid:126)Q and thedirection ˆ x of the applied temperature gradient. Thefunction g ( u ) = κ D,p =0 ( k ) /κ ( T ) is plotted in Fig. 2 forthe case where Λ min = 10 d . In the small k limit, thevalue is g ≈ − u / p > κ D,p ( k ) κ ∞ = (cid:90) ω D dωc ( ω ) (cid:16) ω D ω (cid:17) p g (cid:16) k Λ min (cid:16) ω D ω (cid:17) p (cid:17) , (31)where κ ∞ = C ∞ v Λ min ( T ) / p ≥
3, Eq.(31) diverges in the k → ω limit diverges logarithmically at p = 3,and more severely at higher p . In theory, the divergenceis cut off by finite sample size. This is rarely seen ex-perimentally, since p = 2 “N scattering,” and Akhieserdamping provide alternatives. Finite sample size in-troduces a complication. The boundaries destroy thehomogeneity of the theory. Even phonons with k sig-nificantly larger than k min have Λ > L , and carry heatballistically, unaware of the spatial variation of sampletemperature T ( x ). To first approximation, this gives alower limit ( ω (Λ = L )) below which the integral Eq.(31) q ( p /A) k ( q ) / k ( ) continuum, p=0discrete, p=0discrete, p=2 FIG. 2. Thermal conductivity versus wavevector, computedin Debye approximation for Λ min = 10 d , and plotted versuswavevector 0 ≤ q ≤ π/d . The top curve is the continuummodel, Eqs.(28,30) with constant mean free path Λ = Λ min .The middle curve is for the discrete model, Eq.(32) with con-stant mean free path p = 0. The bottom curve, from Eq.(35),uses a more realistic frequency-dependent mean free path,with p = 2. is cut off. Phonons with ω min < ω < ω (Λ = L ) give bal-listic currents, less than their bulk diffusive contribution,and outside the usual local version of the Fourier law. Inbulk crystals, this contribution can be important at low T , but not at higher T where such phonons are a verysmall minority and the ballistic component is negligible.The Debye model can also be used to find expressionsfor the non-local conductivity κ ( q ) of the periodic slabmodel, analogous to κ ( k ) of the continuum model. Notethat q and k (discrete and continuous wavevector) are thenotational clue indicating the model being solved. Firsttake the unphysical model with constant mean free path( p = 0). Using Eq.(21), the analog of Eq.(28) is κ D,p =0 ( q ) = κ ( T ) h (cid:48) ( q, Λ min ) = κ ( T ) h ( w ) (32)where h (cid:48) ( q, Λ min ) = h ( w ) is h ( w ) = 3 cos( qd/ (cid:20)(cid:90) d cos θ cos θe iqd/ + iw cos θ + (cid:90) − d cos θ cos θe − iqd/ + iw cos θ (cid:21) . (33)Here w = 2 sin( qd/ min /d is the discretized version of u = k Λ min that appears in Eq.(30). Note that, although h depends separately on q and on Λ min , the additional q -dependence beyond that contained in w plays only therole of a fixed parameter, while the w dependence ac-quires additional importance when the mean-free pathΛ Q acquires ω Q -dependence. Performing the d cos θ integral gives h ( w ) = 3 cos( qd/ w Re (cid:110) we iqd/ + ie iqd ln (cid:104) iwe − iq/d (cid:105)(cid:111) . (34)This reduces to Eq.(30) in the small q limit, underthe replacements q → k and w → u . The function h = κ D,p =0 ( q ) /κ ( T ) is also shown in Fig. 2. Up un-til q ≈ . π/d ( q ≈ π/ Λ min ) the discrete and contin-uum versions fall almost equally rapidly with q to < . q = ± π/d . This means,for example, that there is no response to input heating˙ e ( (cid:96) ) = ˙ e exp( ± iπ(cid:96) ) = ( − (cid:96) ˙ e . The reason is that adja-cent junctions have currents J ( (cid:96) + ) = ± ˙ e/
2. Therefore,the slab current (the average of the two adjacent junctioncurrents) is zero.For the more physical case of p >
0, the answer is κ D,p ( q ) κ ∞ = (cid:90) ω D dωc ( ω ) (cid:16) ω D ω (cid:17) p h (cid:16) q Λ min (cid:16) ω D ω (cid:17) p (cid:17) (35)The function κ D,p =2 ( q ) /κ ∞ is shown in the classical limit( c ( ω ) = D ( ω )) in Fig. 2 as the bottom curve. The con-ductivity falls more much rapidly with q , which meansincreased non-locality. This is not surprising. Spatialmemory extends much farther because of the longer meanfree paths. These results are translated back to coordi-nate space in Fig. 3. The p = 0 (constant Λ) resultsfall exponentially (as exp( − md/ Λ)) if the supercell size
N d exceeds Λ sufficiently. The p = 2 case behaves dif-ferently. At a distance md = 8Λ min = 80 d , the value of κ ( m ) is ≈ κ (0) / .
4, falling much more slowly than anexponential.
VI. NUMERICAL RESULTS
Fig. 4 illustrates the computed spatial temperaturevariation for the discrete slab model. The behaviorclosely resembles that found by Zhou et al. in theirclassical MD simulation of GaN. Therefore, I believe thatthe PBE, as extended here to discrete slabs, and mod-eled in RTA and in the Debye approximation, correctlycaptures the physics. The calculations of Fig. 4 usethe classical limit C ( T ) = C ∞ , with mean free pathΛ = Λ min ( ω D /ω ) , Λ min = 20 d , and various total celllengths L ranging from 40 d to 2560 d . Symmetry requires T ( N/ m ) = − T ( N/ − m ), and an inflection pointin T ( (cid:96) ) at (cid:96) = N/
4. At the largest N shown, the dis-tance N/ min , andthe answer for κ eff ( N ) is still 10% lower than the macro-scopic ( N → ∞ ) limit. The temperature profile is ac-curately linear, over the 60 slabs shown, for the threelargest lengths L , but the slopes are 10 to 30% higherthan the macroscopic limit.Zhou et al. invoke a Mattheissen’s rule justifica-tion for extrapolation to L → ∞ , namely the ideathat “boundary scattering” causes 1 /κ to behave like -80 -60 -40 -20 0 20 40 60 80slab number10 -5 -4 -3 -2 -1 k ( m ) / k N= 40, p=0N= 80, p=0N=160, p=0N= 40, p=2N= 80, p=2N=160, p=2
FIG. 3. The non-local thermal conductivity, κ ( m ) definedin Eq.(5), and computed using Eq.(35) for the slab model,Fourier transformed back to coordinate space. The minimummean free path Λ min is 10 d ; the supercell has N slabs of width d , where N is shown in the figure legend. The period N repetitions are shown. Black and red curves use the powerlaw p = 0 , Q = Λ min ( ω D /ω Q ) p . −30 −20 −10 0 10 20 30distance from mid−point−40−30−20−10010203040 t e m pe r a t u r e e xc u r s i on FIG. 4. Temperature versus distance from the midpoint (cid:96) = N/ (cid:96) = 0) and cold ( (cid:96) = N/
2) points. Allcomputations used the discrete slab formula (Eq.(35)) andΛ min = 20 d . The power law 1 /τ Q is ω Q . The cell size N was40 for the steepest curve, then 80, 160, 320, 640, and 2560. / Λ + 1 /L . However, their cell has periodic boundaryconditions and therefore no actual boundary. The dis-crete version of the PBE presented here uses correct sta-tistical theory and incorporates the ring geometry by de-sign. Therefore, it should correctly describe the rate atwhich κ eff ( L ) converges to the L → ∞ limit, and provideguidance for extrapolation. Fig. 5 shows how κ eff ( L )converges as L increases. Models with constant mean -3 -2 -1 L min /L10 -4 -3 -2 -1 k / k e ff ( L ) - p=0p=1p=2 FIG. 5. Rate of convergence of 1 /κ eff ( L ) to its L → ∞ limit, by a logarithmic plot as a function of Λ min /L . Allcomputations used the discrete slab formula (Eq.(35)) andΛ min = 10. From bottom to top, the power law p of 1 /τ Q ∝ ω pQ is p = 0, 1, and 2. The numerical power law of κ /κ − ∝ (1 /L ) s is s ≈ ∞ , 2, and 1/2. free path converge exponentially, as exp( − /L ). Thisis true both for continuum and discrete cases. This is notsurprising, but I have not yet found a simple proof. Thefactor of 4 just relates to the fact that only 1/4 of thering (∆ m = N/
4) is available for decay. Discrete mod-els with mean free paths diverging as 1 /ω ( p = 1) and1 /ω ( p = 2) are also shown in Fig. 5. They apparentlyconverge algebraically, as 1 /L for p = 1 and 1 / √ L for p = 2. These powers were found numerically from theslope of the log-log graphs. I conjecture that the behav-ior is κ bulk − κ eff ( L ) ∝ (1 /L ) (3 − p ) /p . This scaling for p = 2 is in rough accord with the highest temperaturesimulation by Zhou et al. , where convergence was non-linear in 1 /L . A plot of their numerical results versus / √ L instead of 1 /L gives a significantly better straightline. Unfortunately, the extrapolated value then falls be-low 1 /κ = 0, indicating that even larger simulation cellsare needed before extrapolation can be relied on. Resultsfor GaN by Lindsay et al. , using DFT and PBE, agreewell with experiment. They also show that isotopicallypure GaN will have κ (300 K ) ≈ /κ closer to 0 than to the extrapolated values ofZhou et al. . VII. DISCUSSION
Non-local heat transport is seldom explicitly dis-cussed. The non-locality is hidden if the temperaturegradient is uniform. Then the form of the non-locality(contained in κ ( k ) for a homogeneous system) is irrel-evant, since only the k → κ ( q ) (Eq.(35)) has, I believe, negligible discretization er-ror, and correctly includes the nanoscale corrections thatenter when mean free paths are comparable to the simu-lation cell length L in an MD simulation. This situationis hard to avoid for materials with good crystalline orderand weak anharmonicity, such as the GaN simulation ofref. 13.Picturing heat transport as explicitly non-local mayhave benefits. For example, in MD simulations of κ bythe “direct method,” it would be sensible to impose heatin a periodic fashion, ˙ e ( (cid:96) ) = cos(2 π(cid:96)/N ), or ˙ e ( q ) con-taining only the smallest non-zero q = 2 π/N d allowed.This simplifies Eqs.(25,26). More important, it shouldenable extrapolation to the q → T ( (cid:96) ) wouldbe used at all (cid:96) to extract T ( q ), rather than using only a few points of T ( (cid:96) ) to extract a gradient at the midpoint.A future paper on this topic is planned.A number of recent papers formulate theories of heatconductivity in nanoscale systems . Microscopic the-ory is needed, not just to supplement simulation, butmore importantly, to aid experiment in interpretingnanoscale effects. Non-local effects are evident in heattransport by nanoscale samples of good conductors likegraphene. Landauer methods are often preferred,but have limitations when inelastic scattering is present.Meir-Wingreen-type approaches often used to supple-ment Landauer methods in electron transport, and canbe generalized to heat transport . The PBE, which isevidently useful to model non-locality, and includes in-elasticity, can perhaps be exploited in new ways to sim-plify and unify some of these problems. VIII. ACKNOWLEDGEMENTS
This work was suggested by discussions with M.-V.Fernandez-Serra and J. Siebert, whose stimulation isgratefully acknowledged. I also thank D. Broido, G.Chan, Y. Li, K. K. Likharev, J. Liu, A. J. H. McGaughey,S. Ocko, T. Sun, R. M. Wentzcovitch, and an anonymousreferee for helpful input. This work was supported in partby DOE grant No. DE-FG02-08ER46550. ∗ [email protected] G. D. Mahan and F. Claro, Phys. Rev. B , 1963 (1988). C. Dames and G. Chen, Thermal Conductivity of Nanos-tructured Thermoelectric Materials, CRC Handbook,edited by M. Rowe (Taylor and Francis, Boca Raton, FL,2006), Ch. 42. G. Chen,
Nanoscale energy Transport and Conversion (Ox-ford University Press, Oxford, 2005), Ch. 7. Z. M. Zhang,
Nano/Microscale Heat Transfer (McGraw-Hill, New York, 2007), Ch. 5. D. G. Cahill, W. K. Ford, K. E. Goodson, G. D. Mahan,A. Majumdar, H. J. Maris, R. Merlin, and S. R. Phillpot,J. Appl. Phys. , 793 (2003). A. Majumdar, J. Heat Transfer , 7 (1993). J. A. Johnson, A. A. Maznev, J. Cuffe, J. K. Eliason, A.J. Minnich, T. Kehoe, C. M. Sotomayor Torres, G. Chen,and K. A. Nelson, Phys. Rev. Lett. , 025901 (2013). A. A. Maznev, J. A. Johnson, and K. A. Nelson, Phys.Rev. B , 195206 (2011). K. C. Collins, A. A. Maznev, Z. Tian, K. Esfarjani, K. A.Nelson, and G. Chen, J. Appl. Phys. , 104302 (2013). D. N. Payton, M. Rich, and W. M. Visscher, Phys. Rev. , 706 (1967). F. M¨uller-Plathe, J. Chem. Phys. , 6082 (1997). A. Maiti, G. D. Mahan, and S. T. Pantelides, Solid StateCommun. , 517 (1997). X. W. Zhou, S. Aubry, R. E. Jones, A. Greenstein, and P.K. Schelling, Phys. Rev. B , 115201 (2009). B.-Y. Cao and Y.-W. Li, J. Chem. Phys. , 024106(2010). D. P. Sellan, E. S. Landry, J. E. Turney, A. J. H. Mc-Gaughey, and C. H. Amon, Phys. Rev. B , 214305(2010). R. Peierls, Ann. Phys. , 1055 (1929). J. M. Ziman,
Electrons and Phonons (Oxford UniversityPress, London, 1960). A. J. Minnich, Phys. Rev. Lett. , 205901 (2012). C. Hua and A. J. Minnich, Phys. Rev. B , 094302 (2014). W. Li, J. Carrete, N. A. Katcho, and N. Mingo, ComputerPhys. Commun. , 1747 (2014). C. Herring, Phys. Rev. , 954 (1954). A. Ward and D. A. Broido, Phys. Rev. B , 085205(2010). K. Esfarjani, G. Chen, and H. T. Stokes, Phys. Rev. B ,085204 (2011). A. Akhieser, J. Phys. (USSR) , 277 (1939). H. J. Maris, in
Physical Acoustics , edited by W. P. Masonand R. N. Thurston (Academic, New York, 1971), Vol. 8,p. 279. L. Lindsay, D. A. Broido, and T. L. Reinecke, Phys. Rev.Lett. 109, 095901, (2012) G. Lebon, H. Machrafi, M. Grmela, and Ch. Dubois, Proc.R. Soc A , 3241 (2011). S. G. Das and A. Dhar, Eur. Phys. J. B , 372 (2012). K. S¨a¨askilahti, J. Oksanen, and J. Tulkki, Phys. Rev. E , 012128 (2013). F. Yang and C. Dames, Phys. Rev. b , 035437 (2013). R. Landauer, Z. Phys. B - Condensed Matter , 217(1987). D. E. Angelescu, M. C. Cross, and M. L. Roukes, Super-lattices Microstruct. , 673 (1998). N. Mingo, Phys. Rev. B , 113308 (2003). Y. Meir and N. S. Wingreen, Phys. Rev. Lett.68