Analysis of Boundary Conditions for Crystal Defect Atomistic Simulations

Numerical simulations of crystal defects are necessarily restricted to finite computational domains, supplying artificial boundary conditions that emulate the effect of embedding the defect in an effectively infinite crystalline environment. This work develops a rigorous framework within which the accuracy of different types of boundary conditions can be precisely assessed. We formulate the equilibration of crystal defects as variational problems in a discrete energy space and establish qualitatively sharp regularity estimates for minimisers. Using this foundation we then present rigorous error estimates for (i) a truncation method (Dirichlet boundary conditions), (ii) periodic boundary conditions, (iii) boundary conditions from linear elasticity, and (iv) boundary conditions from nonlinear elasticity. Numerical results confirm the sharpness of the analysis.


Introduction
Determining the geometry and energies of defects in crystalline solids is a key problem of computational materials science [47,Ch. 6]. Defects generally distort the host lattice, thus generating long-ranging elastic fields. Since practical schemes necessarily work in small computational domains (for example, "supercells") they cannot explicitly resolve these fields but must employ artificial boundary conditions (periodic boundary conditions appear to be the most common). To assess the accuracy and in particular the cell size effects of such simulations, numerous formal results, numerical explorations, or results for linearised problems can be found in the literature; see for example [3,8,17,27] and references therein for a small representative sample.
The novelty of the present work is that we rigorously establish explicit convergence rates in terms of computational cell size, taking into account the long-ranged elastic fields. Our framework encompasses both point defects and straight dislocation lines. Related results in a PDE context have recently been developed in [5].
The second motivation for our work is the analysis of multiscale methods. Several multiscale methods have been proposed to accelerate crystal defect computations (for example atomistic/continuum coupling [25,30] or QM/MM [4]), and our framework provides a natural set of benchmark problems and a comprehensive analytical substructure for these methods to assess their relative accuracy and efficiency. In particular, it provides a machinery for the optimisation of the (non-trivial) set of approximation parameters in multiscale schemes.
The mathematical analysis of crystalline defects has traditionally focused on dislocations [1,2,18,20] and on electronic structure models [10,11]; however, see [9] for a comprehensive recent review focused on point defects. The results in the present work, in particular the decay estimates on elastic fields, also have a bearing on this literature since they can be used to establish finer information about equilibrium configurations; see for example, [21].

