Continuum limits of discrete isoperimetric problems and Wulff shapes in lattices and quasicrystal tilings

We prove discrete-to-continuum convergence of interaction energies defined on lattices in the Euclidean space (with interactions beyond nearest neighbours) to a crystalline perimeter, and we discuss the possible Wulff shapes obtainable in this way. Exploiting the"multigrid construction"of quasiperiodic tilings (which is an extension of De Bruijn's"pentagrid"construction of Penrose tilings) we adapt the same techniques to also find the macroscopical homogenized perimeter when we microscopically rescale a given quasiperiodic tiling.


Introduction
The question of what crystal shapes are induced by what kind of interactions has preoccupied researches since the beginning of the field of crystallography. Mathematically, the study of crystal shapes has been first put on a firm ground within the continuum theory, starting with the work of Wulff [49], later reformulated and extended by Herring [22] and others [16,31,48]; see also [44] and references therein, for the connection to anisotropic perimeter functionals. In the continuum study, the role of the microscopic structure of the material considered is not modelled explicitly, and one starts by studying a surface energy of the form where φ : R d → [0, +∞] is a 1-homogeneous convex function, E ⊂ R d is a finiteperimeter set, ∂ * E is the reduced boundary of E and H d−1 is the (d−1)-dimensional Hausdorff measure in R d . See the book [33] for details. The optimizer of P φ among unit-volume competitors gives the shape of an ideal crystal with anisotropy φ, and its shape is called the Wulff shape corresponding to φ, see [22,44].
In this work, we focus on the link between discrete energy-minimization models and the minimization giving rise to the Wulff-shape problem. We will think of a discrete crystal to be a fixed-cardinality minimum-energy subset of the vertices of a lattice, or a subset of tiles in a quasiperiodic tiling. The energies that we consider are sums of pairwise interactions that respect the periodic or quasiperiodic structure. We establish compactness and Γ-convergence results in which, when we scale down the lattice as we increase the cardinality of point or tile configurations, the discrete energy functionals converge to a perimeter functional as in (1.1). We then consider the effect of modifying the discrete interaction model at the microscopic scale, on the macroscopic limit Wulff shape obtained in the Γ-limit. This endeavor fits within the general theory of discrete-to-continuum limits for crystals and quasicrystals. See [4] for the triangular lattice, [9] for another approach for quasicrystals and [7] for a homogenization result on the Penrose tiling, and the discussion below for more related results.

Setting and main results for lattice energies
We consider a lattice L = M Z d ⊂ R d , where M ∈ GL(d). We define det L as | det M | whenever L = M Z d (we refer to [6] for a proof that this is a well-defined quantity). We consider configurations X N = {x 1 , . . . , x N } lying in L and we define the energy Here V : L → (−∞, 0] is a fixed potential which may quantify the fall in energy due to the formation of atomic bonds in a crystal, for example. We first consider the case where V vanishes outside of a finite subset N ⊂ L such that span Z N = L (see Section 2.6.1 for more general cases). We will be interested in the surface-type which counts the energy excess due to missing bonds. Indeed the energy (1.2) rewrites as E(X N ) = C E N + F(X N ), where C E = w∈L\{0} V (w) is a "bulk" term independent of the shape of X N . To every configuration X ⊂ L (not necessarily having N points) we associate, denoting by U (1.5) We consider the following convergence: given a sequence (X N ) N ∈N , we say that X N converges to a set E ⊂ R n if E N (X N ) → E locally in measure (also referred to as the "L 1 loc convergence", identifying sets with their characteristic function).
Theorem 1.1 (Gamma convergence for crystals). The functionals N − d−1 d F N Γconverge, with respect to the topology above, to a functional of the form P V := P φ V as in (1.1), where As we discovered after the completion of the preliminary version of this paper, this result had been already proven by Gelli in her PhD thesis [21] in a more general form (see also [8] and [1]). We decided to leave the proof (even if the ideas are very close to those in [21]) because in our simplified case some of the intricacies of the general case are not present, and moreover the argument will be referenced later in the proof of the quasicrystal case given by Theorem 3.3. In order to prove Theorem 1.1 we first show in Section 2.1 that it is sufficient to prove it when the lattice L is Z d . Then in the rest of Section 2 we prove it for L = Z d . We also note that for every finite perimeter set E in R d we have , as proved in Proposition 4.11, so that it is not restrictive to assume that V is symmetric. We also prove the following compactness result, which motivates the chosen convergence. there exists a subsequence X N k and a finite perimeter set E such that E N k (X N k ) → E in L 1 loc .
As a consequence of Theorem 1.1 and Proposition 1.2 we obtain the following. Corollary 1.3. Minimizers of F N converge locally in measure, up to rescaling and possibly a translation, to a finite perimeter set E that minimizes (1.1) for its own volume constraint. More precisely (as explained in Section 4), in the case of the anisotropy (1.6) this Wulff shape coincides with the Minkowski sum of segments given by (1.7) Remark 1.4. (i) In Section 2.6.1, in order to simplify the analysis in the case when N does not span Z d , we also introduce the following convergence, which has been widely used in the literature: given X ⊂ L we define the empirical measure µ N (X) := 1 N x∈X δ x/N 1/d . (1.8) Then by definition X N converges to a set E if the empirical measures µ N (X N ) converge to ½ E weakly as measures. Observe that, when considering subsets of Z d , this convergence is equivalent to the "L 1 loc " considered above, and to state and prove Theorem 1.1 we could equivalently use the functionals (ii) The restriction #X N = N in (1.5) is not necessary, and Theorem 1.1 would be true even without it. However we chose to put it so that the proof of the recovery sequence becomes more precise (we can construct sets with exactly N points), and so that we can talk about minimizers of F N (which without a cardinality constraint would be trivial) and thus state Corollary 1.3.
A particular case of Theorem 1.1 appears in [4] within the study of triangularlattice configurations in the plane. This global convergence result for discrete energy functionals was successively made more quantitative near the minimum in [38], who proved the N 3/4 -law for fluctuations near the minimizer in the hexagonal case (see also [34] for the 3-dimensional case and [10,36] more in general). Results describing the structure of configurations minimizing important discrete functionals in an "unconstrained" setting, i.e. without restricting the configurations to a lattice, are available in very few cases, in dimensions 2, 3,8,24 in different models (see [5,11,12,14,19,35,45,46]), and in these cases too, the optimal limit shape can be shown to coincide with the Wulff shape of the corresponding lattice. Of the above works, note that [5,12] work with a potential V which involves interactions beyond nearest-neighbors, motivating our choice of including general V in practice.
Our proof of the Γ-liminf inequality in Theorem 1.1 is based on splitting the contributions to the energy appearing in (1.3) into contributions from the single edge directions in the support of V , and taking the limit on each one separately with the help of Reshetnyak's theorem. The Γ-limsup inequality is by polyhedral approximation, like the one performed in a special case in [4].
One of the advantages of our method, especially for the Γ-liminf case, is that it has indicated us a strategy for treating the quasiperiodic case, via a relatively nontechnical discussion. We expect that the same strategy can extend to more general quasicrystals and glass-like generalizations to configurations constructed from configurations of hypersurfaces.

Setting and main results in the quasicrystal case
The term quasicrystal refers to a class of generalized lattices that are not periodic, but possess some form of quasiperiodicity. A great interest in these kinds of arrangements arose in crystallography in the 80's (see e.g. [17,26,30]), when it was famously observed by Shechtman [42] that some metal alloys create diffraction patterns with five-fold symmetries that could not be explained by periodic arrangements of atoms. These patterns were then explained exactly by a "quasiperiodic" structure that never repeats but has atomic Fourier transform (thereby indicating some version of periodicity). We refer to [41] for a mathematical introduction to quasicrystals. In the mathematical community the most famous quasiperiodic arrangement is arguably the Penrose tiling, a tiling of the plane created with the use of two kinds of rhombuses as described by De Bruijn in [15]. An algebraic precursor of the idea of a quasicrystal can be traced back to Meyer [37].
Quasicrystals can be satisfactorily modeled by a variety of alternative non-equivalent mathematical definitions, depending on the precise focus of a given model or theory, and we refer to [20,28,29] for a comparison between some (but not all) of the different possible definitions.
The choice of definition which allows to directly connect to the theorems in Section 1.1 is the so-called "multigrid construction" of quasicrystals, introduced as far as we could find in De Bruijn [15] in the second part of his paper, and extended in [20] to the setting considered here. Wulff shapes of quasicrystals have been compellingly characterized in the physics literature for example in [23,24], and thus our work here consists in writing complete proofs of the energy convergence which formalizes [23,24] within the theory of Γ-convergence, and slightly generalizes the results to the full multigrid setup [20]. Amongst other constructions of large classes of quasicrystals, we mention the cutand-project method, also formulated in [15] for the Penrose tiling case. Gähler and Rhyner [20] extended the Penrose description from De Bruijn and proved that tilings by parallelohedra can be constructed by one method if and only if they can be constructed by the other.

Energy and Γ-convergence result in the quasicrystal case
We now define perimeter energies on the space of finite unions of tiles in quasiperiodic tilings. We consider a given quasiperiodic tiling T of R d by parallelotopes ("tiles"), which are produced through the "multigrid construction". This means that the polyhedral complex of the tiling T is dual to the one formed by dissecting R d by a number of families of parallel hyperplanes. We refer to Section 3.1 for the precise definitions, and to Figures 1, and 2 for a simple example. We consider the energy of a set T ⊂ R d which is a union of finitely many tiles from T as the following perimeter functional: where ν is the normal to ∂T , and where w is a nonnegative weight function which is defined on the finitely many possible directions of ν. The crucial point is that this functional can be rewritten in the form for some non-negative potential W (with a sign convention opposite to the crystal case), and where the sum runs among all the possible directions of ν. A more detailed description is given in Section 3.2, in which formula (1.10) is repeated in (3.17) and reexpressed in a dual space in (3.18). This rewriting makes it possible to interpret the perimeter-type functional (1.10) as a superposition of interaction potentials of type (1.3), which we know how to handle. We then define (1.11) and then we have the following analogue of Theorem 1.1.
Theorem 1.5. The functionals F N defined in (1.11) Γ-converge, with respect to the L 1 loc topology, to the functional P W = P φ W , for φ W of the same form as (1.6) (see (3.19b) and the preceding discussion for the definition). Moreover, if W (ν) > 0 for every normal ν to a tile, sequences with equibounded energy are compact in L 1 loc . In particular, the minimizers X N of E W from (1.10), amongst N -tile configurations converge, up to rescaling, to a finite perimeter set E that minimizes the anisotropic perimeter P W .

The search of general Wulff shape constructions
In both our main theorems, we find that the limit anisotropic perimeter functionals correspond to φ which is a sum of terms of the form W (v)| v, · | with W (v) > 0 (possibly take W = −V for the crystal case). What are the possible Wulff shapes corresponding to these perimeters? Surprisingly, this natural question, which is thoroughly investigated in physics papers [23,24], does not seem to be well-studied in the mathematical literature. Therefore, in Section 4 we collect and describe the basic results in this direction and provide a few new examples to illustrate some phenomena. It follows from classical convex geometry (see [40]) that if φ W = φ W 1 ± φ W 2 where W 1 , W 2 are finitely supported potentials and φ W is defined as in (1.6), then the Wulff shape W W is the Minkowski sum/difference of the Wulff shapes of φ W j , j = 1, 2. Therefore for positive finitely supported V the corresponding Wulff shape is a Minkowski sum of segments, sometimes named a zonotope, and we directly have the following: Theorem 1.6 (Wulff shapes under constant-sign potentials). A set W ⊂ R d is obtainable as limit optimal shape from energies as in (1.3) for nonpositive V with finite support, or as in (1.10) for positive W , if and only if W is a zonotope.
It is often claimed in the mathematical literature that this remains the case even for signed finitely supported W . However we point out that (as observed amongst others by [24] in the quasiperiodic case) this is not true. Indeed we can find examples of simple signed W , both coming from lattices and from multigrid quasicrystals, in which the Wulff shape is not a zonotope. This indicates that real-world crystalline shapes such as the pyritohedron or general truncated octahedra, are possible within our model, for signed W .
We leave as an interesting future direction the extension of Γ-convergence results such as Theorems 1.1 and 1.5 to the case of general signed W (with W = −V for Theorem 1.1). We believe that the results remain true whenever W is such that φ W is strictly positive on the unit sphere. The techniques considered here need to be refined for that case and we leave the extension to future work.

Structure of the paper
In Section 2.1 we prove equation (2.5) which allows us to change coordinates and reduce the study of lattice energies to the case of Z d . Section 2 is devoted to the crystal case, with the proof of Theorem 1.1 and of Proposition 1.2, as well as to some degenerate analogues, described in Section 2.6.1. Section 3 is devoted to the quasicrystal case, with the proof of Theorem 3.3. Section 4 is devoted to the study of possible Wulff shapes that can appear as continuum minimizers for the limit energies from our main theorems. It also describes the first steps for the study of Wulff shapes for signed potentials W . Finally, Section 5 includes sketches of some direct generalizations of our results and a short discussion of what seem interesting open directions for future work.
Acknowledgements. The authors wish to thank Maria Stella Gelli for pointing out that Theorem 1.1 was a specialization of previous results [1,8,21]. GDN has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme under grant agreement No 757254 (SINGULARITY). MP is supported by the Chilean Fondecyt Iniciación grant number 11170264 entitled "Sharp asymptotics for large particle systems and topological singularities".

Lattice case
In this Section we prove Theorem 1.1. We first prove in Subsection 2.1 that we can reduce to the case of L = Z d . In Subsection 2.2 we split the energy according to the direction of the bonds appearing in (1.3). We then relate each of these energies to a suitable anisotropic perimeter of a certain set associated to X N . In this way we can rewrite the total energy E N as the superposition of anisotropic perimeters in different directions v (those for which V (v) is non zero) of certain approximations of X N as union of cylinders with axis along v. This will help us deduce in Subsection 2.3 the lim inf inequality from the lower semicontinuity of perimeter-type functionals. Then in Subsection 2.4 we prove the lim sup inequality by approximation with polyhedral sets through a direct construction.

Reduction to the integer lattice
Every lattice L can be expressed as L = M Z d with M ∈ GL(d). Moreover there is a direct correspondence between configurations in L and in Z d : and thus to find the Γ-limit for a general lattice L we can translate the problem in Z d , find the Γ-limit there, and then go back to the original lattice. The following result shows that when translating the problem from Z d to any lattice L, the perimeter functional P V in (1.1) behaves well.
Proof. We apply the area formula [3, Thm. 2.91] to the map M −1 and the (n − 1)rectifiable set ∂ * E: where J ν E (y) ⊥ M −1 is the Jacobian determinant of M −1 restricted to the hyperplane orthogonal to ν E (y). Now we claim that ν E (M −1 y) = M * ν E (y) |M * ν E (y)| . Indeed, choose a basis e 1 , . . . , e d−1 for the tangent space to E, which is equal to ν ⊥ E . Then M e 1 , . . . , M e d−1 is a basis for the tangent space to E, which coincides with ν ⊥ E . Therefore we have 0 = ν E , M e i = M * ν E , e i . As M * ν E is orthogonal to e 1 , . . . , e d−1 , it must be a multiple of ν E (M −1 y), as desired. We thus obtain that (2.1) equals We have We now want to prove that We first claim that J The tangential jacobian J V M of a linear map M with respect to a hyperplane V can be computed in the following way: consider an (n − 1)-dimensional unit cube Q inside V , then J V M = H n−1 (M Q). Consider now a cube expressed as a Minkowski Therefore, using also the fact that Due to the fact that restriction to a subspace and inverse commute for invertible maps, we also find

This proves the claimed identity (2.3). Equation (2.2) thus becomes
This is a sublattice of Z d of dimension d − 1 (see e.g. [6, Ch. VII.2]), lying in the linear subspace (Rv) ⊥ . As a lattice of such subspace, we can select a basis b 1 , . . . , b d−1 and define its fundamental cell We finally consider the sublattice and its fundamental cell Given a translation vector τ ∈ Z d ∩ U v we also define the translated sublattices L v,τ := τ + L v . For fixed v, a volume estimate using the fundamental cells gives that the number of distinct sublattices of the form L v,τ is given by

Splitting of the energy
For every fixed v ∈ Z d \ {0} and every fixed translation vector τ ∈ Z d we define an energy functional by setting for X ⊂ L which yields a decomposition of the total energy as We observe that for every fixed v and τ , the only terms appearing in F v,τ are of the type V (v), so that effectively

Relation with anisotropic perimeter
Now for every given v and τ we associate to X a set E v,τ (X), made of translated copies of U v , and we relate the energy F v,τ (X) to a suitably defined anisotropic perimeter in direction v of E v,τ . More precisely, to every configuration X ⊂ Z d , and to every fixed v ∈ Z d \ {0} and τ ∈ Z d , we associate the following set and its contraction by a factor of N 1 d : Every paralelotope x + U v in the above definition of E v,τ N (X) has exactly one face for which normal vector ν there holds v, ν + > 0, and v = |v|ν for this face. Furthermore, such face has area H d−1 (U L v ⊥ ) and all other faces are parallel to v. Hence the contribution in (2.6) coming from bonds contained in L v,τ can be expressed as an anisotropic area functional in direction v, as F v,τ (X) = P v (E v,τ (X)), where for finite-perimeter E ⊂ R d we define P v (E) as follows, ν(x) being the normal to the essential boundary ∂ * E (cf. [3, Definition 3.60]): For a general measurable E ⊂ R d and v, τ as above we define: This functional automatically satisfies 3)), and will be used in the proof of Γ-convergence.

