Continuum Surface Energy from a Lattice Model

We investigate connections between the continuum and atomistic descriptions of deformable crystals, using certain interesting results from number theory. The energy of a deformed crystal is calculated in the context of a lattice model with general binary interactions in two dimensions. A new bond counting approach is used, which reduces the problem to the lattice point problem of number theory. The main contribution is an explicit formula for the surface energy density as a function of the deformation gradient and boundary normal. The result is valid for a large class of domains, including faceted (polygonal) shapes and regions with piecewise smooth boundaries.


Introduction
This article is concerned with the derivation of continuum surface energy from a standard lattice model, by exploiting results related to certain lattice point problems of number theory, e.g. [BL, BR, Hu, IKM, Pi].
We study the energy of a crystal, modelled as the part of a Bravais lattice L contained in a reference region Ω ⊂ R d , with atoms (elements of Ω ∩ L) interacting through a pair potential ϕ.
The potential may have unrestricted range but must decay fast enough. The crystal is subjected to a smooth deformation y : Ω → R d . The energy under consideration is E{Ω, y} = x∈Ω∩L z∈(Ω∩L)\x ϕ (|y(z) − y(x)|) To approach the continuum limit, one may scale the lattice, i.e., replace L by εL and rescale the potential to ϕ ε = ϕ( · ε ), then study asymptotics of the energy as ε → 0 [BBL, Mo]. Equivalently, one can rescale the region to rΩ and the deformation to y r = ry( · r ), with r = 1/ε, but leave L and ϕ unscaled. The emphasis of the present work is on the dependence of the energy on the geometry Rosakis of the boundary ∂Ω. Our main result is as follows. For the case d = 2, suppose Ω is a convex region with piecewise smooth boundary subject to certain restrictions (spelled out in Proposition 5.2) with outward unit normal n, and that the deformation is homogeneous, y(x) = F x, x ∈ Ω for some 1 F ∈ M 2×2 + . Then the energy satisfies as k → ∞, k ∈ Z. In the right hand side of (1.1), Ω(k) is a suitbaly defined region containing the same lattice points as the dilated region kΩ = {z : z = kx, x ∈ Ω}, namely Ω(k)∩L = kΩ∩L. Our main contribution is the following explicit formula for the surface energy density γ • : M 2×2 + × S 1 → R: Specifically, when Ω is a lattice polygon, Ω(k) is a certain rational polygon containing kΩ, with sides parallel to those of kΩ, and such that Ω(k) ∩ L = kΩ ∩ L. However, Ω(k) is not a dilation of Ω in general. In this case, the main contribution to the o(k) term in (1.1) is an O(1) corner energy that is obtained exactly.
On the other and, if Ω is a smooth C 2 strictly convex region, one may choose Ω(k) = kΩ; moreover the result then holds for any real (not only integer) sequence k → ∞; see Proposition 5.2.
Our results show that for a large class of regions Ω, if one replaces Ω(k) by kΩ, (1.1) does not hold, unless γ • is replaced by another functionγ (see (1.7) below), whose dependence on the normal n involves a dense set of discontinuities. The hypotheses of the standard surface energy minimization theorem (for the two-dimensional case see Dacorogna and Pfister [DP] and Fonseca [Fo] for the three dimensional version) yielding the Wulff shape may not be fulfilled in general. In contrast, γ • (F, n) from (1.2) is Lipschitz in n.
The above formula for the surface energy density is analogous to the well-known Cauchy-Born formula for the stored energy function in the first term of (1.1): In more than one dimension, the first rigorous derivation of continuum energy functions from atomistic models is due to Blanc, Le Bris and Lions [BBL], who study (among other problems) the asymptotics of the energy 2 of a crystal Ω ∩ εL, subject to a prescribed smooth deformation y : Ω → R 3 as ε → 0. The dominant term is the usual elastic energy Ω W (∇y(x))dx with W 1 M 2×2 + is the set of 2 × 2 matrices with positive determinant. 2 This energy is divided by #(Ω ∩ εL) and has rescaled potential ϕ ε .
given by (1.3). The next term, of order ε in Theorem 3 of [BBL], is a surface integral over ∂Ω that involves values of the deformation gradient and unit normal. The form of this surface energy is not explicit and it is not clear to what extent it can be expressed as a function of those variables. Terms of order ε 2 include a volume integral of an explicitly determined higher gradient energy, but also surface terms; the latter are left unspecified.
As shown in one dimension by Mora-Corral [Mo], the higher order terms in the asymptotic expansion of the energy in powers of ε depend on the choice of the sequence of ε → 0. In Theorem 3 of [BBL], this choice is restricted by the hypothesis that there exist a sequence ε = ε k → 0 as k → ∞, such that #(Ω ∩ ε k L) = |Ω|/ε d k (in dimension d). Letting r = 1/ε and scaling Ω instead of the lattice, this means that for some sequence r k → ∞, #(r k Ω ∩ L) = |r k Ω|. (1.4) In the present work we rely on bond counting arguments instead of asymptotics to a large extent.
A byproduct of this approach is an explanation of this sequential dependence issue. Our method hinges on finding, for each w ∈ L, the w-bond number of Ω: For large r, the dominant contribution to N w (rΩ) is #(rΩ ∩ L). Finding the asymptotics of this number as r → ∞ is the lattice point problem of number theory [BR, IKM, Ts]. This reduces to studying the asymptotics of the lattice point remainder R(r) = #(rΩ ∩ L) − |rΩ|, the difference of the two sides of (1.4). In two dimensions, the problem is open for general domains with piecewise C 1 boundary, while even the Gauss circle problem (Ω the unit disk, L = Z 2 ) is not completely settled [Hu]. Through (1.5), the lattice point remainder enters our estimates for the energy E{rΩ, y r }, whose asymptotic form thus depends on the sequence r k through R(r k ). This can be problematic as R is discontinuous and highly oscillatory. In general, the behavior of R(r) depends strongly on the shape of ∂Ω.
For Ω a lattice polygon (one whose vertices are lattice points), R(r) is of same order as the surface energy-R(r) = O(r d−1 ) in dimension d-and can be characterized explicitly; see, e.g., Lemma 2.2 below. For smooth convex domains, as shown by van der Corput [Co], R(r) = O(r 2/3 ), between the orders of the surface and the gradient energy of [BBL], but difficult to characterize [Hu]. Hypothesis (1.4) made by [BBL] is equivalent to existence of a sequence r k such that R(r k ) = 0, thus it eliminates certain undesirable higher order terms from a Riemann sum of the elastic energy. Unfortunately however, it is not known for which choices of domain Ω such a sequence exists. We adopt an alternative approach that avoid making such a hypothesis.
Crystals typically occur in faceted form in their natural state (for instance, the Wulff shape, e.g., [He, DP, Fo]). This is because of surface energetics affecting crystal growth, but also because cleavage fracture creates new surfaces along special crystallographic planes. This means that they can be modelled as crystallographic polyhedra, whose facets inhabit crystallographic planes (that contain a two-dimensional sublattice of L).
In Section 2 we assume that Ω is a lattice polytope, i.e, one whose vertices are lattice points.
This does not sacrifice too much generality over crystallographic polytopes. Indeed, if Ω is crystallographic polytope, then kΩ is a lattice polytope for some k ∈ Z. In addition, there is a lattice If Ω is convex, then Ω ′ = conv{Ω ∩ L}. In view of Theorem 3 of [BBL], one expects that the dominant surface energy term does not involve higher gradients of the deformation. Accordingly, it suffices to assume that the deformation is homogeneous (affine). To keep the geometry simple, we confine our analysis to two dimensions. Unlike [BBL, Mo], initially we do not employ a limit process, but rather a bond counting technique. The computation of the energy is reduced to that of (1.5). We then show that this calculation reduces to a number of lattice point problems. The solution of the latter for lattice polygons is furnished by Pick's Theorem [Pi]. The lattice point remainder R(k) is known exactly (for k ∈ Z) and contributes to the surface energy explicitly, being of the same order.
In Section 3 we compute the energy of polygonal crystals. For an interatomic potential of finite but arbitrary range, we obtain the energy of essentially any convex lattice polygon exactly (Proposition 3.1). This result is not asymptotic and does not suffer from the sequential dependence issue explored in [Mo]. Let the deformation be y(x) = F x, x ∈ Ω. The energy equals the exact sum of the elastic energy Ω W (F )dx plus the surface energy ∂Ω γ ⋄ (F,n)dx, plus the corner energy N i=1 τ (F, n i , n i−1 ), summed over the N vertices of Ω. The surface energy density is explicitly obtained: wheren is a normal to ∂Ω whose components on each facet are the Miller indices (irreducible integers) of the corresponding lattice plane, and ϕ is the interatomic potential. The corner energy τ (F, n i , n i−1 ) is also explicit but more complicated; apart from F , it depends on the two unit normals of the facets meeting at the ith vertex.
For an infinite range potential this result retains only asymptotic validity for a lattice polygon kΩ as k → ∞; the three energies just mentioned are the first three terms of the asymptotic expansion of the energy for large k (Proposition 3.4).
In Section 4, we consider regions with smooth boundaries. Because of its construction based on lattice polygons, the surface energy density (1.6) is only defined for "rational" directions of the surface normal; n = (ν 1 , ν 2 ) ∈ S 1 is called rational if ν 2 /ν 1 is a rational number or ν 1 = 0, irrational otherwise. It is natural to ask how (1.6) can be extended to irrational normals. When Ω is strictly convex and ∂Ω is smooth for example, the normal is irrational almost everywhere on ∂Ω. We start by letting ∂Ω be of class C 2 with positive curvature. The key observation is that the convex hull of all lattice points contained in such an Ω is a lattice polygon. This allows us to use number-theoretic results on the asymptotic properties of such hulls due to Bárány and Larman [BL]; see also the survey [IKM]. Perhaps surprisingly, the surface energy density for smooth strictly convex regions (Proposition 4.1) is different from (1.6). It is given by (1.2), where n is the unit normal to ∂Ω and can take on irrational values. The difference is due to the lattice point remainder [Co, Hu], which is now of lower order than the surface energy. As a result, the asympotic expression for the energy of inflated regions rΩ is sequence-independent; the sequence of r → ∞ is not restricted to be integer but arbitrary.
We then consider more general regions with piecewise C 1 boundary that comprises flat facets as well as curves with positive curvature. For such regions, the surface energy density function, now defined for all n ∈ S 1 , is obtained in Proposition 4.4: with γ ⋄ from (1.6) and γ • from (1.2). The dependence of the surface energy density on the normal is rather pathological. Specifically,γ(F, ·) : S 1 → R is continuous at irrational n, discontinuous at rational n, and almost nowhere differentiable (Proposition 4.6). Because of this, the surface energy density need not satisfy the usual hypotheses of the Wulff theorem (determining the domain that minimizes the surface energy under fixed measure); see e.g. [DP, Fo], but also Remark 5.1.
In Section 5, we resolve the difficulties due to discontinuous dependence of the surface energy on the unit normal. This dependence is due to the behaviour of the lattice point remainder of regions with rational boundary normal. We then alter the region Ω so as to change its measure, but not the lattice points it contains. The goal is that the lattice point remainder of the modified region should be of lower order than the surface energy. For example, if Ω, hence kΩ, is a lattice polygon, translate each side of kΩ outwards by half the distance to the next crystallographic plane with the same normal. This results in a rational polygon Ω(k) that contains the same lattice points as kΩ.
Note, however, that Ω(k) is not a rescaling of Ω in general. The lattice point remainder of Ω(k) is O(1) as k → ∞, of lower order than the surface energy. This allows us to write the latter in the form ∂Ω(k) γ • (F, n)ds. The associated surface energy density γ • , given by (1.2), is Lipschitz continuous in the unit normal. This and additional considerations discussed in Section 5, show that γ • is the appropriate density for the determination of the Wulff shape that minimizes the surface energy ∂Ω γ • (F, n)ds over a suitable class of regions Ω with fixed measure [He, Fo].
A more realistic approach to surface energy would allow for "relaxation" of atomic positions from the macroscopic deformation near the boundary. Such deviations might be determined by minimization of the atomistic energy. This is a formidable problem in the present setting (more than one dimension, general boundary geometry, arbitrary interaction range, nonconvex potentials).
One of the few results in this direction is due to Braides and Cicalese [BC]; they obtain the relaxed surface energy in one dimension using Γ-convergence. The result is not explicit and seems difficult to compare quantitatively with the explicit "constrained" energy of Mora-Corral [Mo]. In two dimensions, Theil [Th] calculates the relaxed surface energy of a crystal with quadratic short range potentials; the result is in the form of a perturbation of the constrained surface energy.
In order to obtain quantitative information on the difference between the relaxed and constrained surface energies, numerical optimization of the atomistic energy was recently performed for a completely unconstrained, Lennard-Jones two-dimensional crystal [Ro]. Atomic positions were allowed to relax from initial positions forming a lattice triangle or hexagon with low Miller-index boundary. The constrained energy was obtained by minimizing over the deformation gradient matrix of a homogeneous deformation that the atoms are constrained to follow. It was found that the difference between the relaxed and constrained surface energies is typically less than three percent (after the appropriate scaling and bulk energy is accounted for). This suggests that in some situations the relaxed and constrained surface energies may be quite close. In analogous one-dimentional computations, the results agree qualitatively with the conclusions of [BC], while the difference between the relaxed and constrained surface energies is less than one percent. Values of this difference computed in three dmensions using density-functional theory for low Miller-index surfaces in various metals are usually less than three percent; see, e.g., [EHF].
Many of the results presented here, in particular expressions (1.6) through (1.7) for the surface energy density, are valid for three-dimensional crystals as well [Ro].

