Journal of the Mechanics and Physics of Solids Nonequilibrium thermomechanics of Gaussian phase packet crystals: Application to the quasistatic quasicontinuum method

The quasicontinuum (QC) method was originally introduced to bridge across length scales by coarse-graining an atomistic ensemble to significantly larger continuum scales at zero temperature, thus overcoming the crucial length-scale limitation of classical atomic-scale simulation techniques while solely relying on atomic-scale input (in the form of interatomic potentials). An associated challenge lies in bridging across time scales to overcome the time- scale limitations of atomistics at finite temperature. To address the biggest challenge, bridging across both length and time scales, only a few techniques exist, and most of those are limited to conditions of constant temperature. Here, we present a new general strategy for the space–time coarsening of an atomistic ensemble, which introduces thermomechanical coupling. Specifically, we evolve the statistics of an atomistic ensemble in phase space over time by applying the Liouville equation to an approximation of the ensemble’s probability distribution (which further admits a variational formulation). To this end, we approximate a crystalline solid as a lattice of lumped correlated Gaussian phase packets occupying atomic lattice sites, and we investigate the resulting quasistatics and dynamics of the system. By definition, phase packets account for the dynamics of crystalline lattices at finite temperature through the statistical variances of atomic momenta and positions. We show that momentum–space correlation allows for an exchange between potential and kinetic contributions to the crystal’s Hamiltonian. Consequently, local adiabatic heating due to atomic site motion is captured. Moreover, in the quasistatic limit, the governing equations reduce to the minimization of thermodynamic potentials (similar to maximum-entropy formulation previously introduced for finite-temperature QC), and they yield the local equation of state, which we derive for isothermal, isobaric, and isentropic conditions. Since our formulation without interatomic correlations precludes irreversible heat transport, we demonstrate its combination with thermal transport models to describe realistic atomic- level processes, and we discuss opportunities for capturing atomic-level thermal transport by including interatomic correlations in the Gaussian phase packet formulation. Overall, our Gaussian phase packet approach offers a promising avenue for finite-temperature nonequilibrium quasicontinuum techniques, which may be combined with thermal transport models and extended to other approximations of the probability distribution as well as to exploit the variational structure. equation, which yields the equations of motion for the phase space parameters. We also show that the interatomic correlations are of fundamental importance for modeling interatomic heat flux. Since such correlations increase the degrees of freedom significantly, they are neglected in subsequent sections while retaining 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 vanishing limit in quasistatics. We show that, in the latter limit, the governing equations define the local thermomechanical equilibrium of the system and yield a rate-independent local thermal equation of state. To capture thermal transport, we adopt Venturini et al.’s linear Onsager kinetics-based model. We also discuss the time scale imposed by the Onsager


A B S T R A C T
The quasicontinuum (QC) method was originally introduced to bridge across length scales by coarse-graining an atomistic ensemble to significantly larger continuum scales at zero temperature, thus overcoming the crucial length-scale limitation of classical atomic-scale simulation techniques while solely relying on atomic-scale input (in the form of interatomic potentials). An associated challenge lies in bridging across time scales to overcome the timescale limitations of atomistics at finite temperature. To address the biggest challenge, bridging across both length and time scales, only a few techniques exist, and most of those are limited to conditions of constant temperature. Here, we present a new general strategy for the space-time coarsening of an atomistic ensemble, which introduces thermomechanical coupling. Specifically, we evolve the statistics of an atomistic ensemble in phase space over time by applying the Liouville equation to an approximation of the ensemble's probability distribution (which further admits a variational formulation). To this end, we approximate a crystalline solid as a lattice of lumped correlated Gaussian phase packets occupying atomic lattice sites, and we investigate the resulting quasistatics and dynamics of the system. By definition, phase packets account for the dynamics of crystalline lattices at finite temperature through the statistical variances of atomic momenta and positions. We show that momentum-space correlation allows for an exchange between potential and kinetic contributions to the crystal's Hamiltonian. Consequently, local adiabatic heating due to atomic site motion is captured. Moreover, in the quasistatic limit, the governing equations reduce to the minimization of thermodynamic potentials (similar to maximum-entropy formulation previously introduced for finite-temperature QC), and they yield the local equation of state, which we derive for isothermal, isobaric, and isentropic conditions. Since our formulation without interatomic correlations precludes irreversible heat transport, we demonstrate its combination with thermal transport models to describe realistic atomiclevel processes, and we discuss opportunities for capturing atomic-level thermal transport by including interatomic correlations in the Gaussian phase packet formulation. Overall, our Gaussian phase packet approach offers a promising avenue for finite-temperature nonequilibrium quasicontinuum techniques, which may be combined with thermal transport models and extended to other approximations of the probability distribution as well as to exploit the variational structure.

