Closing the gap between atomic-scale lattice deformations and continuum elasticity

Crystal lattice deformations can be described microscopically by explicitly accounting for the position of atoms or macroscopically by continuum elasticity. In this work, we report on the description of continuous elastic fields derived from an atomistic representation of crystalline structures that also include features typical of the microscopic scale. Analytic expressions for strain components are obtained from the complex amplitudes of the Fourier modes representing periodic lattice positions, which can be generally provided by atomistic modeling or experiments. The magnitude and phase of these amplitudes, together with the continuous description of strains, are able to characterize crystal rotations, lattice deformations, and dislocations. Moreover, combined with the so-called amplitude expansion of the phase-field crystal model, they provide a suitable tool for bridging microscopic to macroscopic scales. We show that the elastic field in presence of defects, as computed from the complex amplitudes, accurately matches continuum elasticity results that include the contribution of dislocation cores. The rotation field in tilted crystals and polycrystalline systems is also addressed, delivering ready-to-use orientational order parameters. This study enables the in-depth analysis of elasticity effects for macro- and mesoscale systems taking microscopic details into account.

Strains and defect-induced deformations have tremendous effects on the macroscopic properties of single and poly-crystalline materials [1]. These effects have fostered a huge variety of studies for more than a century, starting with the first theories describing the elastic field generated by dislocations in solids [2,3].
Deformation of crystal lattices, although involving changes in the positions of atoms, are crucial for understanding the behavior of systems defined on larger length scales [4]. Continuum mechanics, and associated continuous elastic fields, are very useful for describing elastic effects on the mesoscopic and/or macroscopic scales. In this approach, a continuous representation of the displacement of atoms in a lattice with respect to a reference crystal is employed [5]. This description is useful for either relatively simple distortions, as the one given by pure elastic deformations and rotation, or for deformations induced by the presence of dislocations [5][6][7]. Such an approach, indeed, can be exploited to provide in-depth studies of material properties allowing for direct comparisons with experiments and/or a-priori predictions, as, e.g, for plasticity onset in complex heterostructures [8,9] or elasticity effects on material transport mechanisms and morphological evolution [10,11].
For some applications, however, continuum mechanics is not enough as neglecting the description of atoms leads to a crucial loss of information. For instance, this happens when looking at contributions of the dislocation core to elastic fields [12,13] and, in turn, at dislocation nucleation, motion, and reaction. In these cases, in order to describe material properties by elasticity theory, the elastic fields must be described within mesoscale [14,15] or atomistic approaches. Typically, severe restrictions * marco.salvalaglio@tu-dresden.de are present for these methods in the description of long time-and large length-scales.
An attempt to overcome the timescale limits of atomistic approaches, by focusing on diffusive timescales, lead to the development of the so-called phase-field crystal (PFC) method. It focuses on the dimensionless atomic density field difference, n, which is related to the probability of finding an atom in a certain position in space filtering out vibrations on lattice sites [16][17][18]. It provides good descriptions of elasticity [19] and dislocation dynamics [20] even if it usually requires fine spatial discretizations. This latter limitation is overcome by the complex amplitude expansion (APFC) [21][22][23] of the PFC model for which both long time scales and large length scales can be examined. It consists of a coarse-grained representation of the density n that is expressed by the sum of Fourier modes representing specific lattice symmetries. The slowly-oscillating complex amplitudes of these modes, η j , are then the variables used to characterize the crystalline lattice. APFC has proved useful in the advanced modeling of materials as illustrated in studies of compositional domains [24], binary alloys [25], dislocation dynamics [26], morphology and motion of dislocation networks at grain boundaries [27] and the generic tuning of materials properties [28][29][30]. However, the basic concept of APFC, namely the coarse-graining of an explicit lattice representation by focusing on the complex coefficients of Fourier modes, can be readily applied to any atomistic description as obtained, e.g, from theoretical modeling, atomistic simulations, or experimental imaging.
The first aim of this work is to show how to exploit this representation to derive expressions for continuous elastic fields independent of lattice symmetry and system dimensionality from the complex amplitudes. In practice, we describe how to reconstruct strain and rotation fields from atomic density provided that the complex amplitudes functions of the corresponding Fourier modes are known. The procedure is illustrated by exploiting the results of numerical simulations for strained/tilted crystalline and poly-crystalline systems. To this purpose, we numerically solve the equations of the APFC model. This choice is related to the second general aim of this work. Indeed, by the analysis of the strain field, we also aim at a further assessment of the APFC approach as a reliable coarse-grained method to account for microscopic effects while describing extended defects such as dislocation and grain boundaries. The elastic field in presence of defects, derived from complex amplitudes as computed by APFC simulations, reproduce predictions of continuum mechanics. However, the results also indicate some deviation from continuum mechanics due to the effect of atomistic structure at the dislocation cores, as it has also been observed in other atomistic-modeling approaches [12,13].

