Asymptotic Expansion of the Elastic Far-Field of a Crystalline Defect

Lattice defects in crystalline materials create long-range elastic fields which can be modelled on the atomistic scale using an infinite system of discrete nonlinear force balance equations. Starting with these equations, this work rigorously derives a novel far-field expansion of these fields: The expansion is computable and is expressed as a sum of continuum correctors and discrete multipole terms which decay with increasing algebraic rate as the order of the expansion increases. Truncating the expansion leaves a remainder describing the defect core structure, which is localised in the sense that it decays with an algebraic rate corresponding to the order at which the truncation occurred.


Introduction
The field of continuum solid mechanics has been highly successful in providing robust predictions of material behaviour at a wide range of length-scales. In crystalline materials in particular, it is recognised that the predictions made by the equations of linear elasticity are valid with tolerable errors even when resolving features such as defects whose characteristic size is close to that of the interatomic spacing. This said, when considering processes which involve the genuine thermodynamic, electronic and chemical properties of such defects, the fundamental discreteness of matter becomes crucial and no single continuum model can be sufficient to completely capture the fine detail of a material's behaviour at this scale. Moreover, it is exactly these fine details which determine phenomena such as a material's yield strength and behaviour under cyclic loading, both of which are crucial to understand for engineering applications.
As a result, a range of theoretical techniques have been developed over the last 60 years which seek to predict and compute defect behaviour in crystals, connecting discrete and continuum models of these materials. Broadly, these approaches can be divided into two categories, namely concurrent and sequential modelling strategies. Models in the former class combine discrete and continuum models into a single system which can then be numerically solved simultaneously, while those in the latter class generally involve iteration over separate models acting at different scales. In particular, the last 25 years has seen a great deal of research activity focused on concurrent strategies, with one of the most significant developments in this area being the quasicontinuum method [TOP96] and many variants thereof [LO13]. In contrast, the present work revisits sequential strategies for accurately modelling defects, but with a new perspective that recent progress in the study of multi-scale models has enabled, and our results on the structure of crystalline defects have direct consequences for the design and evaluation of more general models, including both purely atomistic and concurrently coupled models.
Starting from a discrete energy for a material defined on an infinite domain [EOS16], we develop a hierarchy of linear continuum PDE systems which can be efficiently derived and solved numerically. The solutions to these PDE systems form a sequence of smooth predictors which describe the far-field behaviour of the lattice strain around a defect to within arbitrary accuracy, and thus provide increasingly accurate boundary conditions to be used on the discrete model when confined to a finite domain. The key idea behind our approach is to exploit the knowledge that the variation in the strain field generally decays smoothly away from the core of any localised defect, and hence the far-field behaviour is more and more accurately predicted by the solutions of the continuum PDE models we define. These far-field approximate solutions are coupled with the properties of the discrete defect core, encapsulated by the spatial moments of the acting forces and expressed as a multipole expansion. The relative simplicity of these moments provides an elegant, computable way to transfer information from the nonlinear discrete problem to the continuum hierarchy. The coupling moreover is "weak" in the sense that a term in the continuum hierarchy only requires information on the multipole terms of strictly lower order. Indeed all terms are defined and are computable in sequence without concurrent coupling.
Our approach has connections with classical approaches to modelling defects in continuum linear elasticity using the defect dipole tensor [Esh56,NH63] (also known as the elastic dipole tensor or the double force tensor) and to Sinclair's work on atomistic models of fracture [SL72,Sin75]. More recent related work includes that of Trinkle and coworkers, where lattice Green's functions have been used to improve the accuracy of defect computations [Tri08,TT16], along with mathematical developments in our understanding of the regularity of discrete strain fields advanced by the authors and coworkers [HO12,EOS16,BS16,BBO19]. In particular, [BBO19] explores some of the initial ideas of our mathematical strategy in a simplified setting. The approach presented here unifies many of the ideas involved in these previous works into a single framework and expands them systematically to higher orders.
More generally, the powerful structural results we present here serve as a useful tool in any discussion of crystalline defects where high accuracy is required. In particular, we will outline how our results may be used as a foundation for a rigorous numerical analysis of defect algorithms and also provide a path to systematically improve their accuracy.

Methodology
Our starting point is to consider a total energy E(u) for displacements u of an infinite lattice Λ of atoms. We will make the mild assumption that the total energy of a displacement is expressed as a sum of site energies, i.e. contributions to the total energy arising from the environment of each atom. Under the assumption of frame indifference, this energy (and the site energies which make up the total energy) must depend only upon the relative displacements between atoms.
An equilibrium displacementū of the energy E satisfies the force balance equation and it is this infinite system of discrete equations which we study. Since equilibrium displacementsū exhibit decay properties away from the core of point defects and dislocations [EOS16], we can develop these equations around the state of zero displacement to derive approximate equations for the equilibrium displacementū in the far field. This comes in two steps: • As a first step, we can expand the site potential around zero to obtain linear and nonlinear lattice operators, which still depend on finite differences (atomic bonds); and • As a second step, we can use an expansion for the displacement itself to replace finite differences with a gradient and higher derivatives, obtaining continuum PDE approximations to the discrete lattice equations.
Note already, that the latter Taylor expansion requires sufficiently smooth continuum displacements and cannot be applied to the discreteū itself, so some care needs to be taken when applying this strategy. The PDE approximations come in the form of a hierarchy of corrector equations. Crucially, each one of these corrector equations has the same form: the continuum linear elasticity (CLE) equation must be solved but with different right-hand sides (forces) that depend on previous terms in the expansion. Such hierarchies of corrector equations have been explored by other authors before, e.g., see [EM07] for a formal hierarchy of similar corrector equations in the defect-free setting. In the present work, we take great care to define all terms rigorously, ensure sufficient regularity, and provide sharp estimates for all resulting error terms. We then combine these continuum correctors with a multipole series which can be obtained from the moments of linearised residual forces. Overall, this enables us to obtain far-field approximations to arbitrarily high order forū, which is characterised by a discrete remainder whose locality (decay) is precisely controlled.
While our general approach and many of our theoretical results are generic, some significant technical and conceptual challenges come into play when applying the results to specific defects. Specifically, questions regarding the geometry and the relation of a defect state to a reference lattice, as well as the precise properties of the lattice Green's function after adjusting for the geometry. These challenges can sometimes be overcome adhoc in leading order (as in [EOS16] for edge dislocations) but require a more detailed understanding for higher orders. Our main focus here is the general methodology. We will therefore restrict ourselves to two defect types that are geometrically relatively simple. Namely, general point defects and (straight) screw dislocations.

Outline
The paper is organised as follows: In § 2, we develop our general theory independent of specific defects. We present our main theoretical results concerning the decay of displacement fields and their approximation by multipole expansions for linearised lattice models in § 2.2 and then outline the derivation of the continuum corrector PDEs in § 2.3.
In § 3, we apply our methodology to our nonlinear atomistic models of point defects in Theorem 3.1 and screw dislocations in Theorem 3.5, and explain the implications for the convergence of numerical methods which exploit the result. § 4 presents our conclusions, and discusses the outlook for extending and applying these results in the future.
The proofs of our decay estimates for the linearised models are then provided in § 5, a full discussion of the lattice Green's function in § 6, and the proofs of Theorem 3.1 and 3.5 and in § 7.