Introduction
Crystalline solids exhibit physical and chemical transport phenomena across wide ranges of length and time scales. This includes the transport of charges (Butcher, 1986;Ziman, 2001), heat (Ziman, 2001), and mass (Weiner, 2012), as well as mechanical failure. Understanding such phenomena is crucial from both a fundamental scientific standpoint as well as to further advance technologies ranging from solid-state batteries (Kim et al., 2014a) to thermal management systems (Hicks and Dresselhaus, 1993) to failureresistant metallic structural components (Hirth, 1980) -all exposed to complex dynamic conditions varying over time scales of a few microseconds to several hours and length scales of a few nanometers to a meter. Such a variety of length and time scales underlying the transport phenomena calls for simulation techniques that capture the physical processes across all length and time scales involved. While continuum mechanics and related finite element (FE) and 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 various transport phenomena at atomic-level length scales (angstroms to tens of nanometers) and time scales (femto-to nanoseconds) (Tuckerman, 2010). Where a higher level of accuracy is required, methods such as Density Functional Theory (DFT) or the direct computation of Schrödinger's equation have aimed at capturing the quantum coupling of molecular-level physics. All of the aforementioned techniques specialize in the approximate ranges of length and time scales mentioned above. However, each of those techniques poses restrictive assumptions at smaller scales while becoming prohibitively costly at larger scales, hence making scale-bridging techniques attractive (Srivastava and Nemat-Nasser, 2014;van der Giessen et al., 2020). For instance, the state-of-the-art DFT-based multiscale modeling techniques developed by Motamarri et al. (2020) make ab-initio accuracy available to larger atomic ensembles, simulating systems consisting of approximately 4000 atoms with a large computational cost met by massively parallel supercomputers.
Several concurrent scale-bridging techniques have been developed over the past few decades, particularly focusing on multiscale thermomechanical modeling of crystalline materials (cf. Xu and Chen, 2019 for a detailed review of some of the available techniques). The atomistic-to-continuum (AtoC) method developed by Wagner et al. (2008) utilizes a coupling between a predetermined atomistic domain and an overlaid continuum domain, discretized using a suitable FE method. The continuum domain is used for providing a heat bath to the atomistic domain, simultaneously minimizing the difference between the temperature fields in both domains, hence establishing a two-way coupling. In the bridging-domain method (BDM), the pre-determined atomistic and continuum subdomains are coupled by imposing a weak displacement compatibility condition at the intersecting nodes (Belytschko and Xiao, 2003). Nodal mechanical forces are computed using the total Hamiltonian of the system constructed using the atomistic and continuum subdomains as well as the weak compatibility conditions. Chen (2009) developed the concurrent atomistic-continuum (CAC) method in which microscopic balance equations, derived using the theory of Irving and Kirkwood (1950), are solved to determine the thermomechanical deformation of atomic sites within large continuum-scale subdomains. Unlike the AtoC and BDM methods, the CAC method can capture lattice defects, such as dislocations and cracks, in the continuum subdomain. The coupled atomistic/discrete-dislocation (CADD) method of Shilkrot et al. (2002) is another multiscale method that allows the movement of dislocation defects across the atomistic and continuum subdomains. While CAC requires only the interatomic potential as the constitutive input, CADD requires reduced-order continuum constitutive models to determine long-range elastic stress fields. All of the aforementioned methods require a-priori knowledge of atomistic and continuum subdomains, and they solve the thermomechanical deformation in both subdomains using distinct methodologies. This, unfortunately, inhibits a seamless transition from atomistic-to continuum-scale subdomains. Furthermore, finite-temperature variants of all those techniques typically use timeintegration at the scale of atomic vibrations to resolve the available phonon modes in the domain (Chen et al., 2017), thus not allowing for temporal coarsening.
The quasicontinuum (QC) method of Tadmor et al. (1996a) is a concurrent scale-bridging method that aims to solve the problem of spatial upscaling from atomistic to continuum length scales via intermediate mesoscopic scales (Miller et al., 1998). Starting from a standard atomistic setting, a carefully selected set of representative degrees of freedom (associated with repatoms) reduces the computational costs and admits simulating continuum-scale problems by restricting atomistic resolution to where it is in fact needed. Different flavors of QC have been proposed based on the interpolation of forces on repatoms (Knap and Ortiz, 2001), or based on approximating the total Hamiltonian of 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 formulation of  enables a seamless spatial scale-bridging within the QC setup, which does not require a strict separation of (nor apriori knowledge about) atomistic and non-atomistic (coarsened) subdomains within the simulation domain. The capabilities of this fully nonlocal QC technique have been demonstrated, e.g., by large-scale quasistatic total-Lagrangian simulations of dislocation interactions during nano-indentation , nanoscale surface and size effects , and void growth and coalescence (Tembhekar et al., 2017). In this work, we adopt the theoretical nonlocal, energy-based QC formulation of  in a new, updated Lagrangian framework and implementation for spatial upscaling.
While spatial coarse-graining is thus achieved by approximating an atomic ensemble by a subset of atoms, temporal coarsegraining requires approximate modeling due to the unavailability of the trajectories of all the atoms at a given time. Furthermore, the Hamiltonian dynamics of an ensemble of atoms couples length and time scales in the system (Evans and Morriss, 2007), which is why spatial scale-bridging techniques have often been applied to systems at zero temperature or a uniform temperature only (Tadmor et al., 2013). One way to model uniform temperature is to apply a global ergodic assumption for every atomic site, thus yielding space-averaged trajectories of atoms as phase-averaged trajectories. Suitable global thermal equilibrium distributions (such as NVT ensembles Tadmor and Miller, 2011) are used for phase averaging of trajectories. The motion of atoms on these phase-averaged P. Gupta et al. trajectories is governed by the phase averages of interatomic potentials and kinetic energy. Due to thermal softening of interatomic potentials upon phase-averaging, the accessible time-scales in isothermal dynamic simulations increase, hence enabling time coarsening. Kim et al. (2014b) further increased the accessible time limits of Tadmor et al.'s method using hyperdynamics (Voter, 1997) to capture thermally activated rare events, such as atomic-scale mass diffusion. However, most prior work has been restricted to local harmonic approximations of interatomic potentials and isothermal deformation (Tadmor and Miller, 2011;Tadmor et al., 2013) at a uniform temperature. Using Langevin dynamics is an alternative strategy, in which a stochastic thermal forcing is added to 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 (uniform temperature), which is why it is computationally costly. Kulkarni et al. (2008) introduced a variational mean-field approach for modeling non-uniform temperature, which approximates the global distribution function of the ensemble as a product of local entropy-maximizing (or max-ent ) distributions, constraining the local frequency of atomic vibrations by 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 temperature and 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 FE setup of the QC formulation and empirical bulk-thermal conductivity values. Venturini et al. (2014) extended the max-ent approach to 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. Combining Venturini et al.'s max-ent-based approach to include non-uniform interstitial concentrations, Mendez et al. (2018) used an Arrhenius-type master-equation model to achieve diffusive transport in long-term atomistic simulations. In the specific case of isothermal diffusion of impurities in crystals, the maxent framework of Venturini et al. (2014) resembles the diffusive molecular dynamics (DMD) formulation of Li et al. (2011), who used an isotropic Gaussian atomic density cloud to model uncorrelated vibrations of atoms.
We here present a similar nonequilibrium finite-temperature formulation which approximates the global distribution function of the ensemble based on Gaussian Phase Packets (GPP), which is a different ansatz from the max-ent one and exhibits time evolution governed by the Liouville equation. Previously, Ma et al. (1993) studied the dynamic Gaussian Phase Packet (GPP) ansatz as a trajectory-sampling technique for uniform-temperature MD simulations. They approximated the distribution function of each atom in an ensemble as a correlated Gaussian distribution with the covariance matrix as the mean-field or phase space parameter. Such an approximation yields the evolution equations of the covariance matrix by either directly substituting the GPP ansatz into the Liouville equation (strong form) or by using the Frenkel-Dirac-McLachlan variational principle (McLachlan, 1964). The resulting equations -combined with appropriate ergodic assumptions -may be integrated in time 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 relaxing towards a state of thermodynamic equilibrium with uniform temperature. Their work was inspired by the application of GPPs in quantum mechanics by Heller (1975), who used it for calculating parametrized solutions of the Schrödinger's equation.
In this work, we apply the GPP ansatz to an ensemble of atoms to elucidate the nonequilibrium and (local-thermal) equilibrium thermomechanical quasistatics and dynamics of the system. We show that an approximation of interatomic correlations is required for modeling atomistic-level transport phenomena. In an ensemble of atoms, such an approximation increases the degrees of freedom by ( 2 ), which is computationally costly. Therefore, we assume uncorrelated atoms or interatomic independence. Our formulation may be considered as a dynamic extension of the max-ent methods of Kulkarni et al. (2008) and Venturini et al. (2014) as well as the DMD approach of Li et al. (2011), to which our theory reduces in the quasistatic limit. We show that incorporating momentum-displacement correlations in a Gaussian ansatz for the distribution function elucidates the local dynamics of atoms and the energy exchange from kinetic to potential, thus dynamically capturing the thermomechanical deformation. When all atoms are assumed to be uncorrelated (such a solid crystal of independent atoms is also known as an Einstein's solid), the GPP approximation fails to capture thermal transport due to non-uniform thermomechanical deformation and hence requires additional phenomenological modeling to that end. Moreover, to achieve temporal coarsening, our GPP approach highlights the importance of vanishing interatomic correlations, approaching the quasistatic approximation.
To simulate irreversible heat flux, we combine the GPP framework within the quasistatic approximation with Venturini et al.'s linear Onsager kinetics to model local thermal transport. Although our GPP approach under the quasistatic approximation resembles quasistatic max-ent (Kulkarni et al., 2008;Venturini et al., 2014), our correlated Gaussian ansatz highlights the physical significance of assuming vanishing cross-correlations across all degrees of freedom. Furthermore, it emphasizes the need for additional thermodynamic assumptions required for modeling the thermomechanical deformation of the crystal, due to the loss of knowledge of the temporal evolution of the correlations. Our approach furthermore furnishes the quasistatic isothermal max-ent setup with alternative formulations for isobaric and isentropic conditions. The remainder of this study is structured as follows. In Section 2 we review the fundamentals of nonequilibrium statistical mechanics based on the Liouville equation and the GPP approximation. We derive the weak formulation using the variational structure of the Liouville equation, which yields the equations of motion for the phase space parameters. We also show that the interatomic correlations are of fundamental importance for modeling interatomic heat flux. Since such correlations increase the degrees of freedom significantly, they are neglected in subsequent sections while retaining 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 vanishing limit in quasistatics. We show that, in the latter limit, the governing equations define the local thermomechanical equilibrium of the system and yield a rate-independent local thermal equation of state. To capture thermal transport, we adopt Venturini et al.'s linear Onsager kinetics-based model. We also discuss the time scale imposed by the Onsager P. Gupta et al. Fig. 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 the atomic sites is modeled using the linear Onsager kinetics of Venturini et al. (2014). kinetics-based model and its time-stepping constraints. In Section 4 we discuss our QC implementation of the local thermomechanical equilibrium equations combined with thermal transport and demonstrate its use in a new, update-Lagrangian, distributed-memory QC solver for coarse-grained atomistic simulations. Finally, Section 5 concludes our analysis and discusses limitations and future prospects.

Nonequilibrium thermodynamics of Gaussian phase packets
We begin by discussing the nonequilibrium modeling of deformation mechanics of crystalline solids, using the GPP approximation in which atoms are treated as Gaussian clouds, centered at the mean phase space positions of the atoms. We briefly review the application of the Liouville equation for analyzing the statistical evolution of a large ensemble of atoms subject to high-frequency vibrations. Such random vibrations are modeled by assuming local phase space coordinates (positions and momenta) drawn from local Gaussian distributions. We discuss the impact of assuming independent Gaussian distributions for each atom (termed interatomic independence hereafter) and the corresponding inability of the model to capture nonequilibrium thermal transport. Moreover, the independent Gaussian assumption allows us to formulate the dynamical system of equations governing the atoms and the corresponding mean-field parameters. In subsequent sections, we use the insights gained here to formulate isentropic and non-isentropic quasistatic problems of finite-temperature crystal deformation (see Fig. 1 for a schematic description).

Hamiltonian dynamics
The time evolution of an atomic crystal, modeled as an ensemble of classical particles, is fully characterized by the generalized positions = { ( ), = 1, … , } and momenta = { ( ), = 1, … , } of all particles, which evolve in time according to Hamilton's equations, is the Hamiltonian of the system, ( ) represents the potential field acting on all atoms in the system, ( ) = −∇ ( ) is the net force on the th atom, and is the mass of the th atom. In classical MD, Eqs. (2.1) are solved for the trajectories ( ( ), ( ) ) of all atoms over time. As a consequence, such simulations require femtosecond-level time step resolution (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 the phase space coordinate = ( ( ), ( )). P. Gupta et al. Here and in the following, we use ∈ R 6 for brevity to represent the momenta and positions of the particles in three dimensions (3D). It is convenient to introduce the distribution ( , ) such that the quantity is the probability of finding the system of atoms such that the position and momentum of the th atom lie within ( , + d ) and ( , + d ), respectively. The time evolution of the probability distribution ( , ) is governed by the Liouville equation (Evans and Morriss, 2007;Tadmor and Miller, 2011): with the Liouville operator  defined by anḋ= (̇,̇). 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 where {⋅, ⋅} defines the Poisson bracket (Landau and Lifshitz, 2013). We note that solving (2.4) combined with (2.6) at all times is equivalent to solving the equations of motion (2.1) with (2.2) (Evans and Morriss, 2007). In general, solving (2.4) requires a discretization of the 6 -dimensional phase space { ⊆ R 6 ∶ ∈ } and of time for a system of atoms in 3D. To avoid such a general discretization, which is computationally prohibitively expensive, one resorts to physically sensible approximations of the probability distribution ( , ). One such approximation is the GPP ansatz discussed in Section 2.3, 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 derivation of such equations of motion is straightforward if the mean-field parameters can be identified as functions of moments of the random variable . However, a general framework for obtaining such equations of motion can be derived using the variational structure of the Liouville equation, as outlined below.

Variational structure and weak formulation of the Liouville equation
The Liouville equation in (2.4) may be identified as the strong form of the equation of motion of ( , ) in phase space. The weak form of the Liouville equation can be written as where ∈ R is a test function in real space (the denoting a variation in the variational-calculus sense), is some instance of time, {⋅, ⋅} denotes the Poisson brackets, and ⊆ R 6 denotes the phase space. Analogous to the Hamiltonian setting of classical mechanics (Lanczos, 1986), (2.7) can be identified as the first variation of the functional with respect to , where we abbreviated  = ( , ). We refer to [] as the Liouville action. Integrating (2.8) by parts with respect to , we obtain (2.9) The above equation resembles the general definition of the action integral in classical mechanics (cf. Chapter VI in Lanczos, 1986) and analogously defines the Liouville Hamiltonian (, ) as,  (2.12) P. Gupta et al. Substituting (2.11) in (2.8), considering = , and enforcing stationarity with respect to , we obtain the following set of equations of motion for parameters ( ) ( = 1, … , ): The phase average of any function ( ) is denoted by ⟨ ( )⟩ and defined as where ℎ is Planck's constant and the factor ! ℎ 3 is introduced to normalize by the phase space volume (Landau and Lifshitz, 1980). Using (2.14), (2.13) can be rewritten as This system of equations yields the equations of motion of parameters ( ), which parametrize the trajectories of the system governed by (2.4) and can be considered as the classical analogue of the Frenkel-Dirac-McLachlan variational principle used in quantum mechanics (McLachlan, 1964). For parameters ( ), which can be directly written as phase-averaged functions of the random variable (of mathematical form ⟨ ( )⟩ where (⋅) is some smooth function), the equations of motion can be derived using the identity (Evans and Morriss, 2007;Zubarev, 1974) (see Appendix A for a brief discussion). Although not exploited in the following, the variational setting of the Liouville-based approach pursued in this study bears potential with regards to numerical solution schemes, as it allows for a systematic discretization and convergence study of solutions, among others.

Crystal lattice of Gaussian phase packets
The above framework requires an ansatz for the probability distribution. As discussed in Section 1, the mean-field approximationbased temporal-coarsening approaches (Kulkarni et al., 2008;Li et al., 2011;Venturini et al., 2014) start by constructing an ansatz for the distribution ( ) directly at steady state. If the Hamiltonian is constrained as in Venturini et al. (2014), then the distribution function ( ) is only a function of , satisfying  = 0. While these restrictions may admit quasistatic solutions, they do not provide insight into the dynamic evolution nor do they consider the nature of the process which led to a quasistatic solution. In addition, considering only variance constraints on the distribution function (Kulkarni et al., 2008;Li et al., 2011) enforces a steady state of the mean-field parameters, as shown in the section below.
As a new approach, we here start with a multivariate Gaussian phase packet (GPP) ansatz to the distribution function with no assumptions about vanishing correlations, and we systematically discuss the consequences of eliminating interatomic and crosscorrelations of the degrees of freedom. The GPP approximation was first introduced in the context of quantum mechanics by Heller (1975) and used for classical systems by Ma et al. (1993). Both Heller (1975) and Ma et al. (1993) substituted Gaussian distributions into the Liouville equation (2.4) (Schrödinger's equation for quantum systems) to obtain the dynamical evolution of the phase space parameters. Following their idea, we approximate the system-wide probability distribution ( , ) as a multivariate Gaussian distribution, i.e., where is the 6 × 6 covariance matrix composed of the interatomic and momentum-displacement correlations, represents the vector of all atoms' mean positions and momenta, and the partition function ( ) is defined by Phase averaging of using (2.14) with (2.17) confirms that ( ) = ⟨ ( )⟩. We further conclude that can be written as a block-matrix with components such that each block represents the correlation among displacements and momenta of atoms and .
Consequently, instead of time-resolving the positions and momenta of all atoms (as in MD), we may equivalently track the time evolution of the atomic ensemble through the phase space parameters ( ( ), ( ) ) , which includes the mean momenta and positions of all atoms as well as the covariance matrix. Since parameters ( ) and ( ) can be identified as phase functions, application of (2.16) to ( ) and ( ) yields the dynamical equations that govern their time evolution: (2.20) P. Gupta et al. Application of (2.15) to ( ) and ( ) yields identical equations, as derived in Appendix B. Let us further specify the second equation in (2.20). Writing the components of covariance matrix as ) , (2.21) and assuming that all atoms have the same mass without loss of generality, identity (2.16) yields the following time evolution equations for the submatrices identified above: . (2.22d) Eqs. (2.22) govern the evolution of the pairwise momentum and displacement correlations of atoms and , and they must be solved to obtain the interatomic correlations at all times. Eqs. (2.22c)-(2.22d) govern the thermomechanical coupling of the crystal. Note that we may identify the pairwise kinetic tensor ( , ) as a measure of temperature, so that (2.22b) and (2.22c) describe the evolution of the system-wide distribution as a result of unbalanced pairwise virial tensors and kinetic tensors (see Admal and Tadmor, 2010 for the tensor virial theorem). The virial tensor ⟨ ( ) ⊗ ( − )⟩ changes with changing displacement correlations of atoms due to varying extents of atomic vibrations ( , ) , thus coupling (2.22d) with (2.22b) and (2.22c). The right-hand side of (2.22a) resembles the tensor form of the interatomic heat current (Sääskilahti et al., 2015;Lepri et al., 2003) and changes with varying correlation matrices , , thus coupling (2.22a) with (2.22b) and (2.22c). Consequently, the imbalance between virial and kinetic tensors in Eqs. (2.22b) and (2.22c) drives the phase space motion of the system of particles, resulting in the time evolution of ( , ) and ( , ) .
In the following, it will be helpful to identify the entropy of the atomic ensemble as where is Boltzmann's constant, and 0 is a constant for a given system with a fixed number of atoms. The entropy rate of change follows as (2.24) Eqs. (2.20) govern the phase space motion of a system of atoms and contain a total of 6 + 36 ( + 1)∕2 evolution equations, solving which is even more computationally expensive than solving the state-space governing equations of MD. Therefore, as a simplifying assumption Ma et al. (1993) and Heller (1975) assumed the statistical independence of atoms (or of states in the quantum analogue), which implies = for ≠ . (2.25) Although this assumption severely limits the applicability of (2.20), it admits a beneficial analytical treatment of the time evolution problem and allows us to gain new insight into the nature of correlations and, especially, into the quasistatic limit. Therefore, we use (2.20) in the following. As we will demonstrate below, a consequence of the interatomic independence assumption (2.25) is that heat transport cannot be resolved correctly. Consequently, we may apply the phase space evolution Eqs. (2.20) and (2.22) only for isentropic (reversible) finite temperature simulations of quasistatic and dynamic processes.

Independent Gaussian phase packets
Since solving the full system of evolution equations is prohibitively expensive, as discussed above, we apply (2.25) and assume non-zero correlations between the position and momentum of each individual atom, but no cross-correlations between different atoms. To this end, we apply the GPP approximation to a single atom , which yields the multivariate Gaussian distribution of the phase space coordinate as where the phase space parameters ( , ) denote the mean phase space coordinate and variance of the th atom, respectively, defined as (2.27) P. Gupta et al.

Fig. 2.
Schematic illustration of an atomic trajectory by GPP dynamics (on the right) compared to an MD trajectory (on the left) for a single particle. Small-scale motions of the particle are approximated by parameters , , and upon averaging over suitable time intervals. As discussed in Section 3.2, in the quasistatic limit, and are proportional to the local temperature and → 0.
The normalization quantity ( ) may be identified as the single-particle partition function is the 6 × 6 covariance matrix of the multivariate Gaussian and accounts for the variance or uncertainty in the momentum and displacement of the th atom. The assumed interatomic independence eliminates the interatomic correlations as independent variables, thus reducing the total number of equations to 27 + 6 for a system of atoms, which govern the time evolution of the kinetic tensor ( , ) , displacement-correlation tensor ( , ) , and the momentum-displacement-correlation tensor combined with the phase-averaged equations of motion, For simplicity, we further make the spherical distribution assumption that all cross-correlations between different directions of momenta and displacements vanish, hence approximating the atomic clouds in phase space as spherical (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.29), we obtain, where we introduced the three scalar parameters tr ( Eqs. (2.31) are 3 coupled scalar ODEs, which determine the changes in the vibrational widths of atoms in phase space ( and ) and the correlation between the displacement and momentum vibrations of the th atom (see Fig. 2). We note that (2.31) are identical to the equations used by Ma et al. (1993), who used the formulation as an optimization procedure for a system of atoms at a uniform constant temperature.
The physical role of momentum-displacement correlation becomes evident upon applying a time-reversal transformation ↦ − to (2.30) and (2.31), which results in the transformations ( , , . Since the correlation signifies the momentum of the th atom in phase space, governing the time evolution of the thermal coordinate , will be referred to as the thermal momentum hereafter. The dynamics and thermodynamics of crystals modeled via the independent GPP approximation can be summarized via eliminating the thermal and dynamical momenta, yielding for every atom = 1, … , Time reversibility of (2.33) highlights that these evolution equations do not capture nonequilibrium irreversible thermal transport at all scales due to the simplification of (2.22) to (2.29) under the interatomic independence assumption (2.25). As discussed above, this assumption is necessary to keep the number of unknowns tractable for large ensembles.
The reversible entropy fluctuations can be obtained by substituting the independent GPP assumption into (2.24), which gives where the local entropy fluctuations of the th atom are given by ) . (2.35) The above equation shows that the local entropy fluctuatioṅis proportional to the thermal momentum . We will return to relation (2.35) when discussing specific types of interatomic potentials in subsequent sections.
Since the interatomic independence assumption results in an incorrect calculation of the interatomic heat flux in (2.22a), one needs to model irreversible thermal transport, e.g., using the linear kinetic potential framework used by Venturini et al. (2014), as discussed in the following.

Dynamics and quasistatics of independent GPPs
We proceed to analyze the dynamic behavior of the GPP Eqs. (2.33) (under the interatomic independence 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 equations to both harmonic and anharmonic potentials for the atomic interactions within the crystal. Although the GPP-based ansatz for the probability density consists of only quadratic functions in (resembling harmonic interactions), the force ( ) is derived from the potential ( ) and can be anharmonic in general. Such an approximation is classically known as the quasi-harmonic approximation because, at equilibrium, the distribution resembles that of a canonical ensemble of harmonic oscillators, even though the interatomic forces are anharmonic. In the following, we show that, for a quasistatic change in the thermodynamic state (i.e., driving the mean dynamical and thermal momenta and , respectively, to zero), the information about the evolution of is lost (cf. (2.33)). Correspondingly, we may assume a specific nature of the thermodynamic process of interest (e.g., isothermal, isentropic, or isobaric) to determine the change in of the GPP atoms during the quasistatic change in the thermodynamic state. Alternatively, the decay of correlation ( ) can be modeled empirically to obtain the change in , if the nature of the thermodynamic process is unknown.

Harmonic approximation
As a simplified case that admits analytical treatment, we first consider a harmonic approximation of the interaction potential ( ), writing where ∈ R 3×3 is the local harmonic dynamical matrix, 0 is the equilibrium potential between the atoms approximated as independent GPPs, and  ( ) represents the neighborhood of the th atom. Inserting potential (3.1) into the evolution Eqs. (2.33) leads to where tr( ) serves an effective force constant, which depends on the number of immediate neighbors of the th atom. Eqs. (3.3) show that the mean mechanical displacement ⟨ ⟩ of the th atom is decoupled from its thermodynamic displacements and for a harmonic potential field between the atoms. The resulting decoupled thermodynamic equations exhibit the following independent eigenvectors ( 0 , + , − ) and corresponding eigenvalues ( 0 , + , − ): so that the general homogeneous solution is composed of constant and oscillatory components. Coefficients ( 0 , + , − ) are determined by the initial thermodynamic state of each atom. The constant component 0 corresponds to = 0 and ∕ = 2 tr( ). To interpret the three terms within this solution, let us formulate the excess internal energy, which for the harmonic approximation (3.1) becomes .
( 3.7) By insertion into (3.7), it becomes apparent that the constant component 0 with = 2 tr( ) has equal average kinetic and potential energies. This equipartition of energy implies that 0 corresponds to the thermodynamic equilibrium. Consequently, the components ± correspond to oscillations about the equilibrium state 0 with frequencies ± = 2 √ 2 tr( )∕ . Therefore, numerical time integration of the GPP equations unfortunately incurs a similar computational cost as a standard MD simulation of the same system (although an equilibrium solution can be identified, oscillations about the equilibrium occur at frequencies on the same order as atomic vibrations). Thus, even though mean motion and statistical information have been separated, the interatomic independence assumption within the GPP framework prevents significant temporal upscaling. Further note that, due to the decoupling of the thermodynamic Eqs. (3.4) from the dynamic equation of motion (3.2), a harmonic GPP lattice exhibits no thermomechanical coupling and hence displays no heating nor cooling under external stress (i.e., it does not expand due to local heating and vice-versa).
Finally, by substituting the harmonic potential (3.1) into (2.35), we obtain the reversible fluctuations in entropy of atom for a system of atoms in a harmonic field: (3.8)

Anharmonic thermomechanical effects
As the above harmonic approximation of the interatomic potential renders the GPP crystal thermomechanically decoupled, we next study the effects of anharmonicity in the potential. As the simplest possible extension of the harmonic potential (3.1), we consider where ∈ R 3×3 is the local dynamical matrix, and denotes a constant anharmonic third-order tensor. With this potential, the evolution Eqs. (2.33) become and where are the components of , (⋅) denotes the th component of vector (⋅), and represents Kronecker's delta (and we use Einstein's summation convention, implying summation over , , ). Note that the second term in (3.10) couples the mechanical perturbations from the thermodynamic perturbations of atom and its neighbors. Moreover, since Eqs. (3.11) contain products of the thermodynamic variables and with the mean mechanical displacements ⟨ ⟩, the anharmonic potential leads to thermomechanical coupling in the GPP evolution equations.
Due to the apparent harmonic nature of most standard interatomic potentials, the time scale of the GPP Eqs. (3.10) and (3.11), when being applied to common potentials, is comparable to the time scale of atomic vibrations, since the system exhibits eigenfrequencies of 2 √ 2 tr( )∕ for a pure harmonic potential (cf. (3.5)). Consequently, numerical time integration of the GPP equations -in the anharmonic as in the harmonic potential case -incurs computational costs comparable to MD, so that (as discussed previously for the harmonic case) the interatomic independence assumption prevents temporal upscaling.

Quasistatics and thermal equation of state
Using the insight gained from the time evolution Eqs. (3.3) and (3.11), we proceed to study the GPP equations within the quasistatic approximation. The latter yields a system of coupled nonlinear equations, whose solution yields the thermodynamic equilibrium state of the crystal with atoms modeled using the GPP ansatz. In the quasistatic approximation, the GPP Eqs. (2.33) with mean mechanical momentum = 0 and thermal momentum = 0 reduce to the following steady-state equations: which are to be solved for the equilibrium parameters ( , , ) for each atom = 1, … , . Substitution of the quasistatic limits → 0, → 0 in (2.33) yields that, at thermomechanical equilibrium, the solution of Eqs. (3.12) corresponds to the equilibrium mean displacements and the equilibrium displacement-variance of all atoms. Therefore, (3.12) provides only two equations for the three unknowns ( , , ) per atom. Analogous to the loss of information about the evolution of the mean mechanical momentum ( ), the quasistatic limit only states that, at final thermodynamic equilibrium, the thermal momentum ( ) has decayed to 0, and and are related by Eq. (3.12b). (Note that substitution of = 0 in (2.33) yields trivially d ∕ d = 0 at quasistatic thermomechanical equilibrium.) To determine the evolution of ( ) during the thermomechanical relaxation of the system towards equilibrium, a model ( ) is required to complement (3.12). Hence, the quasistatic approximation results in the loss of information about the thermodynamics of the process through which the system is brought to the thermomechanical equilibrium. In the following, we consider different such thermodynamic processes and derive the respective process constraint, which -together with (3.12) -yields the result of a thermomechanical relaxation.
To physically approximate the nature of a thermodynamic process, we assume that the ergodic hypothesis holds for quasistatic processes, i.e., where is a sufficiently large time interval, over which the evolution of the system is assumed quasistatic. In the ergodic limit, the momentum variance becomes where is the local temperature of the th atom. As a consequence, to model an isothermal relaxation, we must keep and -by (3.14) -constant for each atom, so that is known for all atoms and the quasistatic Eqs. (3.12) can then be solved for ( , ). By contrast, isentropic equilibrium parameters ( , , ) can be obtained by keeping fixed for each atom during the relaxation. From (2.23) we know that may be interpreted as dimensionless momentum-variance and displacement-variance entropies, respectively. In the following, it will be convenient to use , and , as the mean free parameters instead of and .

Table 1
Summary of the thermodynamic constraints on for the different assumptions about the thermodynamic process under which the system is brought to equilibrium.

Process
Isentropic Isothermal Isobaric Analogously, isobaric conditions can be derived from the system-averaged Cauchy stress tensor (Admal and Tadmor, 2010) where  is the volume of the system. The average hydrostatic pressure of the system is .
( 3.17) Hence, setting = const. in (3.17) is the isobaric constraint, subject to which Eqs. (3.12) become serves as the momentum variance corrected for pressure , yields the equilibrium parameters ( , , , , ) for a given externally applied pressure . Note that for a system of non-interacting atoms, the quasistatic equations subjected to an external pressure reduce to the ideal gas equation  = . Consequently, (3.18) shows that the quasistatic approximation yields the thermal equation of state of the system accounting for the interatomic potential, which enables the thermomechanical coupling within the crystal. The three constraints on for isentropic, isothermal, and isobaric processes -based on the above discussion -are summarized in Table 1.

Helmholtz free energy minimization
The solution ( , , ) of the local equilibrium relations (3.12) may be re-interpreted as a minimizer of the Helmholtz free energy  (note that ( , , ) denotes the whole set of parameters of all atoms constituting the system). The Helmholtz free energy  is defined as with the internal energy of the system being . (3.20) The definition (3.19) implies the local thermodynamic equilibrium definition which can be verified using (3.15) and (3.20). In addition, minimization of  with respect to the parameter sets and , subject to any of the thermodynamic constraints in Table 1 for updating , yields Eqs. (3.12), i.e., A detailed derivation of (3.22b) is provided in Appendix C. Eqs. (3.12) (subject to a suitable thermodynamic constraint) are identical to the max-ent-based formulations of Kulkarni et al. (2008) and Venturini et al. (2014). Kulkarni et al. (2008) developed the max-ent formulation by enforcing constraints on variances of momenta and displacements of the atoms to obtain an ansatz for ( ), which is a special case of the GPP ansatz with no correlations (see Section 2.3). Venturini et al. (2014) developed the max-ent formulation by generalizing the grandcanonical ensemble to nonequilibrium situations and allowing non-uniform thermodynamic properties among atoms. However, for computational feasibility, Venturini et al. (2014) invoked a trial-Hamiltonian procedure to justify the Gaussian formulation of the distribution function (for single-species cases), thus rendering their ansatz also a special case of our GPP ansatz with no correlations. In the above sections, we have shown that such special case arises only as a result of the quasistatic approximation, which enforces vanishing mechanical and thermal momenta. Consequently, within the quasistatic assumption, the GPP based local-equilibrium Eqs. (3.12) are identical to those of Kulkarni et al. (2008) and Venturini et al. (2014). Moreover, subject to the isothermal constraints, GPP quasistatic Eqs.  ). ( ) Free-standing nano-cube composed of 12 × 12 × 12 FCC unit cells of pure single-crystalline copper, illustrating the spatial variation of due to a varying number of neighbors for atoms in the bulk vs. near the surface, at = 1000 K modeled using the potential of Johnson (1988). Thermal expansion is measured using the volume of the inner domain of the crystal (outlined in red). Atoms inside are marked in green. Fig. 4. Thermal expansion of a free-standing nano-cube (12 × 12 × 12 atomic unit cells) and of an infinite crystal of pure single-crystalline Cu, modeled using the Extended Finnis-Sinclair potential of Dai et al. (2006) and the exponentially decaying potential of Johnson (1988). ( ) Comparison of computed volumetric changes with experimental data from Nix and MacNair (1941) and MD results ( ref = ( = 273 K)), where ( ) is the volume of the inner subdomain of the nano-cube in Fig. 3( ) for the finite-cube calculation and 3 ( ) for the infinite-crystal calculation, at temperature . ( ) Variation of the displacement-variance entropy with temperature for the nano-cube calculation. As the temperature increases, the vibrational kinetic energy increases, resulting in an increase of at thermal equilibrium due to the equipartition of energy. Due to different numbers of interacting atomic neighbors, atoms on the corners, edges, faces, and in the bulk exhibit different values of (see Fig. 3). Note that the value of depends on the interatomic potential (shown results are for potential of Johnson (1988)).

Isothermal validation: thermal expansion of a Cu single-crystal
Within the scope of the present study, we aim to validate the computational implementation of the local-equilibrium Eqs. (3.12) within a finite-temperature, updated-Lagrangian 3D QC framework, which will be detailed and used in Section 4 to perform coarsegrained nonequilibrium thermomechanical simulations. For numerical validation of the isothermal quasistatic setting, we here compute the thermal expansion coefficient of a single-crystal of pure copper (Cu). We use an in-house developed updated-Lagrangian QC library, which is based on a simplicial mesh of tetrahedral elements in 3D and which adopts the seamless coarsening strategy of , so that all regions in a simulation (including the fully-resolved domain) are meshed. For computing phase-space averages, we use third-order multivariable Gaussian quadrature (Stroud, 1971;Kulkarni, 2007). at temperature undergoing the deformation measured by the strain tensor ( ) . The central atom (red) interacts with the full neighborhood and is relaxed isentropically, keeping the rest of the atoms mechanically fixed. Its phase-averaged potential ⟨ ( )⟩ inf models the average phase-averaged potential of an infinite crystal. ( ) Domain of 18 × 18 × 18 FCC unit cells of Cu. The domain is initially relaxed isothermally at the temperature at which the elastic constants are to be obtained. Afterwards, all the atoms are displaced according to the affine deformation (3.23). Blue atoms in the outer-layer are held mechanically fixed. Red atoms, initially bounded within the red-outlined box of size equivalent to 6 × 6 × 6 FCC unit-cells, are relaxed mechanically. All the atoms are kept thermally fixed since the elastic constants correspond to the curvature of the free energy against mechanical deformation only (Venturini, 2011). ⟨ ( )⟩ sub is the phase-averaged potential of the red atoms and ⟨ ( )⟩ box is the phase-averaged potential of all atoms in the domain, both averaged over the respective number of atoms. As ( ) increases, both phase-averaged potentials vary. Curvatures of both ⟨ ( )⟩ sub and ⟨ ( )⟩ box with respect to ( ) at the respective energy minima are used for computing the elastic constants. Fig. 4 illustrates solutions of the isothermal local-equilibrium relations for varying fixed uniform temperature = for all atoms within a perfect Cu single-crystal, modeled using the EAM potentials of Dai et al. (2006) and Johnson (1988). To simulate the expansion of an infinite crystal, we determine the lattice parameter and displacement-variance entropy which minimize the Helmholtz free energy  of an atom under the influence of its full centrosymmetric neighborhood in a periodic lattice (cf. Fig. 3( )). Kulkarni (2007) used the same infinite-crystal formulation, in which the mechanical forces are evaluated using the lattice spacing as the independent variable (see also the appendix of Tembhekar (2018)). To demonstrate the effect of free surfaces, we repeat the same calculation for a free-standing Cu nano-cube containing 12 × 12 × 12 FCC unit cells, subject to the local-equilibrium Eqs. (3.12), isothermal constraint, and free boundaries at various temperature values in our QC solver at full resolution (cf. Fig. 3( )). Results are also included in Fig. 4(a). As temperature increases, the lattice spacing increases. In turn, as the spacing between atoms increases, the displacement-variance entropy also increases (cf. Fig. 4( )). Further evident from that graph, is not uniform within the crystal, owing to varying numbers of interacting neighboring atoms on the corners, edges, and faces as well as within the bulk of the crystal. Sufficiently far away from the free boundaries, each atom is embedded in an approximately perfect infinite crystal and hence characterized by uniform values of the displacement entropy. To minimize the effects of free boundaries, we compute the deviation in volume of a bounding box enveloping a 2 × 2 × 2 volume at the center of the crystal with temperature (outlined in red in Fig. 4( ) at the center of the nano-cube) to compute the thermal expansion coefficient in the free-standing cube case. To validate our results, we further report comparable thermal expansion data obtained from the MD code LAMMPS (Plimpton, 1995), using an NPT ensemble with periodic boundaries (labeled 'MD' in Fig. 4( )).
The comparison shows that the GPP-based formulation (which in this isothermal setting is equivalent to the max-ent formulation) yields reliable data for the chosen potentials with marginal errors up to half the melting temperature. At higher temperature, the GPP results under-predict the thermal expansion -as can be expected -with errors being significantly smaller for the nearestneighbor-based potential of Johnson (1988). We note that the GPP-based results depend on the chosen third-order multivariable Gaussian quadrature for computing phase-space averages. Even though higher-order quadrature has been shown to yield improved accuracy (Tembhekar, 2018) (at considerably increased computational cost), a detailed case study of the impact of the order of quadrature is beyond the present scope.

Isentropic validation: elastic constants of a Cu single-crystal
Going beyond uniform thermal expansion, we proceed to benchmark the GPP-based formulation with isentropic process conditions by calculating the uniaxial and shear components of the linear elastic stiffness tensor ( 11 and 44 , respectively, in Voigt-Kelvin notation, see Reddy, 2007), and the bulk modulus -which we collectively refer to as elastic constants in the following. As a sample system we again choose single-crystalline pure Cu.
Calculation of the aforementioned elastic constants ( 11 , 44 , and ) was previously used as a benchmark by Venturini (2011) for a crystal with impurities, and by Tembhekar (2018), who all used the max-ent formulation (without specifying the isentropic process condition). To compute the elastic constants, we follow a procedure similar to the one by . We consider a domain of 18 × 18 × 18 FCC unit-cells of Cu, subject to the local equilibrium Eqs. (3.12) and the isentropic process P. Gupta et al.  Dai et al. (2006), compared against the experimental data from Chang and Himmel (1966) and Overton and Gaffney (1955).
constraint (see Fig. 5( )), so as to model an isentropic deformation process around the undeformed ground state at a given (initial) temperature . We use the isentropic constraint for relaxing the deformed crystal, because the elastic constants are measured in an adiabatic setting (Overton and Gaffney, 1955). Initially, the domain is relaxed isothermally with free boundaries. After the initial relaxation, atoms are displaced according to, where ( ) = is the engineering strain measure at the th deformation step, is the strain increment, and is the base deformation matrix. For 11 , 44 , and , it is defined as, respectively, (3.24) To model an infinite crystal without size effects imposed by free surfaces, we also consider a domain of 2 × 2 × 2 FCC unit cells of Cu, initialized using the lattice parameter and displacement-variance entropy which minimize the Helmholtz free energy at the given temperature (see Fig. 3( )). The domain is deformed using (3.23), isentropically relaxing the central atom while keeping its neighborhood mechanically fixed by the external deformation (Fig. 5( )). As the atoms are displaced, the free energy of each atom changes due to changing potential energy. For the infinite-crystal model, the atom at the center interacts with its whole neighborhood, which deforms as per (3.23) and exhibits a change in its phase-averaged potential ⟨ ( )⟩ inf . The affine deformation mimics the straining of a perfect infinite crystal. The phase-averaged potential of the central atom, ⟨ ( )⟩ inf , is measured for all increments to compute the elastic constants. When modeling a finite-sized, free-standing crystal, relaxation is performed while holding the atoms in a layer touching the free boundaries of the domain mechanically fixed (blue atoms in Fig. 5( )). The phaseaveraged potential of the whole domain, ⟨ ( )⟩ box , and that of a sub-domain within the bulk of the domain, ⟨ ( )⟩ sub (also averaged over the number of respective type of atoms), are measured for all . We use a 4th degree polynomial fit through the phase-averaged potentials to obtain continuous functions ⟨ ( )⟩ box ( ), ⟨ ( )⟩ sub ( ), and ⟨ ( )⟩ inf ( ), and to compute the respective curvatures at ( ) = , for which ⟨ ( )⟩ sub , ⟨ ( )⟩ box , and ⟨ ( )⟩ inf are minimal. The elastic constants at temperature are obtained from the curvatures via (Venturini, 2011) where ( , ) is the atomic volume at and temperature . We note that, by definition, the change in free energy  ( , ) with mechanical strain measure is used for computing the elastic constants (Venturini, 2011;Phillips and Rob, 2001). Since the change in  ( , ) is identical to the change in the phase-averaged potential ⟨ ( )⟩ for only mechanical deformations, we have used ⟨ ( )⟩ in (3.25). Fig. 6 shows the change of ⟨ ( )⟩ box , ⟨ ( )⟩ sub , and ⟨ ( )⟩ inf with for various temperatures and the elastic constants 11 , 44 , and , obtained using the GPP-based local-equilibrium Eqs. (3.12), and experimental data (Overton and Gaffney, 1955;Chang and Himmel, 1966). Results were computed for decreasing values of to ensure convergence. The reported results correspond to = 0.0005. The computed elastic constants display similar thermal softening as observed in experiments, yet they exhibit an offset at all temperatures since the EAM potential of Dai et al. (2006) was calibrated using the elastic constants at room temperature, which are treated as those at 0 K in the present formulation (see Table 4 in Dai et al. (2006) and Table 3.1 in Kittel (1976)). At 0 K, the values of 11 and 44 obtained from the infinite-crystal model ( 11 = 168.41 GPa, 44 = 75.41 GPa) match exactly the values reported by Dai et al. (2006). Those obtained from the finite-crystal simulation setup ( 11 = 164.34 GPa, 44 = 78.01 GPa using ⟨ ( )⟩ sub and 11 = 154.20 GPa, 44 = 73.12 GPa using ⟨ ( )⟩ box ) deviate from those reported by Dai et al. (2006) due to the size-effects posed by the free surfaces, as expected. The overall behavior captures the thermal softening of pure Cu well.

Linear Onsager kinetics
We reiterate that Eqs. (3.22) are identical to the max-ent framework, as derived and utilized previously (Kulkarni et al., 2008;Ariza et al., 2012;Venturini et al., 2014;Ponga et al., 2015;Tembhekar, 2018). However, the max-ent framework is based on a different motivation and does not rely on the physical insight gained in Sections 2.3, 3.1.1, and 3.1.2 about the dynamics of atoms at long and short time intervals. Moreover, the GPP framework clearly highlights the information loss as a result of the quasistatic approximation (otherwise hidden in the max-ent framework), which also enables the modeling of thermomechanical deformation of crystals under various thermodynamical conditions (e.g., isentropic, isobaric, or isothermal processes). Furthermore, recall that in Section 2.3 we showed that the interatomic independence assumption precludes the framework from capturing any changes in local temperature due to unequal temperatures of neighboring atoms. Such an assumption is also made in the max-ent framework. In reality, a non-uniform temperature distribution arises readily in non-uniformly deformed crystals (differences in local atomic spacings lead to different local temperatures). As both max-ent and the present GPP formulation cannot model heat flux as a consequence of such local temperature differences, we require explicit thermal transport models to capture heat transport. To this end, we here adopt the linear Onsager kinetics model introduced by Venturini et al. (2014), based on a quadratic dissipation potential -as discussed in the following.
The total rate of change of entropy of an atom can be decomposed into a reversible and an irreversible change, i.e., where ,rev is the reversible heat addition and ,irrev is the irreversible change signifying the net influx of heat into the atom due to a non-uniform temperature distribution. In a dynamic system (see Sections 3.1.1 and 3.1.2 ), the reversible change is composed of fluctuations about the equilibrium state due to the local harmonic nature of the interatomic potential, and it is proportional to the thermal momentum . Within the quasistatic approximation, the information about the evolution of ( ) as the system relaxes towards the equilibrium is lost since → 0 is imposed, thus rendering the reversible changes in entropy unknown. Therefore, such a reversible change is imposed implicitly by the thermodynamic constraints (see Table 1). For an isolated system of atoms, the reversible heat exchange vanishes ( rev = 0), and the system relaxes to equilibrium adiabatically, which we here term free relaxation. Note that this free relaxation is not isentropic, since the temperature can be nonuniform as a result of an imposed non-uniform deformation of the crystal, resulting in irreversible thermal transport that increases the entropy. We model such irreversible change ,irrev by adopting the kinetic formulation introduced by Venturini et al. (2014): where is the local, pairwise heat flux, driven by a local pairwise discrete temperature gradient through the kinetic potential . Using the dissipation inequality, Venturini et al. (2014) formulated the discrete temperature gradient as = 1 − 1 . (3.28) Within the linear assumption, the kinetic potential is modeled as, where denotes a pairwise heat transport coefficient (which is treated as an empirical constant in this work), and represents the pairwise average temperature. Eqs. (3.27) and (3.29) yield the entropy rate kinetic equation for a freely relaxing system of atoms as (3.30) P. Gupta et al. To use the model, one needs to calibrate the heat transport coefficients , for which experimentally measured values of for a given arrangement of atoms can be used. To this end, Venturini et al. (2014) derived the following discrete-to-continuum relation between and the thermal conductivity tensor at the location of atom : where is the volume of the atomic unit cell in the crystal. The obtained values of may be interpreted to capture the interatomic heat current and regarded as the intrinsic property of the material. In the model of Ponga and Sun (2018), an empirical parameter (exchange rate) was fitted against theoretically or experimentally characterized thermal conductivity values.
It is important to point out that Eq. (3.31) was derived from first-order approximations of temperature differences between interacting atoms, using Taylor expansions (Venturini et al., 2014). Depending on the thermal boundary conditions and the size of the domain, temperature differences between interacting atoms may not be negligible. For significant temperature differences, (3.31) is inaccurate, hence, we here identify the values of by simulating a differentially heated Cu nanowire of square cross-section up to steady state. Specifically, we fit the kinetic coefficient 0 for Cu, using the thermal conductivity measurements for Cu nanowires obtained by Mehta et al. (2015). To this end, we consider a Cu nanowire of square cross-section of size 145 × 43 × 43Å 3 , with the central axis of the nanowire oriented along the -direction and all edges aligned with the ⟨100⟩-directions (see Fig. 7 ). Atoms in the region < are thermostated at a temperature of 360 K, while the atoms in the region > are thermostated at a temperature of 240 K. All the boundaries are considered as free surfaces and the system is relaxed isentropically, followed by the diffusive step for thermal transport (as explained in Algorithm 1 below). As the system evolves according to (3.30) (using full atomistic resolution everywhere in the nanowire), a uniform macroscopic heat flux is established between < < (see Fig. 7 ). Dividing the nanowire into atomic planes marked by , the macroscopic flux across the plane at = is given by where and ∕ are temperature and entropy generation at plane . From (3.32) the approximate thermal conductivity is obtained as where is the cross-sectional area, and is the length across which the temperature difference is maintained. To find 0 , we do not use Eq. (3.31) since it is valid only for small temperature differences. Instead, we start from a reference guess of 0 = 0.1 nW/K, compute the macroscopic flux , and from it the thermal conductivity via (3.33) and compare it to the experimental value. To obtain an approximate value of 0 for our simulations, we considered = 100 W∕(m ⋅ K) as determined experimentally by Mehta et al. (2015) (cf. Figure 3 in Mehta et al., 2015). The initial guess of 0 is updated using the secant method, until the conductivity P. Gupta et al. value of = 100 W∕(m⋅K) is reached. The obtained value of 0 ≈ 15.92 eV∕(ns⋅K) = 2.55 nW/K is used in all subsequent simulations. We note that the numerical values, however, do not affect the physical modeling framework and are only representative values to be used in subsequent simulations.
In reality, large deformations cause defects which modify the phonon and electron scattering properties of a crystal, due to which 0 may vary with deformation (the same applies to considerably affine straining of a crystal). However, such effects are beyond the scope of the current work. As shown in Fig. 7 , the temperature profile is linear and the macroscopic thermal flux defined by (3.33) is constant at steady state, thus highlighting that the discrete model yields a behavior similar to Fourier's heat conduction law, where material between two isothermal boundaries exhibits a linear temperature distribution and constant heat flux, given that the conductivity is uniform. The model can be extended to deformation-dependent values of 0 -yet this requires a detailed (and presently unavailable) understand of scattering phenomena in strained crystals. Therefore, we here assume 0 = const. as a first-order approximation.
Eqs. (3.12) combined with (3.30) complete the nonequilibrium thermomechanical atomistic model. Every atom is assumed to be in thermal equilibrium at some local temperature , and mechanical and thermal interactions of the atoms with different spacings and different temperatures are modeled using the interatomic potential ( ), the local equation of state (e.g. (3.18)), and the entropy kinetic equation (3.30). For a general system, the assumed thermodynamic process constraints yield the reversible heat exchange (see Table 1). For instance, in an isothermal constraint, , remains constant and , changes with mechanical deformation to satisfy the local equation of state (e.g. (3.18)) at equilibrium, thus changing the net entropy (see (3.15)).
While the nature of the assumed thermodynamic process via which the system relaxes depends on the macroscopic and microscopic boundary conditions, the process is generally assumed quasistatic with respect to the fine-scale vibrations of each atom. However, the thermal transport equation (3.30) introduces a time scale to the problem governing the irreversible evolution of entropy. If we consider a two-atom toy model, consisting of only two atoms with temperatures and (Fig. 8) where we have assumed equal coefficients = 0 for both atoms. Let us further assume the thermomechanical relaxation of ( , ) takes place through an isentropic process following the quasistatic Eqs. which shows that, when using a simple first-order explicit-Euler finite-difference scheme for time integration, the time step is restricted by for numerical stability. Venturini et al. (2014) used 0 = 0.09 nW/K for the bulk region of a silicon nanowire, which yields a maximum allowed time step of 0.23 ps. With higher-order explicit schemes, a larger maximum time step may be obtained, but the maximum allowed time step remains at a few picoseconds. Such restriction arises because the bulk thermal conductivity is used to fit the atomistic parameter , which reduces the length scale as well as the time scale. Such a restriction also highlights the P. Gupta et al. coupling of length and time scales (Evans and Morriss, 2007). For larger temperature differences and larger systems, the nonlinear equation (3.35) yields the following stability limit on time step (see Appendix D for the derivation) (3.38) which is stricter than the linear stability limit in (3.37). In summary, even though we simulate quasistatic behavior at, in principle, considerably larger time scales by the GPP-based formulation, the use of the linear thermal transport model severely restricts the usable time step size in simulations.

Algorithmic implementation
Before applying QC coarse-graining to the GPP-based formulation, let us summarize the proposed thermomechanical transport model for simulating the quasistatic deformation of crystals composed of GPP atoms (see Algorithm 1 for details). Our numerical scheme decomposes each load step as follows: • Step 1: Given the equilibrium parameters for each atom from the (previous) th load step, an external stress/strain is applied to the system at load step + 1.
• Step 2: Quasistatic relaxation, solving (3.12) subject to one of the constraints in Table 1 to obtain the intermediate state ) for each atom . • Step 3: Staggered time stepping that alternates between irreversible updates of the total entropy and quasistatic reversible relaxation steps of all variables, until convergence is achieved. Specifically, the total entropy is updated irreversibly from ( * ) over a time interval , using Eq. (3.30) and explicit forward-Euler updates, and reversibly using the assumed thermodynamic process constraints during the subsequent thermomechanical relaxation (see the two contributions in (3.26)). Time step depends on the external stress/strain rate applied. By definition, the thermomechanical relaxation is assumed quasistatic, hence only slow rates with respect to atomistic vibrations can be modeled. However, the irreversible transport imposes a time-scale restriction via (3.38). Consequently, time integration via suitable time steps must be continued for steps, such that ∑  Table 1.
This completes a single staggered time step. Note that the update ( * * ), +1 → ( * ), +1 corresponds to the reversible entropy update to satisfy (3.15) during the thermomechanical relaxation and depends on the assumed thermodynamic process constraint. In a full quasistatic setting, the transport equation is driven towards a steady state witḣ( * ), = 0, which defines the convergence criterion and hence determines the total number of time steps. , followed by Step 1 until the final load step.
To complete Step 2, we use a combination of the Fast Inertial Relaxation Engine (FIRE) of Bitzek et al. (2006) and nonlinear generalized minimal residual (NGMRES) using the Portable, Extensible Toolkit for Scientific Computation (PETSc) (Brune et al., 2015). We employ a forward-Euler finite-difference scheme to update the entropy due to irreversible transport in Step 3. The steps are explained in detail as pseudocode in Algorithm 1.