Outline
Our approach consists in placing the defect in an infinite crystalline environment, for simplicity say Z d , where d ∈ {2, 3}. Let w : Z d → R m be the unknown displacement of the crystal, then we decompose w = u 0 + u, where the predictor u 0 specifies the far-field boundary condition through the requirement that the corrector u belongs to a discrete energy spaceẆ 1,2 (a canonical discrete variant oḟ H 1 (R d )). We then formulate the condition for w to be an equilibrium configuration asū ∈ arg min E (u) u ∈Ẇ 1,2 , where E (u) is the energy difference between the total displacement w = u 0 + u and the predictor u 0 . It is crucial that u 0 is an "approximate equilibrium" in the far-field, expressed through the requirement that δE (0) ∈ (Ẇ 1,2 ) * . If this condition fails, then inf{E (u) | u ∈Ẇ 1,2 } = −∞. For this reason, we think of u 0 as a predictor and u as a corrector. For the case of dislocations, the choice of u 0 is non-trivial, as the "naive" linear elasticity predictor does not take lattice symmetries correctly into account. In Section 3.1 we present a new construction that remedies this issue.
We will not be concerned with the existence of solutions to (1), which is a difficult problem even for the simplest classes of defects. Uniqueness can never be expected for realistic interatomic potentials.
Assuming that a solution to (1) does exist, we may then analyze its "regularity". More precisely, under a natural stability assumption we estimate the rate at which u (and its discrete gradients of arbitrary order) approach zero, for example, where Du( ) is a finite difference gradient centered at ∈ Z d , r = 0 for point defects and r = 1 for straight dislocation lines.
These regularity estimates then allow us to establish various approximation results. For example, we can estimate the error committed by projecting an infinite lattice displacement field u to a finite domain by truncation. This motivates the formulation of a Galerkin-type approximation scheme for (1) (see Sections 2.3 and 3.4) This is a finite dimensional optimisation problem with dim(Ẇ 1,2 N ) ≈ N , and our framework yields a straightforward proof of: supposeū is a strongly stable (cf. (9)) solution to (1) then, for N sufficiently large, there exists a solutionū N to (3) such that where r = 0 for point defects, r = 1/2 for straight dislocation lines, and ū − u N Ẇ 1,2 is a natural discrete energy norm. Note that N is directly related to the computational cost of solving (3). It is interesting to note that the rate N −1/2 is generic; that is, it is independent of any details of the particular defect. We prove a similar error estimate for periodic boundary conditions in Section 2.4.
In Sections 2.5, 2.6, 3.6, 3.7 we then consider two types of concurrent (or, selfconsistent) boundary conditions that use elasticity models to improve the far field corrector. In these approximate models, we solve a far-field problem concurrently with the atomistic core problem in order to improve the boundary conditions placed on the atomistic core. We also mention that, based on the approach advocated in the present paper, improved a/c schemes with superior optimal convergence rates have recently been developed in [24,36].
Our numerical experiments in Sections 2.7, 3.8 mostly confirm that our analytical predictions are sharp, however, in a preprint [16] we also show cases where they do not capture the full complexity of the convergence behaviour.
Finally, we note that for the sake of brevity this paper only contains the proofs for the decay estimates (2), but not the proofs of the approximation results. Complete proofs, as well as additional numerical tests and some clarifying remarks are contained in the preprint [16].
Restrictions Our analysis is restricted to static equilibria under classical interatomic interaction with finite interaction range. We see no obstacle to include Lennard-Jones type interactions, but this would require finer estimates and a more complex notation. However, we explicitly exclude Coulomb interactions or any electronic structure model and hence also charged defects (see, for example, [9][10][11]17,27]).
As a reference atomistic structure we admit single-species Bravais lattices. Generalising to multi-lattices may require substantial additional work.
We exclude curved line defects, grain or phase boundaries, surfaces or cracks. Moreover, we exclude the case of multiple or infinitely many defects. We hope, however, that our new analytical results on single defects will aid future studies of this setting; for example, see [21] for an analysis of multiple screw dislocations based on the present work.
Notation Notation is introduced throughout the article. Symbols that are used across sections are listed in Appendix B. We only briefly remark on some generic points. The symbol ·, · denotes an abstract duality pairing between a Banach space and its dual. The symbol | · | normally denotes the Euclidean or Frobenius norm, while · denotes an operator norm. C is a generic positive constant that may change from one line of an estimate to the next. C will always remain independent of approximation parameters (such as domain size), of lattice position or of test functions. However, it may depend on the interatomic potential or some fixed displacement or deformation field (for example, on the solution). The dependencies of C will normally be clear from the context, or stated explicitly. To improve readability, we will sometimes replace C with .
For a differentiable function f , ∇ f denotes the Jacobi matrix and ∇ r f = ∇ f ·r the directional derivative. The first and second variations of a functional E ∈ C 2 (X ) are denoted, respectively, by δ E(u), v and If Λ ⊂ R d is a discrete set and u : Λ → R m , ∈ Λ and + ρ ∈ Λ, then we define the finite difference D ρ u( We will normally specify a specific stencil R associated with a site and define Du( u denotes a jth order derivative, and D j u defined recursively by D j u := D D j−1 u denotes the jth order collection of derivatives.

Atomistic Model
We formulate a model for a point defect embedded in a homogeneous lattice. To simplify the presentation, we admit only a finite interaction radius (in reference coordinates) and a smooth interatomic potential. Both are easily lifted, but introduce non-essential technical complications.
Let d ∈ {2, 3} and A ∈ R d×d nonsingular, defining a Bravais lattice AZ d . The reference configuration for the defect is a set Λ ⊂ R d such that, for some For analytical purposes it is convenient to assume the existence of a background mesh, that is, a regular partition T Λ of R d into triangles if d = 2 and tetrahedra if d = 3 whose nodes are the reference sites Λ, and which is homogeneous in For each u : Λ → R m we denote its continuous and piecewise affine interpolant with respect to T Λ by I u. Identifying u = I u we can define the (piecewise constant) gradient ∇u := ∇ I u : R d → R m×d and the spaces of compact and finite-energy displacements, respectively, by . The grey lines indicate the interaction bonds, R , between atoms, for a nearest-neighbour mode, as well as the auxiliary triangulation T Λ It is easy to see [32,34] thatẆ c is dense inẆ 1,2 in the sense that, if u ∈Ẇ 1,2 , then there exist u j ∈Ẇ c such that ∇u j → ∇u strongly in L 2 . Each atom ∈ Λ may interact with a neighbourhood defined by the set of lattice vectors R ⊂ (B r cut ∩(Λ− ))\{0}, where r cut > 0, and we let Du( ) := D R u( ). We assume without loss of generality that if ( , + ρ) is an edge of T Λ , then ρ ∈ R .
This assumption implies, in particular, that ∇u L 2 ≈ |Du| p 2 for any p ∈ ; that it, V should be understood as a site energy difference.) Then the energy functional for compact displacements is given by for u ∈Ẇ c .
We assume throughout, that V is homogeneous outside the defect core, that is, R ≡ R and V ≡ V for all | | R def , and it is point symmetric, Without loss of generality, we also assume that Under these assumptions we can extend the definition of E toẆ 1,2 .
One now uses the fact that the summand in the first group scales quadratically, while δE (0) is a bounded linear functional. The details are presented in Section 5.2.
Lemma 1 implies that the atomistic minimisation problem is meaningful: we seekū where arg min denotes the set of local minimizers. We are not concerned with the existence of solutions to (8), but take the view that this is a property of the lattice and the interatomic potential. We shall assume the existence of a strongly stable equilibriumū ∈Ẇ 1,2 , by which we mean that δE (ū) = 0 and there exists c 0 > 0 such that Since E ∈ C k (Ẇ 1,2 ) and k 2 it is clear that a strongly stable equilibrium is also a solution to (8) (but not vice-versa).

Remark 1.
In [21], (9) is proven rigorously for an anti-plane screw dislocation, under restrictive assumptions on the interatomic potential. We cannot see how one might in general prove such a result, however, in all numerical experiments that we have undertaken to date we observe it a posteriori.

Regularity
Our approximation error results in subsequent sections require estimates on the decay of the elastic fields away from the defect core. These results do not require strong stability of solutions, but only stability of the homogeneous lattice, If (9) holds for any u ∈Ẇ 1,2 , then (10) holds with c A c 0 ; see [16,Section B.2].
Theorem 1. Suppose k 3, that the lattice is stable (10), and that u ∈Ẇ 1,2 is a critical point, δE (u) = 0. Then there exist constants C > 0, u ∞ ∈ R m such that, for 1 j k − 2, and for | | sufficiently large, Proof. The proof for the cases j = 0, 1 is given in Section 6.3. The proof for the case j > 1 is given in Section 6.4.
In what follows we assume k 4, although some results are still true with k = 3.

Clamped Boundary Conditions
The simplest computational scheme to approximately solve (8) is to project the problem to a finite-dimensional subspace. Due to the decay estimates (11) it is reasonable to expect that simply truncating to a finite domain yields a convergent approximation scheme.
We choose a computational domain and solve the finite-dimensional optimisation problem Since dim(Ẇ 0 (Ω R )) < ∞, (13) is computable. Moreover, since it is a pure Galerkin projection of (8) it is relatively straightforward to prove an error estimate.

Theorem 2.
Letū be a strongly stable solution to (8), then there exist C, R 0 > 0 such that, for all R R 0 there exists a strongly stable solutionū 0 R of (13) satisfying Proof (Idea of proof). We construct a truncation operator T R : . Since δE and δ 2 E are continuous it follows that δ 2 E (T Rū ) is positive definite for sufficiently large R and that as R → ∞. The inverse function theorem (IFT) yields the existence of a solution u 0 R , for sufficiently large R, satisfying and consequently also ∇ū 0 R − ∇ū L 2 C ∇T Rū − ∇ū L 2 . Finally, the regularity estimate of Theorem 1 yields the stated rate in terms of R. The complete proof is detailed in the preprint [16].

Remark 2. (Computational cost)
In addition to the assumptions of Theorem 2, assume that R ≈ N 1/d , which is a shape regularity condition for Ω R , then the error estimate (14) reads In particular, if (13) can be solved with linear computational cost, then (15) is an error estimate in terms of the computational cost required to solve the approximate problem.

Periodic Boundary Conditions
For simulating point defects (and often even dislocations), periodic boundary conditions appear to be by far the most popular choice. To implement periodic boundary conditions, let ω R ⊂ R d be connected such that B R ⊂ ω R , and Let Ω R := ω R ∩ Λ be the periodic computational domain and Ω per R := α∈Z d (Bα + Ω R ) the periodically repeated domain (with an infinite lattice of defects). For simplicity, suppose that ω R is compatible with T Λ , that is, there exists a subset T R ⊂ T Λ such that clos(ω R ) = ∪T R . The space of admissible periodic displacements is given bẏ The energy functional for periodic relative displacements u ∈Ẇ per (Ω R ) is given by For this definition to be meaningful, we assume for the remainder of the discussion of periodic boundary conditions that The computational task is to solve the finite-dimensional optimisation problem Proof (Idea of proof). The proof proceeds much in the same manner as for Theorem 2, but some details are more involved due to the fact that (16) is not a Galerkin projection of (8). The main additional difficulty is that the strong convergence ∇T Rū | ω R → ∇ū| ω R does not immediately imply stability of the periodic hessian, that is, To prove this, we consider the limit as R → ∞ (with an arbitrary sequence of associated domains Ω R ) and decompose test functions into a core and a far-field "sufficiently slowly". We then show that stability of δ 2 E (ū) implies positivity of H R v co , v co while stability of the homogeneous lattice (10) implies positivity of H R v ff , v ff . The cross-terms vanish in the limit. Thus, we obtain (18) for sufficiently large R. The details are given in the preprint [16]. 2. Compared with Theorem 2 we now only control the geometry in the computational domain Ω R . We can, however, "post-process" to obtain a global defect geometryv per := Π Rū per R (slightly abusing notation sinceū per R ∈Ẇ 1,2 (Λ)), for which we still get the estimate ∇v per − ∇ū L 2 C R −d/2 .

Boundary Conditions from Linear Elasticity
We consider a scheme where the elastic far-field of the crystal is approximated by linearised lattice elasticity. Our formulation is inspired by classical as well as recent multiscale methods of this type [22,42,43,45], but simplified to allow for a straightforward analysis.
We fix a computational domain Ω R ⊂ Λ such that B R ∩Λ ⊂ Ω R (for R R def ) and we linearise the interaction outside of Ω R This results in a modified approximate energy difference functional Analogously to Lemma 1 it follows that E lin R can be extended by continuity to a functional E lin R ∈ C k (Ẇ 1,2 ). Thus, we aim to compute Remark 4. The optimisation problem (20) is still infinite-dimensional, however, by defining Ω R := Ω R ∪ ∈Ω R and the effective energy functional for any u : Ω R → R m , it can be reduced to an effectively finite-dimensional problem. The reduced energy E red Ω can be computed efficiently employing lattice Green's functions or similar techniques [22,42,43,45]. This process introduces additional approximation errors, which we ignore subsequently. Thus, we only present an analysis of an idealised scheme, as a foundation for further work on more practical variants of (20).

Theorem 4.
Letū be a strongly stable solution to (8), then there exist C, R 0 > 0 such that for all domains Ω R ⊂ Λ with B R ∩ Λ ⊂ Ω R and R R 0 , there exists a strongly stable solution of (20) satisfying Proof (Idea of proof). The computational space is the same as for the full atomistic problem, hence the error is determined by the consistency error committed when we replaced V with V lin in the far-field. This error is readily estimated by a remainder in a Taylor expansion, which immediately implies that After establishing also stability of δ 2 E lin R , which follows from a similar argument we obtain ∇u lin R − ∇ū L 2 Dū 2 4 (Λ\Ω R ) , and employing the decay estimate (11) yields the stated error bound. The details of the proof are given in the preprint [16].

Boundary Conditions from Nonlinear Elasticity
A natural progression is to employ nonlinear elasticity in the far-field instead of linear elasticity, and determine whether this can improve further upon the approximation error. In this context it is only meaningful to employ continuum nonlinear elasticity, since our original atomistic model can already be viewed as a lattice nonlinear elasticity model. This leads us to considering a class of multiscale schemes, atomistic-to-continuum coupling methods (a/c methods), that has received considerable attention in the numerical analysis literature in recent years. We refer to the review article [25] for an introduction and comprehensive references. A key conceptual difference, from an analytical point of view, between a/c methods and the methods we considered until now is that they exploit higher-order regularity, that is, the decay of D 2ū , rather than only decay of Dū. Methods of this kind were pioneered, for example, in [30,40,41,46].
Due to the relative complexity of a/c coupling schemes we shall not present comprehensive results in this section, but instead illustrate how existing error estimates can be reformulated within our framework. This extends previous works such as [31,33,35] and presents a framework for ongoing and future development of a/c methods and their analysis; see for example [23,24,26,36], and references therein, for works in this direction.
We choose an atomistic region Ω a R ⊂ Λ, an interface region , and let I R denote the corresponding nodal interpolation operator. We let R and R c denote the sizes of Ω a R and ω R in the sense that for some c 0 > 0. As space of admissible displacements we definė We consider a/c coupling energy functionals, defined for u ∈Ẇ 0 (T R ), of the form where the various new terms are defined as follows: is the effective volume of T , where vor( ) denotes the Voronoi cell associated with the lattice site ; is an interface potential, which specifies the coupling scheme; is the Cauchy-Born strain energy function, which specifies the continuum model.
With this definition it is again easy to see that E ac R ∈ C k (Ẇ 0 (T R )). We now aim to compute The choice of the interface site-potentials V i is the key component in the formulation of a/c couplings. Many variants of a/c couplings exist that fit within the above framework [25]. In order to demonstrate how to apply our framework to this setting, we shall restrict ourselves to QNL type schemes [14,35,41], but our discussion applies essentially verbatim to other force-consistent energy-based schemes such as [28,38,39]. For other types of a/c couplings the general framework is still applicable; see in particular [24] for a complete analysis of blending-type a/c methods.
As a starting point of our present analysis we assume a result that is proven in various forms in the literature, for example, in [28,33,35]: we assume that there exist η > 0 and c 1 > 0 such that there exists a strongly stable solutionū ac R to (24) satisfying provided that h D 2ū 2 (Λ∩(ω R \B R )) + Dū 2 (Λ\B Rc/2 ) η. Such a result follows from consistency and stability of an a/c scheme and applying of the Inverse Function Theorem along similar lines as in the preceding sections.

Proposition 1.
Letū be a strongly stable solution of (8) and assume that (22) and (25) hold. Further we require that ω R and T R satisfy the following quasi-optimality conditions: Then there exists a constant C, depending on η, c 2 , c 3 , c 4 , and β such that, for R sufficiently large, Proof (Idea of proof). The proof consists in estimating the right-hand side of (25) in terms of R. Note that assuming (26) ensures that the truncation term Dū 2 (Λ\B Rc/2 ) does not dominate the coarse-graining term h D 2ū 2 (Λ∩(ω R \B R )) .
Remark 5. 1. It is interesting to note that an atomistic/continuum coupling is not competitive when compared against coupling to lattice linear elasticity. The primary reason for this is that the loss of interaction symmetry which causes a first-order coupling error at the a/c interface (the finite element error could be further reduced by considering higher order finite elements [13]).
2. Using our framework, the analysis in [13] suggests that one can generically expect the rate R −d−2 for the energy error.
3. To convert (27) into an estimate in terms of computational complexity, we note that, if we also have |h(x)| c 5 (|x|/R) β with β > 1, then the total number of degrees of freedom (in the atomistic and continuum region) is bounded by N dof C R d . The error estimate then reads

Setup We present an examples of an interstitial-type point defect in a two-dimensional triangular lattice
cf. Figure 1. (A second example, with a vacancy defect is presented in the preprint [16].) The reference configuration is given by Λ = AZ 2 ∪ {(1/2, 0)}. For each ∈ Λ, let R( ) denote the set of directions connecting to , defined by the bonds displayed in Figure 1. Then, the site energy is defined by To compute the equilibria we employ a robust preconditioned L-BFGS algorithm specifically designed for large-scale atomistic optimisation problems [37]. It is terminated at an ∞ -residual of 10 −7 .
We exclusively employ hexagonal computational domains. We slightly redefine N , letting it now denote the number of atoms in the inner computational domain, that is, #Ω N in the ATM-DIR, ATM-PER and LIN methods and #Ω a N in the AC method. Then, our analysis predicts the following rates of convergence for both model problems, where the rate N −2 for the energy in the AC case is predicted in [13]. We make some final remarks concerning the LIN and AC methods:

LIN:
For the experiments in this paper, we did not implement a variant based on Green's functions. Instead, we chose as an inner domain Ω N a hexagon of side-length K ≈ N 1/2 within a larger domain of a hexagon of side-length K 3 . This modification of the method does not affect the convergence rates. AC: To generate the finite element mesh we use the construction described in [23,26].

Discussion of Results
The graph of N versus the geometry error and the energy error are presented in Fig. 2. All slopes are as predicted with mild preasymptotic regimes for the ATM-PER and AC methods. The only exception is the energy for the LIN method, which displays a faster decay than predicted. We can offer no explanation at this point.

Dislocations
We now present an atomistic model for dislocations and analogous regularity and approximation results. To avoid excessive duplication we will occasionally build on and reference Section 2. Our presentation also builds on the descriptions in [2,19]. For more general introductions to dislocations, including modeling aspects as well as analytical and computational solution strategies we refer to [7,18].

Atomistic Model
We consider a model for straight dislocation lines obtained by projecting a threedimensional crystal. Briefly, let BZ 3 denote a three-dimensional Bravais lattice, oriented in such a way that the dislocation direction can be chosen parallel to e 3 and the Burgers vector can be chosen as We consider displacements W : BZ 3 → R 3 of the three-dimensional lattice that are independent of the direction of the dislocation direction, that is, e 3 . Thus, we choose a projected reference lattice Λ := AZ 2 := {( 1 , 2 ) | ∈ BZ 3 }, and identify W (X ) = w(X 12 ), where w : Λ → R 3 , and here and throughout we write a 12 = (a 1 , a 2 ) for a vector a ∈ R 3 . It can be readily checked that this projection is again a Bravais lattice.
We may again choose a regular triangulation T Λ satisfying T Λ + ρ = T Λ for all ρ ∈ Λ. Each lattice function v : Λ → R m has an associated P1 interpolant I v : R 2 → R m and we identify ∇v = ∇ I v. Further, we recall the definition of the spacesẆ c ,Ẇ 1,2 from (4).
x 1 } the "branch-cut" (cf. (31)), chosen such that Γ ∩ Λ = ∅. In order to model dislocations the site energy potential must be invariant under lattice slip. Normally, this is a consequence of permutation invariance of the site energy, but here we will formulate a minimal assumption. To that end, we define the slip operator S 0 acting on a displacement w : Λ → R 3 , or w : This operation leaves the three-dimensional atom configuration corresponding to that is, Y S represents only a relabelling of the atoms. Therefore, formally, if V (Dw) is the site energy potential as a function of displacement, then it must by invariant under the map w → S 0 w: In (34) below we will restate this assumption for a restricted class of displacements only, which will allow us to continue to employ the finite range interaction assumption.
Dislocations in an infinite lattice store infinite energy due to their topological singularity. We therefore decompose the total displacement w = u 0 +ū into a farfield predictor u 0 and a finite energy core correctorū ∈Ẇ 1,2 . There is no unique way to specify u 0 , but a natural choice is the continuum elasticity solution: For a function u : R 2 \ Γ → R m that has traces from above and below, we denote these traces, respectively, by where the tensor C is the linearised Cauchy-Born tensor (derived from the interaction potential V ; see Section 7 for more detail). In our analysis we require that applying the slip operator to the predictor map u 0 yields a smooth function in the half-space wherer is defined in Lemma 2 below. That is, we require that S 0 u 0 ∈ C ∞ (Ω Γ ). Except in the pure screw dislocation case (b 12 = 0) u lin does not satisfy this property. The origin of this conundrum is that linearised elasticity assumes infinitesimal displacements, yet we apply it in the large deformation regime near the defect core.
To overcome this technical difficulty, instead of u 0 = u lin , we define the predictor arg(x) denotes the angle in (0, 2π) between b 12 ∝ e 1 and x, and η ∈ C ∞ (R) with η = 0 in (−∞, 0], η = 1 in [1, ∞) and η > 0 in (0, 1). While the distinction between u 0 and u lin is crucial, it arises from a subtle technical issue and could be ignored on a first reading, especially in view of the following lemma. (10), then u lin is well-defined. Forr sufficiently large, ξ :

