Nonequilibrium thermomechanics of Gaussian phase packet crystals: application to the quasistatic quasicontinuum method
NNonequilibrium thermomechanics of Gaussian phase packet crystals:application to the quasistatic quasicontinuum method
Prateek Gupta, Dennis M. Kochmann
Mechanics & Materials Lab, Department of Mechanical and Process EngineeringETH Z¨urich, 8092 Z¨urich, Switzerland
Abstract
The quasicontinuum method was originally introduced to bridge across length scales – from atomistics to signifi-cantly larger continuum scales – thus overcoming a key limitation of classical atomic-scale simulation techniqueswhile solely relying on atomic-scale input (in the form of interatomic potentials). An associated challenge liesin bridging across time scales to overcome the time scale limitations of atomistics. To address the biggest chal-lenge, bridging across both length and time scales, only a few techniques exist, and most of those are limited toconditions of constant temperature. Here, we present a new strategy for the space-time coarsening of an atom-istic ensemble, which introduces thermomechanical coupling. We investigate the quasistatics and dynamics ofa crystalline solid described as a lattice of lumped correlated Gaussian phase packets occupying atomic latticesites. By definition, phase packets account for the dynamics of crystalline lattices at finite temperature throughthe statistical variances of atomic momenta and positions. We show that momentum-space correlation allows foran exchange between potential and kinetic contributions to the crystal’s Hamiltonian. Consequently, local adia-batic heating due to atomic site motion is captured. Moreover, in the quasistatic approximations the governingequations reduce to the minimization of thermodynamic potentials such as Helmholtz free energy (depending onthe fixed variables), and they yield the local equation of state. We further discuss opportunities for describingatomic-level thermal transport using the correlated Gaussian phase packet formulation and the importance ofinteratomic correlations. Such a formulation offers a promising avenue for a finite-temperature non-equilibriumquasicontinuum method that may be combined with thermal transport models.
Keywords:
Quasicontinuum, Multiscale Modeling, Atomistics, Non-Equilibrium Statistical Mechanics,Updated-Lagrangian
1. Introduction
Crystalline solids exhibit physical and chemical transport phenomena across wide ranges of length and timescales. This includes the transport of charges (Butcher, 1986; Ziman, 2001), of heat (Ziman, 2001), and ofmass (Weiner, 2012), as well as mechanical failure. Understanding such phenomena is crucial from both a fun-damental scientific standpoint as well as to further advance technologies ranging form solid-state batteries (Kimet al., 2014) to thermal management systems (Hicks and Dresselhaus, 1993) to failure-resistant metallic struc-tural components (Hirth, 1980) – all exposed to complex dynamic conditions. The variety of length and timescales underlying such transport phenomena mandates the need for simulation techniques that capture the phys-ical processes at all length and time scales involved. While continuum mechanics and related finite elementand phase field methods have been successful at modeling physical processes at relatively large length scales(typically micrometers and above) and time scales (milliseconds and above) (Hirth, 1980; Mendez et al., 2018),molecular statics (MS) and molecular dynamics (MD) have been successful at elucidating the physics of varioustransport phenomena at atomic-level length scales (angstroms to tens of nanometers) and time scales (femto- tonanoseconds) (Tuckerman, 2010). Where a higher level of accuracy is required, methods such as Density Func-tional Theory (DFT) or the direct computation of Schr¨odinger’s equation have aimed at capturing the quantumcoupling of molecular-level physics, tracking particle trajectories at subatomic length scales (on the order of − − − m) and time scales ( − − − s). All of the aforementioned techniques specialize in the ap-proximate ranges of length and time scales mentioned above. However, each of those techniques poses restrictive Email addresses: [email protected] (Prateek Gupta), [email protected] (Dennis M. Kochmann)
Preprint submitted to Journal of the Mechanics and Physics of Solids January 29, 2021 a r X i v : . [ c ond - m a t . m e s - h a ll ] J a n ssumptions at smaller scales while becoming prohibitively costly at larger scales, hence making scale-bridgingtechniques attractive (Srivastava and Nemat-Nasser, 2014; van der Giessen et al., 2020).The quasicontinuum (QC) method of Tadmor et al. (1996a) is one such method that aims to solve the problemof spatial scale-bridging from atomistic to continuum length scales via intermediate mesoscopic scales (Milleret al., 1998). Starting from a standard atomistic setting, a carefully selected set of representative degrees of free-dom ( repatoms ) reduces the computational costs and admits simulating continuum-scale problems by restrictingatomistic resolution to where it is in fact needed. Different flavours of QC have been proposed based on theinterpolation of forces on repatoms (Knap and Ortiz, 2001), or based on approximating the total Hamiltonianof the system using quadrature-like summation and sampling rules (Eidel and Stukowski, 2009; Dobson et al.,2010; Gunzburger and Zhang, 2010; Espanol et al., 2013). For example, the nonlocal energy-based formulationof Amelang et al. (2015) enables a seamless spatial scale-bridging within the QC setup, which does not requirea strict separation of (nor a-priori knowledge about) atomistic and non-atomistic (coarsened) subdomains withinthe simulation domain. The capabilities of this fully nonlocal QC technique have been demonstrated, e.g., bylarge-scale quasistatic total-Lagrangian simulations of dislocation interactions during nano-indentation (Ame-lang et al., 2015), nanoscale surface and size effects (Amelang and Kochmann, 2015), and void growth andcoalescence (Tembhekar et al., 2017). In this work, we adopt the nonlocal energy-based setup of Amelang et al.(2015) in a new, updated Lagrangian setting for spatial upscaling.While spatial coarse-graining is thus achieved by approximating an ensemble of atoms with a subset ofatoms, temporal coarse graining requires approximate modeling due to the unavailability of the trajectories of allthe atoms, at a given time. Furthermore, Hamiltonian dynamics of an ensemble of atoms couples length and timescales in the system (Evans and Morriss, 2007), which is why spatial scale-bridging techniques have often beenapplied to systems at zero temperature or at uniform temperature only (Tadmor et al., 2013). One way to modeluniform temperature is to apply a global ergodic assumption for every atomic site, thus yielding space-averagedtrajectories of atoms as phase-averaged trajectories. Suitable global thermal equilibrium distributions (such asNVT ensembles (Tadmor and Miller, 2011)) are used for phase averaging of trajectories. However, most priorwork has been based on local harmonic approximations of interatomic potentials (Tadmor and Miller, 2011;Tadmor et al., 2013). Another way is to use Langevin dynamics, in which a stochastic thermal forcing is addedto the dynamics of atoms to account for thermal vibrations (Qu et al., 2005; Marian et al., 2009). Unfortunately,such an approach poses a time integration restriction even for systems in thermodynamic equilibrium (uniformtemperature), which is why it is computationally costly. Kulkarni et al. (2008) introduced a variational approachfor modeling non-uniform temperature, which approximates the global distribution function of the ensemble asa product of local entropy maximizing (or max-ent ) distributions, constraining the local frequency of atomicvibrations using the local equipartition of energy of every atom. This local ergodic assumption yields time-averaged trajectories as phase-averaged trajectories and thus enables the definition of non-uniform temperatureand internal energy. However, the interatomic independence or statistically uncorrelated local distributions,inherent in that approach, precludes the transport of energy across the atoms. As a remedy, Kulkarni et al.(2008) modeled thermal transport using the finite-element setup of the QC formulation and empirical bulk-thermal conductivity values. Venturini et al. (2014) extended the max-ent approach to non-uniform temperaturesand non-uniform interstitial concentrations of solute atoms in crystalline solids to also include mass diffusion.Specifically, they modeled transport using linear Onsager kinetics, derived from a local dissipation inequality.Here, we introduce a different formulation based on Gaussian Phase Packets (GPP) from the max-ent distributionapproximation to understand the dynamics of local distribution functions of atoms.Previously, Ma et al. (1993) studied the dynamic Gaussian Phase Packet (GPP) approximation as a tra-jectory sampling technique for uniform-temperature molecular dynamics simulations. They approximated thedistribution function of each atom in an ensemble as a correlated Gaussian distribution with the covariance ma-trix as the mean-field or phase-space parameter. Such an approximation yields the evolution equations of thecovariance matrix by either directly substituting the Liouville equation (strong form) or by using the Frenkel-Dirac-McLachlan variational principle (McLachlan, 1964). The resulting equations may be integrated in time,combined with appropriate ergodic assumptions, to infer the locally averaged physical quantities of the system(such as temperature). However, Ma et al. (1993) applied the formulation on a small system of atoms relaxingtowards thermodynamic equilibrium. Their work was inspired by the application of GPPs in quantum mechanicsby Heller (1975), who used it for calculating parameterized solutions of the Schr¨odinger’s equation.In this work, we use the GPP approximation applied to an ensemble of atoms to elucidate the nonequilibriumand (local-thermal) equilibrium thermomechanical quasistatics and dynamics of the system. We show that an2pproximation of interatomic correlations is required for modeling atomistic-level transport phenomena. In anensemble of N atoms, such an approximation increases the degrees of freedom by O ( N ) , which is computa-tionally costly. Therefore, we employ uncorrelated atoms or interatomic independence for further modeling andapplications. In addition, we combine the GPP framework within the quasistatic approximation with Venturiniet al.’s linear Onsager kinetics to model local thermal transport. In Section 2 we review the fundamentals ofnonequilibrium statistical mechanics based on the Liouville equation and the GPP approximation. We show thatthe interatomic correlations are of fundamental importance for modeling interatomic heat flux. However, suchcorrelations tend to increase the degrees of freedom significantly, which is why they are neglected in this work.However, we retain the phase-space correlation of each atom, which we later identify as the thermal momentum .In Section 3 we discuss the importance of thermal momentum in dynamics and the implications of its vanishinglimit in quasistatics. We show that in the quasistatic limit the equations define local thermomechanical equilib-rium of the system and yield the rate-independent local thermal equation of state. To capture thermal transport,we review Venturini et al.’s linear Onsager kinetics-based model and its application. We also discuss the timescale imposed by the Onsager kinetics-based model and its time stepping constraints. In Section 4 we discussthe QC implementation of the local thermomechanical equilibrium equations combined with thermal transportand demonstrate its use in a new update-Lagrangian distributed-memory QC solver for coarse-grained atomisticsimulations. Finally, in Section 5 we conclude our analysis and discuss limitations and future prospects.
2. Nonequilibrium thermodynamics of Gaussian Phase Packets
In this section, we discuss the nonequilibrium modeling of deformation mechanics of crystalline solids us-ing the GPP approximation in which atoms are treated as Gaussian clouds, centered at the mean phase-spacepositions of the atoms. We briefly review the application of the Liouville equation for analyzing the statisti-cal evolution of a large ensemble of atoms subject to high-frequency vibrations. Such random vibrations aremodeled by assuming local phase-space coordinates (positions and momenta) drawn from local Gaussian dis-tributions. We discuss the impact of assuming independent Gaussian distributions for each atom (termed inter-atomic independence hereafter) and the corresponding inability of the model to capture nonequilibrium thermaltransport. Moreover, the independent Gaussian assumption allows us to formulate the dynamical system ofequations governing the atoms and the corresponding mean-field parameters. In subsequent sections, we use theinsights gained here to formulate isentropic and non-isentropic quasistatic problems of finite-temperature crystaldeformation (see Figure 1 for a schematic description).
Figure 1: Schematic illustrating a typical nonequilibrium QC study of a domain with atomistic and coarsened subdomains.All atomic sites are modeled using the Gaussian Phase Packet (GPP) approximation. The transport of energy among theatomic sites is modeled using the linear Onsager kinetics of Venturini et al. (2014). .1. Hamiltonian dynamics The time evolution of an atomic crystal, modeled as an ensemble of classical particles, is fully characterizedby the generalized positions q = { q i ( t ) , i = 1 , . . . , N } and momenta p = { p i ( t ) , i = 1 , . . . , N } of all N particles, which evolve in time according to the Hamiltonian equations, d p i d t = − ∂ H ∂ q i = −∇ q i V ( q ) = F i ( q ) , (2.1a) d q i d t = ∂ H ∂ p i = p i m i , (2.1b)where H ( p , q ) = N (cid:88) i =1 p i · p i m i + V ( q ) (2.2)is the Hamiltonian of the system, V ( q ) represents the potential field acting on all atoms in the system, and m i is the mass of the i th atom. Equations (2.1) are solved in standard molecular dynamics studies of crystals, inwhich trajectories ( p i ( t ) , q i ( t )) of all atoms are resolved in time. As a consequence, such simulations requirefemtosecond-level time steps (Tuckerman, 2010; Tadmor et al., 1996a,b) and are unable to capture long-time-scale phenomena within computationally feasible times.To this end, we consider the statistical treatment of the Hamiltonian dynamics governed by (2.1) and (2.2),in which the local vibrations of all atoms about their mean positions are modeled as random fluctuations in thephase-space coordinate z = ( p ( t ) , q ( t )) . Here and in the following, we use z ∈ R N for brevity to represent themomenta and positions of the N particles in three dimensions (3D). It is convenient to introduce the distribution f ( z , t ) such that the quantity d P = f ( p , q , t ) N (cid:89) i =1 d p i N (cid:89) i =1 d q i (2.3)is the probability of finding the system of atoms such that the position and momentum of the i th atom lie within ( q i , q i + d q i ) and ( p i , p i + d p i ) , respectively. The probability distribution f ( z , t ) is governed by the Liouvilleequation (Evans and Morriss, 2007; Tadmor and Miller, 2011) ∂f ( z , t ) ∂t + i L f = 0 with f ( z ,
0) = f ( z ) , lim z →∞ f ( z , t ) = 0 , (2.4)with the Liouville operator L defined by i L = ∂∂ z · ˙ z + ˙ z · ∂∂ z (2.5)and ˙ z = ( ˙ p , ˙ q ) . Here and in the following, dots denote time derivatives. Given the equations of motion in (2.1)and the Hamiltonian of the system in (2.2), the Liouville operator in (2.5) becomes i L = ˙ z · ∂∂ z + ∂∂ z · ˙ z = ˙ p · ∂∂ p + ˙ q · ∂∂ q + ∂∂ p · ˙ p + ∂∂ q · ˙ q = − ∂ H ∂ q · ∂∂ p + ∂ H ∂ p · ∂∂ q . (2.6)We note that solution of (2.4) combined with (2.6) at all times is equivalent to solving the equations of mo-tion (2.1) with (2.2) (Evans and Morriss, 2007). In general, such a solution requires a discretization of the N -dimensional phase-space { Γ ⊆ R N : z ∈ Γ } and of time for a system of N atoms in 3D. To avoid such acomputationally intensive discretization, we parametrize f ( z , t ) using the GPP approximation (detailed below),which yields the equations of motion of the mean phase-space coordinates and the respective fluctuation auto-and cross-correlations, categorized as phase-space or mean-field parameters. The GPP approximation was first introduced in the context of quantum mechanics by Heller (1975) andused for classical systems by Ma et al. (1993). Both Heller (1975) and Ma et al. (1993) substituted Gaussiandistributions into the Liouville equation (2.4) (Schr¨odinger’s equation for quantum systems) to obtain the dy-namical evolution of the phase-space parameters. Furthermore, Ma et al. (1993) noted that identical relations4an be obtained using the Frenkel-Dirac-McLachlan variational principle for the Liouville equation, commonlyused for deriving equations of motion of phase-space parameters in quantum mechanics (McLachlan, 1964).Following their idea, we approximate the system-wide probability distribution f ( z , t ) as a multivariate Gaussiandistribution, i.e., f ( z , t ) = 1 Z ( t ) e − ( z − z ( t )) T Σ − ( t )( z − z ( t )) , (2.7)where Σ is the N × N covariance matrix composed of the interatomic and momentum-displacement corre-lations, z represents the vector of all atoms’ mean positions and momenta, and the partition function Z ( t ) isdefined by Z ( t ) = 1 N ! h N (cid:90) R N e − ( z − z ( t )) T Σ − ( t )( z − z ( t )) d z = (2 π ) N N ! h N √ det Σ . (2.8)The phase average of any function A ( z ) is denoted by (cid:104) A ( z ) (cid:105) and defined as (cid:104) A ( z ) (cid:105) = 1 N ! h N (cid:90) R N A ( z ) f ( z , t ) d z , (2.9)where h is the Planck’s constant and the factor N ! h N is the normalizing factor for the phase-space volume (Lan-dau and Lifshitz, 1980). This confirms that z ( t ) = (cid:104) z ( t ) (cid:105) . We further conclude that Σ can be written as ablock-matrix with components Σ ij = 1 N ! h N (cid:90) R N ( z i − z i ) ⊗ ( z j − z j ) f ( z , t ) d z = (cid:104) ( z i − z i ) ⊗ ( z j − z j ) (cid:105) , (2.10)such that each block represents the correlation among displacements and momenta of atoms i and j . Conse-quently, we obtain the phase-space parameters ( z ( t ) , Σ ( t )) . That is, from now on we track the atomic ensemblethrough the mean positions and momenta of all atoms, z ( t ) , and their covariance matrix, Σ ( t ) – rather thantime-resolving the positions and momenta directly.Time evolution equations for the phase-space parameters are obtained from the Liouville equation (2.4) byusing the following identity for the phase average of any phase function A ( z ) (Evans and Morriss, 2007) d (cid:104) A (cid:105) d t = 1 N ! h N (cid:90) R N f ( z , t ) d A d t d z = (cid:28) d A d t (cid:29) (2.11)(see Appendix A for a brief discussion). Application to z ( t ) and Σ ( t ) yields the dynamical equations d z d t = (cid:104) ˙ z (cid:105) , d Σ ij d t = (cid:10)(cid:0) ˙ z i − ˙ z i (cid:1) ⊗ ( z j − z j ) (cid:11) + (cid:10) ( z i − z i ) ⊗ (cid:0) ˙ z j − ˙ z j (cid:1)(cid:11) . (2.12)Let us further specify the second equation. Writing the components of covariance matrix Σ ij as Σ ij = (cid:32) Σ ( p , p ) ij Σ ( p , q ) ij Σ ( q , p ) ij Σ ( q , q ) ij (cid:33) , (2.13)and assuming that all atoms have the same mass m , identity (2.11) yields the following time evolution equations: d Σ ( p , p ) ij d t = (cid:10) F i ( q ) ⊗ (cid:0) p j − p j (cid:1)(cid:11) + (cid:104) ( p i − p i ) ⊗ F j ( q ) (cid:105) , (2.14a) d Σ ( p , q ) ij d t = (cid:10) F i ( q ) ⊗ (cid:0) q j − q j (cid:1)(cid:11) + Σ ( p , p ) ij m , (2.14b) d Σ ( q , p ) ij d t = (cid:104) ( q i − q i ) ⊗ F j (cid:105) + Σ ( p , p ) ij m . (2.14c) d Σ ( q , q ) ij d t = Σ ( p , q ) i + Σ ( q , p ) i m , (2.14d)5quations (2.14) govern the evolution of the pairwise momentum and displacement correlations of atoms i and j , and they must be solved to obtain the interatomic correlations at all times. Equations (2.14c)-(2.14d) gov-ern the thermomechanical coupling of the crystal. Note that we may identify the pairwise kinetic tensor Σ ( p , p ) ij as a measure of temperature, so that (2.14b) and (2.14c) describe the evolution of the system-wide distributionas a result of unbalanced pairwise virial tensors and kinetic tensors (see Admal and Tadmor (2010) for the ten-sor virial theorem). The virial tensor (cid:10) F i ( q ) ⊗ (cid:0) q j − q j (cid:1)(cid:11) changes with changing displacement correlations ofatoms due to varying extents of atomic vibrations Σ ( q , q ) ij , thus coupling (2.14d) with (2.14b) and (2.14c). Theright-hand side of (2.14a) resembles the tensor form of the interatomic heat current (S¨a¨askilahti et al., 2015;Lepri et al., 2003) and changes with varying correlation matrices Σ q , p ij , thus coupling (2.14a) with (2.14b) and(2.14c). Consequently, the imbalance between virial and kinetic tensors in equations (2.14b) and (2.14c) drivesthe phase-space motion of the system of particles, resulting in the time evolution of Σ ( p , p ) ij and Σ ( q , q ) ij .We identify the entropy S of the atomic ensemble as S = − k B (cid:104) ln f (cid:105) = k B (cid:32) N [1 + ln(2 π )] − ln ( N !) + ln (cid:32) √ det Σ h N (cid:33)(cid:33) = S + k B ln (cid:32) √ det Σ h N (cid:33) , (2.15)where k B is Boltzmann’s constant, and S is a constant for a given system with a constant number of atoms. Theentropy rate of change follows as d S d t = k B Σ d(det Σ )d t . (2.16)Overall, equations (2.12) govern the phase-space motion of a system of atoms and contain N + 36 N ( N +1) / equations, solving which is even more computationally expensive than solving the state-space governingequations of MD. As a simplifying assumption, Ma et al. (1993) and Heller (Heller, 1975) assumed the statisticalindependence of atoms (i.e., of states in the quantum analogue), which implies Σ ij = for i (cid:54) = j. (2.17)As shown in the following section, the assumption (2.17) severely limits the applicability of (2.12), since underthis assumption the transport of heat is not correctly resolved in the evolution equations. Consequently, we mayapply the phase-space evolution equations (2.12) and (2.14) only for isentropic (reversible) finite temperaturesimulations of quasistatic and dynamic processes. Since solving the full system of evolution equations is prohibitively expensive, as discussed above, let usapply (2.17) and assume non-zero correlations between the position and momentum of each individual atom, butno cross-correlations between different atoms. To this end, we apply the GPP approximation to a single atom i ,which yields the multivariate Gaussian distribution of the phase-space coordinate z i as f i ( z i , t ) = 1 Z i e − ( z i − z i ) T Σ − i ( z i − z i ) so that f ( z , t ) = N (cid:89) i =1 f i ( z i , t ) , (2.18)where the phase-space parameters ( z i , Σ i ) denote the mean phase-space coordinate and variance of the i th atom,respectively, defined as z i ( t ) = 1 h (cid:90) R f i ( z i , t ) z i d z i = (cid:104) z i (cid:105) , Σ i = (cid:104) z i ⊗ z i (cid:105) . (2.19)The normalization quantity Z i ( t ) may be identified as the single particle partition function Z i ( t ) = 1 h (cid:90) R e − ( z i − z i ) T Σ − i ( z i − z i ) d z i = (cid:18) πh (cid:19) (cid:112) det Σ i . (2.20) Σ i is the × covariance matrix of the multivariate Gaussian and accounts for the variance or uncertainty in themomentum p i and displacement q i of the i th atom.The assumed interatomic independence eliminates the interatomic correlations as independent variables, thusreducing the total number of equations to N +6 N for a system of N atoms, which govern the time evolution of6 igure 2: Illustration of an atomic trajectory by GPP dynamics (on the right) compared to a molecular dynamics trajectory(on the left) for a single particle. Small-scale motions of the particle are approximated by the parameters Ω , Σ and β uponaveraging over suitable time intervals. As discussed in Section 3.2, in the quasistatic limit, Ω and Σ are proportional to thelocal temperature. the kinetic tensor Σ ( p , p ) i , displacement-correlation tensor Σ ( q , q ) , and the momentum-displacement-correlationtensor Σ ( p , q ) i = (cid:16) Σ ( q , p ) i (cid:17) T via d Σ ( p , p ) i d t = (cid:104) F i ( q ) ⊗ ( p i − p i ) (cid:105) + (cid:104) ( p i − p i ) ⊗ F i ( q ) (cid:105) , (2.21a) d Σ ( p , q ) i d t = (cid:104) F i ( q ) ⊗ ( q i − q i ) (cid:105) + Σ ( p , p ) i m , (2.21b) d Σ ( q , q ) i d t = Σ ( p , q ) i + Σ ( q , p ) i m , (2.21c)combined with the phase-averaged equations of motion, d (cid:104) p (cid:105) i d t = (cid:104) F (cid:105) i , d (cid:104) q (cid:105) i d t = (cid:104) p (cid:105) i m . (2.22)For simplicity, we further make the spherical distribution assumption that all cross-correlations between differentdirections of momenta and displacements vanish, hence approximating the atomic clouds in phase-space asspherical (thus eliminating any directional preference of the atomic vibrations). While being a strong assumption,this allows us to reduce the above tensorial evolution equations to scalar ones. Specifically, taking the trace tr ( · ) of (2.21), we obtain, dΩ i d t = 2 (cid:104) F i ( q ) · ( p i − p i ) (cid:105) , (2.23a) dΣ i d t = 2 β i m , (2.23b) d β i d t = (cid:104) F i ( q ) · ( q i − q i ) (cid:105) i m , (2.23c)where we introduced the three scalar parameters tr (cid:16) Σ ( p , p ) i (cid:17) = 3Ω i , tr (cid:16) Σ ( q , q ) i (cid:17) = 3Σ i and tr (cid:16) Σ ( p , q ) i (cid:17) = 3 β i . (2.24)Equations (2.23) are N coupled scalar ODEs, which determine the changes in the vibrational widths of atomsin phase-space ( Ω i and Σ i ) and the correlation β i between the displacement and momentum vibrations of the i th atom (see Figure 2). We note that (2.23) are identical to the equations used by Ma et al. (1993), who used theformulation as an optimization procedure for a system of atoms at a uniform constant temperature.The physical role of momentum-displacement correlation β i becomes evident upon applying a time-reversaltransformation t (cid:55)→ − t to (2.22) and (2.23), which results in the transformations ( q i , Ω i , Σ i ) (cid:55)→ ( q i , Ω i , Σ i ) and ( p i , β i ) (cid:55)→ ( − p i , − β i ) . Consequently, the correlation β i signifies the momentum of the i th atom in phase-space,governing the time evolution of the thermal coordinate Σ i . Hence, β i will be referred to as the thermal momentum i = 1 , . . . , N d (cid:104) q (cid:105) i d t = (cid:104) F (cid:105) i m , (2.25a) d Σ i d t = 2Ω i m + 2 (cid:104) F i ( q ) · ( q i − q i ) (cid:105) m , (2.25b) dΩ i d t = 2 (cid:104) F i ( q ) · ( p i − p i ) (cid:105) . (2.25c)Time reversibility of (2.25) highlights that these evolution equations do not capture nonequilibrium irreversiblethermal transport at all scales due to the simplification of (2.14) to (2.21) under the interatomic independenceassumption (2.17). As mentioned previously, such interatomic independence is necessary to keep a feasiblenumber of unknowns for large ensembles. The reversible entropy fluctuations can be obtained by substitutingthe independent GPP assumption into (2.16), which gives d S d t = k B Σ d det Σ d t = N (cid:88) i =1 (cid:18) k B Σ i d det Σ i d t (cid:19) = N (cid:88) i =1 d S i d t , (2.26)where the local entropy fluctuations of the i th atom are given by d S i d t = k B Σ i d det Σ i d t = k B (cid:0) Ω i Σ i − β i (cid:1) dd t (cid:0) Ω i Σ i − β i (cid:1) = k B (cid:0) Ω i Σ i − β i (cid:1) (cid:18) i β i m (cid:2) (Ω i Σ i ) − β i (cid:3) + 2 (cid:2) Ω i Σ i F i · ( p i − p i ) − β i F i · ( q i − q i ) (cid:3)(cid:19) . (2.27)The above equation shows that the local entropy fluctuation ˙ S i is proportional to the thermal momentum β i . Wewill return to relation (2.27) when discussing specific types of interatomic potentials in subsequent sections.Since the interatomic independence assumption results in an incorrect calculation of the interatomic heat fluxin (2.14a), one needs to model irreversible thermal transport, e.g., using the linear kinetic potential frameworkused by Venturini et al. (2014), as discussed in the following.
3. Dynamics and Quasistatics of independent GPPs
We proceed to analyze the dynamic behavior of the GPP equations (2.25) (under the interatomic indepen-dence assumption) and subsequently deduce the quasistatic behavior as a limit case. For instructive purposes(and because the (quasi-)harmonic assumption plays a frequent role in atomistic analysis), we apply the equa-tions to both harmonic and anharmonic potentials that describe atomic interactions within the crystal. We showthat, for a quasistatic change in the thermodynamic state of a crystal composed of the GPP atoms (i.e., drivingthe mean dynamical and thermal momenta p and β , respectively, to zero), the information of the evolution of Ω is lost (cf. (2.25)). Correspondingly, we may assume a specific nature of the thermodynamic process of interest(e.g., isothermal, isentropic, etc.) to determine the change in Ω of the GPP atoms during the quasistatic changein the thermodynamic state. Alternatively, the decay of correlation β ( t ) can be modeled empirically to obtainthe change in Ω , if the nature of the thermodynamical process is unknown. As a simplified case that admits analytical treatment, we consider a harmonic approximation of the interac-tion potential V ( q ) , writing V ( q ) = V + 12 N (cid:88) i =1 (cid:88) j ∈N ( i ) ( q i − q j ) T K ( q i − q j ) , (3.1)8here K ∈ R × is the local harmonic dynamical matrix, V is the equilibrium potential between the atomsapproximated as independent GPPs, and N ( i ) represents the neighbourhood of the i th atom. Equations (2.25)with (3.1) become d (cid:104) q (cid:105) i d t = − (cid:88) j ∈N ( i ) K (cid:10) q i − q j (cid:11) (3.2)and d Σ i d t = 2Ω i m − n i tr( K ) m Σ i , (3.3a) dΩ i d t = − n i tr( K ) β i , (3.3b)where n i tr( K ) is an effective force constant, which depends on the number of immediate neighbours of the i th atom, denoted by n i . Equations (3.3) show that the mean mechanical displacement (cid:104) q (cid:105) i of atom i is decoupledfrom its thermodynamic displacements Ω i and Σ i for a harmonic potential field between the atoms. The resultingdecoupled thermodynamic equations dΩ i d t = − n i tr( K ) β i , dΣ i d t = 2 β i m , and d β i d t = Ω i m − n i tr( K )Σ i , (3.4)exhibit the following independent eigenvectors ( φ , φ + , φ − ) and corresponding eigenvalues ( ω , ω + , ω − ) : φ = mn i tr( K )10 , φ ± = − mn i tr( K )1 − imω ± , ω = 0 , ω ± = ± i (cid:114) n i tr( K ) m , (3.5)so that a general homogeneous solution Ω i ( t )Σ i ( t ) β i ( t ) = a φ + a + φ + e iω + t + a − φ − e iω − t , (3.6)is composed of constant and oscillatory components. Coefficients ( a , a + , a − ) are determined by the initial ther-modynamic state of each atom. The constant component φ corresponds to β = 0 and Ω i = Σ i / (2 mn i tr( K )) .To interpret the three terms within this solution, let us formulate the excess internal energy, which for the har-monic approximation (3.1) becomes E = (cid:104)H(cid:105) = (cid:104) V ( q ) (cid:105) + N (cid:88) i =1 (cid:10) | p i | (cid:11) m = N (cid:88) i =1 (cid:18) Ω i m + n i tr( K )Σ i (cid:19) . (3.7)By insertion into (3.7), it becomes apparent that the constant component φ with Ω i = 2 mn i tr( K )Σ i hasequal average kinetic and potential energies. This equipartition of energy implies that φ corresponds to thethermodynamic equilibrium. Consequently, the components φ ± correspond to oscillations about the equilibriumstate φ with frequencies ω ± = 2 (cid:112) n i tr( K ) /m . Due to the decoupling of the thermodynamic equations (3.4)from the dynamic equation of motion (3.2), a harmonic GPP lattice exhibits no thermomechanical coupling andhence displays no heating or cooling of the lattice under external stress (and expansion due to local heating).Finally, by substituting the harmonic potential (3.1) into (2.27), we obtain the reversible fluctuations inentropy of atom i for a system of atoms in a harmonic field: d S i d t = 3 k B n i tr( K ) β i ˙ β i (cid:0) Ω i Σ i − β i (cid:1) Ω i Σ i − β i . (3.8) As an harmonic approximation of the interatomic potential renders the GPP crystal thermomechanicallydecoupled, as discussed above, we next study the effects of anharmonicity in the potential. As the simplestpossible extension of the harmonic potential, we now consider V ( q ) = V + 12 N (cid:88) i =1 (cid:88) j ∈N ( i ) ( q i − q j ) T K ( q i − q j ) + 16 N (cid:88) i =1 (cid:88) j ∈N ( i ) ( q i − q j ) T ζ (cid:2) ( q i − q j ) ⊗ ( q i − q j ) (cid:3) , (3.9)9here K ∈ R × is the local dynamical matrix, and ζ denotes a constant anharmonic third-order tensor. Withthis potential, the evolution equations (2.25) become d (cid:104) q (cid:105) i d t = − (cid:88) j ∈N ( i ) K (cid:10) q i − q j (cid:11) − (cid:88) j ∈N ( i ) ζ (cid:10) ( q i − q j ) ⊗ (cid:0) q i − q j (cid:1)(cid:11) , = − (cid:88) j ∈N ( i ) K (cid:10) q i − q j (cid:11) − (cid:88) j ∈N ( i ) ζ (cid:2) (Σ i + Σ j ) I + (cid:10) q i − q j (cid:11) ⊗ (cid:10) q i − q j (cid:11)(cid:3) (3.10)and d Σ i d t = 2Ω i m − m n i tr( K )Σ i + (cid:88) j ∈N ( i ) ζ lmn (cid:10) ( q i − q j ) l ( q i − q j ) m ( q i − (cid:104) q i (cid:105) ) n (cid:11) , = 2Ω i m − i m n i tr( K ) − (cid:88) j ∈N ( i ) ζ lmn (cid:16) δ ml (cid:10) q j (cid:11) n + δ nl (cid:10) q j (cid:11) m (cid:17) , (3.11a) dΩ i d t = − n i tr( K ) β i − (cid:88) j ∈N ( i ) ζ lmn (cid:10) ( q i − q j ) l ( q i − q j ) m ( p i − (cid:104) p i (cid:105) ) n (cid:11) , = − β i n i tr( K ) − (cid:88) j ∈N ( i ) ζ lmn (cid:16) δ ml (cid:10) q j (cid:11) n + δ nl (cid:10) q j (cid:11) m (cid:17) , (3.11b)where ζ lmn are the components of ζ , ( · ) l denotes the l th component of vector ( · ) , and δ represents Kronecker’sdelta (and we use Einstein’s summation convention, implying summation over l, m, n ). Note that the secondterm in (3.10) couples the mechanical perturbations with the thermodynamic perturbations of atom i and itsneighbors. Moreover, since equations (3.11) contain products of the thermodynamic variables Σ and β with themean mechanical displacements (cid:104) q (cid:105) , the anharmonic potential leads to thermomechanical coupling in the GPPevolution equations.Due to the apparent harmonic nature of most standard interatomic potentials, the time scale of the GPPequations (3.10) and (3.11), when being applied to common potentials, is comparable to the time scale ofatomic vibrations, since the system exhibits eigenfrequencies of (cid:112) n i tr( K ) /m for a pure harmonic potential(cf. (3.5)). Consequently, numerical time integration of the GPP equations incurs a similar computational cost asa standard molecular dynamics simulation of an identical system. Thus, even though mean motion and statisticalinformation have been separated, the interatomic independence assumption within the GPP framework preventssignificant temporal upscaling. Using the insight gained from the time evolution equations (3.3) and (3.11), we proceed to study the GPPequations within the quasistatic approximation. The latter yields a system of coupled nonlinear equations, whosesolution yields the thermodynamic state of the crystal composed of the GPP atoms. In the quasistatic approx-imation, the GPP equations (2.25) with mean mechanical momentum p i = 0 and thermal momentum β i = 0 reduce to the following steady-state equations: (cid:104) F i ( q ) (cid:105) = 0 , Ω i m + (cid:104) F i ( q ) · ( q i − q i ) (cid:105) , (3.12)which are to be solved for the equilibrium parameters ( q i , Σ i , Ω i ) for each atom i . As a result, the quasistaticassumption results in the loss of information about the process through which the system is brought to theequilibrium state, since the evolution of β i is unknown during the process. Moreover, equations (3.12) areinsufficient to solve for all three equilibrium phase-state parameters ( q i , Σ i , Ω i ) . Consequently, the quasistaticequations (3.12) can only be solved for each atom if a specific thermodynamic process is assumed and posedas an additional constraint. To elucidate the nature of a thermodynamic process, we assume that the ergodichypothesis holds for quasistatic processes, i.e., (cid:104) A ( z ) (cid:105) = 1 τ (cid:90) τ A ( z ) d t, (3.13)10here τ is a sufficiently large time interval, over which the evolution of the system is assumed quasistatic. In theergodic limit, the momentum variance becomes Ω i = mk B T i , (3.14)where T i is the local temperature of the i th atom. Since Ω i is proportional to the local temperature, the quasistaticequations (3.12) can be solved for the equilibrium parameters ( q i , Σ i ) for an isothermal process by keeping Ω i for all atoms constant. By contrast, isentropic equilibrium parameters ( q i , Σ i , Ω i ) can be obtained by keeping S i fixed for each atom. From (2.15) we know that S i = S ,i + 3 k B ln (cid:18) √ Ω i Σ i h (cid:19) = ˜ S + 3 k B S Σ ,i + 3 k B S Ω ,i = const ., (3.15)where ˜ S = S ,i − k B ln h . Upon using suitable dimensional constants of unit values, parameters S Ω ,i = ln Ω i and S Σ ,i = ln Σ i may be interpreted as dimensionless momentum variance and displacement-varianceentropies, respectively. In the following, it will be convenient to use S Ω ,i and S Σ ,i as the mean free parametersinstead of Ω i and Σ i . Analogously, isobaric conditions can be derived from the system-averaged Cauchy stresstensor (Admal and Tadmor, 2010) σ = − V (cid:88) i (cid:28) p i ⊗ p i m + F i ( q ) ⊗ ( q i − q i ) (cid:29) , (3.16)where V is the volume of the system. The average hydrostatic pressure p of the system is p = − tr( σ )3 = 1 V (cid:88) i (cid:18) k B T i + (cid:104) F i ( q ) · ( q − q i ) (cid:105) (cid:19) . (3.17)Hence, setting p = const . in (3.17) is the isobaric constraint, subject to which equations (3.12) become (cid:104) F i ( q ) (cid:105) = 0 , (cid:18) k B T i − pVN (cid:19) + (cid:104) F i ( q ) · ( q − q i ) (cid:105) . (3.18)Solving these equations, in which Ω i = m (cid:16) k B T i − pVN (cid:17) serves as the momentum variance corrected for pres-sure p , yields the equilibrium parameters ( q i , S Σ i , S Ω i ) for a given externally applied pressure p . Note that fora system of non-interacting atoms, the quasistatic equations subjected to an external pressure p reduce to theideal gas equation pV = N k B T . Consequently, equation (3.18) shows that the quasistatic approximation yieldsthe thermal equation of state of the system accounting for the interatomic potential, which enables the ther-momechanical coupling within the crystal. The three constraints on Ω i for isentropic, isothermal, and isobaricprocesses – based on the above discussion – are summarized in Table 1.For numerical validation, we compute the thermal expansion coefficient and elastic constants (bulk modu-lus, shear modulus, and uniaxial modulus) of a single-crystal of pure copper (Cu). For computing phase-spaceaverages we use third-order multivariable Gaussian quadrature (Stroud, 1971; Kulkarni, 2007). Figure 3 illus-trates solutions of the isothermal equilibrium relations at various levels of fixed uniform temperature T i = T for all atoms within a Cu crystal, modeled using the EAM potentials of Dai et al. (2006) and Johnson (1988).As the temperature T increases, the local separation of atoms increases, thus increasing the volume of the crys-tal. Furthermore, as the spacing between atoms increases, the displacement-variance entropy S Σ also increases(cf. Figure 3 ( b ) ). Further evident from that graph, S Σ is not uniform within the crystal, owing to varying numbers n i of interacting neighboring atoms on the corners, edges, and faces as well as within the bulk of the crystal.Process Isentropic Isothermal IsobaricConstraint Ω i Σ i = const . Ω i = k B mT i Ω i = k B mT i − pV mN Table 1: Summary of the thermodynamic constraints on Ω i for the different assumptions about the thermodynamic processunder which the system is brought to equilibrium. a ) ( b ) ( c ) Figure 3: Thermal expansion of cubes of pure single-crystalline copper, composed of and atomic unit cells, andmodeled using the Extended Finnis-Sinclair potential (Dai et al., 2006) and the exponentially decaying potential of John-son (1988). ( a ) Comparison of computed volumetric changes with the experimental data obtained from Nix and MacNair(1941). ( b ) Variation of the displacement-variance entropy S Σ with temperature. As the temperature increases, the vi-brational kinetic energy increases, resulting in an increase of Σ at thermal equilibrium due to the equipartition of energy. ( b, c ) Due to different numbers of interacting atomic neighbours, atoms on the corners, edges, faces, and in the bulk exhibitdifferent values of S Σ . Note that the value of S Σ depends on the interatomic potential (shown results are for a -unit-cellcube modeled using Johnson’s potential). ( c ) shows the spatial variation of S Σ for the same cube at 1000 K. Figure 4 shows computed elastic constants (bulk modulus, shear modulus, and uniaxial modulus) using thelocal equilibrium relations (3.12) for different fixed, uniform temperatures applied to a pure Cu single-crystalmodeled using the potential of Dai et al. (2006). This already served as a benchmark in Kulkarni et al. (2008),Venturini (2011), and Tembhekar (2018), who all used the max-ent formulation. Since the GPP quasistaticequations (3.12) are identical to those of the max-ent formulation, we obtain identical results. Elastic moduliare computed from the local equilibrium relations with the isentropic constraint, since the elastic constants aremeasured in an adiabatic setting (Overton Jr and Gaffney, 1955). We note that the values of the reported elasticconstants are offset at all temperatures, since Dai et al.’s EAM potential was calibrated using the elastic constantsat finite temperature values, which are treated as those at 0 K in the present formulation.
The solution ( q , S Σ , S Ω ) of the local equilibrium relations (3.12) may be re-interpreted as a minimizer of theHelmholtz free energy F (note that ( q , S Σ , S Ω ) denotes the whole set of parameters of all N atoms constitutingthe system). The Helmholtz free energy F is defined as F ( q , S Σ , S Ω ) = inf S (cid:40) E ( q , S Σ , S ) − (cid:88) i Ω i S i k B m (cid:41) , (3.19)with the internal energy of the system being E ( q , S Σ , S ) = (cid:104)H(cid:105) = (cid:88) i (cid:18) i m + (cid:104) V i ( q ) (cid:105) (cid:19) . (3.20)The definition (3.19) implies the local thermodynamic equilibrium definition Ω i k B m = ∂E∂S i , (3.21)which can be verified using (3.15) and (3.20). In addition, minimization of F with respect to the parameter sets q and S Σ , subject to any of the thermodynamic constraints in Table 1 for updating S Ω , yields equations (3.12),i.e., − ∂ F ∂ q i = 0 = ⇒ (cid:104) F i ( q ) (cid:105) = 0 , (3.22a)and − ∂ F ∂S Σ ,i = 0 = ⇒ i m + (cid:104) F i ( q ) · ( q i − q i ) (cid:105) = 0 . (3.22b)12 igure 4: Elastic constants of a Cu single-crystal (black: uniaxial, red: bulk, blue: shear modulus) across a wide rangeof temperatures evaluated following the procedure outlined by Venturini (2011) and Amelang et al. (2015), using the EAMpotential of Dai et al. (2006), compared against the experimental data from Chang and Himmel (1966) and Overton Jr andGaffney (1955). A detailed derivation of (3.22b) is provided in Appendix B. We point out that equations (3.22) are identical tothe max-ent framework, as derived previously (Kulkarni et al., 2008; Ariza et al., 2012; Venturini et al., 2014;Ponga et al., 2015; Tembhekar, 2018). However, the max-ent based framework is based on a different motivationand does not rely on the physical insight gained in Sections 2.2, 3.1.1, and 3.1.2 about dynamics of atoms atlong and short time intervals. Moreover, the GPP framework highlights the information loss as a result of thequasistatic approximation (otherwise hidden in the max-ent framework) which also enables the modeling of ther-momechanical deformation of crystals under various thermodynamical conditions (e.g., isentropic, isobaric, orisothermal processes). Furthermore, recall that in Section 2.2 we showed that the interatomic independence as-sumption precludes the framework from capturing any changes in local temperature due to unequal temperaturesof neighbouring atoms. Such an assumption is also made in the max-ent framework. In reality, a non-uniformtemperature distribution may easily arise in a non-uniformly deformed crystal lattice. For example, when a crys-tal is deformed by applied mechanical stresses, the local temperature rises as atoms are compressed to smallerinteratomic distances. Therefore, both max-ent and the present GPP formulation require explicit thermal trans-port models in order to capture the latter. To model the transport of heat in a non-uniformly deformed lattice,we here adopt the linear Onsager kinetics model introduced by Venturini et al. (2014), based on the quadraticdissipation potential – as discussed in the following.
The total rate of change of entropy of an atom can be decomposed into a reversible and an irreversiblechange, i.e., d S i d t = q i, rev T i + q i, irrev T i , (3.23)where q i, rev is the reversible heat addition and q i, irrev is the irreversible change signifying the net influx ofheat into the atom due to a non-uniform temperature distribution. In a dynamic system (see Sections 3.1.1and 3.1.2), the reversible change is composed of fluctuations about the equilibrium state due to the local harmonicnature of the interatomic potential, and it is proportional to the thermal momentum β . Within the quasistaticapproximation, the information about the evolution of β ( t ) as the system relaxes towards the equilibrium is lost13ince β → is imposed, thus rendering the reversible changes in entropy unknown. Therefore, such a reversiblechange is imposed implicitly by the thermodynamic constraints (see Table 1).In general, for an isolated system of atoms the reversible heat addition vanishes ( q rev = 0 ), and the systemrelaxes to equilibrium adiabatically, which we here term free relaxation . Note that this free relaxation is notisentropic, since the temperature can be non-uniform as a result of an imposed non-uniform deformation of thecrystal, resulting in irreversible thermal transport that increases the entropy. We model such irreversible change q i, irrev by adopting the kinetic formulation introduced by Venturini et al. (2014): q i, irrev = (cid:88) j (cid:54) = i R ij , R ij = ∂ Ψ ∂P ij , (3.24)where R ij is the local, pairwise heat flux, driven by a local pairwise discrete temperature gradient P ij through thekinetic potential Ψ . Using the dissipation inequality, Venturini et al. (2014) formulated the discrete temperaturegradient as P ij = 1 T i − T j . (3.25)Within the linear assumption, the kinetic potential Ψ is modeled as, Ψ = 12 N (cid:88) i =1 (cid:88) j (cid:54) = i A ij T ij P ij with T ij = 12 ( T i + T j ) , (3.26)where A ij denotes a pairwise heat transport coefficient (which is treated as an empirical constant in this work),and T ij represents the pairwise average temperature. Equations (3.26) and (3.24) yield the entropy rate kineticequation for a freely relaxing system of atoms as d S i d t = 1 T i (cid:88) j (cid:54) = i A ij T ij P ij = 1 T i (cid:88) j (cid:54) = i A ij T ij P ij , (3.27)which yields the following discrete-to-continuum relation between A ij and the thermal conductivity tensor κ i atthe location of atom i (Venturini et al., 2014): κ i = 12 V i N (cid:88) j =1 A ij (cid:0) q i − q j (cid:1) ⊗ (cid:0) q i − q j (cid:1) , (3.28)where V i is the volume of the atomic unit cell in the crystal. Venturini et al. (2014) validated the thermal transportmodel based on linear Onsager kinetics by demonstrating the capturing of size effects of the macroscopic thermalconductivity of Silicon nanowires. Experimentally measured values of κ i for a given arrangement of atoms canbe used to obtain A ij from (3.28). The obtained values of A ij may be interpreted to capture the interatomic heatcurrent and regarded as the intrinsic property of the material.Relations (3.27) (for i = 1 , . . . , N ) define a system of nonlinear equations governing the irreversible changeof entropy of each atom due to thermal transport with its neighboring atoms. In general, thermal transportconsists of nonlinear interactions of phonons across all wavenumbers (Ziman, 2001), which are globally definedin the system. Comparing with (2.14a), equation (3.27) can be regarded as a phenomenological approximationof the interatomic correlations, with A ij as the empirical coefficient. Consequently, we here deal with localdissipative entropy transport, in which atoms correlate only with their interacting neighbors via the transportequation (3.27) (see Figure 5 for an illustration).Equations (3.12) combined with (3.27) complete the nonequilibrium thermomechanical atomistic model.Every atom is assumed to be in thermal equilibrium at some local temperature T i , and mechanical and thermalinteractions of the atoms with different spacing and different temperatures are modeled using the interatomicpotential V ( q ) , the local equation of state (e.g. (3.18)), and the entropy kinetic equation (3.27). While thenature of the thermodynamic process via which the system relaxes depends on the macroscopic and microscopicboundary conditions, the process is generally assumed quasistatic with respect to the fine-scale vibrations ofeach atom. However, the thermal transport equation (3.27) introduces a time scale to the problem governing14 igure 5: Schematic illustration of the local entropy dissipation via equation (3.27) : atom i interacts with its interatomicneighbors, whereby differences in local temperature are responsible for heat flux. the irreversible evolution of entropy. If we consider a two-atom toy model, consisting of only two atoms withtemperatures T i and T j (Figure 5), equations (3.27) reduce to dd t (cid:18) S i S j (cid:19) = (cid:18) A T ij P ij /T i A T ji P ji /T j (cid:19) , (3.29)where we have assumed equal coefficients A ij = A for both atoms. Let us further assume the thermomechanicalrelaxation of ( q , S Σ ) takes place through an isentropic process following the quasistatic equations (3.12). Whensubstituting equation (3.15), equation (3.29) becomes dd t (cid:18) T i T j (cid:19) = 2 A k B (cid:18) T ij P ij T ji P ji (cid:19) . (3.30)The stationary state of equation (3.30) is a state with uniform temperature, T i = T j . When assuming a leading-order perturbation about the stationary state, equation (3.30) yields dd t (cid:18) T i T j (cid:19) ≈ A k B (cid:18) T j − T i T i − T j (cid:19) , (3.31)which shows that, when using a simple first-order explicit-Euler finite-difference scheme for time integration,the time step ∆ t is restricted by ∆ t ≤ k B A (3.32)for numerical stability. Venturini et al. (2014) used A = 0 . nW/K for the bulk region of a silicon nanowire,which yields a maximum allowed time step of . ps. With higher-order explicit schemes, a larger maximumtime step may be obtained, but the restriction remains at a few picoseconds. Such restriction arises because thebulk thermal conductivity κ is used to fit the atomistic parameter A ij , which reduces the length scale as wellas the time scale. Such a restriction also highlights the coupling of length and time scales (Evans and Morriss,2007). For larger temperature differences and larger systems, the nonlinear equation (3.30) yields the followingstability limit on time step ∆ t (see Appendix C for the derivation) ∆ t (cid:18) A k B (cid:19) (cid:88) j T ij T i T j ≤ . (3.33)In the following sections, we will apply the linear kinetic formulation (3.27) to the QC framework to under-stand the deformation of a coarse-grained Cu crystal. As a first step, we fit the kinetic coefficient A for Cu,using the thermal conductivity measurements for Cu nanowires obtained by Mehta et al. (2015). To this end, weconsider a Cu nanowire of square cross-section of of size 145 × ×
43 ˚A , with the central axis of the nanowireoriented along the x -direction and all edges aligned with the (cid:104) (cid:105) -directions (see Figure 6 a ). Atoms in the re-gion x < x l are thermostatted at a temperature of 360 K, while the atoms in the region x > x r are thermostattedat a temperature of 240 K. As the system evolves according to (3.27) (using full atomistic resolution everywhere15 a ) ( b ) Figure 6: Thermal conduction in a Cu nanowire with square cross-section and edges aligned with the (cid:104) (cid:105) -directions. ( a ) Schematic showing the cross-sectional planes used for computing the macroscopic thermal flux. ( b ) Spatial variationof temperature (black solid line) and macroscopic heat flux (red dashed line) along the axis of the nanowire. in the nanowire), a uniform macroscopic heat flux is established between x l < x < x r (see Figure 6 b ). Dividingthe nanowire into atomic planes marked by x i , the macroscopic flux across the plane at x = x i is given by J x,i = i (cid:88) j =1 T j d S j d t , (3.34)where T j and dS j /dt are temperature and entropy generation at plane x j . From (3.34) the approximate thermalconductivity κ is obtained as κ = J x,i S A ∆ x ∆ T , (3.35)where S A is the cross-sectional area, ∆ x is the length across which the temperature difference ∆ T is maintained.To find A , we perform the above simulation with an initial guess of A , then compute the macroscopic flux J x,i and from it the thermal conductivity κ via (3.35). Based on the computed κ , we tune the value of A until the desired value of the effective thermal conductivity κ is obtained. In metals, phonon and electronscattering events contribute to the thermal and electrical resistivities (Ziman, 2001). Due to surface scatteringof electrons in nanowires, the thermal conductivity drops significantly. We therefore use a reduced value ofthe thermal conductivity, κ = 100 W / ( m · K ) , instead of the bulk thermal conductivity of W / ( m · K ) , asexperimentally determined by Mehta et al. (2015), yielding an approximate kinetic constant A ≈ . eV / ( ns · K ) = 2 . nW/K.Before applying the QC coarse-graining, let us summarize the proposed thermomechanical transport modelfor simulating quasistatic deformation of crystals composed of GPP atoms (see Algorithm 1 for details):• Step 1 : Given the equilibrium parameters (cid:16) q ( n ) i , S ( n )Σ i , S ( n )Ω i , S ( n ) i (cid:17) from the (previous) n th load step, anexternal stress/strain is applied to the system at load step n + 1 ,• Step 2 : Quasistatic relaxation, solving (3.12) subject to one of the constraints in Table (1) to obtain theintermediate state (cid:16) q ( ∗ ) i , S ( ∗ )Σ i , S ( ∗ )Ω i , S ( ∗ ) i (cid:17) .• Step 3 : Staggered time stepping that alternates between irreversible updates of the total entropy and qua-sistatic reversible relaxation steps of all variables, until convergence is achieved. Specifically, the totalentropy is updated from S ( ∗ ) i over a time interval δt , using equation (3.27) and explicit forward-Eulerupdates. The interval δt ( n ) depends on the external stress/strain rate applied. By definition, the ther-momechanical relaxation is assumed quasistatic, hence only slow rates with respect to atomistic vibra-tions can be modeled. However, the irreversible transport imposes a time-scale restriction via (3.33).Consequently, time integration via suitable time steps ∆ t k must be continued for K steps, such that (cid:80) Kk =1 ∆ t k = δt . After each update of the entropy from S ( ∗ ) ,ki → S ( ∗ ) ,k +1 i , quasistatic relaxation of16tate (cid:16) q ( ∗ ) ,ki , S ( ∗ ) ,k Σ i , S ( ∗ ) ,k Ω i , S ( ∗ ) ,ki (cid:17) → (cid:16) q ( ∗ ) ,k +1 i , S ( ∗ ) ,k +1Σ i , S ( ∗ ) ,k +1Ω i , S ( ∗ ) ,k +1 i (cid:17) follows by solving (3.12)with a constraint from Table 1. In a full quasistatic setting, the transport equation is driven towards a steadystate with ˙ S ( ∗ ) ,Ki = 0 , which defines the convergence criterion and hence determines the total number oftime steps.• Step 4 : Assignment of the final state as (cid:16) q ( n +1) i , S ( n +1)Σ i , S ( n +1)Ω i , S ( n +1) i (cid:17) = (cid:16) q ( ∗ ) ,Ki , S ( ∗ ) ,K Σ i , S ( ∗ ) ,K Ω i , S ( ∗ ) ,Ki (cid:17) ,followed by Step 1 till the final loadstep. Algorithm 1:
Single load step from n to n + 1 of the quasistatic thermomechanical transport modelfor irreversible deformation of crystals composed of GPP atoms. For a fully quasistatic transportsimulation δt ( n ) → ∞ . Result: (cid:16) q ( n +1) i , S ( n +1)Σ i , S ( n +1)Ω i , S ( n +1) i (cid:17) Input: (cid:16) q ( n ) i , S ( n )Σ i , S ( n )Ω i , S ( n ) i (cid:17) , δt ( n ) , tol; k ← ;quasistatic reversible relaxation of (cid:16) q ( n ) i , S ( n )Σ i , S ( n )Ω i , S ( n ) i (cid:17) to (cid:16) q ( ∗ ) ,ki , S ( ∗ ) ,k Σ i , S ( ∗ ) ,k Ω i , S ( ∗ ) ,k (cid:17) by solving(3.12) with a constraint from Table 1; t ← ;compute ˙ S ( ∗ ) ,ki ; ˙ S i ← ˙ S ( ∗ ) ,ki ; while t < δt ( n ) and (cid:113) N (cid:80) i ˙ S i > tol do compute ∆ t k satisfying the constraint (3.33);irreversible update of S ( ∗ ) ,ki to S ( ∗ ) ,k +1 i using (3.27) and a forward-Euler finite-difference scheme; S ( ∗ ) ,ki ← S ( ∗ ) ,k +1 i ;quasistatic reversible relaxation of (cid:16) q ( ∗ ) ,ki , S ( ∗ ) ,k Σ i , S ( ∗ ) ,k Ω i , S ( ∗ ) ,k (cid:17) to (cid:16) q ( ∗ ) ,k +1 i , S ( ∗ ) ,k +1Σ i , S ( ∗ ) ,k +1Ω i , S ( ∗ ) ,k +1 (cid:17) by solving (3.12) with a constraint from Table 1; k ← k + 1 ; ˙ S i ← ˙ S ( ∗ ) ,ki ; t ← t + ∆ t k ; end (cid:16) q ( n +1) i , S ( n +1)Σ i , S ( n +1)Ω i , S ( n +1) i (cid:17) ← (cid:16) q ( ∗ ) ,ki , S ( ∗ ) ,k Σ i , S ( ∗ ) ,k Ω i , S ( ∗ ) ,k (cid:17)
4. Finite-temperature updated-Lagrangian quasicontinuum framework based on GPP atoms
Having established the GPP framework, we proceed to discuss the application of the thermomechanical andcoupled thermal transport model (equations (3.12) and (3.27), respectively) to an updated-Lagrangian QC formu-lation for coarse-grained simulations. Previous zero- and finite-temperature QC implementations have usuallybeen based on a total-Lagrangian setting (Ariza et al., 2012; Tadmor et al., 2013; Knap and Ortiz, 2001; Ame-lang et al., 2015; Ponga et al., 2015; Kulkarni et al., 2008), in which interpolations are defined and hence atomicneighborhoods computed in an initial reference configuration. Unfortunately, such total-Lagrangian calculationsincur large computational costs and render especially nonlocal QC formulations impractical (Tembhekar et al.,2017), since atomistic neighborhoods change significantly during the course of a simulation, so that the initialmesh used for all QC interpolations increasingly loses its meaning and atoms that form local neighborhoods inthe current configuration may have been considerably farther apart in the reference configuration. Therefore,we here introduce an updated-Lagrangian QC framework to enable efficient simulations involving severe de-formation and atomic rearrangement. Moreover, we adopt the fully-nonlocal energy-based QC formulation ofAmelang et al. (2015), since an energy-based summation rule allows for a consistent definition of the coarse-grained internal energy and all thermodynamic potentials of the system.17 a ) ( b ) Figure 7: Illustration of the the quasicontinuum (QC) framework based on GPP quasistatics (eqs. (3.12) ) combinedwith the linear Onsager kinetics for thermal transport (eq. (3.27) ). ( a ) The thermomechanical transport parameters ( q k , S Σ ,k , S Ω ,k , ˙ S k ) are the degrees of freedom of repatom k . Repatoms are shown as red circles, sampling atoms assmall white circles. All thermodynamic potentials and hence the quasistatic forces and thermal fluxes are computed froma weighted average over a set of sampling atoms (e.g., forces and fluxes for repatom k are governed, among others, bysampling atom α and its atomic neighbors j ). The fully-nonlocal formulation bridges seamlessly and adaptively from fullatomistics to coarse-grained regions. ( b ) Computation of sampling atom weights for the updated Lagrangian implementa-tion using tetrahedral solid angles.
Within the QC approximation, we replace the full atomic ensemble of N GPP atoms (as described in Sec-tion 2.2) by a total of N h (cid:28) N GPP representative atoms ( repatoms for short), each having the thermomechan-ical transport parameters ( q k , S Σ ,k , S Ω ,k , ˙ S k ) as their degrees of freedom (see Figure 7). The position of eachand every atom in the coarse-grained crystal lattice is obtained by interpolation. For an atom at location q hi inthe reference configuration, the thermomechanical transport parameters in the current configuration are obtainedby interpolation from the respective parameters of the repatoms: q hi S h Σ ,i S h Ω ,i ˙ S hi = N h (cid:88) k =1 q k S Σ ,k S Ω ,k ˙ S k N k ( q i ) , (4.1)where the subscript h denotes that the parameters are interpolated from the N h GPP repatoms, and N k ( q hi ) is the shape function/interpolant of repatom q k evaluated at q hi . In the following we use linear interpolation(i.e., constant-strain tetrahedra), while the method is sufficiently general to extend to other types of interpolants.Based on the interpolated parameters from (4.1), the free energy F ( q, S Σ , S Ω ) of the crystal is replaced by theapproximate free energy F h of the QC crystal with F h ( q h , S h Σ , S h Ω ) = N (cid:88) i =1 (cid:18) Ω hi m i − Ω hi S hi k B m i (cid:19) + (cid:104) V ( q ) (cid:105) = N (cid:88) i =1 (cid:18) Ω hi m i − Ω hi S hi k B m i + (cid:104) V i ( q ) (cid:105) (cid:19) , (4.2)where we assumed that the decomposition V ( q ) = (cid:80) i V i ( q ) of the interatomic potential holds. Furthermore,we allow the masses of all atoms to be different, denoting by m i the mass of atom i . Equation (4.2) definesthe free energy of the system accounting for all the atoms N with their thermomechanical parameters evaluatedusing the N h repatoms in a slave-master fashion.To reduce the computational cost stemming from the summation over all N atoms in (4.2), sampling rulesare introduced, which approximate the full sum by a weighted sum over N s (cid:28) N carefully selected samplingatoms (Eidel and Stukowski, 2009; Iyer and Gavini, 2011; Amelang et al., 2015; Tembhekar et al., 2017).Specifically, we adopt the so-called optimal summation rule of Amelang et al. (2015) to sample the free energyat N s sampling atoms. Consequently, the approximate free energy F h is further approximated as F h ( q h , S h Σ , S h Ω ) ≈ N s (cid:88) α =1 w α (cid:18) Ω hα m α − Ω α S α k B m α + (cid:104) V α ( q ) (cid:105) (cid:19) , (4.3)18here w α is the sampling weight of the α th sampling atom.We use the first order summation rule of Amelang et al. (2015), in which all nodes and the centroid of eachsimplex (tetrahedron in 3D) are assigned as sampling atoms. Amelang et al. (2015) computed the sampling atomweights using the geometrical division of the simplices by planes at a distance r from the nodes (Figure 7 ( b ) )and adding the corresponding nodal volume to the respective sampling atom weight, while the rest of the simplexvolume was assigned to the Cauchy-Born-type sampling atom at the centroid. Here, we point out that a simplerweight calculation is possible by considering the spherical triangle generated by the intersection of simplex e with a ball of radius r centered at one of the nodes. Considering the arcs of the spherical triangle subtend angles α , β , and γ at the opposite points, the area of the triangle is given by ( α + β + γ − π ) r . Hence, the approximatevolume of the enclosed region is v eα ≈ r α + β + γ − π ) , (4.4)and w α = ρ e (cid:80) e v eα are the sampling atom weights at nodes, where ρ e is the density of simplex e (expressedas the number of atoms per unit volume). For the centroid sampling atoms, the remaining volume times ρ e isassigned as its sampling weight w α . Since the deformation is affine within each element e , sampling atomweights in coarse-grained regions change negligibly in a typical simulation and are therefore kept constantthroughout our simulations. In the following we will also need a separate set of repatom weights (cid:98) w k , whichwe calculate by lumping the sampling atom weights w k to the repatoms: each repatom receives the weight ofitself (each repatom is a sampling atom) plus one quarter of the Cauchy-Born-type centroidal sampling atomswithin all adjacent elements e : (cid:98) w k = w k + (cid:88) e w e . (4.5)Given the sampling atom weights w α , minimization of the approximate free energy F h ( q h , S h Σ , S h Ω ) givenin equation (4.3) with respect to degrees of freedom ( q k , S Σ ,k ) of the k th repatom yields the local mechanicalequilibrium conditions − ∂ F h ∂ q k ≈ − N s (cid:88) α =1 w α ∂ (cid:104) V α ( q ) (cid:105) ∂ q k = − N s (cid:88) α =1 w α (cid:28) ∂V α ( q ) ∂ q k (cid:29) = 0 (4.6)and the corresponding thermal equilibrium conditions − ∂ F h ∂S Σ ,k ≈ N s (cid:88) α =1 w α (cid:18) α m α ∂S Σ ,α ∂S Σ ,k − ∂ (cid:104) V α ( q ) (cid:105) ∂S Σ ,k (cid:19) = N s (cid:88) α =1 w α (cid:18) α m α ∂S Σ ,α ∂S Σ ,k − (cid:28) ∂V α ( q ) ∂ q · ∂ q ∂S Σ ,k (cid:29)(cid:19) = 0 . (4.7)Substituting the interpolation from (4.1) into (4.6) and (4.7) yields (for repatoms k = 1 , . . . , N h ) − N s (cid:88) α =1 w α (cid:88) j ∈N ( α ) (cid:28) ∂V α ( q ) ∂ q j (cid:29) N k ( q j ) + (cid:28) ∂V α ( q ) ∂ q α (cid:29) N k ( q α ) = 0 (4.8a)and N s (cid:88) α =1 w α α m α N k ( q α ) − (cid:88) j ∈N ( α ) (cid:28) ∂V α ( q ) ∂ q j · (cid:0) q j − q j (cid:1)(cid:29) N k ( q j ) + (cid:28) ∂V α ( q ) ∂ q α · ( q α − q α ) (cid:29) N k ( q α ) = 0 . (4.8b)Note that these are similar to the equilibrium conditions derived by Tembhekar (2018) following Kulkarni et al.’smax-ent formulation, although the max-ent formulation bypasses the thermodynamic relevance of the parame-ters in the dynamic setting. As noted previously in Section 3.2, the local thermal equilibrium equation (4.8b)corresponds to the local equation of state of the system, here providing the equation of state of the coarse-grainedquasicontinuum. Solving equations (4.6) and (4.7) subject to one of the constraints from Table 1 depending uponthe assumption of the thermodynamic process yields the variables ( q k , S Σ ,k , S Ω ,k ) for all repatoms, thus yieldingthe thermodynamically reversible solution for the deformation of the system.To introduce the irreversible thermal conduction of the linear Onsager kinetic model (cf. equation (3.27)),we solve for the entropy rates ˙ S k of all repatoms and evolve the entropy in time for each repatom. To this end,19e notice that the first term in equation (4.8b) represents a thermal force due to the thermal kinetic energy of thesystem. Since that thermal force in our energy-based setting follows a force-based summation rule of Knap andOrtiz (2001), the the entropy rate calculation of the repatom k simplifies to (cid:98) w k d S k d t = N s (cid:88) α =1 w α d S α d t N k ( q α ) = N s (cid:88) α =1 w α T α (cid:88) j ∈N ( α ) A αj T αj P αj N k ( q α ) , (4.9)where (cid:98) w k is the repatom weight. Note that the kinetic potential in equation (3.26) can, in principle, also becoarse-grained analogously to the free energy in equation (4.3). However, the resulting calculation of the entropyrates is computationally costly since it involves summation over all repatoms for each sampling atom calculation,which is why this approach is not pursued here. Equations (4.8a) and (4.8b) combined with the thermodynamicconstraints in Table 1 and the coarse-grained thermal transport model in equation (4.9) yield the solution of ageneric thermomechanical deformation problem subject to a non-uniform temperature distribution and loadingand boundary conditions, as illustrated in Figure 1. We implement the thermomechanical local equilibrium relations (4.8a) and (4.8b) combined with a ther-modynamic constraint from Table 1 and the coarse-grained thermal transport equation (4.9) in an updated-Lagrangian QC setting. The latter is chosen since atoms in regions undergoing large deformation tend to havesignificant neighborhood changes, for which the initial reference configuration loses its meaning in the fully-nonlocal QC formulation, as illustrated by Amelang (2016) and Tembhekar et al. (2017). Consequently, trackingthe interatomic potential neighborhoods in the undeformed configuration incurs high computational costs. Al-ternatively, one could strictly separate between atomistic and coarse-grained regions (as in the local-nonlocalQC method of Tadmor et al. (1996a)), yet even this approach suffers from severe mesh distortion in the coarse-grained regions in case of large deformation, and it furthermore does not easily allow for the automatic trackingof, e.g., lattice defects with full resolution (Tembhekar et al., 2017). It also requires a-priori knowledge aboutwhere full resolution will be required during a simulation. As a remedy, we here deform the mesh with themoving repatoms and we take the deformed configuration from the previous load step as the reference configu-ration for each new load step, thus discarding the initial configuration and continuously updating the referenceconfiguration.For every element e , we store the three initial edge vectors (i.e., three node-to-node vectors forming a right-handed system) in a matrix S e , and the three Bravais lattice vectors indicating the initial atomic arrangementwithin the element in a matrix A e . As the system is relaxed quasistatically under applied loads, all repatomsmove to the deformed configuration (e.g., from load step n = i to n = i + 1 ), thus deforming the edge vectorsof element e from S ei to S ei +1 (and likewise the Bravais basis from A ei to A ei +1 ). The incremental deformationgradient of element e , from step i to i + 1 , can hence be computed from the kinematic relation F ei → i +1 = S ei +1 ( S ei ) − , (4.10) Figure 8: Illustration of the updated-Lagrangian QC implementation at different external load/strain steps denoted by n .The local Bravais basis A n of the highlighted element is shown in blue, with the edge vectors S n . Repatoms are shownas red circles, sampling atoms as small white circles. Deformation gradient F i → i +1 deforms the edge vectors S i to S i +1 and the local Bravais basis A i to A i +1 . A ei +1 = F ei → i +1 A ei . (4.11)Consequently, the integer matrix N , which contains the numbers of lattice vector hops along the element edges,evaluated as N ei = S ei ( A e,i ) − = const ., (4.12)remains constant throughout deformation of a given element e . Moreover, each element edge has a constantnumber vector, denoted by the rows of N ei (see Figure 8). That is, in the updated-Lagrangian setting, thenumber matrix N ei remains constant during deformation. Such conservation of lattice vector hops along theelement edges/faces is particularly useful for adaptive remeshing scenarios, where existing elements may needto be removed and new elements need to be reconnected, with or without changes to the number of lattice sitesused for re-connections. The conservation of lattice vector hops can then be used for computing the Bravaislattice vectors local to new elements. The Bravais lattice vectors are used for calculating the neighborhoods ofthe nodal and centroid sampling atoms belonging to the large elements in the fully nonlocal QC formulation. Thelocal lattice is generated within a threshold radius distance from the sampling atom using those lattice vectors.We use the adaptive neighborhoods calculation strategy of Amelang (2016), which requires larger threshold radiicompared to the interatomic potential cut-off chosen as r th = r cut + r buff , (4.13)where r buff is a buffer distance used for triggering re-computations of neighborhoods, and r cut is the interatomicpotential cut-off. If the maximum relative displacement among the neighbors with respect to a sampling atomexceeds r buff , then neighborhoods of the sampling atom are re-computed (see Amelang (2016) for details).Within the region with atomistic resolution, only nodal sampling atoms have finite weights (close to unity)and hence only their neighborhoods are computed. For such neighborhood calculations Bravais lattice vectors arenot required. Instead, the unique nodes of all elements within the threshold radius of the sampling atom are addedas the neighbors. Consequently, even severely deforming meshes do not require element reconnection/remeshingas long as the deformation stays restricted within the atomistic region, since only nodes of the elements arerequired. Hence, we use meshes with large atomistic regions in the benchmark cases presented below, to restrictthe analysis towards thermodynamics of the deformations. Such simulations do not require adaptive remeshing,the analysis of which is left for future studies. As a benchmark example, we use the updated-Lagrangian QC method discussed above to analyze the effectsof temperature on dislocations and specifically on edge dislocation nucleation under an applied shear stress.We present the analysis for both cases of isothermal and adiabatic conditions (the latter combined with the ( a ) ( b ) Figure 9: Initial conditions of the dislocation dipole setup. ( a ) Initial condition after displacing the atoms according to theisotropic linear elastic displacement field solution. ( b ) Isothermally relaxed state consisting of a single atom void createdby annihilation of the dislocations.Shown are 3D views of the full simulation domain with a magnified view of the fully-resolved central region and a top view. Atoms are colored by the centrosymmetry parameter in arbitrary units and shownbetween threshold values of to . a ) ( b )( c ) ( d ) Figure 10: Comparison of the isothermal nucleation of dislocation dipoles from a single-atom void as obtained fromfully atomistic (left) and QC (right) simulations at varying temperatures in a Cu single-crystal modeled using the EFSpotential (Dai et al., 2006). Shown is the final sheared state of ( a ) snapshot of the atomistic simulation and ( b ) that of theQC simulation. Atoms are colored by the centrosymmetry parameter in arbitrary units,. While the atoms in region A arekept fixed, atoms in region B are allowed to relax. (c,d) The shear stress τ vs. the engineering strain γ is plotted for ( c ) theatomistic simulation and ( d ) the QC simulation. The shear stress is evaluated as the net force in the [110] -direction on theatoms in region A per cross-sectional area in the (111) plane. Faces (112) are periodic. irreversible entropy transport based on linear Onsager kinetics). For both cases, we generate a pair of dislocations(i.e., a dislocation dipole) using the isotropic linear elastic displacement field solutions of edge dislocations withopposite Burgers’ vectors (Nabarro, 1967), superposed linearly in a × × . nm slab of pure single-crystalline Cu, consisting of 125,632 lattice sites and edges oriented along the slip crystallographic directions.Figure 9 ( a ) shows the initial condition generated for the edge dislocation dipole with Burgers’ vectors b = ± a [110] , separated by a distance of ˚A. Displacements are restricted to the (112) plane, while the simulationdomain is set up in 3D with periodic boundary conditions on opposite out-of-plane faces. After initial relaxation,the dislocations annihilate each other due their interacting long-range elastic field, as may be expected at finitetemperature. The result is a line defect in the form of a through-thickness (non-straight) vacancy column (in thefollowing for simplicity referred to as the void ), as shown in Figure 9( b ). We continue loading the simulationdomain in simple shear (moving the top and bottom faces relative to each other), while computing the effectiveapplied shear stress from the atomic forces. At sufficient applied shear, the pre-existing defect will nucleate andemit a dislocation dipole, whose activation energy and behavior depends on temperature. For an assessmentof the accuracy of the QC framework, we carry out both fully atomistic (125,632 atoms) and QC simulations(52,246 repatoms) in isothermal and adiabatic settings. We acknowledge that the QC setup is relatively simpleand there is not yet a significant reduction in the total number of degrees of freedom nor does it involve automaticmesh refinement. Yet, this study presents a simple and instructive example highlighting the efficacy and accuracyof the GPP-based QC formulation introduced in previous sections. In face-centered cubic (FCC) crystals, edge dislocations preferably glide on the close-packed crystallo-graphic { } -planes (Hull and Bacon, 2001). As the void generated due to the initial annihilation of thedislocation dipole is strained under shear deformation, dislocations nucleate from the void at a sufficient level22 igure 11: Comparison of the critical shear stress τ cr and strain γ cr required to nucleate a dislocation dipole from the voidas obtained from QC and from atomistic simulations at various temperatures. The critical strain γ cr is the external shearstrain γ at which τ cr is achieved. of applied shear, propagating in opposite directions, as shown in Figure 10. We apply a shear deformation to allrepatoms in the slab such that, at the n th load step, q ( n ) k, = q ( n − k, + ∆ γ q ( n − k, , (4.14)where indices and refer to the respective components of the mean position vector in the chosen coordinatesystem, and ∆ γ is the applied shear strain increment. As the strain is applied, repatoms in the inner region B(Figure 10) are relaxed while keeping those in the outer region A mechanically fixed to impose the average shearstrain. Note that, due to small deformation in the atomic neighborhoods, the displacement-variance entropy S Σ of repatoms close to the interface between regions A and B changes and, hence, all repatoms in the domainare thermally relaxed assuming an isothermal relaxation (cf. Table 1). While the shear strain is increased, thehorizontal component of the force on all repatoms in region A is computed. The effective shear stress τ on the (111) -plane is computed by normalizing the net horizontal force by the cross-sectional area of the slab. Resultsare shown in Figure 10 ( c ) and ( d ) . Once the stress reaches a critical value, the stress drops as dislocationsnucleate from the void and move to the ends of region B. We observe that the critical stress value decreasesslowly with temperature (see Figure 11). Moreover, the value of the critical stress obtained from a fully atomisticsimulation and the quasicontinuum simulation are within about 6% of each other (see Figure 11), both capturingthe apparent thermal plastic softening of the crystal.Miller and Tadmor (2009) studied a similar 2D scenario with a different crystal orientation, in which thedislocation dipole is stable and the dislocations do not annihilate to form the void. In such a case, the (theoretical)critical shear stress corresponds to the shear stress required to cause dislocation movement in the crystallographicplane. However, in our analysis (which is based on a realistic crystallographic setup) a critical shear stress maybe defined as the shear stress required to nucleate dislocations from the void-like defect. Figure 12: Variation of the position entropy S Σ with temperature inside one of the dislocations nucleated from the void. Thehighlighted atoms are identified using the centrosymmetry parameter (values > are shown). As discussed in Section 3.3,the number of neighbors (and their positions) affects the local interatomic potential of an atom, thus modifying the localvariation of positions. Atoms within the dislocations that are closer than the equilibrium interatomic spacing have smallerposition variance than those that are further apart. igure 13: Shear stress τ on the (111) -plane and dislocation separation distance d vs. applied shear strain γ for isothermaland adiabatic deformations obtained from both atomistic and QC simulations. Note that the differences between isothermalvs. adiabatic data are small, because the temperature increase is not significant. Critical shear stress value deviations arewithin 6% (e.g., isothermal QC: 3.7914 GPa, isothermal atomistics: 3.6389 GPa) and and critical strain values are within12% (isothermal QC: 0.112, isothermal atomistics: 0.100). To simulate the adiabatic activation of dislocations under shear, we repeat the above simulations, now withall repatoms in the domain being relaxed isentropically (cf. Table 1), followed by the thermal transport modelaccording to the steps discussed in Section 3.4 and Algorithm 1. We further assume that the strain-rate issingificantly lower than the time scale imposed by the molecular thermal transport, thus imposing quasistaticconditions for the transport ( δt ( n ) → ∞ in Algorithm 1, see Step 3 at the end of Section 3.4). The initialcondition (again, prepared using the isotropic elastic displacement field solutions) is relaxed isothermally at K, before the adiabatic deformation begins. We compare the adiabatic deformation with the isothermaldeformation of the slab at
K. Figure 13 shows the variation of the shear stress on the (111) -plane withexternal shear strain γ . Due to the mechanical work done by the external shear deformation, the temperature ofthe slab increases slightly, causing apparent softening compared to the isothermal deformation. Figure 14 showsthe spatial variation of temperature as the slab is deformed adiabatically. Before the critical strain, heating causedby local deformation is negligible. As the dislocations are nucleated at the critical strain γ cr , the temperatureof those atoms around the dislocations changes significantly, as shown by the intermediate stage ( q ∗ , T ∗ ) inFigure 14. Further plastic deformation causes increased heating of the slab, particularly due to the restricteddislocation motion beyond the interfaces between regions A and B (Figure 10).As noted above, the critical stress values obtained from QC simulations are within % of those obtained from Figure 14: Local variation of temperature as the slab with void is deformed adiabatically. In the intial stages, the temper-ature rises slowly due to the external deformation. At the critical shear strain γ cr = n cr ∆ γ , dislocations nucleate fromthe void and move along the [110] -direction. The temperature rise of those atoms within the dislocations causes a rapidincrease in the temperature of the slab, as may be expected from the heat generated by work hardening. The temperature ofa few atoms within the dislocations at the intermediate stage n = n ∗ cr before the irreversible transport exceeds the colorbarrange. %. Repatoms on the (110) -and (111) -surfaces, which includes the repatoms on vertices of very large elements and the repatoms in thetransition region between regions A and B, exhibit both mechanical and thermal spurious forces in 3D (Amelanget al., 2015; Amelang and Kochmann, 2015; Tembhekar et al., 2017). These artifacts are expected to be theprimary sources of error in the coarse-graining strategy adopted here. Mechanical spurious forces within theenergy-based fully nonlocal QC setup were discussed in detail in Amelang et al. (2015), and thermal spuriousforces are expected to show an analogous behavior. Therefore, we do not study spurious forces here in detail tomaintain the focus on the thermomechanics of the GPP formulation. Finally, we apply the discussed thermomechanical transport model to the case of nanoindentation into a Cusingle-crystal. While the problem is well studied in 2D using finite-temperature QC implementations (Tadmoret al., 2013), only few QC studies exist studying the finite-temperature effects in 3D. Kulkarni (2007) studied thenanoindentation of a × × unit-cell FCC nearest-neighbour Lennard-Jones crystal, using the Wentzel-Kramers-Brillouin (WKB) approximation, which captures the thermomechanical coupling at comparably lowtemperature only. We here perform nanoindentation simulations of a Cu cube of . µ m side length (approxi-mately 215 × ×
215 unit cells, see Figure 15), modeled by the EAM potential of Dai et al. (2006), underneatha nm spherical indenter modeled using the indenter potential of Kelchner et al. (1998) with a force constantof eV/ ˚A and a maximum displacement of . nm or approximately . times the lattice parameter at K. The crystal consists of approximately 50 million atoms, represented in the QC framework by approximately0.2 million repatoms. The top surface is modeled as a free boundary, while all other boundaries suppress wall-normal displacements, allowing only in-plane motion. Below, we discuss the results for isothermal deformationat T = 0 K ,
300 K ,
600 K and for quasistatic adiabatic deformation, the latter being initially at T = 300 K .Figure 16 shows the variation of the total force on the spherical indenter vs. indentation depth for bothisothermal and quasistatic adiabatic conditions. The force increases nonlinearly with indentation depth, showingthe typical Hertzian-type initial elastic section. With increasing indentation depth, atoms underneath the indentergenerate dislocations and stacking faults, overall creating a complex microstructure, as shown in Figure 17.At the first dislocation activation, the indenter force drops after reaching a critical force. This critical forcedecreases with increasing temperature, indicating plastic softening with increasing temperature. The dislocationsmove towards the boundaries of the atomistic region, gliding in the crystallographic planes of the (111) -family, ( a ) ( b )( c ) Figure 15: Illustration of the undeformed and deformed QC meshes of a . µ m FCC single-crystal of pure Cu. ( a ) Cross-section of the initial, undeformed mesh, ( b, c ) zoomed-in perspective views of the atomistic region (33 × ×
33 unitcells) and the surrounding coarsened regions in the ( b ) undeformed and ( c ) deformed configurations underneath a 5 nmspherical indenter at an indentation depth of . nm. igure 16: Variation of the indenter force with the indenter depth for isothermal ( T = 0 K ,
300 K ,
600 K ) and adiabatic(initially at T = 300 K ) conditions. After an initial elastic region, the curve shows the typical serrated behavior due todislocation activity underneath the indenter. giving way to stacking faults in those planes. As shown in Figure 17, while the initial dislocations maintaintheir structure, the stacking fault structure changes significantly at the same indentation depth as temperatureincreases.Figure 18 shows the spatial temperature distribution and the emergent microstructure during the quasistaticadiabatic deformation of the Cu crystal. With the dislocations local temperature gradients are generated due tolarge gradients in S Σ , which are triggered due to large deviations from a centrosymmetric neighborhood. Thesetemperature gradients are diffused as a result of the thermal transport.We note that, for a thorough quantitative analysis, one may want to obtain results averaged over multiple sim-ulations with initial conditions and/or indenter position slightly perturbed, since the emergence of microstructurebelow the indenter within the highly-symmetric single-crystal is associated with instability and strongly dependson local fluctuations and initial conditions. Such a statistical analysis is deferred to future work.
5. Conclusion and discussion
We have presented a Gaussian phase packets-based (GPP-based) formulation of finite-temperature equilib-rium and nonequilibrium thermomechanics applied to atomistic systems. We have shown that approximating theglobal statistical distribution function with a multivariate Gaussian enables capturing of thermal transport onlyvia interatomic correlations. Due to high computational costs, we have neglected the interatomic correlations,which results in a local GPP approximation of the system. Such a system exhibits reversible dynamics with
Figure 17: Microstructure generated below the 5 nm spherical indenter at an indentation depth of nm for isothermaldeformation of a Cu single-crystal. The generated dislocations move towards the boundaries of the atomistic region,creating stacking faults in the crystallographic planes of the { } -family (shaded in gray). Red repatoms are top surfacerepatoms and blue atoms denote those within the microstructure, identified using the centrosymmetry parameter (values > identify boundary and microstructure atoms only). All repatoms are shown with reduced opacity for comparison ofsize of the microstructure with the atomistic domain. a ) ( b ) Figure 18: Adiabatic deformation of the Cu single-crystal under a spherical indenter. ( a ) Spatial variation of temperaturein the cross-section of the atomistic region. As dislocations and stacking faults are created, local thermal gradients aregenerated which are diffused via the thermal transport model. ( b ) Microstructure generated below the spherical indenterat an indentation depth of nm. The generated dislocations move towards the boundaries of the atomistic region, creatingstacking faults in the crystallographic planes of the { } -family (shaded in gray). Red repatoms are top surface repatomsand blue atoms denote those within the microstructure, identified using the centrosymmetry parameter (values > identifyboundary and microstructure atoms only). Shaded repatoms denote all repatoms with reduced opacity, shown here forcomparison of size of the microstructure with the atomistic domain. thermomechanical coupling, causing local heating and cooling upon movement of atoms relative to local neigh-borhoods. Moreover, in the quasistatic limit we have shown that the equations yield local mechanical and thermalequilibrium conditions, the latter yielding the local equation of state of the atoms based on the interatomic forcefield. To capture the irreversibility due to local thermal transport triggered by the adiabatic heating/cooling ofatoms, we have coupled the quasistatic framework with linear Onsager kinetics. Such a model involves an em-pirical coefficient fitted to obtain approximate bulk conductivity measurements and captures the experimentallyobserved size-effects of the thermal conductivity, as shown by Venturini et al. (2014). Moreover, we have shownthat the time scale imposed by the atomic-scale transport is of the same order as that of atomic vibrations. Whilethe atomic-scale thermal transport imposes a small time scale, as the system reaches a non-uniform steady state,the local heat flux imbalance decreases. Below a tolerance value, the heat transport can be terminated, yielding asteady state solution. Based on the global multivariate Gaussian, interatomic correlations may be fitted to obtaincorrelation functions (akin to interatomic potentials), which can help develop the transport constitutive propertiesof atomistic systems and also advance current understanding of the long-standing nanoscale thermal transportproblem. Finally, we have combined the quasistatic thermomechanical equations based on the local GPP ap-proximation with thermal transport in a high-performance, distributed-memory, updated-Lagrangian 3D QCsolver, which is capable of modeling thermomechanical deformation of large-scale systems by coarse-grainingthe atomistic ensemble in space. Benchmark simulations of dislocation nucleation and nanoindentation underisothermal and adiabatic conditions showed convincing agreement between coarse-grained and fully resolvedatomistic simulations. Since the time integration of the atomic transport can be terminated for small heat fluximbalance (discussed in Algorithm 1), the quasistatic simulations offer significant advantages over traditionalMD studies, which can tackle only high strain rates. The presented methodology of coupling local thermal equi-librium with a surrogate empirical model of thermal transport and spatial coarse-graining (by the QC method)can model deformation of large crystalline systems at mesoscales and at quasistatic loading rates. Due to thetime-scale limitations of MD, a one-to-one comparison of the presented simulations with finite-temperature MDsimulations is prohibitively costly. A detailed analysis of the accuracy of the spatial coarse-graining of the ther-momechanical model presented here, and a comparison with suitable MD simulations. qualifies as a possibleextension of this work. For such comparisons, however, very large scale nonequilibrium molecular dynamics(NEMD) simulations are required at sufficiently slow strain rates. Acknowledgments
The support from the European Research Council (ERC) under the European Union’s Horizon 2020 researchand innovation program (grant agreement no. 770754) is gratefully acknowledged. The authors thank MichaelOrtiz for stimulating discussions. 27 eferences
Admal, N. C., Tadmor, E. B., 2010. A unified interpretation of stress in molecular systems. Journal of Elasticity100 (1-2), 63–143.Amelang, J. S., 2016. A fully-nonlocal energy-based formulation and high-performance realization of the qua-sicontinuum method. Ph.D. thesis.Amelang, J. S., Kochmann, D. M., 2015. Surface effects in nanoscale structures investigated by a fully-nonlocalenergy-based quasicontinuum method. Mechanics of Materials 90, 166 – 184.Amelang, J. S., Venturini, G. N., Kochmann, D. M., 2015. Summation rules for a fully nonlocal energy-basedquasicontinuum method. Journal of the Mechanics and Physics of Solids 82, 378–413.Ariza, M., Romero, I., Ponga, M., Ortiz, M., 2012. Hotqc simulation of nanovoid growth under tension in copper.International Journal of Fracture 174 (1), 75–85.Butcher, P., 1986. The theory of electron transport in crystalline semiconductors. In: Crystalline SemiconductingMaterials and Devices. Springer, pp. 131–194.Chang, Y., Himmel, L., 1966. Temperature dependence of the elastic constants of cu, ag, and au above roomtemperature. Journal of Applied Physics 37 (9), 3567–3572.Dai, X., Kong, Y., Li, J., Liu, B., 2006. Extended finnis–sinclair potential for bcc and fcc metals and alloys.Journal of Physics: Condensed Matter 18 (19), 4527.Dobson, M., Luskin, M., Ortner, C., 2010. Accuracy of quasicontinuum approximations near instabilities. Jour-nal of the Mechanics and Physics of Solids 58 (10), 1741–1757.Eidel, B., Stukowski, A., 2009. A variational formulation of the quasicontinuum method based on energy sam-pling in clusters. Journal of the Mechanics and Physics of Solids 57 (1), 87–108.Espanol, M. I., Kochmann, D. M., Conti, S., Ortiz, M., 2013. A γ -convergence analysis of the quasicontinuummethod. Multiscale Modeling & Simulation 11 (3), 766–794.Evans, D. J., Morriss, G. P., 2007. Statistical Mechanics of Nonequilbrium Liquids. ANU Press.Gunzburger, M., Zhang, Y., 2010. A quadrature-rule type approximation to the quasi-continuum method. Multi-scale Modeling & Simulation 8 (2), 571–590.Heller, E. J., 1975. Time-dependent approach to semiclassical dynamics. The Journal of Chemical Physics 62 (4),1544–1555.Hicks, L., Dresselhaus, M. S., 1993. Thermoelectric figure of merit of a one-dimensional conductor. Physicalreview B 47 (24), 16631.Hirth, J. P., 1980. Effects of hydrogen on the properties of iron and steel. Metallurgical Transactions A 11 (6),861–890.Hull, D., Bacon, D. J., 2001. Introduction to dislocations. Butterworth-Heinemann.Iyer, M., Gavini, V., 2011. A field theoretical approach to the quasi-continuum method. Journal of the Mechanicsand Physics of Solids 59 (8), 1506–1535.Johnson, R., 1988. Analytic nearest-neighbor model for fcc metals. Physical Review B 37 (8), 3924.Kelchner, C. L., Plimpton, S., Hamilton, J., 1998. Dislocation nucleation and defect structure during surfaceindentation. Physical Review B 58 (17), 11085.Kim, S.-P., Datta, D., Shenoy, V. B., 2014. Atomistic mechanisms of phase boundary evolution during initiallithiation of crystalline silicon. The Journal of Physical Chemistry C 118 (31), 17247–17253.28nap, J., Ortiz, M., 2001. An analysis of the quasicontinuum method. Journal of the Mechanics and Physics ofSolids 49 (9), 1899–1923.Kulkarni, Y., 2007. Coarse-graining of atomistic description at finite temperature. Ph.D. thesis, California Insti-tute of Technology.Kulkarni, Y., Knap, J., Ortiz, M., 2008. A variational approach to coarse graining of equilibrium and non-equilibrium atomistic description at finite temperature. Journal of the Mechanics and Physics of Solids 56 (4),1417–1449.Landau, L., Lifshitz, E., 1980. Statistical Physics, Part 1: Volume 5.Lepri, S., Livi, R., Politi, A., 2003. Thermal conduction in classical low-dimensional lattices. Physics Reports377 (1), 1–80.Ma, J., Hsu, D., Straub, J. E., 1993. Approximate solution of the classical liouville equation using gaussianphase packet dynamics: Application to enhanced equilibrium averaging and global optimization. The Journalof Chemical Physics 99 (5), 4024–4035.Marian, J., Venturini, G. N., Hansen, B., Knap, J., Ortiz, M., Campbell, G., 2009. Finite-temperature non-equilibrium quasicontinuum method based on langevin dynamics. Modelling and Simulation in MaterialsScience and Engineering, vol. 18, N/A, December 10, 2001, pp. 015003 18 (LLNL-JRNL-412903).McLachlan, A., 1964. A variational solution of the time-dependent schrodinger equation. Molecular Physics8 (1), 39–44.Mehta, R., Chugh, S., Chen, Z., 2015. Enhanced electrical and thermal conduction in graphene-encapsulatedcopper nanowires. Nano letters 15 (3), 2024–2030.Mendez, J., Ponga, M., Ortiz, M., 2018. Diffusive molecular dynamics simulations of lithiation of siliconnanopillars. Journal of the Mechanics and Physics of Solids 115, 123–141.Miller, R., Tadmor, E., Phillips, R., Ortiz, M., 1998. Quasicontinuum simulation of fracture at the atomic scale.Modelling and Simulation in Materials Science and Engineering 6 (5), 607.Miller, R. E., Tadmor, E. B., 2009. A unified framework and performance benchmark of fourteen multiscaleatomistic/continuum coupling methods. Modelling and Simulation in Materials Science and Engineering17 (5), 053001.Nabarro, F. R. N., 1967. Theory of crystal dislocations. Clarendon Press.Nix, F., MacNair, D., 1941. The thermal expansion of pure metals: copper, gold, aluminum, nickel, and iron.Physical Review 60 (8), 597.Overton Jr, W., Gaffney, J., 1955. Temperature variation of the elastic constants of cubic elements. I. Copper.Physical Review 98 (4), 969.Ponga, M., Ortiz, M., Ariza, M., 2015. Finite-temperature non-equilibrium quasi-continuum analysis ofnanovoid growth in copper at low and high strain rates. Mechanics of Materials 90, 253–267.Qu, S., Shastry, V., Curtin, W., Miller, R. E., 2005. A finite-temperature dynamic coupled atomistic/discretedislocation method. Modelling and Simulation in Materials Science and Engineering 13 (7), 1101.S¨a¨askilahti, K., Oksanen, J., Volz, S., Tulkki, J., 2015. Frequency-dependent phonon mean free path in carbonnanotubes from nonequilibrium molecular dynamics. Physical Review B 91 (11), 115426.Srivastava, A., Nemat-Nasser, S., 2014. On the limit and applicability of dynamic homogenization. Wave Motion51 (7), 1045–1054.Stroud, A. H., 1971. Approximate calculation of multiple integrals. Prentice-Hall.29admor, E. B., Legoll, F., Kim, W., Dupuy, L., Miller, R. E., 2013. Finite-temperature quasi-continuum. AppliedMechanics Reviews 65 (1).Tadmor, E. B., Miller, R. E., 2011. Modeling materials: continuum, atomistic and multiscale techniques. Cam-bridge University Press.Tadmor, E. B., Ortiz, M., Phillips, R., 1996a. Quasicontinuum analysis of defects in solids. Philosophical Mag-azine A 73 (6), 1529–1563.Tadmor, E. B., Phillips, R., Ortiz, M., 1996b. Mixed atomistic and continuum models of deformation in solids.Langmuir 12 (19), 4529–4534.Tembhekar, I., 2018. The fully nonlocal, finite-temperature, adaptive 3d quasicontinuum method for bridgingacross scales. Ph.D. thesis, California Institute of Technology.Tembhekar, I., Amelang, J. S., Munk, L., Kochmann, D. M., 2017. Automatic adaptivity in the fully nonlocalquasicontinuum method for coarse-grained atomistic simulations. International Journal for Numerical Meth-ods in Engineering 110 (9), 878–900.Tuckerman, M., 2010. Statistical mechanics: theory and molecular simulation. Oxford University Press.van der Giessen, E., Schultz, P. A., Bertin, N., Bulatov, V. V., Cai, W., Cs´anyi, G., Foiles, S. M., Geers, M. G.,Gonz´alez, C., H¨utter, M., et al., 2020. Roadmap on multiscale materials modeling. Modelling and Simulationin Materials Science and Engineering 28 (4), 043001.Venturini, G., Wang, K., Romero, I., Ariza, M., Ortiz, M., 2014. Atomistic long-term simulation of heat andmass transport. Journal of the Mechanics and Physics of Solids 73, 242–268.Venturini, G. N., 2011. Topics in multiscale modeling of metals and metallic alloys. Ph.D. thesis, CaliforniaInstitute of Technology.Weiner, J. H., 2012. Statistical mechanics of elasticity. Courier Corporation.Ziman, J. M., 2001. Electrons and phonons: the theory of transport phenomena in solids. Oxford UniversityPress. Appendix A. Time evolution of phase averaged quantities
The time evolution of the phase average of a phase-space quantity A ( z ) can be derived using the representa-tive solution of the Liouville equation, ∂f∂t + i L f = 0 = ⇒ f ( z , t ) = e − i L t f ( z ,
0) = e − i L t f ( z ) , (A.1)where f ( z , is the initial condition and e − i L t is the propagating operator, which transforms the probabilitydistribution initially defined at phase-space coordinate z to the probability distribution at z ( t ) (Evans andMorriss, 2007). Furthermore, the time evolution of the phase-space quantity A ( z ) is given by d A d t = i L A = ⇒ A ( z , t ) = e i L t A ( z ,
0) = e i L t A ( z ) . (A.2)Equations (A.1) and (A.2) reveal that the operators e ± i L t transport the probability distribution f ( z ) and phase-space quantities A ( z ) defined in terms of z of a system of particles, given that z also changes in time as thesystem of particles evolves.Operator i L is self-adjoint and satisfies the property (cid:90) Γ A ( z ) i L B ( z ) d z = (cid:90) Γ ( − i L ) A ( z ) B ( z ) d z (A.3)30or real-valued A and B and Γ ⊆ R N . Using this property, the time evolution of the phase average of aphase-space quantity A ( z ) , defined in terms of z , is obtained from N ! h N d (cid:104) A (cid:105) dt = dd t (cid:90) Γ f ( z , t ) A ( z ) d z = (cid:90) Γ dd t (cid:2) e − i L t f ( z ) (cid:3) A ( z ) d z = (cid:90) Γ e − i L t f ( z ) i L tA ( z ) d z . (A.4)Using equation (A.2), we obtain d (cid:104) A (cid:105) d t = 1 N ! h N (cid:90) Γ f ( z , t ) d A d t d z = (cid:28) d A d t (cid:29) . (A.5)We note that equation (A.5) is obtained using only the evolution equations (A.1), (A.2) and property (A.3),which hold for any Hamiltonian system and thus contain no time-coarsening approximations. Accordingly, time-variational formulations such as the Frenkel-Dirac-McLachlan variational principle (McLachlan, 1964) leads toidentical equations. Appendix B. Quasistatic GPP as Helmholtz free energy minimization
The Helmholtz free energy F as a function of parameter set ( q , S Σ , S Ω ) is defined as F ( q , S Σ , S Ω ) = E ( q , S Σ , S ) − (cid:88) i Ω i S i k B m i (B.1)with the relation Ω i k B m i = ∂E∂S i . (B.2)Minimization of F with respect to the set q yields − ∂ F ∂ q i = 0 = ⇒ − (cid:28) ∂V ( q ) ∂ q i (cid:29) = (cid:104) F i ( q ) (cid:105) = 0 . (B.3)To minimize F with respect to the set S Σ , we consider the relation q i − q i = (cid:112) Σ i x i , (B.4)for some normalized vector x i . From equation (B.4) it follows that ∂ q i ∂S Σ ,i = ∂ q i ∂ √ Σ i ∂ √ Σ i ∂S Σ ,i = (cid:112) Σ i x i = q i − q i . (B.5)Finally, minimization of F with respect to the set S Σ yields the following set of equations: − ∂ F ∂S Σ ,i = 0 = ⇒ i m i − (cid:28) ∂V ( q ) ∂ q i · ∂ q i ∂S Σ ,i (cid:29) = 3Ω i m i + (cid:104) F i ( q ) · ( q i − q i ) (cid:105) = 0 . (B.6)Equations (B.3) and (B.6) are identical to the quasistatic GPP equations (3.12). Appendix C. Time step stability bounds for entropy transport
Applying a forward-Euler explicit time discretization to equation (3.30), we obtain t ( k ) (cid:32)(cid:18) T i T j (cid:19) ( k +1) − (cid:18) T i T j (cid:19) ( k ) (cid:33) = 2 A k B T , ( k ) ij T ( k ) i T ( k ) j (cid:18) T i − T j T j − T i (cid:19) ( k ) , (C.1)where superscript ( k ) implies a quantity evaluated at the k th time step. Rearranging the above equation yields (cid:18) T i T j (cid:19) ( k +1) = T ( k ) (cid:18) T i T j (cid:19) ( k ) , (C.2)31here T ( k ) is the transition matrix at the k th time step, defined by T ( k ) = I + 2 A ∆ t ( k ) k B T , ( k ) ij T ( k ) i T ( k ) j (cid:18) − − (cid:19) . (C.3)For numerical stability, the transition matrix must have eigenvalues with magnitude ≤ , which yields the bound A ∆ t ( k ) k B (cid:32) T , ( k ) ij T ( k ) i T ( k ) j (cid:33) ≤ . (C.4)Applying the above limit to a system in which the i thth