I. AMPLITUDE EXPANSION OF PHASE-FIELD CRYSTAL MODEL
The PFC model accounts for the lattice structure by means of a continuous periodic field n describing the dimensionless atomic probability density and/or the average of atoms position over the vibrational timescale [16][17][18]. It is based on a free-energy functional, F n , that describes a first-order transition between a disordered/liquid phase, where n is constant, and an ordered/crystalline phase, where n is periodic. F n can be written, where ∆B 0 , B x 0 , v and t are parameters as in Ref. [31]. In the crystalline state n can be generally approximated as sum of plane waves as with n 0 (x) the average density, η j (x) the amplitude of each plane wave and k j the reciprocal space vector representing a specific crystal symmetry (see Supp. Info. S1).
In the so-called amplitude expansion of the PFC model (APFC) [21][22][23], the amplitudes are assumed to vary on a length scale much larger than the atomic spacing and are the variables used to describe a given crystalline system. Lattice symmetries are described by means of a fixed set of vectors k j . The complex amplitudes functions η j (x) allow for distortions and rotations of the crystal lattice with respect to a reference state accounted for by k j vec-tors. The free energy, expressed in terms of η j 's, reads with G j ≡ ∇ 2 +2ik j ·∇ and corresponds to a complex polynomial of η j and η * j and depends on the specific crystalline symmetry as reported in Ref. [29] for triangular, body-centered cubic (bcc) and face-centered cubic (fcc) lattices. The evolution laws for η j 's read Simulations reported in the following are performed exploiting the finite element toolbox AMDiS [32,33] with a semi-implicit integration scheme as reported in Ref. [29]. Periodic boundary conditions (PBC) are used for all the boundaries of the simulation domains. To describe a crystalline phase, the parameters entering the free energy are set to favor the crystal phase as follows: B x = 0.98, v = 1/3, t = 1/2 and ∆B = 0.02 [25,27,29].

II. STRAIN AND ROTATION FIELDS FROM COMPLEX AMPLITUDES
The amplitudes functions describing a deformation of the crystal, u, with respect to the reference state corresponding to a specific set of k j entering (2), can be computed as which defines N independent equations. φ j are the real values describing the j-th amplitudes in a relaxed, unrotated crystal. They can be computed by minimization of the energy functional in (3) assuming constant, real amplitudes for each different length of k j . (5) can be rewritten as with ϕ j = arg(η j ) = arctan [Im(η j )/Re(η j )]. In order to determine the components of deformation fields u i from amplitudes, (6) must be inverted, resulting in a system of d equations with d the dimensionality of the system. In 2D, by selecting d = 2 amplitudes labelled by generic indexes l and m, with the two components of the displacement field obtained by index permutations on the group (i, j) = (x, y) and ij the 2D Levi-Civita symbol. In 3D, by selecting d = 3 amplitudes labelled by generic indexes l, m, and n, u i ≡ u 3D i we obtain with the three components of the displacement field obtained by index permutations on the group (i, j, κ) = (x, y, z). Amplitudes must be chosen in order to have a non vanishing denominator of the prefactor entering (8).
Without loss of generality we fix here l = 1, m = 2 and n = 3. For small deformations, the strain tensor ε can be written ε = (1/2)[∇u + (∇u) T ]. The strain components can then be explicitly computed from η j by means of spatial derivatives of (7) (d = 2) or (8) (d = 3). This leads to ∂ϕ j /∂x i terms. Notice that, ϕ j are inherently discontinuous due to their functional form. However, η j are, by definition, continuous complex functions in both their real and imaginary part and Since |η j | 2 = φ 2 0 > 0 almost everywhere in the crystal phase, the terms ∂ϕ j /∂x i and then ε can be readily computed. |η j | 2 only vanishes exactly at the dislocation core position for some indices, j, consistent with continuum elasticity theory. By exploiting the components of the displacement field, rotations of the crystal structure with respect to the reference orientation can also be evaluated as ω = ∇ × u. ω ij corresponds to the rotational angle in the x i -x j -plane. For d = 2, ω ij ≡ ω, describes rotation in the two-dimensional domain. For the sake of readability, the expressions of ε and ω for d = 2, 3, as functions of ∂ϕ j /∂x i , are reported in the Supp. Info. S2.

