where indices $1$ and $2$ refer to the respective components of the mean
position vector in the chosen coordinate system, and $\Delta\gamma$ is the
applied shear strain increment. As the strain is applied, repatoms in the
inner region B (Figure 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 $S_{\Sigma}$ 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
increased, the horizontal component of the force on all repatoms in region A
is computed. The effective shear stress $\tau$ on the $(\overline{1}11)$-plane
is computed by normalizing the net horizontal force by the cross-sectional
area of the slab. Results are shown in Figure 12$(c)$ and $(d)$. 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 Figure 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 Figure
13), both capturing the apparent thermal plastic softening of the crystal.

Figure 16: Local variation of temperature as the slab with void is deformed
adiabatically. In the intial stages, the temperature rises slowly due to the
external deformation. At the critical shear strain
$\gamma_{\mathrm{cr}}=n_{\mathrm{cr}}\Delta\gamma$, 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
$n=n^{*}_{\text{cr}}$ before the irreversible transport exceeds the colorbar
range.

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. However, in our
analysis (which is based on a more realistic crystallographic setup since the
slip planes of the dislocations are close-packed planes of $\\{111\\}$ family)
a critical shear stress may be defined as the shear stress required to
nucleate dislocations from the void-like defect.

####  4.2.2 Adiabatic

To simulate the quasistatic adiabatic activation 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.4 and Algorithm 1. As noted
above, since the boundaries do not allow any heat exchange out of the domain
to the surrounding (thermally insulated), the term _adiabatic_ is used to
describe the setup. We further assume that the strain-rate is singificantly
lower than the time scale imposed by the molecular thermal transport, thus
imposing quasistatic conditions for the transport ($\delta t^{(n)}\to\infty$
in Algorithm 1, see Step 3 at the end of Section 3.4).  As we describe below,
the thermomechanical deformation approaches isothermal conditions since the
quasistatic assumption results in homogenization of the temperature field to a
large extent. 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. Figure 15 shows the
variation of the shear stress on the $(\overline{1}11)$-plane with external
shear strain $\gamma$. 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. Figure 16 shows the spatial
variation of temperature as the slab is deformed adiabatically. Before the
critical strain, heating caused by local deformation is negligible. As the
dislocations are nucleated at the critical strain $\gamma_{\mathrm{cr}}$, the
temperature of those atoms around the dislocations changes significantly, as
shown by the intermediate stage $(\overline{\boldsymbol{q}}^{*},T^{*})$ in
Figure 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 feature of thermomechanical deformation at slow
strain rates was observed by (Ponga et al., 2016) while studying strain rate
effects on nano-void growth in magnesium and (Ponga and Sun, 2018) while
studying 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
(Figure 12).

As noted above, the critical stress values obtained from QC simulations are
within $6$% of those obtained from the fully atomistic simulations.
Furthermore, the critical strain values are within $12$%. Repatoms on the
$(110)$\- and $(\overline{1}11)$-surfaces, which includes the repatoms on
vertices of very large elements and the repatoms in the transition region
between regions A and B, exhibit both mechanical and thermal spurious forces
in 3D (Amelang et al., 2015; Amelang and Kochmann, 2015; 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 Amelang et
al. (2015), 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.

$(a)$$(b)$$(c)$

Figure 17: Illustration of the undeformed and deformed QC meshes of a $0.077$
$\mu$m FCC single-crystal of pure Cu. $(a)$ Cross-section of the initial,
undeformed mesh, $(b,c)$ zoomed-in perspective views of the atomistic region
(33$\times$33$\times$33 unit cells) and the surrounding coarsened regions in
the ($b$) undeformed and $(c)$ deformed configurations underneath a 5 nm
spherical indenter at an indentation depth of $0.75$ nm. Figure 18: Variation
of the indenter force with the indenter depth for isothermal
($T=0~{}\mathrm{K},~{}300~{}\mathrm{K},~{}600~{}\mathrm{K}$) and adiabatic
(initially at $T=300~{}\mathrm{K}$) conditions. After an initial elastic
region, the curve shows the typical serrated behavior due to dislocation
activity underneath the indenter.

###  4.3 Thermal effects on nanoindentation of copper

Figure 19: Microstructure generated below the 5 nm spherical indenter at an
indentation depth of $1$ nm for isothermal deformation of a Cu single-crystal.
The generated dislocations move towards the boundaries of the atomistic
region, creating stacking faults in the crystallographic planes of the
$\\{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 boundary and microstructure
atoms only). All repatoms are shown with reduced opacity for comparison of
size of the microstructure with the atomistic domain.

$(a)$$(b)$

Figure 20: Adiabatic deformation of the Cu single-crystal under a spherical
indenter. $(a)$ Spatial variation of temperature in the cross-section of the
atomistic region. As dislocations and stacking faults are created, local
thermal gradients are generated which are diffused via the thermal transport
model. $(b)$ 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 boundary and microstructure atoms only). Shaded repatoms denote all
repatoms with reduced opacity, shown here for comparison of size of the
microstructure with the atomistic domain.