The Bond Counting Approach
For subsets P , Q of R n , define the Minkowski sum P ⊕ Q = {p + q : p ∈ P, q ∈ Q} and write The lattice is L = Z 2 unless otherwise noted.
Remark 2.1. All of our results can be immeditately adapted to any Bravais Lattice L * by incorporating the linear mapping from L onto L * into the deformation. Then expressions like (1.2) remail valid if L is replaced by L * , provided the the linear mapping from L onto L * has unit Jacobian determinant.
We assume that the reference region Ω ⊂ R 2 is a convex body, or a compact convex set with is the set of all w-bonds of Ω (bonds with bond vector w). We will use the abbreviation The energy of the homogeneous deformation y(x) = F x can be written as The factor of 1/2 occurs since b(x, w) = b(x + w, −w) and the potential ϕ is even in w. Interchanging the order of summation above we obtain Evidently, in order to determine the energy, it suffices to calculate, for each w ∈ L, the w-bond number of Ω, i.e., N w (Ω) = #B w (Ω); see (2.1): Clearly the number of w-bonds "starting" in Ω equals the number of lattice points of Ω: Some of these bonds are not contained in B w (Ω): (2.5) so that S + w is the part of ∂Ω through which w points outwards. Denote by T † w (Ω) the set of all w-bonds that intersect S + w and terminate outside Ω.
(2.6) Some of these bonds "straddle" Ω, that is, have both endpoints outside Ω but intersect ∂Ω; specifically, Then obviously in view of (2.4), As a result, Roughly speaking, the number of w-bonds in Ω equals the number of lattice points in it, minus the number of bonds that traverse the boundary at least once, plus the number of bonds that traverse the boundary twice. The reason for the splitting (2.8) is that each term can be evaluated using results from geometric number theory.
One important case we will consider is when Ω ⊂ R 2 is a convex lattice polygon. In particular, The number of lattice points in Ω, #(Ω ∩ L), is addressed by Pick's Theorem, [Pi, Re, BR], a variant of which is the following Lemma 2.2. Let Ω be a simple closed lattice polygon with facets S i and outward Miller normal Equivalently, letting θ i be the (dihedral) angle between normals of facets meeting at the ith vertex, (2.10) Proof. Pick's Theorem [Pi, Re] states that (since Ω is closed). If two neighbouring lattice points in a facet S i differ bym i ∈ L (with relatively contains both its endpoints. Now the Miller normaln i =m ⊥ i , so that |n i | = |m i | and (2.9) follows. Also, (2.10) is a trivial consequence of (2.9), given that the sum in (2.10) equals 1.
Remark 2.3. The shape of naturally occurring crystals is very often faceted (polyhedral). Thus one might start by assuming that Ω is a polygon, though not necessarily a lattice polygon. In that case If Ω is a crystallographic polygon, so that its facets are contained in crystallographic lines, then its vertices need not be lattice points. However, one can then show that there is some integer k such that kΩ is a lattice polygon.
Remark 2.4. Eq. (2.10) has an interesting interpretation. It exactly equates a discrete quantity (number of atoms in Ω) with a continuum expression: the "volume" integral of a bulk density, plus the "surface" integral of a surface density, plus contributions of corners. We will show in the sequel that both the w-bond number N w (Ω) and the energy admit analogous representations.
The first term in (2.8) is given by (2.9). Turning to the second term, let P i (w) be the parallelogram b 0 ⊕ S i with two parallel sides S i and w + S i if w · n i > 0, P i (w) = ∅ otherwise. Then it is easy to (2.14) In general, P (w) is not convex. However, if one defines then Ω w is a convex lattice polygon, being the Minkowski sum of two such sets. In fact, Also P (w) = Ω w \ Ω, while Ω ⊂ Ω w . This and (2.14) imply The right hand side can be evaluated using Lemma 2.2 for each term. Note that ∂Ω w comprises ∂Ω \ S + w , w + S + w and two w-bonds joining these two pieces. In view of (2.12) the result is where x = (x + |x|)/2 for x ∈ R and n = n i on S i is the unit outward normal to ∂Ω. Here |b 0 |/|w| = gcd(w).
We will show next that for |w| small enough compared to the facets of Ω, Q(w) consists of one or two triangles, each having a vertex at one of the two ends of the simple polygonal line S + w . For example, if Ω = [0, 3] 2 and w = (1, 1), Q(w) consists of the triangle with vertices (0, 3), (1, 4) and (1, 3) and its image under reflection about the (1, 1)-axis. Any b ∈ T ‡ w (Ω) intersects two different facets of ∂Ω by (2.7). Let where v i ∈ Z 2 are the vertices of Ω. The shortest line segment with endpoints on non-adjacent facets has length δ. If |w| < δ, b ∈ T ‡ w (Ω) necessarily intersects two adjacent facets, say S i and S i−1 meeting at some vertex v i , with outward normals n i , n i−1 (where n 0 = n N ). Since both endpoints of b are outside Ω, w ·n i and w ·n i−1 must have opposite signs. Then in case w ·n i > 0 and w ·n i−1 < 0, x + w is in the triangle with vertices v i , v i + w and the intersection of S i and w + S i−1 , which is therefore part of Q(w). If the reverse inequality holds, the triangle with vertices v i , v i + w and the intersection of w + S i and S i−1 is part of Q(w). Regarding lattice point count, both cases reduce to the triangle with base b 0 and sides normal to n i and n i−1 : In addition, the relative interior of the base b(v i , w) of the triangle with endpoints v i , v i + w is also part of Q(w) and contains gcd(w) − 1 lattice points. Consequently, if |w| < δ, ( 2.22) Unfortunately, T (w, n i , n i−1 ) is not a lattice polygon in general, since q need not have integer coordinates and Lemma 2.2 does not apply. Instead, we count the lattice points inside the triangle more directly: Lemma 2.5. Suppose (w · n i )(w · n i−1 ) < 0 and let T = T (w, n i , n i−1 ) ⊂ R 2 be the triangle of (2.21). Let u ∈ Z 2 be such that {u,w} is a lattice basis for Z 2 . Then where for w ∈ Z 2 and unit n, m ∈ R 2 with (w · n)(w · m) < 0, Proof. Let n = n i , m = n i−1 . Since in (2.21) q · n = 0, q = λn ⊥ for some λ ∈ R. Then solving (q − w) · m = 0 for λ gives q as in the second of (2.23). Letw = w/ gcd(w) = (w 1 ,w 2 ) and suppose u = (u 1 , u 2 ) ∈ Z 2 solves u ·w ⊥ = 1, orw 2 u 1 −w 1 u 2 = 1. This is solvable by Bezout's Lemma since gcd(w 1 ,w 2 ) = 1. Then the matrix A = col(u,w) has unit determinant u ·w ⊥ = 1 and integer entries, hence so does A −1 = row(w ⊥ , u ⊥ ). As a result {u,w} is a lattice basis for Z 2 , while the linear transformation with matrix A −1 is lattice invariant . Now T ′ = A −1 T has vertices 0, (0, k) ∈ Z 2 and p = (α, β), where k = gcd(w), (α, β) = (q ·w ⊥ , q · u ⊥ ); (2.25) in general p is not a lattice point. Suppose for the moment that α > 0. Then For x ∈ R let ⌊x⌋ ′ be the greatest integer strictly less than x and ⌈x⌉ ′ the least integer strictly greater than x. Then the number of lattice points on a segment {(x 1 , x 2 ) : x 1 = j, µ < x 2 < ν}, where j ∈ Z and µ < ν ∈ R equals ⌊ν⌋ ′ − ⌈µ⌉ ′ + 1. Hence, Since ⌊x⌋ ′ = ⌈x⌉ − 1 and ⌈x⌉ ′ = ⌊x⌋ + 1, the above reduces to s(α, β, k) in (2.24). It then follows from (2.25) and (2.23) that #(T ′ ∩ L) = N T (w, n, m). The linear transformation with matrix A is lattice invariant and thus #(AT ′ ∩ L) = #(T ′ ∩ L) [BP], while AT ′ = T . In case α < 0, reflect T ′ by replacing α by |α|. If α = 0 thenT ′ = ∅.
This together with (2.22) gives (2.26) To obtain an expression for the w-bond number of Ω, merely substitute (2.9), (2.18) and (2.26) into (2.8) and rearrange. Observe that for a given bond vector w, N w (Ω) is completely determined by the area |Ω|, the lengths |S i | of the facets, and their orientations through the Miller normalsn i : Lemma 2.6. Suppose |w| < δ, cf. (2.20). Then the w-bond number of Ω is given by (2.27) Remark 2.7. The above can readily be written in a form similar to (2.10): as a bulk integral, plus a "surface" integral, plus corner contributions, for suitable normal-dependent densities g and h. See also Remark 2.4.
The present approach of counting bonds has certain similarities with the bond density lemma of Shapeev [Sh].

