## Abstract

A thermomechanical equivalent continuum (TMEC) theory is developed for the deformation of atomistic particle systems at arbitrary size scales and under fully dynamic conditions. This theory allows continuum interpretation of molecular dynamics (MD) model results and derivation of thermomechanical continuum constitutive properties from MD results under conditions of general macroscopically transient thermomechanical deformations, which are not analysed by statistical mechanics. When specialized to the more specific conditions of non-deforming systems in macroscopic equilibrium, this theory yields certain results that are identical to, or consistent with, the results of statistical mechanics. Coupled thermomechanical continuum equations and constitutive behaviour are derived using MD concepts in a time-resolved manner. This theory is a further advancement from the purely mechanical equivalent continuum (EC) theory developed recently. Within the meaning of classical mechanics, the TMEC theory establishes the ultimate atomic origin of coupled thermomechanical deformation phenomena at the continuum level. The analysis is based on the decomposition of atomic particle velocity into a structural deformation part and a thermal oscillation part. On one hand, balance of momentum at the structural level yields fields of stress, body force, traction, mass density and deformation as they appear to a macroscopic observer. On the other hand, balance of momentum for the thermal motions relative to the macroscopically measured structure yields the fields of heat flux and temperature. These quantities are cast in a manner as to conform to the continuum phenomenological equation for heat conduction and generation, yielding scale-sensitive characterizations of specific heat, thermal conductivity and thermal relaxation time. The structural deformation and the thermal conduction processes are coupled because the equations for structural deformation and for heat conduction are two different forms of the same balance of momentum equation at the fully time-resolved atomic level. This coupling occurs through an inertial force term in each equation induced by the other process. For the structural deformation equation, the inertial force term induced by thermal oscillations of atoms gives rise to the phenomenological dependence of deformation on temperature. For the heat equation, the inertial force term induced by structural deformation takes the phenomenological form of a heat source.

## 1. Introduction

Material deformation processes fundamentally occur at the atomic level through atomic motion and atomic structural change. Continuum theories of material deformation at the macroscopic, mesoscopic and microscopic scales are primarily phenomenological descriptions of the atomic events. These theories involve quantities such as stress, body force, traction, displacement, temperature, mass density and internal state variables (ISVs). Molecular dynamics (MD) theories, on the other hand, involve quantities such as interatomic force, external force, atomic mass and atomic displacement. While MD theories provide explicit resolution of particle motions and structures of atomic systems, continuum theories provide approximate characterizations of the atomic events and atomic structural states in an aggregate sense. Continuum models often have far fewer numbers of degrees of freedom (DOF) than MD models. At the atomic level, structural deformation and thermal fluctuation make up the total atomic motion. MD models explicitly track the total displacement. Continuum theories, however, treat these two parts of atomic motion separately and, to a degree, independently. Specifically, while the structural deformation part is explicitly modelled at the continuum level, the thermal oscillation part is only accounted for phenomenologically, in the form of heat energy. Reconciliation of the differences in the two descriptions and, more importantly, the integration of the two frameworks of analyses for cross-scale characterization are a great challenge in physics, material science and mechanics.

One of the current trends in computational materials science is the drive to forge links between different length and time-scales, as reviewed by, e.g. Curtin & Miller (2003), Miller (2003), Tan (2003) and Liu *et al*. (2004*b*). In this endeavour, simultaneous atomistic and continuum investigations (e.g. Buehler *et al*. 2003,, 2004) have been increasingly carried out. At the atomistic scale, models based on quantum mechanical, especially density-functional theories for electronic properties are linked to MD and kinetic Monte Carlo simulations. Subsequent steps include integration of discrete and continuum models and coarse graining to higher scales. In particular, the quasicontinuum (QC) approach (Tadmor *et al*. 1996*a*,*b*; Miller *et al*. 1998; Shenoy, *et al*. 1999*a*,*b*; Knap & Ortiz 2001; Shilkrot *et al*. 2002*a*) and the atomic-scale finite-element method (AFEM; Liu *et al*. 2004*a*) integrate atomistic description into a continuum framework of analysis. Their advantages include accounting for certain atomistic deformation mechanisms, reduction of numbers of DOF and facilitation of multiscale hierarchical analyses of material behaviour. Similarly, other approaches that integrate MD and continuum descriptions (e.g. Kohlhoff *et al*. 1991; Gumbsch 1995; Rudd & Broughton 2000; Ogata *et al*. 2001; Shilkrot ,*et al*. 2002*a*,*b*), Wagner & Liu (2003), Park *et al*. (2004) and Shilkrot *et al*. (2004) also offer multiscale modelling capabilities. A primary feature in many of these models is the simultaneous use of atomistic and continuum regions.

A related but separate issue in multiscale analysis is the development of descriptions of material behaviour from pure atomistic simulations without the use of any phenomenological continuum regions. This question can be stated differently. If MD calculations can be carried out on a satisfactorily large scale, how can we derive continuum constitutive laws from the results of such calculations? This is clearly a quite realistic task as far as the behaviour of nanowires, nanotubes and nanocoils are concerned, since many full-scale MD simulations are already being carried out. Statistical mechanics is a very important theoretical framework that has provided highly successful approaches in this endeavour. For example, the Kubo formula (e.g. Kubo *et al*. 1991) establishes the relation between the atomistic behaviour and the time-independent constitutive properties of materials under macroscopically non-deforming or infinitesimally small deformation conditions. However, it does not pertain to the behaviour of materials undergoing large and macroscopically dynamic (fully non-equilibrium at the macroscopic scale) deformations. Indeed, the theoretical framework for obtaining the constitutive behaviour (mechanical and thermal) of an atomistic system undergoing macroscopically large and dynamic deformations has not been developed so far. Even for quasi-equilibrium conditions, one can argue that statistical mechanics does not maintain strict fidelity to molecular processes in a spatially resolved manner. Specifically, statistical mechanics does not provide a characterization of the non-uniform deformation fields (e.g. the spatial distributions of strain and strain rate). Although stress can be obtained using statistical mechanics under equilibrium or quasi-equilibrium conditions, the work conjugacy between the stress field and the rate-of-deformation field for spatially non-uniform deformations is mostly not discussed in statistical mechanics. As such, statistical mechanics may not be as useful for analysing deformation mechanisms associated with dislocations and other defects involving heterogeneous fields at the nanoscale, since it does not provide full field characterizations of relevant quantities. Also, since statistical approaches do not resolve non-local interactions and discreteness- and non-locality-induced length scales, they may not be perfectly suitable for analysing nanoscale processes where such issues are important.

Ultimately, in the context of multiscale modelling and characterization of material behaviour ‘from the ground up’ (from *ab initio*, first principle, MD, to micro-, meso- and macroscopic continuum models), the formulation of continuum descriptions from the results of MD simulations must be carried out with a high degree of faithfulness, under the conditions of general macroscopically dynamic (non-equilibrium), large, and heterogeneous deformations in which the coupling between material behaviour and deformation must be considered. The important first step in such an endeavour is the continuumization of atomic results. Subsequent steps include scaling, reduction of number of DOF and extraction of phenomenological quantities. The development of continuum form fields should maintain strict equivalence of momentum, work rates, mass and kinetic energy between the discrete fields and the continuum fields at all times and all size scales. Furthermore, the non-locality of atomic interactions and discreteness- and non-locality-induced length scales must also be described in the continuum representation. Another issue is the treatment of the thermal oscillation part and the structural deformation part of atomic motion. They must be delineated as one approaches higher length scales, since continuum theories at higher scales treat these two parts separately, while MD models explicitly track both parts as one entity.