Lemma 2. (i) Suppose that the lattice is stable
Proof. The proof is given in Section 5.3. Statement (ii) implies that the net-Burgers vector of u 0 (and hence of any u 0 + u, u ∈Ẇ 1,2 ) is indeed b. Moreover, the fact that S 0 u 0 ∈ C ∞ (Ω Γ ) will allow us to perform Taylor expansions of finite differences.Statement (iii) indicates that u 0 is an approximate far-field equilibrium, which allows us to use u 0 as a far-field boundary condition (see Lemma 3 below).
In order to keep the analysis as simple as possible we would like to keep the convenient assumption made in the point defect case of a finite interaction range in reference configuration. At first glance this contradicts the invariance of the site energy under lattice slip (29), but we can circumvent this by restricting the admissible corrector displacements. Arguing as in [16, Section B.1] we may choose sufficiently large radiir A ,m A and define Upon choosingm A ,r A sufficiently large, we can ensure that any potential equilibrium solution is contained in A . Thus, the restriction of admissible displacements to A is purely an analytical tool, which ensures that we can treat V as having finite range, despite admitting slip-invariance. For We can now rigorously formulate the assumptions on the site energy potential: such that for each u ∈ A , and w = u 0 + u, the site energy associated with a lattice site is given We assume again that V (0) = 0 (that is, V is the energy difference from the reference lattice) and that R, V are point symmetric (6). We shall assume throughout that V is invariant under lattice slip, reformulating (30) as In addition, to guarantee lattice stability both before and after shift we assume that not only D but also S −1 DS include nearest-neighbour finite differences (or equivalent): The global energy (difference) functional is now defined by Proof (Idea of the proof). The main idea is the same as in the point defect case. The proof that δE (0) ∈Ẇ −1,2 is based on the construction of u 0 in terms of the linear elasticity predictor u lin , which guarantees that u 0 is an "approximate equilibrium" in the far-field. See [20] for a similar proof applied in the simplified context of a screw dislocation. The complete proof (given in Section 5.4) for our general case requires a combination of the proof in [20] and the concept of elastic strain introduced in Section 3.2.
The variational problem for the dislocation case is Since A is open, if a minimiserū exists, then δE (ū) = 0. We call a minimiser strongly stable if, in addition, it satisfies the positivity assumption (9).