Surface Energy of Lattice Polygons
We are now in a position to compute the energy. Consider first a finite-range potential that only involves bonds within a bounded set. Let the bond range R ⊂ L \ {0} be symmetric, so that w ∈ R =⇒ −w ∈ R. Allow the interatomic potential ϕ w (·) to depend explicitly on w, require and define the energy of the homogeneous deformation y(x) = F x, x ∈ Ω, where ϕ w : (0, ∞) → R is not restricted to be regular in any way.
Proposition 3.1. For F ∈ M 2×2 + ,m ∈ Z 2 and unit n, m ∈ R 2 define the stored energy function the surface energy density function and the vertex energy function

5)
where the sector step function and θ(n, m) is the angle between n and m, while N T is defined in Lemma 2.5. Suppose the bond range R is bounded with max w∈R |w| < δ, cf. (2.20). Letn =n i on S i . Then the following expression is exact: Proof. As in the argument leading to (2.3), one can write (3.2) as By the hypothesis on R, Lemma 2.6 holds for all w ∈ R. Multiply (2.27) by ϕ w (|F w|) and sum the result over w ∈ R. Interchange the order of summations, noting that w∈R w ·n i ϕ w (|F w|) = 1 2 w∈R |w ·n i |ϕ w (|F w|) by the symmetry of R and the first of (3.1), also that the sum of the (dihedral) angles between normals of facets meeting at vertices N i=1 θ(n i , n i−1 ) = 1, and finally that summation over w in the sector of R where (w · n)(w · m) < 0 can be replaced by summation over R provided the summand is multiplied by H n,m (w).
Remark 3.2. The above result is not asymptotic but exact, since we have made no use of asymptotics so far. It applies to convex lattice polygons that are arbitrary apart from the restriction that the bond range is smaller than the characteristic size δ of (2.20).
Next we consider infinite-range potentials, where R = L \ {0}. We seek the energy of the kth dilation kΩ of the region Ω, k ∈ Z + . Here we have no choice but to let k be an integer; otherwise kΩ is not a lattice polygon in general. The following will be useful.
For convenience we suppose that the interatomic potential ϕ w (·) = ϕ(·) (does not explicitly depend on w), although this is not essential.
Proposition 3.4. Suppose the interatomic potential ϕ : (0, ∞) → R satisfies the following: for each r 0 > 0 and for some constants C = C(r 0 ) and d > 2, |ϕ(r)| < Cr −(2+d) for r ∈ [r 0 , ∞). (3.8) Proof. Note that δ(kΩ) = kδ(Ω) = kδ in (2.20), so that Lemma 2.6 for kΩ holds provided Split the energy as follows: (3.10) Now it is clear that for any w ∈ L and k ∈ Z + , for some constant C > 0, since all bonds within kΩ start in kΩ and by Lemma 2.2 applied to kΩ (the dominant term in (2.9) would be |kΩ| = k 2 Ω). This provides a bound for the second term in (3.10): where we invoked (3.7), α > 0 is such that |F z| > α|z| for all z ∈ R 2 and we used Lemma 3.3 with ρ = kδ and p = d; and C is a generic constant with possibly different values each time it appears.
The first term in (3.10) is covered by Proposition 3.1 applied to kΩ, since w ∈ R k means |w| < kδ = δ(kΩ). Noting that |kΩ| = k 2 |Ω|, |kS i | = k|S i |, Proposition 3.1 implies where W k , γ k and τ k are given by (3.3), (3.4) and (3.5) with R k in place of R; see (3.9). Recalling that W , γ ⋄ and τ are defined by the same equations with R = L \ {0}, using Lemma 3.3 with M = 2, we may estimate (omitting arguments) We only demonstrate the third of these, the others being easier. Recall that in (3.5), N T is the number of lattice points in the interior of a certain triangle T whose area is bounded above by C|w| 2 , cf. Lemma 2.5. By Pick's Theorem (2.11) (applied to the lattice parallelogram of smallest area A containing T , and having the same base) the area A exceeds N T hence N T < C|w| 2 . Also gcd(w) ≤ |w|, |H n,m | ≤ 1, hence we have from (3.5), C|w| 2 |ϕ(|F w|)| < C w∈Z 2 , |w|>kδ |w| −d < Ck 2−d proceeding as in (3.11). By (3.13), replacing W , γ ⋄ and τ by W k , γ k and τ k in (3.12) produces an error of O(k 2−d ). Combine this with (3.11) and (3.10) to obtain (3.8).

