PHYSICAL REVIEW E 90, 032112 (2014)

Heat conduction in nanoscale materials: A statistical-mechanics derivation of the local heat flux Xiantao Li* Department of Mathematics, The Pennsylvania State University, University Park, Pennsylvania 16802, USA (Received 12 June 2014; published 11 September 2014) We derive a coarse-grained model for heat conduction in nanoscale mechanical systems. Starting with an all-atom description, this approach yields a reduced model, in the form of conservation laws of momentum and energy. The model closure is accomplished by introducing a quasilocal thermodynamic equilibrium, followed by a linear response approximation. Of particular interest is the constitutive relation for the heat flux, which is expressed nonlocally in terms of the spatial and temporal variation of the temperature. Nanowires made of copper and silicon are presented as examples. DOI: 10.1103/PhysRevE.90.032112

PACS number(s): 05.70.Ln, 44.10.+i, 61.46.Bc, 46.15.−x

I. INTRODUCTION

Over the past two decades, there has been increasing interest in modeling energy transport at the microscopic scale. Such effort has been driven by the progress in designing and manufacturing micromechanical and electrical devices, for which thermal conduction properties have great influences on the performance. Many previous studies (see, e.g., [1–14]) have focused on the calculation of the heat conductivity and its dependence on the shape, length, and temperature of the system. A typical setup is to connect the boundaries to two heat baths with different temperature, modeled by stochastic (Langevin) or deterministic (Nos´e-Hoover [15–17]) thermostats. The molecular dynamics (MD) equations are solved in the interior to drive the system to a steady state, at which point the heat flux can be measured to estimate the heat conductivity. More general transient heat conduction problems, however, pose a much greater challenge for direct nonequilibrium molecular dynamics (NEMD) simulations at the atomic scale. For instance, while it is often relatively straightforward to incorporate quantities such as displacement, velocity, temperature, and pressure into MD simulations as constraints, the temperature gradient is very difficult to impose. The temperature gradient that can be imposed is usually on the order of 107 –108 K/m, which is too large to model realistic systems. It is unclear whether results obtained from NEMD can be appropriately extrapolated to the correct regime. Meanwhile, the linear response theory takes advantage of the small spatial or temporal variations by a perturbation around a thermodynamic equilibrium [18]. The most important result from this traditional theory is the Green-Kubo formula for transport coefficients. The derivation is well justified for mechanically driven systems. For a thermally driven system, however, the mathematical derivation usually requires assumptions that are difficult to state or verify. In addition, the Green-Kubo formula determines the heat conductivity in terms of the correlation of the heat flux averaged over the entire system. It does not offer a formula for the heat conductivity at a local point, which is particularly important for a heterogeneous sample. Recently, Petravic and Harrowell [19] generalized the linear response approach and with great mathematical clarity

*

[email protected]

1539-3755/2014/90(3)/032112(10)

derived a formula for the heat conductivity in terms of the fluctuation at the boundary. This paper aims to develop a more general heat conduction model. Motivated by the comments in the review paper [20], we present a coarse-grained description as a more efficient way to model heat conduction problems. Instead of particle trajectories, we choose to work with locally averaged quantities, e.g., average displacement, velocity, and energy. The time evolution of these coarse-level quantities, driven by the inhomogeneity of the temperature, exactly follows a system of conservation laws. In principle, this offers a means to consider thermoelastic processes at a much larger scale. However, to obtain a complete description at this scale, a constitutive relation is needed for the fluxes in the conservation laws. The second part of the coarse-graining procedure is to determine an explicit formula for the local heat fluxes along each crystal plane. We propose a linear response approximation, relating the energy flux to the temperature gradient. It may appear to the readers that this is already a standard approach in statistical physics [18]. The main departure from the traditional approach is a linear response approximation around a quasilocal thermodynamic equilibrium, which allows the temperature to be nonuniform. To obtain the quasilocal equilibrium, we adopt the statistical-operator approach [21,22]. By using nanowires as a specific example, we prove that the quasilocal equilibrium provides an approximate solution to the Liouville equation associated with the MD model. However, the accuracy of the approximation is only guaranteed for a short time. More importantly, the average heat flux with respect to this approximate phase-space density is zero. In order to extend this approximation to a longer time scale and predict the heat flux, we introduce the linear response approximation, which leads to a closed-form expression for the heat flux. An important observation from this model is that the resulting constitutive relation is highly nonlocal: It depends on the temperature gradient everywhere in the system. The nonlocality in time is reflected in a convolutional integral. Namely, the flux is also determined by the history of the temperature. This can be viewed as one generalization of Fourier’s law, typically observed at the macroscopic scale. Although nonlocal models have been suggested, e.g., by Tzou [23], the model derived in this paper is obtained directly from the full MD model instead of empirical observations. Further, the parameters have been shown to have their roots at the atomic scale.

032112-1

©2014 American Physical Society

XIANTAO LI

PHYSICAL REVIEW E 90, 032112 (2014)

In the following, we first present a careful mathematical derivation of the conservation laws, with the hope to clarify the correct molecular expression of the heat flux, which, unfortunately, has been a subject of debate [24]. We then follow the statistical-operator approach to arrive at the quasilocal thermodynamic equilibrium. Section II C presents the linear response approach and the resulting expression for the heat flux.

We let the domain I be one of the cells. Within this cell, we define the average mass, displacement, momentum, and energy  1  mI = mi , q I = mi r i , mI i∈ i∈I I (2)   mi v i , EI = Ei . pI = i∈I

II. MATHEMATICAL DERIVATIONS A. Problem setup