III. DEFORMATIONS INDUCED BY DISLOCATIONS
Let us consider pairs of dislocations in a 2D triangular lattice that form at the interface between layers of different atomic spacing in order to accommodate a misfit strain. As performed in Ref. [29], the dislocations can be described by APFC by setting an initial condition for η j in order to reproduce opposite deformations (namely initial strains) ± ε. This can be done using (5) with in-plane displacement field u(r) defined by u x = ±(a tri /L x )x and u y = 0, where a tri is the distance between maxima of the density as in (2) for triangular symmetry and L x is the size of the computational domain, both evaluated along the x axis. These opposite displacements are imposed as illustrated in Fig. 1(a). The relaxation of such an initial condition until the defect shape is stationary in a square system of linear dimension L i = 40a tri is illustrated in Fig. 1(a). A 2 is shown on the left and the reconstructed density n from (2) is shown on the right. Two pairs of dislocations form within the simulation domain where A 2 is constant in the solid and decreases at defects. Although the APFC approach does not allow for the exact representation of atoms at the dislocation-core, the lattice distortion is well described as shown by the illustration of atomic planes by solid lines, highlighting the presence of one specific defect. The deformation of the lattice can be quantified at the atomic level by evaluating the Burgers vector, b, as shown in Fig. 1(b). Notice that it corresponds to a lattice spacing inx direction, i.e. |b| = 4π/ √ 3. The continuous strain fields can be readily computed by exploiting the approach of Sec. II. ε xx is shown in Fig. 1(c) while, for the sake of completeness, other components are shown in the Supp. Info. S3.
The elastic field as computed from the complex am-plitudes matches very well with continuum elasticity as shown in Fig. 1(d)-(g). The latter can be computed provided that the dislocation character, namely edge, screw or mixed [7], as well as b are known. As detailed in the Supp. Info. S4, the defects in Fig. 1(c) can be modeled as a 2D array of edge dislocations and the strain field according to continuum elasticity (CE) ε ce can be computed as superposition of the elastic field of single dislocations. Figs. 1(d)-(g) show ε ce by solid red lines along 1 , 2 , and 3 as defined in Fig. 1(c). The strain components computed from η j 's are shown by black circles. An almost perfect agreement is found far away from dislocation cores as observed along 2 for ε xx . The same holds for the other cases except for region close to the dislocation cores, namely at l i ∼ ±10a tri where, however, a continuous description of lattice deformations is not well-posed. Deeper insight into the atomistic nature of the dislocation core can be obtained by focusing on the analysis of the elastic field as in Refs. [12,13]. Therein, the components of the strain field in the presence of a straight dislocation, ε d ij , are decomposed as ε d ij = ε ff ij + ε core ij , with ε ff ij a far-field predictor of the elastic field and ε core ij a correction due to the core effects. A decay ε ff ij ∼ r −1 , with r the distance from the dislocation core, is expected in agreement with continuum elasticity theory. The dislocation core effects are characterized by a decay ε core ij ∼ r −2 independently on the lattice symmetry, except for some high-symmetry nominal position of the core where ε core ij ∼ r p with p < −2 [12,13]. ε ce ij , corresponds to an explicit expression of ε ff ij . The contribution of the dislocation core can then be analysed by evaluating ε core ij = ε apfc ij − ε ce ij , with ε apfc ij the elastic field components computed from η j . The scaling of these quantities is illustrated in Fig. 2 where a system as in Fig. 1 is considered with L i = 320a tri . ε apfc xx (red squared) and ε ce xx (green triangles) are shown along the line 1 (see Fig. 1(c)) with r = 0 the position of one of the defects. In the large r limit the strain shows a decay that is very close to ∼ r −1 , but some small deviations are observed. By evaluating the difference, ε apfc xx − ε ce xx , (blue circles) we observe a faster decay that scales as ∼ r q with q = −2.05 ± 0.15 recalling the more localized correction that can be ascribed to the dislocation core. We can then conclude that the elastic field description proposed here not only matches the description of continuum elasticity far from dislocation cores but also includes a correction typical of refined descriptions including atomic-scale effects.
So far only a 2D system has been discussed. However, the evaluation of elastic field can be readily provided also in 3D regardless the lattice symmetry (see Supp. Info. S2). To this purpose we consider a configuration with peculiar 3D features as shown in Fig. 3. An initial condition mimicking layers with opposite, biaxial strain along thex andŷ direction is considered. η j 's are initialized as in Fig. 1 with u y = u x [29] and, without loss of generality we select a crystal having bcc lattice symmetry (by setting k j accordingly, see Supp. Info. S1). The size of the domain is set to L i = 60a bcc , with a bcc the lattice constant for bcc arrangements. Interfaces between layers are (001) planes. APFC simulations account for the formation of networks of dislocations at the interface as shown in Fig. 3(a). In this figure, the grey structure corresponds to the region where A 2 < 0.8 max(A 2 ), i.e. to the defects as it is the region where A 2 significantly decreases. The initial interfaces between layers are higlighted by blue (top) and green (bottom) planes. The atomic structure at the defects as constructed by (2) and the Burgers vector (here |b| = 2 √ 2π) are illustrated in Fig. 3(b) in a small 2D region around the defect lying on the yz-plane Π highlighted in Fig. 3(a). The strain field computed from amplitudes is illustrated by means of ε yz in Fig. 3(c). By accounting for the parallel dislocations forming along both the in-plane directions, the elastic field can be approximated by continuum elasticity. A comparison of the ε yz component as obtained from amplitudes (APFC) computed by the simulations illustrated in Fig. 3(a) and from continuum elasticity (CE) by adapting the equations of Supp. Info. S2 (accounting then for two sets of dislocations oriented alongx and y directions and having perpendicular Burgers vector as in Fig. 3(b)) is reported in Fig. 3(d) showing a general agreement similar to Fig. 1.

