In this work, we present a concurrent atomistic-continuum (CAC) method for modeling and simulation of crystalline materials. The CAC formulation extends the Irving-Kirkwood procedure for deriving transport equations and fluxes for homogenized molecular systems to that for polyatomic crystalline materials by employing a concurrent two-level description of the structure and dynamics of crystals. A multiscale representation of conservation laws is formulated, as a direct consequence of Newton's second law, in terms of instantaneous expressions of unit cell-averaged quantities using the mathematical theory of distributions. Finite element (FE) solutions to the conservation equations, as well as fluxes and temperature in the FE representation, are introduced, followed by numerical examples of the atomic-scale structure of interfaces, dynamics of fracture and dislocations, and phonon thermal transport across grain boundaries. In addition to providing a methodology for concurrent multiscale simulation of transport processes under a single theoretical framework, the CAC formulation can also be used to compute fluxes (stress and heat flux) in atomistic and coarse-grained atomistic simulations.

## I. INTRODUCTION

Solid state physics describes the structure of all crystals in terms of primitive unit cells. These primitive unit cells, each containing an identical group of atoms, are the building blocks of crystals. When added continuously in space, they form the structure of a crystal and completely fill the space the crystal occupies,^{1} cf. Fig. 1.

Although the motion of an atom in a crystal is an irregularly fluctuating function of time, it turns out that the dynamics of a crystal can be most readily described, not in terms of individual atoms, but in terms of traveling waves, named lattice vibrations by Max Born. These vibrational modes are called phonons and are of two types: acoustic and optical, cf. Fig. 2. With the former, atoms in a unit cell move in the same phase, resulting in the deformation of the lattice; with the latter, atoms undergo relative motion within the lattice cells, leaving the lattice unchanged. Thus, the displacements of atoms in a crystal can be described as the sum of continuous lattice displacements and the internal displacements of the atoms relative to the lattice.^{2}

It is seen that the solid state physics description of the structure and dynamics of crystals is based on a concurrent atomistic-continuum (CAC) approach: the structure is continuous at the lattice level but discrete at the atomic level; the lattice deformation is continuous until a defect is generated, cf. Fig. 3, while the internal motion is discrete.

Historically, atomistic and continuum descriptions offer two fundamentally different approaches to our understanding of crystalline matter. From the atomistic viewpoint, matter consists of discrete particles; while from the continuum viewpoint, matter is infinitely divisible. The continuum view leads to formulations of field equations of conservation laws. Supplemented by constitutive relations, such as Hooke's law and Fourier's law, these conservation equations serve as the governing equations in continuum modeling of materials.

This work introduces the concurrent atomistic-continuum (CAC) method that links and unifies the atomistic and continuum descriptions of crystalline materials. The Irving-Kirkwood procedure for linking the molecular description to hydrodynamical equations is reviewed in Sec. II; the CAC formalism is introduced in Sec. III; the finite element (FE) formulation of CAC as well as the fluxes and temperature in the FE representation is outlined in Sec. IV; numerical examples are presented in Sec. V, followed by a summary in Sec. VI. Detailed derivations of the conservation equations using the mathematical theory of distributions for monatomic crystals are included in Appendix A and that for polyatomic crystals in Appendix B.

## II. THE IRVING-KIRKWOOD PROCEDURE

The statistical mechanics formulation to link discrete and continuum descriptions of conservation laws was pioneered by Irving and Kirkwood (IK) in 1950.^{3} In their landmark paper, IK defined local densities in terms of molecular variables using the Dirac delta for “localized density,”^{4} “since mass or momentum of any molecule may be considered as localized at that molecule.”^{3} A local density of a dynamic phase function $A(rk,vk)$, where $rk$ and$vk$ are the position and velocity of the *k*th particle, is defined as an ensemble-averaged point function $a\xaf(x,t)$ as

where the symbol “≡” means “equal by definition.”

Equation (1) is usually interpreted as mapping the phase functions in the 6N-dimensional phase space to the local densities in the 4-dimensions of physical space and time. Based on the local density definition in Eq. (1), IK derived the rate of change for the densities of mass, momentum, and energy in the form of partial differential equations that relate the density of a conserved quantity $a\xaf(x,t)$ and its flux $J\xaf(x,t)$ as

Formulas for fluxes (stress and heat flux) are then obtained through comparing the rate equations with the differential equations of hydrodynamics.

The IK formalism has inspired numerous efforts to link molecular variables to field quantities in classical continuum mechanics.^{5–8} To understand the motivation of a new formalism, we recall the following noted by Irving and Kirkwood^{3} and by Kirkwood.^{9}

Local densities and fluxes in the IK formulation are point functions. “The point functions, although averaged neither over space nor time, satisfy equations that are identical in form to the equations of hydrodynamics for a single component, single phase system.”

^{3}These equations were called “the hydrodynamical-like equations” in the IK paper.^{3}It was noted that “to obtain the hydrodynamical equations themselves it is merely necessary to perform the appropriate space and time averages.”^{3}The densities for mass, momentum, and energy were defined without restriction to a single component single phase system. For fluxes (stress and heat flux), IK stated that “we must impose this restriction” and the flux formulas obtained with their formalism are “only valid for a single component, single phase system.”

^{3}The statistical mechanical theory of transport processes published in a series of 14 papers by Kirkwood and co-workers only considers rigid molecules, ignoring the internal degrees of freedom (DOFs) of the molecules. In the first one of the series, Kirkwood envisioned extension of the formulation to molecules possessing internal degrees of freedom (DOFs).