We have in mind a crystalline system with N atoms, the motion of which is determined by the Newton equation r˙ i = v i , mi v˙ i = −∇ r i U.

(1)

Here mi , r i , and v i are, respectively, the mass, coordinate, and velocity of the ith atom; U = U (r 1 ,r 2 , . . . ,r N ) is the interatomic potential, which we will assume to be empirical and short ranged. Initially, the system is prepared with some random distribution, which will be discussed later. The model based on (1) is known as the MD model. It provides a fundamental description of many important aspects of physics, which however, may not be easy to access directly. Our principal objective is to derive simpler models, in which quantities of our direct interest are the primary variables. A coarse-graining (CG) procedure of this type for MD models has been one of the most important and challenging problems in molecular modeling [25]. Most existing methods (see, e.g., [26–52]) work with the canonical ensemble at a constant temperature and do not apply to the heat conduction problems where the temperature is not uniform. As the first step we partition the system into cells, each of which contains multiple atoms, as demonstrated in Fig. 1. In practice, this can be done by defining supercells, each of which contains multiple primitive cells of the crystal lattice. Alternatively, one may start with a set of points in the domain and perform a Voronoi tessellation [53].

i∈I

As is customary in solid mechanics, we follow the reference coordinate. Thus, the notion i ∈ I simply indicates that the reference position of the ith atom is in the domain I . The reference position is associated with the underlying lattice. Intuitively, q I and pI represent, respectively, the average coordinate and momentum in the cell. In continuum mechanics [54], heat flux is often introduced at the level of energy balance. Therefore, we have also defined the average energy EI . This requires that the total potential energy be written as  Ui , (3) U= i

in which Ui is the energy due to atoms surrounding the ith atom. In (2), the average energy is then given by Ei =

mi v 2i + Ui . 2

(4)

It is worthwhile to point out that this step is motivated by how continuum mechanics models are derived: One starts with a general control volume and then the notions of momentum and heat fluxes will follow from the fundamental conservation laws. However, instead of postulating a constitutive relation, we first express the fluxes in terms of molecular trajectories by following (1) to ensure the microscopic origin. By taking time derivatives, we find that  def d pI = tI = f i, (5) dt i∈ I

with f i = −∇ r i U. To interpret the right-hand side as a flux instead of a body force, we express f i as  f ij , f ij = − f j i . (6) fi = j =i

ΩI Ω

J

FIG. 1. (Color online) Partition of the system: The atoms are grouped into different cells. The interface between the neighboring cells I and J is highlighted. This partition is generated from a Voronoi tessellation.

It might be intuitive that f ij represents the force on the ith atom due to the j th atom; however, this does not imply that the interaction is pairwise since f ij may depend on the position of other atoms. As will be shown below, this decomposition is important for defining the fluxes locally at the cell interfaces. This assumption, along with the assumption of the energy partition (3), deserves much attention even though all the existing empirical potentials have natural decomposition of those forms [55]. We have addressed the issues of consistency in a separate work [56]. Next, with a substitution into (5), one finds cancellations when both i and j are inside I . Therefore, we get  d f ij . (7) pI = dt i∈ ,j ∈ /

032112-2

I

I

HEAT CONDUCTION IN NANOSCALE MATERIALS: A . . .

PHYSICAL REVIEW E 90, 032112 (2014)

This suggests that we define the momentum flux (traction) between two neighboring cells, e.g., I and J , as  f ij . (8) t I,J = i∈I ,j ∈J

Clearly, we have t I,J = −t J,I . In addition, the conservation of momentum is satisfied,  d pI = t I,J . (9) dt J Due to the short-range interaction, this flux term only depends on the coordinate of the atoms near the interface. As an example, for the embedded-atom model (EAM) [57], the potential energy is expressed in the function form   1  U = ϕ(rij ) + E(ρi ), 2 1iN 1j N,j =i 1iN (10)  ρi = ρ(rij ). 1j N,j =i

Here we have adopted the usual notation in molecular simulations: r ij = r i − r j and rij = |r ij |. The function ρ represents the influence of local electron density and because of the nonlinearity of the function E, the interaction is of multibody nature. For this model, the most commonly used decomposition is given by r ij f ij = −ϕ  (rij ) − [E  (ρi ) + E  (ρj )]ρ  (rij ) , rij (11) 1 ϕ(rij ) + E(ρi ), Ui = 2 j =i which clearly satisfies (6). Of more direct interest to our problem is the energy balance. So we calculate the time derivative dtd EI . For the EAM potential, with the force and energy decomposition (11), the energy flux can be expressed as  r ij 1  d ϕ  (rij ) EI = · (v i + v j ) dt 2 i∈ ,j ∈ rij / I I  r ij . (12) − [E  (ρi )v j + E  (ρj )v i ] · rij Similar to (8), we define the local energy flux  r ij 1  qI,J = ϕ  (rij ) · (v i + v j ) 2 i∈ ,j ∈ rij I J  r ij   − [E (ρi )v j + E (ρj )v i ] · . rij

FIG. 2. (Color online) Coarse graining the copper nanowire. The system is aligned in the 100 direction with dimension 46.3 × 14.5 × 14.5 nm3 . For the CG procedure, it is divided into blocks, in each of which the total energy (momentum) is chosen as CG variables. The flux jI +1/2 = qI,I +1 is defined locally at the interface of two blocks and will be computed from (33).

is tempting to use the formula 1  q I,J = f ij · (v i + v j ), 2 i∈ ,j ∈ I