IV. LATTICE ROTATIONS AND POLYCRYSTALLINE SYSTEMS
We focus here on the analysis of deformations and rotations in polycrystalline systems, which involve the evolution of grain boundaries (GB) and their influence on the strain fields.
We study first a simple system made of a straight GB forming between two crystals with a symmetric tilt. In particular, a rectangular domain, L x × L y withx = [10]

andŷ = [01]
, is considered with a straight vertical GB at the center. The relative tilt angle between the two crystals, namely 2θ, is set by initializing the η j functions as with δk j (θ) = k j · R(θ) − k j and R(θ) the counterclockwise rotation matrix. A ±θ tilt is imposed for the left and the right part of the simulation domain respectively (see also Fig. 4(a)). By using PBC a GB with infinite extension is considered. A second GB is also expected, that is shared between the left and right periodic boundary of the simulation domain. L x can be chosen arbitrarily while L y has to match the periodicity of amplitudes alongŷ. Specific details about the simulation setup can be found in Ref. [29].  Fig. 4(a) illustrates the result of relaxing the aforementioned initial condition by means of the A 2 field in a small region at the GB, formed by an array of dislocation, showing three defects as minima of A 2 . ω is reported in Fig. 4(b). A constant value of this field is obtained when moving away from the GB, while a modulation at the defects and at the GB is observed, reproducing the effect of the dislocation on the local orientation of the crystal lattice. The values at which ω saturates in the crystals correspond to ±θ/2 imposed in the initial condition. A measure of the tilt angle is then obtained by exploiting the deformation field. It should be noted that it is not possible to directly obtain the tilt angle by inverting equations (10) due to the phase term. Fig. 4(c) also illustrates the strain component ε xy , calculated from amplitudes η j as discussed in the previous section. In this instance a significant superposition of the strain lobes is obtained, due to the proximity of the dislocations.
The importance of extracting the rotation of grains by means of a scalar field defined everywhere becomes more clear when looking at poly-crystalline systems, as illustrated in Fig. 5. Therein a system with L x × L y = 320a tri × 320a tri is considered. Here the initial condition contains ten crystal seeds, namely circular regions with |η j | > 0 with random positions, random radii ranging from 10π to 20π and integer values of θ ∈ [−5 • , 5 • ]. According to the parameters of the free energy, the crystal phase has the lowest energy and the crystals grow and merge forming a complex network of small-angle GBs made of dislocations. Afterward, these GBs evolve resulting in shrinkage of grains and annihilation of dislocations. This is illustrated in details in Fig. 5 by means of A 2 , ω, and the deformation field in terms of ε xy starting from the merged crystal. This APFC simulation reproduces the typical features of the evolution in polycrystalline systems. In addition, evaluation of deformation fields does not require any data post-processing and solving for η j 's would also allow for the simulation of larger system than atomistic approaches [34]. Faceted GBs are obtained (here in 2D corresponding to closed polygonal chain). Moreover, the motion of dislocations leading to the shrinkage of grains occurs along preferred directions, revealing the accurate description achieved by APFC despite its coarse-grained nature. ω is nearly constant within the grains and varies at the GBs, thus it can be used to identify single grains as they are characterized by different tilts. Moreover, ω accounts for the contribution of single dislocations. Interestingly, the rotation of grains is known to increase when shrinking due to the deformation induced by dislocations which get closer and closer. Values of ω are actually found to increase when grains shrink. This is in agreement with 2D atomistic calculations showing an increase of the interface energy and a rotation of the grains during similar processes [19,35,36]. A three-dimensional system consisting of a single spherical crystal that is rotated with respect to the larger single crystal that it is embedded in is shown in Fig. 6. L i = 50a fcc with a fcc the lattice constant for fcc arrangements. A rotation of 5 • about the [111] direction is set. The initial condition for η j is set by exploiting (10) (see also Ref. [27]). The resulting dislocation network after a first relaxation phase is shown in Fig. 6(a). As in Fig. 3(a), the grey structure corresponds to the dislocations. A closed network of defects forms, having the peculiar structure of defects at twisted grain boundaries when the normal of the spherical inclusion approaches the rotational axis. As expected by classical grain growth theory [37] and 3D simulations [27] this defect network evolves leading to the shrinking of the rotated inclusion. The extraction of the rotation field on the yz-plane Π highlighted in Fig. 6(a) is illustrated in Fig. 6(b). Different stages during the shrinkage are reported. ω yz is zero outside the inclusion and almost constant inside. As expected and observed in Fig. 5(a), a short wavelength modulation is observed at defects. ω yz increases during the shrinkage of the grain pointing out a rotation of the grain [19,35,36]. The corresponding ε yz field is shown in Fig. 5(c).