To achieve the continuumization stated above, an equivalent continuum (EC) representation has recently been developed for discrete atomistic systems under conditions of general dynamic deformation, Zhou & McDowell (2002) and Zhou (2003). This new theoretical framework uses an explicit expression of all atomic DOF. Work-conjugate continuum stress and deformation fields, as well as all other work- and momentum-preserving kinetic quantities and mass distribution are defined. The EC is equivalent to its corresponding MD system in that, at all times, it preserves the linear and angular momenta of the particle system, it conserves the internal and external mechanical work rates and it has an equal amount of kinetic energy and contains the same amount of mass as the particle system. The equivalence between the discrete system and EC holds for the entire system and for volume elements defined by any subset of particles in the system, therefore, averaging and characterization across different length scales are possible and size-scale effects can be explicitly analysed. This theory is a nanoscale mechanical theory, since just like in the MD system, the kinetic energy of the EC includes contributions from both the thermal part and the structural dynamic part of the atomic velocity. In this sense, the EC is a fully faithful continuum form representation of the MD model. Just as MD models explicitly track the absolute particle motion and the total kinetic energy, the EC model does the same in a fully time-resolved manner. The EC is a continuum form interpretation of the results of atomistic calculations. It is different from the QC or mixed models cited previously. Its primary objective is the continuumization of MD solutions, offering an approach for continuum relations to be formulated using the results of pure MD calculations.

The EC theory does not provide a delineation of the thermal behaviour and the macroscopically observed structural deformation of atomic systems. In this paper, this theory is extended to account for the thermal and mechanical processes at continuum scales. The new thermomechanical equivalent continuum (TMEC) theory described herein allows scale-dependent continuum descriptions of material behaviour in the form of coupled thermomechanical processes of deformation and heat conduction to be obtained from molecular dynamics simulations at the atomic level. The development is based on a decomposition of the atomic velocity into a relatively slower-varying structural deformation part and a high-frequency thermal oscillation part. The establishment of equivalence between the discrete system and the continuum system follows the same approach as the development of the mechanical EC theory. In addition to all mechanical field quantities, one important new aspect is the delineation of thermal energy in the form of atomic vibrations (phonons), introduction of temperature, evaluation of specific heat, calculation of thermal conductivity and the assessment of thermal wave relaxation time. The application in this regard of the TMEC theory developed in this paper is the calculation of temperature and the structural deformation under fully dynamic conditions. Specifically, the theory proposes an answer for the question of ‘which part of the atomic motion is heat and which part is macroscopic structural deformation?’ Most analyses in statistical mechanics on the thermal conductivity assume that the system is non-deforming and in macroscopic equilibrium, therefore, there is no macroscopic deformation velocity. On the other hand, if a solid is indeed going through a macroscopic non-uniform deformation, it would be incorrect to calculate the heat flux using the total atomic velocity, since thermal conduction in solids is relative to the lattice, not relative to fixed spatial locations. This paper offers a theoretical consideration regarding issues such as ‘at the atomistic level, how should heat conduction processes be analysed in materials going through large, non-uniform plastic deformations under fully transient conditions?’

## 2. Thermomechanical equivalent continuum

### (a) Fully resolved atomistic configuration and macroscopically perceived configuration

Consider a dynamically deforming system of *N* particles that occupies space *V* and has an envelope of surface *S*, as shown in figure 1. At time *t*, particle *i* (1≤*i*≤*N*) has position **r**_{i}, displacement **u**_{i} and velocity . denotes material time derivative. The interparticle force applied on particle *i* by particle *j* is *f*_{ij}(**r**_{ij}), where is the central distance between particle *i* and particle *j*. Note that Newton's third law requires that **f**_{ij}=−**f**_{ji}. This analysis here admits pair-wise Lennard–Jones type potentials and multi-body interactions models such as the embedded atom method (EAM) potentials (Daw & Baskes 1983), the modified EAM potentials (Baskes 1987; Baskes *et al*. 1989) and other potentials (Carlsson 1990; Gavezzotti 2002).

To develop a scalable equivalent continuum representation of the particle system, we consider a sub-volume element with closed surface *S*^{e} associated with a subset of *M*(≤*N*) particles in the ensemble. Assume that *M*_{S} out of the *M* particles (*M*_{S}≤*M*) are on surface *S*^{e}, therefore defining it. The remaining *M*−*M*_{S} particles are in the interior of *V*^{e} and are considered internal particles for *V*^{e}.

Obviously, MD calculations fully track the instantaneous, absolute position **r**_{i} and absolute velocity . Within the framework of classical mechanics, these MD concepts are on solid ground. We can call *V* (or *V*^{e}) the *fully resolved instantaneous atomistic configuration* of the system (or subsystem).

Any attempt to measure **r**_{i} and experimentally would require an ‘observer’ or instrumentation that has a sufficiently fast temporal resolution and a sufficiently fine spatial resolution, although such an observer or such instrumentation may not exist. At least in most cases, the high-frequency (thermal) part of the atomic velocity cannot be or is not explicitly resolved. For example, the thermal oscillation part of atomic motion is only phenomenologically accounted for as heat in continuum theories. This is both convenient and necessary in most microscopic and macroscopic theories. In order to admit such an analysis in the EC theory and allow a transition to higher-scale continuum theories from MD, a TMEC theory is developed in this section. To achieve this objective, a decomposition of the atomic displacement **u**_{i} into a structural deformation part and a thermal oscillation part is considered. Based on this decomposition, a thermomechanical version of the equivalent continuum theory is developed. Specifically, the structural–thermal decomposition can be written as(2.1)The fact that is the high-frequency part and is the relatively ‘slowly’ changing part dictates that the decomposition be carried out with regard to the time evolution of **u**_{i} and not be solely accomplished in a spatial manner without regard to time history. This decomposition is dependent on both the size scale and time scale of the analysis pursued. Therefore, there is no unique demarcation between the two parts that can be universally applied under all scales. However, some basic features for the decomposition can be outlined. Specifically, this division can be obtained through a Fourier analysis in time or in both time and space. This may be accomplished, for example, by obtaining a spectral representation of in the form of(2.2)In the above expressions, , *ν* is frequency and *ν*_{cutoff} is a cutoff frequency whose choice depends on the time and size scales of analysis. Since thermal oscillations are typically at frequencies in the range of 0.5–50 THz (Born & Huang 1988), and structural deformation occurs at much lower frequencies, *ν*_{cutoff} can be determined on a case-by-case basis for specific systems, conditions and size and time-scales. To put this in perspective, we note that the high-velocity motion of atoms at the structural level (associated with, for example, dislocation or vacancy motion) by itself does not necessarily correspond to high frequencies and, therefore, does not pose a challenge for the separation of the thermal vibrations and the structural motion of atoms using the method proposed here. The key is the frequency of motion, not velocity. The duration of such atomic motion is more important. Any mechanical loading (with shock wave propagation being near the high-rate extreme) would most probably have rise times of least 400–500 ps (Gahagan *et al*. 2000), with most shock waves having rise times of the order of several to hundreds of ns. These give rise to structural deformation frequencies of, at most, 2 GHz, which is several orders of magnitude lower than the frequencies of thermal oscillations of the order of 0.5–50 THz. Challenges might arise for laser-induced processes with durations or rise times of the order of several ps (Branch & Smith 1996) or shorter than 1 ps (corresponding to a frequency of 1 THz). Indeed, ultra-high frequency deformation is intimately related to the thermal motion of atoms. The handling of this issue is a topic worth considering by many. The discussions in this paper are limited to conditions under which the frequencies of structural deformation and those for the thermal motions are clearly distinguishable. This assumption permits the analysis here to be applied to most problems, including most laser-induced shock processes at the ns level.

Inverse transforms would yield the structural and the thermal velocities, i.e.(2.3)Even though the analysis is carried out for each atom separately in a purely temporal manner, it can allow the spatial propagation of lattice waves at different frequencies to be revealed. This is simply because the spatial distribution of the structural part of atomic motions is represented in **u**_{i} for each atom. It is important to point out that the spectral analysis outlined here is, perhaps, only one out of many possible approaches for obtaining the decomposition in equation (2.1). Appendix A provides a discussion on the averaging methods for decomposition frequently used in the literature. The analysis that follows assumes that a decomposition has indeed been achieved in an agreeable manner, which may or may not be the approach proposed in equations (2.2) and (2.3). Regardless of how the decomposition is achieved or how one prefers to achieve the decomposition, the analysis that follows remains the same.