Finally, we apply the discussed 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 studying the finite-temperature effects in 3D. Kulkarni
(2007) studied the nanoindentation of a $32\times 32\times 32$ unit cell FCC
nearest-neighbour 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$ $\mu$m side length (approximately
215$\times$215$\times$215 unit cells, see Figure 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 boundary, while all other
boundaries suppress wall-normal displacements, allowing only in-plane motion.
Below, we discuss the results for isothermal deformation at
$T=0~{}\mathrm{K},~{}300~{}\mathrm{K},~{}600~{}\mathrm{K}$ and for quasistatic
adiabatic deformation, the latter being initially at $T=300~{}\mathrm{K}$.

Figure 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 observed by Ponga et al. (2015) in nano-void growth in a
copper crystal. These PDLs are created along the slip planes of the crystal of
$\\{111\\}$ family, as shown in Figure 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 with increasing temperature. 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
Figure 19, while the initial dislocations maintain their structure, the
stacking fault structure changes significantly at the same indentation depth
as temperature increases.

Figure 20 shows the spatial temperature distribution and the emergent
microstructure during the quasistatic adiabatic deformation of the Cu crystal.
With the dislocations local temperature gradients are generated along the PDLs
generated in the slip planes of $\\{111\\}$ type due to large gradients in
$S_{\Sigma}$, which are triggered due to large deviations from a
centrosymmetric neighborhood (as identified in Figure 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.

##  5 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. We have shown that approximating the global statistical
distribution function with a multivariate Gaussian ansatz enables capturing of
thermal transport only via interatomic correlations. Due to high computational
costs, we have neglected the interatomic correlations, which results in a
local GPP approximation of the system. Such a system exhibits reversible
dynamics with 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. 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 obtain 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 as 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 convincing agreement between coarse-grained and fully
resolved atomistic simulations. Since the time integration of the atomic
transport can be terminated for 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. 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, very large scale nonequilibrium
molecular dynamics (NEMD) simulations are required at sufficiently slow strain
rates.

## Acknowledgments

The support from the European Research Council (ERC) under the European
Union’s Horizon 2020 research and innovation program (grant agreement no.
770754) is gratefully acknowledged. The authors thank Michael Ortiz for
stimulating discussions and Miguel Spinola for aiding with the LAMMPS
simulations for numerical validation.

## References

  * Admal and Tadmor (2010) Admal, N. C., Tadmor, E. B., 2010. A unified interpretation of stress in molecular systems. Journal of Elasticity 100 (1-2), 63–143.
  * Amelang (2016) Amelang, J. S., 2016. A fully-nonlocal energy-based formulation and high-performance realization of the quasicontinuum method. Ph.D. thesis.
  * Amelang and Kochmann (2015) Amelang, J. S., Kochmann, D. M., 2015. Surface effects in nanoscale structures investigated by a fully-nonlocal energy-based quasicontinuum method. Mechanics of Materials 90, 166 – 184.
  * Amelang et al. (2015) Amelang, J. S., Venturini, G. N., Kochmann, D. M., 2015. Summation rules for a fully nonlocal energy-based quasicontinuum method. Journal of the Mechanics and Physics of Solids 82, 378–413.
  * Ariza et al. (2012) Ariza, M., Romero, I., Ponga, M., Ortiz, M., 2012. Hotqc simulation of nanovoid growth under tension in copper. International Journal of Fracture 174 (1), 75–85.
  * Belytschko and Xiao (2003) Belytschko, T., Xiao, S., 2003. Coupling methods for continuum model with molecular model. International Journal for Multiscale Computational Engineering 1 (1).
  * Bitzek et al. (2006) Bitzek, E., Koskinen, P., Gähler, F., Moseler, M., Gumbsch, P., 2006. Structural relaxation made simple. Physical review letters 97 (17), 170201.
  * Brune et al. (2015) Brune, P. R., Knepley, M. G., Smith, B. F., Tu, X., Jan 2015. Composing scalable nonlinear algebraic solvers. SIAM Review 57 (4), 535–565.
URL http://dx.doi.org/10.1137/130936725

  * Butcher (1986) Butcher, P., 1986. The theory of electron transport in crystalline semiconductors. In: Crystalline Semiconducting Materials and Devices. Springer, pp. 131–194.
  * Chang and Himmel (1966) Chang, Y., Himmel, L., 1966. Temperature dependence of the elastic constants of cu, ag, and au above room temperature. Journal of Applied Physics 37 (9), 3567–3572.
  * Chen et al. (2017) Chen, X., Xiong, L., McDowell, D. L., Chen, Y., 2017. Effects of phonons on mobility of dislocations and dislocation arrays. Scripta Materialia 137, 22–26.
  * Chen (2009) Chen, Y., 2009. Reformulation of microscopic balance equations for multiscale materials modeling. The Journal of chemical physics 130 (13), 134706.
  * Dai et al. (2006) Dai, X., Kong, Y., Li, J., Liu, B., 2006. Extended finnis–sinclair potential for bcc and fcc metals and alloys. Journal of Physics: Condensed Matter 18 (19), 4527.
  * Dobson et al. (2010) Dobson, M., Luskin, M., Ortner, C., 2010. Accuracy of quasicontinuum approximations near instabilities. Journal of the Mechanics and Physics of Solids 58 (10), 1741–1757.
  * Eidel and Stukowski (2009) Eidel, B., Stukowski, A., 2009. A variational formulation of the quasicontinuum method based on energy sampling in clusters. Journal of the Mechanics and Physics of Solids 57 (1), 87–108.
  * Espanol et al. (2013) Espanol, M. I., Kochmann, D. M., Conti, S., Ortiz, M., 2013. A $\gamma$-convergence analysis of the quasicontinuum method. Multiscale Modeling & Simulation 11 (3), 766–794.
  * Evans and Morriss (2007) Evans, D. J., Morriss, G. P., 2007. Statistical Mechanics of Nonequilbrium Liquids. ANU Press.
  * Gunzburger and Zhang (2010) Gunzburger, M., Zhang, Y., 2010. A quadrature-rule type approximation to the quasi-continuum method. Multiscale Modeling & Simulation 8 (2), 571–590.
  * Heller (1975) Heller, E. J., 1975. Time-dependent approach to semiclassical dynamics. The Journal of Chemical Physics 62 (4), 1544–1555.
  * Hicks and Dresselhaus (1993) Hicks, L., Dresselhaus, M. S., 1993. Thermoelectric figure of merit of a one-dimensional conductor. Physical review B 47 (24), 16631.
  * Hirth (1980) Hirth, J. P., 1980. Effects of hydrogen on the properties of iron and steel. Metallurgical Transactions A 11 (6), 861–890.
  * Hull and Bacon (2001) Hull, D., Bacon, D. J., 2001. Introduction to dislocations. Butterworth-Heinemann.
  * Irving and Kirkwood (1950) Irving, J., Kirkwood, J. G., 1950. The statistical mechanical theory of transport processes. iv. the equations of hydrodynamics. The Journal of chemical physics 18 (6), 817–829.
  * Iyer and Gavini (2011) Iyer, M., Gavini, V., 2011. A field theoretical approach to the quasi-continuum method. Journal of the Mechanics and Physics of Solids 59 (8), 1506–1535.
  * Johnson (1988) Johnson, R., 1988. Analytic nearest-neighbor model for fcc metals. Physical Review B 37 (8), 3924.
  * Kelchner et al. (1998) Kelchner, C. L., Plimpton, S., Hamilton, J., 1998. Dislocation nucleation and defect structure during surface indentation. Physical Review B 58 (17), 11085\.
  * Kim et al. (2014a) Kim, S.-P., Datta, D., Shenoy, V. B., 2014a. Atomistic mechanisms of phase boundary evolution during initial lithiation of crystalline silicon. The Journal of Physical Chemistry C 118 (31), 17247–17253.
  * Kim et al. (2014b) Kim, W. K., Luskin, M., Perez, D., Voter, A., Tadmor, E. B., 2014b. Hyper-qc: An accelerated finite-temperature quasicontinuum method using hyperdynamics. Journal of the Mechanics and Physics of Solids 63, 94–112.
  * Kittel (1976) Kittel, C., 1976. Introduction to solid state physics.
  * Knap and Ortiz (2001) Knap, J., Ortiz, M., 2001. An analysis of the quasicontinuum method. Journal of the Mechanics and Physics of Solids 49 (9), 1899–1923.
  * Kulkarni (2007) Kulkarni, Y., 2007. Coarse-graining of atomistic description at finite temperature. Ph.D. thesis, California Institute of Technology.
  * Kulkarni et al. (2008) Kulkarni, Y., Knap, J., Ortiz, M., 2008. A variational approach to coarse graining of equilibrium and non-equilibrium atomistic description at finite temperature. Journal of the Mechanics and Physics of Solids 56 (4), 1417–1449.
  * Landau and Lifshitz (1980) Landau, L., Lifshitz, E., 1980. Statistical Physics, Part 1: Volume 5.
  * Lepri et al. (2003) Lepri, S., Livi, R., Politi, A., 2003. Thermal conduction in classical low-dimensional lattices. Physics Reports 377 (1), 1–80.
  * Li et al. (2011) Li, J., Sarkar, S., Cox, W. T., Lenosky, T. J., Bitzek, E., Wang, Y., 2011. Diffusive molecular dynamics and its application to nanoindentation and sintering. Physical Review B 84 (5), 054103.
  * Ma et al. (1993) Ma, J., Hsu, D., Straub, J. E., 1993. Approximate solution of the classical liouville equation using gaussian phase packet dynamics: Application to enhanced equilibrium averaging and global optimization. The Journal of Chemical Physics 99 (5), 4024–4035.
  * Marian et al. (2009) Marian, J., Venturini, G. N., Hansen, B., Knap, J., Ortiz, M., Campbell, G., 2009\. Finite-temperature non-equilibrium quasicontinuum method based on langevin dynamics. Modelling and Simulation in Materials Science and Engineering, vol. 18, N/A, December 10, 2001, pp. 015003 18 (LLNL-JRNL-412903).
  * McLachlan (1964) McLachlan, A., 1964. A variational solution of the time-dependent schrodinger equation. Molecular Physics 8 (1), 39–44.
  * Mehta et al. (2015) Mehta, R., Chugh, S., Chen, Z., 2015. Enhanced electrical and thermal conduction in graphene-encapsulated copper nanowires. Nano letters 15 (3), 2024–2030.
  * Mendez et al. (2018) Mendez, J., Ponga, M., Ortiz, M., 2018. Diffusive molecular dynamics simulations of lithiation of silicon nanopillars. Journal of the Mechanics and Physics of Solids 115, 123–141.
  * Miller et al. (1998) Miller, R., Tadmor, E., Phillips, R., Ortiz, M., 1998. Quasicontinuum simulation of fracture at the atomic scale. Modelling and Simulation in Materials Science and Engineering 6 (5), 607.
  * Miller and Tadmor (2009) Miller, R. E., Tadmor, E. B., 2009. A unified framework and performance benchmark of fourteen multiscale atomistic/continuum coupling methods. Modelling and Simulation in Materials Science and Engineering 17 (5), 053001.
  * Motamarri et al. (2020) Motamarri, P., Das, S., Rudraraju, S., Ghosh, K., Davydov, D., Gavini, V., 2020\. Dft-fe–a massively parallel adaptive finite-element code for large-scale density functional theory calculations. Computer Physics Communications 246, 106853.
  * Nabarro (1967) Nabarro, F. R. N., 1967. Theory of crystal dislocations. Clarendon Press.
  * Nix and MacNair (1941) Nix, F., MacNair, D., 1941. The thermal expansion of pure metals: copper, gold, aluminum, nickel, and iron. Physical Review 60 (8), 597.
  * Overton Jr and Gaffney (1955) Overton Jr, W., Gaffney, J., 1955. Temperature variation of the elastic constants of cubic elements. I. Copper. Physical Review 98 (4), 969.
  * Plimpton (1995) Plimpton, S., 1995. Fast parallel algorithms for short-range molecular dynamics. Journal of computational physics 117 (1), 1–19.
  * Ponga et al. (2015) Ponga, M., Ortiz, M., Ariza, M., 2015. Finite-temperature non-equilibrium quasi-continuum analysis of nanovoid growth in copper at low and high strain rates. Mechanics of Materials 90, 253–267.
  * Ponga et al. (2016) Ponga, M., Ramabathiran, A. A., Bhattacharya, K., Ortiz, M., 2016. Dynamic behavior of nano-voids in magnesium under hydrostatic tensile stress. Modelling and Simulation in Materials Science and Engineering 24 (6), 065003.
  * Ponga and Sun (2018) Ponga, M., Sun, D., 2018. A unified framework for heat and mass transport at the atomic scale. Modelling and Simulation in Materials Science and Engineering 26 (3), 035014.
  * Qu et al. (2005) Qu, S., Shastry, V., Curtin, W., Miller, R. E., 2005. A finite-temperature dynamic coupled atomistic/discrete dislocation method. Modelling and Simulation in Materials Science and Engineering 13 (7), 1101.
  * Reddy (2007) Reddy, J. N., 2007. An introduction to continuum mechanics. Cambridge university press.
  * Sääskilahti et al. (2015) Sääskilahti, K., Oksanen, J., Volz, S., Tulkki, J., 2015. Frequency-dependent phonon mean free path in carbon nanotubes from nonequilibrium molecular dynamics. Physical Review B 91 (11), 115426.
  * Shilkrot et al. (2002) Shilkrot, L., Miller, R. E., Curtin, W., 2002. Coupled atomistic and discrete dislocation plasticity. Physical review letters 89 (2), 025501.
  * Srivastava and Nemat-Nasser (2014) Srivastava, A., Nemat-Nasser, S., 2014. On the limit and applicability of dynamic homogenization. Wave Motion 51 (7), 1045–1054.
  * Stroud (1971) Stroud, A. H., 1971. Approximate calculation of multiple integrals. Prentice-Hall.
  * Tadmor et al. (2013) Tadmor, E. B., Legoll, F., Kim, W., Dupuy, L., Miller, R. E., 2013. Finite-temperature quasi-continuum. Applied Mechanics Reviews 65 (1).
  * Tadmor and Miller (2011) Tadmor, E. B., Miller, R. E., 2011. Modeling materials: continuum, atomistic and multiscale techniques. Cambridge University Press.
  * Tadmor et al. (1996a) Tadmor, E. B., Ortiz, M., Phillips, R., 1996a. Quasicontinuum analysis of defects in solids. Philosophical Magazine A 73 (6), 1529–1563.
  * Tadmor et al. (1996b) Tadmor, E. B., Phillips, R., Ortiz, M., 1996b. Mixed atomistic and continuum models of deformation in solids. Langmuir 12 (19), 4529–4534.
  * Tembhekar (2018) Tembhekar, I., 2018. The fully nonlocal, finite-temperature, adaptive 3d quasicontinuum method for bridging across scales. Ph.D. thesis, California Institute of Technology.
  * Tembhekar et al. (2017) Tembhekar, I., Amelang, J. S., Munk, L., Kochmann, D. M., 2017. Automatic adaptivity in the fully nonlocal quasicontinuum method for coarse-grained atomistic simulations. International Journal for Numerical Methods in Engineering 110 (9), 878–900.
  * Tuckerman (2010) Tuckerman, M., 2010. Statistical mechanics: theory and molecular simulation. Oxford University Press.
  * van der Giessen et al. (2020) van der Giessen, E., Schultz, P. A., Bertin, N., Bulatov, V. V., Cai, W., Csányi, G., Foiles, S. M., Geers, M. G., González, C., Hütter, M., et al., 2020. Roadmap on multiscale materials modeling. Modelling and Simulation in Materials Science and Engineering 28 (4), 043001.
  * Venturini et al. (2014) Venturini, G., Wang, K., Romero, I., Ariza, M., Ortiz, M., 2014. Atomistic long-term simulation of heat and mass transport. Journal of the Mechanics and Physics of Solids 73, 242–268.
  * Venturini (2011) Venturini, G. N., 2011. Topics in multiscale modeling of metals and metallic alloys. Ph.D. thesis, California Institute of Technology.
  * Voter (1997) Voter, A. F., 1997. A method for accelerating the molecular dynamics simulation of infrequent events. The Journal of chemical physics 106 (11), 4665–4677.
  * Wagner et al. (2008) Wagner, G. J., Jones, R., Templeton, J., Parks, M., 2008. An atomistic-to-continuum coupling method for heat transfer in solids. Computer Methods in Applied Mechanics and Engineering 197 (41-42), 3351–3365.
  * Weiner (2012) Weiner, J. H., 2012. Statistical mechanics of elasticity. Courier Corporation.
  * Xu and Chen (2019) Xu, S., Chen, X., 2019. Modeling dislocations and heat conduction in crystalline materials: atomistic/continuum coupling approaches. International Materials Reviews 64 (7), 407–438.
  * Zhang and Curtin (2008) Zhang, F., Curtin, W., 2008. Atomistically informed solute drag in al–mg. Modelling and Simulation in Materials Science and Engineering 16 (5), 055006.
  * Ziman (2001) Ziman, J. M., 2001. Electrons and phonons: the theory of transport phenomena in solids. Oxford University Press.
  * Zubarev (1974) Zubarev, D. N., 1974. Nonequilibrium statistical thermodynamics. Consultants Bureau.
  * Zubarev et al. (1996) Zubarev, D. N., Morozov, V., Ropke, G., 1996. Statistical Mechanics of Nonequilibrium Processes, Volume 1 (See 3527400834): Basic Concepts, Kinetic Theory. Vol. 1. Wiley-VCH.

##  Appendix A Time evolution of phase averaged quantities

The time evolution of the phase average of a phase-space quantity
$A(\boldsymbol{z})$ can be derived using the representative solution of the
Liouville equation,

$\frac{\partial f}{\partial t}+i\mathcal{L}f=0\ \implies\
f(\boldsymbol{z},t)=e^{-i\mathcal{L}t}f(\boldsymbol{z}_{0},0)=e^{-i\mathcal{L}t}f(\boldsymbol{z}),$
(A.1)

where $f(\boldsymbol{z}_{0},0)$ is the initial condition and
$e^{-i\mathcal{L}t}$ is the propagating operator, which transforms the
probability distribution initially defined at phase-space coordinate
$\boldsymbol{z}_{0}$ to the probability distribution at $\boldsymbol{z}(t)$
(Evans and Morriss, 2007; Zubarev et al., 1996). Furthermore, the time
evolution of the phase-space quantity $A(\boldsymbol{z})$ is given by

$\frac{\;\\!\mathrm{d}A}{\;\\!\mathrm{d}t}=i\mathcal{L}A\ \implies\
A(\boldsymbol{z},t)=e^{i\mathcal{L}t}A(\boldsymbol{z}_{0},0)=e^{i\mathcal{L}t}A(\boldsymbol{z}).$
(A.2)

Equations (A.1) and (A.2) reveal that the operators $e^{\pm i\mathcal{L}t}$
transport the probability distribution $f(\boldsymbol{z})$ and phase-space
quantities $A(\boldsymbol{z})$ defined in terms of $\boldsymbol{z}$ of a
system of particles, given that $\boldsymbol{z}$ also changes in time as the
system of particles evolves.  Operator $i\mathcal{L}$ satisfies the property

$\int_{\Gamma}A(\boldsymbol{z})i\mathcal{L}f(\boldsymbol{z})d\boldsymbol{z}=\int_{\Gamma}(-i\mathcal{L})A(\boldsymbol{z})f(\boldsymbol{z})d\boldsymbol{z}$
(A.3)

for real-valued $A$ and $f\to 0$ as $\boldsymbol{z}\to\partial\Gamma$ where
$\partial\Gamma$ is the boundary of $\Gamma$ and
$\Gamma\subseteq\mathbb{R}^{6N}$. For a $\Gamma$ almost covering
$\mathbb{R}^{6N}$, $f(\boldsymbol{z})$ approaches 0 as any component of
momentum approaches infinity or any spatial dimension approaches infinity
(probability of finding classical atoms far away from their mean position must
decay to 0) . Using this property, the time evolution of the phase average of
a phase-space quantity $A(\boldsymbol{z})$, defined in terms of
$\boldsymbol{z}$, is obtained from

$N!\,h^{3N}\frac{\;\\!\mathrm{d}\left\langle{A}\right\rangle}{dt}=\frac{\;\\!\mathrm{d}}{\;\\!\mathrm{d}t}\int_{\Gamma}f(\boldsymbol{z},t)A(\boldsymbol{z})\;\\!\mathrm{d}\boldsymbol{z}=\int_{\Gamma}\frac{\;\\!\mathrm{d}}{\;\\!\mathrm{d}t}\left[e^{-i\mathcal{L}t}f(\boldsymbol{z})\right]A(\boldsymbol{z})\;\\!\mathrm{d}\boldsymbol{z}=\int_{\Gamma}e^{-i\mathcal{L}t}f(\boldsymbol{z})i\mathcal{L}tA(\boldsymbol{z})\;\\!\mathrm{d}\boldsymbol{z}.$
(A.4)

Using equation (A.2), we obtain

$\frac{\;\\!\mathrm{d}\left\langle{A}\right\rangle}{\;\\!\mathrm{d}t}=\frac{1}{N!\,h^{3N}}\int_{\Gamma}f(\boldsymbol{z},t)\frac{\;\\!\mathrm{d}A}{\;\\!\mathrm{d}t}\;\\!\mathrm{d}\boldsymbol{z}=\left\langle{\frac{\;\\!\mathrm{d}A}{\;\\!\mathrm{d}t}}\right\rangle.$
(A.5)

We note that equation (A.5) is obtained using only the evolution equations
(A.1), (A.2) and property (A.3), which hold for any Hamiltonian system and
thus contain no time-coarsening approximations. Accordingly, time-variational
formulations such as the Frenkel-Dirac-McLachlan variational principle
(McLachlan, 1964) leads to identical equations.

##  Appendix B Quasistatic GPP as Helmholtz free energy minimization

The Helmholtz free energy $\mathcal{F}$ as a function of parameter set
$(\overline{\boldsymbol{q}},S_{\Sigma},S_{\Omega})$ is defined as

$\mathcal{F}(\overline{\boldsymbol{q}},S_{\Sigma},S_{\Omega})=E(\overline{\boldsymbol{q}},S_{\Sigma},S)-\sum_{i}\frac{\Omega_{i}S_{i}}{k_{B}m_{i}}$
(B.1)

with the relation

$\frac{\Omega_{i}}{k_{B}m_{i}}=\frac{\partial E}{\partial S_{i}}.$ (B.2)

Minimization of $\mathcal{F}$ with respect to the set
$\overline{\boldsymbol{q}}$ yields

$-\frac{\partial\mathcal{F}}{\partial\overline{\boldsymbol{q}}_{i}}=0\
\implies\ -\left\langle{\frac{\partial
V(\boldsymbol{q})}{\partial\boldsymbol{q}_{i}}}\right\rangle=\left\langle{F_{i}(\boldsymbol{q})}\right\rangle=0.$
(B.3)

To minimize $\mathcal{F}$ with respect to the set $S_{\Sigma}$, we consider
the relation

$\boldsymbol{q}_{i}-\overline{\boldsymbol{q}}_{i}=\sqrt{\Sigma_{i}}\boldsymbol{x}_{i},$
(B.4)

for some normalized vector $\boldsymbol{x}_{i}$. From equation (B.4) it
follows that

$\frac{\partial\boldsymbol{q}_{i}}{\partial
S_{\Sigma,i}}=\frac{\partial\boldsymbol{q}_{i}}{\partial\sqrt{\Sigma_{i}}}\frac{\partial\sqrt{\Sigma_{i}}}{\partial
S_{\Sigma,i}}=\sqrt{\Sigma_{i}}\boldsymbol{x}_{i}=\boldsymbol{q}_{i}-\overline{\boldsymbol{q}}_{i}.$
(B.5)

Finally, minimization of $\mathcal{F}$ with respect to the set $S_{\Sigma}$
yields the following set of equations:

$-\frac{\partial\mathcal{F}}{\partial S_{\Sigma,i}}=0\ \implies\
\frac{3\Omega_{i}}{m_{i}}-\left\langle{\frac{\partial
V(\boldsymbol{q})}{\partial\boldsymbol{q}_{i}}\cdot\frac{\partial\boldsymbol{q}_{i}}{\partial
S_{\Sigma,i}}}\right\rangle=\frac{3\Omega_{i}}{m_{i}}+\left\langle{\boldsymbol{F}_{i}(\boldsymbol{q})\cdot\left(\boldsymbol{q}_{i}-\overline{\boldsymbol{q}}_{i}\right)}\right\rangle=0.$
(B.6)

Equations (B.3) and (B.6) are identical to the quasistatic GPP equations
(3.12b).

##  Appendix C Time step stability bounds for entropy transport

Applying a forward-Euler explicit time discretization to equation (3.33), we
obtain

$\frac{1}{\Delta t^{(k)}}\left(\left(\begin{matrix}T_{i}\\\
T_{j}\end{matrix}\right)^{(k+1)}-\left(\begin{matrix}T_{i}\\\
T_{j}\end{matrix}\right)^{(k)}\right)=\frac{2A_{0}}{3k_{B}}\frac{T^{2,(k)}_{ij}}{T^{(k)}_{i}T^{(k)}_{j}}\left(\begin{matrix}T_{i}-T_{j}\\\
T_{j}-T_{i}\end{matrix}\right)^{(k)},$ (C.1)

where superscript $(k)$ implies a quantity evaluated at the $k^{\text{th}}$
time step. Rearranging the above equation yields

$\left(\begin{matrix}T_{i}\\\
T_{j}\end{matrix}\right)^{(k+1)}=\boldsymbol{T}^{(k)}\left(\begin{matrix}T_{i}\\\
T_{j}\end{matrix}\right)^{(k)},$ (C.2)

where $\boldsymbol{T}^{(k)}$ is the transition matrix at the $k^{\mathrm{th}}$
time step, defined by

$\boldsymbol{T}^{(k)}=\boldsymbol{I}+\frac{2A_{0}\Delta
t^{(k)}}{3k_{B}}\frac{T^{2,(k)}_{ij}}{T^{(k)}_{i}T^{(k)}_{j}}\left(\begin{matrix}1&-1\\\
-1&1\end{matrix}\right).$ (C.3)

For numerical stability, the transition matrix must have eigenvalues with
magnitude $\leq 1$, which yields the bound

$\frac{2A_{0}\Delta
t^{(k)}}{3k_{B}}\left(\frac{T^{2,(k)}_{ij}}{T^{(k)}_{i}T^{(k)}_{j}}\right)\leq
1.$ (C.4)

Applying the above limit to a system in which the $i^{\mathrm{th}}$ atom has
multiple neighbors, we obtain the constraint in equation (3.36).