Electron-transport enhanced molecular dynamics for metals and semi-metals Reese E. Jones, Jeremy A. Templeton, Gregory J. Wagner, David Olmsted, and Normand Modine Abstract In this work we extend classical molecular dynamics by coupling it with an electron transport model known as the two-temperature model. This energy balance between free electrons and phonons was first proposed in 1956 by Kaganov and coworkers but has recently been utilized as a framework for coupling molecular dynamics to a continuum description of electron transport. Using finite element domain decomposition techniques from our previous work as a basis, we develop a coupling scheme that preserves energy and has local control of temperature and energy flux via a Gaussian isokinetic thermostat. Unlike previous work on this subject, we employ an efficient, implicit time integrator for the fast electron transport which enables larger stable time-steps than the explicit schemes commonly used. A number of example simulations are given that validate the method, including Joule heating of a copper nanowire and laser excitation of a suspended carbon nanotube with its ends embedded in a conducting substrate.

1

INTRODUCTION

In modeling non-equilibrium thermal transport in solids, classical molecular dynamics (MD) has the primary strength of explicitly representing phonon modes and the scatterers of phonons such as defects and boundaries. In addition, MD naturally represents the ballistic transport that can accompany reduced dimensionality, the diffusive transport expected at the macroscale and mixtures of the two [1–3]. On the other hand, free electrons and their role in thermal transport are missing. These electronic effects are vital in applications such as laser processing of materials [4–7], designing thermoelectric materials [8–11], estimating heat transport in conducting nanotubes and nanowires [8, 12], and in the performance of electronic devices in general [13]. Predictive models of the phenomenology of the interactions between the charge carriers and the atoms represented in MD exist. For instance, the two temperature model (TTM) [14–18] is essentially an energy balance between phonons and free electrons derived from the second velocity moment of the Boltzmann transport equations (see, e.g., G. Chen’s text [19, Chapter 8] for de1

tails). The TTM equations model the transport and relaxation of hot electrons as they interact with a relatively cold lattice. The use of an incoherent transport model for the carrier electrons, as opposed the coherent/wave-like transport described by ab initio quantum methods, can be justified on sufficiently large scales. For instance, the related drift-diffusion model (DDM) of electron transport is in wide-spread use in semiconductor device modeling [20]. In Section 2 we argue, by means of a Knudsen number based on the ratio of the electron mean free-path to the characteristic length-scale of the system, that the lower limit for these types of models can be on the order of the lattice parameter. As previous authors have pointed out, the analogy of the TTM, DDM and related electron transport models with hydrodynamics affords some of the same treatment, e.g. the use of a Reynolds-like number to define different regimes of transport [21], and the same intuition. Along those lines, in Section 2 we give a dimensional analysis of the timescales relevant to the TTM, as well as a brief exposition of its phenomenology using parameters for a typical metal, semi-metal and semiconductor. As a point of clarification regarding the use of a TTM description of semiconductors, note that there are other multiple temperature models [22–24] which split the thermal energy not between electrons and phonons but between acoustic and optical phonons; however, these models are not the subject of the present work. In addition to the continuum description of electron flow, another concern with using the TTM or any temperature-based model at the nanoscale is the existence of a local “equilibrium” temperature especially in modeling a highly non-equilibrium process. Alternate local, non-equilibrium definitions to the customary kinetic temperature exist [25–28]. However, the existence of an optimum non-equilibrium temperature measure is still a matter of debate and research [29]. Here we adopt a local kinetic temperature for the lattice and are careful to choose time- and length-scales such that its fluctuations are acceptable and the kinetic and potential energies are roughly in equipartition. Due to significant atomistic effects in many of the applications of interest, especially short-pulse laser heating of metals, the continuum TTM has recently been extended with MD coupling schemes to simulate the rapid exchange of thermal energy between electrons and the lattice [30, 31]. In this work, like these previous efforts, we take the TTM partial differential equations (PDEs) as a framework for coupling the electron and phonon mediated energy transport. Here, we use it to govern the evolution of the electron system represented on a finite element (FE) grid and the coupling between this system and the phonon system represented (in part) on an MD lattice. This development builds upon the recent work of Zhigilei and co-workers [30,32–35] and Xu and co-workers [31,36], by adding : (a) a rational development of the coupling algorithm through the principle of energy conservation, (b) an implicit subcycle time integrator for the electron transport that preserves energy, (c) a finite element-based domain decomposition for use in larger scale problems with local regions of interest, and (d) field-based thermostats capable of handling flux and temperature based constraints.1 Sections 3 and 4 give the details of these developments, including the 1 Apparently,

[37] uses version of domain decomposition with their finite difference scheme

2

derivation of a coupling method from energy conservation and Gauss’s principle of least constraint, as in [38]. In Section 5, a simple example simulation is used to validate the method. Next, the extension of the coupled TTM to applications that involving Joule heating of conductors is illustrated with a Cu nanowire. Here we are able to simulate the boundary layers and the steady state of the transport problem. By the use of the TTM as a platform for a full hydrodynamic description of the electron transport, this method differs from the ad hoc approach to Joule heating of MD lattices given in [39, 40]. This Joule heating example also shows our capability to control the end temperatures of the lattice while simultaneously exchanging heat with electron system over the remainder of the lattice. As a last example, we use the FE based domain decomposition to simulate the laser heating of a suspended metallic carbon nanotube (CNT) with its ends embedded in a substrate modeled by a continuum. Here, we employ nonlinear exchange laws not previously utilized in coupled TTM simulation. Finally, the paper is concluded in the last section with a brief summary of results and directions for future work.

2

TWO TEMPERATURE MODEL

In order to design the MD-FE coupling scheme and the time integration method we need to quantify the time- and length-scales of the TTM and understand its phenomenology. As mentioned in the Introduction, the two temperature model [14] is a set of equations ρe ce θ˙e = −∇ · qe − g + ρe re ρp cp θ˙p = −∇ · qp + g + ρp rp

(1)

derived from the second velocity moment of the Boltzmann transport equation for two types of interacting carriers, electrons (e) and phonons (p), with respective temperatures θe and θp . The symbol ρ denotes the respective mass densities and c the heat capacities per unit mass. The electron and phonon-mediated enˆ (θ, ∇θ) that are functions of ergy transport are modeled with heat fluxes q = q their respective temperatures and temperature gradients. The electron-phonon exchange rate g = gˆ(θe − θp ) accounts for the exchange of heat due to electronphonon collisions, which has the obvious symmetry gˆ(θe − θp ) = −ˆ g(θp − θe ). The terms ρe re and ρp rp represent external sources of heat, the latter of which, for instance, can be used to model injection of heat into the electron system by laser irradiation or Joule heating. Joule heating is the work done by the electric field on electron flow which appears explicitly in the derivation of the TTM [19, Chapter 8]. Due to this term the TTM by itself is not closed and typically requires additional balances of electron density and momentum which come from lower moments of the Boltzmann transport equation [20, 23]. Direct but the method is not described in detail.

3

sources in the phonon system are included here for completeness but not considered in the subsequent simulations. Together these equations (1) represent the balance of energy for a continuum transporting heat with two carriers but not undergoing significant deformation. The heat fluxes, qe and qp , result from the ensemble averages of the third moments of their respective velocities in the Boltzmann transport equation combined with the self-scattering terms. They are typically assumed to follow diffusive Fourier behavior, e.g. qp = −kp ∇θp where the conductivity tensor kp is simply kp = kp I for an isotropic material. Also the energy exchange term, which can be highly non-linear2 [35, 41, 42], can be approximated to first order as a linear function g ≈ ge-p (θe − θp ) of the temperature difference (θe − θp ). With these approximations and taking all the material properties to be constant and isotropic, a linear system of coupled, diffusive partial differential equations results ρe ce θ˙e = ke ∇2 θe − ge-p (θe − θp ) + ρe re (2) ρp cp θ˙p = kp ∇2 θp + ge-p (θe − θp ) + ρp rp In addition to the well-known diffusion solutions, see e.g. [43], where the two carriers are isolated by taking ge-p = 0, the system (2) admits other simple analytical solutions. First, a uniform exchange solution, with no spatial dependence, i.e. ∇2 θe = 0 and ∇2 θp = 0, and consequently no diffusion, is given by   t Θ ≡ θe − θp = Θt=0 exp − τe-p (3) ǫ ≡ ρe ce θe + ρp cp θp = ǫt=0 where τe-p =