V. CONCLUSIONS
A continuous description of lattice deformations exploiting the complex amplitudes of Fourier modes representing the periodicity of crystal lattices has been derived. With this approach, deformation fields can be readily computed from an atomistic representation of crystals without any ad-hoc post-processing procedure and independently of system dimensionality and lattice symmetry.
This framework has been shown to achieve its full potential combined with the APFC approach. Other than trivially providing a direct access to amplitudes and allowing for describing either strain or rotated crystals, the APFC model easily allows for large length-scale simulations approaching the ones typical of continuum theories, still retaining essential microscopic features. This is further demonstrated here by showing that APFC encodes dislocation-core contributions to the elastic fields, which are missing in continuum elasticity. In addition, future developments and optimizations of the numerical methods adopted here will further push the capabilities of such an approach. Therefore, a bridging scale theoretical framework for the investigation of crystalline structures, coupling microscopic details with macroscopic properties, is here disclosed.
It is also worth mentioning that the evaluation of deformations and of the local orientation of crystals presented here paves the way to plenty of applications such as microstructural analysis and grain recognition, optimization of numerical method (e.g., orientation-based meshing criterions for APFC) and multiscale modeling of orientation-dependent effects such as, e.g., the coupling of microstructure evolution to magnetic fields.

Supporting Information for:
Closing the gap between atomic-scale lattice deformations and continuum elasticity