Finite-temperature updated-Lagrangian quasicontinuum framework based on GPP atoms
Having established the GPP framework, we proceed to discuss the application of the thermomechanical and coupled thermal transport model (Eqs. (3.12) and (3.30), respectively) to an updated-Lagrangian QC formulation for coarse-grained simulations. Previous zero-and finite-temperature QC implementations have usually been based on a total-Lagrangian setting (Ariza et al., 2012;Tadmor et al., 2013;Knap and Ortiz, 2001;Ponga et al., 2015;Kulkarni et al., 2008), in which interpolations are defined and hence atomic neighborhoods computed in an initial reference configuration. Unfortunately, such total-Lagrangian calculations incur 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 initial mesh used for all QC interpolations increasingly loses its meaning and atoms that form local neighborhoods in the 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 deformation and atomic rearrangement. Moreover, we adopt the fully-nonlocal energybased QC formulation of , since an energy-based summation rule allows for a consistent definition of the coarse-grained internal energy and all thermodynamic potentials of the system.
Within the QC approximation, we replace the full atomic ensemble of GPP atoms (as described in Section 2.3) by a total of ℎ ≪ representative GPP atoms (repatoms for short), each having the thermomechanical transport parameters ( , , , , ,̇) as their degrees of freedom (see Fig. 9). The position of each and every atom in the coarse-grained crystal lattice is obtained by Algorithm 1: Single load step from to + 1 of the quasistatic thermomechanical transport model for irreversible deformation of crystals composed of GPP atoms. For a fully quasistatic transport simulation ( ) → ∞ and the thermal gradients are diffused till the RMS entropy generation rate is higher than some tolerance .
interpolation. For an atom at location ℎ in the reference configuration, the thermomechanical transport parameters in the current configuration are obtained by interpolation from the respective parameters of the repatoms: where the superscript ℎ indicates that a parameter is interpolated from the ℎ repatoms, and ( ℎ ) is the shape function/interpolant of repatom evaluated at ℎ . In our updated-Lagrangian setting, the respective previous load step serves as the reference configuration used to define the mesh and the above interpolation. In the following we use linear interpolation (i.e., constantstrain tetrahedra in 3D), while the method is sufficiently general to extend to other types of interpolants. Based on the interpolated parameters from (4.1), the free energy  ( , , ) of the crystal is replaced by the approximate free energy  ℎ of the QC crystal with where we assumed that the decomposition ( ) = ∑ ( ) of the interatomic potential holds. Furthermore, we allow the masses of all atoms to be different, denoting by the mass of atom . Eq. (4.2) defines the free energy of the system accounting for all atoms with their thermomechanical parameters evaluated using the ℎ repatoms in a slave-master fashion.
To reduce the computational cost stemming from the summation over all atoms in (4.2), sampling rules are introduced, which approximate the full sum by a weighted sum over ≪ carefully selected sampling atoms (Eidel and Stukowski, 2009;Iyer and Gavini, 2011;Tembhekar et al., 2017): where is the sampling weight of the th sampling atom. Specifically, we adopt the first-order optimal summation rule of , in which all nodes and the centroid of each simplex (tetrahedron in 3D) are assigned as sampling atoms.  computed the sampling atom weights using the geometrical division of the simplices by planes at a distance from the nodes (Fig. 9( )) and adding the corresponding nodal volume to the respective sampling atom weight, while the rest of the simplex volume was assigned to the Cauchy-Born-type sampling atom at the centroid. Here, we point out that a simpler weight calculation is possible by considering the spherical triangle generated by the intersection of simplex with a ball of radius 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 ( + + − ) 2 . Hence, the approximate volume of the enclosed region is and = ∑ are the sampling atom weights at nodes, where is the density of simplex (expressed as the number of atoms per unit volume). For the centroid sampling atoms, the remaining volume times is assigned as its sampling weight . Since the deformation is affine within each element owing to the chosen linear interpolation, sampling atom weights in coarse-grained regions change negligibly in a typical simulation and are therefore kept constant throughout our simulations. In the following we will also need a separate set of repatom weightŝ, which we calculate by lumping the sampling atom weights to the repatoms: each repatom receives the weight of itself (each repatom is a sampling atom) plus one quarter of the Cauchy-Born-type centroidal sampling atoms within all adjacent elements : (4.5) Given the sampling atom weights , minimization of the approximate free energy (4.3) with respect to degrees of freedom ( , , ) of the th repatom yields the local mechanical equilibrium conditions and the corresponding thermal equilibrium conditions Substituting the interpolation from (4.1) into (4.6) and (4.7) yields (for repatoms = 1, … , ℎ ) Note that these are similar to the equilibrium conditions derived by Tembhekar (2018)   system, here providing the equation of state of the coarse-grained quasicontinuum. Solving Eqs. (4.6) and (4.7) subject to one of the constraints from Table 1 depending upon the assumption of the thermodynamic process yields variables ( , , , , ) for all repatoms, thus yielding the thermodynamically reversible solution for the deformation of the system. To introduce seamless coarse-graining of the linear Onsager kinetic model for irreversible thermal conduction governed by Eq. (3.30), we solve for the entropy rateṡof all repatoms and evolve the entropy in time for each repatom. To this end, we notice that the first term in Eq. (4.8b) represents a thermal force due to the thermal kinetic energy of the system. Since that thermal force in our energy-based setting follows a force-based summation rule (Knap and Ortiz, 2001), the entropy rate calculation of the repatom simplifies tô wherêis the repatom weight. Note that the kinetic potential in (3.29) can, in principle, also be coarse-grained analogously to the free energy in (4.3). However, the resulting calculation of the entropy rates is computationally costly since it involves summation over all repatoms for each sampling atom calculation, which is why this approach is not pursued here. Eqs. (4.8a) and (4.8b) combined with the thermodynamic constraints in Table 1 and the coarse-grained thermal transport model in (4.9) yield the solution of a generic thermomechanical deformation problem subject to a non-uniform temperature distribution and loading and boundary conditions, as illustrated in Fig. 1. Convergence of the force-based summation rule was analyzed by Knap and Ortiz (2001) for repatom forces. Hence, Eq. (4.9) converges to (3.30) as the coarse-grained mesh is refined down to atomistic resolution (weightŝ and approach unity, and the dependence on all sampling atoms excluding the repatom vanishes). Consequently, even in the coarse-grained regions, Eq. (4.9) approximates the atomistic thermal transport model governed by (3.30). As shown in Section 3.6, the linear Onsager kinetic model approximates Fourier's law-type thermal transport for a relaxed crystal, since the temperature reaches a linear distribution at steady state and the macroscopic flux reaches a constant value for differentially heated boundaries (neglecting any deformation dependence of the thermal conductivity). Since the deformation in large coarse-grained elements in a QC simulation is expected to be small, the coarse-grained equation (4.9) is expected to approximate Fourier's law-type thermal transport sufficiently well.