ρe c e ρp c p . ge-p (ρe ce + ρp cp )

(4)

It describes the uniform relaxation of an initial temperature difference between the two systems such that the temperature difference follows a simple exponential decay with the total energy and its density ǫ remaining constant. Second, a steady-state (θ˙e = 0 and θ˙p = 0), quasi-one dimensional solution on a finite domain3 is given by       L L−x x Θ ≡ θe − θp = csch Θx=0 sinh + Θx=L sinh ℓe-p ℓe-p ℓe-p (5)
x ¯ ¯ ¯ ¯ θ ≡ ke θe + kp θp = θx=0 + θx=L − θx=0 L where s ke kp ℓe-p = . (6) ge-p (ke + kp ) 2 More

discussion of this subject is given in Section 5.3 and the Appendix. solution for the difference temperature (5a) on a semi-infinite domain can be recovered by taking the limit L → ∞. 3 The

4

Here, the temperature difference is again exponential and the conductivityweighted sum θ¯ is linear in space. This solution shows how boundary layers and non-uniform heating can occur in a conductor near reservoirs of “hot” electrons.

2.1

Time and length-scales

By examining three fundamental solutions of (2), namely isolated phonon diffusion, isolated electron diffusion, and uniform exchange, equation (3), three fundamental timescales can be identified: ρp c p ∗ 2 (7) L phonon diffusion τp = kp ρe c e ∗ 2 (8) L electron diffusion τe = ke ρe c e ρp c p electron-phonon exchange τe-p = (9) ge-p (ρe ce + ρp cp ) where L∗ is a characteristic length-scale in the specific problem. Note that in the atomistically-coupled problem a fourth timescale, that of atomic vibration, will become important since it determines the stable time-step of the MD evolution and is tied to the period of thermal fluctuations. This issue will be revisited in Section 4. In order to determine representative values of these timescales, we estimated the relevant properties for a metal, Cu; a semi-metal, graphite/CNT; and a semiconductor, Si, which are enumerated in Table 1 with the details of the calculation given in the Appendix. Though a temperature does not formally enter into the dimensional analysis, most of the parameters are temperature dependent so reference temperatures θe = 1000 K, i.e. moderately hot electrons, and θp = 300 K, i.e. room temperature phonons, were chosen. Clearly, the trend shown in Table 1 for the electronic heat capacity ρe ce , electron-mediated thermal conductivity ke and the electron-phonon coupling parameter ge-p is mainly due to the variation in the number of free electrons. The resulting characteristic time- and length-scales are given in Table q 2, where the length-scale ℓe-p was

e will be discussed in the following defined in (6) and the length-scale ℓe = gke-p section. In a small system, i.e. nano-devices, the electron diffusion τe is typically the fastest timescale. In larger systems, i.e. macroscopic, the exchange timescale τe-p is typically the fastest since it does not scale with L∗ . For nano-devices, the slowest timescale is typically phonon diffusion τp in metals but can be the exchange timescale τe-p in semiconductors where the free electron density is much smaller. All are much larger than the characteristic time-step of MD which is a fraction of an atomic vibrational period and usually less than a femtosecond. One assumption particular to MD used in these estimates is the Dulong-Petit law [44, Chapter 22] 3κB ρp c p = (10) Vα which represents the classical, high-temperature limit of heat capacity. Here Vα is the tributary volume of an atom, i.e. the volume of a unit cell divided by

5

Material Cu C Si

Debye temperature [K] 345 2500∗ 645

heat capacity [J/K-m3 ] ρe c e ρp c p 96600 3.52 × 106 1500 9.40 × 105 0.587 2.07 × 106

conductivity [J/K-m-s] ke kp 244.3 11.2 2.44 100∗ 0.0244 120

coupling coef. [J/K-m3 -s] ge-p 2.6 × 1017 9.0 × 1014 1.0 × 1011

Table 1: Electron and phonon material parameters for copper, graphite/carbon nanotube, and doped silicon, at θp = 300 K, θe = 1000 K. Note that values marked with an asterisk are not well-defined and therefore approximate.

Material Cu C Si

τp [ps] 31.4 0.940 1.72

τe [ps] 0.0395 0.0614 0.002040

τe-p [ps] 0.362 1.66 0.111

ℓe-p [nm] 6.41 51.5 494.

ℓe [nm] 30.7 52.1 494.

Table 2: Time-scales, exchange length-scale and dimensionless parameters for L∗ = 10.0 nm

the number of atoms per unit cell and κB is the Boltzmann constant. Another relevant assumption is the choice of an appropriate length-scale of resolution L∗ which affects both diffusion timescales. At least three separate physical criteria emerge in the coupled TTM: the exchange length, the lower limit of the diffusion description, and the volume size needed to reasonably define a local temperature. Unlike the latter two, the exchange length ℓe-p provides an upper limit on the characteristic length-scale if boundary layers are to be resolved. From a physical perspective, carrier relaxation/equilibration can impact device performance, and thus it is important numerically that the FE mesh be fine enough to resolve the relevant length scales while still being coarse enough so that the energy transported by the electrons can be accurately described by continuum equations. Here, Table 2 shows that ℓe-p ranges on the order of 5 − 50 nm (or a few to tens of lattice parameters in length) for metals and semi-metals. On the other-hand, the lower limit of the diffusion model can be estimated using electron mean free path data from [45] and a Knudsen number, defined as the ratio of L∗ to the mean free path of the carrier, must be larger than one.4 The data in [45] is for very high energy (temperature) electrons, the lowest values are 0.64 nm for Cu and 7.6 nm for Si at 200 eV (2.3 × 106 K), but the data is very linear with temperature. Extrapolating to room temperature, the mean free path for Cu is 0.42 nm and for Si 0.47 nm both of which are on the order of the lattice parameter for the respective materials. Lastly, in order to define a local quasi-equilibrium temperature such that the temperature 4 Discussion of the continuum approximation with respect to carrier density is postponed until the following section.

6

fluctuations are less than an order of magnitude of the mean, the volume given by L∗ 3 must contain at least a few hundred atoms in practice. For Cu there are approximately 80 atoms per cubic nm and for Si there are about 25 atoms per cubic nm. So, using these three competing criteria an appropriate length-scale of resolution can be estimated to be in the range 1 − 10 nm.

2.2

Dimensional and asymptotic analysis

In order to determine the scaling of boundary layers and transients we consider a non-dimensional version of the linear TTM equations (2): d τ τ τ (θe − θp ) + re θe = ∇2 θe − dt τe τg,e τr,e τ τ τ d (θe − θp ) + rp θp = ∇2 θp + dt τp τg,p τr,p

(11)

where a length-scale L is used to rescale x → L1 x and a timescale τ is used to rescale t → τ1 t. In addition, temperature has been rescaled by a reference temperature θ and the external energy sources by R. Clearly, ∇2 → L−2 ∇2 d d → τ1 dt . The timescales τe and τp have been defined previously in (7b) and dt and (7a), respectively. Note that both of these timescales are proportional to L2 and therefore are very sensitive to the characteristic size of the system. The ρ c ce new timescales are defined simply as: τg,e = ρgee-p , τg,p = gpe-pp , τr,e = creeθ , and c θ

τr,p = rpp . τ First we examine balance for electrons. Define parameter ǫ = τg,e , which is e q ke small as long as L ≫ ℓe = ge-p (see Table 2), and choose the timescale to be τ = 1ǫ τe . Then if we omit the external source for simplicity, we arrive at ǫ2

d θe = ǫ∇2 θe − (θe − θp ) dt

(12)

To first order, the two temperatures are in equilibrium θe = θp ; however, (12) is a singular perturbation problem. Since ǫ multiplies the highest order derivative, length needs to be rescaled near the boundaries to satisfy all the boundary conditions. For a boundary at x1 = 0, the rescaling x → ξ = { √1ǫ x1 , x2 , x3 } leads to 0 = ∇2ξ θe − (θe − θp ) (13)

Given non-trivial boundary conditions at ξ1 = 0 and the matching condition θe = θp√as ξ → ∞ in the far field, we get a boundary layer with characteristic length ǫL = ℓe-p . (Note unlike the solution (5) which concerned the difference in the two temperatures and defined ℓe-p , the length-scale ℓe is only germane to boundary layers in the electron temperature profile.) Likewise, rescaling time shows that transients in the electron system are significant at times on the order of ǫ2 τ = τg,e . Finally, we can expand θe as a perturbation series θe = θe,0 + ǫθe,1 + O(ǫ2 ) . 7