Models and Notation
The Atomistic Energy and General Notation. Our results are concerned with modelling of crystalline defects, i.e. local regions of non-uniform atom arrangements embedded in a homogeneous host crystal. In this section we will start with the homogeneous setting itself. In § 3 we will then look at the description of defect configurations and discuss how the homogeneous results can be applied there.
The homogeneous crystal is described by a Bravais lattice Λ := AZ d as the reference, where d ≥ 2 and A ∈ R d×d is non-singular, and displacements u : Λ → R N , where we allow N = d in order to model a range of scenarios (for example, in some models for pure screw dislocations, we have d = 2, N = 3, while in anti-plane shear d = 2, N = 1).
We denote discrete differences by D ρ u(ℓ) := u(ℓ + ρ) − u(ℓ) for ℓ, ρ ∈ Λ. Later on we will also look at higher discrete differences which we denote by Next we look at an interaction neighbourhood R ⊂ Λ\{0}. We assume throughout that R is finite, R = −R, and that it spans the lattice span Z R = Λ. Based on R we define the discrete difference stencil Du(ℓ) := D R u(ℓ) := (D ρ u(ℓ)) ρ∈R . Again, later we will consider higher discrete differences and in particular apply D k times for which we use the simple notation D k u = D . . . Du. With the discrete difference stencil we can formally write down the energy of u as where V : R N ×R → R is the site energy. We will assume throughout that V ∈ C K (R N ×R ) for some K and satisfies the natural and very mild symmetry assumption Note that this is only a formal definition of the energy as the sum might not converge. All these quantities might also need to be adjusted to inhomogeneous generalisations E def , Λ def , R ℓ , and V ℓ to allow for a desired defect structure. Both of these aspects will be discussed in detail in § 3. In particular, we will make E precise and establish differentiability properties for any u in the discrete energy space H 1 which is defined by For future reference we also define the dense subspace H c ⊂ H 1 of displacements with compact support, Variations of E can then be written as More generally, we use the notation T [a 1 , . . . , a k ] for a multi-linear operator T . For symmetric operators we will shorten the notation further and write We denote the Hessian at zero by It will be convenient to interpret these objects as linear functionals, belonging to (H 1 ) * , acting on the last test function. We will often use a point wise representation based on the ℓ 2 scalar product. For example, where Div A = − ρ∈R D −ρ A ·ρ is the discrete divergence for a matrix field A : Λ → R N ×R . In this notation we can also write down the force equilibrium equations δE(u) = 0 in the pointwise form 0 = δE(u)(ℓ) = − Div ∇V (Du(ℓ)) .

Lattice Stability and The Green's Function
We assume throughout that the Hamiltonian H = δ 2 E(0) is stable and will equivalently call the lattice stable (see [HO12]), which by definition is the case if and only if there exists a c 0 > 0 such that For stable operators H there exists a lattice Green's function G : Λ → R N ×N such that H[Ge k ](ℓ) = e k δ ℓ,0 for all 1 ≤ k ≤ N and such that see [EOS16]. Here we used the notation |ℓ| 0 := |ℓ| + 2 to write decay rates in the discrete setting in a more compact form which we will do throughout. We will also often just write G k := Ge k .
The Cauchy-Born Continuum Model Our atomistic lattice model naturally gives rise to corresponding continuum model based on the Cauchy-Born rule. In the continuum setting, one has an energy of the form The energy density W is given by the Cauchy-Born rule for any M ∈ R N ×d , where c vol = |det A| > 0 is the volume of a lattice cell. We will later see in more detail the usefulness and limitations of the nonlinear continuum model for defect problems. We will also make heavy use of the linearised continuum problem for our corrector equations. These are given through the continuum Hamiltonian H C = δ 2 E C (0). In our pointwise notation the equilibrium equations is where C = ∇ 2 W (0). These are the standard continuum linear elasticity (CLE) equations. The lattice stability (4) in particular also implies the Legendre-Hadamard stability of C (see [HO12,BS16]), so that (7) is elliptic and allows for a continuum Green's function or fundamental solution G C : Notation for Tensors To work with various higher order tensor products throughout, we establish a precise but compact notation: Given a k-tuple of vectors in R d , σ = (σ (1) , . . . , σ (k) ) ∈ (R d ) k we denote their k-fold tensor product The vector space spanned by these tensor products is denoted (R d ) ⊗k , and it is easy to see this space is isomorphic to R d k . We also write when considering the k-fold tensor product of a single vector v ∈ R d .
Let S k denote the usual symmetric group of all permutations which act on the integers {1, . . . , k}. This action can be extended to k-tuples and tensor products by defining π(σ) := (σ (π(1)) , . . . , σ (π(k)) ) and π(σ ⊗ ) := σ (π(1)) ⊗ · · · ⊗ σ (π(k)) for any π ∈ S k and σ ∈ (R d ) k . For any σ ∈ (R d ) k , we define the symmetric tensor product by The space spanned by these symmetric tensors is then denoted by (R d ) ⊙k , and is a vector subspace of (R d ) ⊗k . The natural scalar product on (R d ) ⊗k and (R d ) ⊙k is denoted by A : B for A, B ∈ (R d ) ⊗k and, as usual, is defined to be the linear extension of In particular, for u :

General Results for the Linearised Equation
At the most fundamental level, our results concern the characterisation of the far-field behaviour of lattice displacements u : Λ → R N that are close to equilibrium in the far-field. More precisely, given a stable Hamiltonian H as introduced in the previous section we will characterise the decay of a general lattice displacement u provided that the (linearised) residual forces f (ℓ) := H[u](ℓ) decay sufficiently rapidly as |ℓ| → ∞. In § 2.3 we will then show how to use these results for the linearised operator to obtain characterisations of the far-field behaviour of equilibrium displacements in our full nonlinear interaction model.
With this definition we have the following result. Remark 2.2. If H[u] satisfies the assumptions for some α < −1 instead of α ∈ N 0 , then there is no logarithmic factor needed in the result, i.e., (10) The same holds true for the Theorems 2.4 and 2.6 below, if α < −1.
Remark 2.3. When initially reading both this theorem and the theorems in §2 and §3, we suggest to ignore the logarithmic terms and focus only on the algebraic rates. However, the treatment of the logarithmic terms is an important aspect of both our theorems and proofs since they appear to be intrinsic to the expansion, and not due to suboptimal estimates.
Suppose now that we have a general elastic field u with non-vanishing moments (8) but still fast decay of the residual forces. Then, we can decompose u into a truncated multipole expansion -higher order derivatives of the lattice Green's function defined in § 2.1 -corresponding to the nonvanishing moments and a far-field remainder that exhibits the improved decay established in Theorem 2.1. This idea is made precise in the next result.
and the remainder decays as It is convenient to have a continuum reformulation of this multipole expansion, to avoid having to work with the discrete Green's function and its discrete derivatives. Towards that end we exploit the connections between continuum and discrete Green's functions and derive higher order continuum approximations of the discrete Green's function. Specifically, in § 6 we will construct a sequence of continuum kernels with the following properties.
Theorem 2.5. There are unique kernels G n ∈ C ∞ (R d \{0}; R N ×N ) such that for all ℓ ∈ Λ and j, p ∈ N 0 and such that the G n are positively homogeneous of degree (2 − 2n − d) if n ≥ 1 or d ≥ 3, while in the case n = 0, d = 2 we have G 0 (ℓ) = A log|ℓ| + ϕ(ℓ), where A ∈ R N ×N and ϕ is 0-homogeneous. Furthermore, we find G 0 = G C .
The G n , n > 0, are higher order corrections resolving the atomisticcontinuum error. We will give precise definitions of all the G n in § 6. In addition we want to point out that the G n are practically computable via Fourier methods. A formula for that is also given in § 6.
Returning to the multipole expansion, if p = 1, 2 in Theorem 2.4, then one can replace the lattice Green's function G with the continuum Green's function G C . However, for a higher order continuum description, higher order G n need to be used. If we also use Taylor expansions to get actual derivatives we get a pure continuum expansion.

