From atomistic systems to linearized continuum models for elastic materials with voids

We study an atomistic model that describes the microscopic formation of material voids inside elastically stressed solids under an additional curvature regularization at the discrete level. Using a discrete-to-continuum analysis, by means of a recent geometric rigidity result in variable domains [27] and {\Gamma}-convergence tools, we rigorously derive effective linearized continuum models for elastically stressed solids with material voids in three-dimensional elasticity.


Introduction
In the last years, the understanding of stress driven rearrangement instabilities (SDRI) has attracted huge interest both from the physics and the mathematics community [5,35,36,37,45]. These morphological instabilities of interfaces are generated by the competition between elastic bulk and surface energies, including many different phenomena such as brittle fracture, formation of material voids inside elastically stressed solids, or hetero-epitaxial growth of elastic thin films.
From a mathematical point of view, the common feature of functionals describing SDRI is the presence of both stored elastic bulk and surface energies. In the static setting, problems arise concerning existence, regularity, and stability of equilibrium configurations obtained by energy minimization. In the framework of linearized elasticity, these issues are by now mostly dealt with in dimension two [24,25,37,38] and only recently in dimension three [17,27]. In this work, we focus on the formation of material voids inside elastically stressed solids. At the continuum level and in the regime of linear elasticity, the variational formulation consists in considering functionals defined on pairs of function-set, namely where E ⊂ Ω represents the (sufficiently regular) void set within an elastic body with reference configuration Ω ⊂ R d , and u is the corresponding elastic displacement field. The first part of the functional represents the elastic energy depending on the linear strain e(u) := 1 2 (∇u + (∇u) T ), where C : R d×d → R d×d denotes the fourth-order positive semi-definite tensor of elasticity coefficients, whose kernel contains the subspace of skew-symmetric matrices. The surface energy depends on a possibly anisotropic density ϕ evaluated at the outer unit normal ν E to ∂E ∩ Ω.
From an analytical point of view, in order to guarantee existence of equilibrium configurations for Dirichlet boundary value problems, it is necessary to consider an effective relaxation of the energy given in (1.1). Indeed, even though for fixed E the functional F (·, E) is weakly lower semicontinuous in H 1 (Ω\ E; R d ) and, for fixed u, F (u, ·) can be regarded as a lower semicontinuous functional on sets of finite perimeter with respect to the L 1 -convergence of sets, the energy defined on pairs (u, E) is not lower semicontinuous. This is due to the fact that, along a sequence, the voids may collapse in the limit into a discontinuity of the limiting displacement, which lies in the space of generalized functions of bounded deformation (GSBD) [21]. The relaxation has to take this phenomenon into account, in particular collapsed surfaces need to be counted twice in the relaxed energy. Another interesting and challenging question is to establish a link between the relaxation of (1.1) and the corresponding model of nonlinear elasticity via a simultaneous linearization-relaxation procedure. In that situation, the linear elastic energy is replaced by an integral functional of a nonlinear elastic energy density of the full deformation gradient ∇y (where y : Ω \ E → R d represents the deformation), which is rotationally invariant and is minimized on SO(d). This analysis has been carried out in the physically relevant dimensions d = 2, 3 under the addition of a curvature regularization term in [27], where the main novelty is a quantitative geometric rigidity result for variable domains (see Theorem 2.1 therein).
The scope of this paper is to derive the relaxation of the model in (1.1) as an effective continuum limit of physically relevant atomistic models. We will present here only the case d = 3, since the case d = 2 is completely analogous, and some arguments along the proofs can actually even be simplified. The passage from atomistic to continuum models (see [7]) via Γ-convergence [6,20] has been carried out in various contexts, including elasticity [1,8,10,11,44], fracture [2,9,19,29,30,31,41,43], or more general problems containing free discontinuities [4,39,42].
The analysis in some of these contributions, see [9,31,39], shows that for standard mass-spring models (e.g., governed by Lennard-Jones-type interactions) the effective bulk and surface energies are of the same order only in the case of infinitesimal elastic strains. In this case, the discrete-tocontinuum limits are coupled with a simultaneous passage from nonlinear to linear elasticity. We follow the approach of Schmidt [44], which can be motivated by regrouping the interactions of a mass-spring model, cf. for instance [10]. More precisely, given the reference set Ω, we denote by ε > 0 the lattice spacing of the underlying reference lattice Z ε (Ω) := εZ 3 ∩ Ω and by δ ε > 0 the linearization parameter depending on the lattice spacing, which represents the typical order of the elastic strain and satisfies δ ε → 0 as ε → 0. Here, we would like to point out that our analysis extends to general Bravais lattices. For expository reasons, however, we only consider the case of Z 3 . To a given pair (y, E), where the lattice points E ⊂ Z ε (Ω) represent the presence of voids inside the material at the microscopic level and y : Z ε (Ω) \ E → R 3 represents its deformation, we associate an elastic energy F el ε (y, E) := δ −2 ε i∈Zε(Ω)\E ε 3 W el ε,cell (i, y, E) , (1.2) where the cell energy W el ε,cell (i, y, E) (specified more precisely in Subsection 2.2) represents the elastic energy needed in order to locally deform the neighborhood of i ∈ Z ε (Ω) \ E through the deformation y. Additionally, the presence of voids induces an extra energy contribution of discreteperimeter type due to the missing bonds at the boundary of the configuration, namely, where V ⊂ Z 3 describes the set of neighbors of a lattice point (possibly containing also next-nearest and next-next-nearest neighbors) and c ξ > 0 is the energy cost "per unit area" in the reference lattice with respect to the neighboring point i + εξ. This energy keeps track of the local valence of each atom in E, i.e., the number of atoms missing in order to complete its neighborhood.
In [39], the discrete-to-continuum limit of the energy given by the sum of the two terms in (1.2) and (1.3) has been derived in the two-dimensional setting of elastic thin films, which from a modeling viewpoint, in the continuum setting, corresponds to the simplifying assumption that in (1.1) the void E is given by a supergraph of a function. This analysis also includes the reasoning why the elastic strain needs to be of order δ ε ∼ √ ε, as then (1.2) and (1.3) are of the same order in ε, see also [9,31]. The main ingredient for the proof in [39] is a (by now classical) rigidity estimate, see [33,44], which allows to infer that configurations y with bounded elastic energy have to be close in a suitable quantitative sense to a global rotation on the set Z ε (Ω) \ E. This fact then permits to linearize around the respective rotation and in this way to obtain a linearized elastic energy in the limit. However, the analogous rigidity result in dimension d = 3, or for the more general problem of elastic materials with voids, is missing due to the possibly highly complex geometry of the interfaces between bulk and void. In [27], this issue has been circumvented in a continuum setting by considering an additional curvature regularization term in the energy fuctional, allowing for a novel, refined rigidity estimate (cf. Theorem 2.1 therein). In order to employ this result in the present paper for our discrete-to-continuum analysis, we add an additional discrete curvature term in the functional, of the form The specific form of the curvature term is introduced in Subsection 2.3. It can be understood as a suitable discretization of an integral of the second fundamental form of surfaces on a mesoscopic scale η ε > 0, satisfying ε ≪ η ε ≪ 1. Moreover, we justify the use of this term by presenting an example in the framework of the Embedded-Atom Model (EAM), see Subsection 2.5. This is a semi-empirical, many-atom potential aiming at describing the atomic structure of metallic systems by including a nonlocal electronic correction [22,23]. The total energy of a pair (y, E) is then given by where the terms are introduced in (1.2)-(1.4), respectively. Under suitable scaling assumptions on the curvature energy and under prescribed Dirichlet boundary conditions, we show that the energies (F ε ) ε>0 Γ-converge with respect to an appropriate topology to the relaxation of (1.1). We refer to Subsection 2.4 for the detailed statements, and in particular to (2.14) for the definition of the relaxed energy. Let us highlight that, as an additional feature of possibly independent interest, our effective continuum and linearized limit shows the validity of the so-called Cauchy-Born rule [15,34], namely that, loosely speaking, individual atoms follow the macroscopic deformation gradient. Mathematically, this corresponds to the fact that suitably defined discrete gradients of the atomic displacements reduce to classical gradients in the continuum limit [44]. The proof follows a semi-discrete approach. First of all, to each pair (y, E) we associate a piecewise constant interpolation of y on ε-cubes and a continuum set Q ε (E) being the union of ε-cubes, i.e., Q ε (E) = i∈E (i + ε[− 1 2 , 1 2 ) 3 ). In order to apply the results obtained in [27], it is necessary to modify and smoothen the continuum void set Q ε (E), see Lemmata 3.4-3.5. More precisely, we first replace Q ε (E) by a set E ηε , being a union of η ε -cubes, where η ε > 0 is the mesoscopic scale related to the curvature term (1.4). Afterwards, we smoothen the set E ηε in a suitable way to pass to a smooth setÊ ε . The core of our proof lies in the fact that the specific form of (1.4) ensures that our modifications on the η ε -scale are energetically convenient and lead to sharp estimates on the perimeter and curvature energy ofÊ ε . After having associated to the pair (y, E) an interpolation and a smooth continuum set, we are in the position to apply the compactness result for materials with voids in the continuum setting, see [27,Proposition 3.1]. We point out that this result yields compactness for rescaled displacement fields u := δ −1 ε (y − id) in the space GSBD [21], and that its proof fundamentally relies on the rigidity estimate in [27,Theorem 2.1]. At this point, a further step consists in proving the validity of the Cauchy-Born rule in the setting of GSBD functions. As the analogous result for Sobolev functions [44] is not directly applicable, beforehand we need to perform a delicate approximation of GSBD functions by Sobolev functions, using recent blow-up techniques for GSBD functions [16,14,18,28] and generalizing them to a discrete setting.
With the compactness result and the Cauchy-Born rule in GSBD at hand, the Γ-liminf inequality for the elastic energy follows from standard arguments by employing the strategy devised in [44]. For the surface energy, it is essential to perform the above modifications carefully, in order to have a sharp estimate on the continuum (anisotropic) surface energy in terms of the discrete perimeter and discrete curvature energy. This, together with the lower bound for the elastic energy, allows to conclude the Γ-liminf inequality. The Γ-limsup is performed via a density argument, which allows to reduce the problem to constructing recovery sequences for smooth functions and sets.
The paper is structured as follows. After fixing some notation at the end of the Introduction, we present our model and the main results in Section 2. In particular, we introduce the various terms of the discrete energy functional in Subsections 2.1-2.3. Our main results are stated in Subsection 2.4 and their proofs are the content of Sections 3-6. In particular, in Section 3 we prove the compactness result, while in Section 4 we address the Cauchy-Born rule for discrete symmetrized gradients. Since this result can be of independent interest when dealing with discreteto-continuum problems for elastic materials with surface discontinuities, we formulate and prove it in any dimension. Sections 5-6 are devoted to the proofs of the Γ-liminf and the Γ-limsup inequality, respectively.
Notation. We close the Introduction with some basic notation. Let d ∈ N, d ≥ 2, denote the dimension of the ambient Euclidean space. For the most part of the paper we will focus on the case d = 3, except for Section 4, the content of which is for any dimension d. Given an open bounded Lipschitz set Ω ⊂ R d , we denote by M(Ω) the collection of all measurable subsets of Ω. By A(Ω) we indicate the collection of open subsets of Ω and by A reg (Ω) the subfamily consisting of those open sets E ⊂ Ω such that ∂E ∩ Ω is a (d − 1)-dimensional C 2 -submanifold of R d . Manifolds and functions of C 2 -regularity will be called regular in the following. Given A ∈ M(Ω), we denote by int(A) its interior, by A its closure, and by A c := Ω \ A its complement in Ω. The diameter of A is denoted by diam(A). Moreover, for r > 0 we let (1.5) (1.7) By P : R d×2 d → R d×2 d we define the orthogonal projection onto the space orthogonal to R d×2 d s,t . Let ε > 0 represent the interatomic distance between neighboring atoms in the lattice reference configuration. For A ∈ M(Ω) we set Z ε (A) := εZ d ∩ A. For any i ∈ εZ d , we defineî := i + ε( 1 2 , . . . , 1 2 ). Moreover, for every i ∈ Z ε (Ω) and any set of lattice points E ⊂ Z ε (Ω) (which in the sequel will represent the void set at the microscopic level), we set (1.8) We also define where, for every i ∈ R d and ρ > 0, Q ρ (i) := i+ρ[− 1 2 , 1 2 ) d is the half-open cube in R d of sidelength ρ centered at i. Given also λ > 0 and a cube Q := Q ρ (i), we set Q λρ := Q λρ (i). When ρ = 1, we omit the dependence on the sidelength and simply denote Q(i) := Q 1 (i). For a set of finite perimeter A ⊂ R d [3], we also define the family of cubes with respect to the shifted lattice ε (1.10) For simplicity, we writeQ ε :=Q ε,R d . We also introduce cubes of sidelength 2ε, namelŷ The following fact will be used several times in the sequel: for every η > 0 and i ∈ ηZ d , it holds that Finally, unless stated otherwise, 0 < c ≤ 1 ≤ C will denote generic dimensional constants, whose values will be allowed to vary from line to line.