Surface Energy for More General Boundaries
We examine the surface energy density function γ ⋄ in (3.4) more closely, paying attention to its dependence on the surface normal. Due to its construction, γ ⋄ (F, ·) :M → R is defined only for "rational directions", that is, on the set of Miller normals M = {n :n = (ν 1 , ν 2 ) ∈ Z 2 , gcd(ν 1 , ν 2 ) = 1}. The first term (involving the sum) reduces to a function of the unit normal n, and trivially admits a unique continuous extension onto the whole of the unit circle S 1 . There is no such extension for the second term. Define the rational and irrational direction sets as respectively, whereM is defined in (4.1). Thus a vector is rational (irrational) if the tangent of the angle it makes with the usual basis vectors is rational (irrational). Since facets of lattice polygons have rational normals, the surface energy density γ ⋄ is defined only for such directions. Note that for each n ∈ S 1 R there is a uniquen =n(n) ∈M withn/|n| = n. The question arises as to how one can extend the definition ofγ ⋄ (F, n) = γ ⋄ (F, |n(n)|n), n ∈ S 1 R , to the whole of S 1 . This is related to another question: what is the surface energy when ∂Ω is smooth, for example ∂Ω = S 1 ? It turns out that this question can be answered, at least partially, using the present approach. The basic idea is that even if ∂Ω is not polygonal, but smooth, the convex hull of all lattice points inside Ω is a convex lattice polygon.