The Full Far-Field Expansion
Theorem 2.4 lays out a path on how to construct good far-field approximations for a solutionū of the atomistic equations δE(ū) = 0. Instead of looking at the solutions directly it suffices to construct an approximate solution of the equationsû such that the remainder r defined byū =:û + r has small linearized forces H[r] in the far field. For this approach to be useful, it is desirable thatû is both easy to understand analytically and practically computable.
Our goal is to construct smooth continuum approximations through the addition of successive corrector terms u C i up to the desired order. We writē u =û p + r p = u C 0 + u C 1 + · · · + u C p + r p and aim to achieve |H[r p ]| |ℓ| −d−p−1 0 , so that with Theorem 2.4 we havē ; that is, the remainder is highly localised around the defect core.
The precise statements for both point defects and screw dislocations are given in § 3. The full rigorous construction of the u C i is given in the proofs in § 7. We do however want to formally outline this construction here.
The first step is a Taylor expansion of the energy around the lattice or, more precisely, the potential V around zero, giving Separating out the linear terms and inserting the ansatzū =û p + r p = u 0 + u 1 + · · · + u p + r p gives The next step then is to Taylor expand the discrete differences in H[u i ] and δ k+1 E(0)[û p ] k leading to continuum differential operators. In particular, the leading order term for H[u i ] is c vol H C [u i ] followed by higher order differential operators. Formally, u i is of the order |ℓ| 2−d−i 0 with each derivative or discrete difference adding one order of decay. This allows us to group all the resulting terms into S i , based on their order of decay. We thus obtain We can therefore use the PDE to define u C i . When we look at all the details of this construction in § 7, we will see that the u i as used on the right hand side of (15) has to include the multipole terms of that order, which we can write as In particular, it is then important to know that with the help of Theorem 2.6 we can use the their smooth continuum variant instead of the discrete one. The S i in (15) are in general non-linear and higher order differential operators. Crucially though, S i only depends on the terms u 0 , . . . , u i−1 and not u i itself. That means that the equation defining u i is always the same second order, elliptic continuum linear elasticity equation With u C i defined this way, most of the terms on the right hand side of (14) cancel out and a precise estimate of the higher order errors gives the desired estimate for H[r p−1 ].
In § 7, we give a precise definition of the S i in the corrector equation (15). The grouping of the terms into different S i depends on the dimension d. As an example, for d = 2, the first three S i are given by is the Cauchy-Born energy density for A ∈ R N ×d and H SG is a linear differential operator describing a strain-gradient term in linear elasticity. It is defined by If on the other hand d = 3, then

Far Field Expansion for Crystalline Defects
We now demonstrate how our general structural results can be directly applied to obtain precise characterisations of the discrete elastic far-fields surrounding crystalline defects. Our results apply directly to points defects and screw dislocations; edge and mixed mode dislocations require additional ideas due to their more non-trivial lattice topology and will therefore not be discussed here. In addition, we briefly outline how these characterisations give rise to novel algorithms for simulating such defects.

Point defects
We consider point defects first. We briefly review the setting of [EOS16] to motivate the formulation of our main result in this context. First, we assume that the point defect has a reference configuration Λ def ⊂ R d , d ≥ 2 which is locally finite and homogeneous outside some defect radius R def , meaning Let R ℓ denote a finite interaction range for each site x ∈ Λ def and assume that there is a family of site energies V ℓ ∈ C K (R dR ℓ ), ℓ ∈ Λ def . Moreover, assume N = d and that there is a homogeneous interaction range R and site energy V ∈ C K (R dR ) for all sites of Λ, and that The potential energy under a displacements u : Λ def → R d and u : Λ → R d are then, respectively, given by With these definitions E def is then well defined on H c and have a unique Allowing Λ def = Λ in B R def admits defects such as vacancies and interstitials, while allowing inhomogeneity of V ℓ admits impurities and foreign interstitials.
A point defect can be thought of as a finite-energy equilibrium of E, that is, equilibrium displacementsū def ∈ H 1 (Λ def ) such that A notationally convenient approach is to simply projectū def to the homogeneous lattice. That is, we defineū : This is of course only one of many possible projections, which we have made only for the sake of notational convenience. Our subsequent results are essentially independent of how this projection is performed. Most importantly, because we haveū =ū def outside the defect core we obtain that for |ℓ| large enough. Here δ ℓ (ℓ ′ ) := δ ℓℓ ′ . With a small amount of additional work one can in fact show that This motivates the setup for our next result.
and such that the u C i satisfy the PDEs in equation Remark 3.2. In particular, up to order p = d, we have the pure multipole expansion Effects from the nonlinearity and higher order derivatives are only noticeable in terms beyond that.
Remark 3.3. In the point defect case it is likely possible to reduce the number of logarithms in the estimate somewhat for all orders. We do not want to explore this in detail but want to point out the log-free estimate for low orders which directly follows from (23) and estimates on the lattice Green's function based on Theorem 2.5. To be precise, for p < d we havē

Screw dislocations
Now let us consider screw dislocations. Again, our modelling follows the setup in [EOS16] and [BBO19]. We consider a straight screw dislocation with periodic behaviour along the dislocation line so that we can project to the lattice to a two-dimensional lattice on the normal plane to describe the behaviour. Hence we have d = 2 and N = 3 meaning u : Λ ⊂ R 2 → R 3 , though N is left arbitrary in the following to include for example the case in [BBO19] where N = 1.
Again we have a finite interaction range R and a site energy V ∈ C K (R dR ). The potential energy is then formally given by However that sum will usually not converge so we follow [EOS16] and consider the energy differences instead, where u CLE is the continuum linear elasticity solution.
More precisely, u CLE solves where b ∈ R 3 , b e 3 , is the Burgers vector of the screw dislocation,x ∈ R 2 is the reference position of the dislocation core and Γ := {x ∈ R 2 : x 2 = x 2 , x 1 ≥x 1 } a branch-cut chosen such that Γ ∩ Λ = ∅. We want to point out that the precise positioning ofx is not crucial and does not have physical meaning as the difference between two shifted solutions is in the energy space H 1 . Equation (29) was missed in [EOS16] but is in fact crucial for the results there to be true. It encodes the assumption that the system has zero net force and thus avoids spurious solutions of the type g(x) = u CLE (x) + G C (x −x). However, the standard construction of a solution u CLE , which can be found, e.g., in the book by Hirth and Lothe [HL82], already takes it into account. Therefore, the results of [EOS16] and later works that build on it remain correct provided such a solution u CLE is employed.
The following observation links it to the atomistic setting.
Therefore (29) is equivalent to either of these sums vanishing.
Indeed, the property δE(u)(ℓ) = 0 is heavily used in [EOS16] and we will use their results here.
With the definition (24) the energy E is then well defined on u CLE + H c and has a unique continuous extension to u CLE + H 1 . And with V ∈ C K we also find E ∈ C K , see [EOS16, Lemma 3]. (30) and such that the u C i satisfy the PDEs (59) and u C 0 = u CLE . Furthermore, the remainder r p+1 satisfies the estimate Remark 3.6. Contrary to the point defect none of these terms are expected to vanish in general, except for a few special cases which are explored in [BBO19]. In particular, the regularity assumption cannot be weakened as in the point defect case. Indeed, our general theory without looking at any special cases requires K ≥ J +2+⌊ p d−1 ⌋. As far as our proof goes the number of logarithms is optimal for d = 2, though probably not for higher dimension. We also expect this to be generic for the theorem itself as indeed the u C i will (in general) contain higher and higher logarithmic terms. However, in special cases these logarithmic terms in the u C i do not necessarily always appear, as explored in [BBO19].