which can be easily derived for pair potentials. However, this is in general inconsistent with the definition of local energy (4). In fact, the definition of the energy flux has received a great deal of attention recently [58–60]. In particular, the issue of the consistency with the energy balance (14) has been addressed [56]. We have also derived a consistent expression of the fluxes for the Tersoff potential [61], another popular model used in practice. The expressions, however, are too lengthy to be included here. In light of the heavy details, we have presented the results in a separate work [56]. What has been derived so far is merely the conservation laws associated with the averaged quantities. Although they do represent the dynamics at a much coarser level, Eqs. (9) and (14) are not closed, i.e., there is no constitutive relation. The next section is devoted to a closure procedure, which will connect the fluxes to the CG quantities. B. Quasilocal equilibrium

In order to best describe the ideas, let us consider a nanowire system, similar to the one shown in Fig. 2. We divide the system into cells along the longitudinal direction. In this case, the flux between the cells I and I +1 are tI,I +1 and qI,I +1 . The energy balance is simply expressed as d EI (t) = qI,I +1 (t) − qI −1,I (t). (15) dt Ideally, we would like to take the ensemble averages of the conservation laws (9) and (15) and only focus on the dynamics of the averaged quantities. For this purpose, we need to consider the probability distribution associated with the MD model. Let us introduce the notation r = (r 1 ,r 2 , . . . ,r N ), v = (v 1 ,v 2 , . . . ,v N ).

(13)

Again we have qI,J = −qJ,I and the conservation of energy  d EI (t) = qI,J . (14) dt J The summation on the right-hand side is over all the neighboring cells. Remark 1. A correct expression of the energy flux is obviously crucial for studying heat conduction processes. It

J

(16)

The probability density ρ(r,v,t) for the MD model (1) is governed by the Liouville equation [62,63] ρt + Lρ = 0,

ρ(·,0) = ρ0 (·).

(17)

Here L is the Liouville operator and the dot in parentheses indicates the phase-space variables, i.e., ρ(·,t) = ρ(r,v,0). The many-particle Louiville equation (17), which in principle is equivalent to the MD model, provides the full description of the system. The passage from molecular trajectories to the

032112-3

XIANTAO LI

PHYSICAL REVIEW E 90, 032112 (2014)

average Aˆ of a quantity A can be expressed as  ˆ = A(t) A(r,v)ρ(r,v,t)d r dv.

Notice that βI =

(18)

R6N

It is in general not practical to directly solve the Liouville equation numerically. However, as we will show, it is possible to find approximate solutions parametrized by a few CG variables. In particular, we seek a local eq thermodynamic equilibrium, denoted by ρloc and referred to as quasilocal equilibrium. This approach, motivated by the nonequilibrium statistical-operator method for quantum statistical physics [21,22], is particularly useful for systems eq far from equilibrium. The first basic requirement for ρloc is that it is non-negative and its integral is one:  eq ρloc d r dv = 1. (19)

R6N

 Eˆ I (t) =

R6N

(20) eq EI ρloc (r,v,t)d r

dv.

These conditions make certain that the approximate density yields the corresponding averages. Under these constraints, we seek a particle distribution using the maximum entropy principle. More specifically, for a fixed t, we consider the functional  S[ρ] = − ρ(r,v,t) ln ρ(r,v,t)d r dv − (L − 1)  × −



 ρ(r,v,t)d r dv − 1

I

−





 p

LI ·  LE I

pI (r,v)ρ(r,v,t)d r dv − pˆ I  ˆ EI (r,v)ρ(r,v,t)d r dv − EI . (21)

I

The first term in S[ρ] is the entropy and the following terms are penalty terms that enforce the constraints with L, L p , and LE being the Lagrange multipliers. In this paper we consider only the problem of heat conduction and we neglect the influence of the average displacement and velocity. In this case, the variational approach yields    exp − I βI (t)EI eq ρloc (r,v,t) = . (22) Z(t) Here Z(t) plays a role similar to the partition function and βI is determined from the conditions  eq Eˆ I (t) = ρloc (r,v,t)EI d r dv. (23) Substitution of (22) into (21) yields the entropy  S= βI EI + ln Z.

(25)

It is immediately recognized as the inverse temperature in the I th cell: βI = kB1TI . Furthermore, we have ∂S = 0. ∂βI Therefore, the entropy S is only a function of Eˆ I . The quasilocal equilibrium (22) is a generalization of the canonical ensemble and allows the temperature to be inhomogeneous. In this case, the heat capacity becomes a tensor and is given by 1 ∂ Eˆ I = (EI − Eˆ I )(EJ − Eˆ J ). ∂TJ kB TJ2

R6N

The second requirement is that the ensemble averages of eq the CG variables with respect to ρloc are consistent with the observed values, here denoted by pˆ I and Eˆ I . Namely, we have  eq pˆ I (t) = pI ρloc (r,v,t)d r dv,

∂S . ∂ Eˆ I

(26)

Remark 2. Notice that in (2) the average mass is independent of the phase variables. Therefore, it has not been enforced: It can be computed a priori. Further, we have shown only the implication of the constraints from the energy. Together with the momentum constraints, the quasilocal equilibrium takes the form 

  pˆ 1 eq ρloc (r,v,t) = βI EI + I · pI . (27) exp − Z(t) mI I In the following discussion we assume that the average momentum can be neglected so that we can focus on the heat conduction. In this case, the energy flux q coincides with the heat flux since there is no convection. Remark 3. We have chosen to use the inverse temperature β. This is similar to the continuum mechanics formulation in [54]. However, the usual temperature T clearly can be used as well. The main departure from the usual Gibbs distribution is that the probability density (22) allows the temperature to be inhomogeneous, which makes heat conduction possible. It can be directly verified that when the averages of the CG variables are smooth, i.e., the spatial and temporal variations eq are small, the quasilocal equilibrium ρloc is an approximate solution of the Liouville equation (17). Later this will be illustrated through a one-dimensional example. C. Nonequilibrium heat flux