Proposition 4.1.
Let Ω ⊂ R 2 be strictly convex and ∂Ω be C 2 with positive curvature. Suppose ϕ is as in Proposition 3.4, but with d > 3. Define the reduced surface energy density γ Then for any sequence r = r k → ∞ as k → ∞ (r k ∈ R + , k ∈ Z + ), where n : ∂Ω → S 1 is the unit outward normal to ∂Ω.

Remark 4.2.
This asympotic result for inflated regions rΩ, is sequence-independent; that is, the sequence of r → ∞ is not restricted to be integer but is arbitrary. This occurs because the lattice point remainder R(r) = O(r 2/3 ) [Co, Hu] is of lower order than the surface energy. In contrast, the surface energy for lattice polygons, or the more general regions considered in Proposition 4.4 depends on the sequence of dilation factors. The dependence on the dilation sequence is thoroughly studied in one dimension in [Mo].
Proof. For each r > 0, let Ω r = conv(rΩ ∩ L). Then Ω r ⊂ rΩ is a convex lattice polygon, while rΩ ∩ L = Ω r ∩ L. Hence, in view of (2.2), where y(x) = F x for x ∈ rΩ. The calculation of E{Ω r , y} proceeds as above with one exception.
Since Ω r ⊂ rΩ, and rΩ \ Ω r contains no lattice points, (2.19) and (2.16) imply Q w (Ω r ) ⊂ Q w (rΩ). (4.7) Let q, q ′ ⊂ ∂(rΩ) be the two points of ∂(rΩ) where the tangent vector is w, and B rρ ⊂ rΩ be a disk with ∂B rρ tangent to ∂(rΩ) at q, where ρ is the smallest radius of curvature of ∂Ω. Also let B ′ rρ ⊂ rΩ be a similar disk tangent to ∂(rΩ) at q ′ . Then for r large enough it is easy to see that The connected component ofQ w (B rρ ) containing q is contained inside an isosceles triangle with base a w-bond (with length |w|), and height the distance from the base middle to the intersection of the two circles ∂B rρ and w + ∂B rρ ; these are tangent to the base at its endpoints. The triangle height is thus bounded by C(r)|w|, where C(r) approaches zero for large r. A crude but sufficient upper bound of the lattice point count of this set, hence also of the right hand of (4.8), is C|w| 2 , with C independent of r. In view of of (4.7), #T ‡ w (Ω r ) is also bounded by C|w| 2 . This estimate replaces the sum over vertices (second sum) in (2.27). Since w∈L\{0} |w| p ϕ(|F w|) are absolutely convergent for p = 0, 1, 2 as one infers from Lemma 3.3, it follows that Here we have used (4.2) and (4.4), then (2.10), in which the last term (sum) equals 1, together with the fact rΩ ∩ L = Ω r ∩ L. We turn to ∂Ωr γ • (F, n)ds. Recalling (4.4), a typical term involves ∂Ωr |w · n|ds = 2|Proj w ⊥ ∂Ω r ||w|, (4.10) |Proj w ⊥ ∂Ω r | being the length of the projection of ∂Ω r onto a line perpendicular to w. This follows after splitting ∂Ω r into two pieces, over which w · n is ≥ 0 and ≤ 0, and using the Divergence Theorem on each. Next, we show that (4.11) where the constant C is independent of r > 1 and w. There are lattice points z − and z + ∈ ∂Ω r ∩ L, such that Ω r lies entirely between lattice lines l − , l + with normal w ⊥ and containing z − , z + , respectively. Consider the part of ∂(rΩ) that lies outside the strip bounded by l + and l − . It consists of two disjoint arcs, one to the "right" of l + and the other to the "left" of l − . The length of the projections of these two arcs onto the w ⊥ axis equals the difference in (4.11). Let c + be the arc to the right of l + (with endpoints in l + ). Let s be the region bounded by c + and l + . The only lattice points it contains are in l + . This is true since rΩ \ Ω r is free of lattice points. By the strict convexity of rΩ, there is a unique q ∈ c + where the normal to c + is w ⊥ . Consider the osculating circle of c + at q. Let s ′ ⊂ s be the portion of the osculating disc contained in s; it is a circular segment whose height (in the direction w ⊥ ) equals the thickness of s (the length of its projection onto a line along w ⊥ ). The radius of the circle is rρ for some ρ > 0. There are two possibilities. Either s ′ lies between l + and the next lattice line l ′ with normal w ⊥ to the right of l + , or it extends beyond l ′ to the right.
In the first case the height of the segment s ′ is 1/|w|, the distance between adjacent lattice lines with normal w ⊥ . In the second case, let s ′′ be the portion of s ′ to the right of l ′ . Then s ′′ is also a circular segment and free of lattice points. Suppose its chord length is c and height is h. Since the radius of the circular arc is rρ, we have h 2 − 2rρh + c 2 /4 = 0. Solving this for h/(rρ) and using the inequality 1 − √ 1 − x < x for 0 < x < 1 yields h < c 2 /(4rρ). Now since the circular segment s ′′ is free of lattice points and its chord is in l ′ , the chord length c < |w| ≤ |w| (since the distance between adjacent lattice points in l ′ is |w|.) Hence h < |w| 2 /(4rρ). The total height of the larger circular segment s ′ is h + 1/|w| which is thus bounded by C|w| 2 for r ≥ 1. The thickness of s in the direction normal to w is the same as this height. This shows (4.11).
According to Proposition 4.1, when ∂Ω is smooth and strictly convex, so that the normal vector is irrational almost everywhere on ∂Ω, the surface energy density is given by (4.4); in contrast, for lattice polygons (with rational normal a.e. on ∂Ω), the surface energy density is given by (4.2). This suggests that we combine the two expressions in defining a surface energy density for all values of the unit normal. That will allow us to treat a more general case with Ω a (not necessarily strictly) convex body. We do place some restrictions on ∂Ω: flat parts of ∂Ω must be lattice segments (with rational normals). Corners have to be lattice points.
We now state the main result of this section: Proposition 4.4. Assume that Ω is a convex body with ∂Ω Lipschitz, and that there is a finite set is a C 2 curve and one of the following two alternatives holds: Suppose ϕ is as in Proposition 3.4, but with d > 3. Define the extended surface energy densitŷ γ(F, ·) : S 1 → R as follows: with γ • defined in (4.4) and S 1 R , S 1 I defined in (4.3). Then as k → ∞, k ∈ Z + , where n : ∂Ω → S 1 is the unit outward normal to ∂Ω.
Proof. We now choose r = k ∈ Z + and let Ω k = conv(kΩ∩L). The part of the proof of Proposition 4.1 prior to (4.9) is easily adapted to the present setting, so that once again, as k → ∞, with γ ⋄ as in (4.2), Let ∂Ω f be the union of those S i that are straight segments and ∂Ω c the union of the S i with positive curvature, so that ∂Ω = ∂Ω f ∩ ∂Ω c . By hypothesis, for k ∈ Z + we have kv i ∈ ∂(kΩ) ∩ L, hence Our hypotheses regarding ∂Ω c , specifically alternative (i), ensure that n ∈ S 1 I a.e. on k∂Ω c , while (ii) implies that n ∈ S 1 R a.e. on k∂Ω f . Using (4.15), rewrite the above as It remains to show thatR(k) = O(k 2/3 ) as k → ∞, k ∈ Z + . Let i ∈ J c , so that S i satisfies alternative (i) in the statement of Proposition 4.4. Let S i k be the portion of ∂Ω c k between kv i and kv i+1 , i.e., terminating at these two points and containing no other kv j . Let the strictly convex body D i be such that ∂D i = Γ i . Let G i k be the bounded region whose boundary is kS i ∪ S i k ; this is well defined since both curves terminate at kv i and kv i+1 .
in view of (4.14) applied to D i for r = k ∈ Z + .
Next, note that S i k ⊂ ∂D i k . As a result, ds < Ck 2/3 (4.21) by Lemma 4.3 with D = D i .
Next, we turn to the difference of the last two integrals in (4.19). Recalling (4.4), we write this as follows: i∈Jc w∈L\{0} where n is the outward unit normal to k∂Ω and ∂Ω k in the first two integrals, whileñ is outward where the estimate follows from (4.12) by replacing Ω of Proposition 4.1 by D i ; the constant C is independent of k. Since the sum w∈L\{0} |w| 3 ϕ(|F w|) converges absolutely by hypothesis, so does the double sum in the previous equation; therefore This together with (4.20) and (4.21) shows thatR(k) = O(k 2/3 ). The normal is irrational a.e. on ∂Ω c . Consequently k∂Ωc γ • ds = k∂Ωcγ ds = k ∂Ωcγ ds, and (4.16) follows from (4.18), since (4.6) holds.