Remark 6.
One can also formulate anti-plane models for pure screw dislocations by restricting A to displacements of the form u = (0, 0, u 3 ) and also computing a predictor of the form u lin = (0, 0, (u lin ) 3 ). Note also that for pure screw dislocations, (33) is ignored. In the anti-plane case we may also choose A =Ẇ 1,2 since only slip-invariance in anti-plane direction is required, that is, the topology of the projected two-dimensional lattice remains unchanged.
To define in-plane models for pure edge dislocations one restricts A to displacements of the form u = (u 1 , u 2 , 0). The predictor u lin does not simplify in this case.
All our results carry over trivially to these simplified models.

Remark 7. The definition of the reference solution with branch-cut
In this case the predictor solution u 0 would be replaced with S 0 u 0 . Let the resulting energy functional be denoted by It is straightforward to see that, if δE (ū) = 0, then δE S (Sū) = 0 as well. This observation means, that in certain arguments, an estimate onū in the left half-space where no branch-cut is present immediately yields the corresponding estimate on Sū in the right half-space as well.

Remark 8.
Another source of arbitrariness comes from the precise definition of the predictor u 0 , for example, through the choice of the dislocation core position x or the choice of smearing function η. Indeed, more generally, arbitrary smooth modifications to u 0 are allowed as long as they do not significantly change the farfield behaviour. While such changes to the predictor u 0 affect the resulting corrector u (the solution of (37)), the total displacement u 0 +ū remains unchanged in the sense that, if u 0 = u 0 + w 0 is a modified predictor, thenū =ū − w 0 is again a solution of (37).

Elastic Strain
The transformation u 0 → S 0 u 0 produces a map that is smooth in Ω Γ , and which generates the same atomistic configuration. It is therefore natural to define the elastic strains The analogous definition for the corrector displacement u is otherwise.
The slip invariance condition (34) can now be rewritten as Linearity of S and hence ofD implies δV (D(u 0 + u)), Dv = δV (e +Du),Dv , and so forth.