We note that thermal oscillations are usually associated with no net momentum at the system level. Specifically, for systems in macroscopic equilibrium(2.4)where 〈 〉 denotes averaging over a sufficient period of time.

Under the decomposition in equation (2.1), the position of particle *i* can be written as(2.5)This decomposition of atomic positions defines the difference between the *actual configuration* (*V*^{e} and *S*^{e}) and the *macroscopically perceived* configuration ( and ) of an arbitrary volume element. An illustration of the macroscopically perceived configuration and the actual configuration is given in figure 2. While the pure mechanical EC theory in Zhou & McDowell (2002) and Zhou (2003) is formulated over the actual configuration, the TMEC here must be formulated over the macroscopically perceived configuration. The reason is obvious: macroscopic continuum theories are written relative to the macroscopically observed size and shape of a system ( and ).

The force on particle *i* due to atoms or agents that are *external* to the system under consideration is . The total force on *i* is(2.6)Here, the summation is over those particles inside the system of *N* particles that interact directly with particle *i*. It is worthwhile to point out that, due to non-local interactions, external force can exist for particles both in the interior of *V* and on surface *S*. Note that the concepts of internal and external forces are specific to the particular sub-volume of considered. So, in general, and , except for . For example, for the volume element *V*^{e} illustrated in figure 1, the forces between atoms 1 and 2 (**f**_{12}, **f**_{21}) as well as those between atoms 3 and 4 (**f**_{34}, **f**_{43}) are internal; while the forces between atoms 5 and 6 (**f**_{56}, **f**_{65}) and those between atoms 7 and 8 (**f**_{78}, **f**_{87}) are external. Note that all these forces are internal when the system is considered as a whole. This point will become clearer shortly.

Note that balance of momentum must be satisfied by a system as a whole and by any portion of the system. For a given portion of a system, internal forces and external forces are fundamentally different and play different roles. Therefore, if two systems are to be equivalent and balance of momentum is to be satisfied at any size scale, their internal force work rate, external force work rate and inertial work rate at all size scales must be equal.

To delineate the internal forces and external forces for volume element systematically, we use a local (capital) index ‘*I*’ for the atoms associated with *V*^{e}. With this notation, each atom inside has two indices, one is the local index *I* (1≤*I*≤*M*) and the other is its global index *i* (1≤*i*≤*N*) as the atom is also part of the complete system . Since a unique correspondence between *I* and its global counterpart *i* can be established, we use them interchangeably here for convenience of discussion. Under this notation, ‘*j*≠*I*’ and ‘*j*=*I*’ should be interpreted, respectively, as ‘*j*≠*i*’ and ‘*j*=*i*’. For example, *j*=*I* should be read as ‘the particle with global index *j* and the particle with local index *I* (and therefore global index *i*) are the same particle’.

These notations are used here to delineate the total force on atom *I* in exerted by other atoms *also inside* (either in the interior of or on surface ) and the total force exerted by atoms and agents *outside* . A distinction must be made between these internal and external interactions. Note that the total force on *I* is , and(2.7)Here, is the fraction of the atomic bond that is spatially within element . It pertains to the bond between atoms *I* and *J* that are both inside . In general, when atoms are randomly distributed (as in amorphous materials), is determined by the dihedral angle of the element as a fraction of the sum of such angles (≤360°) of all elements associated with the particular bond. Specifically, , with being the dihedral angle in element ‘*e*’ associated with the bond between atoms *I* and *J*, and *q* being the number of elements connected to the bond. Consider, for example, the BCC lattice in figure 3. The bond between atoms 1 and 3 is shared by four tetrahedral cells (each of the four tetrahedral cells is considered a volume element with *M*=4). Therefore, for the tetrahedral element defined by atoms 1, 2, 3 and 4 (and for each of the other three cells) for this bond. For bonds on surface (both atoms of the bond are on ), the sum of such angles is <360°.

### (b) Two perspectives of atomic particle motion

The equation of motion for atom *i* is(2.8)There are two distinctive yet related perspectives towards this equation. These two aspects can be illustrated using two different observers. An illustration of this issue is given in figure 4. The first observer is stationary and has a macroscopic resolution for atomic motion. This observer uses a slower-responding oscilloscope that is insensitive to motions at frequencies above *ν*_{cutoff}. Consequently, is not explicitly resolved or is ‘invisible’ at the size and time-scales analysed by this macroscopic observer. To understand the motion as it appears to the macroscopic observer, we rewrite equation (2.8) as(2.9)where is the high-frequency inertial force . Since is not directly measurable for the macroscopic observer, is not directly measurable. The macroscopically observed trajectory is (see figure 4)(2.10)However, the observer can infer the effect of and its existence from the deviation of the average path (which he or she actually measures) from what is predicted by(2.11)In the above equation, and are, respectively, the internal and external forces on atom *i* when it is on the macroscopically perceived trajectory in equation (2.10). represents the solution to the above equation, which in general differs from the displacement solutions to equations (2.8) and (2.9). Equation (2.11) consists of only kinetic and kinematic quantities perceivable to the observer. It describes the trajectory if thermal oscillations simply did not occur ( at all times). It is important to point out that equation (2.11) is cited here for comparison and perspective only. It is not used to develop any equation of consequence in this paper, since it is not an actual equation for the physics here. Rather, it is a construed equation one can use for comparison with equation (2.10) in order to determine if or thermal oscillations actually exist. We must note that, although equation (2.10) describes the ‘perceived’ deformation of an atomic system at longer time and higher size scales, it is equation (2.9) that governs the actual deformation on which the ‘perception’ is based. The perceived state is important since it defines the volume and surface that form the configuration required by a macroscopic analysis. The analysis of macroscopic deformation (balance of momentum and balance of energy), however, must use equation (2.9).

The second perspective towards equation (2.8) is based on the observation of an observer travelling with the structural velocity of an atom. Assume that this observer uses an oscilloscope that can resolve motions at all frequencies. To this observer, the particle motion then appears to be governed by(2.12)where is the inertial force due to the structural deformation. Compared with , is slow-changing. Since is invisible to this moving observer, he or she can only infer the effect (and, therefore, the existence) of by noticing the deviation of from what is predicted by(2.13)In the above equation, and are, respectively, the internal and external forces on atom *i* calculated using the atomic position as perceived by the observer moving with . represents the solution to the above equation, which in general, differs from the displacement solutions to equations (2.8), (2.9), (2.11) and (2.12). Equation (2.13) consists of only kinetic and kinematic quantities perceivable to the moving observer. It describes the trajectory of the thermal oscillations if the structural deformation simply did not occur ( at all times). It is important to point out that, just like equation (2.11), equation (2.13) is cited here for comparison and perspective only. It is not used to develop any equation of consequence in this paper. It will become clear later that equations (2.9) and (2.12) govern, respectively, the macroscopically observed deformation and the thermal behaviour of a material system. These two equations yield the coupled mechanical and thermal equations at various scales, respectively. We will also see that the inertial forces and give rise to the coupling between the mechanical and the thermal processes. The fact that equation (2.12) is written from the perspective of an observer moving with a particle's macroscopic velocity indicates that heat conduction analysed here is relative to material mass at the macroscopically perceived positions . Therefore, a material time derivative should be used later in the continuum representation of internal energy. This also requires the heat flux to be defined relative to macroscopically perceived, current, atomic positions (as opposed to **r**_{i} or the original positions **r**_{i0}).

### (c) Structural deformation: balance laws and field quantities