^{9}

## III. THE CAC FORMALISM

The CAC formalism is an extension of the IK procedure for homogeneous hydrodynamic systems to a two-level structural description of crystals that includes the atomic DOFs within each lattice cell.^{10,11} With such an extension, the CAC formulation is applicable to polyatomic crystalline materials and can be used to solve atomic trajectories or quantify instantaneous fluxes, conferring full characteristics of a concurrent atomistic-continuum method.

### A. Ensemble-averaged local densities

For a polyatomic crystal of $Nl$ unit cells, with each unit cell containing $N\alpha $ atoms, the density of a phase function $A(rk\xi ,vk\xi )$ at point ** r** in the physical space can be expressed as

where $rk\xi $ and $vk\xi $ are the position vector and velocity vector of the ξth atom in the *k*th unit cell, and the symbol ⟨ ⟩ denotes ensemble averaging that involves $NlNa$ fold integrals of the phase function ** A** and a probability distribution function in the phase space.

Recall that the structure and deformation of a crystal are continuous in space at the lattice level but are discrete at the atomic level. It is, therefore, consistent with solid state physics to describe the crystal structure in terms of the positions of lattice cells and the relative positions of atoms inside the lattice cells. For this purpose, we resolve the sum in the definition of the local density in Eq. (3) into contributions from groups of atoms within the same lattice cell and express the locations relative to the unit cell for those atoms, i.e.,

where $rk$ is the position of *k*th unit cell and $\Delta rk\xi $ is the relative position of ξth atom in the unit cell *k*.

The density in Eq. (4) is expressed in terms of the direct product of two deltas, with the two field variables ** x** and

**representing a distinction of large- and small-scale variations of the variable**

*y***in Eq. (3). Note that the ensemble average in Eq. (4) also involves $NlNa$ fold integrals and, although expressed in terms of different phase variables, the values of a local density defined in Eqs. (3) and (4) can be mathematically equivalent. However, different from that in Eq. (2), a conservation equation for the density defined in Eq. (4) relates a conserved quantity to two fluxes, $J\xaf1$ and $J\xaf2$, owing to the use of large- and small-scale variations in the latter, i.e.,**

*r*### B. Space and time-averaged local densities

Evans and Morriss argued that, since conservation laws hold instantaneously, conservation equations and fluxes should be definable without ensemble averaging.^{12} This has been demonstrated recently for monatomic crystals.^{13–16} To derive the conservation equations and flux expressions that hold instantaneously based on the two-level structural description of crystals, we consider a unit cell of volume *V* located at point ** x** in the physical space; within the unit cell, there are $N\alpha $ volume elements, each containing one atom and having volume

*V*satisfying $\u2211\alpha V\alpha =V$. Since the distribution of the unit cells is continuous and homogeneous in space, we may define a local density per unit-cell volume as

_{α}where $a\alpha (x,t)=\u2211kA(rk\alpha ,vk\alpha )\delta V(x\u2212rk)$ is the contribution of the *α*th atom to $a(x,t)$ and can be expressed as

and

Since $\u2211\alpha =1N\alpha \delta V\alpha (y\u2212\Delta rk\xi )=1$, one can see that Eq. (7) satisfies Eq. (6), and $a(x,y,t)=a\alpha (x,t)$ holds for all $y\u2208V\alpha $.

Note that the box functions defined in Eqs. (8) and (9) have a jump discontinuity at the enclosing surfaces *V* and *V _{α}*, respectively. These functions do not have derivatives in the classical sense. It thus requires a generalization of the concept of a function. The corresponding mathematical theory is known as the theory of distributions.

^{17}This theory will be used throughout this work as the mathematical tool to derive the conservation laws for physical quantities described by nondifferentiable functions.

In atomistic simulations, the equation of motion is solved in discrete time steps. Accordingly, local densities should be further averaged over a time-interval. The contribution of the *α*th atom to the mass density *ρ _{α}*, momentum density $p\alpha (p\alpha =\rho \alpha v\alpha )$, and energy density $E\alpha (E\alpha =\rho \alpha e\alpha )$ per unit volume can then be defined as

where $\Phi k\xi $ is the potential energy of atom *kξ*; $\delta \xafV$ denotes the average of *δ _{V}* over a time-interval

*Δt*, i.e., $\delta \xafV(rk\u2212x)=1\Delta t\u222b0\Delta t\delta V[rk(t+\tau )\u2212x]d\tau $. The time evolution of a conserved quantity, i.e., the conservation equation, can then be expressed in terms of a lattice-level flux

**as**

*J*or in terms of atomic-level fluxes *J*_{1} and *J*_{2} as

where ** n** is the outward unit normal vector to the surface element,

*J*_{1}is the part of flux that flows across the enclosing surface ∂

*V*, and

*J*_{2}is the part of flux that flows back and forth inside

*V*and hence only cross the bounding surface ∂

*V*, as shown in Fig. 4.

_{α}It is noted that by integrating Eq. (5) over *V* and invoking the divergence theorem, one obtains an integral equation in the form of Eq. (14). However, Eq. (5) is not valid instantaneously and it requires the fluxes to be differentiable. By contrast, the conservation equation in Eq. (14) holds instantaneously in the distributional sense and the fluxes are not required to be continuous functions.

### C. Conservation equations and fluxes