Regularity
The regularity of the predictor u 0 is already stated in Lemma 2. We now state the regularity of the correctorū. It is interesting to note that the regularity of the dislocation correctorū is, up to log factors, identical to the regularity of the displacement field in the point defect case, which indicates that the dislocation problem is computationally no more demanding than the point defect problem. Indeed, this will be confirmed in Section 3.4.
In the pure screw case where b 12 = 0 we simply have D =D.

Clamped Boundary Conditions
To extend clamped boundary conditions to the dislocation problem, we prescribe the displacement to be the predictor displacement outside some finite computational domain Ω R ⊂ Λ. Thus, we may think of these boundary conditions as asynchronous continuum linearised elasticity boundary conditions. This amounts to choosing a corrector displacement space analogous toẆ 0 (Ω R ) in the point defect case, and the associated finite-dimensional optimisation problem reads Theorem 6. Letū be a strongly stable solution to (37), then there exist C, Proof. See the preprint [16].

Periodic Boundary Conditions
It is possible to extend periodic boundary conditions to the dislocation case by considering a periodic array of dislocations with alternating signs. In practise the computational domain then contains a dipole or a quadrupole. It then becomes necessary to estimate image effects, for which our regularity results are still useful, but which requires substantial additional work. Hence, we postpone the analysis of periodic boundary conditions for dislocation to future work, but refer to [8] for an interesting discussion of these issues.

Boundary Conditions from Linear Elasticity
We now extend the lattice linear elasticity boundary conditions to the dislocation case. The linearisation argument (19) should now be carried out for the full displacement w = u 0 + u, and reads but this is invalid whenever the interaction neighbourhood crosses the slip plane Γ . Instead, we must first transform the finite difference stencils as follows: recall the definition of Ω Γ from (32) and the definition of elastic strain e andDu from (38) and (39), then we definẽ According to Lemma 2 and Theorem 5, if u =ū, then |D 0 w( )| = O(| | −1 ), hence we may linearize with respect to this transformed finite different stencil. Using the slip invariance condition (34), we obtain and we therefore define the energy difference functional where V lin is the same as in the point defect case, where Ω R ⊂ Λ is the "inner" computational domain. It follows from minor modifications of the proof of Lemma 3 that E lin R can be extended by continuity to a functional E lin R ∈ C k (A ). Thus, we aim to compute Theorem 7. Letū be a strongly stable solution to (37), then there exist C, R 0 > 0 such that for all domains Ω R ⊂ Λ with B R ∩ Λ ⊂ Ω R and R R 0 , there exists a strongly stable solution of (47) satisfying Proof (Idea of proof). The proof is similar to the point defect case, the main additional step to take into account being that the linearisation is with respect to the full displacement u 0 +ū. Since ∇u 0 ∼ |x| −1 it therefore follows that the linearisation error at site is only of order O(| | −2 ), while in the point defect case it was of order O(| | −2d ). This accounts for the reduced convergence rate. See the preprint [16] for the complete proof.

Remark 10.
The key difference between the schemes (44) and (47) is that the former employs a precomputed continuum linear elasticity boundary condition while the latter computes a lattice linear elasticity boundary condition on the fly. It is therefore interesting to note that, for dislocations, solving the relatively complex exterior problem yields almost no qualitative improvement over the basic Dirichlet scheme (44). Indeed, if the cost of solving the exterior problem is taken into account as well, then the scheme (47) may in practice become more expensive than (44).
The main advantage of (47) appears to be that the boundary condition need not be computed beforehand, but could be computed "on the fly". We speculate that this can give a substantially improved prefactor when the dislocation core is spread out, for example, in the case of partials.

Boundary Conditions from Nonlinear Elasticity for Screw Dislocations
The formulation of a/c coupling methods for general dislocations is not straightforward. We therefore consider only the case of pure screw dislocations and postpone the general case to future work. Thus, we assume that b = e 3 , and in this case, only the invariance of V in the normal direction is relevant: We set up the computational domain and approximation space as in Section 2.6. To define the energy functional, we first construct a modified interpolant that takes into account the discontinuity of the full displacement across the slip plane, similarly to the elastic strain used in Section 3.6, where I R is the nodal interpolation with respect to T R . With this definition, the energy difference functional is given by where V i , W, v eff T are defined as in Section 2.6. We seek to compute We again let R and R c be the sizes of Ω a R and ω R , and assume that there exists η > 0 and c 1 > 0 such that there exists a strongly stable solutionū ac R to (50) satisfying provided that hD 2 (u 0 +ū) 2 (Λ∩(ω R \B R )) + Dū 2 (Λ\B Rc/2 ) η. Letū be a strongly stable solution of (37) and assume that (51) and (52) hold. Further we require that ω R and T R satisfy the following quasi-optimality conditions: Then there exist R 0 , C depending on η, c 2 , c 3 , c 4 , and p, such that for all R R 0 there exists a strongly stable solutionū ac R to (50) satisfying ∇ū ac Proof. See the preprint [16].

Setup
We consider the anti-plane deformation model of a screw dislocation in a BCC crystal from [20], the main difference being that we admit nearest neighbour many-body interactions instead of only pair interactions. Thus, we only give a brief outline of the model setup. The choice of dislocation type is motivated by the fact that the linearised elasticity solution is readily available. Briefly, let BZ 3 = Z 3 ∪ (Z 3 + (1/2, 1/2, 1/2) T ) denote a BCC crystal, then both the dislocation core and Burgers vector point in the (1, 1, 1) T direction. Upon rotating and possibly dilating, the projection AZ 2 of the BCC crystal is a triangular lattice, hence we again assume (28). The linear elasticity predictor is now given by u lin (x) = 1 2π arg(x −x), where we assumed that the Burgers vector is b = (0, 0, 1) T andx is the centre of the dislocation core.
Let the unknown for the anti-plane model, the displacement in e 3 direction, be denoted by z( ) := y 3 ( ), then we use the EAM-type site potential with φ(r ) = ψ(r ) = sin 2 (πr ) and G(s) = 1 2 s 2 . The 1-periodicity of φ, ψ emulates the fact that displacing a line of atoms by a full Burgers vector leaves the energy invariant.
We apply again the remaining remarks in Section 2.7.1.

Discussion of Results
We choosex = (1/3, 1/(2 √ 3)) T , which places the dislocation core slightly off-centre to avoid spurious effects of a high-symmetry anti-plane setting. Two additional tests (including the high-symmetry case) are presented in the preprint [16].
The results are shown in Fig. 4. We observe precisely the predicted rates of convergence. However, it is worth noting that although the asymptotic rates for ATM, LIN and AC are identical (up to log-factors), the prefactor varies by an order of magnitude. The "dip" in the energy error for the LIN method is likely due to a change in sign of the error.
The prefactor is a crucial piece of information about the accuracy of computational schemes that our analysis does not readily reveal. Ideally, one would like to establish estimates of the form Dū where C * and p can be given explicitly, however much finer context-sensitive estimates would be required to achieve this.