To analyse the behaviour of the system associated with macroscopically observable structural deformation , we consider a variation of this velocity and invoke the principle of virtual work for an arbitrary element *as it appears to the macroscopic observer*. Under this condition, the equivalence can be stated as(2.14)where and are the volume and surface of the element defined by the structural part of the atomic position vectors , as illustrated in figure 4. is defined on and . is the fraction of atom *I* that is attributed to element . For example, atom 1 in figure 3 is shared by 24 tetrahedral elements, therefore, for each element. For periodic and amorphous structures alike, can be defined through , with being the solid angle (three-dimensional) or angle (two-dimensional) subtended by an element and *q* being the number of elements connected to an atom. Note that the delineation of and is relative to with defined through equation (2.7). External forces on atoms *in the interior* of contribute only to the body force density , therefore, factor *κ _{I}* is always taken as unity (

*κ*=1) for atoms in the interior of . The external forces exerted on atoms on surface are considered to contribute solely to surface traction , therefore,

_{I}*κ*=0 for atoms on . This partition is somewhat arbitrary and, indeed, any choice 0≤

_{I}*κ*≤1 for surface atoms (along with

_{I}*κ*=1 for interior atoms) will allow external work rate to be preserved. However, the choice of

_{I}*κ*=0 for surface atoms (along with

_{I}*κ*=1 for interior atoms) has the clear advantage of yielding zero body force density as non-local external forces approach zero. This outcome is consistent with local continuum theories.

_{I}The high-frequency term in the above equation is(2.15)The association of with the external force in equation (2.14), rather than , may appear to be somewhat arbitrary at the first glance. However, it is strictly required by Newton's third law and by the necessity of maintaining balance of momentum at all size scales. It is also required by the stipulation that full work rate equivalence be maintained between the EC and the molecular system. Furthermore, it is the only treatment that allows consistent handling of internal, external and inertial work rates between the explicitly resolved particle motion as perceived by the observer and an equivalent continuum representation based on the macroscopically perceived particle motion. In particular, we note that Newton's third law requires that the sum of forces internal to a system, or a portion of a system, vanish. Otherwise, the system or the portion of that system would accelerate by its own effects and create energy and momentum. Here, the high-frequency inertial force term cannot be associated with the internal force because, in general,(2.16)In contrast, since are internal forces, Newton's third law implies that(2.17)

The component equations resulting from equation (2.14) are(2.18)where is defined relative to . The break-up of equation (2.14) into the above component equations or the requirement of strict term-to-term satisfaction of equation (2.14) fully reflects the dynamics of the system. It ensures the satisfaction of balance of momentum and full dynamic equivalence between the TMEC and the particle system at any size scale. The associations of internal forces to internal stress only and external forces to body force and surface traction only are strictly required by balance of momentum and conservation of energy. Otherwise, if one allows internal forces to generate external traction and/or body forces, a system (or part of a system) would accelerate and gain momentum and energy by the effect of its internal forces alone. This would violate balance of momentum and balance of energy. Similarly, if a part of the external forces is used to generate internal stress, part of the momentum and work from the external forces would vanish, again causing violation of balance of momentum and balance of energy.

The relations in equation (2.18) provide the mechanisms for obtaining the deformation part of the equivalent continuum equations. The following analysis is based on these relations. It is important to note that the thermal velocity term in equations (2.14) and (2.18) accounts for the effect of the inertial force associated with the thermal fluctuations of the atoms on structural motion. It acts like an invisible external force. Clearly, it is not associated with a neighbouring (or ‘interacting’) atom, since the neighbouring atom can have totally different thermal velocity and thermal acceleration. The inertial forces of two atoms are not equal in magnitude and opposite in direction. As a matter of fact, no part of the thermal velocities of any two neighbouring atoms can be regarded as being coupled or out of mutual interactions. The inertial force term cannot be treated together with the internal (interactive) forces or stress. Ultimately, this is one of the reasons why thermal motions of atoms do not give rise to a contribution in terms of atomic velocity to stress directly.

The continuum description of the structural deformation is obtained via(2.19)where represents continuum positions inside just like * x* represents positions inside

*V*

^{(e)}. is the shape function and should be interpreted as . The superscript ‘(

*e*)’ is omitted for brevity. Details regarding shape functions can be found in most finite-element method texts. Here, note that is defined using the structural part of the atomic positions over . The corresponding gradient of the virtual velocity field is(2.20)Here, are gradients of the shape functions and ‘⊗’ denotes the tensor product of two vectors. Since(with ( )

^{T}denoting the transpose of a tensor), the virtual work rate of the stress tensor with respect to the virtual rate of deformation is given byBecause are completely arbitrary and independent DOF, the first relation of equation (2.18) leads to the equations for the stress over in the form of(2.21)Note that one of the basic requirements for shape functions is ; therefore, . In general, for an element with

*M*atoms, equation (2.21) yields 3

*M*−6 independent equations. Since the number of independent components in is six, the problem of finding a

*constant*,

*spatially non-varying*,

*average*is overspecified for any choice of

*M*greater than four. Consequently, it is impossible to find a

*constant*,

*spatially non-varying*,

*average*work-equivalent stress over the volume associated with an arbitrary subset of the particle ensemble. It important to point out that this is

*not*to say that a work-conjugate stress field cannot be found on an arbitrary size scale in general. Of course, properly formulated

*spatially varying*can always be found for an element at any scale. More discussions on this follow in the next paragraph. Because of this reason, in general, an average (constant) stress cannot be associated with a work-conjugate deformation field. Although it may be desirable to do so, such a task is not possible because of the disparate number of DOF for the discrete particle subset and the fixed dimensional order of the stress tensor. Parity in the DOF and the order of the stress tensor occurs (six equations and six unknown stress components) only for the simplest three-dimensional cell, the tetrahedron that is associated with four particles. In two dimensions, triangles associated with three atoms are the only possible choice (three equations and three independent unknown stress components). This particular level of continuum characterization is very useful because it fully recognizes the effects of heterogeneities and steep gradients at the scale of individual atoms. This process of establishing dynamical equivalence between the continuum and MD formulations, although performed at interatomic scales, is also very important because it yields the lowest scale continuum fields that may be subsequently subjected to various treatments of continuum averaging, including those of statistical mechanics.

It is important to point out that for *M*>4 in three dimensions and *M*>3 in two dimensions, equation (2.21) requires, in general, spatially varying stress fields with proper choice of integration points in for the relevant numbers of unknowns. Although somewhat more computationally involved, such a pursuit would be quite useful and important since it permits scaling and allows size effects to be quantified through the variation of the size of .

To obtain the traction over the surface area of , consider a surface element defined by *L* particles. The continuum virtual velocity over is(2.22)where is the appropriate shape function over . Substitution of the above equation into the second equation of equation (2.18) yieldswhere is the fraction of that can be attributed to , since may be only a portion of and particle *I* may be on the boundary of (shared by the rest of ). can be defined through , with *q* being the number of surface areas connected to atom *I*. Again, the arbitrariness of requires that , with *I*=1,2,…,*L*. The solution to the above yields the corresponding traction field as(2.23)where are vector solutions of the linear system of equations in the form(2.24)In the above equations, , with *I*, *J*=1,2,…,*L*.

The equations for the body force density resulting from the third equation in equation (2.18) are(2.25)where are the vector solutions of(2.26)In the above relations, , with *I*, *J*=1,2,…,*M*. Similarly, the thermal oscillation-induced ‘body force’ density are(2.27)where are the vector solutions of(2.28)

The mass density for is(2.29)where (*K*=1,2,…,*M*) are solutions of(2.30)with and . Equation (2.30) and the requirement that be independent of yields(2.31)The replacement of in equations (2.14) and (2.18) with the actual structural deformation part of the velocity yields(2.32)

Direct term-to-term equivalence is maintained here, establishing conservation of energy and equivalence of work rates from the perspective of the macroscopic observer. The term associated with thermal oscillations of atoms is(2.33)This term represents the coupling or rate of energy exchange between the structural mode of the deformation and the thermal mode of the deformation as it appears to the macroscopic observer. A positive value of this term indicates energy input from the thermal mode into the mechanical mode. Such cases arise, for example, in thermally driven deformations. Although the thermal force is invisible to the macroscopic observer, its effect on structural deformation, and perhaps more importantly, the energy exchange the above term represents, can be clearly measured and plays an important role in microscopic and macroscopic phenomena. Examples include thermoelastic dissipation, thermoplastic dissipation, electromagnetically driven heat generation and deformation, thermal shock and deformation driven by temperature changes.