Averaging a local density over a unit cell of volume V leads to the integral form of conservation laws that equate the rate of the change of conserved quantities in the unit cell to their fluxes across the enclosing surface of the unit cell. Following from Eq. (10), the equation of conservation of mass can be derived in the form of Eq. (14) as

where $v=x\u02d9$, $v\alpha =p\alpha /\rho \alpha $, and $\Delta v\alpha =v\alpha \u2212v$.

The momentum conservation law can also be derived, using the momentum density defined in Eq. (11), from Newton's second law, as

where *t** _{α}* and $\tau \alpha $ are the momentum flux (stress) that cross ∂

*V*and ∂

*V*, respectively, as a result of atomic interaction and motion involving atom

_{α}*α*. The total stress vector

**on a surface element**

*t**A*of area

^{n}*A*and normal

**, located at position**

*n***in the unit cell at**

*y***, from all the atoms can then be expressed as**

*x*where $tpot$ is the potential part of ** t** that measures the interaction force per unit area transmitted across

*A*, i.e.,

^{n}and $tkin$is the kinetic part of ** t** that measures the flow of momentum per unit area and time across

*A*, i.e.,

^{n}Here, $v~k=vk\u2212v$, $\Delta v~k\alpha =\Delta vk\alpha \u2212\Delta v\alpha $, $\u222brl\eta rk\xi \delta \xafAn(\phi \u2212x)d\phi $, defined in Eq. (A3), and $\delta \xafAn(y\u2212\Delta rk\alpha )$, defined in Eq. (A4), represent the potential and kinetic flux as a line-plane intersection problem in space and time, respectively.

Similarly, denoting *q** _{α}* and

*j**as the parts of heat flux that cross ∂*

_{α}*V*and ∂

*V*, respectively, the time rate of change of the energy density

_{α}*E*, can be derived as

_{α}The total heat flux is then given by

where ** e^{i}** (

*i*= 1, 2, 3) are the orthonormal basis at

**,**

*x*and

The link between atomistic and continuum descriptions at this point has fundamentally departed from the IK formalism. As distinct from classical continuum mechanics, this new formulation describes the structure and properties of a crystalline system as a continuous function in ** x** at the lattice level, cf. Figs. 3 and 5, while discrete in

*α*or

**at the atomic level; fluxes (stress and heat flux) are composed of components resulting from the motion and deformation of lattice cells (denoted as**

*y***and**

*t***, respectively) and the rearrangement of atoms within the cells (denoted as**

*q***τ**and

**j**, respectively). Most importantly, instead of linking molecular variables to classical continuum mechanics equations, the CAC formalism leads to a new representation of conservation laws and fluxes.

### D. Temperature

An important difference between atomistic and classical continuum theories is the definition of temperature; in the former, temperature is a derived quantity, while in the latter, temperature is a fundamental state variable. Specifically, the abiding definition of temperature in classical molecular simulations is the kinetic temperature defined by the kinetic theory in terms of the mean kinetic energy of the random motion of atoms, while that in quantum calculations is based on the mean energy of phonons.

The kinetic theory provides a microscopic description of temperature through the classical mechanics of particles. The equipartition theorem of the kinetic theory asserts that the average kinetic energy of a particle in an equilibrium gaseous system with temperature *T* is $12mvk\alpha 2=32kBT$. Away from thermal equilibrium, the particle velocity is usually replaced by $v~k\alpha $, the difference between particle velocity and the velocity field. Such a temperature definition is consistent with the physical picture of an ideal-gas thermometer in which the velocity is measured relative to the co-moving frame of the kinetic thermometer.^{18} The nonequilibrium temperature is expressed^{19} as

This temperature definition, however, has unambiguous meaning only in thermal equilibrium or an approximate local thermal equilibrium condition such as that in a steady state. To elucidate this point, we decompose the kinetic energy density into two parts,

where $k\alpha 1$ is the kinetic energy density represented by the velocity field and $k\alpha 2$ is the kinetic energy due to the velocity difference between particle velocity and velocity field, usually interpreted as the random motion of particles. Although the definition of total kinetic energy is unambiguous, the velocity field can only be uniquely defined when it is zero (for thermal equilibrium) or a constant (for steady state). This is because the velocity field depends on length and time scales, i.e., the averaging, of the linear momentum. Thus, regardless of whether it is in an atomistic or a multiscale simulation, temperature can only be well defined at thermal equilibrium or steady state. It is a concept describing the state of a system rather than the dynamics of an atom.

Differing from the classical concept, the quantum description of temperature is expressed in terms of phonons. There are more phonons at high temperature and fewer phonons at low temperature. The average phonon number is determined by the equilibrium temperature *T* through the Planck constant *ħ* and the Boltzmann constant $kB$. This relationship is called the Bose-Einstein distribution, i.e.,

where $n(\kappa ,\nu )$ is the phonon number for branch *ν* and wave vector ** κ** in thermal equilibrium. The total kinetic energy of an equilibrium system can then be expressed in terms of all the available phonon modes in the systems, i.e.,