Conclusion
We have introduced a flexible analytical framework to study the effect of embedding a defect in an infinite crystalline environment. Our main analytical results are (1) the formulation of equilibration as a variational problem in a discrete energy space; and (2) a qualitatively sharp regularity theory for minimisers.
These results are generally useful for the analysis of crystalline defects, however, our own primary motivation was to provide a foundation for the analysis of atomistic multi-scale simulation methods, which in this context can be thought of as different means to produce boundary conditions for an atomistic defect core simulation. To demonstrate the applicability of our framework we analyzed simple variants of some of the most commonly employed schemes: Dirichlet boundary conditions, periodic boundary conditions, far-field approximation via linearised lattice elasticity and via nonlinear continuum elasticity (Cauchy-Born, atomistic-tocontinuum coupling). In parallel works [12,23,24,36] this framework has already been exploited resulting in new and improved formulations of atomistic/continuum and quantum/atomistic coupling schemes.

Proofs: The Energy Difference Functionals
This section is concerned with proofs for Lemmas 1 and 3 which state that the energy E can be understood as a smooth functional on the energy space, that is, E ∈ C k (Ẇ 1,2 ) in the point defect case and E ∈ C k (A ) in the dislocation case.

Conversion to Divergence Form
We begin by establishing an auxiliary result that allows us to convert pointwise forces into divergence form without sacrificing fundamental decay properties. If f has compact support, then g can be chosen to have compact support as well.
Since p > d, we obtain that [ f (n) ] p → 0. Moreover, since ∈Z d f (n) ( ) = 0 for all n it follows that f (n) 1 → 0. Further, (58) implies and hence the series ∞ n=0 g (n+1) − g (n) converges. Let g( ) := lim n→∞ g (n) ( ), then (57) implies that g satisfies the identity in (55), and the bound on g p−1 implies the inequality in (55). It remains to note that if f = f ( ) = 0 outside the region | | ∞ L for some L, then f (n) , g (n) , and hence g, are also zero outside this region.
To show the first inequality in (58), we fix = 0, express f + ( ) through f ( ), and estimate The second inequality in (58) is based on the following two estimates: where we denote again (f ,g) := C d ( f, g) and Δg := (g −g)·e d . The first estimate follows from arguments similar to the above. The second estimate, for = (¯ , d ) with d 0, is proved in the following calculation: where we used that for d 0, |(¯ , d −1)| ∞ 1 and the fact that x − 1 2 1 3 (x + 1 2 ) for any x 1. For d > 0 this estimate is obtained in a similar way.
The analogous estimates hold for applications of C d−1 , . . . , C 1 and combining these yields the second inequality in (58).
In addition, if d = 2, under the assumptions of Section 3.1, there existsg : If f has compact support, then g andg can be chosen to have compact support as well.
Proof. One only needs to notice that the assumptions that the operators D andD contain nearest-neighbor finite differences (cf. (7) and (35)) allow to use Lemma 4 to construct the needed g andg.

Proof of Lemma 1
The proof relies on two prerequisites.
Proof. For a very similar argument that can be followed almost verbatim see [34], hence we only give a brief idea of the proof.
Since 1], and . It follows that where C u depends only on |Du| ∞ . In particular, Using similar lines of argument one can prove that F ∈ C k (A ).