Updated-Lagrangian QC implementation
We implement the thermomechanical local equilibrium relations (4.8a) and (4.8b) combined with a thermodynamic 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 have significant 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, tracking the interatomic potential neighborhoods in the undeformed configuration incurs high computational costs. Alternatively, one could strictly separate between atomistic and coarse-grained regions (as in the local-nonlocal QC 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 tracking of, e.g., lattice defects with full resolution (Tembhekar et al., 2017). It also requires a-priori knowledge about where full resolution will be required during a simulation. As a remedy, we here deform the mesh with the moving repatoms and we take the deformed configuration from the previous load step as the reference configuration for each new load step, thus discarding the initial configuration and continuously updating the reference configuration.
For every element , we store the three initial edge vectors (i.e., three node-to-node vectors forming a right-handed system) in a matrix 0 , and the three Bravais lattice vectors indicating the initial atomic arrangement within the element in a matrix 0 . As the system is relaxed quasistatically under applied loads, all repatoms move to the deformed configuration (e.g., from load step = to = + 1), thus deforming the edge vectors of element from to +1 (and likewise the Bravais basis from to +1 , see Fig. 10). The incremental deformation gradient of element , from step to + 1, can hence be computed from the kinematic relation (4.10) which assumes an affine deformation within the element due to the chosen linear interpolation (see Fig. 10). As the element deforms, the lattice vectors also deform in an affine manner: Consequently, the integer matrix , which contains the numbers of lattice vector hops along the element edges, evaluated as = ( , ) −1 = const., (4.12) remains constant throughout deformation of a given element . Moreover, each element edge has a constant number vector, denoted by the rows of (see Fig. 10). That is, in the updated-Lagrangian setting, the number matrix remains constant during deformation. Such conservation of lattice vector hops along the element edges/faces is particularly useful for adaptive remeshing scenarios, where existing elements may need to be removed and new elements need to be reconnected, with or without changes to the number of lattice sites used for re-connections. The conservation of lattice vector hops can then be used for computing the Bravais lattice vectors local to new elements. The Bravais lattice vectors are used for calculating the neighborhoods of the nodal and centroid sampling atoms belonging to the large elements in the fully nonlocal QC formulation. The local lattice is generated within a threshold radius distance from the sampling atom using those lattice vectors. We adopt the adaptive neighborhood calculation strategy of Amelang (2016), which requires larger threshold radii compared to the interatomic potential cut-off chosen as where buff is a buffer distance used for triggering re-computations of neighborhoods, and cut is the interatomic potential cut-off. If the maximum relative displacement among the neighbors with respect to a sampling atom exceeds 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 are not required. Instead, the unique nodes of all elements within the threshold radius of the sampling atom are added as the neighbors. Consequently, even severely deforming meshes do not require element reconnection/remeshing as long as the deformation remains restricted within the atomistic region, since only nodes of the elements are required. Hence, we use meshes with large atomistic regions in the benchmark cases presented below, to restrict the analysis towards thermodynamics of the deformations. Such simulations do not require adaptive remeshing, the analysis of which is left for future studies.

Thermal effects on the shear activation of dislocations
As a benchmark example, we use the updated-Lagrangian QC method discussed above to analyze the effects of 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 constraints (the latter combined with the irreversible entropy transport based on linear Onsager kinetics). The adiabatic constraint here signifies that the simulation domain is thermally isolated from the surroundings and there is no heat exchange between the domain boundaries and the surroundings. For both cases, we generate a pair of dislocations (i.e., a dislocation  Fig. 12. Comparison of the isothermal nucleation of dislocation dipoles from a single-atom void, as obtained from fully resolved (left) and QC (right) simulations at varying temperatures in a Cu single-crystal modeled using the EFS potential (Dai et al., 2006). Shown is the final sheared state of ( ) a snapshot of the fully-resolved simulation and ( ) that of the QC simulation. Atoms are colored by the centrosymmetry parameter in arbitrary units. While the atoms in region A are kept fixed, atoms in region B are allowed to relax. (c,d) The shear stress vs. the engineering strain is plotted for ( ) the fully-resolved simulation and ( ) the QC simulation. The shear stress is evaluated as the net force in the [110]-direction on the atoms in region A per cross-sectional area in the (111) plane. Faces (112) are periodic. dipole) using the isotropic linear elastic displacement field solutions of edge dislocations with opposite Burgers' vectors (Nabarro, 1967), given by (4.14b) superposed linearly onto a 32×25×1.8 slab of pure single-crystalline Cu, consisting of 125,632 lattice sites, whose edges are oriented along the slip crystallographic directions. In (4.14b), ′ 1 and ′ 2 denote the coordinates along the [110]-and [111]-axes, respectively, with respect to the dislocation centers; 1 and 2 are the displacements along these axes, and = ± √ 2 is the magnitude of the Burgers' vector. Fig. 11( ) shows the thus-generated initial configuration of the edge dislocation dipole with Burgers' vectors = ± 2 [110]. separated by a horizontal distance of 80 Å. The centers of the dislocations are separated by a distance of 80 Å in the (horizontal) [110]-direction and a small distance of 1 × 10 −9 Å in the (vertical) [111]-direction. Imposing such displacement fields causes the slip planes of the two dislocations to separate in 2 and leaves a line of atoms on the 2 = 0-plane at the left end of the domain with zero displacements. These atoms are removed from the simulation domain, thus creating void at the edge of the domain. Displacements are restricted to the (112) plane, while the simulation domain is set up in 3D with periodic boundary conditions on opposite out-of-plane faces. After initial relaxation, the dislocations annihilate each other due to their interacting long-range elastic field. The result is a line defect in the form of a through-thickness (non-straight) vacancy column (in the following for simplicity referred to as a void), as shown in Fig. 11( ). This void was created in the initial configuration, when the non-displaced atoms were removed from the simulation domain, and it diffuses to the center of the domain during the initial relaxation. We note that this is a direct consequence of separating the slip planes of the two dislocations in our simulation. If the slip planes are identical, then the dislocations annihilate and form a perfect single-crystal. We deliberately create the void in this fashion, since it ensures that Fig. 13. Comparison of the critical shear stress cr and strain cr required to nucleate a dislocation dipole from the void as obtained from QC and from fully-resolved simulations at various temperatures. The critical strain cr is the external shear strain at which cr is achieved. only two dislocations of opposite orientation are nucleated when the domain is externally sheared. After relaxation, we load the simulation domain in simple shear (moving the top and bottom faces relative to each other), while computing the effective applied shear stress from the atomic forces. Periodic boundary conditions are imposed on ( 112 ) -surfaces, while the rest of the boundaries are included within region A, which is mechanically fixed during relaxation. At sufficient applied shear, the void will nucleate and emit a dislocation dipole, whose activation energy and behavior depend on temperature. For an assessment of the accuracy of the QC framework, we carry out both fully resolved (125,632 atoms) and QC simulations (52,246 repatoms) in isothermal and adiabatic settings. QC simulations are performed on a mesh generated by coarse-graining in the 2 -direction. All three lattice vectors are expanded by a factor of 4 in the coarse-grained region. The atomistic region extends fully in the 1 -and 3 -directions and up to ±51 Å in the 2 -direction. Coarse-graining is performed only in the 2 -direction to prevent the dislocations colliding with the atomistic and coarse-grained subdomain interface. We acknowledge that the QC setup is relatively simple and there is not yet a significant reduction in the total number of degrees of freedom nor does it involve automatic mesh refinement. Yet, this study presents a simple and instructive example highlighting the efficacy and accuracy of the GPP-based QC formulation introduced in previous sections.

Isothermal
In face-centered cubic (FCC) crystals, edge dislocations preferably glide on the close-packed crystallographic {111}-planes (Hull and Bacon, 2001). As the initial void is strained under shear deformation, dislocations nucleate from the void at a sufficient level of applied shear, propagating in opposite directions, as shown in Fig. 12. We apply a shear deformation to all repatoms in the slab such that, at the th load step, where indices 1 and 2 refer to the respective components of the mean position vector in the chosen coordinate system, and is the applied shear strain increment. As the strain is applied, repatoms in the inner region B (Fig. 12) are relaxed, while keeping those in the outer region A mechanically fixed to impose the average shear strain. Note that, due to small deformation in the atomic neighborhoods, the displacement-variance entropy of repatoms close to the interface between regions A and B changes and, hence, all repatoms in the domain are thermally relaxed assuming an isothermal relaxation (cf. Table 1). While the shear strain is P. Gupta et al. Fig. 15. Shear stress on the (111)-plane and dislocation separation distance vs. applied shear strain for isothermal and adiabatic deformations, obtained from both atomistic and QC simulations. Note that the differences between isothermal vs. adiabatic data are small, because the temperature increase is not significant in this example. Critical shear stress value deviations are within 6% (e.g., isothermal QC: 3.7914 GPa, isothermal atomistics: 3.6389 GPa) and critical strain values are within 12% (isothermal QC: 0.112, isothermal atomistics: 0.100). Fig. 16. Local variation of temperature as the slab with void is deformed adiabatically. In the initial stages, the temperature rises slowly due to the external deformation. At the critical shear strain cr = cr , dislocations nucleate from the void and move along the [110]-direction. The temperature rise of those atoms within the dislocations causes a rapid increase in the temperature of the slab, as may be expected from the heat generated by work hardening. The temperature of a few atoms within the dislocations at the intermediate stage = * cr before the irreversible transport exceeds the color bar range.
increased, the horizontal 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. Results are shown in Fig. 12( ) and ( ). Once the stress reaches a critical value, the stress drops as dislocations nucleate from the void and move to the ends of region B. We observe that the critical stress value decreases slowly with temperature (see Fig. 13). Moreover, the value of the critical stress obtained from a fully atomistic simulation and the quasicontinuum simulation are within about 6% of each other (see Fig. 13), both capturing the apparent thermal plastic softening typical for Cu (Shim et al., 2016). Atoms around the nucleated dislocations have severely deformed neighborhoods, which causes the local position entropy to vary within the dislocation cores (see Fig. 14). Miller and Tadmor (2009) studied a similar 2D scenario with a different crystal orientation, in which the dislocation 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 crystallographic plane. By contrast, in our analysis (which is based on a more realistic crystallographic setup since the slip planes of the dislocations are close-packed planes of the {111}-family), the critical shear stress is defined as the shear stress required to nucleate dislocations from the void-like defect.

Adiabatic
To simulate the quasistatic adiabatic nucleation of dislocations under shear, we repeat the above simulations, now with all repatoms in the domain being relaxed isentropically (cf. Table 1), followed by the thermal transport model according to the steps discussed in Section 3.6 and Algorithm 1. As noted above, the term adiabatic refers to a thermally insulated scenario, in which the domain boundaries do not allow for heat exchange with the surroundings. We further assume that the strain rate is significantly lower than the time scale imposed by the molecular thermal transport, thus imposing quasistatic conditions for the transport ( ( ) → ∞ in Algorithm 1, see Step 3 at the end of Section 3.6). The initial condition (again, prepared using the isotropic elastic displacement field solutions) is relaxed isothermally at 300 K, before the adiabatic deformation begins. We compare the adiabatic deformation with the isothermal deformation of the slab at 300 K. Fig. 15 shows the variation of the shear stress on the (111)-plane with external shear strain . Due to the mechanical work done by the external shear deformation, the temperature of the slab increases slightly, causing apparent softening compared to the isothermal deformation. Fig. 16 shows the spatial variation of temperature as the slab is P. Gupta et al. deformed adiabatically. Before the critical strain, heating caused by local deformation is negligible. As the dislocations are nucleated at the critical strain cr , the temperature of those atoms around the dislocations changes significantly, as shown by the intermediate stage ( * , * ) in Fig. 16. Due to quasistatic thermomechanical deformation, the temperature field is homogenized to within 1 K, even after dislocation nucleation from the void. Such close-to-isothermal deformation at slow strain rates was previously observed by Ponga et al. (2016) (who studied strain-rate effects on nano-void growth in magnesium) and by Ponga and Sun (2018) (who studied thermomechanical deformation of carbon nanotubes.) Further plastic deformation causes increased heating of the slab, particularly due to the restricted dislocation motion beyond the interfaces between regions A and B (Fig. 12). As noted above, the critical stress values obtained from QC simulations are within 6% of those obtained from the fully resolved simulations. Furthermore, the critical strain values are within 12%. Repatoms on the (110)-and (111)-surfaces, which includes the repatoms on vertices of large elements and the repatoms in the transition region between regions A and B, exhibit both mechanical and thermal spurious forces in 3D Tembhekar et al., 2017). These artifacts are expected to be the primary sources of error in the coarse-graining strategy adopted here. Mechanical spurious forces within the energy-based fully nonlocal QC setup were discussed in detail in , and thermal spurious forces are expected to show an analogous behavior. Therefore, we do not study spurious forces here in detail to maintain the focus on the thermomechanics of the GPP formulation.