(14)

and substitute this series into (12) and collect terms. The next term after θe,0 = θp in the expansion is θe,1 = ∇2x θp ; so, to second order, θe = θp + ǫ∇2x θp + O(ǫ2 )

(15)

away from boundary layers and transients. Now examine phonon balance. Using same time scaling as in the electron transport analysis, d τe τe 2 ∇ θe + (θe − θp ) θp = dt ǫτp ǫτp,g is obtained. Substitution of the expansion (15) leads directly to   d τe τe ∇2 θp + O(ǫ2 ) + θp = dt ǫτp τg,p

(16)

(17)

Note that because typically τp ≫ τe , the first coefficient on the right hand side, τe ǫτp , is not necessarily a leading order term dominating the second coefficient τe τg,p . In fact, if we unscale (17), we get the dimensional equation ρp c p

d θp = (kp + ke ) ∇2 θp dt

(18)

and the result that the coupling with electrons effectively provides an enhanced conductivity for the phonons.

2.3

Metallic vs. semiconductor devices

As mentioned previously, the main difference between metals and semiconductors with respect to heat transport is the number of free electrons. On a phenomenological level, the high density and mobility of free electrons in metals leads to the electrons being the dominant carrier in energy transport. In this case, the strong electron-phonon exchange simply equilibrates the local energy between the electron and phonon systems. Also, the relatively small heat capacity of the electrons leads to the dominant balance in the electron energy transport equation being diffusion versus exchange in the absence of an external source. This also makes free transient solutions of the electron temperature somewhat trivial due to the rapid relaxation of the electrons to the phonon temperature. In semi-conductors, there are not obvious balances between pairs of dominant terms and the two systems are only weakly coupled due to the relatively small value of ge-p . With an external source r˜e applied to the electrons to balance diffusion, 0 = d θp dt

=

∇2 θe + (θe − θp ) + re

(19)

∇2 θp

(20)

8

we see that heating of the phonons due to exchange is a secondary effect, at least at large length-scales. However, in practice the exchange term will be weak until the source drives the electron temperature to a much higher value than the phonon temperature. In this case, the phonon temperature will essentially remain insulated from the electron temperature, but it will act as a drag on the electron temperature to balance the source. In this mode, the two temperature model serves to correctly compute the electron temperature by accounting for the phonon drag in a powered device. Given the computational requirements of resolving boundary layers and other features on the scale of ℓe-p with MD for semi-conductors, we will limit the remainder of this work to metallic materials with powered semi-conductors deferred to future work. It should be noted that another significant difficulty in modeling semi-conductors with an electronphonon temperature model is the free electron density is so low in these materials that it calls into question the continuum approximation. As one final note on employing the TTM to model electronic devices, metallic CNTs, given their high Debye temperature, have a non-linear electron-phonon exchange behavior g ∼ (θe −θp )5 that is characteristic of high temperature Allen theory [41]. This non-linearity, which is not manifest in noble metals and semiconductors, has a significant effect on the qualitative conclusions drawn from Table 2, especially at higher temperature differences. In fact, these materials can act like semiconductors at low temperature differences and like metals at high temperature differences.

3

MD-FE COUPLING

As mentioned in the Introduction, an atomistic material description is attractive relative to a continuum formulation for detailed representation of phonons, as well as modeling deformation and failure processes. Given the nature of MD, the electron system still needs a continuum method to model the transport of energy mediated by the flow of free electrons and thus there is a clear need to couple a FE representation of the missing physical processes to MD. The basic premise of the proposed coupling method is (a) to replace the phonon energy transport equation (1b) in regions of interest Ωmd with MD and keep the phonon transport equation on the remaining domain Ω \ Ωmd as an efficient surrogate, and (b) use the electron energy transport equation (1b) to model the electron transport everywhere in Ω. To effect this we embed the MD lattice occupying Ωmd ⊂ Ω in a FE grid covering Ω, see the schematic in Figure 1a, which will represent the continuum transport equations (1). Note that in the following I enumerates nodes, and α enumerates atoms. The coupling between the MD and FE regions we propose is predicated on the principle of energy conservation. The total energy should be conserved for a closed system or, more generally, the time rate of change of the total system energy should be equal to the change due to external sources: E˙ ext = E˙ md + E˙ fe

9

(21)

(a) top view

Ω Ωmd

(b) side view, grids separated for clarity Θe g q.n

q.n

Θp

Figure 1: Finite element (FE) and molecular dynamics (MD) domains, top figure (a). The bottom figure (b) shows heat transfer between the electron and phonon systems, upper and lower grid respectively, and between the FE and MD representations of the phonon system, not shaded and shaded regions respectively.

The left hand side of (21) is the sum of all external sources, namely Z Z (n · qe + n · qp ) dA (ρe re + ρp rp ) dV + E˙ ext =

(22)

∂Ω

Ω

where n is the outward unit normal (to the particular region of integration) and the boundary fluxes n · qe and n · qp are computed over the entire boundary. On the subsets of the boundary Γqp and Γqe these fluxes are prescribed and on the remainder of the respective boundaries Dirichlet boundary conditions are applied to the individual fields θp and θe . The MD component Pof the total energy change is simply the sum of the changes in the kinetic K = α 21 mα vα · vα and potential Φ = Φ({xα }) energies X
(23) mα vα · v˙ α − fαmd · vα E˙ md = K˙ + Φ˙ = α

where the interatomic force is given by fαmd = −∂xα Φ and xα , vα are atomic position and velocity, respectively. In general, this term is non-zero due to the presence of a thermostat (which will be discussed later in this section).5 The 5 In [38, Section 5], E md is time averaged to minimize the effects of thermal fluctuations on the empirical temperature. Here, for simplicity, we omit this complication, leaving a discussion of its applicability to the end of Section 4.

10

FE component of the energy change is Z Z fe ˙ ˙ ρp cp θ˙p dV E = ρe ce θe dV + Ω\Ωmd Ω Z Z Z Z (ρp rp + g) dV (24) ∇ · qp dV + (ρe re − g) dV + ∇ · qe dV + = Ω\Ωmd Ω Ω\Ωmd Ω Z Z Z Z Z g dV ρp rp dV − ρe re dV + n · qp dA + n · qe dA + = ∂Ω

Ωmd

Ω\Ωmd

Ω

∂(Ω\Ωmd )

where the phonon system is represented only on Ω\Ωmd and the electron system over the whole domain Ω. Note that the exchange terms cancel over the FE domain, Ω \ Ωmd . A reduced form of the total balance (21) Z Z X
(ρp rp + g) dV (25) n · qp dA − mα vα · v˙ α − fαmd · vα − 0= α

∂(Ωmd )

Ωmd

results from substitution of (22), (23), and (24) into (21). A local energy balance can be postulated on the FE length-scale using a partition of unity NI E˙ I = E˙ Imd + E˙ Ife − E˙ Iext = 0 Z X
md ˙ = NIα vα · mα vα − fα − α

∂Ωmd

NI n · qp dA −

Z

NI (ρp rp + g) dV

Ωmd

(26)

where, in particular, NI is chosen to be the usual FE shape function at node I. With compact support only a subset of the total number of atoms NA contribute to EI and hence, as discussed in Section 2, the FE length-scale is bounded from below by a simulation dependent length-scale at which a local temperature and local thermodynamic equilibrium are surmised to exist, as manifest by energy fluctuations and mean trends. This localized energy balance (26) is a constraint that leads to both a thermostat on the MD lattice and FE temperature evolution equations. There are a wide variety of thermostats available to control the temperature of a system, most notably the commonly used Nos´e-Hoover [46] and Langevin thermostats [47]. Here we are concerned with local and instantaneous control of non-equilibrium processes, and in that mode we adopt a local variant of the Gaussian isokinetic thermostat [48] first proposed in [38]. This treatment is similar to that in [32] and in contrast to the Langevin algorithm employed in [49]. As pointed out in [38], the instantaneous action of the Gaussian isokinetic thermostat has the advantage over the deterministic Nos´e-Hoover and stochastic Langevin thermostats in non-equilibrium applications since both of the latter thermostats have an implicit timescale that gets convolved with the physical time-scales of the dynamic process. This advantage comes at the expense of not generating the canonical ensemble in systems smaller than the thermodynamic limit that are truly in equilibrium. As we will see, the Gaussian isokinetic

11