Lemma 6. Under the conditions of Lemma
Proof. Let v ∈Ẇ c , then we can write the first variation in the form where f ( ) is given in terms of the V ,ρ ; the precise form is unimportant. Point symmetry of the lattice implies that f ( where the inequality u − u(0) 2 (Λ∩B R def +r cut ) ∇u L 2 (B R def +r cut ) follows from the fact that only finite-dimensional subspaces are involved, and for these it is enough to see that for any u such that the right-hand side vanishes, the left-hand side must vanish as well. But this is immediate. This completes the proof.
For u ∈Ẇ c , which according to the two foregoing Lemmas is continuous with respect to thė W 1,2 -topology and thus has a unique extension toẆ 1,2 . Since the first term is C k and the second is linear and bounded, the result E ∈ C k follows as well. This completes the proof of Lemma 1.

Proof of Lemma 2 (Properties of the Dislocation Predictor)
We begin by analyzing the auxiliary deformation map ξ defined in (33) in more detail. To simplify the notation let ζ := ξ −1 throughout this section.

Lemma 7. (a) Ifr is sufficiently large, then
can be continuously extended to the half-space Ω Γ = {x 1 >r + b 1 }, and after this extension we have ζ S ∈ C ∞ (Ω Γ ).
is clearly injective, it follows x 1 = x 1 as well.
(b) The map ξ leaves the x 2 coordinate unchanged and only shifts the x 1 coordinate by a number between 0 and b 1 . Thus, forr /4 > |b 1 |, the statement clearly follows.
In general the proof proceeds by induction. Suppose the result is true for all α with |α| 1 m.
Finally, the coefficient functions (∂ s ζ i − δ is ) are readily seen to also satisfy the same regularity and decay as stated for g α,β with any |β| 1 = |α| 1 . This concludes the proof.
We can now proceed to the proof of Lemma 2. Proof of (i): u 0 is well-defined. The elasticities tensor C is derived from the interaction potential and due to the lattice stability assumption (10) satisfies the strong Legendre-Hadamard condition (see Section 7 for more detail). The existence of a solution to (31) then follows from [18, Section 13-3, Equations 13-78], with ∇u lin ∈ C ∞ (R 2 \ {x}) and and |∇ j u lin (x + x)| |x| − j . Lemma 7 immediately implies that u 0 is also well-defined. This completes the proof of Lemma 2 (i).
Proof of (ii) Let x ∈ Γ ∩ Ω Γ , then For derivatives of arbitrary order, the result is an immediate consequence of (59) and of Lemma 7(c).
Proof of (iii): This statement is again an immediate consequence of (59).

Proof of Lemma 3
The main idea of the proof is the same as in the point defect case, Section 5.2. For u ∈Ẇ c we write It is an analogous argument as in the point defect case to show that F ∈ C k (A ).
To prove that δE (0) is a bounded linear functional, we first use (41) to rewrite it in the form Next, we convert it to a force-displacement formulation, by generalising summation by parts to incompatible gradientsD.
Proof. We let k ∈ Λ and u ∈Ẇ 1,2 , u( ) := δ k . Then we form the expression and show that it vanishes. This result is geometrically evident, but could also be proved by a direct (yet tedious) calculation whose details we omit.
We can now deduce that To prove that δE (0) is bounded we must establish decay of f . For future reference, we establish a more general result than needed for this proof.
Lemma 10. Let f be given by (61), and 0 j k − 2, then there exists C such that Proof. Throughout this proof we will implicitly assume that | | is sufficiently large so that the defect core Br +|b| (x) does not affect the computation. We first consider the case j = 0. Case 1: left halfspace: We first consider the simplified situation when 1 < x 1 , that is we can simply replaceD ≡ D throughout. We will see below that a generalisation to 1 >x 1 is straightforward.
To estimate the first group we expand We have therefore shown (62) for the case j = 0, when lies in the left half-space. Case 2: right halfspace: To treat the case 1 >x 1 , | | sufficiently large, we first rewrite Since S 0 u 0 is smooth in a neighbourhood of | | (even if that neighbourhood crosses the branch-cut), we can now repeat the foregoing argument to deduce again that |S f ( )| | | −3 as well (cf. Remark 7). But since S represents an O(1) shift, this immediately implies also that | f ( )| | | −3 . This completes the proof of (62).
Proof for the case j > 0: To prove higher-order decay, assume again at first that 1 <x 1 and consider τ ∈ R j , j 1, then (2) .

An analogous Taylor expansion as above yields
The term f (2) is readily estimated by multiple applications of the discrete product rule, from which we obtain that | f (2) The generalisation to the case 1 >x 1 is again analogous to above, due to the fact thatD From this point, the argument continues verbatim to the case 1 <x 1 .

Proofs: Regularity
In this section we prove the regularity results, Theorems 1 and 5.

First-Order Residual for Point Defects
Assume, first, that we are in the setting of the point defect case, Section 2.1. To motivate the subsequent analysis we first convert the first-order criticality condition δE (ū) = 0 for (8).
Finally, to state the first auxiliary result, we recall from Section 2.1 the definition of the interpolant I u for discrete displacements u : Λ → R d , which provides point values I u( ) for all ∈ AZ d . Lemma 11. (First-order residual for point defects) Under the assumptions of Theorem 1 there exists g : AZ d → (R m ) R and R 1 , C > 0 such that Proof. Let u ≡ Iū. We rewrite the residual Hu, v as The first group can be written as and where we note that g 1 ( ) is a linearisation remainder and hence |g 1 ( )| |Du( )| 2 for | | sufficiently large. The second group is the residual of the exact solution after projection to the homogeneous lattice AZ d . Writing this group in "force-displacement" format, we observe that f ( ) = ρ∈R D −ρ V ,ρ (Du( )) has zero mean as well as compact support due to symmetry of the lattice. Because of the mean zero condition, we can write it in the form f, v = g 2 , Dv where g 2 also has compact support (cf. Corollary 1).
Finally, the third group vanishes identically, which can for example be seen by summation by parts. Setting g = g 1 + g 2 this completes the proof.

The Lattice Green's Function
To obtain estimates onū and its derivatives from (66) we now analyse the lattice Green's function (inverse of H ). The following results are widely expected but we could not find rigorous statements in the literature in the generality that we require here.
Using translation and inversion symmetry of the lattice, the homogeneous finite difference operator H defined in (65) can be rewritten in the form where R := {ρ + ς | ρ, ς ∈ R} \ {0} and A ρ ∈ R d×d . (Written in terms of V ,ρς , A ρ = ς,τ ∈R,ς −τ =ρ V ,ς τ . Alternatively, one can define A ρ = −2 ∂ 2 Hu,u ∂u(0)∂u(ρ) and arrive at the same result; cf. [19,Lemma 3.4].) Since the Green's function estimates hold for general operators of the form (69) we recall the associated stability for some γ > 0. Next, we recall the definitions of the semi-discrete Fourier transform and its inverse, and where B ⊂ R d is the first Brillouin zone. As usual, the above formulas are wellformed for u ∈ 1 (AZ d ; R m ) andû ∈ L 1 (B; R m ), and are otherwise extended by continuity. Transforming (69) to Fourier space, we get Lattice stability (70) can equivalently be written asĤ (k) γ |k| 2 Id. Thus, if (70) holds, then the lattice Green's function can be defined by We now state a sharp decay estimate for G.

Lemma 12. Let H be a homogeneous finite difference operator of the form (69)
satisfying the lattice stability condition (70), and let G be the associated lattice Green's function. Then, for any ρ ∈ R j , j > 0, or j = 0 if d = 3, there exists a constant C such that Proof. The strategy of the proof is to compare the lattice Green's function with a continuum Green's function.
Letη(k) ∈ C ∞ c (B), withη(k) = 1 in a neighbourhood of the origin. Then, it is easy to see that its inverse (whole-space) Fourier transform η := F −1 [η] ∈ C ∞ (R d ) with superalgebraic decay. From this and (73) it is easy to deduce that where C = C( j, H ) and α ∈ R j is the multi-index defined in the statement of the theorem.
Step 2: Comparison of Green's Functions: Our aim now is to prove that which implies the stated result. (In fact, it is a stronger statement.) We write The explicit representations ofĜ andĜ make it straightforward to show that (one employs the fact thatĜ −1 −Ĝ −1 has a power series starting with quartic terms) Thus, if d − 1 + j is even and we choose 2n := d − 1 + j, then we obtain that Δ n (Ĝ −Ĝ) p α (k) ∈ L 1 (B), which implies that which is the desired result (75).
If d−1+ j is odd, then we can deduce (75) from the result for a larger multi-index α = (α, ρ ) of length j . Namely, fix ∈ AZ d and choose ρ a nearest-neighbour direction pointing away from the origin, then from which (75) easily follows.

Decay Estimates for Du, Point Defect Case
At the end of this section we prove Theorem 1 for the cases j = 0, 1. In preparation we first prove a more general technical result.

Lemma 13.
Let H be a homogeneous finite difference operator of the form (69) satisfying the stability condition (70). Let u ∈Ẇ 1,2 (AZ d ) satisfy p d and h ∈ 2 (AZ d ). Then, for any ρ ∈ R, there exists C 0 such that, for | | 2, Proof. Recall the definition of the Green's function G from Section 6.2 and its decay estimates stated in Lemma 12. Then, for all ∈ AZ d , it holds that and hence, for all σ ∈ R, Applying Lemma 12 and the assumption (76), we obtain For r > 0, let us define w(r ) := sup ∈AZ d , | | r |Du( )|. Our goal is to prove that there exists a constant C > 0 such that where The proof of (78) is divided into two steps.
Step 1: We shall prove that there exists a constant C > 0 and η : R + → R + , η(r ) −→ 0 as r → +∞, such that for all r > 0 large enough, Step 1a: Let us first establish that, for all | | 2r , we have We split the summation into |k| r and |k| > r . We shall write |k| r instead of k∈AZ d ,|k| r , and so forth.
Step 1b: Let us now consider the remaining group in (77), which we must again estimate for all | | 2r .
Step 2: Let us define v(r ) := r d z(r ) w(r ) for all r > 0. We shall prove that v is bounded on R + , which implies the desired result. Multiplying (79) There exists r 0 > 0 such that, for all r > r 0 , Cη(r ) 1 2 . This implies that, for all r > r 0 , Denoting F := sup r r 0 v(r ) and reasoning by induction, we obtain that, for all r > r 0 , where N (r ) C log(2+r ). Finally, the above inequality implies that v(r ) C + F and thus v is bounded on R + . This implies (78) and thus completes the proof of the lemma.

Decay Estimates for Higher Derivatives, Point Defect Case
From Section 6.3 we now know that |Du( )| C| | −d for | | sufficiently large, and more generally we can hope to, inductively, obtain that |D i u( )| C| | 1−d−i . Using this induction hypothesis we next establish additional estimates on the righthand side g in (66).
as well, which gives a first crude estimate for the decay. Exploiting this observation, the proofs of the higher-order decay estimates take a somewhat simpler form, as they need to address the nonlinearity.
where g is defined in (66).
Let α = (α 1 , . . . , α j ) ∈ R j be a multi-index. For any "proper subset" α = (α i ) i∈I , I {1, . . . , j}, we have by the assumptions made in the statement of the lemma that Thus, applying the discrete product rule (85), we obtain, for τ ∈ R s , s 2, Using, moreover, the estimates for τ ∈ R j+1 , we can conclude that This, together with (84), completes the proof.
To complete the proof of Theorem 1 we need a final auxiliary lemma that estimates decay for a linear problem.

Lemma 15. Let H be a homogeneous finite difference operator of the form (69)
satisfying the stability condition (70). Let u ∈Ẇ 1,2 (AZ d ) satisfy where p > d and j 0. Then, for i = 1, . . . , j and ρ ∈ R i , there exists C > 0 such that Proof. The proof is a straightforward application of the decay estimates for the Green's function. For the sake of brevity, we shall only carry out the details for the case j = 2. This will reveal immediately how to proceed for j > 2.
For all ∈ AZ d , ς, ς ∈ R, we have We again split the summation over |k| | |/2 =: r and |k| > r . In the set |k| > r the estimate is a simplified version (due to the absence of the nonlinearity) of Step 1b in the proof of Lemma 13, which yields |k|>r ρ∈R In the set |k| < r , we carry out a summation by parts, where χ r,ς (k), ν r,ς (k) ∈ {−1, 0, 1}. To see this, consider two discrete functions a, b and the characteristic function χ(k) = 1 if |k| r and χ(k) = 0 otherwise. Then, This establishes the claim that the coefficients χ r,ς , ν r,ς belong indeed to {−1, 0, 1}. The summation over |k| r + |ς | can be bounded using a simplified variant of the estimates in Step 1a of the proof of Lemma 13 and the decay assumption for g. This yields The "boundary terms" in (89) (second group on the right-hand side) are estimated by Thus, in summary, we can conclude that The only modification for the case j > 2 is that j − 1 summation by part steps are required instead of a single one. This completes the proof of Lemma 15.
Proof (Theorem 1, Case j 2). The statement of Theorem 1, Case j 2, is an immediate corollary of Lemmas 14 and 15.

Proof of Theorem 5, Case j = 1
We now adapt the arguments of the foregoing sections to the dislocation case. Remembering that Du 0 ( ) → 0 as | | → ∞ we begin by recalling the definitions of e =D 0 u 0 andDu from Section 3.2, noting that |e( )| | | −1 .
Let u :=ū, v ∈Ẇ c and | | sufficiently large, then (41) Upon defining the linear operator we obtain that We can now generalise Lemma 11 as follows.
is sufficiently large. Applying the bounds for η 1 and G again, we therefore obtain that Thus, we can estimate To summarize the proof up to this point, we have shown that, if | | is sufficiently large and if B 3| |/4 ( ) ∩ Γ = ∅, then Reflection argument: We now extend the argument to the case when B 3| |/4 ( )∩ Γ = ∅. According to Remark 7, we have δE S (Su) = 0 (recall that in the definition of E S we have replaced u 0 with S 0 u 0 ). This new problem is structurally identical to δE (u) = 0, except that the branch-cut Γ is now replaced with Γ S . Therefore, it follows that (95) holds, but u replaced with Su and for all 1 >x 1 , | | sufficiently large. It is now immediate to see that we can replace DSu with S −1 DSu =Du without changing the estimate. Thus we obtain In the preprint [16] we also give a direct argument for (96) using purely algebraic manipulations. Conclusion: We now consider arbitrary . We rewrite (96) in a way that allows us to apply the argument similar of Step 2 in the proof of Lemma 13. We begin by noting that Fix ε > 0, then there exists r 0 > 0 such that D u 2 (B s 2 ( )) ε, whenever | | r 0 .
We can now apply the argument of Step 2 in the proof of Lemma 13 to obtain that w(r ) r −1 and hence |Du( )| | | −1 .
Having established a preliminary pointwise decay estimate onDū, we now apply a boot-strapping technique to obtain an optimal bound. Proof (Theorem 5, Case j = 1). In view of Remark 7 (cf. part (2) in the proof of Lemma 17) we may assume, without loss of generality, that belongs to the left half-plane, that is, 1 <x 1 . We again define v and η 1  To estimate the first group we note that T 1 = g,Dv , hence we can employ the residual estimates from Lemma 16. Combining Lemma 16 with Lemma 17 we have |g(k)| |k| −2 , which readily yields Here, we used the observation that due to the fact that D ρ Sw(k) = D ρ w(k ), where |ρ − ρ | + |k − k | 1.

Proof of Theorem 5, Case j > 1
In view of case j = 1 and also of Lemma 15(b) it is natural to conjecture that Suppose that we have proven this for i = 1, . . . , j − 1. Then the triangle inequality immediately yields |D j u( )| | | − j log | |, which is of course sub-optimal, but it allows us again to apply a bootstrapping argument. In the dislocation case, this requires two steps, corresponding to cases (a) and (b) of the following lemma.
The term T 2 can be estimated analogously as in the proof of the case j = 1 in Section 6.5, noting that by the same argument as used there, |D D ρ G(k − )| |k − | − j−1 . Thus, one obtains The term T 1 : First, we split The second term is readily estimated, using |Dv(k)| | − k| − j−1 , by To estimate S 1 we first notice that, provided that | | is chosen sufficiently large, this sum only involves values of g, v away from Γ , that is,D ≡ D and we can write S 1 = |k− | | |/2 g(k), Dv(k) .
We are now in a position to mimic the argument of Lemma 15 almost verbatim, only having to take care to take into account the slower decay of g. Namely, according to the hypothesis stated at the beginning of the present proof, and employing Lemma 18 we have |D i g(k)| |k| −i−2 log r |k|. This in turn yields an additional log-factor in the estimate |S 1 | | | − j−1 log r +1 | |.
Conclusion: Arguing initially with r = 1, we obtain from the preceding arguments that |D j u( )| | | − j−1 log 2 | |. This initial estimate implies that, at the beginning of the proof, we may in fact choose r = 0, and therefore, we even obtain the improved bound |D j u( )| | | − j−1 log | |. Recalling that we assumed (without loss of generality) 1 <x 1 , so that in fact we have |D j u( )| | | − j−1 log | |, this completes the proof.

A.1 Cauchy-Born Model
Consider a Bravais lattice AZ d with site potential V : (R m ) R → R ∪ {+∞}. Consider the homogeneous continuous displacement field u : R d → R m , u(x) = Fx for some F ∈ R m×d . Then interpreting u as an atomistic configuration, the energy per unit undeformed volume in the deformed configuration u is If u, u 0 : R d → R m are both "smooth" (that is, |∇ 2 u(x)|, |∇ 2 u 0 (x)| 1), then is a good approximation to the atomistic energy-difference ∈AZ d V (Du( ))−V (Du 0 ( )). The potential W : R m×d → R ∪ {+∞} is called the Cauchy-Born strain energy function. Detailed analyses of the Cauchy-Born model are presented in [6,15,34]. In these references it is shown that both the Cauchy-Born energy and its first variation are second-order consistent with atomistic model, and resulting error estimates are derived.

A.2 Linearised Elasticity
A continuum linear elasticity model that is consistent with the atomistic description can be obtained by expanding the Cauchy-Born strain energy function W to second order: where we employed summation convention. Let C jβ iα := ∂ F iα F jβ W (0), then employing cancellation of the linear terms, we obtain the linearised energy-difference functional and the associated equilibrium equation is (This equation becomes non-trivial when supplied with boundary conditions or an external potential, either or both arising from the presence of a defect.) If the lattice AZ d is stable in the sense that, for some γ > 0, (cf. (9), (70)) then the tensor C satisfies the Legendre-Hadamard condition and hence the linear elasticity equations are well-posed in a suitable function space setting [19,44]. Finally, we remark that, the linear elasticity model can also be obtained by first deriving a quadratic expansion of the atomistic energy and then taking the long-wavelength limit (continuum limit). This yields the relationship between the continuum Green's function and the lattice Green's function exploited in the proof of Lemma 12.