Liminf inequality
We now put together two facts: • By Reshetnyak's Theorem [3,Thm. 2.38], functionals as in (2.8) are lower semicontinuous with respect to the L 1 loc convergence of sets, since the 1homogeneous extension of the integrand from (2.8) is a convex function.
From these two facts, given any sequence This proves the lim inf inequality.

Limsup inequality
in such a way that they have positive scalar product with ν, is equal to . Proof. We start by proving (2.11). If ν, v ≤ 0 then no bonds [x, x + v), x ∈ L v,τ cut the hyperplane ν ⊥ in a direction making positive scalar product with ν, and also (2.11) gives zero contribution, thus the thesis follows. From now on we concentrate on the case ν, v > 0. In this case an edge [x, x + v) as above intersects ν ⊥ if and only if it intersects it while having positive scalar product with ν.
We pass to a model case for clarity first: There exists an invertible affine map x → Ax − τ which sends L v,τ to Z d and sends the v + τ to e 1 . The union of all edges of the form [x, x + v), x ∈ L v,τ is sent by this map to the countable union of lines R × Z d−1 and the hyperplane ν ⊥ is sent to a hyperplane transverse to all these lines since ν ⊥ ∦ v. The hyperplane Aν ⊥ − τ then meets each line R × Z d−1 exactly once. Coming back to the original coordinates, we find that the number of intersections of a set in ν ⊥ with segments [x, x + v), x ∈ L v,τ is equal to the number of intersections with the lines L v ⊥ + vR + τ . In the case of the set Int[P, L v,τ ] this number is also equal to the number of cylinder sets (U ⊥ v + τ + x) + Rv appearing in the union (2.10), since each such cylinder set contains exactly one such line. We now claim that for each y ∈ R d there holds We can apply the above for y = x + τ corresponding to the ters in (2.10) to obtain (2.11) by (d − 1)-dimensional volume comparison. To prove (2.13), we proceed by an elementary reasoning, although faster alternative proofs are possible. Note that it suffices to prove the case y = 0. Then observe that ν, v + = ν, v in our case. If we cut the infinite cylinder U ⊥ v + Rv into equal parallelotopes by hyperplanes orthogonal to ν passing through the equally spaced points Z v |v| , then the average number of parallelotopes per unit length in the direction v/|v| is the same as if we had cut it by hyperplanes orthogonal to v, thus the cut-out parallelotopes have equal volumes in the two cases. Since the volume of a parallelotope is equal to its basis area times the height relative to that basis, we get (2.13) directly.
We now prove (2.12). To do this, we note that any cylinder U v ⊥ + y + Rv that meets ∂P must also be included in its R-neighborhood for R larger or equal to the diameter of the intersection of U v ⊥ + Rv with ν ⊥ . The latter is bounded by by the same reasoning as in the first part of the proof. Then as P is a convex polytope, we can use the Hausdorff measure of ∂P to control the volume of its R-neighborhood.
We now pass to the construction of the recovery sequence. For a given measurable set E we need to construct sets We first prove the statement under the further assumption that ♯N < ∞.
Step 1. Approximating E by polyhedral sets of volume 1. By a classical approximation result we approximate E by polyhedral sets, more precisely we find a sequence of polyhedral sets E j such that where P is the standard perimeter. By Reshetnyak's theorem [3, Theorem 2.38] the same convergence holds for any anisotropic perimeter functional such as P V . Dilating by a factor converging to 1 we can also impose that |E j | = 1.
From now E is assumed to be a fixed polyhedral set of volume 1 and we work with its rescaling N 1 d E which has volume N .
Step 2. Approximation of polyhedral sets by discrete sets Y N . For large N ∈ N we will consider the discrete sets We have that 14) in which N E only depends on E (precisely, it depends on the rate of convergence of the limit in the definition of H d−1 (∂E)), and C d only depends on d.
We can then subtract or add a set of at most k N : From now on we focus on approximating the value of F(Y N ), and we will use the notations of Section 2.2.2. If we show that F v,τ (Y N ) satisfies the good bounds for all v ∈ N , τ ∈ U v , then we can sum the bounds and use the triangle inequality for approximating F(Y N ).
Step 3. Approximating F v,τ (Y N ) by the contribution of interiors of faces. We fix v ∈ N , τ ∈ U v . We denote by P a face of the polytope N 1 d E. We consider the discretization of P adapted to F v,τ (Y N ) given by Int[P, L v,τ ] from (2.10). For each face P of N 1 d E, the total number of bonds congruent to v cut by the interior of P in L v ⊥ is controlled via Lemma 2.3 and is given by (2.11).
As P meets each of the parallel lines τ + x + Rv, x ∈ U ⊥ v at most once, the number of bonds cut by P is the same as the number of lines of this type cut by P . Also, note that each cylindrical set of the form contains exactly one of the above lines, and it is not counted within the contributions of Int[P, We bound from above the number of x ∈ L v,τ for which this happens as follows: where in the last step we used the fact that projections decrease Hausdorff measure.
Summing up our reasoning so far, we show that −F v,τ (Y N ) can be approximated by the sum of one V (v) for each one of the lines that meet the interiors of all faces P , and the error terms are bounded with the help of (2.16) and (2.11). Note that, in order to obtain the precise contribution in F v,τ , we multiply all terms by the negative factor V (v), which reverses the inequalities. We get: We use now (2.12) and the definition of P v from Section 2.2.3, to get that Step 4. Conclusion of the proof. We may re-express the sums in ( Then we find, using also the previous estimate (2.15), , as desired.

Compactness
To prove compactness we use that the energy (1.4), and apply standard compactness results for finite perimeter sets. We thus obtain compactness in the topology of local convergence in measure, which is the natural one to expect (there are sequences with equibounded energy that loose mass at infinity).
Remark 2.4. In the two-dimensional case, a connectedness assumption for the sequence of sets is sufficient to at least imply compactness up to translations, because perimeter controls diameter for connected sets. This is the requirement considered for instance in [4]. However in higher dimension this is not sufficient anymore, as is clear by considering a set with a long "tentacle".
We prove compactness in the case when span Z N = Z d . For a discussion about what happens if the span is not Z d see Subsection 2.6.1. In order to control the perimeter of E N (X N ) with the energy F(X N ) we need the following combinatorial lemma.
for some constant K depending on V , N and d only.
Proof. Since E N (X N ) is a union of cubes, the contributions to its perimeter come from (d − 1)-dimensional faces each of which is orthogonal to some canonical basis vector e i , and separate a point x ∈ X N and a point x ± e i ∈ X N , for some basis vector e i . Therefore By the assumption that span Z N = Z d , we can choose a basis b = (v 1 , . . . , v d ) of Z d made of vectors in N . This means that for every standard basis vector e i we can write e i = d j=1 a j i v j for some integer coefficients a j i . Repeatedly using the rule Putting everything together the conclusion follows with K = 2dA/c.
for some constant K not depending on N . The claimed compactness now follows from the standard compactness of finite perimeter sets, see e.g. [3, Theorem 3.39].

Sparse or lower dimensional lattices
In this section we consider the case that a canonical ambient lattice L ⊂ R d is given, but N := {v ∈ R d : V (v) = 0} does not span the whole of L. As usual, we restrict to the case L = Z d , as our problem is affine-invariant, and we consider the case span Z N = Z d . The discussion of the Γ-limit of our perimeter functional follows analogous principles in these cases, but the rescalings that we use depend on dim span R N as well. We start by dealing with the case that this dimension is < d.

Lower dimensional structures
In the case where W := span R N is a proper subspace of R d (say dim W = k < d) the Γ-convergence result is still true, but the L 1 loc convergence is not the natural one anymore. This is clear considering a potential with V (±e 1 ) = −1 and zero otherwise, where the optimal structures are 1-dimensional segments in direction e 1 , and therefore under the scaling (1.4) the sets E N (X) are not compact (they converge weakly to zero when seen as measures as in (1.8)). However it is clear that in this case, after a rescaling by a factor N −1 , the optimal structures converge weakly to segments. As we shall now sketch, there is a natural topology where compactness holds even when dim W < d: the topology of "L 1 loc -convergence on every slice parallel to W ". Actually for simplicity we will consider the empirical measures associated to X N and use the weak convergence of measures. We first partition Z d according to subspaces parallel to W : we consider It is clear that T W λ leaves P(W ) invariant. The replacement for the sets (1.4) is given by the empirical measures The measures µ W N (X) are supported on P(W ) and have total mass 1 if #X = N . It is clear that under the energy (1.2) the only interaction is among points living in the same subspace parallel to W , while different spaces do not interact. We define the following notion of convergence on the space of empirical measures. Definition 2.6 (Convergence on slices). Let W := span R N , dim W = k. A sequence of µ j of measures on P(W ) converges to µ if for every τ ∈ τ (W ) the measures µ j (τ + W ) converge to µ (τ + W ).
We then define the space of slices with finite perimeter, given by all sets of the form where we identify E τ with a (k-dimensional) finite perimeter set in W . On this space we define the functional where P V is defined by (1.1) and (1.6) (and where again we identify E ∩ Z with a finite perimeter set inside Z). We then have the following result. We omit the proof since it is a direct consequence of Theorem 1.1 and Proposition 1.2 applied to every slice. • Compactness: every sequence µ N such that sup N E(µ W N ) < ∞ admits a subsequence that converges in the topology of Definition 2.6 to a measure of the form µ = 1 det(span Z N ) H k E, for some E as in (2.21).
• Γ-convergence: the functionals (1.5) converge to 1 det L P W V (see (2.22)) under the same topology, where we identify a set E as in (2.21) with the measure H k E.

Sparse lattices
We finally briefly mention the case of sparse lattices, that is when The energy (1.2) of a finite X ⊂ Z d splits as a sum of energies of coset subconfigurations X ∩ (span Z N + x), x ∈ C, which do not interact. In this case the statement of the Γ convergence has to be modified to take into account the possibility that each such X ∩ (span Z N + x), when rescaled, converges to a different finite perimeter set. The limit space is thus composed of finite superpositions of finite perimeter sets, and the Γ-limit is a sum of energies of the type (1.1) (or (2.22) in case of lower dimensional structures) among such a decomposition. We do not write the detailed results and we limit to observe here that they can be deduced in case of need with little modifications, and with the help of the observations above, from the main results.

Quasicrystal case
The purpose of this Section is to extend the results of Section 2 to nearest-neighbor interaction energies for subsets of quasicrystal tilings.

Multigrid construction of quasicrystals
We here introduce the definition and notation for multigrids and for the associated dual tiling.

Multigrids in the "dual space"
We start with the construction of multigrids, which we imagine to live in a space "dual" to the one in which our quasicrystal tiling will live.
For g ∈ R d , γ ∈ R we define a hyperplane grid as follows: If G ⊂ R d is a finite set and to each g ∈ G we associate a number γ g ∈ R, we call the multigrid with normals G and translations γ the collection of hyperplanes We will sometimes avoid the mention of γ, G when they are clear from the context, and will write H g := H(g, γ g ) and H g,k := H(g, k, γ g ) in that case. With this convention, for g ∈ G and k ∈ Z we also define the slabs We will assume in the following that every multigrid we consider is in general position, that is, we suppose that any d hyperplanes {H kg,g : g ∈ J}, with J ⊂ G, ♯J = d must intersect at a single point, which we denote by x(J, k J ), where k J := (k g ) g∈J .
Moreover we will assume that no more than d hperplanes intersect simultaneously.

Quasicrystal tiling in the "primal space"
We next associate to a multigrid as in (3.2) a tiling of R d by parallelotopes, which will constitute our quasicrystal in the "primal" space. In particular to every intersection point in the multigrid we associate a parallelotope P (x) of the tiling.
Convention. We will interpret J ⊂ G not just as a set of vectors g but as an ordered set. This allows to denote without ambiguity det(J) := det(g 1 , . . . , g d ), if J = {g 1 , . . . , g d }.
To do this, we fix once and for all an ordering of G and for each subset J ⊂ G we will always order its elements in increasing order with respect to this ordering.
Before explaining the construction we need to introduce a further datum, namely a map G ∋ g → g ∈ R d (the collection of g will give the possible directions of the 1-dimensional edges of the parallelotopes P (x) cf. Figure 2), on which we require the following compatibility condition: Note that this implies that g → g is bijective. For the constrution it would be equivalent to require the product to be always negative, the important thing being that the sign does not depend on g 1 , . . . , g d . We observe that a standard choice to produce Penrose tilings in the plane is to set g = g ⊥ .
The vectors g, in conjunction with the multigrid M (G, γ) above, define a tiling of R d as follows. To every x = x(J, k J ) we associate the unique vector k = k G ∈ Z G which "extends" k J and satisfies the admissibility condition and we consider the parallelotope P (x) ⊂ R d defined as a Minkowski sum as follows: The condition (3.5) ensures that the union of all such tiles is a tiling of R d . This is claimed in [20, p. 270] for dimensions up to 3, where the authors say that they believe the same criterion holds in general dimension but did not check this in detail.
In the case of Penrose tilings, a proof is also present in [15], for very specific choices g, g. We provide a proof in the general case in Proposition 3.2.
Notations. We denote the set of all possible vertices x(J, k J ) as in (3.4) by X = X (G, γ) and the set of tiles P (x) defined as in (3.7) associated to x ∈ X (G, γ) by T = T (G, γ, G).

Bounded distortion between primal and dual space configurations
For X, Y ⊂ R d a map φ : X → Y has bounded distortion (or is a BD map) if sup x∈X |φ(x) − x| < +∞. We say that Y is BD to X if there exists a BD bijection φ : X → Y . Following the above notation, we will show that the bijection that relates the points x(J, k J ) as in (3.6) to the centers of the corresponding parallelotopes (3.7) from the associated tiling, is a BD map, up to an affine distortion: Lemma 3.1. Let X and T be the admissible multigrid points and tiles associated to M (G, γ) and to G as above. Let A : R d → R d be the affine map defined by and let φ : X → R d be the map associating to a vertex x ∈ X the center of the corresponding parallelotope P (x), namely Moreover, assuming the condition (3.5), the map A is invertible.
Proof. From the compatibility condition (3.5) and the Cauchy-Binet formula for the determinant we immediately obtain the invertibility of A. Indeed we can write where G and G are the d × (♯G) matrices with columns respectively running through ( g |g| ) g∈G and ( g |g| ) g∈G (in the same order with respect to the mapping ∼). Now Cauchy-Binet formula gives that where J runs through all subsets of {1, . . . , ♯G} with d elements and where G J stands for the d × d minor of G obtained considering only the columns corresponding to J. The condition (3.5) implies that every term in the right hand side of the above expression has the same sign (the sign of the determinant does not change if we multiply every column by a positive factor), and therefore we obtain the invertibility of A.
Next, if x satisfies (3.6) then k g = ⌊ 1 |g| x, g |g| − γ g ⌋. As a consequence, denoting by {·} the fractional part, we have (3.10) The maximum modulus of the right hand side of (3.10) is achieved at one of the vertices of the convex polytope on the right, and this leads directly to The conclusion follows.
The previous result allows us to transfer the problem from the dual to the primal space with a finite error. Roughly speaking, since we rescale down the space to find the Wulff shape, this finite error asymptotically vanishes in the rescaling, and this implies that the Wulff shape in the dual and in the primal space will be the same, up to an affine map.
We next prove the following: The beginning of the proof is standard within the theory of polyhedral complexes, however we include it for the sake of keeping the presentation self-contained and more transparent.
Proof. Let K be the polyhedral complex generated by the hyperplane arrangement {H k,g : k ∈ Z G }. We follow a standard construction in order to fix an explicit realization of the dual complex of K in the dual space. To do this, first consider the barycentric subdivision K 1 of K, i.e. the rectilinear simplicial complex with one vertex at the barycenter of each k-cell of K, for each k = 0, . . . , d, and whose d-cells are formed by the barycenters of f 0 , . . . , f d for all choices of facets f 0 , . . . , f d of K such that f 0 f 1 · · · f d . Now the d-cell of the dual complex K * of K corresponding to vertex x can be identified with the union of cells of K 1 which contain x; these cells of K 1 form the so-called "star" of x, denoted Star(x). Lowerdimensional cells of K * can then be obtained by intersecting d-cells.
Due to the assumption on our multigrid constants γ, which is chosen so that no more than d hyperplanes H g,k intersect (see the end of Section 3.1.1), the star in K 1 of each vertex x of K is a polyhedral cell complex isomorphic to the one given by the barycentric subdivision of the facet complex of parallelotope P (x), denoted by K 1 (P (x)). This fact can be proved by induction on d, together with the fact that the complex with cells P (x) is combinatorially dual to the one generated by the hyperplanes H g,k and we leave the details to the reader.
For each x vertex of K, from the cell-complex map φ x : Star(x) → K 1 (P (x)) we can define a mapφ x : Star(x) → P (x) which sends cells to cells, respects cell inclusions and is affine on the interiors of cells of Star(x) of any dimension. Asφ x is a piecewise affine bijection, it thus has Jacobian of constant sign. Furthermore, this sign is equal to the sign of det(g : g ∈ J) det( g : g ∈ J), which is constant by assumption (3.5).
We now glue theφ x to find a continuous piecewise affine mapφ : R d → R d . For this map to be continuous, we need to verify a compatibility condition. Let x = x ′ ∈ X be such that Star(x) ∩ Star(x ′ ) = ∅. Equivalently, x, x ′ are vertices of the same cell of K. In this case, let j be the lowest dimension of a cell of K containing both x and x ′ . This means that we find an index k G ∈ Z G such that the slabs S g,kg , g ∈ G as in (3.3) all contain x, x ′ either in their interior or in their boundaries, and there are precisely j choices of g ∈ G such that x, x ′ are in the same boundary component of S g,kg . From definition (3.7) it follows that P (x) ∩ P (x ′ ) have a common (d − j)dimensional facet. By further inspection of indices, we see that then both φ x and φ x ′ send Star(x) ∩ Star(x ′ ) to the barycentric subdivision of this common facet, and thusφ x andφ x ′ coincide on (Star(x) ∩ Star(x ′ )), as desired.
By estimating bilipschitz constants, we find thatφ is at bounded distance from φ from Lemma 3.1, and this shows thatφ is at bounded distance from the affine bijective map A from the same lemma, and thus it is surjective. Thusφ is a piecewise affine surjective map of R d such that for each x ∈ X the restrictionφ| Star(x) is bilipschitz with its image and (by (3.5)) preserves orientation.
If we had a continuous covering from R d to R d , it would have to be injective because R d is simply connected. For ourφ, we cannot ensure that it is a homeomorphism locally near the boundary of Star(x) for x ∈ X , however the degree-theoretic proof based on preserved orientations still works, as follows. If we take a large cycle Σ = ∂M in K around y ∈ R d , we have that the winding degree satisfies deg(φ, y, ∂M ) = 1, where deg(φ, y, ∂M ) is defined as the winding number ofφ (y ′ )−φ(y) |φ(y ′ )−φ(y)| around the origin as y ′ ∈ ∂M . The fact that this degree is 1 for ∂M far a way from y follows by approximation, from the fact thatφ is a bounded distortion of an affine bijective map. Thus by degree theory (for a proof see [32, Thm. 4]) we have for any ♯{φ −1 (y)} ≤ deg(φ, y, ∂M ) = 1, and thus bijectivity ofφ follows, and in particular the P (x), x ∈ X give a tiling, as desired.

Multigrid lines and parallelotope rails
We will denote the normalized directions of lines obtained by intersecting d − 1 hyperplanes of M (G, γ) as follows: (3.11) Due to condition (3.5), in fact for each line Rv + c, v ∈ V, there exists exactly one choice of J ′ v := {g 1 , . . . , g d−1 } as in (3.11) and these vectors are linearly independent, therefore the normalisation condition on v in (3.11) is well-posed. Note that the requirement that there exists c such that d−1 j=1 H 0,g j = Rv + c is equivalent to v ⊥ g j for j = 1, . . . , d − 1.
We note the following further properties: 1. If v ∈ V has uniquely associated hyperplanes directions J ′ v = {g 1 , . . . , g d−1 } ⊂ G as in (3.11) (with indices respecting the induced ordering fixed on G) then to v ∈ V we associate in the primal space the vector v ∈ R d determined by the conditions v ⊥ g j for j = 1, . . . , d − 1, | v| = 1 and det( g 1 , . . . , g d−1 , v) > 0 Again this is a good definition because, as seen before, the mapping G ∋ g → g is bijective and because we assumed the condition (3.5). It follows that { v : v ∈ V} are the normal vectors of the faces of parallelotopes P (k) as in (3.7), where v is normal to the faces generated by g 1 , . . . , g d−1 . Also, all normal vectors to the faces of P (k) are generated in this way, up to a sign.
2. With the above notation, the successive intersections of Rv with the hyperplanes H g,k , g ∈ G \ {g j } d−1 j=1 correspond to "parallelotope rails", i.e. chains of parallelotopes T = P (x) ∈ T (see (3.7)), in which neighboring parallelotopes have in common exactly the face with normal v.
3. Such rails were described for d = 2, 3 in [23,24]. We will also call a dual rail the set of collinear points in the dual space corresponding to tiles forming a rail.
Note that dual rail directions are in direct correspondence with the vectors v ∈ V, however in the primal space the vector v in general is not giving the direction of the parallelotope rail corresponding to v ∈ V: this direction is instead parallel to G G T v.

Density of multigrid sublattices and of rails
For each J ⊂ G of cardinality ♯J = d, we form a nondegenerate translated lattice Λ J ⊂ R d by intersecting d-ples of planes from the families corresponding to g ∈ J (cf. equation (3.1)): (3.12) in which {p J } = ∩ g∈J H g,0 is determined by the translation numbers γ g , g ∈ J from the definition of our multigrid, and the notaton Λ * := {z ∈ R d : z, p ∈ Z, ∀p ∈ Λ} denotes the dual lattice of Λ. In particular where like in (3.9), G J is the matrix whose columns are ( g |g| ) g∈J . We will also use later the following geometric consideration: if v is one of the multigrid directions (3.11) that generate Λ J , then v ⊥ intersects the multigrid lines parallel to v in a sublattice Λ J ′ v (with unit cell denoted by U Λ J ′ v ) which is related to Span Z g |g| 2 : g ∈ J ′ v by a formula like (3.12). Therefore, analogously to (3.13), for each J ⊃ J ′ v , ♯J = d, there holds: (3.14) Given a fixed rail direction v, for every possible choice of J such that J ⊃ J ′ v , the lattice Λ J has a generator λ v,J v parallel to v, with λ v,J > 0. Furthermore we have which shows that the leftmost quantity does not depend on J, but only on v. We will use this fact in the proof of the lim inf inequality.
Furthermore, we find that if P is a polygon in the hyperplane ν ⊥ , where ν ∈ R d is a unit vector, then the number of dual rails that meet P per unit area of P (and in the limit of "large P ") equals .
where J ′ v = {g 1 , . . . , g d−1 } ⊂ G denotes the set of vectors as in the definition (3.11) of v ∈ V. The above equality follows from the fact that v ⊥ g for g ∈ J ′ v and |v| = 1 by the definitions of J ′ , V, and is the multigrid analogue of (2.11).
The vertices of the multigrid are the union of all Λ J as above, and since the traslations γ g are such that no more than d planes from the multigrid meet simultaneously, then this union is disjoint.

Energy functional in primal and dual spaces
Let T be our quasiperiodic tiling of R d . The parallelotopes from T have edges which are translations of the segments [0, g] for g ∈ G. The (d − 1)-dimensional facets of elements of T have H d−1 -measures equal to | det( g : g ∈ J ′ )| 1 , where J ′ ⊂ G, ♯J ′ = d − 1, indexes their set of edges.
The same information as determined by the mention of J ′ can be equivalently encoded in terms of the associated rail directions v J ′ , i.e. using the observations and notation of Section 3.1.4. In this case we define a potential on the multigrid directions, W : V → R + which associates to v ∈ V defined in (3.11), the weighted surface area given by where J ′ v ⊂ G is related to v ∈ V via (3.11), and w( v) is the weight which we associate to faces with normal direction v, which is the "primal space" vector corresponding to v.
Note that here we use the opposite sign convention for W compared to the case of lattices (in which V was nonpositive), which seems more natural in this case.
With notation (3.16) we then define the energy of a set T ⊂ R d , which is a union of finitely many tiles, as follows: (3.17)

Energy in the dual space
To a finite union T of tiles from T there corresponds a finite set of vertices X ⊂ X and to faces of tiles of T there correspond edges between vertices in X . We can rewrite the energy (3.17) in terms of the dual space elements as where (x, y) is the open segment in R d with endpoints x, y. To prove the last equality in (3.18), note that for each v, in (3.17) we are counting pairs {x, y} corresponding to adjacent tiles sharing a face with normal vector v, only one of which is in T , which is a reformulation of (3.17). The rewriting (3.18) corresponds to considering a particle interaction between vertices of the multigrid, where only pairs of first neighbours (with respect to the natural graph structure of the multigrid given by its 1-skeleton) are taken into account.

Rescaled energy functionals and Γ-convergence statement
We now recall the definition of the functionals that allow us to formulate our Γconvergence result: • for each N ∈ N we define otherwise.
• As a limit of F N we find the desired perimeter functional: (3.19a) where, with A = G G T as in Lemma 3.1, Λ J as in (3.12), and U Λ J ′ v as in formula (3.14), we have Our main result then is the following.
Theorem 3.3. The functionals F N Γ-converge, with respect to the L 1 loc topology, to the functional P W . In particular, the minimizers X N of E W from (3.18), amongst N -tile configurations converge, in a weak sense, up to rescaling, to a finite-perimeter set E that minimizes P W .

Density of subtilings
We now perform a "statistical analysis" of the tiling, counting how many tiles of every kind appear on average in a fixed portion of the space, proving a formula for the asymptotic density of each subtiling. We then apply these estimates to prove the main result of this section, namely Proposition 3.8, which gives a correspondence between convergence of a sequence in the primal space and convergence of suitable sublattices in the dual space. Given the whole set of tiles T and J ⊂ G, ♯J = d, we denote by T J the family (subtiling) composed by all tiles corresponding to points of the sublattice Λ J , namely We denote by T J = T ∈T J T the corresponding union. Given a subset X N ⊂ X and the associated union of tiles T N = x∈X N P (x), we define X J N := X N ∩ Λ J and T J N := T N ∩ T J . Lemma 3.4 (Criterion for measure convergence). Let E, E N ⊂ R n be sets of finite measure. Then E N → E in measure if and only if both of the following hold: Proof. We only prove that (i) and (ii) imply the convergence, because the other implication is trivial. Assumption (ii) implies the convergence of the L 1 -L ∞ pairing ½ E N , ½ F → ½ E , ½ F when F is a ball. Thanks to Vitali's covering theorem, and using also assumption (i), we extend this convergence to any measurable F . Moreover, given any L ∞ nonnegative function g, we can write Since the functions ½ E N ∆E are positive, this implies strong convergence to 0, which is equivalent to |E N ∆E| → 0.
The usefulness of the previous criterion stems from the fact that in order to prove convergence in measure, it is sufficient to prove an estimate for the volume on balls, without regard to where the set is located inside the balls. Using the decomposition of the multigrid in sublattices we will be able to precisely estimate the number of points of a given configuration inside balls, and using the bounded distortion property between primal and dual space we will transfer this information back and forth.
Remark 3.5. We also mention the following result due to Visintin [47], that we could have used in the conclusion of the previous lemma: if a sequence u N ∈ L 1 (R n , R m ) is converging weakly to u, and if u(x) is an extremal point of the closed convex hull co({u N (x)} N ) for a.e. x, then u N → u strongly in L 1 . In our case the assumption on the convex hull is verified since all functions are characteristic functions.
Recall the definition of the matrices G, G and G J , G J given below (3.9), and of the affine map A with linear part G G T given in (3.8).
In particular the subtiling T J has asymptotic density ρ J . Moreover, for every measurable F with finite measure, we have Proof. Let D be the diameter of the tiles in T J . We have that The second term is bounded by |B R | − |B R−D | = O(R d−1 ). In the first term, every tile is entirely contained in B R and therefore contributes with |T | = | det( g : g ∈ J)| = (Π g∈J |g|) det G J to the sum. We thus obtain To estimate the cardinality we make use of the bounded distortion property given by Lemma 3.1. It follows that Here we have used that, thanks to the bounded distortion (Lemma 3.1), the tile P (x) associated to x is at most at a bounded finite distance from Ax, and therefore the error obtained counting the points x instead of the tiles P (x) is of order O(R d−1 ) (that corresponds to the number of points in a finite neighbourhood of ∂B R ). Putting together the last two equations we obtain (3.20).
The last statement about measurable sets follows from a rescaling and an application of Vitali's covering theorem.
In the previous lemma we proved the existence of a density for every subtiling. In the following lemma we prove a more general statement: whenever a sequence E N is converging in measure to a set E, then every "restricted subtiling" E N ∩ (N −1/d T J ) is uniformly spread inside E N , with the same density as the global one ρ J .
Lemma 3.7. Let E be a measurable set of finite measure. If E N → E in measure then for every measurable set F we have for ρ J as defined in (3.20).
Proof. From the convergence in measure we deduce that |E N ∩ (E ∩ F )| → |E ∩ F | for any F . Moreover and the last term converges to J ρ J |E ∩ F | = |E ∩ F | by Lemma 3.6. Since the previous inequality holds term by term, we obtain that every single term has to converge to its upper bound, that is Since the limit of Proposition 3.8 (L 1 -correspondence between primal and dual). Let T N = x∈X N P (x) be a sequence of tiled sets, and suppose that N −1/d T N → E in measure for some set E of finite measure. Then, setting We make use of the characterization of convergence in measure given by Lemma 3.4. The aim is thus to prove that, for every ball B, We first rewrite the left-hand side in terms of T J N , using the bounded distortion property. We have which converges to | det Λ J | | det( g: g∈J)| |E ∩ AB| = | det Λ J | | det( g: g∈J)| | det A| |A −1 E ∩ B| by Lemma 3.7. Recalling (3.13) and the definition of ρ J (see (3.20)) we obtain We have thus proved (3.22), and by Lemma 3.4 we reach our conclusion.

Compactness
Compactness of sequences of tile unions T N with equibounded F N -energy as defined in (1.11) is a direct consequence of compactness of finite perimeter sets, see [3,Theorem 3.39]. Indeed the assumption that W (v) > 0 when v is a normal to a face in the tiling directly implies that the functional F N bounds the perimeter of T N .

Liminf inequality
In order to prove the lim inf inequality we will analyze every sublattice Λ J separately, and for each we will prove the lower bound (3.29) concerning the energy of the bonds in a fixed direction v ∈ V. To put these bounds together and prove the full lim inf inequality, we will need the following combinatorial lemma. Let us start by defining the edge-perimeter sets EP v , EP J,v as follows: for S ⊂ Λ J we define and for X ⊂ X we let EP v (X) := {{x, y} ⊂ X : P (x), P (y) share a face, x ∈ X, y / ∈ X, x − y v} . Lemma 3.9. Let X ⊂ X be a finite set and let v ∈ V. Then Proof. Consider a multigrid direction v ∈ V and note that there are ♯G − d + 1 distinct choices of J such that for some λ J > 0 the vector λ J v is a generator of Λ J . These J are obtained by adding any new vector from G to the (d − 1)-ple J ′ v , associated to v as in the definition of V.
We claim that the image of φ v,w contains EP J,v (X ∩ Λ J ) ∩ {{x ′ , y ′ } : x ′ , y ′ ∈ ℓ}. Indeed, let {x ′ , y ′ } belong to the latter set, and let's say that x ′ ∈ X ∩ Λ J and Then [x ′ , y ′ ] is subdivided by the multigrid points X into a concatenation of segments [x, y] ⊂ ℓ such that the first segment has starting point in X and the last one has end point outside X. Therefore there exists at least one segment in the concatenation that has one end in X and the other outside X, i.e. it is an element of E(X). It follows from the definition of φ v,w that φ v,w ({x, y}) = {x ′ , y ′ } in this case, proving the claim. More precisely, By taking a union over w ∈ G \ J v , the above claim implies that the multimap given by (3.27) Taking the union of (3.27) over the multigrid lines of the form ℓ = Rv + c, and observing that φ v is actually (♯G − d + 1)-to-one because distinct Λ J 's are disjoint, the bound (3.25) follows.

Proof of the liminf inequality
We consider T N which is a union of N tiles from T , and such that The corresponding X N ⊂ X in the dual space can be partitioned along the multigrid lattices, into the subsets X J N := X N ∩ Λ J . Recall the definition of the auxiliary sets E J N := x∈X J N (x + U Λ J ), first introduced in (3.21). By Proposition 3.8 we have the L 1 -convergence and we may then use the same setup as in the lattice case, Section 2.3, to obtain where we used equation (3.15). Summing over all J that contain J ′ v and using Lemma 3.9, we obtain Taking now the lim inf of the last sum and using (3.29) we obtain the sum of ♯G − d + 1 = ♯{J ⊂ G : ♯J = d, J ′ v ⊂ J} equal terms, which cancels the factor 1 #G−d+1 . We are left with: Finally we also sum among all v ∈ V, using that by definition (see (3.18)) (3.31) Finally, using the change of variables formula given by Proposition 2.1 for M = A and with V replaced by W • A −1 , we obtain exactly the lim inf inequality involving the functional defined by (3.19) and (3.19b).

Limsup inequality
The proof follows the overall strategy from Section 2.4, with a few additions. The statement that we prove, analogous to the one proved in Section 2.4, is the following.
For a given measurable set E there exists N 0 ∈ N depending only on E, T , X such that for N ∈ N,

Strategy comparison to the crystal case:
The main difference to the case of lattices treated in Section 2.4 is that the sets T N do not have volume proportional to N , as the tiles from T do not all have the same volume. On the other hand, the set X in the dual space is a union of lattices, due to the multigrid construction. This implies sharp volume bounds for each lattice, analogous to the ones from Section 2.4.
Due to the above, we set up our desired approximation via steps analogous to Steps 1-4 from Section 2.4, extended to the multigrid setting. This produces sets of tile centers X N , and we use the Bounded Distortion Lemma 3.1 to determine distortion bounds in an extra step 5, in order to show that our final tile set T N dual to X N satisfies In what follows, we will describe the main changes required for the quasicrystal proof, referring to the steps from Section 2.4 for details.
Step 0. Reduction to the case of A = G G = Id in Lemma 3.1. By the affine invariance of our functionals, we can rescale the multigrid by (G G) −1 and reduce to the above mentioned case. From now on we thus assume that G G = Id. This allows to simplify the notation and focus on more essential ideas.
Step 1. Approximating E by polyhedral steps of volume 1. This step is exactly the same as in Section 2.4.
Step 2. Approximation of polyhedral sets by discrete sets from the multigrid. Note that X is the disjoint union of the lattices {Λ J : J ⊂ G, ♯J = d}, in which Λ J has density 1/| det Λ J |. We then consider a rescaling of E (now assumed to be polyhedral and of volume 1) of cardinality close to N by correcting for the density of X , which is the sum of densities of Λ J : Condition (2.14) for the multigrid setting follows by applying the bound for lattices to all the lattices Λ J , with J ⊂ G, ♯J = d, and using the triangle inequality. We obtain a version of (2.14) with a constant depending on X due to this: As a consequence of (3.33), we can continue and possibly add or remove a cluster of multigrid points obtaining from Y N a set X N satisfying (2.15) with constants depending on X now, concluding this step.
Steps 3-4. Approximating E W (T N ) by the contribution of interiors of faces. Our decomposition of the energy will now use the rails from Section 3.1.5 to decompose the energy E W into contributions coming from all families of rails, which replaces the role of the contributions from parallel rays in Step 3 of Section 2.4. Note that ) (see the paragraph preceding (3.14) for this notation). With these substitutions we reach a quasicrystal analogue of (2.16), (2.17) and (2.18), with the only difference being that C ′ d now depends on the multigrid data T , and that in the quasicrystal case P v (N 1 d E) has to be defined as the contribution of v ∈ V to P W , i.e. the contribution of the summand in (3.19b) corresponding to v only. With this we conclude Steps 3-4. We note here a weakened simplified version of the final estimate analogous to (2.19), in which C 1 , C 2 only depend on T , X , E, W (see (2.19) for a more precise dependence in the lattice case): Step 5. Passage from the dual to the primal space. Note that we have E W (X N ) = E W (T N ), and the volume of tiles in T which correspond to the symmetric difference between (3.33). Thus in order Remark 3.11. Although not used here, it is interesting to note that by the Laczkovic characterization of bounded displacement [27], the BD property shown in Lemma 3.1 is actually equivalent to a property of the type (3.35). The original work of [27] treats the case of cubic tiles, while substitution tilings are treated in [43] based on the same basic idea, and probably the multigrid analogue can be treated similarly too, based on Hall's marriage lemma. We only prove the implication from Lemma 3.1 to (3.35) here.
Proof of Lemma 3.10: We denote E N := N We note that if C diam is the maximum diameter of a tile from T , and C BD := sup x∈X |ψ(x) − Ax| with the notation of Lemma 3.1, then As the constants C diam , C BD depend only on G, G, it follows that for N sufficiently large (depending on G, G and E) such that the bound (3.35) holds, with C depending only on G, G.

Conclusion of proof of Theorem 3.3
Recall that we have already proved the Γ-liminf inequality corresponding to the statement of Theorem 3.3 in Section 3.5.
For the Γ-limsup part, we use the X N as recovery sequence, and we showed in

Building complex Wulff shapes from simpler ones
In this section we briefly discuss the relation between the Wulff shapes of two potentials and the Wulff shape of their sum. We then show in Proposition 4.11 that for all our purposes it is sufficient to consider symmetric potentials, that is those satisfying W (w) = W (−w) for every w, since the energy is invariant under symmetrization of W .

Basic properties
We first relate the functional to equivalent formulations. Instead of φ = φ W as before, we consider first the case of a general convex positively 1-homogeneous φ : R d → [0, +∞). We recall (see [44] for a proof) that the Wulff shape (i.e. a shape homothetic to any minimizer of ∂ * E φ(ν E )dH d−1 under fixed volume constraint) is then given by the subdifferential ∂ − φ(0) at the origin: Remark 4.1. In case W φ as defined in (4.1) has zero volume, then the isoperimetric problem is not well posed, however the definition (4.1) still is the natural extension by continuity, of the definition for the non-degenerate case. Note that the Wulff shape W φ obtained in this case is the optimal shape for the case discussed in Section 2.6.1.
For any nonnegative φ we have that W φ is compact, convex and contains the origin. Moreover, the origin is an interior point of W φ if W φ has nonzero volume. Viceversa if K is a convex set containing the origin then we can define φ K (x) := sup{ x, y : y ∈ K}, and note that then φ W φ = φ and W φ K = K. If φ(x) = φ(−x) and φ(x) > 0 for x = 0, then φ defines a norm on R d , which is dual to the norm whose unit ball is K, and which has the polar body K * of K as unit ball. Common terminologies refer to φ K as the surface tension function of K or as the support function of K.

Sums and differences of potentials and of Wulff shapes
Given two potentials φ 1 and φ 2 a natural question to ask is what is the Wulff shape associated to φ 1 + φ 2 . As shown in [40,Thm. 1.7.5], the answer is the Minkowski sum of the Wulff shapes of φ 1 and φ 2 , i.e. we have that We note that a special case in which (4.2) gives important information, is that of , in other words W φ is a zonotope, i.e. a Minkowski sum of segments. Equivalently, a zonotope is a convex polytope all of whose 2-dimensional faces are centrally symmetric [40,Thm. 3.5.1]. More generally, a zonoid is a convex body which is the Hausdorff-distance-limit of zonotopes. A zonoid Z has in general support function φ Z representable as a superposition (4.3) in which µ Z is a positive measure. Using 1-homogeneity, we may further impose that µ Z be supported on the unit sphere S d−1 ⊂ R d , in which case µ Z is uniquely determined by Z (see [40,Thm. 3.5.3]).

Recall that Minkowski subtraction is defined by
, or A − B can be characterized as the maximal convex set C such that C + B ⊂ A. As a consequence, φ A−B is the convexification of φ A − φ B , i.e. the largest positively 1-homogeneous convex function bounded above by φ A − φ B . As already noted e.g. in [13, Lem. 1(b)] with a different notation, we then have that Therefore, for convex 1-homogeneous φ 1 , φ 2 such that φ 1 , φ 2 , φ 1 − φ 2 ≥ 0, we also have the following expression via Minkowski subtraction: Note that if φ(x) = φ(−x), j = 1, 2 and W φ is 2-dimensional, then it is centrally symmetric, and is a zonoid. Centrally symmetric Wulff shapes of dimensions ≥ 3 are generally not zonoids: see the examples below, in which φ is of the form (4.5) below.
The important case for our work is that of φ of the special form 5) for finite N ⊂ R d , and where V : N → R. We can always split the sum into two sums φ = φ + − φ − corresponding to the decomposition N = N − ∪ N + , so that sign(V (v)) = ± on N ± , then by (4.4) W φ = W φ + − W φ − . This shows that for φ as in (4.5), the Wulff shape W φ is the Minkowski difference of two zonotopes: this is in general not a zonotope, as shown by a few examples below.
Question: Are all convex centrally symmetric polyhedra expressable as the difference of two zonotopes? (This is trivially true in dimension 2, thus the question is really interesting for dimensions 3 and higher.)