Thermal effects on nanoindentation of copper
Finally, we apply the thermomechanical transport model to the case of nanoindentation into a Cu single-crystal. While the problem is well studied in 2D using finite-temperature QC implementations (Tadmor et al., 2013), only few QC studies exist that study finite-temperature effects in 3D. Kulkarni (2007) studied nanoindentation of a 32 × 32 × 32 unit cell FCC nearest-neighbor Lennard-Jones crystal, using the Wentzel-Kramers-Brillouin (WKB) approximation, which captures the thermomechanical coupling at comparably low temperature only. We here perform nanoindentation simulations of a Cu cube of 0.077 μm side length (approximately 215 × 215 × 215 unit cells, see Fig. 17), modeled by the EAM potential of Dai et al. (2006), underneath a 5 nm spherical indenter modeled using the indenter potential of Kelchner et al. (1998) with a force constant of 1000 eV/Å 3 and a maximum displacement of 1.26 nm or approximately 3.5 times the lattice parameter at 0 K. The crystal consists of approximately 50 million atoms, represented in the QC framework by approximately 0.2 million repatoms. The top surface is modeled as a free surface, while all other boundaries suppress wall-normal displacements, allowing only in-plane motion. Below, we discuss the results for isothermal deformation at = 0 K, 300 K, 600 K and for quasistatic adiabatic deformation, the latter being initially at = 300 K. Fig. 18 shows the variation of the total force on the spherical indenter vs. indentation depth for both isothermal and quasistatic adiabatic conditions. The force increases nonlinearly with indentation depth, showing the typical Hertzian-type initial elastic section. With increasing indentation depth, atoms underneath the indenter generate dislocations and stacking faults, overall creating a complex microstructure consisting of prismatic dislocation loops (PDLs), as also observed by Ponga et al. (2015) in nano-void growth in Cu. These PDLs are created in P. Gupta et al. Fig. 18. Variation of the indenter force with the indenter depth for isothermal ( = 0 K, 300 K, 600 K) and adiabatic (initially at = 300 K) conditions. After an initial elastic regime, the curve shows the typical serrated behavior due to dislocation activity underneath the indenter.  As dislocations and stacking faults are created, local thermal gradients are generated which are diffused via the thermal transport model. ( ) Microstructure generated below the spherical indenter at an indentation depth of 1 nm. The generated dislocations move towards the boundaries of the atomistic region, creating stacking faults in the crystallographic planes of the {111}-family (shaded in gray). Red repatoms are top-surface repatoms and blue atoms denote those within the microstructure, identified using the centrosymmetry parameter (values > 5 identify surface and microstructure atoms only). Shaded repatoms denote all repatoms with reduced opacity, shown here for comparison of the size of the microstructure with the atomistic domain. the {111} slip planes of the crystal, as shown in Fig. 19. At the first dislocation activation, the indenter force drops after reaching a critical force. This critical force decreases with increasing temperature, indicating plastic softening. The dislocations move towards the boundaries of the atomistic region, gliding in the crystallographic planes of the {111}-family, giving way to stacking faults in those planes. As shown in Fig. 19, while the initial dislocations maintain their structure, the stacking fault structure changes significantly at the same indentation depth as temperature increases. Fig. 20 shows the spatial temperature distribution and the emergent microstructure during the quasistatic adiabatic deformation of the Cu single-crystal. With the dislocations, local temperature gradients are generated along the PDLs generated in the {111} slip planes due to large gradients in , which are triggered due to large deviations from a centrosymmetric neighborhood (as identified in Fig. 16). These temperature 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 simulations with initial conditions and/or indenter position slightly perturbed, since the emergence of microstructure below the indenter within the highly-symmetric single-crystal is associated with instability and strongly depends on local fluctuations and initial conditions. Such a statistical analysis is deferred to future work.