The quasilocal equilibrium provides a means to close the momentum and energy balance equations (9) and (14): Instead of computing the ensemble averages with respect to ρ(·,t), we eq compute them based on ρloc (·,t). The corresponding average is denoted by ·QL . This is quite similar to the closure procedure in kinetic theory, where a local Maxwellian is used to derive hydrodynamic equations [64]. For the current problem, we obtain the energy equation d EI QL = qI,I +1 QL − qI −1,I QL . dt

(28)

To close the energy equation, it remains to find the average of the heat flux. With direct calculations, one finds that (24)

I

032112-4

qI,I +1 QL = 0 ∀ I.

(29)

HEAT CONDUCTION IN NANOSCALE MATERIALS: A . . .

PHYSICAL REVIEW E 90, 032112 (2014)

eq

This is because the average velocity with respect to ρloc for each atom is zero. At this level, there is no energy transport. Conceptually, this is because heat conduction is an irreversible phenomenon and does not arise from this part of the closure procedure. In the statistical-operator approach [21,22], the terms that can be deeq rived from ρloc are referred to as reversible forces. Meanwhile, a common idea to derive dissipative forces is to seek a highereq order approximation, i.e., ρ ≈ ρloc + ρ  , where ρ  represents a correction. Numerous forms have been suggested [18,21,22], but they are difficult to justify. Further, the well established techniques in kinetic theory, including the Hilbert and Chapman-Enskog expansions [64], do not apply here because of the absence of a collision term and small parameters. Here we propose a different approach, motivated by the classical linear response theory [18]. We observe from direct calculations that ∂ eq eq eq ρ + Lρloc = ζ (r,v,t)ρloc . (30) ∂t loc where  ζ (r,v,t) = − {[βI +1 (t) − βI (t)]qI,I +1 I

+ [EI − Eˆ I (t)]β˙I (t)}.

(31)

When the temperature is smooth in space and time, the function ζ is small. This can be justified on the ground that the real temperature gradient, velocity gradient, and strain eq rate are small at the atomic scale. Therefore, ρloc satisfies a modified equation with a small right-hand side. Based on this observation, one can use techniques in the theory of differential eq equations and show that ρ(·,t) ≈ ρloc (·,t), at least for short time scales. Next we assume that initially the system is at a quasilocal equilibrium. This is often interpreted as a good preparation of the system. We express the solution of the Liouville equation (17) using the semigroup operator [63] eq

ρloc (·,t) = e−tL ρloc (·,0) + eq



eq

0

J

J

¯

κI,J (t) = qI,I +1 (t),qJ,J +1 , ηI,J (t) = qI,I +1 (t),EJ . (38) It is also reasonable to assume the stationarity in space, i.e., κI,J (t,s) = κI −J (t − s) and ηI,J (t,s) = ηI −J (t − s). These functions appear in the convolutional integral in (36) and will also be referred to as kernel functions. Higher-order approximations will involve the correlation of three or more quantities. Here κI,J and ηI,J are not independent parameters. In particular, from (15) we deduce that η˙ I,J (t) = −κI,J (t) + κI,J −1 (t).

eq



+

ηI,J (t − s)βJ (s)|t0

 J

=− +

(34)

(35) eq

[−κI,J (t − s) + κI,J −1 (t − s)]βJ (s)ds

0



ηI,J (t)βJ (0)

 J

qI,I +1 (t − s)EJ ρloc (r,v,s)d r dv.

t

J

The second term can now be interpreted as a correction, embodying the nonequilibrium contribution. This constitutes the basis for deriving the dissipative force for transport phenomena such as heat conduction, which we now discuss. We compute the average heat flux using (34). It is useful for this calculation to define the correlation functions  eq κI,J (t,s) = qI,I +1 (t − s)qJ,J +1 ρloc (r,v,s)d r dv,

R6N

0

=

0



(39)

Using this relation, one can use integration by parts and rewrite the second term in (36) as  t ηI,J (t − s)β˙J (s)ds

(33)

ηI,J (t,s) =

(36)

e−βH , (37) Z where β¯ is the reference temperature, e.g., the average temperature in the entire system. We regard this approximation as a linear response approximation and can simply write the correlation functions as

e−(t−s)L ζ (·,s)ρloc (·,s)ds.

R6N

ηI,J (t,s)β˙J (s)ds. 0

ρ eq =

(32)

Combining the equations, we find that  t eq eq ρ(·,t) = ρloc (·,t) − e−(t−s)L ζ (·,s)ρloc (·,s)ds.

t

Here we neglected the effect of the boundary. Boundary conditions will be discussed in the next section. So far, the formula (36) is exact and no approximation has been made. The expressions (35) would simplify if we replace eq ρloc by the canonical ensemble

J t

0



+

J

ρ(·,t) = e−tL ρloc (·,0). Further, from (30) we have

Then the average heat flux is given by  t κI,J (t,s)[βJ +1 (s) − βJ (s)]ds qI,I +1  = −

t

κI,J (t − s)[βJ +1 (s) − βJ (s)]ds.

(40)

0

To some extent, the CG model we have derived shares some similarity with nonlocal diffusion models, e.g., [65–69]. However, the models proposed here are derived based on the full MD models instead of empirical observations. D. Thermal boundary condition

In our previous discussion we have neglected the boundary condition. Typically, the temperature at the boundary is maintained using thermostating techniques. As an explicit example we consider the Nos´e-Hoover thermostat at the real