Remark 4.5.
It is interesting that in cases where the normal is rational on a subset of ∂Ω of positive measure, the dilation factors are required to be integers. In contrast, the result of Proposition 4.1 (where the normal is irrational almost everywhere on ∂Ω) is independent of the sequence of dilation factors. In one dimension it is known [Mo] that the coefficients in the asymptotic expansion of the energy depend on this sequence. It should be kept in mind that there is no counterpart in one dimension of an irrational surface, which is purely a higher-dimensional occurrence. The reason for the difference between the rational and irrational cases is the order of the lattice point remainder term.
Proofs of the Wulff theorem associated with surface energy minimisation [Fo] over domains with given measure typically rely on continuity of the surface energy density with respect to the unit normal; see [DP] and Remark 5.1 for a weaker alternative. Perhaps surprisingly, the extended surface energy densityγ(F, ·) : S 1 → R exhibits a dense set of discontinuities as we show next.
Proposition 4.6. Suppose ϕ is as in Proposition 3.4 and fix F ∈ M 2×2 + . Then (ii)γ(F, ·) : S 1 → R defined in (4.15) is continuous at n ∈ S 1 I , discontinuous at n ∈ S 1 R and differentiable at most on a subset of S 1 I of measure zero.
Proof. Arrange the elements of L \ {0} in a sequence: {w j }, j = 1, 2, . . . , such that |w j+1 | ≥ |w j |, and define g j (n) = (−1/4)ϕ(|F w j |)|w j · n| for n ∈ S 1 . Then clearly g j : S 1 → R is Lipschitz on S 1 and (formally for the moment) γ • (F, n) = ∞ j=1 g j (n). Now since |g j | ≤ M j = |ϕ(|F w j |)| |w j | on S 1 and the series ∞ j=1 M j = w∈L\{0} |ϕ(|F w|)| |w| converges in view of Lemma 3.3, then G k (n) = k j=1 g j (n) converge uniformly as k → ∞ to γ • (F, n) on S 1 by the Weierstrass M test. Since n → |w · n|, n ∈ S 1 is Lipschitz with constant |w|, the Lipschitz constant of G k is bounded The uniform convergence of the G k together with the uniform bound on their Lipschitz constants guarantee that the limit function γ • (F, ·) is also Lipschitz on S 1 and (i) holds.
To show (ii), consider the function In other words, letting n = (ν 1 , ν 2 ) ∈ S 1 , otherwise. (4.22) Then one hasγ (4.23) By (i), it suffices to prove that h is continuous at irrational n and discontinuous at rational n to show the continuity part of (ii). In fact, h is very similar to the Thomae function T (x) = 1/q for x = p/q, p, q coprime integers (x rational), and zero for x irrational; see e.g. Proposition 4.1 in [Sa]. Adapting these results to the h is trivial in view of (4.22). Thus h is continuous at irrational n and discontinuous at rational n and so isγ(F, ·). Also h is nowhere differentiable by a simple adaptation of Proposition 6.1, [Sa]. Since by part (i) γ • is Lipschitz, it is differentiable a.e. on S 1 by the Rademacher theorem. Thenγ(F, ·) fails a.e. to be differentiable by (4.23). Also it is not differentiable at rational n as it is not continuous there.