Remark 4.3.
A statement that we found in several instances in the convex analysis literature (see e.g. [40,Cor. 3.5.7], which is based on [39]), is that if W φ is a polytope with φ as in (4.5), then it is necessarily a zonotope: this is false as shown by the examples below. The main correction required in [40] seems to be in the proof of [40,Lem. 3.5.6], when it is mentioned that two integrals over sets B, B ′ with respect to a weight ρ are equal: this can be false for general (possibly concentrating) ρ, such as ρ with atoms on B \ B ′ , as is our case. It is important to emphasize that the above theorem and proof remain valid for regular convex bodies, and counterexamples are possible only the case of polytopes and non-smooth convex bodies.
Example 4.4 (Quasicrystal non-zonotope Wulff shapes). One of the main new ideas in [24] was to allow negative weights, to produce new limiting isoperimetric shapes in quasicrystals. In particular, seemingly for the first time in [24], a possible choice of N and of a signed weight V over N is given, such that if φ as in (4.5), then W φ is a regular dodecahedron. As mentioned above, such W φ is not a zonotope however, importantly, it was observed in actual quasicrystalline materials. The choice of [24] is to take N = N + ∪ N − where N − are the directions of vertices of a regular icosahedron centered at the origin, and N + are the unit vectors v pointing towards the midpoints of the sides of the same icosahedron. Then they set V = −1 on N − and V (v) = 5/6 on N + . Also icosahedral Wulff shapes can be obtained with φ as in (4.5), see [24, Fig. 6].
Remark 4.5. Related to the above example, note that any choice also gives a dodecahedral W φ . Here the constant 5/6 is taken directly from [24], whereas the lower bound on c + follows by requiring φ V > 0. The allowed interval for c + is given by the conditions that φ > 0 and that only the faces with normals N − give supporting planes of W φ . Example 4.6 (Crystalline non-zonotope Wulff shapes). A. It is worth mentioning that non-zonotope crystal shapes can in theory be formed, and are actually observed in nature, also in the presence of true lattice-like microscopic structure, such as in pyrite, which forms a (non-regular) dodecahedral Wulff shape, also known as a "pyritohedron", and having vertices (±1, ±1, ±1) and cyclic permutations of 0, ± 3 2 , ± 3 4 .
This shape can be obtained with φ in the form (4.5), for example with We can set V + (±e j ) = 2/3 and V + (v) = 4 21 for the remaining vertices in N + , and then the Wulff shape of φ V + is the zonohedron (recall that by definition a zonohedron is a 3-dimensional zonotope) with set of edge vectors equal to those of the above pyritohedron. Note that W φ V + then has 12 decagonal sides corresponding to pentagonal pyritohedron sides, as well as 30 rectangular faces and 20 hexagonal faces. If now V − is constant on N − , the Wulff shape of φ V − is a zonohedron with sides equal to the normals to pyritohedron faces. If the constant value of V − is large enough, the Minkowski subtraction of W φ V − from W φ V + , which gives W φ V + −φ V − , has by Remark 4.2 the effect of removing suitable sides of the decagon facets from W φ V + , and thus allows to obtain the pyritohedron as As a simpler situation realizable in Z 3 , and such that the Wulff shape W φ is an octahedron (again not a zonotope, thus not realizable with positive V ), we can take φ as in (4.5) and N = N − ∪ N + where N + = {cyclic permutations of (0, ±1, ±1)} , N − = {±e 1 , ±e 2 , ±e 3 }, and again we take V = V + − V − = c + 1 N + − 1 N − . Note that N + are the nearestneighbors in a face-centered-cubic (FCC) lattice, and the Wulff shape corresponding to V + is a truncated cuboctahedron of sidelength c + 2 √ 2, while the Wulff shape corresponding to V − is a cube of sidelength 2. Following the previous remark as well as Remark 4.2 we find that values c + ∈ 1 4 , 1 2 , give an octahedral Wulff shape.