032112-5

XIANTAO LI

PHYSICAL REVIEW E 90, 032112 (2014)

fi − ζL v i , mi

M  2 ˙ζL = 1 mi v i − 3MkB TL . QL i=1

(41)

M  i=1

 ∂ ∂ ∂ ρ+ ρ + ζ˙L ρ − 3Mζρ. v˙i · ∂ xi ∂v ∂ζ i L i=1 M

p˙ i ·

0.02

0

1

2

3

4

5

6

time

Here we have assumed that the first M atoms are in the heat bath on the left with temperature TL . The added variable ζL is to control the temperature. The Liouville operator associated with the Nos´e-Hoover dynamics is given by [16] Lρ =

L=44.80 nm L=22.40 nm

−0.02 0 0.01

kernel function

x˙ i = v i , v˙ i = −

0.04

kernel function

time scale [16]. In particular, for the atoms at the left boundary, we follow the equations

L=44.80 nm L=22.40 nm

0.005 0 −0.005 −0.01 0

1

2

3

4

5

6

time

(42)

The energy for the last cell on the right is defined in a similar way. Then, with a direct calculation, one can show that the local equilibrium defined in (22) still satisfies (30). As a result, the formula for the heat flux (36) still holds. III. EXAMPLES

In this section we present results from several numerical tests. In these experiments we work with nanowires made from copper or silicon. It is clear from our derivation that the theory also applies to other systems, e.g., graphite sheets; however, we choose the quasi-one-dimensional systems because of the simple forms of the energy flux. A. Copper nanowire

The first example is a copper nanowire, as shown in Fig. 2. We align the system in the horizontal direction 100. We choose the reference temperature to be T = 100 K and at this temperature, the lattice spacing is 0.3621 nm. We also choose the embedded-atom potential [70]. The cross section of the wire has a square shape with a width of 1.811 nm. We use the consistent expression for the energy flux (13) in all of the following tests. We will use the standard unit in molecular simulations. Namely, we use angstroms as the length unit and eV as the energy unit (from the EAM [70]). Such selection will determine the units of other quantities. In particular, the corresponding unit for time is 0.081 155 ps and the unit for the correlation function κ is 5.9179 × 10−6 W2 . We first compute the correlation functions with respect to ρ eq in (37), i.e., the canonical ensemble. In the computation we use the Nos´e-Hoover thermostat technique [15–17] to equilibrate the system to the desired temperature. Then the

simulation is continued for 5 × 106 steps to generate enough samples. Numerous tests were carried out for copper wires of different length. Figure 3 shows the correlation functions computed from these tests. We can see that even when the wires are of different length, the correlation functions are in good agreement, at least within the measured time period. This seems to suggest that the correlation functions are insensitive to the length of the system. Next we compare the computed correlation functions to the correlation functions for the bulk, which is modeled by periodic boundary conditions in all directions. As one can observe from Fig. 4, these correlation functions are clearly different. The difference can be attributed to the special

0.04

kernel function

0

FIG. 3. (Color online) Comparison of the correlation functions of the heat flux: the kernel function κI,I (t) (top) and the kernel function κI,I +1 (t) (bottom). The unit for time is 0.081 155 ps and the unit for the kernel function κ is 5.9179 × 10−6 W2 .

Nanowire Bulk 0.02

0

−0.02 0

0.5

1

1.5

2

2.5

3

3.5

4

time 0.01

kernel function

Similarly, a heat bath can be introduced at the boundary on the right to maintain the system at a different temperature, e.g., TR . In the presence of the heat bath at the boundary, we define the energy in the first cell as

 1 1 (43) E0 = Vi + mi v 2i + QL ζL2 . 2 2 i∈

Nanowire Bulk

0.005 0 −0.005 −0.01 0

0.5

1

1.5

2

2.5

3

3.5

4

time

FIG. 4. (Color online) Comparison of the correlation function of the heat flux in the nanowire and in the bulk: the kernel function κI,I (t) (top) and the kernel function κI,I +1 (t) (bottom). The unit for time is 0.081 155 ps and the unit for the kernel function κ is 5.9179 × 10−6 W2 .

032112-6

PHYSICAL REVIEW E 90, 032112 (2014)

kernel function

HEAT CONDUCTION IN NANOSCALE MATERIALS: A . . .

kernel function

0.02

0.01

0

0.04

L=39.53 nm L=19.57 nm

0.02 0 −0.02 0

0.5

1

1.5

2

2.5

time −0.01

kernel function

−3

−5

I−J

0

5 0

1

2

3

4

5

6

time

4

x 10

L=39.53 nm L=19.57 nm

2 0 −2 −4 0

0.5

1

1.5

2

2.5

time

FIG. 5. (Color online) Plot of the correlation functions κI,J (t). The unit for time is 0.081 155 ps and the unit for the kernel function κ is 5.9179 × 10−6 W2 .

geometry of the system since the four sides of the nanowire are free surfaces and the surface to volume ratio is large. Another observation from Figs. 3 and 4 is that the peak of the correlation function κI,J (t) will be further delayed as |J − I | becomes larger. Figure 5, where we show several correlation functions in the same plot, clearly confirms this observation. The behavior is reminiscent of a wave source starting at the middle and propagating in both directions. This might be related to what was observed in the energy diffusion [71]. B. Silicon nanowire