Recall that the number of phonon modes is equal to the number of the particle degrees of freedom of the system and that the total kinetic energy per phonon mode in the classical limit is equal to *k _{B}T/*2

*.*The difference between the quantum and classical temperatures can thus be quantified through equating the phonon and classical descriptions of the total kinetic energies,

^{20}i.e.,

For a given quantum temperature, the corresponding classical temperature can be calculated using Eq. (28). Figure 6 presents such calculations for LiF and SrTiO_{3}. As can be seen from Fig. 6, the two temperatures are different at low temperature, but the difference becomes negligible as the temperature increases. At 0 K quantum temperature, LiF and SrTiO_{3} have a zero-point energy equivalent to 237.1 K and 279.7 K classical temperature, respectively, and hence there are corresponding motions at 0 K, i.e., the zero-point vibrations.

## IV. FINITE ELEMENT FORMULATION

As shown in Fig. 5, a local density such as the mass density in a polyatomic system is not a continuous function at the atomic scale. By decomposing the atomic position into the position of the lattice cell, ** x**, and its internal position within the unit cell at

**, we are able to express a density**

*x*

*a**(*

_{α}

*x**,t*) or

**(**

*a***,**

*x, y**t*) at an equilibrium or steady state system as a homogeneous and continuous function in

**. Consequently, the field equations of the conservation laws can be discretized in**

*x***and solved using continuum-based numerical methods such as the finite element (FE) method.**

*x*### A. Temperature in FE formulation

As shown in Eq. (25), the decomposition of the total kinetic energy into $k\alpha 1$ and $k\alpha 2$ depends on the length scale of the velocity field. In particular, in a FE model $k\alpha 1$ is represented by FE nodal velocities, in which the displacement field in an element is interpolated with FE shape functions $S\xi (x)$, i.e.,

where $U\xi \alpha (t)$ is the displacement vector of *α*th atom embedded within the ξth node of the element; the contribution of the *α*th atom to the kinetic energy in an element of volume *V _{e}* is thus given by

The kinetic energy of an element can also be expressed in terms of phonons using the phonon dispersion relations of the FE model.^{21} However, a FE model employing the usual trilinear shape functions cuts off phonons whose wavelengths are smaller than the element size. The motion of the FEs thus only represents a subset of the phonons of the system. This means that the FE mesh, i.e., the FE displacement approximation, uniquely determines the decomposition of total kinetic energy into $k\alpha 1$ and $k\alpha 2$, where $k\alpha 1$ is represented by FE nodal velocities, while $k\alpha 2$ is related to $k\alpha 1$ through the temperature-dependent phonon number $n(\kappa ,\nu )$.

A special case is for systems at high-temperature thermal equilibrium, in which every phonon mode or every particle DOF has the same energy $kBT/2$. The total kinetic energy of an 8-node element containing $nl=Ve/V$ unit cells is thus given by

with

and

### B. Finite element equations

Substituting the conservation equations of mass into the linear momentum equation, we obtain

Denoting the sum of the two surface integrals of kinetic stresses as $f\alpha T(x,t)$, noting that the averaged $f\alpha T(x,t)$ over a long-time duration can be linked to $k\alpha 2$ and hence the local temperature, we have

where $\u2207xT$ and *λ* can be determined from $\u2207x(k\alpha 1+k\alpha 2)$, in which $\u2207xk\alpha 1$ can be calculated based on FE nodal velocities, and $\u2207xk\alpha 2$ can be determined from $\u2207xk\alpha 1$ if we assume local thermal equilibrium and employ the assumption used in the linearized phonon Boltzmann equation, i.e.,$\u2207xnk(x,t)=\u2207xnk0(x,t)$,^{22} where $nk0$ and $nk$ are the equilibrium and nonequilibrium phonon mode distribution for wave vector *k*, respectively.

Equation (34) can be solved using the finite element method, with the weighted residual of Eq. (34) over an element being given by

This is the weak (Galerkin) form of Eq. (34). The integral can be numerically evaluated using Gaussian quadrature, and the equation can be solved for equilibrium or nonequilibrium processes.

For systems at thermal equilibrium,$\u2207xT=0$, $f\alpha T$ is a fluctuating function of time with a zero mean everywhere except at the system's boundary if thermal expansion is constrained. It can thus be modeled as a constant force applied at the boundary. When atomic-scale fluctuations are important to the simulated phenomena such as thermal activation of dislocation motion or phase transition, $f\alpha T$ can be modeled as a periodic or randomly fluctuating force in time with the frequencies of the phonons that are cut off by the FE shape functions.

For nonequilibrium processes, $\u2207xT$ is a constant within an element if the usual linear interpolation is used. $f\alpha T$ can then be modeled as a body force with a constant mean in each element but with periodic or random fluctuation in *t*. It is noted that, to model it as a random force, $f\alpha T$ needs to be partitioned into a fluctuating part and a damping part, with the latter to be determined by and to balance the fluctuation force, according to the fluctuation-dissipation theorem.^{23}

### C. Fluxes in FEA solutions

Similar to that in MD,^{13–15,24} fluxes across a surface element in a CAC model can be computed through calculations of the interaction forces and atomic motions that cross the surface element. Not only can the stress and heat flux at any given surface element be calculated, their spatial distributions can also be plotted as contours for the entire model.

With the usual trilinear shape functions to approximate the displacement field in FE, maximum fluxes occur at the nodes of the elements. Thus, only fluxes at the nodes and a few other locations on the boundaries of the elements are needed for the flux contours, such as the unit cells shown in Fig. 7. In regions of structural disorders or discontinuities, such as those near a dislocation core or crack tip, atomic-scale fluxes can be calculated using Eqs. (18) and (19) and plotted as a contour at atomic resolution.

## V. NUMERICAL EXAMPLES

Currently, there are two public versions of CAC codes, the CAC package in LAMMPS and PyCAC. The general procedure for a CAC simulation can be summarized as follows:

Select an interatomic potential(s) for the system;

Discretize the system based on its morphological length scales, using atomic resolution for atomic-scale structural disorders such as grain boundaries (GBs), while using FEs for single crystal regions with the mesh size determined by, e.g., the smallest possible dislocation spacing;

Relax the system with initial temperature. One may need to systematically change the initial configurations to obtain the dislocation or GB structures that have the minimal energy;

Apply loading and boundary conditions;

Run simulation to solve for and to visualize atomic trajectories, which are either atomically resolved or interpolated through FE nodal displacements, for observations of dynamic processes or emerging phenomena, using software such as Ovito for atomic visualization or Tecplot for FE visualization; and

Calculate fluxes (stress or heat flux) or other transport properties.

The accuracy and efficiency of the CAC method have been tested through one-to-one comparisons with MD in space- and time-resolved simulations of phenomena such as crack initiation and branching,^{26,27} phase transitions,^{28} dislocation nucleation,^{29–32} loop formations^{33–36} interactions with interfaces and other defects,^{37–43} as well as phonon-dislocation interaction^{44,45} and phonon-GB interaction.^{46} In general, for an element containing *N* by *N* by *N* unit cells, the ratio of the number of mathematical operations in CAC to that in MD scales as *O*(*N*^{−3}).^{42} This means that CAC can be very efficient by using large elements, e.g., for polycrystalline materials with micrometer-sized grains or for crystals with low dislocation density, i.e., large dislocation spacing.

To illustrate the applicability of CAC, in this section, we present three sets of prior simulation results. These simulations were performed at low temperature without considering the fluctuations from phonons whose wavelengths are smaller than the element size.

### A. Atomic-scale structures of interfaces

CAC has been employed to study the structure and energy of grain boundaries (GBs) and their dynamic interaction with cracks and dislocations in SrTiO_{3} polycrystals.^{42,43} In these studies, GBs are atomically resolved, while regions away from GBs are coarse-grained. Simulation results show that both the energy and atomic displacements of SrTiO_{3} GBs from CAC agree well with those that from MD if the width of the atomically resolved region for the GBs is larger than 12 Å.^{42} The ghost forces reported at the atomic-continuum interface in many multiscale models are not present in the CAC models.^{42} The atomic-scale GB structures from CAC^{42} are shown in Fig. 8 to compare well with those from experimental results.^{47}

Figure 9 presents CAC simulation results of the misfit dislocation networks in the phase interfaces in heteroepitaxial thin films. CAC is shown to have reproduced the experimentally observed dislocation networks in PbTe/PbSe (001) interfaces.^{48} The dislocation spacing in the system as a function of the epilayer thickness obtained by CAC is also in good agreement with the experimental measurements. In addition, the critical thickness above which misfit dislocations are observed to form in the heteroepitaxial thin film compares very well with experimental observations.^{49}

### B. Fracture and dislocations

Since the only constitutive relation in CAC is the nonlocal force field, and CAC is of an integral form, continuity between FEs is not required. Consequently, nucleation and propagation of dislocations and cracks can be simulated in the FE region via sliding and separation between elements. Simulation results can be output in terms of elements or in terms of atomic trajectories. This makes CAC a convenient tool for the simulation and visualization of the dynamic processes of defect nucleation, interaction, and propagation.

Figure 10 presents CAC simulation results for the loading of a notched solid normal to the notch depth. The simulations have captured the dynamics processes of crack initiation, propagation, and branching. They also reveal that the interactions of the stress waves that emitted from the propagating crack with those that are reflected back by the boundaries of the specimen lead to a sharp increase of the crack speed, which in turn triggers the dynamic instability and hence the branching of cracks in micrometer-scale specimens.^{26} In addition, we have measured the critical stress intensity factor and observed a crack speed near the speed of Rayleigh wave before crack branching.

CAC has also been used to characterize the complex dynamics of fast moving dislocations, including the energy intensities and the wavelengths of phonons emitted from sonic dislocations, as well as the velocity-dependent stress fluctuations around the core of moving dislocations. Figure 11 presents snapshots comparing dislocation motion simulated using CAC with a uniform mesh and that using MD, with the number of DOFs of the CAC model being 1% of that of MD. The dislocation velocity is measured to be ∼2900 m/s in both CAC and MD simulations. Mach cones in a V-shaped pattern of the phonon wave-fronts are observed in the wake of the sonic dislocations. Analysis of simulation results based on a wavelet transform shows that the faster a dislocation is moving, the longer the emitted phonon wavelength.^{29–32}

Figure 12 presents CAC simulation results of the sequential slip transfer of *a*_{0} ̸2 [110](111) dislocations across a Σ3 (111) coherent twin boundary (CTB) in Ni^{56} under ostensibly quasistatic strain rate conditions. Simulation results shows that CAC can serve as an effective tool for the study of slip transfer reactions of sequences of incoming and outgoing long-range dislocation pileups of mixed dislocation characters. Five embedded atom method (EAM) potentials^{51–55} were compared; results demonstrate the uncertainty in computed dislocation-interface reactions associated with the use of a variety of interatomic potentials and suggest that the applicability of dislocation/GB interaction criteria in the literature for slip transfer may have limitations.

### C. Phonon thermal transport

CAC solves for approximate atomic trajectories. Collective motion of atoms, i.e., phonons, and the nucleation and propagation of defects, as well as phonon scattering by defects and other phonons, naturally emerges in the simulation. It thus provides a unified treatment for phonon transport and defect dynamics, as both are represented by the same atomic motion, without the need to assume the nature of phonon transport or scattering, or the mechanism for defect nucleation, propagation, and interactions.

Figure 13 presents CAC simulation results of a heat pulse propagating across GBs.^{46} A phonon representation of heat pulses composed of spatio-temporal Gaussian wave packets obeying the Bose-Einstein distribution^{57} is used to mimic the coherent lattice excitation by ultrashort lasers. As can be seen from Fig. 13, the simulation has captured the following:

The phenomenon of phonon focusing “caustics,”

^{58}manifested as the concentration of energy flux in the sixfold symmetric focusing pattern shown in Fig. 13(a);The Kapitza resistance, shown in Fig. 13(d) as the discontinuities in the kinetic energy at the GBs;

Simultaneous ballistic and diffusive phonon transport across GBs, and a co-existence of coherent-incoherent phonon scattering by the GBs. These are clearly shown in Figs. 13(b)–13(d) as both wavelike transport of the propagating heat pulse and the breaking of the 6-fold symmetry of the phonon focusing pattern due to diffusive transport and incoherent scattering by the GBs; and

Local reconstruction of the GB structures, as shown in Fig. 13(e), as a result of the phonon-GB interaction during the heat pulse propagation.

## VI. SUMMARY AND DISCUSSIONS

We have presented the CAC method for simulation of transport processes in crystalline materials. The formalism extends the Irving-Kirkwood statistical mechanical theory of transport processes for single-phase, single-component molecular systems to a concurrently coupled atomic-continuum description of polyatomic crystals. Conservation equations are derived in terms of instantaneous expressions of cell-averaged quantities using the mathematical theory of distributions.^{59} Fluxes are then obtained that quantify the flow of conserved quantities across the lattice cells as well as those that flow back and forth within the cells, as direct consequences of the local density definition and Newton's second law.

Basic features of the CAC formulation, in comparison with classical continuum mechanics and molecular dynamics, may be summarized as follows:

As in classical continuum mechanics, the CAC governing equations are field equations, and hence can be discretized and solved using the finite element method. However, unlike the situation with continuum mechanics, each element node in CAC represents a unit cell that further contains a group of atoms; the atomic interaction and atomic-level crystal structure are thus fully built into the formulation, with no additional constitutive equations other than an interatomic potential being needed.

By including the internal DOFs of each lattice cell associated with atomic movement in the formulation, CAC can resolve the atomic structure of crystalline materials and reproduce both acoustic and optical phonon modes. It also can measure fluxes at the atomic scale, such as stress near a dislocation core or heat flux across an interface or a defect.

Since CAC is derived in a bottom-up manner from the atomistic, it becomes identical to the underlying atomistic model when discretized at the lattice level. CAC thus has the full applicability of MD in predicting materials behavior without the need for

*a priori*assumption for mechanisms. It can also go beyond MD to scale up for mesoscale simulations by utilizing the continuous structure and deformation in regions where structure disorders or other discontinuities are absent or where atomic resolution is not needed.CAC is thus naturally a concurrent multiscale tool that resolves details to full atomic resolution at regions of interest, while coarse-graining to thousands of atoms per finite element elsewhere; it admits the propagation of dislocations, addresses stacking fault structures in the FE regions, and provides a unified treatment for material microstructures, defects, and phonon transport.

A comparison of representative coarse grained (CG) methods can be found in Ref. 60, and a detailed review of dynamic multiscale methods can be found in Ref. 61. CAC is currently the only concurrently coupled atomistic and continuum method whose continuum description is not based on the classical continuum mechanics equations but a reformulated field equation of conservation laws. As distinct from other CG and multiscale methods, the CAC formulation links and unifies atomistic and continuum description of transport equations. There are two unique capabilities of the CAC formulation:

It can serve as a multiscale theory for concurrent atomistic-continuum simulation of polyatomic crystalline materials under a single set of governing equations, thereby expanding the current atomic-interaction-based simulation capabilities from the nanoscale to the mesoscale crystalline materials.

It provides consistent formulas for calculations of continuum quantities such as stress and heat flux in atomistic or CG atomistic simulations of heterogeneous materials, noting that the majority of flux formulas in the literature or currently implemented in MD software are only applicable to homogenized, single-phase, single-component materials.

Irving and Kirkwood recognized that the classical hydrodynamics is a single-scale theory as a result of a single-scale materials description.^{3} Beyond the single-scale continuum mechanics or hydrodynamics, we have the generalized continuum mechanics (GCM) that advocate and employ a concurrent two-level structural description of materials.^{62–64} Well-established GCM theories include the Micromorphic theory,^{65} the Micropolar theory,^{66} and the Microstructural theory.^{67} Unlike CAC, these GCMs were formulated via a top down approach.^{68} Although the link between molecular dynamics and GCM has been attempted,^{69–72} the popular GCMs are two-level continuum-continuum formulations,^{60} lacking a description for the discrete structure or motion at the atomic scale.

We have also presented three sets of simulation results in this tutorial to illustrate the applicability of CAC in reproducing atomic-scale structures of interfaces, dynamic processes of dislocations and fracture, and phonon transport across GBs. These are low-temperature simulations, in which thermal fluctuations with wavelength smaller than the size of the elements in the FE regions are ignored.

For simulation of defect dynamics at finite temperature, the FE mesh size in CAC is determined by the density of defects as well as by the dominant wavelengths of phonons that interact with the defects. In regions where atomic-scale accuracy is necessary to capture interface and defect reconstruction, full atomic resolution is needed.

For low-temperature transport processes in which phonon-phonon scattering has a negligible effect, the usual FE trilinear shape functions may be replaced by those that combine the trilinear shape functions and Bloch wave functions consisting of phonons whose wavelengths are cut off by the shape function. This allows short-wavelength phonons from an atomically resolved region to traverse a coarse-grained FE region and enter another atomically resolved region^{73} so as to model phonon scattering by dislocations, impurities, GB or phase interfaces with all possible phonons.

For finite-temperature phonon transport when high-frequency short-wavelength phonons are dominated by thermally resistive^{74} or diffusive processes, the surface integrals of the kinetic momentum fluxes, i.e., $f\alpha T$, in the momentum equation can be modeled as a body force with a constant mean but fluctuating in time in terms of the frequencies of the missing phonons. This is intended to reproduce the phonon density of state of the underlying atomistic model, thereby approximating the effect of the thermal fluctuations by the missing phonons. Phonon transport represented by the motion and deformation of the elements, whether it is diffusive, ballistic, or simultaneously ballistic and diffusive, will naturally emerge in the simulation as the consequence of the conservation laws, the interatomic potential, and the initial temperature.

Finally, we would like to note that, while the conservation equations and expressions of fluxes and temperature in the CAC formulation are identical to those in an atomistic model of particles, the FE solutions to the conservation equations are approximate numerical solutions. Consequently, simulated phenomena and calculated fluxes using the FE approximations feature a trade-off of efficiency and accuracy. There are thus coarse-graining errors that may increase with the FE mesh size, depending on the nature of the phenomena and the morphological length scales of the materials. Advanced algorithms that can improve the accuracy and efficiency of the simulation, such as adaptive mesh refinement and advanced multiple time step algorithm, are needed. Also, the current CAC formulation is based on the traditional interatomic potentials with well-defined site energies. Future work will be needed to extend the formulation to quantum-accurate, machine-learned potentials.

## ACKNOWLEDGMENTS

This material is based upon research supported by the U.S. Department of Energy (DOE), Basic Energy Sciences under Award No. DE-SC0006539. The contributions of D.L.M. and Y.C. are partially supported by the National Science Foundation (NSF) under Grant Nos. CMMI 1761553 and 1761512, respectively. We thank Zexi Zheng for generating the data plotted in Fig. 6.

### APPENDIX A: CONSERVATION EQUATIONS FOR MONATOMIC CRYSTALS

The classical formulations of boundary value problems are based on the assumption that their solutions are smooth (i.e., with continuous derivatives) and that the equation is satisfied at every point in the domain of the problems. For problems that involve point mass, point charge, or other discreteness or discontinuities, the smoothness assumption breaks down. The mathematical tool for such problems is the theory of distributions, also known as the theory of generalized functions, developed by Laurent Schwartz in the 1950s.^{17} One generalized function frequently used in science and engineering is the Dirac delta-function introduced by Paul Dirac in his research in quantum mechanics and is employed in most statistical mechanics formulations to link a particle description to a field representation. The theory of distributions enables a rigorous solution to the mathematical difficulties in linking discrete and continuum descriptions that have been bypassed by *ad hoc* constructions or by heuristic arguments.

Appendixes A and B present formal derivations of conservation laws based on the mathematical theory of distributions. A more detailed discussion of the integral form of conservation laws and atomistic formulas for fluxes for monatomic crystals can be found in Refs. 13–16. Due to additional quantities such as temperature discussed in this work, there are notational changes in this paper.

#### 1. Local density definition and the balance equation of linear momentum for monatomic crystals

We define the volumetric density of mass, linear momentum, and energy over a unit cell and a time-interval as

where $\delta \xafV$ is the time step-averaged box function $\delta V$ defined in (8) and is used to serve the role of averaging the IK point-function local density over a unit cell *V* and a time-interval *Δt*.

The derivative of $\delta V$ does not exist everywhere as a classical function, but it exists as a distribution. A distribution *f( x)* is not required to have a value at any point

**, it is defined by its action on a test function**

*x**φ(*that infinitely many times differentiable and has a bounded support. For example, the Dirac delta-function

**x**)*δ(*is defined as (

**x**)*δ, φ) = φ(0).*In this sense, the derivative of

*δ*is the surface distribution that acts on a test function

_{V}*φ*as $(\u2202\u2202t\delta V(x\u2212rk,\phi (x,t))=\u222b\u222b\u25ef\u2202V\u2061\phi vk\u22c5ndS$. Any smooth function

*g(*can be viewed as a distribution whose action on a test function

**x**)*φ(x)*is defined as the integral of the product

*g(*over the whole range of the variable

**x**) φ(**x**)

*x**.*Furthermore, any distribution

*f(*is proved to be the limit of a sequence of smooth functions

**x**)*f*in the sense that the numerical sequence

_{ɛ}*(f*converges to the number

_{ɛ}, φ)*(f, φ)*for any test function as the regularization parameter

*ɛ*approaches 0. In this sense, the derivative of

*δ*can formally be written as a definite integral of a

_{V}*δ*-function, that is, as the limit of the corresponding integral in which the

*δ*-function is replaced by a regularized distribution

*δ*that is a smooth function and converges to

_{ɛ}*δ*as $\epsilon \u21920$ in the distributional sense,

^{59}i.e.,

For the rest of the appendixes, any formal integral of a *δ*-function, like the one in the right-hand side of [Eq. (A2)], is understood as a distribution defined by such a limit. Any relation between distributions obtained in this way is proved to be independent of the choice of regularization.^{59}

The internal force density is also a volume density. An important fact is that the total internal force in a volume element *V* is equal to the total interaction force transmitted across the enclosing surface ∂*V*, which can be described as a line-plane intersection problem. We denote an oriented surface element that is centered at ** x** with normal

**and area**

*n**A*as

*A**; the intersection of the line segment*

^{n}

*L**and surface element*

_{kl}**can be expressed as**

*A*^{n}It is noted that the time-interval average of $\u222bLkl\delta An(x\u2212\phi )d\phi $ over a time step is itself.

The crossing of an atom through a surface element during a time step *Δt* can also be expressed as the intersection of the path of the atom $rk$ and a surface element ** A^{n}** with area

*A*and normal

**as**

*n*Expressing the total internal force in the volume element *V* as the total interaction forces crossing the bounding surface ∂*V*, we obtain the internal force density in terms of site energy $\Phi k$, according to Eq. (A3), as

The time rate of change of the linear momentum within *V* can be derived, as a result of Newton's second law applying to the particles in *V*, using Eqs. (A5) and (A2), as

where the stress vector ** t** on the surface element of normal

**is expressed as a line-plane intersection problem, cf. Eqs. (A3) and (A4), as**

*n*where $Fkl\u225c\u2212(\u2202\Phi l/\u2202rk\u2212\u2202\Phi k/\u2202rl)$ is a measure of the strength of the interaction between atom *k* and atom *l*.

#### 2. Conservation equation of energy for monatomic crystals

Similarly, the time rate of change of the energy density can also be expressed in terms of atomic variables as

The first term in the right-hand side of Eq. (A8) can be expressed as

and the second term, according to Eq. (A2), as

Substituting Eqs. (A9) and (A10) into Eq. (A8), we obtain the integral form of the energy conservation law as

and the heat flux vector in the indicial notation as

### APPENDIX B: CONSERVATION EQUATIONS FOR POLYATOMIC CRYSTALS

#### 1. Balance equation of linear momentum for polyatomic crystals

The rate of change of the density of linear momentum defined in Eq. (11) can be expressed as

The first term in the right-hand side of Eq. (B1) can be derived according to Newton's second law and, in the absence of a body force, it is equal to the internal force density $f\alpha int(x)$. Note that the translational symmetry of potential energies requires that

This leads to the expression of the total force acting on atom *kξ* due to all the other atoms as

Figure 14 shows that the total interaction force on atom *α* in the unit cell at point ** x** can be decomposed into the interaction forces on atom

*α*due to the interactions between atoms in different unit cells and that between atoms in the same unit cell. This then leads to the definitions of surface tractions, $t\alpha pot$ and $\tau \alpha pot$. The internal force density $f\alpha int(x)$ can then be written, according to Eq. (B3), using the line-plane intersection defined in Eq. (A3), as

where

The second term in Eq. (B1) can be expressed according to Eq. (A2), noting that all the cross terms involving the velocity differences vanish over a time-interval average and $\u2211\alpha =1N\alpha \delta V\alpha (y\u2212\Delta rk\xi )=1$, as

and the third term as

with

Substituting Eqs. (B4), (B7), and (B8) into Eq. (B11), we obtain the conservation equation of linear momentum as

where

The total stress vector on a surface element at point ** y** in the unit cell located at

**due to the interaction and motion of all the atoms is then expressed as a sum of potential and kinetic parts as**

*x*where

and

The stress tensor can then be expressed in terms of the stress vector as

#### 2. Conservation equation of energy for polyatomic crystals

The rate of change of the energy density can be written as

With

the first term in the right-hand side of Eq. (B17) can be expressed, using Eq. (B18) and the line-plane intersection theorem, as

The second term can be derived according to the distributional derivative of the box function expressed in Eq. (A2) as

and the third term as

This leads to

The total heat flux is thus expressed as

with

and

Irving and Kirkwood noted that a definition for the site energies of particles is required for definition of energy density as well as for the derivation of a field representation of energy conservation law in terms of particle variables.^{3} The conservation equations and the expressions for fluxes derived in the appendixes are valid for any analytical many-body potentials with well-defined site energies. A detailed comparison and discussion of the fluxes defined as a surface density presented in Appendix A and the fluxes formulas expressed as a volume density can be found in Ref. 13. Specifically, analytical and MD simulation results presented in Ref. 13 show that replacing the Dirac δ with a finite-sized volume weighting function changes the fundamental nature of fluxes as a surface density. Consequently, the dynamic balance between the rate of change of the total momentum or energy in a volume element and the fluxes across the surface boundary cannot be established, leading to the failure of the volume averaged flux formulas in satisfying the momentum and energy conservation laws as well as typical transport boundary conditions. The equivalence between flux formulas expressed as a pair density derived using the kinetic theory^{75,76} and that expressed as a single density derived using statistical mechanics,^{3} as ensemble averaged point functions or as surface-averages, can be found in Ref. 16.

## REFERENCES

*Multiscale Materials Modeling for Nanomechanics*, edited by C. Weinberger and G. Tucker (Springer, Cham, Switzerland) Vol. 245, pp. 261–296.

_{3}symmetric tilt grain boundaries

_{3}Al