It is important to point out that the ‘macroscopic’ representation of equation (2.32) differs from the regular continuum balance of energy relation in the form of(2.34)Although this relation is quite similar in form to that in equation (2.32), they have somewhat different meanings depending on the scale. The quantities in equation (2.34) (denoted by an underscore) are *phenomenological* quantities developed for approximate descriptions of the atomic counterparts in equation (2.32). They may not be based on strict, full equivalence of momentum, energy, work rates and mass at the atomic level. Since they are higher-scale concepts with much fewer numbers of DOF, they do not provide resolution of atomic scale behaviour as do the quantities in equation (2.32) (denoted by an overbar). The following correspondence can be outlined based on physical mechanisms, without necessarily implying equality (equality holds if all the atomic DOF are used in the continuum representation):(2.35)

Of particular interest is the identification of the high-frequency, thermal coupling term with the stress work rate in the first correspondence (rather than with the body force term in the second correspondence) in (2.35). This identification reflects the continuum phenomenological description of thermoelastic dissipation and thermoplastic dissipation through the stress work, instead of through a ‘body-force-like’ term like that in equation (2.32). Although this continuum phenomenological treatment is convenient and allows material behaviour to be described with various degrees of accuracy at different length and time-scales, the theoretical analysis here shows that it is fundamentally incompatible with the requirement of full work rate and momentum equivalence between the continuum description and the atomic reality. The form of such a treatment does not allow full momentum and work equivalence to be obtained, even if the model is assigned the same number of DOF as the discrete atomic system. From the viewpoint of multiscale material behaviour characterization from the ‘ground up’ through molecular dynamics simulations, this association of thermal dissipation and thermo-structural coupling with stress work rate can be used to formulate proper higher-scale, lower-DOF, phenomenological constitutive laws. For example, it can be defined that(2.36)

Therefore,(2.37)It must be pointed out that equation (2.36) does not give a unique definition for . More importantly, it is not possible to choose such that can be made to satisfy balance of momentum (equation (2.14) or (2.21)). Furthermore, this phenomenological association of thermal dissipation with the stress work rate cannot be construed as to imply dependence of stress on thermal velocity.

### (d) Thermal oscillations: balance of energy and thermal fields

The analysis above is carried out from the perspective of a higher-scale observer with focus on the structural part of the atomic deformation . We now focus our analysis on the thermal part of the deformation , from the perspective of an observer moving at velocity with full resolution for thermal motions.

Equation (2.12) and the principle of virtual work over with respect to yield (after replacement of by in the end)(2.38)In the above relation,(2.39)Also, term-to-term equality is maintained and(2.40)where is defined over through(2.41)Also, is calculated through equations (2.27) and (2.28) with , and being replaced by , and , respectively. represents work coupling between the thermal mode of deformation and the structural mode of deformation as it appears to an observer moving with the structural velocity . Specifically, it is the rate of energy input into the thermal mode of atomic motion.

Equation (2.38) represents balance of energy as it appears to the observer travelling with the structural velocity of an atom in the discrete description or to an observer travelling with velocity of the (continuumized) mass point occupying * x* at time

*t*in the equivalent continuum description. This equation is the atomic description of the thermal process in a coupled thermomechanical deformation event. The continuum phenomenological characterization is(2.42)where is the rate of phenomenological density of heat source (thermal energy generated in the continuum per unit volume and per unit time), is the phenomenological heat flux (thermal energy flowing across unit surface area per unit time),

*e*

_{T}is the density of the kinetic thermal energy or the part of the internal energy associated with the thermal velocity of atoms, and

*e*

_{S}is the density of the potential part of the internal energy associated with interatomic force interactions. Specifically,

*e*

_{S}is due to the thermal displacement away from the macroscopically observed position ; therefore, if all atoms happen to be at their macroscopically observed positions ,

*e*

_{S}=0. Since equation (2.42) is phenomenological in nature and involves much fewer DOF than equation (2.38), full and direct equivalence between the atomic equation and the macroscopic phenomenological equation may not be established. The following correspondences can be outlined based on physical mechanisms, without necessarily implying equality (equality holds if all the atomic DOF are used in the continuum representation):(2.43)

Note that Reynolds transport theorem (Malvern 1969) specifies that

Equations (2.38) and (2.43) allow phenomenological forms of heat equations at various scales to be formed based on the results of MD calculations. Specifically, the second relation above can be used to obtain characterizations of the thermal behaviour of materials. To this end, we define the temperature at an atomic position through(2.44)where *k* is the Boltzmann constant. We note that the standard definition of the *average* temperature for several atoms is only given in the statistical sense over a volume of a system with *M* atoms throughThe extension of this definition to individual atoms here is consistent with the statistical interpretations of temperature for aggregates of atoms. The statistical consistency is in the sense that the temperature of an aggregate of atoms is the average of the temperatures of the individual atoms. Some may argue that thermal energy (therefore temperature) can only be measured (primarily) at higher scales in a statistically averaged, phenomenological manner and it is difficult to measure the temperature of individual atoms and, therefore, the concept of the temperature of an individual atom is difficult to justify. A new and different perspective can be offered here. The advancement in scanning tunnelling microscopy (STM) and atomic force microscopy (AFM) has brought to reality the ability to resolve and manipulate the positions of individual atoms, e.g. Alivisatos *et al*. (1996) and Fishlock *et al*. (2000). The thermal motions of individual atoms play a significant role in such processes and can be analysed. If the thermal velocity of an atom can be clearly defined, the kinetic energy (thermal energy) associated with it can be defined and quantified as well. The atomic temperature defined in equation (2.44) has relevance, application and scientific justification for nanoscience of individual atoms and small clusters of atoms. This extension of definition of temperature to an individual atom allows thermal analyses to be carried out at the nanoscale. We also note that this definition fully lends itself to interpretations in a statistical sense for higher-scale systems. Under this condition, the temperature field inside can be expressed as(2.45)It is specifically noted that the concept of temperature, just like the concept of thermal velocity, is inherently scale-dependent. The equivalence of thermal energy at an arbitrary size scale of allows specific heat *c*^{(e)} to be defined through(2.46)where is the kinematic part of the internal (thermal) energy per unit volume inside and is the potential part of the total thermal energy per unit volume inside , such that the last line of equations in equation (2.43) are followed and that *e*_{S}=0 if atoms associated with are at their corresponding macroscopically observed positions . To illustrate the consistency of equation (2.46) with the definition of specific heat in statistical mechanics, we consider a system in macroscopic equilibrium with a uniform thermomechanical state ( is non-changing) for which for all values of *I* in . Here, ‘〈 〉’ denotes proper averaging over time. Under such conditions, equations (2.44) and (2.45) yield . Since thermal vibrations of atoms around their macroscopically observed positions have very small amplitudes at low to moderate temperatures, the thermal motions are effectively linear oscillators. Consequently, there is an equal partition of the total internal energy into the kinetic part and the structural part in the statistical sense, i.e. . Therefore, equation (2.46) yields *c*^{(e)}=3*k*/*m*, which is the specific heat per unit amount of mass, giving rise to a total heat capacity of *C _{V}*=(3

*k*/

*m*)(

*m*)(

*N*)=3

*kN*for the system (

*N*is the total number of particles in ). This is the Dulong–Petit law in statistical mechanics (Born & Huang 1988; Toda

*et al*. 1991). Obviously, when specialized to systems in macroscopic equilibrium, the definition of specific heat in equation (2.46) is consistent with statistical mechanics' definition of where

*E*is the total energy of the system.

The governing equation for the heat flux inside is the local version of equation (2.42) in the form of(2.47)This equation and the boundary requirement that(2.48)over allow **q**^{(e)} to be uniquely determined. The determination of **q**^{(e)} at all times yields .