S1. RECIPROCAL SPACE VECTORS FOR DIFFERENT SYMMETRIES
In the APFC approach the crystal lattice symmetry is fixed by considering a proper set of N vectors k j to be used in (2) and in the following. For triangular symmetry N = 3 and reciprocal-space vectors are with k 0 = 1. For bcc symmetry N = 6 and the reciprocal-space vectors are with k 0 = √ 2/2. For fcc symmetry N = 7 and the reciprocal-space vectors are with k 0 = √ 3/3.

S2. STRAIN-FIELD AND ROTATION-FIELD COMPONENTS FROM AMPLITUDES
In Sec. II the derivation of the displacement field u and, in turn, of the strain ε and rotation ω tensors is discussed. The result is that they are functions of k j vectors and ∂ϕ j /∂x i with ϕ j = arctan [Im(η j )/Re(η j )].
For 2D systems, the components of ε are given by For 3D systems, ε ij are given by

S3. 2D MULTILAYER: STRAIN-FIELD COMPONENTS
FIG. S1: Strain field components for the 2D multilayer system illustrated in Fig. 1 in the main text. They are computed from ηj with l = 1 and m = 2.
In Sec. III (Fig. 1), the strain field is illustrated as a two-dimensional color map in terms of ε xx only, while comparisons with the elastic field computed from continuous elasticity (see also Sec. S4) are shown for all the components of ε. Fig. S1 shows all the components of such a system. Moreover, it also shows the directions i along which the elastic fields are compared in Fig. 1.

S4. DISLOCATION ELASTIC FIELD FROM CONTINUUM ELASTICITY THEORY
The elastic field of a dislocation can be described by continuum elasticity theory, provided that the dislocation character, namely edge, screw or mixed [i], as well as the Burgers vector b are known. For d = 2, the Burgers vector as in Fig. 1(c) is actually intrinsically perpendicular to the ideal dislocation line, that is, in turn, perpendicular to the two-dimensional domain. Therefore, dislocations observed therein have edge character. The analytic stress field σ d of an edge dislocation with core in x c = (0, 0) and b aligned to the x axis is here considered [i]. In order to remove the singularity inherently present at the dislocation core, the regularization introduced in Ref. [ii] is considered (and here reported in the assumption of infinite, straight dislocation). The stress field components read σ xx /σ 0 = −y 3ζ 2 + 3x 2 + y 2 , σ yy /σ 0 = −y ζ 2 − x 2 + y 2 , σ zz /σ 0 = −2νy 2ζ 2 + x 2 + y 2 , σ xy /σ 0 = x ζ 2 + x 2 − y 2 σ xz = σ yz = 0 E the Young modulus, and ν the Poisson ratio. ζ is the parameter controlling the regularization of the elastic field at the core, here arbitrarily set to ζ = |b|/2. The system illustrated in Fig. 1 consists of a portion of an infinite 2D array of dislocations, as imposed by PBC, with the same Burgers vector when moving along the interface between layers, and opposite when moving in the direction perpendicular to the interface (see also the orientation of dislocations in Fig. 1). The total stress field, σ tot is then obtained by superposing the field σ d originated from each defect in such an array, accounting for the proper orientation of b and a proper shift of x c . The strain field is computed by exploiting Hooke's law, i.e., ε ce = [(1 + ν)/E]σ tot − (ν/E)tr(σ tot )I (S10)