A Continuous Surface Energy Density
There are two issues associated with the surface energy densityγ. The first issue is the lack of continuity ofγ(F, ·). This suggests that the surface energy minimisation problem, that of minimising the integral ∂Ωγ (F, n)ds over a suitable class of regions Ω with |Ω| fixed, may actually be ill posed.
Remark 5.1. The standard hypothesis for surface energy minimisation in three dimensions is continuity ofγ(F, ·) [Fo]. However, in two dimensions, as shown by Dacorogna and Pfister [DP], lower semicontinuity ofγ(F, ·) suffices. It is easy to show from (4.15) and Proposition 4.6 thatγ(F, ·) is indeed lower semicontinuous provided W (F ) < 0. The latter inequality is not unreasonable; for example, it is satisfied for values of F near the minimum of W (F ), when the latter is given by (1.3) with ϕ a standard Lennard-Jones potential.
The second issue is that the surface energy minimisation problem with densityγ(F, ·) is not physically appropriate, since fixing |Ω| is not the same as fixing the total mass, or equivalently, the number #(Ω ∩ L) of lattice points of Ω. If the minimisation were over the class of lattice polygons with fixed lattice point number, the appropriate constraint would fix |Ω| + ∂Ω 1/(2|n|)ds instead of |Ω|, by virtue of Lemma 2.2. For a lattice polygon, the lattice point remainder R(k) = #(kΩ ∩ L) − |kΩ| can be written as (5.1) Rosakis using Lemma 2.2. It seems that R is implicated in both issues raised above. Being O(k), it contributes to the surface energy and gives rise to the term 1 2|n| W (F ) in (4.15), (4.23), which is the one responsible for the lack of continuity ofγ. Also, surface energy minimisation over domains with fixed measure would seem to make physical sense only if their lattice point remainder #(Ω∩L)−|Ω| vanishes, so that constraining |Ω| fixes the lattice point number, hence the mass. One way to ensure this might be to seek a sequence of dilation factors r k ∈ R satisfying condition (1.4) imposed by [BBL], i.e., R(r k ) = 0. It is not clear for what choices of Ω this is possible, and we modify this approach in two ways.
First, we relax the condition R(r k ) = 0 and require instead that there is a sequence r k such that so that the lattice point remainder is of lower order than the surface energy, which is O(r k ). This is satisfied for the smooth regions with positive boundary curvature of Section 4, where the fact that R(r) = O(r 2/3 ) for any real sequence r → ∞ was exploited in proving Proposition 4.1. As a result, the density γ • in (4.5) is continuous in the unit normal by extension to the whole of S 1 ; see Proposition 4.6 (i).
Second, in case Ω is a lattice polygon, or the "mixed" region of Proposition 4.4, we rewrite the energy in terms of an "equivalent" region Ω(k) containing the same lattice points as the scaled region kΩ. Accordingly, from (2.2) it is clear that E{kΩ, y} = E{Ω(k), y}. Observe that, given the set of atoms that are within a convex region kΩ, there is some freedom in choosing an alternative convex region Ω(k) containing precisely the same atoms. By choosing Ω(k) in a specific way, we can ensure that the lattice point remainder of Ω(k) is of lower order than the surface energy. For lattice polygons this can be done as follows. The "interplanar" distance between adjacent parallel lattice lines with Miller normaln is 1/|n|. If Ω is a lattice polygon, construct Ω ′ by moving each side with Miller normaln i of ∂Ω outward by 1/(2|n i |), half the interplanar distance. Then extend the translated sides, so that they once more intersect in the same order as before. Thus Ω ′ is a rational polygon [BR] (not a lattice polygon) that contains the same atoms as Ω, with sides parallel to those of Ω and vertex angles the same as those of Ω. In general though, it is not a dilation of Ω, although Ω ⊂ Ω ′ . Performing the same operation on kΩ for each k ∈ Z yields Ω(k). Since the layers added to kΩ have measure equal to k ∂Ω 1/(2|n|)ds to dominant order, it follows from (5.1) that the lattice point remainder #(Ω(k) ∩ L) − |Ω(k)| = o(k), so that (5.2) is satisfied for r k = k ∈ Z. Writing the energy in terms of the modified region Ω(k), one arrives at the following representation:  The O(1) term is a correction due to intersection, in the neighbourhood of corners, of layers corresponding to adjacent sides, since directions and thicknesses of layers are k-independent. The second equality above follows from (2.10) of Lemma 2.2. The O(1) terms in (5.5) are actually constant (depend only on Ω and not on k) as is easily shown. This establishes the middle assertion in (5.3). Since the distance of adjacent lattice lines with normaln i is 1/|n i |, the added layers (whose thickness is half that distance) contain no new lattice points; thus the first assertion of (5.3) holds true, while the last is trivial. The first of (5.3) ensures that E{kΩ, y} = E{Ω(k), y}. Now (5.4) follows immediately from Proposition 3.4, (5.5) and the definitions (3.4) and (4.4).