Γ-limits for signed potentials V
While it is interesting to note that non-zonohedral Wulff shapes can occur, as mentioned in [24] and explained in Section 4.2, another relevant question in the framework of the current paper is whether these shapes can occur as Γ-limits of discrete perimeter functionals like in Theorems 1.1 and 3.3. In this section we show that this is possible but only under additional hypotheses, without which the convergence result can be false. In particular the study at the discrete level presents new difficulties which are invisible at the continuum level. On the other hand, note the following pathological example: say that for N a square, X N = K N ∩ Span Z (N + ), where K N is a square of sidelength √ 2N with edges parallel to the vectors of N + (thus K N is equal to a rescaling of the continuum Wulff shape corresponding to weight φ V for c > 1). In this case, in the "bulk" of K N , all available N − -bonds participate to its V -perimeter. For large N the number of negative contributions is O(N ) while possible positive contributions are O(N 1 2 ), corresponding to bonds from the "boundary layer" of K N . This means that the discrete V -weighted perimeter of X N is −CN + O(N 1 2 ) for some constant C depending on V . This holds even for c > 1 very large, and exhibits an essential difference compared to the continuum perimeter φ V (which as mentioned above, is positive for c > 1). The above phenomenon can be exploited to prove that the Γ-convergence is also false for the above V : if to the X N from the above paragraph we add all the missing points of K N −N 1 ∩ Z 2 for some fixed large N 1 , it is not hard to see that the so-obtained set X N has cardinality N ≤ 2N and discrete V -weighted perimeter bounded below by for large enough N 1 . For the above to hold, we need to choose N 1 depending on V only, and thus we can fix it independently of N . Informally speaking, the negative contributions in a N 1 -neighborhood of ∂K N can always compensate the positive contributions in a √ 2-neighborhood of ∂K N . On the other hand, if N 1 is fixed then where K 0 is a cube of sidelenght √ 2 and sides parallel to vectors from N + , which is a rescaling the Wulff shape of φ V . The φ V -perimeter of K 0 is positive and the rescaled V -weighted perimeters of X N are negative, showing that the Γ-convergence result of Theorem 1.1 is false for all our V with c > 1.
Note that the pathological example (4.7) can be modified to hold in any dimension and in a variety of situations. We formulate the following condition on V , depending on a parameter ǫ > 0: This condition forbids the pathology of Example (4.7), at least for the case of finite N , and it may help for an extension of the Γ-convergence result of Theorem 1.1 to signed V , if satisfied. We leave this investigation for future work.
Note that Example (4.7) produces non-Γ-converging discrete energies mainly based on the structure of N + , N − and not as much on the further degrees of freedom of the weights V + , V − . In the general case this phenomenon is confirmed by the following Proposition 4.8. We recall that if N ⊂ R d is a subset, we say that a set A ⊂ R d is N -connected, if for any x, y ∈ A there exists an N -path in A from x to y, i.e. there exists n ∈ N and points x = x 0 , x 1 , . . . , x n = y ∈ A such that x i − x i−1 ∈ N + for all i. With this definition, we have the following: Proof. The proof is in a similar spirit to the one of Lemma 2.5, so we refer to that proof more details.
To each pair of N -neighbors x, y ∈ L we assign an N + -path in L from x to y, denoted p x,y . To do this, it is sufficient to arbitrarily choose p 0,v for all v ∈ N − , define p 0,v = {0, v} for v ∈ N + , then extend the definition by translation to all pairs of N -neighbors x, x + v. Since N is finite, we find by a simple covering argument that there exists C 1 > 0 depending on our choice of paths p 0,v , v ∈ N − only, such that max v∈N + max z∈L ♯{p x,y : {z, z + v} is a step of p x,y } ≤ C 1 .
If X ⊂ L and we have x ∈ X, x + v / ∈ X, v ∈ N , then p x,x+v has at least one edge {x i−1 , x i } such that x i−1 ∈ X, x i / ∈ X. It follows that Thus if we choose C V = C 1 inf v∈V + V (v) the thesis follows.
Remark 4.9. While the choice of constant C V from Proposition 4.8 which follows from the proof method is in general far from optimal, there exists such C V such that (4.6) holds and the Wulff shape W W for W as in Proposition 4.8, is not a zonotope. For example, part B of Example 4.6 satisfies the hypotheses of Proposition 4.8.
The ensued potential C V V + − V − for very large C V > 1 creates as Wulff shape a truncated octahedron whose hexagonal sides have alternating very long and very short edges, and thus are not centrally symmetric, so that the corresponding W W is not a zonotope.
Remark 4.10. The same considerations as in the rest of this subsection, in particular the ones from Remark 4.9, can be applied to the case of quasicrystals under signed potentials V with similar conclusions, but the condition from Proposition 4.8 can only justify some of the non-zonohedral Wulff shape limits corresponding to the φ V -Wulff shapes examples presented in [24], and not all of them. For a complete theory a sharper condition on V for the extension of Theorems 1.1 and 3.3 to signed V is needed. We plan to address the study of these sharp versions to future work.