thermostat also provides a compatible framework for incorporating control of heat flux in the manner of Ikeshoji and Hafskjold [50] in specified regions. In the MD region, Gauss’s principle of least constraint can be written as X X min mα (v˙ α − v˙ ∗ )2 + (27) λI h˙ I v˙ α

α

α

I

∗ and provides a variational framework for the thermostat control. In (27), v˙ α = 1 md are the unconstrained accelerations and λI are the Lagrange multipliers mα fα used to enforce the coupling constraints. In our formulation, the nodal nonholonomic constraint h˙ I is the localized energy balance h˙ I = E˙ I = 0, equation (26), or a local temperature constraint T˙I (vα ) = (θ˙p )I between the rate of change of the temperature TI , determined by the MD state, and the empirical temperature (θp )I , scaled by an appropriate dimensional factor. Specifically, the temperature constraint takes the form: ! X    X 1 3κ B h˙ I = ρp cp T˙I − (θ˙p )I NIα Vα = T˙I − (θ˙p )I = 0 (28) NIα 2 2 α α

Here we have employed the Dulong-Petit relation (10) and a representative volume given by the sum of the shape function NIα = NI (xα ) evaluated at the atomic position xα times the atomic tributary volume Vα . Up to this point we have not given a prescription for the atomic temperature TI and, here, the fact that Gauss’s principle relies on the variation of the time derivatives of the atomic velocities necessitates a kinetic definition. As in Wagner et al. [38, Eqns. (9-13)], we use a local kinetic temperature defined by 1 X ˜ NIα mα vα · vα (29) TI = 3κB α ˜Iα = Here, N

P 1 NIα β NIβ

is a normalized version of the shape function at node

I that results from a lumped least squares minimization described in [38]. As a consequence of the form of both types of constraints, the associated EulerLagrange equations are 1X ∗ λI NIα mα vα = 0 (30) mα (v˙ α − v˙ α )+ 2 I

and h˙ I = 0. A field-based Gaussian isokinetic thermostat follows:6 P v˙ α = m1α fαmd − 12 I NIα λI vα    P 3κB ˙ md  ( θ ) f − N p I  α α Iα 2 P R  1 − RdΩmd NI n · qp dA α,J NIα (mα vα · vα )NJα λJ =  2  − Ωmd NI (ρp rp + g) dV

(31) for I ∈ Nθp

(32) for I ∈ N \ Nθp

6 Note that in the case of a flux-based control, the right-hand side of Eqn. (32) and λ could be split into their exchange and FE-MD conduction parts and each solved separately. This could be advantageous since the two components act on different time-scales.

12

where Nθp is a set of nodes where temperatures are prescribed and on the remainder the energy balance E˙ I = 0 is enforced. The per-atom drag P force customarily employed in kinetic thermostats can be defined as fαλ = 21 I NIα λI vα . Now we turn to the governing equations of the FE electron and phonon temperatures, which are Z Z XZ
NI q¯e dA ∇NI · qe + NI (ρe re − g) dV − NI ρe ce NJ dV (θ˙e )J = Ω

J

XZ J

Γqe

Ω

NI ρp cp NJ dV (θ˙p )J =2

Ω

X α

+

Z

NIα vα · fαmd + fαλ

Ω\Ωmd





∇NI · qp + NI (ρp rp + g) dV −

Z

Γqp

NI q¯p dA (33)

and involve the drag force fαλ and all the sources of heat external to the MD system. Here, q¯e and q¯p represent prescribed surface fluxes of heat. Equation (33a) is simply the Galerkin projection of (1a). Equation (33b) is a Galerkin projection of (1b) on Ω \ ΩMD with the substitution of Z X X
(34) NI ρp cp θ˙p dV ≈ 2 NIα mα vα · v˙ α = 2 NIα vα · fαmd + fαλ Ωmd

α

α

as in [38, Section 4].7 Unlike in our previous work [38], in (33b) we explicitly remove all internal contribution due to MD and do not rely on the “atomic” quadrature implicit in (34) canceling the native quadrature in elements fullyfilled with the lattice. In the special case where the system is isolated and Ωmd = Ω, the equations (33) and (31) reduce to Z XZ
∇NI · qe + NI (ρe re − g) dV NI ρe ce NJ dV (θ˙e )J = J

and

Ω

Ω

J

XZ

NI ρp cp NJ dV (θ˙p )J =2

Ω

X α

(35)

NIα vα · fαmd + fαλ

Z 3κB X NI g dV NIα Tα NJα λJ = 2 Ω


(36)

α,J

where Tα = 3κ1B mα vα · vα . Since equation (35b) is essentially a time derivative of the definition of kinetic temperature (29) given the assumption the classical Dulong-Petit limit (10), it can be replaced with the simple equation (θp )I = TI . This reduced system, with a lumped solution of (36), is similar to that employed in [31, 32]. 7 A further reduction was made in [38, Eqn. 29] by canceling the boundary flux term with half the component of the thermostat power term due to conduction, which is possible in the absence of electron-phonon exchange which occurs across the whole domain.

13

In summary, Figure 1b illustrates the flow of heat between the two systems of carriers due to exchange and between the MD lattice and FE mesh due to boundary conduction. The method is comprised of : (a) the control of the MD system, equations (31) and (32) , (b) the evolution equations for the FE fields (33), and (c) the time integration scheme, which will be described in the next Section.

4

TIME INTEGRATION

As Section 2 made clear, the electron balance occurs at much faster timescales that the corresponding phonon transport. This widely recognized fact motivates the design of time integration schemes for TTM-like PDEs such that the electron system is integrated with a smaller time-step that the phonon system, see, for example, the so-called “gallop-trot” method [22]. In addition to the disparity in timescales, the expense of calculating the interatomic forces fαmd relative to the operations involved in the FE-based equations motivates special treatment of the fast electron diffusion timescale. Specifically, the proposed scheme allows the MD system to evolve with its normal time-step based on the vibrational frequencies of the atomic lattice. Another guiding principle in constructing the integration scheme was the re-use of existing MD integration schemes and previous coupling method [38]. After spatial discretization of the continuum part of the coupled TTM, but before temporal discretization, we have three sets of coupled ordinary differential equations (ODEs), (31) and (33), as well as an algebraic equation, (32): θ˙ e = Fe (θ e , θp ) θ˙ p = Fp (θe , θp ) + Fmd (x, v, λ)
1 v˙ = fmd (x) + fλ (v, λ) m λ = λ(x, v, θ˙ p )

(37)

where bold symbols have been introduced for vectors such that x, v, fmd , fλ contain the atomic state for atoms α = 1, ..., NA (position, velocity, force due to the potential and force due to the thermostat, respectively) and θe , θp , λ contain the nodal information ( empirical electron and phonon temperatures, and drag coefficient, in order) for nodes I = 1, ..., NN . The function Fe represents the right-hand side of (33a), Fp represents the FE part of the right-hand side of (33b), and Fmd represents the MD part of the right-hand side of (33b). First, we solve (37d) for λ, construct fλ and use a Verlet half step 1 ˜ n+1/2 = vn + ∆t fλn v 2

(38)

to update the atomic velocities. Second, a third-order Gear [51] predictor is

14

used update the phonon temperature ˜ n+1 =θn+1/2 + ∆t θ˙ n + 1 ∆t2 θ ¨n θ p p p p 2 n n ˜θ˙ n+1 = ¨ θ˙ + ∆t θ p

p

¨ ˜θn+1 = p

p n ¨ θp

(39)

