{'fraction_non_alphanumeric': 0.052041654521847426, 'fraction_numerical': 0.014335303002199316, 'mean_word_length': 3.9546409347512146, 'pattern_counts': {'":': 1, '<': 0, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 32, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The role of the patch test in 2D atomistic-to-continuum coupling methods

For a general class of atomistic-to-continuum coupling methods, coupling multi-body interatomic potentials with a P1-finite element discretisation of Cauchy--Born nonlinear elasticity, this paper adresses the question whether patch test consistency (or, absence of ghost forces) implies a first-order error estimate. In two dimensions it is shown that this is indeed true under the following additional technical assumptions: (i) an energy consistency condition, (ii) locality of the interface correction, (iii) volumetric scaling of the interface correction, and (iv) connectedness of the atomistic region. The extent to which these assumptions are necessary is discussed in detail.


Introduction
Defects in crystalline materials interact through their elastic fields far beyond their atomic neighbourhoods. An accurate computation of such defects requires the use of atomistic models; however, the size of the atomistic systems that are often required to accurately represent the elastic far-field makes atomistic models infeasible or, at the very least, grossly inefficient. Indeed, atomistic accuracy is often only required in a small neighbourhood of the defect, while the elastic far field may be approximated using an appropriate continuum elasticity model.
Atomistic-to-continuum coupling methods (a/c methods) aim to exploit this fact by retaining atomistic models only in small neighbourhoods of defects, and coupling these neighbourhoods to finite element discretisations of continuum elasticity models; see Figure 1(c,d). By employing a coarse discretisation of the continuum model, such a process can achieve a considerable reduction in computational complexity, however, some of the first a/c methods suffered from the so-called "ghost force problem": While homogeneous deformations are equilibria of both the pure atomistic and the pure continuum model, they are not equilibria of certain a/c models [37,50,35] due to spurious forces -the "ghost forces" -that can arise at the interface between the atomistic and continuum regions.
Much of the literature on a/c methods has focused on constructing a/c methods that did not exhibit, or reduced the effect of the "ghost forces" [50,51,52,15,4,20,10,33,49,23,56,22]. A straightforward solution was the introduction of force-based (i.e, non-conservative) methods [25,50,51,10,20,33]. The construction of accurate energy-based coupling mechanisms turned out to be more challenging. Several creative approaches providing partial solutions to the problem were suggested [52,15,49,22], however, no general solution exists so far.
The inconsistency of early a/c methods is reminiscent of the inconsistency problems encountered in the early history of finite element methods. A simple criterion to test consistency of finite element methods is the patch test introduced by Irons et al [5]; see also [53,6]. The "ghost force problem" discussed above is precisely the failure of such a patch test. It is well known that, in general, the patch test is neither necessary nor sufficient for convergence of finite element methods; see, e.g., [53,6] where several variants of patch tests are discussed. The same is true for a/c methods: It was shown in [31] that a particular flavour of force-based a/c coupling typically has a consistency error of nearly 100%, even though it does pass the patch test.
(A notable exception is the work of Lu and Ming on force-based hybrid methods [30]. However, the techniques used therein require large overlaps and cannot accommodate sharp interfaces.) Only specific methods are analyzed; the question whether absence of "ghost forces" (or, patch test consistency, as we shall call it) in general implies satisfactory accuracy has neither been posed nor addressed so far. The purpose of the present paper is to fill precisely this gap. After introducing a general atomistic model and a general class of abstract a/c methods, and establishing the necessary analytical framework, it will be shown in Theorem 10, which is the main result of the paper, that in two dimensions patch test consistency together with additional technical assumptions implies first-order consistency of energy-based a/c methods.
1.1. Outline. §2 gives a detailed introduction to the construction of a/c methods. This section also develops a new notation that is well suited for the analysis of 2D and 3D a/c methods. In §2.3. 3 we give a precise statement of the patch test consistency condition. §3 contains a general framework for the a priori error analysis of a/c methods in W 1,pnorms, similar to an error analysis of Galerkin methods with variational crimes. Several new technical results are presented in this section, such as the introduction of an oscillation operator to measure local smoothness of discrete functions ( §3.2.1), an interpolation error estimate for piecewise affine functions ( §3. 2.2), and making precise the assumption made in much of the a/c numerical analysis literature that it is sufficient to estimate the modelling error without coarsening ( §3.5). The main ingredient left open in this analysis is a stability assumption, which requires a significant amount of additional work and is beyond the scope of this paper.
§4 presents two 1D examples for modelling error estimates, which motivate the importance of patch test consistency, and discusses the modelling error estimates that can at best be expected if a method is patch test consistent.
§5 introduces the two main auxiliary results used in the 2D analysis of §6: (i) Shapeev's bond density lemma, which allows the translation between bond integrals and volume integrals; and (ii) a representation theorem for discrete divergence-free P 0 -tensor fields.
Finally, §6 establishes the main result of this paper, Theorem 10: if an a/c method is patch test consistent and satisfies other natural technical assumptions, then it is also first-order consistent. The proof depends on a novel construction of stress tensors for atomistic models, related to the virial stress [2] (generalizing the 1D construction in [32,33]), and a corresponding construction for the stress tensor associated with the a/c energy. Moreover, we discuss in detail to what extent the technical conditions of Theorem 10 are required. For example, we show that, if the atomistic region is finite (i.e., completely surrounded by the continuum region) then the condition of global energy consistency already follows from patch test consistency.
1.2. Sketch of the main result. In this section we discuss the main result in non-rigorous terms. Let E a be an atomistic energy functional and let E ac be an a/c energy functional, which uses the atomistic description in part of the computational domain, and couples it to a finite element discretisation of Cauchy-Born nonlinear elasticity (cf. §2.3.2 for its definition), with suitable interface treatment. Suppose that the following conditions are satisfied: (i) E ac is patch test consistent: Every homogeneous deformation is a critical point of E ac (cf. § 1).
(ii) E ac is globally energy consistent: E ac is exact for homogeneous deformations (cf. § 2.3.2). (iii) Locality and scaling of the interface correction: The interface correction has (roughly) the same interaction range as the atomistic model, and has volumetric scaling (cf. §2.3.6). (iv) The atomistic region Ω a is connected.
(v) Stability: For some p ∈ [1, ∞], and for deformations y in a neighbourhood of the atomistic solution, the second variation δ 2 E ac (y) is stable, when understood as a linear operator from (discrete variants of) W 1,p to W −1,p (cf. §3. 3). (vi) The finite element mesh in the continuum region is shape regular. [7] Under conditions (v) and (vi), we show in §3 that the error between the atomistic solution y a and the a/c solution y ac can be bounded by ∇y a − ∇y ac L p (Ω) E model where Ω is the computational domain, Ω c the continuum region, h a local mesh size function, E ext h is the consistency error for the treatment of external forces, and E model ε is the modelling error, which we describe next. Since y a is not a continuous deformation, but a deformation of a discrete lattice, the various terms appearing above need to be interpreted with care. This is the focus of §3. For example, in the rigorous version of (1), we will replace ∇ 2 y a by an oscillation of a suitably defined gradient.
The step outlined above does not distinguish different variants of a/c methods. The error introduced by the coupling mechanism and the continuum model is contained in the modelling error E model where W −1,p ε is a suitably defined negative Sobolev norm on an atomistic grid. If E ac is not patch test consistent, then, typically, E model ε = O(1). However, if conditions (i)-(iv) hold, and if the problem is set in either one or two space dimensions, then we will prove in §6 that where ε is the atomistic spacing and Ω i an interface region. A rigorous statement of (2) is given in Theorem 10, which is the main result of the paper.
1.3. Basic notational conventions. Vectors are denoted by lower case roman symbols, x, y, a, b, and so forth. Matrices are denoted by capital sans serif symbols A, B, F, G, and so forth. We will not distinguish between row and column vectors. If two vectors are multiplied, then it will be specified whether the operation is the dot product or the tensor product: for a, b ∈ R k , we define where, in a ⊗ b, i is the row index and j the column index. When matrix fields take the role of stress tensors, we will also call them tensor fields and use greek letters σ, Σ, and so forth.
The Euclidean norm of a vector and the Frobenius norm of a matrix are denoted by | · |. The p -norms, p ∈ [1, ∞], of a vector or matrix are denoted by | · | p , and sometimes by · p . For ε > 0, the weighted p -norms, on an index set S , are defined as a p ε (S ) := ε 1/p a p (S ) . The topological dual of a vector space U is denoted by U * , with duality pairing ·, · .
If A ⊂ R d is a measurable set then |A| denotes its d-dimensional volume. The symbol − A denotes |A| −1 A , provided that |A| > 0. If A has a well-defined area or length, then these are denoted, respectively, by area(A) and length(A).
For a measurable function f : Partial derivatives with respect to a variable x j , say, are denoted by ∂/∂x j . The Jacobi matrix of a differentiable function f : A → R k is denoted by ∂f . The symbol ∂ will also be redefined in some contexts, but used with essentially the same meaning as here.
When f is a deformation or a displacement, then we will also write ∇f = ∂f . For r ∈ R d , the uni-directional derivative is denoted by whenever this limit exists. The symbol D r denotes a finite difference operator, which will be defined in §2.1.2. Additional notation will be defined throughout. A list of symbols, with references to their definitions is given in §B.

Introduction to Atomistic/Continuum Model Coupling
In this section we introduce a general multi-body interaction model with periodic boundary conditions. This choice of boundary condition is crucial since the non-locality of the atomistic interactions makes an analysis of Dirichlet or Neumann boundaries challenging. The periodic boundary conditions can also be understood as "artificial boundary conditions" for infinite crystals.
Next, we describe the construction of energy-based a/c methods that couple the atomistic interaction potential with a P 1 -finite element discretization of the Cauchy-Born continuum model. We motivate the patch test ("ghost forces"), and why an interface correction is required to obtain accurate coupling schemes.
2.1.1. Periodic deformations. Let d ∈ {1, 2, 3} denote the space dimension. We will, in subsequent sections, restrict our analysis to d ∈ {1, 2}, however, the introduction to a/c coupling methods is independent of the dimension.
For some N ∈ N, and ε := 1/N , we define the periodic reference cell The space of 2Z d -periodic displacements of L # := εZ d is given by and the space of deformations with prescribed macroscopic strain A by Y is a linear space, while Y A is an affine subspace of Y .
For future reference we define notation that extends sets periodically: For any set A ⊂ R d we define A # := ξ∈2Z d (ξ + A). This notation is consistent with the definition of L # . If A is a family of sets, then We assume that the stored elastic energy of a deformation y ∈ Y is given in the form where is a multi-body interaction potential. Under these conditions, E a ∈ C 2 (Y ). The scaling of the lattice, of the finite difference operator, and of the energy were chosen to highlight the natural connection between molecular mechanics and continuum mechanics. For example, ε d resembles an integral (or, a Riemann sum), while D r resembles a directional derivative. It should be stressed, however, that ε is a fixed parameter of the problem, which is small but does not tend to zero.
The formulation (3) includes all commonly employed classical interatomic potentials (see, e.g., [19,47]): pair potentials such as the Lennard-Jones or Morse potential, bond-angle potentials, embedded atom potentials, bond-order potentials, or any combination of the former, provided that they have a finite interaction range. The generality of the interaction potential also includes effective potentials obtained for in-plane or anti-plane deformations of 3D crystals.
No major difficulties should be expected in generalizing the analysis to infinite interaction ranges, provided the interaction strength decays sufficiently fast. A generalization of the analysis to genuine long-range interactions such as Coulomb interactions is not obvious.

2.1.3.
Assumptions on the interaction potential. For g = (g r ) r∈R ∈ (R d ) R we denote the first and second partial derivatives of V at g, respectively, by Throughout this work we assume the following global bound: there exist constants M a r,s ≥ 0, r, s ∈ R, such that sup where · denotes the 2 -operator norm. This assumption contradicts realistic interaction models and is made to simplify the notation; see §2.2.1 for further discussion of this issue.
In the subsequent analysis we will, in fact, never make direct use of the second partial derivatives ∂ r,s V directly, but only use the resulting Lipschitz property, The proof is straightforward. For future reference, we also define the constant 2.1.4. The variational problem. Let P a ∈ C 2 (Y ; R) be the potential of external forces modelling, for example, a substrate or an indenter. As explained in §2.2.3, one may also use such a potential to model simple point defects such as vacancies or impurities. Given a potential of external forces P a and a macroscopic strain A, we consider the problem of finding local minimizers of the total energy E tot a := E a + P a in Y A , in short, where argmin denotes the set of local minimizers. If y a solves (7), then y a is a critical point of E tot a , that is, where, for a functional E ∈ C 1 (Y ), we define the first variation of E at y, as If E ∈ C 2 (Y ) then the second variation is defined analogously as The same notation will be used for functionals defined on different spaces.
Fix A ∈ R d×d and u ∈ U , then where we have used the fact that u is periodic. Above, and throughout, we use the notation AR = (Ar) r∈R = D R y A (x).

2.2.
Remarks on the atomistic model.

2.2.1.
Invertibility of deformations. In §2.1.2 and §2.1.3 we have assumed that E a is twice differentiable at all deformations y ∈ Y , and that the second partial derivatives of the interaction potential are globally bounded. However, realistic interaction potentials V take the value +∞ if two atoms occupy the same position in space, and hence can only be differentiable at deformations that are one-to-one (i.e., "true" deformations). With only minor additional technicalities such potentials can be admitted in the analysis. The global bound (4) would then be replaced by a local bound and certain explicit bounds on D R y; see, e.g., [40,41].

2.2.2.
Reference cutoff. Another aspect of the atomistic energy (3), which makes it inappropriate for realistic applications is that the interaction potential V has a cut-off radius in the reference configuration. In atomistic models, atoms are unconstrained in their position and hence two atoms that are far apart in the reference configuration may be arbitrarily close, and hence interact, in the deformed configuration.
The reference cutoff in (3) is assumed only for the sake of brevity of the notation. One can, similarly as discussed in §2.2.1, take a more general form of the interaction potential that does not suffer from this drawback and make suitable assumptions on deformations under consideration that control the interaction neighbourhood.

2.2.3.
Modelling crystal defects. Some simple crystal defects can be modelled via the potential of the external forces, P a . The simplest example is an impurity, where a single atom is replaced with an atom from a different species. For (3) this means that the interaction potential is changed from V (D R y(x)) to V mod (x; D R y(x)) in a neighbourhood of the impurity. Alternatively, one may keep the original form of E a and define Similarly, a vacancy can be modelled by simply removing all interactions with a given atom. This would yield a difficulty with the "unused" degrees of freedom for the position of the vacancy atom, which could simply be removed from the system [41]. An interstitial (an additional atom) can also be modelled fairly easily, but one would need to augment the variable y by additional degrees of freedom for the position of the interstitial atom.
Dislocations, which possibly represent the most important class of crystal defects are, in general, more difficult to describe. In the atomistic minimization problem (7) they simply represent a special class of local minimizers, however, in the coupled atomistic/continuum models we discuss below most classes of dislocations less straightforward to embed; however, see [54,34] for simple examples.

2.3.
Construction of a/c coupling methods. The atomistic model problem (7) is a finitedimensional optimisation problem and is therefore, in principle, solvable using standard optimisation algorithms. However, typical applications where atomistic models are employed require of the order 10 9 to 10 12 atoms or more [35,37]. It is therefore desirable to construct computationally efficient coarse grained models.
2.3.1. Galerkin projection. To motivate the idea of atomistic-to-continuum coupling we consider a crystal with a localized defect. Figure 1(a) shows a deformed 2D crystal with an impurity that repels its neighbouring atoms, causing a large local deformation. We observe that, except in a small neighbourhood of the defect, the atoms are arranged as a "smooth" deformation of the reference lattice L # . It is therefore possible to approximate the atomistic configurations from a low-dimensional subspace constructed, for example, using a P 1 -finite element method.
Let Ω := (−1, 1] d , and let T h be a regular [7] triangulation ofΩ, with vertices belonging to L # , that can be extended periodically to a regular triangulation T # h of R d . We make the convention that elements T ∈ T h are closed sets. For T ∈ T h we define h T := diam(T ), and We define the P 1 finite element space and we denote the spaces of piecewise affine displacements and deformations, respectively, by (c) Visualization of the QCE method described in §2.3.2. The blue shaded region is the set Ω qce c , the red atoms inside the white region are the set L qce a (both after deformation).
(d) Visualization of an interface correction as described in §2.3.4. The blue shaded regions is the continuum region Ω c , the green shaded region is the interface region Ω i , and the white region is the atomistic region Ω a . The mesh in Ω i ∪ Ω a is chosen so that it coincides with T ε (cf. §3.1).
For future reference, let I h : Y → Y h denote the nodal interpolation operator. We note Let F # h denote the set of closed edges of the extended triangulation T # h , and let F h denote the set of all edges f ∈ F # h such that area(f ∩ Ω) = 0. Finally, we introduce the spaces of piecewise constant functions P 0 (T h ), P 0 (T # h ), and P # 0 (T h ), defined in a similar manner. The Galerkin approximation of (7) is the coarse-grained minimization problem y a,h ∈ argmin where we recall that E tot a = E a + P a .
For example, if we choose the triangulation T h as in Figure 1(b), then we obtain full atomistic resolution in the neighbourhood of the defect, while considerably reducing the overall number of degrees of freedom by coarsening the mesh away from the defect. In this way it is possible to obtain highly accurate approximations to nontrivial atomistic configurations. Under suitable technical assumptions it is not too difficult to employ the classical techniques of finite element error analysis in this context. Such analyses, including a posteriori grid generation, are given in [28,29,42].

Remark 1 (Regularity of Atomistic Solutions).
In order for the Galerkin projection method, or the subsequent a/c coupling methods, to be accurate we require "regularity" of atomistic solutions. Such a regularity theory does not exist at present, however, most numerical experiments performed on atomistic models for simple lattices indicate "smoothness" of atomistic deformations away from defects. The situation would be different for so-called multi-lattices, which require a homogenisation step and represent a far greater challenge.

2.3.2.
Continuum region & Cauchy-Born approximation. The Galerkin projection (10) reduces the number of degrees of freedom considerably, however, the complexity of computing E a | Y h is not reduced in the same manner. Due to the non-locality of the atomistic interaction E a | Y h cannot be evaluated as easily as in the case of finite element methods for continuum mechanics. Several attempts have been made to use quadrature ideas to approximate E a and render (10) computationally efficient [18,21,24], however, it was shown in [31] that these approximations yield unacceptable consistency errors.
An alternative approach, proposed in [37], is to keep the full atomistic description for atomistically fine elements, while employing the Cauchy-Born approximation for coarser elements as well as an interface region. Following the terminology of [12] we call the resulting method the QCE method (the original energy-based quasicontinuum method).
To formulate this method we choose a set L qce a ⊂ L of atoms that we wish to treat atomistically (the red atoms in Figure 1
With this notation, the QCE energy functional is defined, for y h ∈ Y h , as where W : R d×d → R is the Cauchy-Born stored energy function, Note that W (F) is the energy of a single atom in the Bravais lattice FL # . One may readily check that the complexity of evaluating E qce (or its derivatives) is of the order O(#T h ), that is, of the same order as the number of degrees of freedom. Moreover, it is easy to see that E qce (y A ) = E a (y A ) for all A ∈ R d×d . For future reference, we give a formal definition of this property.

Definition 1 (Global Energy Consistency). We say that an energy functional
More generally one can show that E qce (y h ) is a good approximation to E a (y h ) if ∇y h varies only moderately. Despite these facts, it turns out, as we discuss in §2.3.3 and §4.2, that minimizers of E qce are poor approximations to minimizers of E a .

Remark 2 (The Cauchy-Born Model).
The Cauchy-Born model is a standard continuum model, for large deformations of single crystals. In the absence of defects, solutions of a pure Cauchy-Born model (no atomistic region) can provide excellent approximations to solutions of the atomistic model (7). For example, in [17] it is shown that, for smooth and small dead load external forces, and under realistic stability assumptions on E a , there exist solutions y a of (7) and y c of the Cauchy-Born model, such that where C depends on higher partial derivatives of V and on the regularity of y c . Hence, in the regime of "smooth elastic" deformations, the Cauchy-Born model can be considered an excellent approximation to the atomistic model (3).

The patch test.
The patch test is often employed in the theory of finite element methods [53,6,5] as a simple test for consistency. The test also plays an important role in the design of a/c methods.

Definition 2 (Patch Test Consistency). We say that an energy functional
The terminology "patch test consistency" is motivated by Proposition 1, where we have shown that the exact energy E a does satisfy the patch test (14).
However, the QCE energy E qce is not patch test consistent [50]. This result will be reviewed for a one-dimensional model problem in §4.2, where it will also be shown how failure of the patch test affects the consistency error.
Remark 3. In most of the a/c coupling literature the patch test is stated as the condition that ∂E (y h ) ∂y h (p) y h =y F = 0 for all finite element nodes p.
It is straightforward to see that this condition is equivalent to the variational formulation given in (14).
In the present paper we will focus only on energy-based a/c methods that are patch test consistent, i.e., that remove the "ghost forces" altogether. None of the methods presently available in the literature have resolved this problem in its full generality, however, [15,49,52,22] present several interesting approaches and partial solutions. In the following, we present a generalization of the geometrically consistent coupling method [15], which is the most general approach but leaves some questions concerning its practical construction open.
We assume that all atoms belonging to L \int(Ω c ) are vertices of T h , and we define two sets L gc a , L gc i such that L gc a ∪ L gc i = L \ int(Ω c ), and L gc a ∩ L gc i = ∅. For each x ∈ L gc i , we define a modified interaction potential V (x; ·) ∈ C 2 ((R d ) R ), and we define the a/c energy functional as For the modified potential V one takes a general ansatz with several free parameters, which are then fitted to remove or minimize the ghost force. For example, following the ideas of the quasinonlocal coupling method [52] and the geometrically consistent coupling method [15] one may define The constants C x,r,s can then be determined analytically as in [15,44], or, as proposed in [45], numerically in a preprocessing step. The 2D numerical experiments performed in [45] suggest that it is always possible to determine parameters C x,r,s such that E ac becomes patch test consistent, however, a proof of this fact is still missing. The purpose of the present work is to investigate the question whether patch test consistency is in fact a sufficient condition for first-order consistency of an a/c coupling method. If this would turn out to be false in general, then it would be necessary to develop new approaches for constructing accurate a/c methods.

2.3.5.
General assumptions on the interface correction. For the subsequent analysis we assume an even more general form of the a/c functional than (15).
and we define the continuum, interface, and atomistic regions Ω c := ∪T c h , Ω i := ∪T i h , and Ω a := ∪T a h . (Note that Ω c , Ω i , Ω a are closed sets.) Next, we define the set of all nodes L a ⊂ L that interact with the atomistic region: where the ordered pair (x, x ) ∈ L # × L # is called a bond; here, and throughout, the symbol (x, x ) is also identified with the segment conv{x, x } (in particular, it is closed). To avoid interaction between L a and Ω c , we assume throughout that Next, we define the set of interface bonds that is, the interface functional E i is given as a function of the finite differences D r y(x) of bonds (x, x + εr) that are contained in the interface region. Note also the volumetric scaling ε d .
With this notation, we set The interface functional E i (y h ) specifies the different variants of a/c methods. It is easy to see that the functionals E qce and E gc , discussed above, fit this framework (in the case of E qce we have to drop the assumption (16)).
If we define the total a/c energy as E tot a := E ac + P ac , where P ac is a suitable a/c approximation to P a , the a/c approximation to (7) is If y ac solves (19), then it is a critical point of E tot ac = E ac + P ac : 2.3.6. The locality and scaling conditions. We define notation for first and second partial derivatives of E i as follows: We extend the definition periodically: if b ∈ B i and ξ ∈ 2Z d then ∂ ξ+b E i := ∂ b E i , and we make a similar definition for the second partial derivatives. In our analysis we will require two crucial properties on E i , which we call the locality and scaling conditions: The locality condition implies that the same bonds interact through E i as in the atomistic model. This condition can be weakened, by requiring that only bonds within an O(ε) distance interact, however, such a more general condition would add additional notational complexity.
In the scaling condition we assume that there exist constants M i r,s ≥ 0, r, s ∈ R, such that This condition effectively yields an O(1) Lipschitz bound for δE i in the function spaces we will use. The scaling aspect enters through an implicit assumption on the magnitude of the constants M i r,s , namely, we will assume throughout that the constant is of the same order of magnitude as the constant M a defined in (6).
Remark 4. The factor 1 2 for bonds on the boundary of the interface region is not strictly necessary, since it can be removed by simply replacing M i r,s with 2M i r,s , however, if stated as above it makes the statements of the results in §6 slightly sharper, and moreover simplifies the argument in (68).
The necessity of this factor is related to the fact that we allow E i to depend on bonds that lie on the boundary of Ω # i ; this is made clear in Proposition 13 where we construct a stress function for E ac . Note that if we did not allow E i to depend on these boundary bonds, then it would in fact be impossible to construct patch test consistent a/c methods for non-flat a/c interfaces.

A Framework for the a priori Error Analysis of a/c Methods
When analyzing the error of a numerical method, one should first of all determine the main quantities of interest. For a/c methods, one is usually interested in energy differences between homogeneous lattices and lattices with defects, or critical loads at which defects form or move (i.e., bifurcation points). Since the present paper is mostly theoretical, we will simply focus on the error in the deformation gradient. We note, however, that many aspects of this analysis are crucial ingredients for the analysis of energy differences (see, e.g., [41]) and would usually also enter an analysis of bifurcation points.
We assume from now on that d ∈ {1, 2}. To execute the abstract framework of this section also in 3D, several technical tools as well as the central consistency result need to be developed first.
3.1. Discrete and continuous functions. In the following analysis it will be important to extend the a/c functional E ac to all functions y ∈ Y . To that end, we first define piecewise affine interpolants with respect to an atomistic mesh T ε .
We take a subdivision of the scaled unit cube ε(0, 1) d into d-simplices (in 1D the interval ε(0, 1); in 2D two symmetric triangles; compare with the triangulation of the atomistic region in Figure 1(d)), which we extend periodically to a triangulation T # ε of R 2 with vertex set L # . The restriction of T # ε to Ω is denoted by T ε . Each discrete function v : L # → R k will from now on be identified in a canonical way with its continuous piecewise affine interpolant v ∈ P 1 (T # ε ) k . For future reference we denote the sets of edges of T ε and T # ε , corresponding to the definitions of F h and F # h in §2.3.1, by F ε and F # ε .
3.1.1. Ambiguity of continuous interpolants. If y h ∈ Y h then y h can also be interpreted as a member of Y and therefore has two, possibly different, continuous interpolants. To distinguish them, we make the convention that the symbol y h always denotes the interpolant in P 1 (T h ) d , while a symbol y always denotes the interpolant in P 1 (T ε ) d . If we wish to evaluate the P 1 (T ε ) dinterpolant of a function y h ∈ Y h then we will write I ε y h . To compare a P 1 (T ε ) d -interpolant with a P 1 (T h ) d -interpolant, we use the following lemma. In 1D the result is easy to establish; in 2D it depends on a technical tool that we introduce in §5.1. The proof is given in the appendix.
where we recall that we have defined ∇v L p = |∇v| p L p .

3.1.2.
Extension of the a/c energy. Before we extend the a/c energy E ac we make one last technical assumption, which considerably simplifies the subsequent analysis. We shall assume from now on, that which is not uncommon, then this is no restriction.
With these conventions the a/c energy E ac defined in (18) can be defined canonically for functionsỹ = y + u h ∈ Y + U h by the same formula: We also assume that P ac has a suitable extension to Y . It should be stressed that for general 3.2. Measuring smoothness; An interpolation error estimate. The three main ingredients in the a priori error analysis of Galerkin-like approximations are (i) consistency, (ii) stability, and (iii) an interpolation error estimate. We begin by establishing the latter. To that end, we first need to find a convenient measure of smoothness for discrete functions y ∈ Y .
3.2.1. Measuring smoothness in terms of local oscillation. There are several possibilities to measure the "smoothness" of a discrete function. The most obvious is possibly the use of higher order finite differences, e.g., D e i D e j y(x This last observation motivates the idea to measure smoothness of y by the local oscillation of ∇y. We define the oscillation operator, for measurable sets ω ⊂ R d , and for y ∈ Y , as osc(∇y; ω) := ess sup The sets ω that arise naturally in our analysis will always have O(ε) diameter, which is the reason for the ε −1 -scaling in the definition of osc.
Note, in particular, that if y were twice differentiable, and if diam(ω) ≤ Cε, then we would obtain which further illustrates that the oscillation operator is a reasonable replacement for ∇ 2 y to measure the local smoothness of a piecewise affine function.

3.2.2.
Interpolation error estimate. The smoothness measure we defined in §3.2.1 yields a simple proof of an interpolation error estimate; see Appendix A.
then there exists a constant C I , which depends only on the shape regularity of T h , such that, for all y ∈ Y , p ∈ [1, ∞), 3.3. The stability assumption. The stability of a/c methods relies, firstly, on the stability of atomistic models as well as their Cauchy-Born approximations. It requires a thorough understanding of the physics of a model and in particular more specific information about the interaction potential. Since the focus of the present work is the consistency of a/c methods, we will formulate stability as an assumption.
The simplest notion of stability one may use, which is also closely connected to local minimality, is coercivity of the second variation: where c 0 > 0, and y h ∈ Y h is a suitable deformation in a neighbourhood of the atomistic solution y a , e.g., y h = I h y a . The choice of norm is motivated by the fact that the Cauchy-Born model, and hence the atomistic model, are closely related to second order elliptic differential equations.
Examples of sharp stability estimates for a/c methods in 1D can be found in [13,40,27,26,55]. For pair interactions in 2D the stability of Shapeev's method [49] is established in [41].
More generally, for some p ∈ [1, ∞], p = p/(p − 1), we may assume an inf-sup condition of the form for some constant c 0 > 0. The condition (27) is usually difficult to prove, especially for p = 2, and may indeed be false in general. We will only use it to demonstrate how such a stability result motivates consistency estimates in different negative norms. Examples of 1D inf-sup stability estimates for a/c methods can be found in [10,14,32,42].

3.4.
Outline of an a priori error analysis. The following outline of an a priori error analysis depends on a stability assumption that we will not prove. Moreover, since it primarily serves to motivate the consistency problem, and since a rigorous derivation would be more involved without yielding much additional insight, some steps will be kept vague. Most of these steps are easily made rigorous; the main assumption we make below, which is in fact very difficult to justify rigorously, is that I h y a and y ac are "sufficiently close". See [32,40,42,55] for similar analyses in 1D where all steps are rigorously justified, and [43,41] for a similar semi-rigorous framework, where a proof of this step is replaced by an assumption. Let y a satisfy (8), y ac satisfy (20), and suppose that the stability assumption (27) holds with y h = I h y a . Let e h := y ac − I h y a . Moreover, suppose that ∇I h y a − ∇y ac L ∞ is sufficiently small so that the following approximation can be made precise: Taking the supremum over all v h ∈ U h , and invoking the inf-sup condition (27) and the criticality condition (20), we obtain where we define We split the consistency error δ(E ac + P ac )(I h y a ) W −1,p h into three separate contributions: where we have used (8), and the extension of E ac and P ac for all deformations y ∈ Y constructed in §3.1.
The coarsening error, E coarse h , can be bounded by Lipschitz estimates for δ(E ac + P ac ) and an interpolation error estimate. Using our assumption that T a h ∪ T i h ⊂ T ε it is not difficult to derive Lipschitz estimates of the form where M a is a Lipschitz constant for ∂W (cf. (12)), and M Pac is a Lipschitz constant for δP ac . Combining (28), and (29), the inequality ∇y a − ∇y ac L p ≤ ∇y a − ∇I h y a L p + ∇e h L p , and the interpolation error estimate of Lemma 3, we arrive at the following basic error estimate where The consistency error for the external forces, E ext h , depends on the form of P a and P ac and cannot be discussed at this level of abstraction. The modelling error, E model h , is the focus of the remainder of the present paper.

Remark 5 (Choice of Splitting).
Suppose, for simpliciy, that P ac = P a = 0. In a typical finite element error analysis of continuum mechanics problems one would usually choose a different splitting of the consistency error: This splitting was used in the analysis in [41] and led to a suboptimal estimate of the modelling error, since it still contains some coarsening error.
3.5. The consistency problem. The main step that remains in obtaining an a priori error estimate from (30) is the estimation of the modelling error Most of the numerical analysis literature on a/c methods estimates this modelling error only for the case when T h = T ε . In 1D it is easy to see that this is sufficient, since I ε v h = v h in that case; see also [43]. The following lemma provides the main technical step to explain why it is also sufficient in 2D to consider the case T h = T ε . Its proof uses arguments similar to those in the a posteriori error analysis of continuum finite element methods and is given in Appendix A. for all u ∈ U , u h ∈ U h , where σ ∈ P # 0 (T ε ) d×d ; then there exists a universal constant C M such that, for all p ∈ [1, ∞), The estimate (31), together with Lemma 2, implies the following theorem, where we use the notation Theorem 5. Suppose that T a h ∪ T i h ⊂ T ε and that y ∈ Y ; then, for all p ∈ [1, ∞), with corresponding statement for p = ∞.
Proof. Since E a (I ε y h ) uses only point values of I ε y h , which are the same as for y h , we have Using this fact, we can estimate Due to the assumption that T a h ∪ T c h ⊂ T ε , the first group can be estimated using Lemma 4, with σ = ∂W (∇y), which yields the first term in (33).
Using Lemma 2, the second group can be estimated by Taking the supremum over all u h ∈ U h with ∇u h L p = 1 yields the stated result.
Applying Theorem 5 to the modelling error E model h , defined in (28), we obtain that where Even though E model ε is essentially an upper bound for E model h , it is usually easier to estimate. The consistency problem is to prove a sharp upper bound on E model ε .
In §4 we will discuss two simple 1D examples to determine what can be expected in more general situations. In Theorem 10 we will prove that for an a/c method that is patch test consistent, and satisfies various other technical conditions, one obtains where C is a constant that is independent of y a , but does depend on the interface width, and ω T ⊂ Ω c ∪ Ω i is the interaction neighbourhood defined in (43). Combined with (34) and (30), and using the fact that ε ≤ h T , and that ω T ⊃ ω c T for T ∈ T c ε , this bound yields where c 2 is a constant that is independent of y a . This estimate closely resembles a typical first-order a priori error estimate for a continuum mechanics finite element approximation; see also the interpretation given in §1.2.
It should be stressed again that (35) is not a rigorous error estimate, but depends on various assumptions made in the forgoing discussion, most prominently, the stability assumption (27), and the assumption that ∇I h y a − ∇y ac L ∞ is "sufficiently small". Remark 6. The locality of the patches ω T is crucial. If diam(ω T ) is not of the order O(ε), then it is possible that osc(∇y a ; ω T ) 1 even if y a is globally smooth; see also §6.4.4.

Examples in 1D
In the present section we review the consistency analyses of specific a/c methods to point out the main features and to motivate what may be expected in the general case. Throughout this section we assume that d = 1, R = {±1, ±2}, and that V is given by , where φ 1 , φ 2 ∈ C 2,1 (R) are, respectively, the first and second neighour interaction potentials, which are assumed to be symmetric about the origin. We assume that their derivatives φ i and φ i have global Lipschitz constants m i and m i .
For the 1D analysis it is convenient to write x n = nε, v n = v(x n ), and to write all interactions in terms of the backward difference operator v n = v n − v n−1 ε .
With this notation the atomistic energy can now be rewritten in the form where we note that y n + y n+1 = ε −1 (y n+1 − y n−1 ) describes a second neighbour bond. For future reference we also define the second and third finite differences It is also worth pointing out that v n = ∇v(s) for all s ∈ (x n−1 , x n ).

4.1.
Consistency of the QNL method. We begin with a modelling error analysis for the quasinonlocal coupling method (QNL method) of Shimokawa et al [52]. The geometrically consistent coupling scheme [15] and the method proposed by Shapeev [49] reduce to the same method for 1D second neighbour interactions.
The following presentation follows largely [40], where the QNL method is defined as follows: Let N a = {−K, . . . , K} for some K ≥ 1 and N c = {−N + 1, . . . , N } \ N a ; then, for y ∈ Y , the QNL energy is defined by We observe that we have not modified the first neighbour interactions, but have "split" the non-local second neighbour interactions into local first neighbour interactions in the continuum region.
It is straightforward to rewrite E qnl in the form specified in (18), with and a suitably defined interface functional E i ; however, the form (37) is more convenient for the analysis.
The following modelling error estimate was first established in [40, Thm. 3.1]. Dobson and Luskin [12] treated a quadratic interaction case, using entirely different analytical techniques that gave an even more detailed analysis of the error; Ming and Yang [36] used related methods as [40, Thm. 3.1], but gave a qualitatively less precise estimate of consistency error. An extension of the result to linear finite range pair interactions is given in [27].
Note also that it is shown in [12] that the consistency error of the QNL method in p ε -norms is of the order O(1), that is, the usage of negative norms cannot be avoided.
We will discuss the estimate in detail in §4. 3. The proof of the following result serves as a first guidance on how one may approach proofs of consistency of a/c methods in more general situations.
Proposition 6 (Consistency of the QNL Method). Let y ∈ Y ; then Proof. Throughout the proof we will make use of the fact that the boundary conditions are periodic without comment, treating the boundary as if it belonged to the "interior" of the continuum region.
After inserting the two bounds into (38), and applying several weighted Hölder inequalities one obtains the stated estimate.

4.2.
Inconsistency of the QCE method. As in the previous section, let N a = {−K, . . . , K}, N c = {−N + 1, . . . , N } \ N a , and L a = εN a ; then, under the assumptions and notation set out at the beginning of §4, the QCE energy functional defined in (11) reads The following result is a variant of [43,Thm. 3.2]. Previous analyses of the QCE method [11,12,36] computed the consistency error contribution due to the "ghost forces" explicitly rather than estimating their W −1,p ε -residual contribution.
To obtain the opposite estimate for y = y A , a brief computation gives then ∇u L p (−1,1) = u p ε = 1, and we obtain that

4.3.
Discussion. This discussion of the 1D consistency error estimates largely follows the discussions in [43,40]. We have estimated the modelling errors for two prototypical a/c methods. We see that the leading order terms in the upper bounds are O(ε) and O(ε 1/p ) for the QNL and QCE methods, respectively. However, a much finer distinction should be made.
First, we note that both methods reduce to the Cauchy-Born approximation in the continuum region, and the corresponding contributions are all of second order (see also Remark 2; note also that this requires point symmetry of V , which we have not assumed in general in this paper).
Second, we see that the QCE method (and only the QCE method) has a zeroth-order term Gε 1/p in the interface region. This term occurs since the QCE method is not patch test consistent, that is, homogeneous deformations are not equilibria of the QCE model: The origin and effect of these "ghost forces" are discussed in more detail in [50,35,11,12,36].
We should call this term zeroth order for several reasons: Firstly, it is clearly of zeroth order if p = ∞, in which case the consistency error is related to the error in the W 1,∞ -norm. Secondly, the parameter ε is a constant of the problem and does not tend to zero. As a matter of fact, the accuracy of an a/c method should be related the smoothness of the solution (as opposed to the atomistic scale), and the term Gε 1/p is independent of the magnitude of y in the interface region. The scaling ε 1/p relates only to the width of the interface region.
Finally, it is worth remarking on the first-order consistency term in the interface region for the QNL method. The reason this term is of first order as opposed to second order is the loss of symmetry that is introduced by changing the interaction law at the interface between the atomistic and continuum regions. A recent result of Dobson [8] shows that no a/c method coupling an interatomic potential to the Cauchy-Born continuum model can achieve better than first-order accuracy in the interface region.
Note also, that the second finite differences y n can in fact be written in terms of the oscillation operator: In our analysis in §6, we will ignore the possibility of proving a modelling error estimate that is of second order in the continuum region, but we will be satisfied with an estimate that is globally of first order.

5.1.
The bond density lemma. The bond density lemma is a tool that allows a transition between integrals over bonds, and volume integrals. It was first derived in [49] for the construction of patch test consistent a/c methods for pair potential interactions. Related asymptotic results were used previously in Γ-convergence analyses of atomistic models [3]; the achievement of Shapeev [49] was to obtain a formula that is exact for any triangle.
Before we formulate the result we introduce some notation. Let T ⊂ R 2 be a triangle with vertices belonging to L # . We define the characteristic function χ T : R 2 → R by where B t (x) denotes the closed ball with centre x and radius t. Note that χ T = 1 in int(T ), and χ T = 1/2 on the edges of T . The value of χ T on the corners is not of importance.
Let x, x ∈ R 2 and let ϕ be a function that is measurable on the line segment (x, x ); then we define the line integral, or bond integral, Lemma 8 (Bond Density Lemma ( [49], Lemma 2)). Let T ⊂ R 2 be a non-degenerate triangle with vertices belonging to L # = εZ 2 , and let r ∈ Z 2 ; then As a first application of the bond density lemma we present a proof of Lemma 2 in the appendix.
Remark 7. In the above form, the bond density lemma is false in 3D, which is one of the reasons why the present work is restricted to 2D. Moreover, the condition that the vertices of T belong to lattice sites is also necessary. This is related to the assumption that the vertices of T h belong to L # , however, this is not crucial and could be removed with some additional work.

5.2.
Discrete divergence-free tensor fields in 2D. Our second auxiliary result concerns representations of discrete divergence-free P 0 -tensor fields. Following the construction given by Polthier and Preuß [48], we will give a proof for the periodic setting. This proof also serves to motivate a crucial argument in §6.3.3. Since we will only use the atomistic finite element mesh T ε from now on, we will formulate everything in terms of this mesh. However, all results hold for general periodic triangulations.

5.2.1.
The Crouzeix-Raviart finite element space. The representation of discrete divergencefree tensor fields requires the use of the non-conforming Crouzeix-Raviart finite element space. Recall from §3.1, 2.3.1 the definition of the sets of edges F ε and F # ε . The Crouzeix-Raviart finite element space over T # ε is defined as The degrees of freedom for functions w ∈ N 1 (T # ε ) are the values at edge midpoints, w(q f ), f ∈ F # ε , and the corresponding nodal basis functions are denoted by ζ f . The Crouzeix-Raviart finite element space of periodic functions is defined as The periodic nodal basis functions are defined, for f ∈ F ε , by ζ # f = ξ∈2Z 2 ζ ξ+f .

Path integrals.
For two edges f, f ∈ F # ε , let Γ f,f denote the set of all piecewise affine paths from q f to q f , crossing element edges only in edge midpoints; see Figure 2(a) for an example.
For any piecewise constant vector field σ ∈ P 0 (T # ε ) 2 and for any path γ ∈ Γ f,f , γ = {x(t) | 0 ≤ t ≤ T }, we denote the standard path integral by γ σ · dx =  For piecewise constant tensor fields σ ∈ P 0 (T # ε ) k×2 we define the path integral as Since functions w ∈ N 1 (T # ε ) k have piecewise constant gradients ∇w, and since they are continuous in edge midpoints, it is easy to see that

5.2.3.
Discrete divergence-free tensor fields. The following lemma characterizes discrete divergencefree tensor fields.
Lemma 9. A tensor field σ ∈ P # 0 (T ε ) k×2 satisfies Ω σ : ∇u dx = 0 ∀u ∈ P # 1 (T ε ) k if and only if there exist a constant σ 0 ∈ R k×2 and a function w ∈ N # 1 (T ε ) k such that Sketch of the proof. The reverse direction, that any tensor field of the form σ = σ 0 + ∇wJ has zero discrete divergence, can be checked using a straightforward calculation, using (41).
Step 0. Outline: To simplify the notation, we define α = σJ T and, without loss of generality, we assume that k = 1. Initially, we treat α as a piecewise constant tensor field on all of R 2 , ignoring periodicity. We construct w ∈ N 1 (T # ε ) by explicitly specifying w(q f ) for all edges f ∈ F # ε . We will then show in the last step of the proof that w can be written as the sum of an affine function and a periodic function.
Step 1. Construction of w: Fix a starting edgef ∈ F ε and define w(qf ) = 0. For any edge f ∈ F # ε let γ ∈ Γf ,f and define w(q f ) via the path integral We need to show that this definition is independent of the path.
Let q be a vertex of the triangulation T # ε and let ϕ q be the corresponding nodal basis function with support ω q ; then a fairly straightforward calculation (see, e.g., [48] for the details) shows that where ν ± f are the two unit normals to f , and γ q is the piecewise affine path through edge midpoints circling q; cf. Figure 2(b). Note that the rotation J in the definition of α comes from the fact that tangent vectors are rotated normal vectors.
Since all closed piecewise affine paths can be written as a sum over paths γ q , this implies that γ α · dx = 0 for all closed piecewise affine paths γ, and in particular that the definition of w(q f ) is independent of the choice of path, that is, w is well-defined. Step which immediately implies that ∇w(T ) = α(T ).

A General Consistency Result in 2D
We are now finally in a position to make precise the statement that patch test consistent a/c methods are first-order consistent. Motivated by the example of the QNL consistency result (ignoring, as discussed in §4.3, the second-order consistency of the Cauchy-Born approximation), we would like to prove a result of the form δE ac (y) − δE a (y) W −1,p ε ≤ Cε ∇ 2 y L p (Ωc∪Ω i ) .
As explained in §3.2.1 we will use the oscillation operator (24) to replace the undefined second derivative.
For each element T ∈ T # ε , let ω a T be the interaction neighbourhood of T in the atomistic model, and ω T its union with the set ω c T , restricted to the continuum and interface regions, Note that assumption (16) implies that ω a T ∪ ω c T ⊂ ω T for all T ∈ T c ε .

Theorem 10 (First-order Consistency of Patch Test Consistent a/c Methods).
Suppose that E ac is patch test consistent ( §2.3.3) and globally energy consistent ( §2.3.2), that the locality condition (21) and the scaling condition (22) hold, that int(Ω a ) is connected, that T a h ∪ T i h ⊂ T ε , and that (16) holds. Then, for any y ∈ Y , we have where the oscillation measure osc is defined in (24), the interaction neighbourhood ω T is defined in (43), and the prefactors M T are defined as follows: In (45), the constants M a and M i are defined in (6) and (23), and the "interface width" width(Ω i ) is given by Outline of the proof. We will construct "stress functions" Σ a (y), Σ ac (y) ∈ P 0 (T ε ) 2×2 such that δE a (y), u = Ω Σ a (y) : ∇u dx, and δE ac (y); u = Ω Σ ac (y) : ∇u dx ∀u ∈ U .
If we could prove an estimate of the form Σ a (y; T ) − Σ ac (y; T ) ε osc(∇y; ω T ), then the result would follow immediately. It turns out that this is not possible. Instead, we will use the fact that E ac is globally energy consistent and patch test consistent to construct a correction (cf. Corollary 15, Lemma 17, and §6.2.3) Σ ac (y) = Σ ac (y) − ∇ψ(y)J, which still represents the first variation δE ac (cf. Lemma 9), and which has the property that In addition, we show that Σ ac (y; T ) = Σ a (y; T ) for all T ∈ T a ε . Lipschitz estimates for Σ a and Σ ac , and careful modifications of the argument in the interface region, yield the following result (cf. Lemma 19): and, in particular, which yields the stated first order consistency estimate.

Remark 8.
In 1D, a similar result can be proven using a similar framework but with significantly reduced technicalities. Note, in particular, that the 1D analogue of Lemma 9 is Ω σ · ∇u dx = 0 ∀ u ∈ U if and only if σ is constant.
Hence, the corrector function ψ(F, ·) defined in §6.2.2 is always identically equal to zero, which removes the interface width dependence from the modelling error (cf. §6.3.3). Hence, if d = 1, we obtain (44) again but with modified prefactors Moreover, since ψ ≡ 0 in 1D, and since symmetries are more easily exploited, it is not too difficult in 1D to prove second order consistency in the continuum region.
6.1. The atomistic stress function. A natural "weak" representation of δE a (y), y ∈ Y , is given by Using bond integrals we rewrite this in a form that will be useful for our subsequent analysis. We will then use the bond density lemma whenever we need to transition between bond integrals and volume integrals. This process yields a notion of stress for atomistic models, which is related to the virial stress (see [2] for a recent reference; this connection will be discussed in detail elsewhere). A variant of this result for pair interactions in 1D was developed in [32]. where the stress function Σ a (y) ∈ P # 0 (T ε ) 2×2 is defined as follows: Proof. For the sake of brevity we will write V x,r = ∂ r V (D R y(x)).
Recall that T # = ξ∈2Z 2 (ξ + T ), for T ∈ T ε . It is easy to see that {χ T # | T ∈ T ε } is a partition of unity for R 2 . Hence, using the identity we can rewrite (48) as Interchanging the order of summation and using the fact that ∇ r z = (∇z)r holds db-a.e. (note that, if the bond is aligned with an element edge, then ∇ r z is continuous across that edge) yields The term ∇z(T # )r = ∇z(T )r can be taken outside the bond integral, and hence, employing the identity a · (Gr) = (a ⊗ r) : G, yields Finally, we use the fact that D R y is 2Z 2 -periodic to deduce that The terminology "stress function" for Σ a is motivated by the fact that Σ a (y) takes precisely the same role as the first Piola-Kirchhoff stress tensor in the continuum theory of elasticity.
In the next lemma we prove two useful properties of the atomistic stress function Σ a . We show that Σ a = ∂W under locally homogeneous deformations and give a quantitative estimate for the discrepancy between Σ a and ∂W . Moreover, we have the estimate where M a is defined in (6) and ω a T is defined in (42).
Proof. Part 1: Since [∂ r V (FR) ⊗ r] is independent of x, we can apply the bond density lemma to the sum in curly brackets, to deduce that Recalling that W (F) = V (FR), it can be easily checked that the sum on the right-hand side of the second equality equals ∂W (F). Part 2: Let F = ∇y(T ), V x,r = ∂ r V (D R y(x)), and V F,r = ∂ r V (FR). From part 1 we obtain that We use the Lipschitz property 5 to estimate Since (x, x + εs) ⊂ ω a T , and recalling that F = ∇y(T ), we can further estimate max x ∈(x,x+εs) We combine (55) with (54) and insert the resulting estimate into (53) to arrive at An application of the bond density lemma, and referring to the definition of M a in (6), yields the stated result.
6.2. The a/c stress function. We wish to derive a similar representation of δE ac in terms of a stress function Σ ac , as we did in §6.1 for δE a . A straightforward calculation along the same lines as the proof of Proposition 11, recalling first the definition of the partial derivative ∂ b E i from §2.3.6, yields the following result.
where χ i T is a modified characteristic function: Proof. Employing again the notation V x,r = ∂ r V (D R y(x)), the functional δE ac (y) can be written as The first term on the right-hand side of (56) gives rise to the definition of Σ ac (y; T ) for T ∈ T c ε . (Note that (16) guarantees that the bonds in the second group in (56) do not contribute to Ω c .) After the same calculation as in the proof of Proposition 11, the second group on the righthand side of (56) gives the definition of Σ ac (y; T ) for T ∈ T a ε , as well as the first group in the definition of Σ ac (y; T ) for T ∈ T i ε . The crucial modification to the previous argument is that the modified characteristic functions χ i T , T ∈ T i ε form a partition of unity for Ω i (except at a finite number of points, which do not contribute to bond integrals). Therefore, performing again a similar calculation as in the proof of Proposition 11 to "distribute" the third group on the right-hand side of (56) between interface elements only, we obtain the second group in the definition of Σ ac (y; T ) for T ∈ T i ε . Note that if we hadn't made the modification to the characteristic function, then Σ ac (y; T ), T ∈ T a ε ∪ T c ε would contain contributions from E i .
Since any discrete divergence-free tensor field may be added to Σ ac and still yield a valid representation of δE ac , it is not surprising that, in general, Σ ac does not have the necessary property that Σ ac (y F ; T ) = ∂W (F) for all F ∈ R 2×2 . This can already be observed in the nearestneighbour, flat interface constructions in [44]. Hence, we need to construct an alternative stress function Σ ac representing δE ac that does have the desired properties. This construction will be undertaken in the remainder of this section, using the representation of discrete divergence-free vector fields as gradients of Crouzeix-Raviart functions discussed in §5.2.
6.2.1. Consequences of global energy consistency. Recall from (13) that a functional E ∈ C 1 (Y ) is called globally energy consistent if E (y F ) = E a (y F ) for all matrices F ∈ R 2×2 . The following lemma establishes a simple but crucial consequence of this property.
Proof. From the assumption of global energy consistency, and (12), we obtain that Since y F + ty G = y F+tG this implies that If we apply the foregoing lemma to an a/c functional E ac we obtain the following corollary.
Corollary 15. Suppose that E ac is globally energy consistent, then Proof. This result follows simply from Lemma 14 and the fact that, for all G ∈ R 2×2 , δE ac (y F ), y G = Ω Σ ac (y F ) : G dx = Ω Σ ac (y F ) dx : G.

Consequences of patch test consistency.
The examples given in [44] show that patch test consistency does not necessarily imply that Σ ac (y F ; T ) = ∂W (F). In the current paragraph we characterise the discrepancy between Σ ac (y F ; T ) and ∂W (F). First, we show that the test functions u h ∈ U h in the patch test (14) may be replaced by arbitrary displacements u ∈ U .
Lemma 16. Suppose that E ac is patch test consistent and that T a h ∪ T i h ⊂ T ε ; then we also have Proof. Fix u ∈ U ; then, using the assumption that T a h ∪ T i h ⊂ T ε , we have Since Σ ac (y F ) = ∂W (F) is constant in Ω c , and since I h u = u in Ω i ∪ Ω a , we can therefore deduce that The penultimate equality requires some justification, but follows quite easily from the particular form of E ac assumed in (18) and the assumption that Lemma 17. Suppose that E ac is patch test consistent and globally energy consistent and that where J is a rotation matrix defined in Lemma 9. Moreover, if int(Ω a ) is connected, then we may choose ψ(F) = 0 in Ω a .
Using global energy consistency of E ac and Corollary 15 we obtain that If Ω ∇ψ dx = 0, then Σ 0 = ∂W (F) and hence the result follows.
To prove this, we integrate by parts separately in each element: where, in the last equality, we used the fact that f (ψ + − ψ − ) ds = 0 for all edges f , since ψ is continuous in the edge midpoints. Finally, since Σ ac (y F ; T ) = ∂W (F) for all T ∈ T a ε (cf. Lemma 13 and Lemma 12), it follows that ∇ψ = 0 in Ω a . Hence, if int(Ω a ) is connected, then we can shift ψ by a constant so that ψ = 0 in Ω a . 6.2.3. The modified a/c stress function. We wish to construct a modified a/c stress function Σ ac that can be used to represent δE ac , and satisfies the crucial property that Σ ac (y F ; T ) = ∂W (F) for all F ∈ R 2×2 , T ∈ T i ε . To this end, we generalize the Crouzeix-Raviart function ψ(F) defined in Lemma 17 to arbitrary deformations y ∈ Y . Since we will use later on that the modified functionψ vanishes in Ω a , we require from now on that int(Ω a ) is connected so that we can choose ψ(F, ·) = 0 in Ω a for all F ∈ R 2×2 .
For each y ∈ Y and each face f ∈ F # ε , f = T 1 ∩ T 2 , we define the patch ω f = (T 1 ∪ T 2 ) \ Ω # a , and the deformation gradient averages Note that ω f was defined in such a way that is of no importance, and could have been replaced by any other value. With this notation, and recalling the definitions of the edge midpoints q f and the periodic nodal basis functions ζ # f from §5.2.1, we can definê Note, in particular, thatψ(y F ) = ψ(F) for all F ∈ R 2×2 . It is therefore natural to define the modified stress function In the following lemma we establish some elementary properties of Σ ac .
Lemma 18. Suppose that E ac is energy and patch test consistent, that int(Ω a ) is connected, and that T a h ∪ T i h ⊂ T ε ; then the modified a/c stress function Σ ac , defined in (59), has the following properties: Proof. To prove (60) let z = y B + u for some B ∈ R 2×2 and u ∈ U ; then Ω Σ ac (y) : ∇z dx = δE ac (y), z − Ω ∇ψJ : B + ∇u dx.
Since Ω (∇ψJ) : ∇u dx = 0 by Lemma 9, and since Ω ∇ψ dx = 0 (see the proof of Lemma 17), the representation (60) follows. Property (61) follows from Lemma 17 and the fact thatψ(y F ) = ψ(F): Property (62) follows from the fact that we constructedψ(y) to be zero in Ω a for all y ∈ Y , and from Proposition 13, which states that Σ ac (y; T ) = Σ a (y; T ) for all T ∈ T a ε .
6.3. The Lipschitz property. The final remaining ingredient for the proof of first-order consistency, is a Lipschitz property for Σ ac , similar to the Lipschitz property (52) of Σ a . In order to ensure that there are no modelling error contributions from the atomistic region it turns out to be most convenient to work directly with the stress difference R(y; T ) := Σ ac (y; T ) − Σ a (y; T ).
From (62) we immediately obtain that In the remainder of the section we will estimate R(y; T ) for T ∈ T c ε ∪ T i ε . To motivate the following result we note that, from Lemma 12 and from (61) we see that hence, a suitable Lipschitz estimate for R yields the following result.
Lemma 19. Suppose that all conditions of Lemma 18 hold and, in addition, that E i satisfies the locality and scaling conditions (21) and (22); then where M T is defined in (45).
The proof of this central lemma is split over the following paragraphs.
where, in the last inequality, we used (52) and the fact that ω a T ⊂ ω T for T ∈ T c ε . We still need to estimate ∇ψ(y; T ), which we postpone until §6.3.3.

6.3.2.
Estimates in the interface region. Let T ∈ T i ε , and let V x,r = ∂ r V (D R y(x)), then R(y; T ) = Σ ac (y; T ) − Σ a (y; T ) − ∇ψ(y; which, after combining the first and third group, becomes We will again postpone the estimation of ∇ψ(T ) to §6.3.3, and focus on the terms R (1) (y; T ) and R (2) (y; T ). Let F = ∇y(T ). Using the locality and scaling conditions (21) and (22), we can estimate In the transition from the first to the second line we have used the fact that on bonds that lie on the boundary of Ω # i the constant M i r,s is replaced by 1 2 M i r,s , which effectively replaces χ i T by χ T . Bounding |D s y(x) − Fs| by the local oscillation, and applying the bond density lemma (see Lemma 12 for a similar calculation), we deduce that Following closely the proof of (52), we obtain a similar estimate for R (2) : |r||s|M a r,s osc(∇y; ω T ) = εM a osc(∇y; ω T ).
Note that it is enough to measure the oscillation over ω T (which does not intersect with Ω a ), since (x, x + εs) ∩ Ω a = ∅ for all x ∈ L # \ L # a , s ∈ R. Combining (69) and (70), and using the fact that R(y F ; T ) = 0, we conclude that where width(Ω i ) is defined in (46).
• Specific a/c methods: Suppose that the GCC method described in §2.3.4 is used to construct a patch test consistent coupling scheme, with parameters C x,r,s . Suppose, moreover, that Ω a has two connected components, Ω 1 and Ω 2 , each of which have a portion of the boundary with the same orientation (say, normal e 1 ), as displayed in Figure 3. It is then reasonable to assume that the parameters C x,r,s have the same value in those parts of the interface surrounding Ω 1 and Ω 2 , which would imply that Moreover, since Σ ac (y F ) = ∂W (F) in the continuum region, we would obtain that where γ 2 denotes the curve γ 2 with reversed orientation. This shows that it is possible to choose ψ(F; ·) = 0 in both components of Ω a , and as a consequence, Theorem 10 would remain true. A related issue is the dependence of the modelling error estimate (44) on width(Ω i ), which comes solely from the Lipschitz estimate on F → ψ(F; ·); cf. §6.3.3. Hence, a better understanding of this function may also allow a finer analysis of this undesirable dependence. 6.4.3. The global energy consistency condition. Global energy consistency is a natural and convenient condition that yields the important intermediate result (Corollary 15) that Note also that (73) implies E ac (y F ) = E a (y F ) + c for all F ∈ R 2×2 , where c is a fixed constant that is independent of F; that is, (73) is practically equivalent to global energy consistency.
In some important situations patch test consistency already implies (73). The following result gives such a result for finite atomistic regions.
Proof. According to Lemmas 13 and 17 we have Let G ∈ R 2×2 and choose any u ∈ U such that ∇u = G in Ω i ∪ Ω a ; this is possible due to the assumption that Ω a ∪Ω i ⊂ int(Ω). Integrating by parts twice, letting ν denote the unit outward normal to Ω i ∪ Ω a , and noting that the portions of the surface integrals along ∂Ω cancel each other out, yields Since E ac is patch test consistent, the last term vanishes, and hence the result follows.
It turned out to be difficult to devise a counterexample, which clearly demonstrates that absence (73) can yield an inconsistent method. A more thorough investigation of this condition is still required.
From the definition it is obvious that J(y F ) = 0 for all F ∈ R 2×2 . Moreover, using summation by parts along the two lines L i −,+ ∪ L i +,+ and L i −,− ∪ L i +,− , it is easy to check that the patch test (14) holds. Finally, J satisfies the scaling condition, ∂ (x,x+εr) ∂ (x ,x +εr) J(y) = 2 I, if r = e 1 and x, x ∈ L i a,b , a, b ∈ {+, −}, 0, otherwise.
However, J clearly violates the locality condition.
We may think of J as an a/c functional for E a = 0, or, alternatively, as an additional contribution that can be added to any a/c functional whose interface satisfies L i −,+ ∪ L i −,− ∪ L i +,+ ∪ L i +,− ⊂ int(Ω # i ). Let y ∈ Y A be "smooth" but not affine in the upper half plane {x ∈ L | x 2 ≥ 0}, and let y = y A on the lower interface L i −,− ∪ L i +,− ; then, testing δJ (y) with the unique displacement u ∈ U such that u(x) = y(x) − Ax for x ∈ {(0, 1/2), (1, 1/2)}; This final estimate is scaled like a surface integral and is clearly a zeroth order term if y is smooth but y(0, 1/2) = y(1, 1/2). This shows that the locality condition (21) (or a variant thereof) is indeed necessary to obtain a first-order consistency estimate.
6.4.5. The scaling condition. It is fairly clear that the modelling error estimate can be arbitrarily large without the scaling condition (22). We nevertheless briefly discuss a simple example with a natural interpretation.
Using a similar argument as in the previous paragraph, we define a functional J ∈ C 2 (Y ), J = ε 2 J, J(y) = β where β > 0 is a constant, and where L i ± = x ∈ L x 2 = ±1/2 . It is easy to see that J is patch test consistent, that J(y F ) = 0 for all F ∈ R 2×2 , and that it satisfies the locality condition.
Let y ∈ Y such that y(x) = Ax for x 2 = −1/2, then testing δJ (y) with the unique displacement u ∈ U such that u(x) = y(x) − Ax, for x ∈ L i + , D e 2 u(x) = 0, for all x ∈ L , If y is "smooth" but not affine, then the term in square brackets is of the order O(1). By choosing β arbitrarily large, the modelling error can be made arbitrarily large as well. In particular, the choice β = 1/ε would give a seemingly natural surface scaling to the interface functional, and in this case we would obtain an O(1) modelling error.

Conclusion
A fairly complete consistency analysis of general patch test consistent a/c coupling methods in (one and) two space dimensions was developed in this paper. The main result is the first order modelling error estimate, Theorem 10. The main undesirable condition is the assumption that int(Ω a ) is connected. To remove this assumption a finer analysis of the corrector functions ψ(F, ·) defined in §6.2.2 is required. At this point one cannot exclude the possibility that a/c methods exist for which this assumption is in fact necessary.
Many open problems remain to be answered. First and foremost, one ought to answer the question whether a/c methods satisfying all the conditions of Theorem 10 always exist. In [45], we present a general construction (a variant on the geometrically consistent coupling method [15]) that appears to work in practise, however, we have no proof of this fact so far. Indeed, if it should turn out that in certain cases the "ghost forces" cannot be completely removed, then an extension of Theorem 10 estimating the contribution of the "ghost force" to the modelling error is highly desirable since such a result would provide the correct quantity that needs to be minimized. It is by no means clear that minimizing the "ghost force" itself is the best possible target. A similar analysis would also be useful for estimating the modelling error of blending methods [55].
It should be conceptually straightforward, though technically more demanding, to generalize all results to higher order finite element methods in the continuum region, however, it would then also be desirable to obtain the second order modelling error estimate in the continuum region. Such a result seems difficult to obtain without a more detailed understanding of the corrector functions ψ(F, ·).
An immediate question is whether a variant of the main result is still valid in 3D. This is by no means clear at this point. From a technical point of view, we require generalizations of the two main technical tools: the bond density lemma ( §5.1) and the characterisation of discrete divergence-free P 0 -tensor fields ( §5.2). While the bond density lemma as stated in this paper is false in 3D, one can establish variants that are potentially usefull for a 3D analysis (work in progress). Generalising the explicit construction of §5.2 is entirely open at this point.
Another important and difficult question is the extension to multi-lattices where the Cauchy-Born model is obtained through a homogenization procedure [1,9].
Finally, it should be stressed, that Theorem 10 is a general abstract result, and as such can undoubtedly be improved upon when a specific coupling method is analyzed. It may be possible for specific methods to obtain more information about the corrector functions ψ(F, ·), and hence obtain a better estimate on the dependence of the modelling error on the interface width. For example, the consistency analysis in [41] requires no corrector functions at all, and in the consistency analysis of nearest-neighbour interactions [45] the corrector function vanishes in the continuum region. The proof of Theorem 10 may, however, serve as a general guidance for modelling error estimates in specific cases.
Finally, the stability of a/c methods in 2D/3D is largely open at this point.
Using the periodic bond density lemma, and the fact that {χ T # | T ∈ T ε } is a partition of unity for R 2 , we have Let ϕ q be the nodal basis function associated with the point value at a vertex q, φ f the nodal basis function associated with the normal derivative at an edge f , and let Φ q,α be the nodal basis function associated with the α-component of the derivative ∇z(q).
Note in particular, that the gradients of these nodal basis functions are scale invariant, that is, where C is a fixed constant that is independent of ε.