Accelerated Convergence of Cell Problems
An immediate application of the defect expansions of Theorems 3.1 and 3.5 is that they suggest a novel family of numerical schemes that exploit these expansions to accelerate the simulation of crystalline defects. Here, we will only sketch one such scheme, but leave a more detailed analysis for future work.
Consider the equilibration of a point defect or a screw dislocation near the origin as in Theorems 3.1 and 3.5. We define a family of restricted displacement spaces where atoms are clamped in their reference configurations outside a ball with radius R. Then we can approximate (22), (30) by the Galerkin projection Under suitable stability conditions it is then shown in [EOS16] that for R sufficiently large, where p ∈ {0, 1}. This convergence is an almost immediate corollary of the decay estimate |Dr 1 (ℓ)| |ℓ| −d 0 log p |ℓ| 0 . (For energy minima, [EOS16] can be applied directly while for saddle points the analysis of [BDO20] can be readily adapted.) Our aim now is to accelerate this relatively slow convergence by providing an improved far-field boundary condition.
The overarching principle is to 1. replace the naive far-field predictor with the higher-order predictor 2. and to enlarge the admissible corrector space with the multipole moments That is, the corrector displacement is now parametrised by its values in the computational domain B R ∩ Λ and by the coefficients of the multipole terms.
We can then consider the pure Galerkin approximation scheme: The arguments of [EOS16] leading to (34) are generic Galerkin approximation arguments, leveraging the strong stability condition. They can be followed verbatim up to the intermediate result (Céa's Lemma) The existence ofv R is implicitly guaranteed through an application of the inverse function theorem, due to the fact that the right-hand side in this estimate approaches zero as R → ∞. To estimate the right-hand side we can insert the exact tensors b i,k from the solution representation of Theorems 3.1, 3.5 into v R , in order to obtain Dū − Dû p − Dv R = Dr p+1 − Dw R , where r p+1 is the core remainder term, and hence We can now define w R to be a suitable truncation of r p+1 to the computational domain B R . The details are given in [EOS16, Thm. 2] and immediately yield the following result.
Theorem 3.7. Suppose thatū is a strongly stable solution of (22) or (30); that is, there exists a stability constant c 0 > 0 such that then, for R sufficiently large, there also exists a solutionū R ∈ U Remark 3.8. The scheme (35) cannot be implemented as is since the energy difference functional cannot be evaluated for a displacement with infinite range. However, this highly idealised scheme is of immense theoretical value in that it highlight what could potentially be achieved if this challenge can be overcome. Any practical scheme will necessarily have to engage in the approximate evaluation of the multipole tensors b (i) , for which there are several promising possibilities that we will explore in separate works.
A second challenge for practical implementations is the fast and accurate evaluation of the higher-order far-field predictorû p . All of these approximations require suitable controlled approximation to E, somewhat analogous to quadrature rules or other kinds of variational crimes in the classical numerical analysis context.

Conclusions and Outlook
The main result of the present paper is the fact that the elastic field surrounding a defect in a crystalline solid may be represented to within arbitrary accuracy with three "low-dimensional" ingredients: 1. a series of continuum fields specified through PDEs; 2. a series of multipole moments; and 3. a highly localised discrete core correction.
More specifically, we have shown that by increasing the accuracy of components (1) and (2), the core correction (3) becomes increasingly local. While there is a certain amount of interaction between the components (1) and (2), there is no coupled problem that needs to be solved at any point. Indeed, both series are obtained sequentially order by order and the PDE defining the term of a given order in (1) only depends on lower-order terms of the multipole expansion, but not on the multipole term of the same order.
Our presentation here is restricted to simple lattices and a limited class of defects. Generalisations do require additional technical difficulties to be overcome, but there appears to be no fundamental limitation to extend the method and the results to multi-lattices and a range of other defects in some form. To conclude, we briefly discuss some of these possibilities and limitations.
• Edge and mixed dislocations: Edge and mixed dislocations are technically more challenging as they create a mismatch that affects the two-dimensional reference lattice. To leading order it suffices to correct the CLE solution with an ad-hoc transformation u 0 = u CLE • ξ −1 , see [EOS16], though the analysis also becomes a bit more technical still. For higher orders however more care needs to be taken, not just in the choice of ξ but also in the effect such transformations have on the PDEs. Furthermore, many arguments have additional technical complications due to the need of slip operators to describe the elastic strain.
• Cracks: The full extension of our results to crack geometries appears to be considerably more challenging, as the homogeneous lattice is no longer a particularly good global reference. Thus already the discussion of Green's function is significantly more involved [BHO19]. Formally, we still expect our overall strategy to apply and it is interesting to note that due to different orders of decay the first higher order corrector is already needed to even define the model in the first place rather than "just" improve on it.
• Energy differences for defect transitions: The precise characterisation of the far-field strain in terms of defect continuum fields and a multipole expansion suggests that some level of cancellation in energy differences, e.g. between a saddle point and energy minimum as observed in [BDO20], could be precisely tracked and characterised. Moreover, such results may then also explain improved convergence rates of numerical schemes for energy differences that are often seen in practice.
• Convergence of numerical schemes: A consequence of our analysis, with direct practical value is the construction of improved approximate cell problems that leverage the explicit low-dimensional structure of defect fields that we identified. We have given a hint at how this might be achieved in Theorem 3.7 but much additional work is needed to formulate practical schemes along these lines.
The same line of work can also lead to robust new numerical schemes and analysis of existing schemes for the defect dipole tensor specifically (also called the elastic dipole tensor). Such schemes are of important and ongoing interest in defect physics (e.g., [NMP + 16], [DM18]). In particular, our approach to these terms developed here naturally includes the anisotropic case as well as extensions to higher multipole tensors.
• Higher-order dislocation dynamics models: A further consequence is that we hope that the expansion of the far-field strain we have obtained here allows us to go beyond traditional dislocation dynamics approaches, which rely upon the leading-order CLE description of these defects. By using the structure of our expansion to studying the effect of applied stress fields on defect cores, we can provide more detailed atomistic input into such models. This suggests a route to better connect dislocation dynamics and atomistic approaches, bridging the scale and language gaps between these two simulation methodologies.
• Dynamics and statistical mechanics: Statistical mechanics models, such as free energies or transition rates could in principle benefit from an analysis within our new framework. For example, in the harmonic approximation, the analysis of [BDO20] could be taken as a starting point. It is far less clear whether more nonlinear models could also benefit, and it appears certain that finding similar coarse-grained descriptions of full dynamics of a crystalline far-field would require very different ideas.

Proofs -Decay Estimates
In this section we want to prove Theorem 2.1 and Theorem 2.4. But first, let us cite the following lemma.
Proof. For α = 0 this is [EOS16, Corollary 1]. However the addition of logarithmic terms is trivial. Indeed, one can construct g in the exact same way and can easily carry the logarithmic term through by including them in the weighted norms used in the proof.
The decay estimate for the lowest order found in [EOS16] is indeed based on Lemma 5.1 which then allows for a partial summation in the Green's function representation sum of the remainder given in equation (37). The key idea for the higher order decay estimates is that Lemma 5.1 can actually be extended to higher orders based on vanishing higher order moments as long as one only tries to write symmetric parts in divergence form; see Proposition 5.4 below. We will then use this higher order divergence form with more precise higher order partial summation in specific parts of the lattice in the Green's function representation sum of the remainder given in equation (37) to arrive at the new decay estimates. We begin by establishing two further auxiliary results.
Lemma 5.2. Given a linearly independent set of vectors S ⊂ R d (which must necessarily have #S = k ≤ d), the set of tensors {σ ⊙ : σ ∈ S k } is also linearly independent in (R d ) ⊗k . Furthermore, σ ⊙ = ρ ⊙ with σ, ρ ∈ S k , if and only if ρ = π(σ) for some permutation π ∈ S k .
Proof. Although the proof is straightforward, and likely well-known, we present it for the sake of convenience: Define a scalar product (·, ·) S on R d for which S forms part of an orthonormal basis. This induces a scalar product on the space of tensors by multi-linear extension of (σ ⊗ , ρ ⊗ ) S := k j=1 (σ (j) , ρ (j) ) S .
Consider the scalar product If ρ = π(σ) for all π ∈ S k , then for each π ∈ S k , there exists an index i such that σ (i) = ρ (π(i)) , and consequently (σ (i) , ρ (π(i)) ) S = 0. This entails that each of the products summed in the expression above is zero. Since ρ ⊙ and σ ⊙ are both non-zero and orthogonal in this inner product, they cannot be equal.
The same argument entails that if σ ⊙ = ρ ⊙ , there must exist π ∈ S k such that ρ = π(σ). On the other hand, if ρ = π(σ), we see that We deduce that σ ⊙ and ρ ⊙ are either identical (if and only if σ is a permutation of ρ) or mutually orthogonal in the S-inner product, which implies the stated result.
Lemma 5.3. Let σ ∈ Λ p , then Proof. This identity follows from two observations: First, for a polynomial r in ℓ of degree at most j − 2. Secondly, for any such polynomial. A simple induction then shows the result.
We can now turn towards to a crucial result converting a discrete force field into higher order divergence form. We will slightly abuse notation in the following and let the symmetric part sym A of a tensor A ∈ R N ⊗ (R S ) ⊗k denote the symmetrical part in the later indices only, namely sym(A) := 1 k! π∈S k A ·,π(σ) . For higher order tensor fields we will follow the usual convention from the continuum that Div always applies to the last component.
Proof. p = 1 is covered by Lemma 5.1. By induction, assume the statement is true for a p ∈ N and we now have q−p > d as well as p+1 vanishing moments. In particular, we have already constructed the desired f (0) , . . . , f (p) . Now | sym f (p) (ℓ)| |ℓ| −(q−p) 0 and q − p > d.
We claim that for all j = 0, . . . , p. This is clear for j = 0. By induction, if it is true for j < p, then Note that the decay of f (j+1) is needed not just for the sums to exist but also for the partial summation to be true. This establishes (36). Applying Lemma 5.3, with j = p, and using that I p = 0 we obtain According to Lemma 5.2, the set of the tensors sym σ ⊗ is linearly independent. Additionally, ·,ρ as ρ = π(σ) for some permutation π. Therefore, ℓ sym f (p) (ℓ) = 0 and we can apply Lemma 5.1 to find f (p+1) with the desired properties.
To prepare for the Proof of Theorem 2.1, we fix some u ∈ H 1 and write f (0) := H[u]. We extend the lattice Green's function approach developed in [EOS16] to estimate u. The Green's function satisfies for all u ∈ H c . As the right hand side is not invariant under adding a constant to u, this cannot directly translate to general u ∈ H 1 , but the situation looks better for derivatives.
Lemma 5.5. Let u ∈ H 1 and assume that |f (0) (ℓ)| = |H[u]| |ℓ| −γ 0 for some γ > 1. Then, for all ρ = (ρ 1 , ..., ρ j ) ∈ R j , j ≥ 1, we have Proof. Due to the decay assumption on f (0) the sum converges absolutely. Furthermore, for u ∈ H c the statement is clearly true. The right hand side is a well-defined, continuous, linear functional on H 1 . The result is straightforward when d ≥ 3 or j ≥ 2 as in this case the left hand side is also a continuous, linear functional because |D ρ G k | ∈ ℓ 2 (Λ) if and only if 4 < d + 2j.
To include the case where d = 2 and j = 1 we have to be a bit more careful. Note that, for u ∈ H c z∈Λ As the right hand side is now a well-defined, continuous, linear functional on H 1 , we find for all u ∈ H 1 , that If we now have a u ∈ H 1 where additionally |f (0) | |ℓ| −γ 0 , then the sum z∈Λ f (0) (z)D ρ G k (ℓ − z) converges (absolutely). Consider a smooth cutoff function η R , such that η R : R d → [0, 1] satisfies η R (z) = 1 for |z| ≤ R and η R (z) = 0 for |z| ≥ 2R, as well as |∇ k η R | R −k for k = 1, ..., j. With that, we see as the remaining term can be estimated by That is, (37) holds for all u ∈ H 1 with |f (0) (ℓ)| |ℓ| −γ 0 .
Proof of Theorem 2.1. We now use the Green's function representation (37) to estimate D ρ u(ℓ). Let us consider |f (0) | |ℓ| −d−p 0 log α |ℓ| 0 where we include both the case α ∈ N 0 and α < −1 to include Remark 2.2. Then let us assume there are p vanishing moments.
Although our argument is related to the lower order equivalent in [EOS16], it is technically more complex as both the higher derivatives and the higher z = ℓ z = 0 order divergence form lead to partial summations which, crucially, have to be performed on only specific and separate parts of the lattice. We will therefore provide full details.
This can be done in a very clean way by splitting the sum (37) into four regions; see Figure 5 for a visualisation: Region 1 is the far field, where |z| and |ℓ − z| are comparable. Region 2 is the intermediary area where |z|, |ℓ|, and |ℓ−z| are all comparable. And region 3 and 4 are the areas around z = 0 and z = ℓ, where either |z| can be small but |ℓ| and |ℓ − z| are comparable or |z − ℓ| can be small but |ℓ| and |z| are comparable. Inserting estimates for the residual f (0) and the lattice Green's function G and using this split of the sum indeed gives sharp estimates in absolute value. However, as we will show below, this is only a good estimate if both p = 0 and j = 1, 2. If either p ≥ 1 or j ≥ 3 then the sum in (37) exhibits significant large scale cancellation effects. To get sharp estimates in these cases, we will remove these cancelling terms via separate partial summations in region 3 and 4. The required discrete derivatives are directly available in the case j ≥ 3 or are obtained with the help of Proposition 5.4. To avoid discrete boundary terms, we will not split the four regions sharply but use smooth cutoff functions. The boundary terms in the partial summation are then spread out and can be treated like terms in region 2.
We first estimate the far-field term where the logarithmic term is estimated trivially by log α |z| 0 ≤ log α |ℓ| 0 for negative α, while for α ∈ N 0 the estimate instead follows from partial integrations of the resulting one-dimensional radial integral ∞ |ℓ| r 1−j−d−p log α r dr. The intermediary area is even more direct as we can just estimate the functions uniformly and multiply by the number of lattice points in the area Next, T 3 can be estimated by Hence, we have |T 3 | |ℓ| 2−d−j 0 if either p > 0 or α < −1. If, on the other hand, p = 0 and α ∈ N 0 , then we obtain |T 3 | |ℓ| 2−d−j 0 log α+1 |ℓ|. Finally, for T 4 the estimate is Putting all four estimates together, this completes the proof in the special case where both p = 0 and j = 1, 2.
For p ≥ 1, we need a better estimate on T 3 . We choose f (m) according to Proposition 5.4. We claim that for 0 ≤ m ≤ p Let us prove (40) by induction over m. Clearly, it is true for m = 0 since the second term on the right-hand side is identical to the left-hand side. Given its validity for a m, with m + 1 ≤ p, we now employ Proposition 5.4 and summation by parts to obtain The last term is concentrated in the annulus where D τ η is non-zero and can be estimated by which completes the proof of (40). Using (40) for k = p, we find Therefore, |T 3 | |ℓ| 2−j−d−p 0 for α < −1 and |T 3 | |ℓ| 2−j−d−p 0 log α+1 |ℓ| 0 for α ∈ N 0 . This finishes the proof for j = 1, 2 and arbitrary p.
Finally, for j ≥ 3, we need a better estimate for T 4 . In this case, let us split the difference directions as ρ = (σ, τ ), with σ ∈ R 2 . Then, Note that if at least one discrete derivative falls on η ℓ the estimate is similar to T 2 as |z|, |ℓ|, |ℓ − z| are then comparable for all non-zero terms. We thus get Note that we freely used the estimates of |D j−2 S f (0) | for |D j−2 −S f (0) |, as S spans the lattice and thus these estimates are equivalent.
Theorem 2.4 is almost an immediate consequence of Theorem 2.1. The only additional ingredient we need for the proof is the following lemma.
Proof. For any b one calculates Note that H[D ρ G k ] = D ρ δ 0 e k . In particular, there are no summability problems as this expression has compact support. With that in mind we can use Lemma 5.3 to calculate Therefore, it is sufficient to show that the linear map is a bijection. By dimensionality (recall that S is chosen to be a basis), this is equivalent to T being one-to-one which follows from Lemma 5.2.
Proof of Theorem 2.4. k) : D i S G k , and choose b according to Lemma 5.6. As the moments are linear, the result then follows directly from Theorem 2.1.
Note that even though the choice of b in Lemma 5.6 is unique, that does not necessarily mean that the choice of b for Theorem 2.4 is unique. Indeed, it can happen that for certain b the sum N k=1 b (i,k) : D i S G k exhibits sufficient cancellation to be part of the error estimate, for example if the coefficients b correspond to a discrete approximation of the CLE equation.
Decay estimates for the far-field continuum equations We also want to have decay estimates for the equations (15) introduced in § 2.3. These can be obtained along similar lines to the discrete case. However, we are purely interested in the far-field which simplifies things a bit and avoids a few additional complications specific to the continuum case.
Remark 5.8. It is important to point out that the additional logarithm in the result is a generic aspect of the problem and not just a limitation of our proof. To see that consider α = 0 and p = 1. Specifically, let us take a f with f = div g, |g(ℓ)| |ℓ| −d , and |f (ℓ)| |ℓ| −d−1 . Then our proof below actually shows Even for something as simple as g(ℓ) := E 11 |ℓ| −d one directly gets a logarithm from the integral. That means the logarithm term naturally comes from a summing of the stresses (for p = 1) around the defect core.
Of course there are special cases where there is no logarithmic term appearing. For p = 1, α = 0 that would be any setting with enough cancellation such that sup For example, the first corrector equation in [BBO19] satisfies ∂B R (0) g(z) dz = 0 and indeed there was no logarithmic term.
Sketch of the Proof. As the equation only needs to be satisfied outside of B R 0 , we can assume without loss of generality that f ∈ C ∞ (R d ; R N ), by multiplying f with an appropriate cutoff function. Our ultimate goal will be to define u := f * G C and follow the same estimates as in the discrete case. The main problem in adapting the discrete proof to the continuum is the need to replicate Lemma 5.1 and Lemma 5.4 which allow conversion to appropriate divergence form for p ≥ 1, and we sketch out the route to the analogues of these results here. We follow the general idea of the proof in [EOS16], defining f 0 := f , f n+1 (ℓ) := 3 d f n (3ℓ), and By direct computation it can be checked that and using the bound on f assumed in the statement and p ≥ 1, one concludes that f n → 0 uniformly in any R d \B δ (0). Similarly, we may show that the sequence g n is uniformly summable over the same sets, and we define g := ∞ n=0 g n . Using the relation (41), we find f = − div g in R d \{0}. It also easy to check that g ∈ C ∞ (R d \{0}) and that g satisfies the decay |g(ℓ)| |ℓ| −d−p+1 log α |ℓ| away from zero. In the continuum setting, a possible consequence of this construction is that g may have a singularity at ℓ = 0. However, since we are only interested in the far-field behaviour, we can simply remove any singularity by redefining g on the interior of B R 0 (0) such that g ∈ C ∞ (R d ; R N ×d ), and then redefine f := div g. This operation does not affect the behaviour outside B R 0 (0), and we have therefore recovered the equivalent of Lemma 5.1.
A similar construction and inductive argument to that used in the proof of Lemma 5.4 allows us to generate an extension of |ℓ| −d−p+k log α |ℓ| for ℓ bounded away from zero. To avoid generating a singularity at ℓ = 0 in this case, again we redefine the functions in B R 0 (0) starting at the highest order f (p) .
With this ability to transform f to (higher-order) divergence form, one can now follow the proof of Theorem 2.1 almost verbatim to estimate the derivatives of u = f * G C using the split into the four different contributions defined in (39). Note that the singularity of G C is not a problem as G C and its first derivative are locally integrable everywhere, and we only need to take further derivatives of the Green's function in the estimate of T 3 , which is the part away from that singularity where G C is smooth. Indeed, in the estimate of T 4 , we can use the fact that we have assumed estimates for derivatives of f of arbitrary order and only put one derivative on the Green's function instead of two.

The lattice Green's function and series expansion
In this section, we define the lattice Green's function, and construct a series expansion in terms of continuum kernels proving Theorem 2.5. The general strategy is to expand the inverse of the discrete Fourier multiplier corresponding to the lattice operator H as a series of homogeneous terms about k = 0, and then use tools developed by Morrey in [MJ66] to provide integral representations of these terms in real space which allow us to explicitly characterise the decay and even homogeneity of each term. In the first parts we will also follow some ideas that were developed in [MR02] for the simpler scalar case where N = 1.

Fourier transforms
We begin by providing definitions of the Fourier transforms we deploy here. For f : Λ → R n , we define the Fourier transform F a [f ] : B → R n , where B is the Brillouin zone of the lattice, via The inverse is given by Note that |B| = (2π) d c vol , where c vol is the volume of a fundamental cell of Λ. We note the following properties.
if f decays fast enough such that the sum converges absolutely.
We also use the standard continuum Fourier transform defined for f ∈ L 1 (R d ; R N ) and extended to the space of tempered distributions S ′ (R d ; R N ) by duality. As usual the inverse is given by

Fourier multiplier series manipulation
Next, we formally carry out a series development of the inverse of the discrete lattice operator Fourier multiplier, which we will subsequently analyse in detail.
As before consider the Hamiltonian H = δ 2 E(0) given pointwise by and we write C σρ = (C σiρj ) ij ∈ R N ×N . C satisfies the symmetry C iρkσ = C kσiρ as a second derivative and C iρkσ = C i(−ρ)k(−σ) due to equation (2). Then (see [HO12,BS16]) H has a representation as a Fourier multiplier where the matrices A ρ are implicitly defined through the second identity. As before we assume lattice stability (4). In reciprocal space this entails that there exists c 0 > 0 such that where Id is the N × N identity matrix. We note that H ∈ C ∞ per (B), and by virtue of lattice ellipticity, H is strictly positive definite except when k = 0. This entails that the matrix inverse H −1 (k) exists everywhere in B except at k = 0, and moreover, We define and re-summing, we may express H −1 as a series of terms with increasing homogeneity: Explicitly, A 2n (k) is defined by the finite sum This means A 2n is 2n-homogeneous. In cases where the first term in the expansion of H 2 is a multiple of the identity, i.e. H 2 (k) = h 2 (k) Id, the series simplifies to and yet further simplification can be made if H is a multiple of the identity matrix, e.g., in the scalar setting N = 1, which is the case considered in [MR02].

The lattice Green's function and its development at infinity
With the previous series development complete, we may introduce and analyse the lattice Green's function, which acts as an inverse to the lattice operator H. If d ≥ 3, we define the lattice Green's function to be and if d = 2, we set where γ is the Euler-Mascheroni constant. To demonstrate that the limit in the latter definition exists, we take 0 < δ < ε and calculate In the above argument, we have used thatĤ −1 − A −2 is bounded in a neighbourhood of k = 0, and the fact that |1−e iσℓr | = O(|ℓ|r) for all suitably small r. Rearranging, we have established a uniform bound independent of δ, and hence the limit exists. We now use the series expansion (42) to define a series of functions G n which approximate G. Note that A 2n−2 is locally integrable and defines a tempered distribution unless both n = 0 and d = 2. Except in this special case, we therefore define In the special case where n = 1 and d = 2 we define in the distributional sense, i.e.
With this definition a simple direct calculation also shows that and indeed G C = G 0 .

Alternative representation of G n
We now connect the definitions made above to alternative representations considered by Morrey in [MJ66, chap. 6.2]. These alternative representations will provide us with the ability to directly deduce regularity, decay, and even homogeneity. Furthermore, this representation is promising for computational uses as only a finite surface integral in Fourier space is required for its evaluation.
We begin by setting P = ⌈ d+2n−1 2 ⌉, and h = 2P + 2 − 2n − d, which we note satisfies h ≥ 1. Then, we define for ℓ = 0, where ∆ is the Laplacian, and We remark that these functions satisfy J ′ l (w) = iJ l+1 (w). We also use the following definition from [MJ66, Def. 6.1.4]. We directly want to point out though, that for the purpose of our results here we are only interested in the cases s < 0 and s = 0. So either ϕ is positively homogeneous of a negative degree, or ϕ(ℓ) = ϕ 1 (ℓ) + A log|ℓ|, where A is a constant and ϕ 1 is 0-homogeneous. With this definition in place, we state the following result.
Lemma 6.2. G M n is well-defined, smooth for ℓ = 0, and G M n is essentially Proof. This is part of [MJ66, Thm. 6.2.1]. Note that A 2n being matrixvalued does not change the argument.
In our setting, the only case that might fail to be positively homogeneous is G M 0 for d = 2. As pointed out above, in that case essentially homogeneous means that G M 0 (ℓ) = ϕ(ℓ) + C 1 log|ℓ|, where ϕ is positively 0-homogeneous and C 1 ∈ R N ×N .
We will now show that G n and G M n are in fact identical. Let us first give a short motivation for G M n . The basic idea of the definition (44) is to use the homogeneity of A 2n−2 to isolate the radial component in the Fourier integral and then integrate it out. Naïvely, this approach leads to a integral of the form ∞ 0 r l e irw dr which does not exist for any l and real w. The idea behind the P Laplacians is to ensure that the resulting l is non-negative. If one then adds a small imaginary part to w, the integral is indeed finite and can be computed. This leads to the J l . Lemma 6.3. On R d \{0} the distribution G n is represented by a function which we will also call G n and we have G n = G M n on R d \{0}.
Proof. As discussed in the motivation directly before the Lemma, we first let ℑ(w) > 0 and l ≥ 0. Then Take any δ > 0, and let [i, w] be the line segment in the complex plane, for which we note that z ∈ [i, w] implies ℑ(z) > 0. Define which is easily seen to satisfy (zϕ 0 (z)) ′ = e z −1 z . Using the definition of J 0 and the definition and property of ϕ 0 just noted, we then find where we define In (45), the first term is the desired integral for l = −1, the second term is a renormalization of the singular part for r close to 0, and the third term is of lesser importance and will go to zero at the end. We can simplify p δ 0 (w) further, as the exponential integral E 1 satisfies for z ∈ C\R ≤0 , where γ is the Euler-Mascheroni constant, and [z, ∞] is any contour in C\R ≤0 connecting z to positive real infinity. It follows that

Now, defining
for h ≥ 0 we may calculate inductively that This calculation was restricted to ℑ(w) > 0 to ensure that the integrals converge. But for h ≥ 1 that is true even for real w. So we now can take the limit ℑ(w) → 0 and see that for all w ∈ R, δ > 0, and h ≥ 1, we have that Using the definition in (44), we have obtained the following representation for G M n : We now inspect the three integrals in turn. Clearly, for the first in the sense of tempered distributions. For the second integral, as p δ h is a polynomial of degree h = 2P + 2 − 2n − d, we find (−∆) P p δ h (ℓσ) = 0 except when d = 2 and n = 0, in which case (−∆) P p δ h (ℓσ) = γ + log δ. The integrand in the final integral contains the factor (−∆) P δ(iℓσ) h+1 ϕ h (δiℓσ), which goes to 0 uniformly on compact sets as δ → 0. Since δ was arbitrary, we may pass to this limit, and we do indeed find G m = G M m on R d \{0}.
After having constructed the kernels G n and having established their properties, we can now prove the expansion of the lattice Green's function.
Proof of Theorem 2.5. We have already established the regularity and homogeneity properties of the G n . In the following, fix j ≥ 0 and ρ = (ρ 1 , . . . , ρ j ) ∈ R j . We then still need to show that for all ℓ, p. We claim it is in fact sufficient to show the result holds for any p but where the error decays with two orders less, i.e.
for all ℓ. If this estimate holds, then using equation (47) with p + 1 and the triangle inequality entails instead. If 2n − 2 > −d, note that A 2n η is integrable and has support in B, therefore as well as In all cases we find with as long as m is small enough to ensure that r m is integrable. The smoothness of η, and the properties ofĤ −1 and A 2n discussed in §6.2 ensure that r m is smooth away from k = 0 and is bounded by C|k| 2p−2m+j in a neighbourhood of k = 0. It follows that r m ∈ L 1 if 2m ≤ 2p + j + (d − 1) and therefore (48) implies For G 1−η n , we observe that D α A 2n−2 (1 − η) is integrable for all sufficiently large |α|, and therefore G 1−η n has super-algebraic decay at infinity. The same therefore also holds true for any D ρ G 1−η n . The only remaining claim now is the uniqueness. Assume we have two such series of kernels (G n ) n and (G n ) n . By induction, for p ∈ N 0 let us assume we already know that G n =G n for all n < p, we need to show that G p =G p .
First exclude the case where both p = 0 and d = 2. Then we know that G p andG p are positively homogeneous of order (2 − 2p − d). Using both estimates of order p without any derivatives, we obtain for all ℓ ∈ Λ. Combined with the homogeneity, we thus have Given any ε > 0 we can therefore find an R > 0 such that for all x ∈ { ℓ |ℓ| : ℓ ∈ Λ, |ℓ| > R}. As this set is dense in ∂B 1 (0) and both functions are continuous, (49) holds for all x ∈ ∂B 1 (0). As ε > 0 was arbitrary we thus have G p =G p on ∂B 1 (0) and by homogeneity on all of R d \{0}.
In the case that p = 0 and d = 2 the argument is only slightly more complex. We write G 0 = A log|ℓ| + ϕ andG 0 =Ã log|ℓ| +φ with A,Ã ∈ R N ×N and ϕ,φ both 0-homogeneous. Then for all ℓ ∈ Λ. As ϕ andφ are bounded, we can divide by log|ℓ| and send |ℓ| → ∞ to see that A =Ã. The identity ϕ =φ then follows the same way as before.
We can now rewrite the truncated multipole expansion of Theorem 2.4 in continuum terms.
Proof. This is an immediate consequence of the expansion of the lattice Green's function from Theorem 2.5 and a Taylor expansion of the discrete differences in D ρ G n .
As a consequence we obtain Theorem 2.6.
Proof of Theorem 2.6. The result follows by combining Theorem 2.4 and Lemma 6.4.

Proofs -Far Field Expansions for Crystalline Defects
In this section we prove the results of § 3. Before we get to the main results, let us start with a preliminary proof.
Proof of Proposition 3.4. δE(u)(ℓ) is summable according to [EOS16,Lemma 10]. The same then holds true for H[u](ℓ). The idea of the proof is to use a cutoff η R and sum by parts so that only the far away annulus B 2R \ B R contributes. Here, the linearisation and discretisation errors are small and we can go from the nonlinear forces to the linearised forces and all the way to the continuum forces. So, consider a smooth cutoff function η R , such that η R : R d → [0, 1] satisfies η R (z) = 1 for |z| ≤ R and η R (z) = 0 for |z| ≥ 2R, as well as |∇η R | R −1 . We then have as the linearisation error sums to O(R −1 ). Furthermore, as the discretisation error is also O(R −1 ). Now let us come to the main topic of this section, the proofs of Theorem 3.1 and Theorem 3.5.
At the outset, we recall from the setting discussed in §2 and §3 that we assume the site potentials (and hence the potential energy functional) are of class C K , the number of derivatives we want to estimate is J ≥ 2 and the dimension of the problem is d. We now prove the expansion and the error estimates by induction over p ≥ 0 as long as in the point defect case and in the case of the screw dislocation, where we always have d = 2. We also note that all results only need to be proven for |ℓ| ≥ R for some sufficiently large R.

The case p = 0
The case p = 0 is a consequence of results in [EOS16] for both the point defect case and the case of screw dislocation. To be precise, it follows from the assumptions of both Theorem 3.1 and Theorem 3.5 that J ≤ K − 2.
In this setting, we may apply Theorem 1 in [EOS16], which states that the solutionū in the point defect cases satisfies In the screw dislocation case, we may similarly apply Theorem 5 in [EOS16], implying that In order to combine and summarise these results, we define u 0 := 0 for the point defect, u CLE for the dislocation, so that (51) and (52) imply 7.2 The case p > 0 We now proceed inductively from the case p = 0.
Taylor expansion. We follow the construction that we formally motivated in § 2.3 and expand the solutionū by performing a Taylor expansion of the equilibrium equation for the forces for all ℓ sufficiently large, where g = 0 in the screw dislocations case and g has compact support for point defects.
Recall from §3 that E inherits the same regularity as the homogeneous site potential V , and is therefore C K . For such a variation of the energy we use the pointwise notation for the corresponding ℓ 2 representative. Specifically, that means A Taylor expansion around the homogeneous lattice u = 0 up to order K ≤ K − 2 gives 0 = Div g + δE(ū)(ℓ) (55) We will determine the precise orderK for this expansion as part of our arguments later. We note that we already have explicit expressions for the first two terms in the sum: A homogeneous (Bravais) lattice is always in equilibrium, δE(0) = 0, and the second term is δ 2 E(0)[ū] = H[ū]. For convenience, we will define T (1) p+1 to be the integral error term and also include the compactly supported Div g, i.e.

T
(1) so that we may rewrite (55) as Expansion of solution. Next, we make the inductive assumption that we have already obtained u i for 0 ≤ i ≤ p − 1, and that we can further decompose which are respectively continuum-and lattice-valued functions, so that we may writeū where r p is some remainder. These functions will be assumed inductively to satisfy certain properties which we now make precise: • We will assume that the former terms, u C i ∈ C ∞ (R d \B R 0 (0)), i ≤ p−1, have been constructed to solve a sequence of linear elliptic PDEs, and satisfy decay estimates We note that S i are forcing terms which we determine as part of our argument.
• We assume that the latter terms in (57) have been constructed to take the form for some constant tensors b (i,k) . In particular, H[u MP i ](ℓ) = 0 for |ℓ| large enough, and we assume u MP 0 = 0.
In light of (61), whenever useful, we may also rewrite the multipole contributions in terms of smooth continuum kernels by applying Lemma 6.4. Indeed, by combining the terms of equal homogeneity, Lemma 6.4 implies that we may write for all j ∈ N 0 , and we have the additional error estimate for all j ∈ N 0 . Using the functions u CMP , we therefore have the representationū Using the inductive assumptions made above, we may substitute (58) into the left-hand side of (56), and write The goal is now to split the remainder r p into the next continuum contribution u C p and a new intermediary term s p , and find the correct equation for u C p so the residual becomes even smaller in the far-field. Then we can apply Theorem 2.4 to s p to split this term yet further into multipole terms up to order p and a new remainder r p+1 with the desired improved decay.
Taylor expansion error. With the assumptions above, we turn back to the Taylor expansion (55), and establish the choice of order,K. Our aim is to fully resolve all orders of decay up to and including |ℓ| −d−p These orders of decay come from estimating a k-fold tensor product of Dū with the decay mentioned, and one further additional power arising from the discrete divergence. In particular, to ensure that the order of decay matches or exceeds that of |ℓ| −d−p 0 , we requireK ≥ p d−1 + 1, or in the case of a point defectK ≥ p−1 d + 1. Indeed, this is satisfied by assumption on K, if we setK := K − J − 1. We also note that the expression of T Proof. For ℓ large enough, − Div g(ℓ) = 0, hence we have In the case of a screw dislocation we can use the estimates (53) and (54) for u 0 and r 1 to see that |D jū | |ℓ| −j 0 for 1 ≤ j ≤ J. By distributing the outer discrete derivatives over the tensor product and using this estimate, we deduce that This proves the result in this case asK = K − J − 1 ≥ p + 1 and d = 2. 0 collects all the terms that contain contribution from s p or w MP p . To estimate derivatives of w MP p , we can use (64). For s p we can combine (62) and (67) to obtain for j = 1, . . . , J. This allows us to estimate the second error term.
Proof. We note that T (2) p+1 can be written as a sum of differences, For the k = 2 term, we use (68) and apply the discrete product rule to obtain We note that the shifts ρ appearing in the first estimate which arise from the discrete product rule remain below some finite radius R = R(J). Using the estimates for w MP p and s p gives the desired overall estimate. For k ≥ 3 a similar argument provides the same estimate, as any additional factor just improves the estimate by at least |ℓ| −1 0 log|ℓ| 0 1.
Continuum approximation. We want to construct a sequence of PDEs to resolve (69) to sufficient accuracy and thus describe the far-field behaviour ofū. To that end, we note that a finite difference of any smooth function v can be Taylor expanded to give for some θ ∈ [0, 1]. By iterating this process of Taylor expansion, a similar series approximation can be constructed for higher-order finite differences.
To start, we may use this expansion to approximate the action of H on a smooth function. We note that H was defined by (3) as The symmetry assumption (2) implies ∇ 2 V (0) ρρ ′ = ∇ 2 V (0) (−ρ)(−ρ ′ ) . Therefore all the odd orders cancel and using (71) we get that there exist tensors C 2n ∈ R N 2 ×d 2n for n ≥ 1 such that for any smooth function v, where T (3,1) collects the Taylor error terms. In coordinates, these tensors are Where we used that based on the definition of the Cauchy-Born energy density W in (6).
Construction of u C p . We recall that so far, we have Taylor expanded the site potential and substituted a truncated seriesū p on the right-hand side of (55) to obtain (69), i.e. We note that the remaining series on the left and right-hand sides of this equation involve only the action of discrete difference operators on smooth continuum functions, and so we may now perform a discrete-to-continuum approximation to handle these terms.
On the left, we insert (73) from the discussion of H C as the leading order into the full expansion (72) to see that there exist tensors C n ∈ R N 2 ×d n for n ≥ 3 such that for a smooth function v, We can perform a similar construction for the higher-order terms in the Taylor expansion on the right-hand side of (69), so that in general, there exist tensors C j in the natural tensor space which is isomorph to R N k+1 ×d n for each j = (j 1 , ..., j k ) that satisfies 1 ≤ j m ≤ M k − 1 and k m=1 j m = n for some n ≥ k + 1. This allows us to approximate δ k+1 E(0)[v ⊗k ](ℓ) as Each of the tensors C j is a sum of tensor products formed of components of ∇ k+1 V (0) and a n-fold tensor products of vectors from R. The final term T (3,k) (v) = T (3,k) ∇v(ℓ), . . . , ∇ M k v(ℓ), ∇ M k +1 v collects all the terms that include at least one Taylor error for v. And overall, after inserting the specific v, we can combine these errors to and for notational convenience in the following argument, defineũ i := u C i + u CMP i for i < p andũ p := u C p , which satisfies ∇ jũ i ≤ C j |ℓ| 2−d−j−i 0 log i |ℓ| 0 for all j ∈ N 0 while overall p i=0ũ i =ū p . Let us consider multi-indices i = (i 1 , . . . , k) for choosing multiple u i in a nonlinear term and j = (j 1 , . . . , j k ) for the number of derivatives on each term. We write |i| 1 = m i m for the sum of the components.
Note that as long as we choose the M k large enough, the definition of the S q , q ≤ p, in (78) is independent of their precise choice. Any additional terms from choosing M k even larger only appear in T (4) p+1 . Lemma 7.3. The following estimates holds: for all j ∈ N 0 .
In summary, we may use the definition of S q and the remainder terms T In view of the definition of u C i as a solution to (59), it follows that this equation reduces to the remainder equation which we study in the following section.
Series remainder estimate. Note that (82) is not a linear equation for s p , as s p also appears on the right hand side in T (2) p+1 . In the leading order estimates in [EOS16] this was resolved through techniques coming from regularity theory which we rely on for p = 0. For the higher orders, i.e. p > 0, we can use the (suboptimal) estimates that are available from the previous order. That turns out to be sufficient as the appearance of s p in T Conclusion. To conclude the induction step we still have to look at that the possible addition of lower order terms from (83) as they change the multipole terms from one iteration to the next. Instead of discussing under what conditions these terms vanish, we allow for them to be non-zero and consider the more detailed consequences. Indeed, this is fine as long as u MP 0 = 0 remains true as this is part of the results and as long as u CMP i , i = 1, . . . , p − 1, remain unchanged which is important since they are part of the continuum equations (59). The first just follows from the fact that I 0 [v] = 0 due to the divergence form of T p+1 combined with sufficient decay. For the second part, let us use the continuum expansion of the multipole terms to write We know that v CMP i is (2 − d − i)-homogeneous. Combining the decay estimate (70) for s p with the decay estimates for r p+1 andw MP p+1 , we find for |ℓ| large enough. This is only compatible with the homogeneity if v CMP i = 0 for i = 1, . . . , p − 1. That ends the induction step and thus the proof of Theorem 3.1 and Theorem 3.5.
As a last step we see that for a point defect, for any d, we have u C 0 = 0. We then also find u C 1 = · · · = u C d = 0 as S i = 0 for 1 ≤ i ≤ d. The first non-trivial terms arises when u CMP 1 appears in the nonlinearity, namely S d+1 = 1 2 div ∇ 3 W (0) ∇u CMP 1 2 .