As in [38], the integration scheme for the empirical phonon temperature was deliberately chosen to be higher order than the MD integrator so as not to degrade its accuracy. Third, we integrate electron transport (37a) using a subcycled variant of an implicit Euler scheme m n+i/m ˜ n+1 ) + (1 − α)Fe (θ n+(i−1)/m , θ ˜ n+1 ) (θe − θn+(i−1)/m ) = αFe (θn+i/m ,θ e e p e p ∆t (40) with i the sub-cycle index, m the number of sub-cycles and a parameter α ∈ [0, 1]. Using a Beam and Warming scheme outlined in [52, 53], we take a Taylor n+i/m series expansion of Fe and truncate at first order in ∆θe = (θe )n+i/m − (θe )n+(i−1)/m to obtain Fn+i/m = Fn+(i−1)/m + ∂θ e Fe ∆θe + O(∆θ 2p ) e e then substitute this into (40) :   ) 1 − α∆t ∂θ e Fe ∆θ e = ∆t Fe (θn+(i−1)/m e

(41)

(42)

We solve this system using a finite difference approximation of ∂θ e Fe ∆θe with the generalized minimum residual algorithm (GMRES) [54]. We employed a backward Euler integrator, α = 1, as opposed to the higher order accurate trapezoidal rule, α = 1/2, since that the latter produced oscillatory albeit stable solutions. A major advantage of the implicit scheme over the commonly used explicit sub-cycled algorithm is that the implicit integrator can produce solutions that approach the quasi-static limit even when the electron heat capacity is relatively negligible. Also, relative to explicit Euler, an implicit Euler scheme allows for larger stable time-steps in the presence of nonlinearity, with the usual trade-off with accuracy or computing cost. Fourth, the usual Verlet integration scheme is applied to the MD data vn+1/2

=

xn+1

=

˜ n+1 v

=

1 ˜ n+1/2 + ∆t fmd (xn ) v 2 xn + ∆t vn+1/2 1 vn+1/2 + ∆t fmd (xn+1 ) 2

(43) (44) (45)

Fifth, we solve (37d) again for λ, construct fλ and use a Verlet half step 1 ˜ n+1 + ∆t fλn+1 vn+1 = v 2 15

(46)

to update the atomic velocities. Sixth and last, the Gear corrector for the phonon temperature evolution is applied: n+1 ˜ n+1 , θn+1 ) + Fmd (xn+1 , vn+1 , λn+1 ) θ˙ p = Fp (θ p e ˙ n+1 ˙ n+1 n+1 ˜ R = θp − θp ˜ n+1 + 5 Rn+1 θ n+1 =θ p p 12 n+1 1 ˙ n+1 ¨ ˜ Rn+1 θ p = θp + ∆t2

(47)

It is possible to employ the same causal time filter with an exponential kernel described in [38, Section 5] to the atomic temperatures resulting from (23) to offset the fluctuations caused by elements with relatively small support and therefore few atoms to average over. However, unlike [38] where the characteristic time of the filter was only limited by the phonon relaxation time, the fastest time-scales that needs to be resolved with the TTM are those associated with the electron processes. By the use of an another perturbation solution, we can determine how fast fluctuations on the scale of atomic vibration affects the phonon and electron temperature dynamics. In this analysis we are only interested in the effects of coupling on the electron evolution, so consider (12) but ignore diffusion and the external source. In this case, choose τ = τg,e to obtain the non-dimensional equation for the electron temperature d θe = − (θe − θp ) dt

(48)

and note, in the present context, θp is also varying on a fast timescale, τv . To treat this case, we perform a multi-timescale analysis by defining a small parameter ǫ = ττv and second non-dimensional timescale t2 = 1ǫ t. We assume both temperatures depend on both timescales, θe (t, t2 ) and θp (t, t2 ), and clearly ∂ ∂ d dt → ∂t + ǫ ∂t2 . Similar to the analysis in Section 2.2, we expand the electron temperature in powers of ǫ θe (t, t2 ) = θe,0 (t, t2 ) + ǫθe,1 (t, t2 ) + O(ǫ2 ) then substitute this expansion into (48) and collect terms. The order result in ∂ θe,0 = 0 ∂t2

(49) 1 ǫ

terms (50)

which has the solution θe,0 = θe,0 (t). The order ǫ terms give ∂ ∂ θe,1 = −θe,0 − θe,0 + θp ∂t2 ∂t which has the solution   Z t2 ∂ θp (t, s) ds + C(t) θe,1 = −θe,0 − θe,0 t2 + ∂t 0 16

(51)

(52)

where C(t) is a function of t only. Since the first two quantities on right hand side are secular terms that increase without bound as t2 gets large, stability requires that    Z t2 1 ∂ lim (53) θp (t, s) ds = 0 −θe,0 − θe,0 t2 + t2 →∞ t2 ∂t 0 Rt Defining the time average hθp i = limt2 →∞ t12 0 2 θp (t, s) ds leads to recasting (48) as ∂ θe,0 = − (θe,0 − hθp i) . (54) ∂t This implies that, to first order, the equation for the electron temperature can be taken as the original (48) but using a time filtered value for θp . In practice we use an exponential kernel in this causal time average to give it the property of fading memory [38, Section 5]. Given this analysis, filtering the phonon temperature over timescales less than the intrinsic timescales of the phonon system, i.e. τp and τg,p , which are typically considerably longer than that of the electron system, has utility in the context of the TTM and should not significantly affect the dynamics of the electron transport. An example of the application of this technique will be given in Section 5.3.

5

SIMULATIONS

In order to test the method and illustrate its applications we constructed three model problems using a variety of material systems. We assumed that the empirical potentials were insensitive to variations in local free electron density and electron temperature over the range we simulated, and consequently used published values for their parameterization. The details of how the remaining empirical parameters were estimated are given in the Appendix.

5.1

Relaxation of a Gaussian profile

To verify the methodology we simulated a fictitious material based on LennardJones (LJ) Argon with constant thermal properties ce-p = ρe ce /ρp cp = 1/2, ke /kp = 8 so that de-p = 16 and τe /τp = 1/16. The exchange coefficient ge-p was chosen such that τe-p /τp = 3/4. A film of material 24 unit cells square and 6 unit cells deep, see Figure 2, was given an initial (Gaussian) electronic temperature distribution θe (x) = T0 (3 exp(−(x21 + x22 )/252 ) + 5) and a uniform phonon temperature θp = T0 , where T0 = 20 K. To attain this initial condition, the atomic lattice was thermalized using the local rescale thermostat described in [38, Section 6.4]. Only the center 12 × 12 unit cells of the domain contained “real” atoms but sufficient fixed ghosts were used outside this region to insulate these atoms from free surface effects. These zero temperature ghost layers do not affect the temperature of the atoms undergoing dynamics other than providing soft walls for phonon reflection. The FE grid of 12 × 12 × 1 elements covered 17

the whole domain of 12 × 12 × 6 unit cells and both the lattice and the mesh were prescribed to be periodic in the third dimension. The sequence of temperature profiles in Figure 3 shows the relaxation of the peaked electronic temperature profile to equilibrium through a combination of diffusion and exchange. Note that the atomic temperature at the boundary of the MD region does not exactly track the empirical phonon temperature at every instant due to the nature of the energy preserving flux constraint at the MD-FE boundary. Instead, given that the system is closed we expect overall energy conservation. Figure 4 shows the total energy is essentially constant while the two systems exchange energy. Looking at the difference of the average temperatures, Z 1 θ¯e-p = θ¯e − θ¯p = (θe − θp ) dV , (55) V

in Figure 5 we see the equilibrium state is in accordance with the ratio of the ρ c θ +ρ c θ heat capacities θ¯ss = e ρee cee +ρpp cpp p . Referring to the uniform exchange solution (3), the the average difference temperature θ¯e-p is governed by 1 ¯ θe-p θ¯˙e-p = τe-p

(56)

which has the exact solution θ¯e-p (t) = θ¯e-p (0) exp(−t/τe-p ). Figure 5 shows that the decay is exponential, consistent with this solution, up to the limit where MD driven fluctuations dominate. Similarly, the standard deviation of the electron temperature distribution is related to its diffusion constant, but a quantitative validation is confounded by the finite domain effects. Nevertheless, Figure 6 shows in a qualitative fashion that variance of the electron temperature profile decrease more rapidly than the average temperatures converge, which is in accordance with the relative size of the electron diffusion and exchange time-scales. In real materials these time-scales can be separated by orders of magnitude, thus reinforcing the fact that free transients in the electron temperature typically have little effect on the phonons.

5.2

Joule heating of a nanowire

In this example we create a Cu nanowire 360 unit cells (130.68 [nm]) long with cross-section 4 lattice units square using the EAM potential from [55]. Unlike the first example, here we employ the usual temperature dependencies for ce = γθe , with γ = 96.6 [J/m3 K 2 ] and ke = Lσθe , with Lσ = 107 [1/Ωm]. Like the model in [56], we use an exchange law that is linear in temperature difference with g = 2.6 × 1017 [W/m3 K]. A FE grid 40 elements long and one in cross-section covers the Cu lattice. Since there are no elements lacking atoms and therefore atomic information, there is no need to specify the phonon properties. We did prescribe the electron and phonon temperatures at the ends of the nanowire, θe = 400 K and θp = 300 K, to be distinct and constant due to contact with large reservoirs. Simple Joule heating was simulated by assuming an uniform external source to electrons ρe re = −∇Ψ · Je = 3.2 × 10−4 [W/nm3 ], i.e. the 18

Figure 2: Configuration of the Ar lattice (colored by the phonon temperature) embedded in a FE grid (colored by the electron temperature). Note the 5 layers of ghosts (blue color) are fixed and therefore have zero kinetic temperature.

electric field −∇Ψ and current Je are assumed to be uniform. Furthermore, we assume this source turns on and off instantly due to lower relaxation time for electrons’ momentum versus their energy. The system is heated for 47.5 ps with this uniform external source and then allowed to relax. Figure 7 shows the average temperature vs time. As expected from the analysis of Section 2.2, the response of the electron system to the source is nearly instantaneous but sufficiently resolved by our integration algorithm over the first 0.25/ps. After this transient, electron-phonon exchange becomes dominant and the temperature rate of rise is nearly equal for the two systems. Once the source is turned off, the electron system experiences another rapid transient. Afterward, the two system average temperatures approach each other and their respective boundary temperatures on the timescale of exchange. Figure 8 shows temperature profiles during the heating phase and the relaxation phase, as well as the non-equilibrium steady state where the boundary layer is clearly resolved. In this case the kinetic temperature defined by the atomic velocities is identical to the empirical phonon temperature at all nodes. During heating both the electron and phonon temperature profiles are basically parabolic as expected. Steady state occurs when the conduction into the electron system from the boundaries balances the losses to the boundaries into the phonon system, which are both equal to the total exchange between the two systems.

19

8

electron phonon atomic

7

TEMPERATURE

6 5 4 3 2 1 0 -0.6

-0.4

-0.2 0 0.2 X-COORDINATE

0.4

0.6

Figure 3: Sequence of temperature profiles through a cross-section aligned with the x-axis. The black arrows indicate the progression of profiles with time.

5.3

Laser heating of a carbon nanotube

In this example a metallic (8,8) armchair CNT, 12.6 nm long, is suspended in air by embedding its ends in solid graphite Figure 9. The CNT is modeled with the Tersoff potential [57, 58] and the graphite substrate is modeled with a continuum with the same thermal properties as the CNT. The exposed surface of the reservoirs and the tube are insulated by the air so that no heat crosses those boundaries and the remaining surfaces of the reservoirs are fixed at a constant temperature of 300 K. As in the previous example, the temperature dependency of the heat capacity is ce = γθe with γ = 1.5 [J/m3 K 2 ] and electronic heat conduction is ke = Lσθe with Lσ = 2.443 × 10−3 [W/mK 2 ]. In addition to temperature dependence of the electron heat capacity and conductivity, the measured form of the exchange for CNTs [42] is highly non-linear in temperature, g = h(θe − θp )5 with h = 3.7 × 104 [W/m3 K 5 ]. The fact that the CNT lattice is not space filling creates no algorithmic difficulties since its contribution directly and fully determines the phonon temperature in the regions where there are atoms. In these elements all the effects of the phonon constitutive model are removed, as mentioned in Section 3. Also, to offset the larger fluctuations associated with shape functions with few atoms in their support we employ a time-filter, with characteristic time-scale τ = 0.01 ps, discussed in Section 4 and in [38, Section 5]. The electron system of the CNT interacts with a focused radiation source that has a power input of ρe re = 1.6 × 10−12[W/m3 ] exp(−(x21 + x22 )/(0.1[nm])2 ) which translates to 7.4[J/ps] exp(−(x21 + x22 )/(0.1[nm])2 ) per atom given a tributary volume Vα = (0.36)3 [nm3 ]. We turn on this localized source for 50 ps and then allow the system to relax. The sequence of temperature profiles along the

20

0.3 0.25

ENERGY

0.2

electron phonon atomic total

0.15 0.1 0.05 0 0

2

4

6

8

10

12

14

16

18

TIME

Figure 4: Evolution of energy. The atomic thermal energy is calculated for the lattice, phonon energy is calculated in Ω \ Ωmd and the electron energy is integrated over the full domain Ω.

axis of the tube in Figure 10 shows very localized electron and relatively diffuse phonon temperatures in correspondence with their diffusivities. These profiles through the axis of the CNT extend into FE regions without atoms. In the reservoir regions we see a distinct change in slope due to the reservoirs’ higher thermal mass, especially for the phonons. As the experiments [1] demonstrate, we expect mixed ballistic/diffusive transport in the CNT, which is modeled entirely by the MD. This mixed harmonic/enharmonic transport which must translate to purely diffusive heat flux at the CNT-reservoir boundary, given the nature of the coupling. The large scale oscillations that start to become apparent at about t = 40 ps in Figure 11 indicate that the input energy to electrons eventually excites a strong fundamental mode resonance which can be directly observed in Figure 12. Note that the fact that the lattice is not coincident with the mesh is handled simply by using Lagrangian map of atoms to elements, as in all the previous examples.

6

CONCLUSION

The proposed work improves upon the few existing efforts to add electronic thermal transport effects to molecular dynamics by fully coupling the two temperature model to the lattice using domain decomposition in a rational framework based on the principle of energy conservation and Gauss’s principle. We also extend the regimes of application of coupled TTM algorithms with new electronphonon exchange models and the resolution of boundary layers. By coupling classical MD to FE-based models of the missing physics, we enable the simula-

21

DIFFERENCE AVERAGE TEMPERATURE

5.5 AVERAGE TEMPERATURE

5 4.5 4 3.5

electron phonon steady solution

3 2.5 2 1.5 1 0

2

4

6

8

10

12

14

16

18

10 1 0.1 0.01 0.001 data solution

1e-04 1e-05 1e-06 1e-07 1e-08 1e-09 0

TIME

2

4

6

8

10

12

14

16

18

TIME

Figure 5: Evolution of average temperature (left) and logarithm of the difference in the average temperature (right).

tion of a broad range of physical phenomena from the rapid exchange of heat between the electron and phonon carriers in a lattice through current-induced thermal failure of nanowires. Furthermore, the tight coupling between the MD and FE paradigms efficiently utilizes the inherent strengths of each in modeling energy and charge transport, as well as deformation processes. To enhance the current algorithm to fully account for thermo-mechanical processes we will need to split the velocity of the atoms into a coarse scale representable on the FE grid and a fine scale. The coarse velocity will participate in a continuum momentum balance as well as contribute to the stress power term missing in (1b). This coarse scale velocity would also be subtracted from the atomic velocities to obtain the fine scale velocities upon which the kinetic temperature (29) would be based. In addition, a means for handling the simultaneous transmission of heat and momentum will need to be adopted, see, e.g., [59, 60]. As discussed in Section 2, one issue with using the TTM as a model of electron-phonon energy transport is that the TTM is not closed. It requires other hydrodynamic balances or a constitutive law for the power term associated with the work of the electrochemical potential on the flow electrons. In this paper we assumed a reasonable form for this term in a simple Joule heating example simulation. In the sequel, we intend to use the electron momentum balance and continuity to resolve the local variations in electron density and current. Ultimately, the application of these methods to the simulation of devices with nanoscale features such as defects and material boundaries should further the development and performance of this important technology.

Acknowledgments This work was funded by the Laboratory Directed Research and Development program at Sandia National Laboratories and its support is gratefully acknowledged. Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United States Department of Energy under

22

10 TEMPERATURE VARIANCE

1 0.1 0.01 0.001 1e-04 1e-05 1e-06

electron phonon

1e-07 0

2

4

6

8

10

12

14

16

18

TIME

Figure 6: Evolution of standard deviation of the electron and phonon temperatures over the full domain Ω.

contract DE-ACO4-94AL85000.

References [1] Wang J, Wang JS. Carbon nanotube thermal transport: Ballistic to diffusive. Appl. Phys. Lett. 2006; 88(11):111 909. [2] Chen R, Hochbaum AI, Murphy P, Moore J, Yang P, Majumdar A. Thermal conductance of thin silicon nanowires. Phys. Rev. Lett. 2008; 101(10):105 501. [3] Chang C, Okawa D, Garcia H, Majumdar A, Zettl A. Breakdown of Fourier’s law in nanotube thermal conductors. Phys. Rev. Lett. 2008; 101(7):075 903. [4] Chrisey DB, Hubler GK. Pulsed Laser Deposition of Thin Films. Wiley Interscience: New York, 1994. [5] Baerle D. Laser Processing and Chemistry. Springer: Berlin, 2000. [6] Sundaram SK, Mazur E. Inducing and probing non-thermal transitions in semiconductors using femtosecond laser pulses. Nature Matls. 2002; 1(4):217–224. [7] Wang Y, Xu X, Zheng L. Molecular dynamics simulation of ultrafast laser ablation of fused silica film. Appl. Phys. A 2008; 92(4):849–852.

23

550 AVERAGE TEMPERATURE [K]

AVERAGE TEMPERATURE [K]

1000 900 800 700 600 500 400

ELECTRON TEMPERATURE PHONON TEMPERATURE

300 0

20

40

60 TIME [ps]

80

500

450 ELECTRON TEMPERATURE PHONON TEMPERATURE 400

350

300

100

120

0

0.5

1 TIME [ps]

1.5

2

Figure 7: Evolution of temperature and electron temperature for the full Cu system, full evolution (left) and close-up of initial transient (right). 1200

1200

ELECTRON TEMPERATURE PHONON TEMPERATURE MD TEMPERATURE

1100

1000 TEMPERATURE [K]

1000 TEMPERATURE [K]

ELECTRON TEMPERATURE PHONON TEMPERATURE MD TEMPERATURE

1100

900 800 700 600

900 800 700 600

500

500

400

400

300

300 0

20

40 60 80 100 X-COORDINATE [nm]

120

140

steady state

0

20

40 60 80 100 X-COORDINATE [nm]

120

140

Figure 8: Cu temperature profiles for heating phase (left) and relaxation phase (right). The black arrows indicate the progression of time.

[8] Kim P, Shi L, Majumdar A, McEuen P. Thermal transport measurements of individual multiwalled nanotubes. Phys. Rev. Lett. 2001; 87(21):2155 021– 2155 024. [9] Yang B, Liu WL, Liu JL, Wang KL, Chen G. Measurements of anisotropic thermoelectric properties in superlattices. Appl. Phys. Lett. 2002; 81(19):3588–3590. [10] Li D, Wu Y, Fan R, Yang P, Majumdar A. Thermal conductivity of si/sige superlattice nanowires. Appl. Phys. Lett. 2003; 83(15):3186–8. [11] Snyder GJ, Toberer ES. Complex thermoelectric materials. Nature Matls. 2008; 7(2):105–114. [12] Ou M, Yang T, Harutyunyan S, Chen Y, Chen C, Lai S. Electrical and thermal transport in single nickel nanowire. Appl. Phys. Lett. 2008; 92(6):063 101–1–3.

24

1800

18000

1600

16000

ELECTRON TEMPERATURE [K]

TEMPERATURE [K]

Figure 9: Metallic CNT embedded in a FE mesh showing phonon temperature near the beginning of the heating phase, t = 10 ps

1400 1200 1000 800 600 400 200

14000 12000 10000 8000 6000 4000 2000 0

0

5

10

15 20 25 30 X-COORDINATE [nm]

35

40

0

5

10

15 20 25 30 X-COORDINATE [nm]

35

40

Figure 10: Sequence of temperature (left) and electron temperature (right) profiles along the axis of the CNT.

[13] Majumdar A, Reddy P. Role of electron-phonon coupling in thermal conductance of metal-nonmetal interfaces. Appl. Phys. Lett. 2004; 84(23):4768–4770. [14] Kaganov MI, Lifshits IM, Tanatarov LV. Relaxation between electrons and the crystal lattice. Zhurnal Eksperimental’noi i Teoreticheskoi Fiziki 1956; 31(2(8)):232 – 237. [15] Anisimov SI, Kapeliov BL, Perelman TL. Emission of electrons from the surface of metals induced by ultrashort laser pulses. Zhurnal Eksperimental’noi i Teoreticheskoi Fiziki 1974; 66(2):776 – 81. [16] Qiu TQ, Tien CL. Heat transfer mechanisms during short-pulse laser heating of metals. J. Heat Transfer 1993; 115(4):835–841. [17] Melikyan A, Minassian H, A Guerra I, Wu W. On the theory of relaxation of electrons excited by femtosecond laser pulses in thin metallic films. Appl. Phys. B 1999; B68(3):411 – 14.

25

1400

TEMPERATURE

LEFT RESERVOIR TUBE 1200 RIGHT RESERVOIR 1000 800 600 400 200 0

10

20

30

40

50

60

70

80

90

100

TIME [ps]

Figure 11: Evolution of average temperatures of the two explicitly modeled reservoirs and the CNT.

[18] Jiang L, Tsai HL. Improved two-temperature model and its application in ultrashort laser heating of metal films. J. Heat Transfer 2005; 127(10):1167 – 1173. [19] Chen G. Nanoscale energy transport and conversion. Oxford: Oxford,UK, 2005. [20] Ghione G, Benvenuti A. Discretization schemes for high-frequency semiconductor device models. IEEE Trans. on Antennas and Propagation 1997; 45(3):443 – 56. [21] Majumdar A, Fushinobu K, Hijikata K. Effect of gate voltage on hotelectron and hot-phonon interaction and transport in a submicrometer transistor. J. Appl. Phys. 1995; 77(12):6686–94. [22] Fushinobu K, Majumdar A, Hijikata K. Heat generation and transport in submicron semiconductor devices. ASME J. Heat Transfer 1995; 117(1):25–31. [23] Lai J, Majumdar A. Concurrent thermal and electrical modeling of submicrometer silicon devices. J. Appl. Phys. 1996; 79(9):7353 – 7361. [24] Ong ZY, Pop E. A two-temperature model of narrow-body silicon transistors under steady state and transient operation. Proceedings of the 3rd Energy Nanotechnology International Conference, Amer. Soc. Mech. Eng.: New York, NY, 2009. [25] Casas-Vazquez J, Jou D. Nonequilibrium temperature versus localequilibrium temperature. Physical Review E (Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics) 1994; 49(2):1040 – 8. 26

Figure 12: Fundamental mode excited by focused irradiation. The atoms and the mesh are both colored by the phonon temperature.

[26] Baranyai A. Temperature of nonequilibrium steady-state systems. Physical Review E (Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics) 2000; 62(5):5989 – 97. [27] Bertin E, Dauchot O, Droz M. Temperature in nonequilibrium systems with conserved energy. Physical Review Letters 2004; 93(23):230 601 – 4. [28] Braga C, Travis KP. A configurational temperature Nose-Hoover thermostat. Journal of Chemical Physics 2005; 123(13):134 101 – 1. [29] Casas-Vazquez J, Jou D. Temperature in non-equilibrium states: a review of open problems and current proposals. Reports on Progress in Physics NOV 2003; 66(11):1937 – 2023. [30] Ivanov DS, Zhigilei LV. Combined atomistic-continuum model for simulation of laser interaction with metals: Application in the calculation of melting thresholds in Ni targets of varying thickness. Appl. Phys. A 2004; 79(4-6):977 – 981. [31] Xu XF, Cheng CR, Chowdhury IH. Molecular dynamics study of phase change mechanisms during femtosecond laser ablation. ASME J. Heat Transfer 2004; 126(5):727 – 734. [32] Ivanov DS, Zhigilei LV. Combined atomistic-continuum modeling of shortpulse laser melting and disintegration of metal films. Phys. Rev. B 2003; 68(6):064 114 –. [33] Koci L, Bringa E, Ivanov D, Hawreliak J, McNaney J, Higginbotham A, Zhigilei L, Belonoshko A, Remington B, Ahuja R. Simulation of shock-induced melting of Ni using molecular dynamics coupled to a twotemperature model. Phys. Rev. B 2006; 74(1):12 101–1–4. [34] Duff WH, Zhigilei LV. Computational study of cooling rates and recrystallization kinetics in short pulse laser quenching of metal targets. 8th International Conference on Laser Ablation, vol. 59, Hess W, Herman P, Bauerle 27

D, Koinuma H (eds.), IOP Publishing Ltd: Bristol, England, 2007; 413– 417. [35] Lin Z, Zhigilei LV, Celli V. Electron-phonon coupling and electron heat capacity of metals under conditions of strong electron-phonon nonequilibrium. Phys. Rev. B 2008; 77(7):075 133. [36] Cheng CR, Xu XF. Mechanisms of decomposition of metal during femtosecond laser ablation. Phys. Rev. B 2005; 72(16):165 415 –. [37] Lin Z, Johnson RA, Zhigilei LV. Computational study of the generation of crystal defects in a bcc metal target irradiated by short laser pulses. Phys. Rev. B 2008; 77(21):214 108. [38] Wagner GJ, Jones RE, Templeton JA, Parks ML. An atomistic-tocontinuum coupling method for heat transfer in solids. Comp. Meth. Appl. Mech. Engin. 2008; 197(41-42):3351 – 65. [39] Padgett CW, Brenner DW. A continuum-atomistic method for incorporating Joule heating into classical molecular dynamics simulations. Mole. Sim. 2005; 31(11):749 – 757. [40] Irving DL, L D, Padgett CW, Brenner DW. Coupled molecular dynamics/continuum simulations of Joule heating and melting of isolated copper aluminum asperity contacts. Mod. Simul. Mater. Sci. Eng. 2009; 17(1):015 004–+. [41] Allen PB. Theory of thermal relaxation of electrons in metals. Phys. Rev. Lett. 1987; 59(13):1460 – 3. [42] Hertel T, Fasel R, Moos G. Charge-carrier dynamics in single-wall carbon nanotube bundles: a time-domain study. Appl. Physics A 2002; 75(4):449 – 465. [43] John F. Partial Differential Equations. 4th edn., Springer: New York,NY, 1982. [44] Ashcroft NW, Mermin ND. Solid state physics. Brooks/Cole: London, 1976. [45] Tanuma S, Powell CJ, Penn D. Calculations of electron inelastic mean free paths for 31 materials. Surface and Interface Analysis 1988; 11(11):577– 589. [46] Hoover WG. Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A 1985; 31:1695–1697. [47] van Gunsteren W, Berendsen H. Algorithms for brownian dynamics. Molecular Physics 1982; 45(3):637–47.

28

[48] Evans DJ, Hoover WG, Failor BH, Moran B, Ladd AJC. Nonequilibrium molecular dynamics via Gauss’s principle of least constraint. Phys. Rev. A 1983; 28(2):1016 – 21. [49] Duffy D, Rutherford A. Including the effects of electronic stopping and electron-ion interactions in radiation damage simulations. J. Phys.-Cond. Mat. 2007; 19(1):11 pp. [50] Ikeshoji T, Hafskjold B. Non-equilibrium molecular dynamics calculation of heat conduction in liquid and through liquid-gas interface. Mol. Phys. 1994; 81(2):251 – 61. [51] Gear C. Numerical Initial Value Problems in Ordinary Differential Equations. Prentice-Hall: Englewood Cliffs, New Jersey, 1971. [52] Beam R, Warming R. An implicit finite-difference algorithm for hyperbolic systems in conservation-law form. J. Comp. Phys. 1976; 22(1):87–110. [53] Lowrie RB. A comparison of implicit time integration methods for nonlinear relaxation and diffusion. J. Comp. Phys. 2004; 196(2):566–590. [54] Saad Y, Schultz MH. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J. Sci. Stat. Comp. 1986; 7(3):856–869. [55] Foiles S, Baskes M, Daw M. Embedded-atom-method functions for the fcc metals Cu, Ag, Au, Ni, Pd, Pt, and their alloys. Phys. Rev. B 1986; 33(12):7983–91. [56] Eesley GL. Generation of nonequilibrium electron and lattice temperatures in copper by picosecond laser pulses. Phys. Rev. B 1986; 33(4):2144 – 51. [57] Tersoff J. new empirical-approach for the structure and energy of covalent systems. Phys. Rev. B 1988; 37(12):6991–7000. [58] Tersoff J. modeling solid-state chemistry - interatomic potentials for multicomponent systems. Phys. Rev. B 1989; 39(8):5566–5568. [59] Xiao S, Belytschko T. A bridging domain method for coupling continua with molecular dynamics. Computer Methods in Applied Mechanics and Engineering 2004; 193(17-20):1645–69. [60] Liu X, Li S. Nonequilibrium multiscale computational model. Journal of Chemical Physics 2007; 126(12):124 105. [61] Zimmerman JA, Klein PA, III EBW. Coupling and communicating between atomistic and continuum simulation methodologies. Multiscaling in Molecular and Continuum Mechanics: Interaction of Time and Size from Macro to Nano. Springer: New York, NY, 2007; 439–455.

29

[62] Lukes JR, Zhong HL. Thermal conductivity of individual single-wall carbon nanotubes. ASME J. Heat Transfer 2007; 129(6):705 – 16, doi: 10.1115/1.2717242. [63] Schelling PK, Phillpot SR, Keblinski P. Comparison of atomic-level simulation methods for computing thermal conductivity. Phys. Rev. B 2002; 65(14):144 306–12. [64] Daligault J, Mozyrsky D. Ion dynamics and energy relaxation rate in nonequilibrium electron-ion systems. Phys. Rev. E 2007; 75(2, 2):026 402. There are five basic parameters to determine the two heat capacities, the two conductivity functions and the exchange function for the TTM and only the three that concern the electrons are needed if the system is completely filled with atoms. 1. The heat capacity per volume of the phonon system in the classical limit is given by the Dulong-Petit law: ρp c p =

3kB Vα

(A.1)

where 1/Vα is the atom density i.e. the number of atoms per unit cell divided by the volume of the unit cell. Note that this rule can be problematic for non-bulk lattices like carbon nanotubes and for temperatures below the Debye temperature where the measured heat capacity is typically lower than the Dulong-Petit value due to quantum effects. 2. The heat capacity per volume of the electron system, ρe ce = γθe , can be estimated from the free electron density. Assuming the free electrons behave like a kinetic gas : ρe c e =

2 π 2 ne kB θe = γ(ne )θe 2 EF

(A.2)

For metals, the number of electrons ne is approximately the number of valence electrons per atom times the number of atoms. For semiconductors ne is approximately the doping concentration. The Fermi energy EF can be estimated with ne (A.3) EF = EC + kB θe ln NS and NS = 2



2πmkB θe h2

 23

(A.4)

For Cu γ ≈ 96.6 [J/m3K 2 ] [56], for CNT/graphite γ ≈ 1.5 [J/m3 K 2 ] [42], and for Si γ ≈ 0.00059 [J/m3K 2 ]; the variation being mainly a factor of the difference in the density of free electrons ne .

30

3. The thermal conductivity of the phonon system can be estimated via simulation, in particular MD gives the appropriate classical estimates of kp . For Cu kp ≈ 11.6 [W/mK] (at 300K, [61] using a direct method), for C kp ≈ 100.0 [W/mK] (at 500K, [62] using the Green-Kubo method), and for Si kp ≈ 120.0 [W/mK] (at 500K, [63] using a direct method with extrapolation to bulk). Note the measured values of the total conductivity are about 400 [W/mK] for Cu, 300 − 10, 000 [W/mK] for CNTs and about 130 [W/mK] for Si. 4. The thermal conductivity of the electron system can be estimated from the Franz-Wiedemann law ke = Lσθe = k0

θe θp

(A.5)

One caveat regarding this estimate is the electron-electron collisions must be elastic, see [44, Chapter 16]. Here the Lorentz constant is L = 2.443 × 10−8 [W Ω/K 2 ]. The electrical conductivity for Cu is σ ≈ 107 [1/Ωm], for graphite is σ ≈ 104 − 107 [1/Ωm], and for Si σ = µne ≈ 10−1 − 105 [1/Ωm] for Si, where µ is the carrier mobility. For CNTs a value of 105 [1/Ωm] was chosen as representative, and likewise for Si a middle of the range value was chosen 103 [1/Ωm]. 5. The electron-phonon coupling coefficient can be measured using pumpprobe experiments, see, e.g., [42, 56], or estimated from theory, see, e.g., [35, 41, 64]. In general, the Allen theory [41] predicts that the exchange function follows a linear relationship g ∼ (θe − θp ) at low electron temperatures, which is the model commonly used for metals, and g ∼ (θe − θp )5 at high temperatures. For Cu the measured coefficient for this linear model is ge-p = 2.6 × 1017 [W/m3 K] [56]. For semiconductors, G. Chen [19, Chapter 8] provides an empirical formula for the linear coefficient ge-p ≈ 10−12 [W/m3 K] ne which gives a value of ge-p = 1.0×1011 [W/m3 K] for doped Si. Carbon nanotubes apparently follow the high temperature rule [42]. 144ζ(5)kB γ λ (θe − θp )5 = h(θe − θp )5 (A.6) g= 2 π~ θD where ζ(5) ≈ 1.0369, θD is the Debye temperature, λ is a measured decay rate, and ultimately the coefficient h = 3.7 × 104 [W/m3 K 5 ] for CNTs.

31