The thermal fields of temperature, specific heat and heat flux carry information about the intrinsic thermal constitutive behaviour of the atomic system. They can be used to formulate phenomenological continuum constitutive descriptions of the thermal behaviour of material systems, using MD calculations. Proper forms of thermal constitutive relations should be chosen within the context of specific systems, time scale, size scale and the conditions involved. To illustrate how this can be achieved, we recognize the wave nature of heat conduction at the atomic level and consider the commonly used phenomenological relation between heat flux and temperature gradient in the form of(2.49)where * τ* is the tensor of relaxation times and

**K**^{(e)}is thermal conductivity.

*is positive definite and, for simplicity, can be taken as isotropic, i.e.*

**τ***=*

**τ***τ*, with

**I***being the second order identity tensor.*

**I***τ*is positive, reflecting the fact that the speeds of thermal waves (Joseph & Preziosi 1989) through atomic structures are finite.

**K**^{(e)}is also positive definite, reflecting the fact that heat flows from high temperature regions to low temperature regions. Equation (2.49) leads to heat equations that are hyperbolic in nature, accounting for the fact that ultimately at the atomic level, the conduction of heat is the propagation of high frequency mechanical (thermal oscillation) waves through the atomic structure. This is especially the case at short time-scales in applications such as heating induced by fast laser pulses, e.g. Hector

*et al*. (1992) and Kosinski & Frischmuth (2001). At higher size and time-scales, the wave nature of the thermal process is smeared out and diffusion characteristics dominate, since sufficient time is available for relaxation to occur fully. Alternatively, if the time-scale of analysis is orders of magnitude longer than

*τ*, a transition to the limiting case with

*τ*→0 is justified, resulting in a parabolic heat equation commonly used at higher scales.

The determinations of the relaxation time tensor * τ* and thermal conductivity

**K**^{(e)}may require two separate calculations. First, the determination of

**K**^{(e)}should, preferably, be carried out under steady-state thermal conditions for which , obviating the need for knowledge of

*. After the conductivity is determined, equation (2.49) can then be used to determine*

**τ***, based on the results of at least one separate MD calculation involving transient thermal responses. This should involve the choice of the nine components of*

**τ***through, for example, curve-fitting, for the best description of the MD data. Materials with symmetries have fewer numbers of independent components for*

**τ***and*

**τ***.*

**K**The specific heat and the thermal conductivity defined here are not thermodynamic quantities. They are fully dynamic, instantaneous measures for the thermal behaviour of atomic systems at various size scales. These definitions do not reflect statistical averaging, which is an intrinsic thermodynamic concept. Rather, they are quantities describing the thermal behaviour of atomic systems with explicit resolutions in time and in structure, yet reckoned in the forms of standard (statistical) specific heat and thermal conductivity as to conform to the continuum heat equation. Naturally, these dynamic quantities show oscillations and dependence on size, atomic structure and composition of lattice waves at the atomic level. It is reasonable to expect that as the size scale of analysis is increased and as proper spatial and temporal averages are taken, these quantities will approach the microscopically or macroscopically measured values. Such analyses constitute future development. The quantification and approach introduced here provides an alternative to statistical mechanics for reaching microscopic, mesoscopic and macroscopic continuum scales from a discrete molecular dynamics framework.

## 3. Example: thermomechanical tensile deformation of a rectangular lattice

Although the solutions of dynamic deformation of atomic systems in general are complex and require numerical treatment, a simple example can be analysed here to illustrate the perspectives that the novel TMEC has brought about. Some of the insights obtained here are revealing in terms of the form of stress, scale dependence and thermal–mechanical coupling.

Consider the uniform tension of a rectangular lattice under external forces **f**^{ext} in figure 5. For simplicity, we assume the deformation is strictly uniaxial. Thus, all atomic displacements and velocities are in the horizontal direction. Furthermore, we assume that the deformation is uniform in the vertical direction, therefore, * u_{i}*(

*,*

**x***t*), and are only functions of time and of and are independent of and . Although such a situation is idealized and cannot be easily brought about in laboratory, it is perfectly allowed and has all the fundamental attributes of a dynamically deforming atomistic particle system, including the satisfaction of equation (2.4). It should also be pointed out at the outset of this analysis that although a time-resolved analysis is carried out at the scale of a unit lattice cell here, interpretation of the thermal behaviour discussed has practical meaning primarily in the sense that the results be viewed as time averages over time. Full consistency of the thermal behaviour with higher-scale characterizations occurs only when the size scale of analysis is much larger than what is the case here, and when a proper time average is taken. However, this example allows the issue of continuum representation of the MD system and the solution procedure to be illustrated.

As shown in figure 5, the deformed lattice has dimensions *a* and *b* in the horizontal and vertical directions, respectively. For dimensional consistency, we can regard the lattice as having unit thickness perpendicular to the plane shown, therefore, the atoms are actually cylinders having unit length perpendicular to the plane of analysis. The material is homogeneous and the mass of each atom is *m*. For simplicity, we assume the cutoff radius *R*_{c} is such that *a*<*R*_{c}<2*a* and *b*<*R*_{c}<2*b*; therefore, non-local interactions do not occur. Under the decomposition of equation (2.1), the atomic displacement has a uniform structural deformation part and a thermal oscillation part. Consequently,(3.1)Since the deformation is assumed to be strictly one-dimensional, and . The calculation for the TMEC fields here uses two-dimensional shape functions for rectangular elements, e.g. Reddy (1993). For the rectangular element defined by atoms 1, 2, 3 and 4 in figure 5, the shapes functions are(3.2)

The interatomic forces vary with time (i.e. *f*_{12}=*f*_{12}(*a*,*b*), *f*_{13}=*f*_{13}(*a*,*b*) and *f*_{32}=*f*_{32}(*a*,*b*)) while they are uniform in the *y*-direction since the deformation is uniform in that direction. All relevant continuum fields are summarized in table 1. Both the results for the macroscopic TMEC representation and for the mechanical EC representation with full atomic displacement resolution are given. The stress, traction, body force and deformation fields are consistent with macroscopic expectations. It is important to note the difference and relationship between the stress representations at the macroscopic level and at the fully resolved atomic level. Both representations have similar forms. While the stress (traction and body force as well) at the fully resolved atomic level is evaluated relative to the instantaneous, true positions * r_{I}* of the atoms which define

*V*

^{(e)}, the stress (traction and body force as well) at the macroscopic level is evaluated relative to the macroscopically perceived (structural) positions of the atoms which define . Both evaluations use the full instantaneous forces

*on atoms. Specifically,(3.3)For the volume element defined by atoms 1, 2, 3 and 4 in figure 5, for*

**f**_{I}

**f**_{12}and

**f**_{32}since these forces are shared with neighbouring cells. However, for

**f**_{13}since this force is fully within the element. These coincide with the Cauchy stress expressions

**σ**^{(e)}and in table 1.

The heat flux component in the *x*-direction for the element analysed in table 1 is(3.4)where *C* is an arbitrary constant resulting from the integration and must be determined by equation (2.48) which takes the form of(3.5)Clearly, in equation (3.4) measures the flow rate of energy in the form of high-frequency mechanical waves through the element . This flow is apparently driven by the rate at which the thermal motion of the atoms on the boundary can extract high-frequency mechanical work from or impart high-frequency mechanical work to the surroundings. Of course, the heat flow so driven by the external excitation on the boundary are effected through an imbalance inside between the rate of change of the kinetic and potential energies associated with the thermal motions of atoms and the rate of work of the body force and inertial force over the thermal oscillations. To put it differently, heat flow across is a manifestation, in the form of energy exchange, of the interaction of the high frequency motion of itself with its surroundings.

At the macroscopic scale, a homogenous single crystal is expected to have a spatially uniform thermal conductivity * K*. The mathematical forms of

*T*

^{(e)},

*c*