Conclusion and discussion
We have presented a Gaussian phase packets-based (GPP-based) formulation of finite-temperature equilibrium and nonequilibrium thermomechanics applied to atomistic systems. Application of the Liouville equation to an ansatz of the probability distribution yields equations of motion for all phase-space parameters, thus aiming for temporal upscaling by separating the slow mean motion of atoms from their statistical position and momentum variances at finite temperature. For a general parametrization (or a sequence of ansätze) of the Hamiltonian systems, we have presented the variational structure and weak formulation of the Liouville equation, which alternatively yields the phase-space parameter evolution equations. As a special choice, we have shown that approximating the global statistical distribution function by a multivariate Gaussian ansatz can capture thermal transport at the atomic scale only via interatomic correlations. Due to high computational costs, we have neglected such interatomic correlations, which results in a local GPP approximation of the system. Such a system exhibits reversible dynamics with thermomechanical coupling, causing local heating and cooling upon movement of atoms relative to local neighborhoods. Moreover, in the quasistatic limit, we have shown that the equations yield local mechanical and thermal equilibrium conditions, the latter yielding the local equation of state of the atoms based on the interatomic force field. We have shown that the thus-obtained equilibrium relations are identical to those of the maximum-entropy (max-ent ) ansatz (even though the latter postulates a probability distribution in the quasistatic limit). However, unlike max-ent, our approach yields the full dynamic evolution equations, thereby allowing us to gain insight into the thermomechanical evolution in phase space. We furthermore derived different constraint conditions to be applied to the equilibrium relations for modeling isothermal, isobaric, and isentropic processes. The setup has been validated by computing thermal expansion and elastic constants of a copper single-crystal.
To capture the irreversibility due to local thermal transport triggered by the adiabatic heating/cooling of atoms, we have coupled the quasistatic framework with linear Onsager kinetics. Such a model involves an empirical coefficient fitted to match approximate bulk conductivity measurements and captures the experimentally observed size-effects of the thermal conductivity, as shown by Venturini et al. (2014). Moreover, we have shown that the time scale imposed by the atomic-scale transport is approximately 100 times that of atomic vibrations. While the 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 a steady-state solution. Based on the global multivariate Gaussian ansatz, interatomic correlations may be fitted to obtain correlation functions (akin to interatomic potentials), which can help develop the transport constitutive properties of atomistic systems and also advance current understanding of the long-standing nanoscale thermal transport problem.
Finally, we have combined the quasistatic thermomechanical equations based on the local GPP approximation with thermal transport in a high-performance, distributed-memory, updated-Lagrangian 3D QC solver, which is capable of modeling thermomechanical deformation of large-scale systems by coarse-graining the atomistic ensemble in space. Benchmark simulations of dislocation nucleation and nanoindentation under isothermal and adiabatic constraints showed reasonable agreement between coarse-grained and fully resolved simulations. Since the time integration of the atomic transport can be terminated for a small heat flux imbalance (discussed in Algorithm 1), the quasistatic simulations offer significant advantages over traditional MD studies, which can tackle only high strain rates. The presented methodology of coupling local thermal equilibrium 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 the time-scale limitations of MD, a one-to-one comparison of the presented simulations with finite-temperature MD simulations is prohibitively costly (except for equilibrium properties like the presented thermal expansion coefficient). A detailed analysis of the accuracy of the spatial coarse-graining of the thermomechanical model presented here, and a comparison with suitable MD simulations qualifies as a possible extension of this work. For such comparisons, however, large-scale nonequilibrium molecular dynamics (NEMD) simulations are required at sufficiently slow strain rates, which go beyond the scope of this study.  Applying the above limit to a system in which the th atom has multiple neighbors, we obtain the constraint in Eq. (3.38).