Case (b):
Suppose Ω is a smooth region as in Proposition 4.1. Then choose Ω(k) = kΩ, to that (5.3) follows from [Hu] and note that (5.4) is the same as (4.5) with k = r ∈ R + .

Rosakis
Case (c): Let Ω comply with Proposition 4.4 . For each k let Ω(k) be the set obtained by moving only the flat sides kS i ⊂ ∂Ω f , i ∈ J f of kΩ outwards by 1/(2|n i |) (and discarding portions of the added layers that lie outside the curves Γ j near the endpoints where S i join curved sides of ∂Ω).
Remark 5.3. Provided that a sequence of dilation factors r k exists such that (1.4) holds (as assumed in Theorem 3 of [BBL]), it is possible to modify the results of the previous section to show that for the types of regions considered here, E{r k Ω, y} = r 2 k Ω W (F )dx + r k ∂Ω γ • (F, n)ds + O(1).
Thus the non-explicit surface energy density of Theorem 3 in [BBL] is now determined to be γ • defined in (4.4). The order of the error is due to the vanishing of the lattice remainder R(r k ). Apparently, it does not seem to be known for which regions Ω such a sequence of dilation factors r k exists.
We are thus led to the construction of the regions Ω(k) of Proposition 5.2, which are not dilations of Ω of the form rΩ, since they involve different translations of different facets.
Remark 5.4. Proposition 5.2 indicates that the appropriate problem of surface energy minimisation over regions of fixed mass involves minimising ∂Ω ′ γ • (F, n)ds over a suitable class of domains Ω ′ with |Ω ′ | fixed. The integrand γ • is now Lipschitz continuous in the unit normal as guaranteed by Proposition 4.6. Thus γ • can be used to determine the Wulff shape of the crystal. We must remark, however, that while (5.4) has the aforementioned advantages as regards surface energy minimisation, it is not appropriate as an asymptotic series in k, since the domains of integration depend on the latter variable. The appropriate asymptotic series remains (4.16).