Setting and main results
2.1. Admissible deformations in the discrete setting and discrete gradients. Adopting the notation that was introduced in the previous subsection, we denote the void set at a microscopic level by E ⊂ Z ε (Ω). We consider the associated collection of atoms Z ε (Ω) \ E and corresponding (discrete) deformations y : Z ε (Ω) \ E → R d . We start by introducing admissible deformations satisfying Dirichlet boundary values in a specific sense. In order to define (thickened) Dirichlet boundary values for discrete mappings, we introduce some further terminology and notation. First of all, we make the following geometrical assumption on the Dirichlet boundary ∂ D Ω: there exists a decomposition ∂Ω = ∂ D Ω ∪ ∂ N Ω ∪ N with (the outermost boundary has to be understood in the relative sense with respect to ∂Ω), and there existσ > 0 small enough and x 0 ∈ R d such that for all σ ∈ (0,σ) it holds that where O σ,x0 (x) := x 0 + (1 − σ)(x − x 0 ). (These assumptions are related to a suitable density result for displacements, see [17,Lemma 5.7].) In order to fix boundary values, we introduce an open set U ⊃ Ω such that U and U \ Ω are Lipschitz sets and U ∩ ∂Ω = ∂ D Ω.
The definition in (ii) is for definiteness only. Given y ∈ Y ε (y 0 , ∂ D Ω, E), recalling the notation in (1.6), we define the discrete gradient of y at each point i ∈ Z ε (R d ) as Moreover, we setē(y)(i) := P (∇y(i)) for the orthogonal projection of ∇y onto the orthogonal complement of R d×2 d s,t , cf. (1.7). In the following, when no confusion arises, we will also identify ∇y andē(y) with their piecewise constant interpolation, being equal to ∇y(i) andē(y)(i) respectively, on each Q ε (î). We now focus on the case d = 3 and proceed by introducing the discrete energy. An example fitting to our assumptions is given at the end of the section in Subsection 2.5.