The setup of the system is similar to the copper wire in the previous section. The equilibrium lattice spacing at 100 K is a0 = 0.5499 nm. The system is generated from rectangular unit cells, each of which contains eight atoms. The entire sample is of rectangular shape with a width of 2.1996 nm. For the simulations we use the Tersoff potential [61], which essentially is a multibody potential. The molecular unit for the time scale is 0.029 109 ps and the unit for the correlation functions is 3.0296 × 10−5 W2 . We first conduct equilibrium simulations at the temperature T = 100 K and sample the correlation functions. The results are shown in Fig. 6, where we collected results from two nanowires of different length. Again, the observed agreement suggests that these correlation functions are independent of the wire length. Figure 7 is a collection of kernel functions shown in the same plot. Similar to what was observed from the copper nanowire, the peaks of the kernel functions exhibit a wavepropagation type of behavior. In all the previous tests we computed the correlations at a thermodynamic equilibrium, given by ρ eq in (37) in Sec. II C. To examine this approximation, we have also computed the correlation functions around a nonuniform state with constant temperature gradient. In order to bring the system to such a state, we start with a nanowire of length 39.53 nm and apply thermostats at both boundaries. The temperature on the left and right is maintained at TL = 99K and TR = 101K,

FIG. 6. (Color online) Comparison of the correlation functions of the heat flux for two nanowires of different length: the kernel function κI,I (t) (top) and the kernel function κI,I +1 (t) (bottom). The time scale is 0.081 155 ps and the unit for the kernel function κ is 5.9179 × 10−6 W2 .

respectively. This part of the computation is similar to previous studies on the dimension dependence of thermal conductivity and the methods used here have been motivated by the work [4]. Once the system is brought to the steady state with constant temperature gradient, we compute the correlation function using the block-averaging procedure [17]. Good agreement can be observed from Fig. 8, which suggests that this approximation is quite reasonable. From this simulation we have also observed that the correlation functions are almost translationally invariant, in the sense that κI,J = κI +K,J +K , with large discrepancies observed near the heat bath at the boundaries. This is illustrated in Fig. 9.

FIG. 7. (Color online) Plot of the kernel functions κI,J (t). The time scale is 0.081 155 ps and the unit for the kernel function κ is 5.9179 × 10−6 W2 .

032112-7

XIANTAO LI

PHYSICAL REVIEW E 90, 032112 (2014)

kernel function

0.04 0.03 0.02 0.01 0 −0.01 −0.02

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

time −3

kernel function

5

x 10

0

−5 0

0.5

1

1.5

2

time

FIG. 8. (Color online) Correlation functions κI,I (t) (top) and eq κI,I +1 (t) (bottom) with respect to ρ eq (solid lines) and ρloc (dashed lines). The time scale is 0.081 155 ps and the unit for the kernel function κ is 5.9179 × 10−6 W2 . IV. SUMMARY AND DISCUSSION

This work was motivated by many observations that conventional heat conduction models, represented by the Fourier law, are inadequate to model heat conduction at the microscopic scale. On the other hand, direct molecular dynamics simulations tend to be too expensive and the results contain substantial fluctuations. In this paper we have presented a coarse-graining approach, aiming to describe energy transport at an intermediate scale. We derived the balance equations governing the time evolution of average momentum and energy, which is a discrete analog of continuum mechanics models, but with a generalized constitutive relation. We emphasize three points: (a) The balance laws are derived

directly from the molecular dynamics model and therefore it is guaranteed that we obtain the correct expressions of the fluxes; (b) the constitutive equations, more specifically, the energy flux, are obtained from a linear response approach around a quasilocal equilibrium distribution, which allows nonuniform temperature distribution; and (c) unlike continuum mechanics models, typically expressed in terms of a local spatial gradient, the constitutive law found in this work exhibits a nonlocal dependence on the temperature variations in both space and time. Meanwhile, much further effort is needed toward understanding the fundamental issues in thermal transport in nanoscale materials. The most well known outstanding problem is the dimension dependence of heat conductivity, which for nanomechanical systems has been extensively studied [2–5,11,71–75]. Let us look at the model (36) from this perspective, by naively taking t → +∞. If the temperature gradient, given by Tx , remains steady in time, we have β˙I = 0. Further, βI +1 − βI = hβx = −

kernel function

0.03

κ4,4

0.02

1,1

Direct calculations show that ηI,J (0) = 0, mainly because the heat flux is a linear function of the particle velocity, while the energy is quadratic in v. From (39) and assuming that lim ηI,J (t) = 0,

t→+∞

we find that



+∞

0.01 0 −0.01 0.5

1

1.5

2

1.5

2

time kernel function

−3

4

x 10

κ

4,5

2

κ

1,2

0 −2 −4 0

0.5

1

 κI,J (t)dt =

+∞

κI,J −1 (t)dt.

(46)

0

Continuing with this argument, one can further deduce that the integral is independent of J . Therefore, the expression (45) would be reduced to  +∞ L qI,I +1  = − T κI,I (t)dt. (47) x kB T 2 0

κ

−0.02 0

(44)

where h is the length of each block. We assume that there are M blocks in the system, so the length is given by L = Mh. Substitution of (44) into (36) yields M  +∞  h qI,I +1  = − Tx κI,J (t)dt. (45) kB T 2 J =1 0

0 0.04

h Tx , kB T 2

time

FIG. 9. (Color online) Correlation functions κI,I (t) (top) and eq κI,I +1 (t) (bottom) with respect to ρ eq (solid lines) and ρloc (dashed lines). The time scale is 0.081 155 ps and the unit for the kernel function κ is 5.9179 × 10−6 W2 .