^{(e)}and in table 1 for do not necessarily conform to such macroscopic continuum phenomenological expectations, since they show spatial variations at the lattice level. More important, they are fully dynamic solutions which account for the effects of the transient tensile deformation at the macroscopic level. As shown below, when these solutions are specialized to the conditions of a non-deforming plate in macroscopic equilibrium, consistency with the results of statistical mechanics is obtained. The perspective obtained here is threefold.

First, explicit, time-resolved balance of momentum and conservation of energy is maintained in the TMEC analysis at all times. The result on the heat flux is quite illustrative. The significance of the results here should be more understood in the form of proper averages in time. For example, the temperature, specific heat, thermal conductivity and relaxation time in table 1 are more meaningful if interpreted as(3.6)As in experiments, the determinations of *c*^{(e)} and should be carried out under a state of constant and uniform temperature and a steady state with a temperature gradient, respectively. Specifically, at the uniform and constant temperature state, and 〈*c*^{(e)}〉=3*k*/*m*. At the state of a steady and uniform temperature gradient in the *x*-direction,and , yielding a thermal conductivity in the *x*-direction of . The determination of *τ*, on the other hand, requires a transient impulse at one end. These expressions are consistent with what is given by statistical mechanics.

Second, the atomic structure is explicitly accounted for. Under such conditions, the governing equations for the thermal process are then cast in the form of the macroscopic heat equation. The results at the nanoscopic scale naturally reflect the nature of the atomic structure. Consistency with higher length-scale continuum expectations should emerge as one approaches higher scales. This entails the use of progressively large volume elements. In the process, the lattice-level variations of these quantities should average out and the long range order should become apparent and dominant.

Third, thermal oscillations and structural deformations of atomic structures consist of waves of different frequencies (wavelengths). The wave spectrum is another source giving rise to time- and scale-dependence of responses. Macroscopic response can only be observed when the analysis is carried out at spatial and temporal scales sufficient to represent all components of lattice waves. The unit cell analysis above is at a scale well below many significant component wavelengths in the lattice. While this particular analysis at this particular scale shown in this example does not necessarily represent the higher scale behaviour, it indeed characterizes, rather faithfully and accurately, the atomic interactions at the lattice scale, albeit in a continuum form used at higher scales. The use of such a continuum form at all scales allows comparison, analysis of scaling effects and transition from one scale to another.

## 4. Thermoelastic and thermoplastic dissipation

The issue of obtaining a continuum characterization of plasticity due to atomic rearrangement or configuration change should not be obscured by the conservative (nonlinear elastic) nature of interatomic potentials used in MD. The reversibility of atomic reordering and restorability of broken atomic bonds must not be perceived as to prevent MD models to account for plasticity and dissipation in the continuum sense. MD and continuum theories have vastly different resolutions for defects or atomic structure changes. While MD models fully resolve atomic processes, continuum formulations typically employ far too few state variables to be able to capture the dynamic state of the atomic system. Note that defect generation through atomic rearrangement (e.g. dislocations) is the source of plasticity or dissipation at the continuum level, even though such defects can be reversed at the atomic level *if the right conditions are made available*. Such conditions to return the defect structure to any particular state are very hard to come by, due to the large number of DOF at the atomic level. They practically do not occur under deformations which higher-scale continuum theories are formulated to describe. In other words, macroscopic plasticity theories describe processes that involve changes in defect structures not explicitly accounted for by macroscopic strain fields (phenomenological, with much fewer numbers of DOF compared with MD models) alone. Such processes do not return a body's defect structure to its original state even though the macroscopic strain returns to zero, giving rise to the continuum phenomenological dissipation. A cycle of macroscopic plastic strain must not be equated with the return of atomic structures at the MD level. Realistic MD models at the size and time-scales of processes considered by continuum plasticity theories should capture evolution and irreversibility of defect structures. The equivalent continuum is a faithful representation of the MD system, in terms of work rates, momentum, mass, as well as energy. It naturally characterizes dissipation embedded in the MD models. One of the issues in the grand challenge of multiscale modelling of material behaviour is to devise continuum representations of MD results at different time- and length-scales with appropriately reduced numbers of DOF. One requirement for such developments is accurate reflection of dissipation embedded in the MD results and the EC and TMEC representations developed here.

The decomposition of the atomic level rate of deformation * D* and into an elastic part and a plastic part in the forms of(4.1)is an issue somewhat separate from the formulation of the EC and TMEC theories discussed here. Further developments are needed in order for insights to be gained in this regard. Although it is an independent issue in its own right, it could be handled together with the need to formulate higher scale thermoplastic characterizations of the phenomenological nature. Such an endeavour might include the reduction of the number of DOF and the use of lower resolution state variables including a phenomenological rate of deformation of the form(4.2)

In such a scenario, the first correspondence of equation (2.35) and the third correspondence in equation (2.43) may be invoked to estimate the fraction of stress work converted to heat as(4.3)

This quantity includes the contributions of both thermoelastic dissipation and thermoplastic dissipation. It is more meaningful when interpreted in a proper time-averaged sense, similar to the consideration in equation (3.6). The thermoelastic dissipation rate (*β*^{elas}) is relatively small compared with the thermoplastic rate. To quantify these two dissipation rates, two calculations can be carried out. The first would involve only elastic deformation of an atomic system through which the thermal elastic rate (*β*^{elas}) can be determined. The second calculation should be carried out under the conditions of interest, involving elastic and plastic deformations. The plastic dissipation rate can be determined through *β*^{plas}=*β*−*β*^{elas}. The determination of these quantities are useful in formulating a higher-scale phenomenological decomposition in the form of equation (4.2) when has a much smaller number of DOF than . This is because one basic requirement for the phenomenological plastic rate of deformation is that it yields the fraction of plastic work rate through(4.4)Here, represents the fraction of the higher-scale phenomenological plastic work converted to heat, as measured by, for example, Taylor & Quinney (1937) and Rosakis *et al*. (2000) for metals. Clearly, when *V*^{(e)} is sufficiently large and an average over a sufficiently long period of time is taken, equation (4.4) and knowledge of can be used to formulate a more realistic decomposition represented by equation (4.2). Conversely, if structural considerations dictate the decomposition in equation (4.2) at a particular scale, equation (4.4) and *β*^{plas} can be used to predict through MD calculations.

In addition to thermoelastic and thermoplastic dissipations, equations (2.42) and (2.43) fully account for heat generation through other mechanisms, such as electromagnetic heating. Note that the heat source term depends on body force density , which accounts for the effect of . Of course, high-frequency electromagnetic forces responsible for electromagnetic heating contribute to .

## 5. Discussion

### (a) The atomic origin of coupled thermomechanical behaviour

The TMEC theory for dynamically deforming atomistic systems developed in §2 is an extension of the pure mechanical theory of equivalent continuum described in Zhou (2003). This novel and fundamental theory establishes, within the meaning of classical mechanics, the ultimate atomic origin of coupled thermomechanical deformation processes. It provides a framework for the formulation of coupled thermomechanical continuum descriptions of material behaviour at any spatial and temporal scale ‘from the ground up’, using molecular dynamics simulations. It should be pointed out that the analysis applies to all classical molecular dynamics models regardless of, for example, whether a thermostat is used or whether a macroscopic constant temperature (or pressure) is maintained. This theory is based on a decomposition of the atomic particle velocity into a relatively slower-varying structural deformation part and a high-frequency thermal oscillation part. The continuumization of the discrete molecular dynamics model is carried out with strict adherence to balance of momentum, balance of energy and conservation of mass, in a fully dynamic and time-resolved manner. The establishment of equivalence between the discrete system and the continuum system follows the same approach used for the pure mechanical theory. This approach uses the dynamic principle of virtual work. Since the decomposition of atomic velocity into a structural deformation part and a thermal oscillation part is intrinsically dependent on both size- and time-scales, the development has followed a general framework of analysis, allowing scale-sensitive characterizations of both the structural deformation as it appears to a macroscopic observer and the thermal oscillation as it appear to an observer moving at the structural velocity of an atom. On one hand, balance of momentum at the structural level yields fields of stress, body force, traction, mass density and deformation that maintain full dynamic equivalence between the discrete system and continuum system. This equivalence includes (i) preservation of linear and angular momenta, (ii) conservation of internal, external and inertial work rates and (iii) conservation of mass. On the other hand, balance of momentum in terms of the thermal motions yields the fields of heat flux and temperature. These quantities have been cast in a manner as to conform to the continuum phenomenological equation for heat conduction and heat generation, allowing the specific heat, thermal conductivity and relaxation time at different scales to be quantified. The coupling of the thermal equation of heat conduction and the mechanical equation of structural deformation, as is phenomenologically known at higher continuum scales, occurs because the balance equations at the structural level and the thermal level are only two different forms of the same balance of momentum equation at the fully time-resolved atomic level. This coupling occurs through an inertial force term in each of the two equations ( in equation (2.32) and in equation (2.38)) induced by the other process. For the structural deformation equation, the inertial force term induced by thermal oscillations of atoms gives rise to the phenomenological dependence of deformation on temperature. This point can be seen from the first correspondence in equation (2.35). For the heat equation, the inertial force term induced by structural deformation takes the phenomenological form of a heat source, see equation (2.43). Equations (2.32) and (2.38) represent conservation of energy as it appears to, respectively, the macroscopic observer (who sees ) and the travelling observer (who sees ). Taken together, these two equations specify conservation of energy for the atomic system in the fully time-resolved sense. In this sense, they are equivalent to eqn. (5.3) of Zhou (2003). In summary, the TMEC theory developed here allows scale-dependent continuum descriptions of material behaviour in the form of coupled thermomechanical processes of deformation and heat flow to be obtained from purely mechanical molecular dynamics simulations at the atomic level.

The framework here fully recognizes the spatial scale dependence of deformation. The choice of element *V*^{(e)} is totally arbitrary, allowing size dependence of atomic behaviour to be analysed. Specifically, *V*^{(e)} can be as small as the tetrahedron defined by four neighbouring atoms and as large as the whole system *V*. Length-scale effects due to lattice spacing, thermal and structural wavelengths and structural features such as voids and grains can be explicitly analysed as the size of *V*^{(e)} is varied. Such analyses, of course, carry a very significant computational cost.

The theoretical approach also fully admits the analysis of time-scale influences. Specifically, the decomposition of atomic motion in equation (2.1) is highly dependent on time resolution. At the macroscopic time-scale, high frequency components of atomic oscillations are not explicitly resolvable and are therefore included in . As the time resolution is increased, relatively lower frequency components of atomic oscillations become explicitly resolved and are included in , causing to decrease. In the limit of an analysis with full resolution for atomic motion, and . Consequently, the pure mechanical theory of equivalent continuum in Zhou (2003) is recovered. Of course, temporal and spatial resolutions are often related and cannot be fully separated. Temporal resolutions are higher at small scales. For example, at the size scale of nanowires or electronic components, structural waves along the axial direction may need to be analysed, causing the long wavelength (of the order of 1–2 nm) components of atomic motion to be included in . Such wave components may be regarded as thermal waves at the macroscopic level.

### (b) Constitutive behaviour

It is important to point out that the TMEC representation developed here and any phenomenological characterizations of it at higher scales (not the topic of this paper) carry constitutive information of the material analysed through the thermal and mechanical fields. Specifically, the constitutive behaviour is implied in the histories of deformation. Therefore, the continuum fields obtained here at various time- and size-scales can be used to develop constitutive characterizations at the relevant scales. In this sense, the theoretical framework here is highly relevant to the multiscale analysis of material behaviour. We note that the nonlinearity of the interatomic potential is embedded in the thermal and mechanical fields and their histories. For the mechanical part, the stress evolves nonlinearly with finite deformation. For the thermal part, the nonlinearity may manifest in the form of the temperature dependence of specific heat and thermal conductivity even at equilibrium.

### (c) Relation to statistical mechanics

The time- and structure-resolved approach of analysis here is fully dynamic in nature. The field quantities, especially the thermal quantities, can be expected to approach macroscopically measured quantities when proper averages in both time and space are obtained at the appropriate levels. This approach to reach higher scales from atomic models is distinct from thermodynamic approaches or statistical mechanics. The differences are (i) the TMEC theory is fully dynamic and does not make any assumption regarding a system's proximity to macroscopic equilibrium or the lack of it; (ii) the analysis is carried out on the current (deformed) configuration of a system, therefore, large deformations and defect activities (e.g. dislocation motion) can be analysed; (iii) it permits the analysis of non-equilibrium thermomechanical behaviour during large and transient deformations; (iv) non-uniform full-field characterizations of all relevant quantities (e.g. stress, strain, strain rate, temperature and mass density) are developed, admitting effects of fluctuations, non-local interactions and scale-dependence; and (v) it maintains full equivalence between the discrete molecular system and the continuum representation in terms of momentum, energy and work conjugacy of relevant continuum fields. Another way to look at the difference between statistical mechanics and the TMEC theory here is to note that, while the former captures the behaviour at the ensemble level by focusing on macro state averages (average stress, average temperature), the latter takes a microcanonical approach by developing a full field representation and by capturing local fluctuations in field quantities. The TMEC fields are piecewise continuous, admitting large fluctuations at small size scales. This is an intrinsic feature of the dynamic equivalence at interatomic scales which allows representation of the effects of atomic scale material heterogeneity with relatively high fidelity. The resolution of such interatomic features is important for problems involving steep gradients of the fields and heterogeneities, such as dislocations, interstitial atoms and vacancies. Of course, the fields obtained here can lend itself to continuum treatments, including averaging which constitutes a separate but important task having also to do with transition to higher length scales. In short, the TMEC theory is a novel alternative equivalent continuum approach which, through strict fidelity to the atomic system in terms of full dynamic equivalence and through extension of the statistical concept of heat to individual atoms or small clusters of atoms, offers a framework for explicit analysis of nanoscopic deformation and transport mechanisms in a manner not necessarily possible with traditional statistical means. In particular, the approach may prove important in settings involving short time, transient responses and nanoelectronic devices.

It is also important to point out that the TMEC fields can be used to obtain the appropriate free energy and entropy in a manner consistent with statistical mechanics under relevant conditions. For example, the Helmholtz free energy *F* of can be defined in the same way as in statistical mechanics through (Born & Huang 1988)(5.1)where *E* is the internal energy, *T* is the temperature evaluated from equation (2.45), is the entropy, *U* is the potential energy of the subsystem in evaluated using the macroscopically observed atomic positions , *h* is Planck's constant and is the frequency of the free vibration of the *α*th degree of freedom (*α*=1, 2 or 3) of atom *I*. To evaluate the incremental relation , note that(5.2)

### (d) Other mechanisms for heat conduction

Finally, it is important to point out that the starting point of the theoretical development in this paper is molecular dynamics. The framework used here naturally inherits the assumptions and limitations of molecular dynamics. In particular, we note that the effect of free electrons on heat conduction in metals and alloys is not explicitly accounted for in standard MD models. Reflecting this characteristic of MD, the heat conduction analysed in the TMEC theory developed in this paper concerns solely the contribution of atomic vibrations (or phonons) to thermal conductivity. In order to account for the contributions of free electrons (metals, alloys and semiconductors; e.g. Berman 1976), of electron/phonon interactions, and of bipolar carriers consisting of electron/hole pairs at high temperatures (semiconductors; e.g. Bhandari & Rowe 1988) to thermal conductivity, modifications to the MD framework are required, if a framework within classical mechanics is to be used. This would constitute a significant extension of molecular dynamics which is possible and, perhaps, worthwhile.

## Acknowledgments

Support from the NSF through CAREER award number CMS9984289, from the NASA Langley Research Center through grant NAG-1-02054, and from an AFOSR MURI at Georgia Institute of Technology is gratefully acknowledged.

## Footnotes

As this paper exceeds the maximum length normally permitted, the authors have agreed to contribute to production costs.

- Received May 19, 2004.
- Accepted February 7, 2005.

- © 2005 The Royal Society