Reduction to the symmetric case
In this section we show that for the study of Wulff shapes, it is no restriction to assume that the support function φ is an even function. The case of functionals coming from a potential V : N → R where N ⊂ R d is finite or countable and v∈N |v||V (v)| < ∞ can be reformulated in a form generalizing (4.3), in which instead of the absolute value we consider the positive part of the integrand, and we consider signed finite measures with at most countable support. Indeed, given a potential V as above, we observe that Note that our general class of V includes the case of finite N and constant-sign V from in (1.6). Also note that every φ in the form (4.3) can be put in the form (4.8), just considering the symmetrized measure µ = µ sym Z = 1 2 (µ Z + (−Id) # µ Z ) (which does not change the value in (4.3)).
As we now show, in (4.7) we may always replace " v, x + " by "| v, x |" and consider neighbor sets N = −N and even V only. Indeed, for every finite perimeter set E in R d we have P V (E) = P V sym (E) (4.9) where V sym (v) = 1 2 V (v) + V (−v) . This equality is due to the area formula and to the fact that, for every direction e ∈ S d−1 that we fix, the number of times that the line Re "enters" E is equal to the number of times it "exits" E. Note that (4.9) could fail in the case of relative perimeter functionals P φ (E, Ω), where we compute the perimeter of E relative to a fixed open subset Ω of R d , or when the potential V is not a superposition of terms of type ν, v + . However these cases are not considered here, and the two formulations (4.7) and (4.8) are completely equivalent in our setting. For this reason we will assume from now on that V (and thus φ V ) is symmetric under reflection with respect to the origin. Let µ = (−Id) # µ. Then for every finite perimeter set E ⊂ R d we have Proof. We first prove the result for bounded smooth sets E and then conclude by a density argument. Assume then that E is a bounded, finite perimeter set with smooth boundary. Fix a direction e ∈ S d−1 . We decompose ∂E with respect to the direction e, defining ∂ + e E = {x ∈ ∂E : ν E (x), e > 0} ∂ − e E = {x ∈ ∂E : ν E (x), e < 0} ∂ 0 e E = {x ∈ ∂E : ν E (x), e = 0}.
In particular both N ± e (E, y) are finite for H d−1 -a.e. y ∈ π e ⊥ (E). Now we claim that for H d−1 -a.e. y ∈ π e ⊥ (E) we have N + e (E, y) = N − e (E, y). Indeed, by the area formula we also have π e ⊥ (∂ 0 e E) = 0, which means that for H n−1 -a.e. y ∈ π e ⊥ (E) the line π −1 e ⊥ (y) intersects ∂E always where the normal has nonzero scalar product with e. As observed above, both cardinalities are finite for H d−1 -a.e. y, and by Fubini for H d−1 -a.e. y the 1-dimensional measure of E ∩ π e ⊥ (y) is finite. It follows that E ∩ π −1 e ⊥ (E) is a finite union of bounded intervals, hence we conclude that N + e (E, y) = N − e (E, y) for H d−1 -a.e. y as wanted. We thus obtain This proves the result for smooth, bounded sets.
To prove the general case, let E be a finite perimeter set. Then there exists a sequence E j of bounded, smooth sets such that E j → E in L 1 and P (E j ) → P (E) as j → ∞. By Reshetnyak's theorem [3, Theorem 2.38] the same convergence holds for the perimeters defined by φ µ and φ µ . As a consequence P φµ (E) = lim j→∞ P φµ (E j ) = lim j→∞ P φ µ (E j ) = P φ µ (E).

Final remarks
We collect here some direct generalizations, and some open questions and open directions that would require an essential extra effort compared to the present work, but which could be attacked in future work.

Γ-limits for general signed potentials
The main question that we leave open is to devise a precise, if possible necessary and sufficient, condition on general signed potentials V , so that the Γ-limit result of Theorem 1.1 (or the more general versions of Theorem 1.1 present in [1,8,21] mentioned before) still hold.
The following stability condition on the energy of finite configurations X is a natural hypothesis to impose, when modelling crystalline Wulff shapes: It is natural to study the Γ-limits analogous to our setting, but for general signed V , under the above hypothesis, or under a suitable strengthening of this hypothesis.

Longer range interactions in quasicrystals
There are several possible ways to extend the range of applicability of Theorem 3.3 beyond nearest-neighbor interactions between tiles of the tessellation induced by the multigrid construction of quasicrystals: 1. In the dual space, associate an energy W (y − x) depending on differences of positions of multigrid points x, y belonging to a given fixed lattice, given by any d-ple of hyperplane families from the multigrid construction.
2. Define the interaction W (x, y) between two tiles in the primal space as a function of the sequence of rail directions required to reach y starting from x. For example, in the case of next-nearest-neighbor interactions, let W (x, y) be nonzero only if tile y is a neighbor of x or a neighbor of a neighbor of x. Restricted to neighbors we define W as before, and if there exists z such that x, z belong to a multigrid rail corresponding to multiindex J ′ and z, y are neighbors corresponding to multiindex J ′ , then define W (x, y) as a function depending only on the tuple (J ′ , J ′ ), and independent on the particular choices of x, y.
3. Consider a potential W : R d → [0, +∞) and if x, y are tile centers from the multigrid tiling, define the interaction of the corresponding tiles to be W (y−x).
The first option is the only one in which a direct extension of the setup of Theorem 3.3 will prove a corresponding Γ-convengence result, however we do not develop the theory in this direction because the results would be notationally more involved and we do not have a direct application at hand. The remaining options instead seem better physically motivated by a model of quasicrystalline materials, however require an essential effort beyond the setup of Theorem 3.3 considered here. The comparison between these latter cases seems suited a thorough study in future work.

Curved multigrids
The construction based on merely straight (hyperplane) multigrids that we pursue here can be perturbed, replacing families of hyperplanes by families of curves with richer allowed behavior, which can be well-approximated by hyperplanes at large scale. One then constructs a tiling from the graph of intersections between (d − 1)-ples of surfaces corresponding to different families, by an algorithm similar to Section 3.1.2. The extension of the proof of Theorem 3.3 to this setup seems to us an interesting open direction. On the one hand we mention the use of curved multigrids in modelling (see e.g. [25,Sec. 2.3.2]) and, on the other hand, on the mathematical side, it is interesting to find outthe correct oscillation bounds and metrics for the multigrid hypersurfaces. More precisely, it would be interesting to pursue the determination of natural sufficient regularity conditions on the families of surfaces that we would use, under which extensions of the result of Section 3.1.3 hold true.