This is beginning to look like the Green-Kubo formula and one might conclude that the heat conductivity is linearly proportional to the length of the system. Although there is numerical evidence that seems to suggest this pattern [76], the observations from most of the existing works indicate a different form: κ ∼ Lb with b = 1, especially when the length is larger than the phonon mean free path [71,76]. It has also been suggested by Dhar [20] that instead of taking t → +∞, one should introduce a cutoff time tc that is proportional to the length of the system. We feel that the nonlocal model (36) will be useful in defining local heat conductivity and understanding the dependence on the system geometry; however, further quantitative analysis is needed. Further, the current derivation is based on the classical mechanics of Newton’s equations. Therefore, no quantum effect is considered. This may introduce a large error at low

032112-8

HEAT CONDUCTION IN NANOSCALE MATERIALS: A . . .

PHYSICAL REVIEW E 90, 032112 (2014)

temperature. A more useful approach should start with a quantum or semi-quantum-mechanical model, such as those in [77–82], and derive a coarse-grained description. The current approach still seems promising in this case. In fact, the statistical-operator approach was originally applied to quantum-mechanical systems. The practical difficulty is that the computation for the correlation functions might be too great. Finally, it is our hope to turn the coarse-grained models into a useful computational tool to simulate full thermoelastic behavior at the nanoscale. For this purpose, several

issues concerning the implementation of the coarse-grained models need to be addressed. In particular, one needs a more practical tool to obtain the correlation functions and a reliable numerical method for solving the integro-differential equations.

[1] S. Gill, Z. Jia, B. Leimkuhler, and A. Cocks, Phys. Rev. B 73, 184304 (2006). [2] K. Jolley and S. Gill, J. Comput. Phys. 228, 7412 (2009). [3] J. Che, T. C¸aˇgın, W. Deng, and W. A. Goddard III, J. Chem. Phys. 113, 6888 (2000). [4] J. Chen, G. Zhang, and B. Li, J. Phys. Soc. Jpn. 79, 074604 (2010). [5] D. Donadio and G. Galli, Phys. Rev. Lett. 102, 195901 (2009). [6] A. S. Henry and G. Chen, J. Comput. Theor. Nanosci. 5, 1 (2008). [7] A. J. H. McGaughey and M. Ka, Adv. Heat Transfer 39, 169 (2006). [8] P. K. Schelling, S. R. Phillpot, and P. Keblinski, Phys. Rev. B 65, 144306 (2002). [9] S. G. Volz and G. Chen, Phys. Rev. B 61, 2651 (2000). [10] J.-S. Wang, Phys. Rev. Lett. 99, 160601 (2007). [11] L. Wang, B. Hu, and B. Li, Phys. Rev. E 86, 040101 (2012). [12] A. J. H. McGaughey and M. Kaviany, Int. J. Heat Mass Transfer 47, 1783 (2004). [13] E. Pop, D. Mann, Q. Wang, K. Goodson, and H. Dai, Nano Lett. 6, 96 (2006). [14] A. N. Volkova and L. V. Zhigileib, Appl. Phys. Lett. 101, 043113 (2012). [15] S. Nos´e, J. Chem. Phys. 81, 511 (1984). [16] W. G. Hoover, Phys. Rev. A 31, 1695 (1985). [17] D. Frenkel and B. Smit, Understanding Molecular Simulation: From Algorithms to Applications, 2nd ed. (Academic, New York, 2002). [18] M. Toda, R. Kubo, and N. Hashitsume, Statistical Physics II. Nonquilibrium Statistical Mechanics (Springer, Berlin, 1983). [19] J. Petravic and P. Harrowell, Phys. Rev. E 71, 061201 (2005). [20] A. Dhar, Adv. Phys. 57, 457 (2008). [21] R. Luzzi and A. R. Vasconcellos, Fortschr. Phys. 38, 887 (1990). [22] D. N. Zubarev, Fortschr. Phys. 18, 125 (1970). [23] D. Y. Tzou, Int. J. Heat Mass Transfer 54, 475 (2011). [24] M. H. Khadem and A. P. Wemhoff, Comput. Mater. Sci. 69, 428 (2013). [25] A. Leach, Molecular Modeling: Principles and Applications (Prentice Hall, Englewood Cliffs, 2001). [26] R. E. Rudd and J. Q. Broughton, Phys. Rev. B 72, 144104 (2005). [27] R. E. Rudd and J. Q. Broughton, Phys. Rev. B 58, R5893(R) (1998). [28] A. J. Chorin and P. Stinis, Commun. Appl. Math. Comput. Sci. 1, 1 (2006). [29] S. Izvekov and G. A. Voth, J. Chem. Phys. 125, 151101 (2006).