2.2.
The elastic cell energy. Following [44], the elastic energy of y is defined in terms of a discrete elastic cell energy associated to each lattice point i ∈ Z ε (Ω) \ E. In particular, given where W el ε,cell denotes the cell energy specified below in (2.4) and δ ε > 0 describes the typical order of the elastic strain. The factor ε 3 corresponds to a bulk scaling of the elastic energy. We suppose that δ ε → 0 as ε → 0, i.e., in the discrete-to-continuum limit, we simultaneously pass to the limit of infinitesimal strains. The most natural choice for δ ε in discrete energies featuring bulk and surface terms is δ ε ∼ √ ε, see e.g. mass-spring models with interaction potentials of Lennard-Jones type as mentioned in the Introduction.
The main assumption is that the elastic cell energy can be split into a bulk and a surface part: for each i ∈ εZ 3 \ E the cell energy is defined by (iii) W el bulk is C 3 -regular in a neighborhood ofSO(3) and the Hessian Q bulk := D 2 W el bulk (Z) at the discrete identity Z is positive-definite on the orthogonal complement of the subspace R 3×8 s,t , see (1.7); (iv) W el bulk grows at infinity at least quadratically on the orthogonal complement of the subspace spanned by infinitesimal translations, i.e., surf depends only on the positions of atoms related to I ε (i, E) in the reference configuration, i.e., W el surf (I ε (i, E), F ), for F ∈ R 3×8 , depends on the second variable only through the columns F l , l ∈ I ε (i, E); (vi) The discrete identity Z is an equilibrium configuration also for W el surf . In particular, we assume that there exists a constant C > 0 such that for every I ⊂ {1, . . . , 8} it holds that W el surf (I, ·) ≤ CW el bulk (·) in a neighborhood ofSO (3).
Note that the assumptions (i) and (iii) imply that Q bulk (v, . . . , v) = 0 and Q bulk (Az 1 , . . . , Az 8 ) = 0 for all v ∈ R 3 and A ∈ R 3×3 skew . Assumption (v) ensures that the atomic positions at lattice points in (Z ε (R 3 ) \ Z ε (U )) ∪ E do not affect the elastic energy, i.e., (2.1)(ii) is indeed for definiteness only. The last condition in (vi) was referred to as compatibility condition between the surface and the bulk elastic energy density in [44, Definition 2.2].
2.3. The surface energy. Let V ⊂ Z 3 \ {0} be a finite set of vectors describing the neighbors of a point on the lattice. We assume that An example of such a set of neighbors is V = 3 l=1 µ l e l : µ l ∈ {−1, 0, 1} . This would correspond to a model taking into account nearest, next-nearest, and next-next-nearest neighbor interactions.
For any E ⊂ Z ε (Ω) and A ∈ M(R 3 ), we consider an anisotropic discrete perimeter energy, consisting of pair-interactions among neighboring atoms subordinate to the set V . More precisely, we define where U is the auxiliary set introduced in Definition 2.1 in order to incorporate Dirichlet boundary conditions for the discrete deformations. Here, c ξ > 0 and c ξ = c −ξ for all ξ ∈ V . Since E ⊂ Z ε (Ω), it is clear by the definition that F per ε (E, U ) = F per ε (E, Ω). Thus, when A = U , we omit the dependence on U in the second variable. Note that positive contributions correspond to the situation that an atom at position i + εξ ∈ Z ε (U ) \ E in the reference configuration misses a neighbor i ∈ E lying in the void set E. The factor ε 2 corresponds to the surface scaling.
We complement the discrete perimeter energy with an additional discrete curvature energy, associated to each pair (E, A). This should be thought of as a discrete version of the continuum curvature energy considered in [27,Equation (2.1)]. It will allow us to implement the compactness result [27, Proposition 3.1], obtained through the rigidity result in Theorem 2.1 therein. The discrete curvature energy will be considered at a mesoscale η ε > 0, satisfying η ε /ε ∈ N, η ε → 0, and η ε /ε → +∞ as ε → 0 . (2.6) More details why η ε indeed needs to be chosen as a mesoscale are given below (2.8). Before introducing the specific form of the energy, we will first discuss which configurations of atoms will be considered as flat, in the sense that they will not pay any curvature energy.
Definition 2.2 (Local flatness). Let i ∈ εZ 3 and η ∈ [η ε , +∞). We say that E ⊂ Z ε (Ω) is locally flat in Q η (i) (with respect to U ) if and only if for all Q := Q 2ε (k) ∈Q ε,2 , k ∈ εZ 3 , such that Q η (i) ∩ Q = ∅, one of the following is true (see Figure 1):  We point out that locally flat sets in a cube Q η (i), η ∈ [η ε , +∞), i ∈ εZ 3 , are locally simple coordinate laminates. We refer to Figure 2 for an illustration, Definition 3.1 below for an exact definition, and Lemma 3.2 for a proof of this fact. Due to the fact that η/ε ≥ η ε /ε → +∞ as ε → 0, the number of these laminates inside the cube Q η (i) can possibly tend to infinity, and hence there can be an arbitrarily large number of possibilities for the geometry of the locally flat set E inside Q η (i). In view of this fact, we will introduce sets that have a bounded number of possibilities for their geometry inside Q η (i), independently of ε.  2.3 (Cubic sets). Let i ∈ εZ 3 and η ∈ [η ε , +∞). We say that E ⊂ Z ε (Ω) is cubic in Q η (i) (with respect to U ) if and only if there exist i 0 ∈ Z ε (Q ηε ) and {k 1 , . . . , k n } ⊂ i 0 + η ε Z 3 (possibly empty) such that For an illustration of cubic sets we refer to Figure 3. Note that cubic sets can be locally flat but might also not. Definition 2.4 (Flatness). Let i ∈ εZ 3 and η ∈ [η ε , +∞). We say that E ⊂ Z ε (Ω) is flat in Q η (i) (with respect to U ) if and only if it is cubic and locally flat in Q η (i).
It is elementary to check that flat sets in Q η (i) are a halfspace intersected with Z ε (Q η (i)), see the last example in Figure 3, i.e., in such a case, there exists x 0 ∈ Q η (i) and ν ∈ {±e 1 , ±e 2 , ±e 3 We now introduce a discrete version of the curvature regularization energy that we used in [27], expressed through a curvature cell energy at the mesoscale η ε . In particular, let q ∈ [2, +∞) be fixed from now on. Given any E ⊂ Z ε (Ω) and A ∈ M(R 3 ), we define We again omit the dependence on the set in the second variable for A = U . Our structural assumptions on the curvature cell energy are the following: (i) (Flatness) W curv ε,cell (i, E) = 0 if E is flat in Q ηε (i) ; (ii) (Lower bound for local non-flatness) There exist c > 0 and γ ε > 0, γ ε → 0 as ε → 0, such that W curv ε,cell (i, E) ≥ cγ ε η −1−q ε if E is not locally flat in Q ηε (i) ; (iii) (Upper bound for cubic sets) There exists Let us comment on the scaling of the energy in (ii) and (iii): roughly speaking, local non-flatness corresponds to an energy per atom of order γ ε η −1−q ε , where γ ε represents a curvature regularization parameter. In this case, in a cube Q ηε (i), i ∈ εZ 3 , the overall energy F curv ε (E, Q ηε (i)) is typically of the order γ ε η 2−q ε , since #Z ε (Q ηε (i)) = η 3 ε /ε 3 . Up to the prefactor γ ε , this corresponds to the integral of q-th power of the second fundamental form of a round sphere in R 3 with radius η ε . This effectively relates our choice of discrete curvature to classical continuum curvature notions.
Remark 2.5. Note that, since E ⊂ Z ε (Ω), the Definitions 2.2-2.4 as well as the curvature energy F curv ε (E, U ) are independent of the choice of the auxiliary set U for ε > 0 sufficiently small (and thus η ε small accordingly, see (2.6)).
Remark 2.6 (Different assumptions). For later purposes in the example of Subsection 2.5, let us comment that our results are still valid if in (ii) of the above assumptions, the cubes Q ηε (i) are replaced by smaller cubes Q η (i) for η ∈ [η ε /2, η ε ] with η/ε ∈ N. Our choice of the exact sidelength η ε is only for expository reasons, in order to formulate the above assumptions in a simpler way.
We also assume that the linearization parameter δ ε , the mesoscale η ε , and the curvature regularization parameter γ ε are related to each other via the rates The last condition is [27,Equation (3.4)], and is required in order to apply the compactness result [27, Proposition 3.1]. The first two conditions, which relate the mesoscale on which the curvature energy is defined and the curvature regularization parameter, are necessary in order to suitably modify the void sets at ε-scale to geometrically and energetically more convenient sets at the η εscale. More precisely, the first condition allows for sharp bounds on the perimeter and curvature energy, as well as the cardinality of the modified sets in terms of the energies and cardinality of the original ones, see Lemmata 3.4 and 3.5. The second condition is needed to ensure that the curvature energy for recovery sequences is negligible.
Our goal is to identify the effective continuum limit of the above discrete energies as ε → 0. It turns out that the limiting energy functional coincides with the relaxed, linearized model for material voids in elastically stressed solids studied in [17]. First of all, we introduce the continuum surface energy density ϕ : for ν ∈ R 3 , where V was defined before (2.5). Given u ∈ GSBD 2 (Ω) (see [21] for the definition and details on this function space) and a set of finite perimeter E ∈ M(Ω), we first define the boundary energy term by which is non-trivial if the void goes up to the Dirichlet part of the boundary, or the mapping u does not satisfy the imposed boundary conditions. Here, ν E denotes the measure-theoretical outer unit normal to ∂ * E, ν Ω denotes the outer unit normal to ∂Ω, and tr(u) indicates the trace of u at ∂Ω, which is well-defined for functions in GSBD 2 (Ω), see [21,Section 5].
The limiting linearized elastic energy for pairs (u, E) as above is defined by where e(u) := sym(∇u) is the approximate symmetric gradient of u and Q bulk := D 2 W el bulk (Z), see Assumption (iii) in Subsection 2.2. The limiting surface energy is given by where ν u is the approximate unit normal vector to J u . The total limiting energy is defined by In order to formulate the main result of this paper, we proceed with the definition of convergence for sequences of discrete displacements and void subsets of lattice points. To that end, given E ⊂ Z ε (Ω) and a discrete displacement u ∈ U ε (u 0 , ∂ D Ω, E), we denote byū the piecewise constant interpolation of u, being equal to For the next definition, we denote by L 0 (Ω; R 3 ) the space of measurable maps in Ω (equipped with the topology of convergence in measure) and we recall also the notation introduced in (1.9).
Note that this convergence of (χ Qε(Eε)∩Ω ,ū ε ) ε>0 is referred to as τ -convergence in [27], see the paragraph below [27, Proposition 3.1]. There, also the necessity of the set ω u is discussed, which here is related to the fact that Z ε (U ) might be disconnected into connected components by E ε (in a graph-theoretical sense), and on the components not intersecting the Dirichlet boundary, the behavior of the displacement cannot be controlled.
After these preparations, we now present the main results of the paper.
(i) Then, there exists a subsequence (not relabeled) converging to a pair (u, E) in the d-sense, where u ∈ GSBD 2 (Ω) and E ∈ M(Ω) is a set of finite perimeter.
(ii) Let ω u be the set of finite perimeter in Definition 2.7. There exist sets of finite perimeter and The property stated in (ii) above can be understood as a manifestation of the Cauchy-Born rule for symmetrized gradients: for a given linear macroscopic deformation gradient F , the corresponding symmetrized discrete gradient is given by sym(F )Z. Note that this rule does not enter as an assumption, but is rather a consequence of our analysis. Theorem 2.9 (Discrete-to-continuum Γ-convergence). Under the above assumptions, the sequence of functionals (F ε ) ε>0 Γ-converges to F 0 as ε → 0, with respect to the convergence d.
The proofs of the above theorems will be given in the next sections. Concerning the elastic energy, we will follow the strategy devised in [44], combined with the compactness result in [27, Proposition 3.1], which uses the novel rigidity estimate in [27, Theorem 2.1]. A further ingredient is the Cauchy-Born rule for symmetrized discrete gradients, see Section 4. For the surface part, our idea relies on replacing the discrete void set by a continuum representative, and then using the result in the continuum version, i.e., [27,Theorem 3.2].
Finally, we want to stress that the form (2.10) of the surface energy density is justified by the curvature energy term: the latter ensures that energetically convenient configurations will locally be only along coordinate directions. Thus, no microscopic relaxation takes place and it suffices to calculate the energy density for half-spaces with coordinate vectors as outer unit normal, see the following figure.
ν ε η ε Figure 4. Two microscopic approximations of an interface with normal ν = (−1/ √ 2, 1/ √ 2, 0) to illustrate that no microscopic relaxation takes place. The one on the left, although in general having less perimeter energy at the ε-level compared to the right, has high curvature energy and therefore is not energetically convenient in our model. Note that for specific interaction ranges V , for the one on the right, the interactions in direction ξ = (1, 1, 0) have a positive contribution to the energy, while for the one on the left they do not. The sidelength of the large cubes is ηε, so that the configuration on the right is cubic inside the cube.

2.5.
Example. We close this section with an example of a model with nearest-neighbors and nextnearest-neighbors atomic interactions to which our main results apply. For the elastic cell energy, given K 1 , K 2 > 0, we consider the bulk cell energy The nonnegative term χ is supposed to be nonzero only for discrete gradients which are not locally orientation preserving. In particular, χ is zero in a neighborhood ofSO(3) = SO(3)Z and larger than some The surface cell energy W el surf is chosen appropriately. In [44,Subsection 5.3], it is discussed that this cell energy satisfies all assumptions given in Subsection 2.2.
For the discrete perimeter energy, we consider V = {±e k : k = 1, 2, 3}∪{±e k ±e l : 1 ≤ k < l ≤ 3} and choose any c ξ = c −ξ > 0 for ξ ∈ V . The sum of the elastic and the perimeter energy can be understood as an idealization of interactions of Lennard-Jones-type, cf. [9]. Finally, we show that we can choose a curvature cell energy in the framework of the Embedded Atom Model (EAM) satisfying the conditions (i)-(iii) below (2.7).
Atomic positions induce electronic-cloud distributions. In the EAM, this is modeled by a multibody interaction term of the following form: Here, ρ : R + → R + models the electron-cloud contribution of an atom placed at j on an atom placed at i. The sum ρ i describes the cumulative effect on the atom placed at i of the electronic clouds related to all other atoms. Eventually, the function Ψ ε : R + → R + describes the energy needed to embed an atom at position i in the host electron gas created by the other atoms. For the sake of simplifying the exposition, we prefer to consider an example for d = 2 and for Ω = R 2 . We consider With this particular choice, which corresponds to a two-dimensional EAM where only the cardinality of the nearest-neighbors and next-nearest-neighbors are taken into account, we observe that we can find G ε : In fact, with our choice of ρ we have that Hence, there is a one-to-one correspondence between G ε and Ψ ε . We now choose and define Then, recalling that η ε /ε ∈ N, by changing the order of summation, since for fixed This shows that W curv ε,cell can be chosen such that F curv ε coincides with Φ ε given in (2.17).
Proof. Condition (i) follows from the fact that flat configurations inside a square Q ηε (i) are locally half-spaces intersected with Z ε (Q ηε (i)) (see Picture 2, as well as Lemma 3.2), which corresponds to neighborhood cardinality (3, 2) (at the flat interface) or (4, 4) (inside the set).
Finally, we show (ii) for η ε − 4ε in place of η ε . To this end, in view of (2.18) and (2.19), it By contradiction we assume that The main step of the proof is then to prove that (up to a single rotation of a multiple of π/2 independently of j) for all j ∈ Z ε (Q ηε (i)) \ E, see also Figure 5, where the two possible configurations (up to rotation) are depicted. Indeed, (2.21) would then contradict the assumption that E is not locally flat in Q ηε−4ε (i).
j j  Figure 6. The different cases of atoms j having low energy but having an atom nearby with high energy, illustrating cases (a)-(c) above, respectively.

Compactness
The goal of this section is to prove Theorem 2.8(i). Before we come to the actual proof, we present some auxiliary lemmata which will be essential for the following. Let η ε > 0 be the mesoscale introduced in Subsection 2.3, which satisfies (2.6) and (2.8). Let Ω ⊂ U be bounded open Lipschitz sets in R 3 , as introduced before Definition 2.1, and recall also Remark 2.5. Definition 3.1. Let A ∈ M(U ) and let E ⊂ Z ε (Ω). We say that E is a simple coordinate laminate in A if there exists ν ∈ {e 1 , e 2 , e 3 } and a function h : εZ → {0, 1} such that for every i ∈ Z ε (U ), Proof. We divide the proof into three steps.
Step 1.(Laminate structure of the "interface cubes") Assume that E ∩ Q η (i) = ∅, since otherwise we can choose h ≡ 0 and ν = e 1 . Let Q : Then, by Definition 2.2, up to rotation, we may assume that We show that for all q : We prove this property inductively. First, assume that for l : Indeed, assume by contradiction that this was not the case. Then, since E is locally flat in Q η (i) and {l, l + εe 2 } ⊂ E ∩ Q 2ε (l) by (3.1), one of the following three possibilities holds true: However, due to the local flatness of E in Q η (i) and (3.1), (a) and (b) are not possible, since in both of these cases the local flatness of E with respect to Q would be violated. Therefore, (c) applies, i.e., (3.3) holds true. The same argument shows that for m : Arguing inductively in exactly the same manner, we indeed get (3.2).
Step 2.(One normal to the interface) Next, we show the following, up to a rotation of a multiple of π/2 with respect to an axis passing through a point in the shifted lattice ε(Z 3 + ( 1 2 , 1 2 , 1 2 )) in the direction of one of the coordinate vectors: for all n ∈ εZ 3 such that forQ : Without restriction we assume by contradiction that for j = 1, 2 there exist cubesQ j : E ∩Q 2 = {n 2 , n 2 + εe 1 , n 2 + εe 3 , n 2 + ε(e 1 + e 3 )} .
Step 3.(Conclusions) We claim that (up to a possible rotation of a multiple of π/2 as before) This finishes the proof of the lemma.
In what follows, we will use the following fact several times: This follows from (2.6). As a first step towards the proof of the compactness result, we replace a set E ⊂ Z ε (Ω) by a set that essentially can only vary its shape at the mesoscale η ε . This is formulated in a local version in the next lemma.
Lemma 3.4 (η ε -scale replacement). There exists a constant C = C({c ξ } ξ∈V ) > 0 such that for every E ⊂ Z ε (Ω), the set E ηε ⊂ Z ε (Ω) given by (3.7) satisfies the following. For every open subset U ⊂⊂ U with Lipschitz boundary, there exists ε 0 = ε 0 ( U , U ) ∈ (0, 1), so that for every ε ∈ (0, ε 0 ] the following three properties are satisfied: where ω(ε) → 0 as ε → 0. (ii) The discrete curvature energy of E ηε , see (2.7), satisfies the estimate Proof. In the proof, C > 0 denotes a generic constant that possibly depends on {c ξ } ξ∈V and may vary from line to line. We will assume without loss of generality that E ηε ∩ U = ∅, otherwise there is trivially nothing to prove. In the following, we suppose that ε 0 = ε 0 ( U , U ) ∈ (0, 1) is chosen sufficiently small so that all our subsequent statements are valid for every ε ∈ (0, ε 0 ]. In particular, we assume that (3.8) holds true. We divide the proof into several steps. We first prove (3.9) by distinguishing two different types of cubes (depending on the local geometry of E) and showing the estimate for each one of them separately (Step 1). Then we prove (3.10) ( Step 2) and finally (3.11) ( Step 3). The two subfamilies of cubes in Q ηε, U are defined as Step 1.(Perimeter Estimate) In order to obtain estimate (3.9), we treat the good and bad cubes separately. We first start with the case of bad cubes. Step does not satisfy any of the properties (i)-(iii) of Definition 2.2. We actually have that k ∈ Z ε Q 2ηε (i) , as long as we take ε 0 > 0 small enough. Now, for all j ∈ Z ε Q ηε 4 (k) we have that Q 2ε (k) ⊂ Q ηε 2 (j) and therefore, due to assumption (ii) of the curvature energy (2.7) (see also Remark 2.6), we get (3.14) Notice that, since j ∈ Z ε Q ηε 4 (k) and k ∈ Z ε Q 2ηε (i) , we get that j ∈ Z ε Q 3ηε (i) . This, together with the simple fact that #Z ε (Q ηε 4 (k)) ≥ C η 3 ε ε 3 (for an absolute constant C > 0) and (3.14), allows us to deduce that On the other hand, we claim that Indeed, due to the assumptions on the set of neighbors V , see the beginning of Subsection 2.3, we have #V ≤ 27 and therefore for ε 0 ∈ (0, 1) sufficiently small, for every ε ∈ (0, ε 0 ] (and hence also for η ε > 0 small accordingly) we can estimate where we used the fact that 0 < ε ≤ η ε , and that by the construction of the set E ηε , missing neighbors can only occur at an ε-layer of atoms around ∂Q ηε (i). The last constant C > 0 depends on {c ξ } ξ∈V here. This implies (3.16). Noting that, due to (2.8), we have that γ ε η −q ε → +∞ as ε → 0, and using (3.15) together with (3.16), we obtain (3.13) for ω(ε) := Cγ −1 ε η q ε . This concludes Step 1.1. Next, we treat the case of good cubes. Step By (2.5) and (3.7) we can assume without restriction that E ∩ Q ηε (i) = ∅. Additionally, due to the definitions in (3.12), E is a simple coordinate laminate in Q 3 2 ηε (i) with respect to a unit normal vector ν ∈ {e 1 , e 2 , e 3 }. We assume without loss of generality that ν = e 3 . Since E ∩ Q ηε (i) = ∅, we have that E ηε ∩ Q ηε (i) = Z ε (Q ηε (i)). As E is a simple laminate in Q 3 2 ηε (i) with respect to e 3 , also E ηε is a simple laminate in Q 3 2 ηε (i) with respect to e 3 . Thus, we have that in the formula (2.5) for F per ε (E ηε , Q ηε (i)), the addends ξ ∈ V yielding a positive contribution satisfy ξ 3 = 0. Without restriction we can consider the case ξ 3 > 0. Let j ∈ Z ε (Q ηε (i)) and ξ ∈ V be such that Let now k ∈ E ∩ Q ηε (i) be the point with maximal e 3 -coordinate among all other points in E ∩Q ηε (i) in the same e 3 -column as j, i.e., k 1 = j 1 , k 2 = j 2 , and m 3 ≤ k 3 ≤ j 3 for all m : Therefore, due to the choice of k, we again have that k + εξ / ∈ E. Note that in the previous procedure, since ξ ∞ ≤ 1, for each ξ ∈ V we have that each k ∈ E ∩ Q ηε (i) (with k + εξ / ∈ E) is chosen at most once for a point j ∈ Z ε (Q ηε (i)) ∩ E ηε with j + εξ / ∈ E ηε . Therefore, in (2.5), summing over all j for E ηε and the corresponding k for E, we obtain (3.17).
Step 2.(Curvature Estimate) Next, we prove (3.10). This follows from the property that for every Q ηε (i) ∈ Q ηε, U we have for an absolute constant C > 0. Indeed, once the above inequality is shown, by summing over all Q ηε (i) ∈ Q ηε, U , (3.10) follows from (1.12), (3.8), and (3.18). We now show (3.18). To this end, let Q ηε (i) ∈ Q ηε, U and assume without loss of generality that E ηε is not flat in Q 2ηε (i), since otherwise the statement is trivial by (2.7) and assumption (i) on W curv ε,cell . Since E ηε is cubic in Q 2ηε (i), this also implies that E ηε is not locally flat in Q 2ηε (i), see Definition 2.4. We can also assume that E is not locally flat in Q 5 2 ηε (i), since otherwise E ηε would be flat in Q 5 2 ηε (i), and again the statement would be trivial. Consequently, as E ηε is not locally flat in Q 5 2 ηε (i), we have that, up to rotation, there exists ν ∈ R 3 with ν i ∈ {0, 1} for all i = 1, 2, 3, such that E ηε ∩ Q 2ε (i + ηε 2 ν) does not satisfy any of the conditions (i)-(iii) of Definition 2.2. Moreover, since E is not locally flat in Q 5 2 ηε (i), there exists Q 2ε (k) ∈Q ε,2 such that Q 5 2 ηε (i) ∩ Q 2ε (k) = ∅, but in Q 2ε (k) the set E does not satisfy any of the conditions (i)-(iii) of Definition 2.2. We can now proceed similarly to Step 1.1: for all j ∈ Z ε (Q ηε 4 (k)) we have that Q 2ε (k) ⊂ Q ηε 2 (j) ⊂ Q 3ηε (i) and therefore, due to property (ii) below (2.7) (see also Remark 2.6), we have As before in Step 1.1, along with the fact that #Z ε (Q ηε On the other hand, due to (3.7), we have that E ηε is cubic in Q ηε (j) for all j ∈ Z ε (Q ηε (i)), and therefore by property (iii) below (2.7) and the fact that #Z ε (Q ηε (i) (3.20) Hence, (3.19) and (3.20) imply (3.18). This concludes Step 2.
For a setÊ ⊂ U of finite perimeter, we will also denote by its anistropic perimeter in A ⊂ U with respect to the density ϕ defined in (2.10). Again, if A = U , we omit the dependence on A in the notation. Moreover, for a regular setÊ ∈ A reg (U ), we denote by A the second fundamental form of ∂Ê ∩ U , i.e., |A| = κ 2 1 + κ 2 2 , where κ 1 and κ 2 are the principal curvatures of ∂Ê ∩ U . As a second preliminary step, we replace the set E ηε by a regular set, whose perimeter and curvature energy can be controlled by the discrete perimeter and curvature energy of E ηε , respectively. This is formulated, again in a local fashion, in the following lemma. (i) For the anisotropic perimeter ofÊ with density ϕ, we have

23)
where ω(ε) → 0 as ε → 0. (ii) For the curvature energy ofÊ, we have (3.24) (iii) For the difference in volumes, we have Proof. We divide the proof into several steps. In the first step, we explain how to reduce the problem to the case that E is an η ε -scale set in the sense of (3.7). In Step 2 we show how to construct the setÊ for η ε = 1. In Step 3 we discuss the shape of the interface of the smooth replacement set and derive the estimates in the case that η ε = 1. In Step 4 we perform a scaling argument and derive the corresponding estimates for general η ε > 0, which allows us to conclude.
Step 1.(Reduction to η ε -scale sets) First of all, without loss of generality we can assume that E is an η ε -scale set in the sense of (3.7). To this end, let U ′ ⊂⊂ U be an open Lipschitz set such that U ⊂⊂ U ′ . Let ε 0 > 0 be small enough such that Lemma 3.4 is applicable for U ′ in place of U . We apply Lemma 3.4 for U ′ to obtain E ηε such that E ⊂ E ηε and (3.9)-(3.11) hold true for E ηε . We perform the construction presented in Steps 2-4 below for E ηε in place of E and for U ′ in place of U to get (3.23)-(3.26) for E ηε in place of E. Then, we see that (3.23)-(3.26) also hold for E due to (3.9)-(3.11), E ⊂ E ηε , and . Therefore, from now on we can assume that E is an η ε -scale set.
Step 2.(Construction of the set at scale η ε = 1) Let E ⊂ Z 3 be a finite set and (according to the notation we introduced in (1.9) and the comments thereafter) let Q(E) := i∈E Q(i) . Let (ζ σ ) σ>0 be a family of smooth symmetric mollifiers in R 3 , i.e., there exists a radially symmetric function ζ ∈ C ∞ c (R 3 ; [0, +∞)) such that For every σ > 0 consider the function f σ (x) := χ Q(E) * ζ σ (x) and for t ∈ (0, 1) consider the super-level set E σ,t defined as E σ,t := {x ∈ R 3 : f σ (x) > t}. (Here, for notational convenience, we have suppressed the dependence of f σ on E.) Our goal is to choose c 0 > 0, σ 0 > 0, and t 0 ∈ (0, 1) independently of E such that (i) ∂E σ0,t0 is a regular 2-dimensional manifold; To see these properties, we define the N 0 := 2 27 different sets (E k ) k=1,...,N0 by E k := j∈I k Q(j), I k ⊂ {−1, 0, 1} 3 . We note that for each i ∈ Z 3 we find k such that In other words, there are only N 0 different possibilities for the shape of for all x ∈ Q 3 2 (i). Now let N k be the set of critical values of f k σ0 , which by the Morse-Sard Lemma (cf. Lemma 13.15 in [40]) is a set of L 1 -measure zero. In particular, for N := N0 k=1 N k we have that L 1 (N ) = 0. Thus, we can choose t 0 ∈ (0, 1) \ N such that (i) is satisfied. In order to obtain (iii), due to (3.28), it suffices to choose t 0 and c 0 > 0 small enough depending on ζ σ0 such that (iii) is satisfied for all E k , k ∈ {1, . . . , N 0 }, in place of E. Then, as E σ0,t0 ⊃ Q(E), (ii) follows from the fact that supp ζ σ0 ⊂⊂ B σ0 .
Step 3.(Surface energy estimate in the case η ε = 1) Similarly to the proof of Lemma 3.4, we will consider the two families of good and bad cubes with respect to E, namely Let us denote by the collection of boundary cubes. For a cube Q := Q(i) ∈ ∂Q, let us denote by T Q the (nonempty) collection of those faces of Q that also belong to ∂Q(E). For each T ∈ T Q , let ν T ∈ {±e i } i=1,2,3 be the constant outward pointing unit normal to ∂Q(E) at the face T . For the set E σ0,t0 we denote its outward unit normal at a point x ∈ ∂E σ0,t0 by ν σ0,t0 (x).
Step 3.1.(Good cubes) Let Q = Q(i) ∈ Q good E ∩ ∂Q. Fix a boundary face T ∈ T Q . Provided that σ 0 > 0 is small enough, we can check that E σ0,t0 is also a coordinate laminate in (T ) 2σ0 , (cf. (1.5)), in the sense that ν σ0,t0 (x) = ν T for every x ∈ ∂E σ0,t0 ∩ (T ) 2σ0 . Indeed, suppose without restriction (after possibly a rotation by a multiple of π 2 with respect to a coordinate axis and a reflection) that ν T = e 3 . Then, since ∂E σ0,t0 = {f σ0 = t 0 } and ∇f σ0 is nowhere vanishing on ∂E σ0,t0 (this is again a consequence of the choice of σ 0 and the Morse-Sard Lemma), for every x ∈ ∂E σ0,t0 ∩ (T ) 2σ0 and for l = 1, 2, 3, we have where ν l denotes the l-th component of the outer unit normal to ∂(Q(E) ∩ B σ0 (x)), which satisfies ). Due to (3.27), we have ζ σ0 (x − y) = 0 on ∂B σ0 (x). On the other hand, as a consequence of the laminate structure of E in Q 3 (i), for σ 0 > 0 small enough, we get that ∂(Q(E)) ∩ B σ0 (x) is either empty or is contained in the union of T and the other 8 faces with outer normal e 3 of neighboring cubes to Q(i) that share either an edge or a vertex with T . (Note that these neighboring cubes also belong to Q(E).) Therefore, we have that ν l = δ 3l on ∂(Q(E)) ∩ B σ0 (x), where δ 3l denotes the Kronecker delta. In view of (3.31), this implies For good boundary cubes Q ∈ Q good E ∩ ∂Q, the laminate structure implies that #T Q ∈ {1, 2} and if #T Q = 2, then two opposite faces are in T Q . Therefore, the outer unit normal vector to ∂Q(E) of boundary faces is well-defined up to sign, and denoted by ν Q . We define (3.32), the fact that ν T = ±ν Q , and property (ii) established in Step 2, we have Hence, we obtain while, by the flat laminate structure of E σ0,t0 in (Q) vert 2σ0 , it is clear that Step 3.2.(Bad cubes) Let Q = Q(i) ∈ Q bad E . As noted in (3.28), there are at most N 0 different possibilities for the local shape of E. In particular, there exists a universal constant C > 0 (only depending on σ 0 , t 0 , and these N 0 different configurations (E k ) k=1,...,N0 ) so that for every Q ∈ Q bad E , We also note that ∂E σ0,t0 is covered by the union of ∂E σ0,t0 ∩ (Q) vert 2σ0 , Q ∈ Q good E ∩ ∂Q, together with ∂E σ0,t0 ∩ Q 3 2 (i), Q(i) ∈ Q bad E .
Step 4.(η ε -scale construction) We now use the previous steps in order to construct a regular setÊ from a set E ⊂ Z ε (Ω). By Step 1, we can assume that E is an η ε -scale set in the sense of (3.7). Now, let E 1 := η −1 ε E, which is a finite subset of Z 3 . We then choose (E 1 ) σ0,t0 as constructed in Step 2 and defineÊ byÊ := η ε (E 1 ) σ0,t0 ∩ U . By standard properties of scaling of sets, we have Q ε (E) = η ε Q(E 1 ) (recall (1.9)), as well as The last property directly implies (3.26), and we also get Q ε (E) ∩ U ⊂Ê ⊂ U . Now, denote by where we have used the same notation as in (3.29), (3.30) for the families Q good E1 , Q bad E1 , ∂Q (with respect to E 1 ), and for notational convenience we have suppressed the dependence on the original set E in the notation for the families Q good ηε, U and Q bad ηε, U . By (3.33), (3.34), and (3.35), we have On the other hand, due to (3.36) and elementary scaling properties of H 2 and the second fundamental form, for Q := Q ηε (i) ∈ Q bad ηε, U , we have Finally, by the simple observation we made at the end of Step 3 and by scaling, we get Step 4.1.(Proof of (3.23)) By (3.41) we can estimate We start with the good cubes. We use the same notation T Q and ν T as before in where we also used (3.8) and the fact that η ε /ε ∈ N, see (2.6). Now, let Q ∈ Q bad ηε, U . By the negation of Lemma 3.2 there exists k ∈ εZ 3 such that for Q 2ε (k) ∈Q ε,2 , Q 3ηε (i) ∩ Q 2ε (k) = ∅ and E ∩ Q 2ε (k) does not satisfy any of the properties (i)-(iii) of Definition 2.2. By (2.6), we actually have that k ∈ Z ε Q 7 2 ηε (i) , provided that we take ε 0 > 0 small enough. Arguing as in the proof of Lemma 3.4, for all j ∈ Z ε Q ηε 4 (k) we have that Q 2ε (k) ⊂ Q ηε 2 (j) and therefore, due to property (ii) below (2.7), we get Notice that, since j ∈ Z ε Q ηε 4 (k) and k ∈ Z ε Q 7 2 ηε (i) , we get that j ∈ Z ε Q 4ηε (i) . This, together with the fact that #Z ε (Q ηε 4 (k)) ≥ C η 3 ε ε 3 for an absolute constant C > 0, leads to which proves (3.24), where the last step follows as in (3.45).
Step 4.3.(Proof of (3.25)) By the construction of the setÊ we have that where we recall the notation in (3.37), and set ∂Q ηε, U := {Q ∈ ∂Q ηε : Q ∩ U = ∅}. By the fact that we have reduced to the case that E is an η ε -scale set in the sense of (3.7), we have F per ε (E, Q 3ηε ) ≥ cη 2 ε for all Q ∈ ∂Q ηε, U for some constant c > 0 depending on {c ξ } ξ∈V . Therefore, by (1.12) and (3.8) we obtain where C > 0 only depends on {c ξ } ξ∈V . This finishes the proof of (3.25).
We now have all the necessary ingredients to prove the first part of our compactness theorem. The second part on the validity of the Cauchy-Born rule for discrete symmetric gradients, i.e., (2.16), will be postponed to the next section. LetÊ ε ∈ A reg (U ) be the smooth set obtained from E ε through Lemma 3.5, applied for (U, U m ) in the place of (U, U ) here. This lemma guarantees that Q ε (E ε ) ∩ U ⊂Ê ε ⊂ U , and that we can find a decreasing sequence (ε m ) m∈N , ε m := ε(U m , U ) ∈ (0, 1) such that for every m ∈ N and every ε ∈ (0, ε m ],
By u ε we denote the piecewise affine interpolation of u ε with respect to a standard triangulation subordinate to the lattice εZ 3 , see also Subsection 3.2 in [44] for more details on the definition of the interpolation. By construction, it is clear that u ε | U\Êε ∈ H 1 (U \Ê ε ; R 3 ). By (1.10), (2.6), and (3.48)(iv) we find that Q ⊂ U \ Q ε (E ε ) for every Q ∈Q ε,Um\Êε . Moreover, for each such Q, we get due to [44, Lemma 3.6(i)]. Thus, by our structural assumptions on the discrete elastic energy (cf. Subsection 2.2), which are the same ones as in [44], by [ where the constant C > 0 is independent of m.

The validity of the Cauchy-Born rule for discrete symmetric gradients
In view of our compactness result [27, Proposition 3.1], in particular (3.50) above, the aim of this section is to show that, loosely speaking, weak L 2 -convergence of the symmetric gradients of the piecewise affine interpolations of (u ε ) ε>0 (outside the void set) implies the weak L 2 -convergence of the discrete symmetric gradients (ē(u ε )) ε>0 (outside a slightly modified void set) towards the same continuum limit. As already mentioned, since we believe that this result can be useful also in other discrete-to-continuum problems for elastic materials with surface discontinuities, we phrase the statement in any dimension.
For its formulation, we adopt all the notations introduced in the last subsection of the Introduction, as well as Subsection 2.1. Moreover, as before, for a discrete map u : εZ d → R d we denote by u its piecewise affine interpolation according to a triangulation subordinate to the lattice εZ d , constructed as in [44,Section 3.2]. The collection of d-dimensional simplices used in this construction will be denoted by S ε , and we also define where we recall the notation in and below (1.10). Moreover, note that by construction all dsimplices in S ε have the same L d -measure.
Then, there exist sets of finite perimeter where Z is defined in (1.6).
The second part of the compactness result is a direct consequence of Theorem 4.1.
Proof of Theorem 2.8 (ii). Having in mind (3.48) and (3.50), and the definition ω u :=ω u ∩ Ω, we consider the sets (W ε ) ε>0 and W given by 3) The first condition in (H.1) is clearly satisfied. The second one follows by the construction of the sets. Indeed, by min S 2 ϕ > 0, (3.50)(iv), (3.52), and the fact that sup ε>0 F ε (u ε , E ε ) < +∞, we get sup ε>0 H 2 (∂E * ε ∩ Ω) < +∞. Moreover, we have sup ε>0 H 2 (∂ * ωε u ) < +∞, see before ( The main point is that (H.3) implies (4.2)(iii), up to a modification of the void set which is negligible in the limit (see (4.2)(ii)). A property of this form has been established by Schmidt in [44, Theorem 2.6]: given discrete displacements u ε : εZ d → R d and corresponding piecewise affine interpolations u ε such that u ε ⇀ u weakly in H 1 (Ω; R d ) for some u ∈ H 1 (Ω; R d ), then it holds that ∇u ε ⇀ ∇uZ weakly in L 2 loc (Ω; R d×2 d ) . The proof particularly uses integration by parts for Sobolev functions and cannot be reproduced in our setting with voids. Therefore, a delicate blow-up argument and an approximation of GSBD 2 functions by Sobolev functions needs to be employed to reduce our problem to the setting in (4.4) for the corresponding symmetric gradients. The remainder of this section is devoted to the proof of Theorem 4.1. The reader is invited to skip the following subsections on first reading and to proceed directly with the lower bound in Section 5.  χ Ω\Vεē (u ε ) L 2 ( Ω) < +∞ . In particular, there exists ξ ∈ L 2 loc (Ω; R d×2 d ) such that, up to subsequence (not relabeled), χ Ω\Vεē (u ε ) ⇀ ξ weakly in L 2 loc (Ω; R d×2 d ) . (4.6) The core idea for the proof of the lemma is that discrete gradients can be related to the gradients of the corresponding piecewise affine interpolations, in the sense that there exists an absolute constant C * > 0 such that for all Q ∈Q ε it holds that see [44,Lemma 3.5(i)]. (Note that at this point we are dealing with discrete displacements defined on the whole lattice εZ d , and hence [44, Lemma 3.5(ii),(iii)] do not appear in our setting.) As this result needs to be applied on cubes of sidelength ε, we need to regularize the void sets W ε by introducing V ε , which are, loosely speaking, ε-cubic approximations of W ε . Then, (4.5) follows from (H.3), (4.7), and Korn's inequality on each cube. As Lemma 4.2 already shows (4.2)(i),(ii), it will remain to prove (4.2)(iii). To this end, in view of (4.6), we need to show that Since by (H.1) and (4.2)(ii) we have that χ Ω\Vε → χ Ω\W strongly in L 2 (Ω), it suffices to identify the weak limit ξ for L d -a.e. Lebesgue point x 0 ∈ Ω \ W of u. It is not restrictive to further suppose that the approximate gradient of u at x 0 exists. In fact, recall that for any u ∈ GSBD 2 (Ω), for We can now perform a standard blow-up procedure around such a point x 0 . To this end, we introduce some further notation. Given a sequence (ρ ε ) ε>0 ⊂ (0, +∞), and settingε := ε/ρ ε > 0, consider the rescaled discrete maps ζ ε :εZ d → R d , defined via where x 0,ε ∈ εZ d are chosen such that, for an absolute constant C > 0, there holds The piecewise affine interpolations of ζ ε (again according to [44,Section 3.2]) will be denoted by ζ ε . We mention again that we identifyē(u ε ) andē(ζ ε ) with their piecewise constant interpolations (being constant on each cell in the shifted lattice), i.e., for every i ∈εZ d (using again the notation i := i +ε( 1 2 , . . . , 1 2 )), we haveē(ζ ε )| Qε(î) :=ē(ζ ε )(i), where the fact that ε =ερ ε and (2.2) imply thatē As in the notations in the Introduction, we set V c ε := Ω \ V ε , and we define (Here and in the following, the superscript bl indicates rescaled sets and quantities in the blow-up procedure.) With this at hand, we have the following.
Properties (i)-(v) essentially follow by standard blow-up arguments for GSBD 2 -functions around points where the approximate gradient exists. For convenience of the reader, we will include the main arguments in the proof below. Eventually, property (vi) essentially follows from (i) and (4.11) by a compactness and rescaling argument.
The remaining step is the following lemma. (4.14) The key idea in the proof of the lemma is to approximate the functions ( z ε ) ε>0 by a sequence of Sobolev functions (v ε ) ε>0 also converging to ∇u(x 0 )y. This approximation is achieved by a deep result in the theory of GSBD 2 -functions, namely a Korn inequality for functions with small jump set. From this sequence, we define discrete maps w ε (later denoted by v η ε , but simplified here for convenience) by a sampling argument. We also introduce the corresponding piecewise affine interpolations w ε and we prove that w ε are close to z ε , e(w ε ) are close toē(ζ ε ), and w ε are bounded in H 1 (B 1 ; R d ) independently of ε > 0. In particular, ∇ w ε ⇀ ∇u(x 0 ) weakly in L 2 . We then conclude by applying (4.4) on the sequence of Sobolev functions ( w ε ) ε>0 which allows us to deduce ∇w ε → ∇u(x 0 )Z.
With these preliminary technical lemmata at hand, Theorem 4.1 follows.

4.2.
Proofs of technical lemmata. In this subsection we collect the proofs of the lemmata.
Step 1.(Construction) We start by introducing some extra notation. Recall the definition of the simplices in (4.1) and the comment below it. In particular, n d := #S ε,Q is a purely dimensional constant, Q = S∈Sε,Q S, and ∇ u ε | S is constant for all S ∈ S ε,Q and all Q ∈Q ε . Recalling also the notation introduced in (1.10), let us set where /d , and C d > 0 denotes the relative isoperimetric constant of the d-dimensional cube (see [3, (3.43)] for a version stated on balls instead of cubes). We now define (4.16) Clearly, by construction and the definition of β d > 0 we get (4.2)(i) since for each Q ∈Q bad,1 and for each S ∈ S ε,Q it holds that The key point of the decomposition (4.15) and the modification in (4.16) lies in the fact that and that (4.2)(ii) holds.
and the claim follows. Regarding the second assertion in (4.2)(ii), by construction we have where H d−1 (N ) = 0, see [40,Theorem 16.3]. Therefore, noting that for Q ∈Q bad,2 , due to [40,Theorem 16.3] (again up to a set of zero In view of (4.19), we can apply the the geometric Lemma 4.5 stated at the end of the section for s = 1 2 on P = Q \ W ε for Q ∈Q bad,2 . By this, and by (H.1) and (4.18) we conclude where we used that ∂ * (Q \ W ε ) ∩ int(Q) = ∂ * W ε ∩ int(Q) by construction. Thus, the second assertion in (4.2)(ii) follows.
Here, we used that the affine interpolation r Q coincides with r Q . This along with (H.3) and (4.21) implies (4.5). Eventually, (4.6) follows from weak compactness and (4.5), along with a suitable diagonal argument.
Proof of Lemma 4.3.
Step 1. We again let V c ε := Ω \ V ε . Choose a subsequence (not relabeled) such that (4.6) holds. Without loss of generality, we can assume that x 0 ∈ Ω∩W 0 (cf. [3,Definition 3.60] for the definition of density points), that x 0 is a Lebesgue point of both e(u) and ξ, and that the approximate gradient ∇u(x 0 ) exists. We claim that Since x 0 is a Lebesgue point of ξ, we get the first property in (4.22)(i). The second one follows analogously, by using (H.3) and (4.2)(ii) instead of (4.6), by x 0 ∈ Ω ∩ W 0 , as well as the fact that x 0 is a Lebesgue point of e(u). Moreover, (4.22)(iii) follows from the definition of the approximate gradient ∇u(x 0 ), see (4.8), which we assumed to exist at x 0 . For (4.22)(ii), we notice that (4.5) implies that for the family of Radon measures defined by µ bulk ε := χ V c ε |ē(u ε )| 2 dL d we have that sup ε∈(0,ε0) µ bulk ε ( Ω) < +∞. Hence, there exists a finite positive Radon measure µ bulk on Ω such that (up to a subsequence, not relabeled) µ bulk ε * ⇀ µ bulk in the sense of measures, as ε → 0. In particular, by the Radon-Nikodym theorem we can choose x 0 with the additional property that Moreover, since µ bulk is finite, except for a countable set of ρ ∈ (0, 1), we have µ bulk (∂B 2ρ (x 0 )) = 0. This along with weak convergence of measures shows the first part of (4.22)(ii). The proof of the second part is the same, by using V ε ⊃ W ε and (H.3) in place of (4.5). The proof of (4.22)(iv) is similar: by (4.2)(ii), up to a subsequence (not relabeled) there exists a finite positive Radon measure µ surf on Ω such that the measures H d−1 | ∂ * Vε * ⇀ µ surf on Ω in the sense of measures, as ε → 0. Notice now that for L d -a.e. x 0 we have that and we can thus assume without restriction that x 0 satisfies also (4.23). Indeed, suppose by contradiction that there exists a Borel set C ⊂ Ω with L d (C) > 0 and t > 0 such that Then, as a consequence of [3, Theorem 2.56], we would get µ surf (C) ≥ tH d−1 (C), implying that µ surf (C) = +∞. But this is a contradiction to µ surf being a finite measure. We now derive (4.22)(iv) from (4.23) by using that µ surf (∂B 2ρ (x 0 )) = 0 up to a countable number of values for ρ.
In particular, in view of (4.13)(iv), we have that Since we already know that z ε → ∇u(x 0 )y in measure on B 1 , we can infer from [13, Theorem 1.1, (1.5b)] that for a suitable subsequence we indeed have i.e., (4.13)(v) holds. In a similar fashion, we get . Therefore, to see (4.13)(vi), it suffices to check that ffl B1 F (y) dy = ξ(x 0 ). We can argue as follows. By (4.10) we get By this, along with the fact that ε/ρ ε → 0 as ε → 0, and (4.13)(ii), for ε > 0 small we can estimate This concludes the proof.
Proof of Lemma 4.4. Let x 0 and (ρ ε ) ε>0 be as in Lemma 4.3, and recall thatε := ε/ρ ε → 0 as ε → 0. In the following, we choose ε 0 > 0 such that B 2ρε 0 (x 0 ) ⊂⊂ Ω, and always suppose that ε ∈ (0, ε 0 ) without further notice. In view of the blow-up properties established in Lemma 4.3, we introduce the blow-up versions of our lattice, the cubic decomposition, and the d-dimensional simplices (cf. (4.1)) as Since (4.2)(i) is translation and scaling invariant, we also have that In order to localize our estimates, we fix from now on a small constant δ ∈ (0, 1), with δ > Cε, which will eventually be sent to 0. (Here C > 0 is a sufficiently large dimensional constant, so that all our subsequent estimates and set inclusions are true.) Step 1.(Application of the Korn-Poincaré inequality in GSBD 2 ) We claim that we can find an absolute constant C > 0, sets of finite perimeter (ω ε ) ε>0 ⊂ M(B 1 ), and Sobolev functions (v ε ) ε>0 ⊂ H 1 (B 1 ; R d ) such that v ε = z ε on B 1 \ ω ε , and Moreover, we can find sets of finite perimeter To see this, we first apply the Korn-Poincaré inequality for functions in GSBD 2 (B 1 ) with small jump set in the version of [12,Theorem 1.2]: there exists a dimensional constant C > 0 such that for every ε > 0 there exists a set of finite perimeter ω ε ⊂ B 1 with 29) and Now, by the classical Korn-Poincaré inequality, there exists an infinitesimal rigid motion r ε (y) = By (4.29) the first part of (4.27) follows. Now we claim that Assuming for the moment that (4.31) holds, this implies r ε H 1 (B1) ≤ C( e( z ε ) L 2 (B1) + 1) and thus, by the triangle inequality and (4.30), also the last inequality in (4.27) follows. Let us now show (4.31). We define and we note by (4.13)(iii),(iv) and (4.29) that By (4.30), the fact that v ε = z ε on B 1 \ ω ε , and (4.32) we can estimate i.e., r ε L 2 (B1\Kε) ≤ C( e( z ε ) L 2 (B1) + 1), where the last constant C > 0 depends additionally on ∇u(x 0 ). Next, note that if y ∈ B 1 \ K ε , then also −y ∈ B 1 \ K ε , and thereforê B1\Kε y dy = 0 .
Combining the last two observations, we get In view of (4.33), (4.31) follows. Next, let us construct the sets (T ε ε>0 . We define Let us first prove that In that respect, note that by construction, for every S ∈ S bl ε there exists M ε S ∈ R d×d such that . Therefore, for every S ∈ S bl ε \ S bl ωε such that S ∩ B 1−δ = ∅ and L d (S ∩ V c,bl ε ) > 0, we use the fact that δ > 0 is fixed large enough with respect toε (so that S ⊂ B 1 ) and (4.26) The last inequality, together with (4.27) and the fact that v ε = z ε on B 1 \ ω ε , yields (4.35). Moreover, by the definition of the sets (T ε ) ε>0 , it is clear that which concludes the proof of (4.28) and of this step.
Step 2.(Construction of a good sample for the discrete maps) In this step we define discrete mappings associated to the Sobolev functions v ε . To this end, apart from the parameter δ ∈ (0, 1) that we fixed in the beginning of the previous step (recall that δ ≫ε), we also fix η ∈ (0, 1 2 ). Eventually, in Step 5 of the proof, we will pass to the limits δ, η → 0. From now on, v ε shall denote a fixed representative of the corresponding L 1 -equivalence class. In this step, we show that there exists τ η ε ∈ Q ηε (0) such that 36) and that for the maps v η ε : where v η ε denotes the piecewise affine interpolation of v η ε , subordinate to the lattice Z bl ε . To this end, let us consider the set G η ε ⊂ Q ηε (0) defined as for which we claim that Indeed, consider the function g ε : Q ηε (0) → R, defined by On the one hand, since On the other hand, for η ∈ (0, 1 2 ), by exchanging summation and integration (note that the sum in the definition of the function g ε is finite), we havê Combining the last two inequalities, we arrive at (4.38). Now, let τ ∈ G η ε and S ∈ S bl ε , say S = conv{y 1 , . . . , y d+1 }, where (y l ) l=1,...,d+1 ⊂ Z bl ε . Let l ∈ {1, . . . , d}. By the theory of slicing of Sobolev functions (cf. [3,Proposition 3.105]), for L d -a.e. τ ∈ G η ε we have that v ε | [y l +τ,y l+1 +τ ] ∈ H 1 ([y l + τ, y l+1 + τ ]; R d ), and d dt v ε y l + τ + t(y l+1 − y l ) = ∇v ε y l + τ + t(y l+1 − y l ) (y l+1 − y l ) for a.e. t ∈ (0, 1) . (4.39) Note that |y l+1 − y l | ≤ Cε, that the unit vectors y l+1 −y l |y l+1 −y l | l=1,...,d form a basis of R d , and that |A| ≤ C for a purely dimensional constant C > 0, where A ∈ R d×d denotes here the corresponding transition matrix from the canonical basis {e 1 , . . . , e d } of R d to this basis. Let us also denote by v ε (· + τ ) the piecewise affine function with v ε (y l + τ ) := v ε (y l + τ ) for l = 1, . . . , d + 1, so that by construction, ∇ v ε | S (· + τ ) is constant on S ∈ S bl ε . Then, using Jensen's inequality and Fubini's Theorem, (4.39), and recalling the definition in (1.5), we can estimatê Summing over all simplices S ∈ S bl ε intersecting B 1−δ (and using that δ ≫ε), we infer Here, we used that for every S ∈ S bl ε there is only a bounded number of S ′ ∈ S bl ε , S ′ = S, such that (S)ε ∩ (S ′ )ε = ∅. Thus, the constant C > 0 in the last inequality is purely dimensional. Using the fact that L d (G η ε ) ≥ 1 2 (ηε) d , we can choose τ η ε ∈ G η ε such that This finishes the proof of (4.37). As by definition of G η ε also (4.36) holds, the step is concluded.
Step 5.(Weak convergence of χ V c,bl For the second integral, we can similarly use Hölder's inequality, (4.41), (4.47), (4.13)(vi), and the fact that the domain of integration is nonzero only for x ∈ B 1−2δ , to estimate i∈Lε\L good εˆQε (î)∩B 1−2δ Then, the desired property (4.14) follows by (4.13)(vi), by sending η → 0, the arbitrariness of δ, the boundedness of χ V c,bl ε e(ζ ε ) − e(u)(x 0 )Z in L 2 , and the fact that the set We finish this section with the proof a geometric fact that we used in the proof of Lemma 4.2.
Lemma 4.5. Let s ∈ (0, 1) and let Q ⊂ R d be a cube. Let P ⊂ Q be a set of finite perimeter with L d (P ) ≤ sL d (Q). Then, there exists a constant C > 0 only depending on s and d such that Proof. Without loss of generality (after possibly a rigid motion and a rescaling) we can assume that Q is the (open) unit cube, i.e., Q = (− 1 2 , 1 2 ) d . Thus, we suppose that L d (P ) ≤ s. Moreover, we can assume without restriction that for a small constant δ > 0 depending only on s, to be specified along the proof. Indeed, otherwise the desired inequality clearly holds for a constant depending on δ, as H d−1 (∂ * P ∩ ∂Q) ≤ 2d. We show the statement by induction over the dimension d. The statement for d = 2 can be deduced from [32,Lemma 4.6] as follows: suppose without restriction that P is connected (more precisely, indecomposable, see e.g. [3,Example 4.18]), as otherwise we argue separately for all indecomposable components. By [32,Lemma 4.6], there exists δ 0 ∈ (0, 1) sufficiently small such that, if H 1 (∂ * P ∩ Q) ≤ δ 0 , then for an absolute constant C 2 > 0, where diam(P ) := ess sup{|x − y| : x, y ∈ P }. (There, the proof has been performed for s = 1 2 , but it holds true for any s ∈ (0, 1), provided that δ 0 is chosen sufficiently small depending on s.) Since Q is a cube, we get H 1 (∂ * P ∩ ∂Q) ≤ Cdiam(P ), and since L 2 (P ) ≤ s by assumption, we deduce This concludes the proof for d = 2. Let us now assume for the inductive hypothesis that the statement holds true for some d ≥ 2, and let us prove it also for d + 1. We will only show that where x := (x ′ , x d+1 ) ∈ Q, since the surface measure in the two remaining faces of ∂Q can be controlled by repeating the argument after a rotation of the cube. We claim that L d+1 (P ) ≤ s implies that for s ′ := s+1 2 ∈ (0, 1), In fact, we define δ := δ(s) := s ′ − s = 1−s 2 , and suppose that (4.52) holds. Then, for a.e. − 1 ≤ δ , since otherwise by the coarea formula we find ) > δ , which contradicts our assumption (4.52). This observation along with L d+1 (P ) ≤ s implies for This shows (4.54). By slicing properties of sets of finite perimeter we have that for a.e. t ∈ (− 1 2 , 1 2 ) the set P t := P ∩Q d t has finite perimeter. Thus, by the inductive hypothesis, we can find a constant C > 0 depending on d and s ′ such that for a.e. t ∈ (− 1 2 , 1 2 ). By applying the coarea formula (cf. [40, (18.25)] for g := χ Q and slicing direction e d+1 ), denoting by ν the generalized normal to ∂ * P , and using (4.55), we can estimatê In a similar fashion, the coarea formula also implies that These two estimates along with (4.55) conclude the proof of (4.53).

Lower Bound
In this section we return to the case d = 3 and focus on the proof of the Γ-liminf inequality of Theorem 2.9. In particular, we aim at proving the following result.
We split the proof of this proposition into two separate lemmata, which prove the Γ-liminf inequality for the surface part and the elastic part of the energy separately. It is not restrictive to assume that lim inf ε→0 F ε (u ε , E ε ) < +∞, since otherwise there is nothing to prove. For the next lemmata, recall the definitions in (2.3) and (2.11)-(2.13).
Recalling the definition of the elastic energy in (2.12), and exhausting Ω with compactly supported open subsets Ω ⊂⊂ Ω, the statement follows.

Upper Bound
This final section is devoted to the proof of the Γ-limsup inequality.
Proposition 6.1. Let u = χ Ω\E u ∈ GSBD 2 (Ω) and E ∈ M(Ω) be a set of finite perimeter. Then, there exists (u ε , E ε ) ε>0 , E ε ⊂ Z ε (Ω) and u ε ∈ U ε (u 0 , ∂ D Ω, E ε ), such that (u ε , E ε ) d → (u, E) in the sense of Definition 2.7 and lim sup ε→0 F ε (u ε , E ε ) ≤ F 0 (u, E) . (6.1) Proof. We proceed in several steps. We first reduce to showing (6.1) for pairs of regular sets and displacement maps. Then, we replace these sets with polyhedral sets whose outer unit normal (at the points where this is defined) is a coordinate vector. Subsequently, we perform the approximation on the discrete level. The last step shows the convergence of the discrete energies to the respective continuum one.
Step 1.(Density of regular sets and displacements) By a density result, we can assume that ∂E∩Ω ∈ C ∞ and that u ∈ W 1,∞ (Ω \ E; R 3 ). Indeed, by [17,Theorem 2.2], for each E ∈ M(Ω) with H 2 (∂ * E ∩ Ω) < +∞ and each u = χ Ω\E u ∈ GSBD 2 (Ω), there exists a sequence of sets (E k ) k∈N with ∂E k ∩Ω ∈ C ∞ , χ E k → χ E in L 1 (Ω) and a sequence (u k ) k∈N with u k | Ω\E k ∈ W 1,∞ (Ω\E k ; R 3 ), u k | E k = 0, and tr(u k ) = tr(u 0 ) on ∂ D Ω \ E k , such that u k → u in measure on Ω and Strictly speaking, [17,Theorem 2.2] only ensures that ∂E k is Lipschitz and u k is in H 1 , but in the proof it is shown that ∂E k can be chosen of class C ∞ , see [17,Proposition 5.4], and u k can be approximated by smooth functions by standard approximation arguments. From now on, we therefore assume that ∂E ∩ Ω ∈ C ∞ and that u ∈ W 1,∞ (Ω \ E; R 3 ).
Step 2.(Reduction to polyhedral void sets with coordinate normals) By the regularity of E, the continuity of ϕ (cf. (2.10)), and a covering argument, we can reduce to the case of ∂E being polyhedral. Indeed, up to a localization argument, we can consider a sequence (E k ) k∈N ⊂⊂ E with polyhedral boundary and, for each k ∈ N, extend u to a function u k ∈ W 1,∞ (Ω \ E k ; R 3 ) such that convergence in (6.2) holds. As tr(u) = tr(u 0 ) on ∂ D Ω\E and E k ⊂⊂ Ω, this can be achieved in such a way that tr(u k ) = tr(u 0 ) on ∂ D Ω. Assuming now that E ⊂⊂ Ω has polyhedral boundary and tr(u) = tr(u 0 ) on ∂ D Ω, we prove that we can further restrict to the case of ν E ∈ {±e 1 , ±e 2 , ±e 3 }.
Step 4.(Energy convergence) The upper bound for the elastic energy, namely lim sup follows from the compatibility of W el surf and W el bulk (cf. assumption (vi) in the definition of the discrete elastic cell energy in Subsection 2.2) and the assumed a priori Lipschitz bound on u, see