[30] M. A. Katsoulakis, A. J. Majda, and D. G. Vlachos, Proc. Natl. Acad. Sci. U.S.A. 100, 782 (2003). [31] X. Li, Int. J. Numer. Method. Eng. 83, 986 (2010). [32] X. Li (unpublished). [33] M. Baaden and S. J. Marrink, Curr. Opin. Struct. Biol. 23, 878 (2013). [34] H. Gohlke and M. Thorpe, Biophys. J. 91, 2115 (2006). [35] P. A. Golubkov and P. Ren, J. Chem. Phys. 125, 064103 (2006). [36] A. Gramada and P. E. Bourne, Comput. Phys. Commun. 182, 1455 (2011). [37] S. O. Nielsen, R. E. Bulo, P. B. Moore, and B. Ensing, Phys. Chem. Chem. Phys. 12, 12401 (2010). [38] W. G. Noid, J. Chem. Phys. 139, 090901 (2013). [39] W. G. Noid, J.-W. Chu, G. S. Ayton, V. Krishna, S. Izvekov, G. A. Voth, A. Das, and H. C. Andersen, J. Chem. Phys. 128, 244114 (2008). [40] P. Poulain, A. Saladin, B. Hartmann, and C. Pravost, J. Comput. Chem. 29, 2582 (2008). [41] M. Praprotnik, L. D. Site, and K. Kremer, Annu. Rev. Phys. Chem. 59, 545 (2008). [42] S. Riniker, J. R. Allison, and W. F. van Gunsteren, Phys. Chem. Chem. Phys. 14, 12423 (2012). [43] J. F. Rudzinski and W. G. Noid, J. Phys. Chem. B 116, 8621 (2012). [44] Q. Shi, P. Liu, and G. A. Voth, J. Phys. Chem. B 112, 16230 (2008). [45] M. Stepanova, Phys. Rev. E 76, 051918 (2007). [46] Z. Zhang, L. Lu, W. G. Noid, V. Krishna, J. Pfaendtner, and G. A. Voth, Biophys. J. 95, 5073 (2008). [47] P. Espa˜nol, in Novel Methods in Soft Matter Simulations, edited by M. Karttunen, I. Vattulainen, and A. Lukkarinen, Lecture Notes in Physics Vol. 640 (Springer, Berlin, 2004), pp. 69–115. [48] C. Hij´on, P. Espa˜nol, E. Vanden-Eijnden, and R. DelgadoBuscalioni, Faraday Discuss. 144, 301 (2010). [49] C. Hij´on, M. Serrano, and P. Espa˜nol, J. Chem. Phys. 125, 204101 (2006). [50] D. Kauzlari´c, P. Espa˜nol, A. Greiner, and S. Succi, J. Chem. Phys. 137, 234103 (2012). [51] D. Kauzlari´c, P. Espa˜nol, A. Greiner, and S. Succi, Macromol. Theor. Simul. 20, 526 (2011). [52] D. Kauzlari´c, J. T. Meier, P. Espa˜nol, S. Succi, A. Greiner, and J. G. Korvink, J. Chem. Phys. 134, 064106 (2011). [53] Q. Du, V. Faber, and M. Gunzburger, SIAM Rev. 41, 637 (1999).

ACKNOWLEDGMENTS

This work used the Extreme Science and Engineering Discovery Environment, which is supported by National Science Foundation Grant No. ACI-1053575.

032112-9

XIANTAO LI

PHYSICAL REVIEW E 90, 032112 (2014)

[54] W. Nowacki, Dynamic Problems of Thermoelasticity (Springer, Berlin, 1975). [55] N. C. Admal and E. B. Tadmor, J. Elasticity 100, 63 (2010). [56] X. Wu and X. Li, arXiv:1406.0925. [57] M. S. Daw and M. I. Baskes, Phys. Rev. B 29, 6443 (1984). [58] Y. Chen, J. Chem. Phys. 124, 054113 (2006). [59] A. Guajardo-Cu´ellar, D. B. Go, and M. Sen, J. Chem. Phys. 132, 104111 (2010). [60] E. B. Webb, J. A. Zimmerman, and S. C. Seel, Math. Mech. Solids 13, 221 (2008). [61] J. Tersoff, Phys. Rev. B 37, 6991 (1988). [62] R. Balescu, Equilibrium and Nonequilibrium Statistical Mechanics (Wiley, New York, 1976). [63] D. J. Evans and G. P. Morriss, Statistical Mechanics of Nonequilibrium Liquids (Academic, New York, 2008). [64] C. Cercignani, R. Illner, and M. Pulvirenti, The Mathematical Theory of Dilute Gases (Springer, New York, 1994). [65] N. Burch and R. Lehoucq, Int. J. Multiscale Comput. Eng. 9, 661 (2011). [66] M. D’Elia and M. Gunzburger, Comput. Math. Appl. 66, 1245 (2013). [67] Q. Du, M. Gunzburger, R. B. Lehoucq, and K. Zhou, SIAM Rev. 54, 667 (2012). [68] Y. Z. Povstenko, J. Thermal Stresses 28, 83 (2004).

[69] D. Y. Tzou, J. Thermophys. Heat Transfer 9, 686 (1995). [70] Y. Mishin, M. J. Mehl, D. A. Papaconstantopoulos, A. F. Voter, and J. D. Kress, Phys. Rev. B 63, 224106 (2001). [71] G. Zhang, S. Liu, and B. Li (unpublished). [72] C. W. Chang, D. Okawa, H. Garcia, A. Majumdar, and A. Zettl, Phys. Rev. Lett. 101, 075903 (2008). [73] Y. Chen, D. Li, J. Yang, Y. Wu, J. R. Lukes, and A. Majumdar, Physica B 349, 270 (2004). [74] S.-C. Wang, X.-G. Liang, X.-H. Xu, and T. Ohara, J. Appl. Phys. 105, 014316 (2009). [75] G. Zhang and B. Li, J. Chem. Phys. 123, 114714 (2005). [76] N. Yang, G. Zhang, and B. Li, Nano Today 5, 85 (2010). [77] H. Dammak, Y. Chalopin, M. Laroche, M. Hayoun, and J.-J. Greffet, Phys. Rev. Lett. 103, 190601 (2009). [78] J. L¨u and J.-S. Wang, J. Phys.: Condens. Matter 21, 025503 (2009). [79] G. Y. Panasyuk, G. A. Levin, and K. L. Yerkes, Phys. Rev. E 86, 021116 (2012). [80] A. V. Savin, Y. A. Kosevich, and A. Cantarero, Phys. Rev. B 86, 064305 (2012). [81] J.-S. Wang, X. Ni, and J.-W. Jiang, Phys. Rev. B 80, 224302 (2009). [82] L.-A. Wu and D. Segal, Phys. Rev. E 83, 051114 (2011).

032112-10