Long range order in atomistic models for solids

The emergence of long-range order at low temperatures in atomistic systems with continuous symmetry is a fundamental, yet poorly understood phenomenon in Physics. To address this challenge we study a discrete microscopic model for an elastic crystal with dislocations in three dimensions, previously introduced by Ariza and Ortiz. The model is rich enough to support some realistic features of three-dimensional dislocation theory, most notably grains and the Read-Shockley law for grain boundaries, which we rigorously derive in a simple, explicit, geometry. We analyze the model at positive temperatures, in terms of a Gibbs distribution with energy function given by the Ariza-Ortiz Hamiltonian plus a contribution from the dislocation cores. Our main result is that the model exhibits long range positional order at low temperatures. The proof is based on the tools of discrete exterior calculus, together with cluster expansion techniques.


Introduction
The derivation of the low temperature properties of crystalline solids, starting from a microscopic, atomistic, model, represents a formidable challenge both for theoreticians and practitioners. Realistic atomistic models for solids are characterized by the invariance under the Euclidean symmetries of translations and rotations, which are supposedly broken at low temperatures, as the very existence of crystals in nature witnesses. Unfortunately, from a mathematical point of view, our understanding of the phenomenon of continuous symmetry breaking is still quite limited and, as a consequence, the mathematical theory A. Giuliani of crystalline solids is still in a primitive stage. Even at zero temperature, there are only limited results on the ground state structure of the system: in particular, there are only few, highly simplified, atomistic models for which one can rigorously prove that the ground state is periodic [8,12,14,41]. Even less is known at positive temperatures where most rigorous results are restricted to lattice systems, e.g. [1,20,25]. A notable exception is [5] which establishes the existence of orientational order in a particle system without lattice structure.
Heuristically, we expect that the low energy physics of crystalline materials is dominated by dislocations defects, which interact among each other via an electrostatic-like interaction, and by the formation of grains, which correspond to portions of the crystal with some fixed rotation relative to a background orientation. The grain boundaries are collections of dislocations that are geometrically necessary to connect differently oriented lattices. Remarkably, even though isolated dislocations interact among each other via a Coulomb-like interaction, the energy of a grain appears to scale like the size of its boundary. For a recent mathematical account of this phenomenon, see [31].
There is significant literature on continuum theories for dislocations, see [26] for a starting point. Typically dislocations are represented as closed loops, the energy of a single dislocation loop is proportional to its length, [23]. Discrete dislocation line dynamics represent a very popular simulation technique for studying plasticity since the early 1990s, see e.g. [10] and [27] for a recent account of mathematical results. Continuum models for dislocation configurations have been studied successfully within the framework of Γ-convergence, see e.g. [11,21,22,23]. However, very few results are available on the microscopic derivation of effective continuum theories for dislocations or grain boundaries, see [13,31].
Note that macroscopic effects like plasticity or grain boundary motion are strongly temperature dependent: therefore, it is of particular interest to develop a thermal theory of dislocations, including an equilibrium theory based on the Gibbs distribution.
In this paper, we consider a simple atomistic model for crystalline solids, previously introduced by Ariza and Ortiz [2]. The Ariza-Ortiz model, even if highly simplified, possesses some realistic features expected in real solids, which make it a good starting point for a quantitative understanding of the effects of dislocations and of the formations of grain boundaries. In particular, it has been used to perform discrete dislocation calculations of defects and grain boundaries in graphene, see [3,4,33]. The Ariza-Ortiz model is a discrete model where the interaction energy depends not only on the positions of the particles, but also on the bond structure, see eq.(2.3) below for its precise definition. The model shares some analogies with the Villain model for rotators, in that the energy satisfies an exact additive decomposition property, which allows us to distinguish clearly the elastic ('spin wave') degrees of freedom, and those associated with dislocation defects, see also [30,36,42]. The simplicity of the model allows us to derive sharp estimates on the energy of the grains, on the one hand, and to rigorously characterize key properties of the equilibrium distribution of dislocations at positive temperatures, on the other.
Concerning the kinematics of the Ariza-Ortiz model, we confirm that it supports polycrystalline configuration with energy cost bounded from above by the size of the grain boundary (Theorem 3.2). We also derive sharp asymptotic bounds, albeit in a simpler two-dimensional setting (Theorem 3.3). Our results confirm that the energy density of grain boundaries for small angles is consistent with the Read-Shockley law [38] γ(θ) = θ(c 0 − c 1 log θ) + o(θ), 0 < θ 1, where γ(θ) > 0 is the grain boundary energy density and θ is the orientation difference. See also [31], where the authors establish an upper bound consistent with the Read-Shockley law.
Concerning positive temperatures, we introduce a Gibbs distribution with energy function given by the Ariza-Ortiz Hamiltonian plus a contribution from the dislocation cores. Our main result is that for low temperatures the system exhibits positional long-range order (Theorem 3.1). In particular, this implies that polycrystalline configurations have low probability. To the best knowledge of the authors these are the first rigorous results on dislocations configurations at positive temperature in a microscopic, atomistic, model. See also [6], where similar results have been recently obtained in the context of a related mesoscopic model for crystalline solids. The proof of long-range order is based on the strategy developed in [19,29] for the three-dimensional XY model and other lattice models with Abelian continuous symmetry. The key steps consist in: first, a reduction of the model to an effective model for the dislocation defects, interacting via a tensorial analogue of the electrostatic force; second, a cluster expansion treatment of the latter. The computation of the Green function characterizing the effective interaction among dislocations requires some care, in that the derivation must be compatible with the underlying symmetries of the system, most notably linearized rotational symmetry. This is the key novel feature of the Ariza-Ortiz model, compared to other 'scalar' models treated previously. In this part, we take advantage of the tools of exterior discrete calculus, some aspects of which we briefly review below, for the reader's convenience.
The paper is organized as follows. In Sect.2 we define the Ariza-Ortiz model and discuss its symmetries. In Sect.3 we state our main results, first on the existence of longrange order at positive, low enough, temperatures, then on the energy scaling of grains and grain boundaries. In Sect. 4 we review a few selected aspects of exterior discrete calculus, required in the proofs of our main results. In Sect.5, we prove Theorem 3.1 on long-range positional order. In Sect.6, we prove Theorem 3.3 on the asymptotic computation of the energy of a grain and derive the Read-Shockley law. Finally, in the appendices we collect a few technical results, including the explicit definition of the lattice cellular complex for the face centered cubic lattice, and the asymptotic computation of the correlation decay in the 'spin wave approximation'.

The Ariza-Ortiz model
Let L ⊂ R 3 be the face-centered cubic (FCC) lattice, i.e. L = {n 1 b 1 + n 2 b 2 + n 3 b 3 : n ∈ where: the displacement u : L → R 3 satisfies Dirichlet boundary conditions on Λ, u(x) = 0 if x ∈ Λ; σ : {(x, y) ∈ L 2 : x ∼ y} → L assigns a lattice-valued slip to each nearest neighbor pair (x, y) and it also satisfies Dirichlet boundary conditions on Λ, that is, σ(x, y) = 0 if x, y ∈ Λ. We assume that σ(x, y) = −σ(y, x), so that the energy associated with a nearest neighbor pair (x, y) is independent of the orientation. Moreover, we let x∼y be the sum over the unordered pairs of nearest neighbor sites. The interpretation of σ is that it accounts for crystallographic slip where atoms are being displaced in the direction of the Burger's vector across the slip plane. The deformed configuration is given by the collection of points x + u(x) with x ∈ Λ. The functional H AO accounts for the elastic energy which is caused by the displacement u in the presence of the slip field σ. It should be interpreted as the quadratic approximation of a more complex, non-linear energy. The model has been introduced in [2], which we refer to for details about its microscopic interpretation. See also Appendix A for a heuristic derivation of the model and a discussion about its microscopic meaning.
Note that we study the Ariza-Ortiz model in setting of the FCC lattice because it represents the only simple 3-dimensional lattice, involving only nearest neighbor interactions, satisfying a rigidity estimate a' la Korn, that is, x∼y [(u(x) − u(y)) · (x − y)] 2 x∼y |u(x) − u(y)| 2 , see eq. (5.13) below.

Symmetries of the Ariza-Ortiz model
Consider the infinite volume version of the Ariza-Ortiz energy (2.3), obtained by assuming that u and σ, rather than satisfying Dirichlet boundary conditions on Λ, decay sufficiently fast at infinity so that the infinite sum involved in the definition of the energy makes sense. Such infinite volume Ariza-Ortiz energy is invariant under three different types of symmetry transformations: The presence of the first and third symmetry is a direct consequence of the 'gradient structure' of the Ariza-Ortiz energy, that is, of the fact that it depends on u, σ only upon the combination du − σ. Invariance under linearized rotations is an approximation of the invariance under rotations: . The invariance of the Ariza-Ortiz energy under linearized rotations is a consequence of the observation that (du(x, y) + S (y − x)) · (y − x) = du(x, y) · (y − x) for any skew-symmetric matrix S.
Previously studied models such as the Villain XY model, see, e.g., [18,19], are invariant under the analogues of the first and the third symmetries, but in that context there is no analogue of the second symmetry, which is, instead, a distinctive feature of microscopic models of elasticity. There are significant consequences resulting from the invariance of linearized rotation, most notably the existence of grains, cf. Theorem 3.2. Note that, in a finite box Λ with Dirichlet boundary conditions, the first and second symmetries are broken. On the contrary, the third symmetry is also present in finite volume, provided that v is chosen to satisfy Dirichlet boundary conditions like u. Physically, gauge invariance corresponds to the possibility of conveniently re-labelling the atoms and, correspondingly, of re-defining the nearest neighbours, without any energy cost. Mathematically, gauge invariance implies that the energy only depends on the dislocation part of σ, defined in the following section.
tions on Λ is gauge invariant in the sense that for each v : L → L that satisfies Dirichlet boundary conditions on Λ. To remove this degeneracy we say that two slip fields σ and σ are equivalent if dσ = dσ , with The function q = dσ is called the dislocation part of σ. Note that, if σ satisfies Dirichlet boundary conditions, then also q does, i.e., q(x 1 , x 2 , x 3 ) = 0 for x 1 , x 2 , x 3 ∈ Λ c . A discussion of the link between slip fields without dislocations (dσ = 0) and the existence of v : L → R 3 such that dv = σ can be found in Sect. 4.
The field dσ assigns to each triangular face f , identified with a 3-cycle of nearest neighbor sites, a current flowing orthogonally to f , in the direction induced by the orientation of f . Typically q = dσ is decomposed into a sum of dislocation lines, i.e., q = j q j , where the supports of the q j are the maximal connected components of supp q. Each of these q j can be thought of as a current loop. It will be shown in Sect. 4 that dq = 0, where dq is the discrete analogue of the curl of q: it is a function defined on the elementary cells of L that, on each cell, equals the sum of the values of q on the faces of the cell, with the appropriate orientation. In terms of the current loop representation of q, this curl-free condition means that the current loops are closed.
Denoting by S the set of representatives of non-equivalent slip-fields satisfying Dirichlet boundary conditions (i.e., vanishing on edges contained in Λ c ), we are now in a position to define the expectation of a gauge-invariant observable ϕ (i.e., ϕ(u, σ) = ϕ(u+v, σ+dv) for any v : L → L supported in Λ) with respect to the Boltzmann-Gibbs distribution by and the integral runs over R 3|Λ| (recall that u(x) ∈ R 3 , for x ∈ L, and u(x) ≡ 0 if x ∈ Λ c ). The function W represents the energy contribution of the dislocation cores and has the form where w is even and the sum runs over the unordered set of faces of L (i.e., the set of 3-cycles of nearest neighbor sites in L, modulo their orientation). We assume that for some positive constant w 0 . We remark that the purely additive structure of the core energy, Eq.(3.2), is assumed here just for simplicity: our proofs could be adapted to the case of correlated energies, provided their correlation decays to zero sufficiently fast at large distances, but we prefer to stick to the assumption of exact additivity here, in order to keep technicalities to a minimum. Note that the condition that ϕ, like H AO , is invariant under gauge transformations is a natural requirement: essentially, we are saying that slip fields differing by exact forms are physically un-distinguishable.
We will be specifically interested in the following observables: for x, y ∈ L and v 0 ∈ L * (the dual of L, whose basis vectors m 1 , m 2 , m 3 are defined by the conditions b i · m j = 2πδ i,j , see (C.1)), we let and we let be the corresponding two-point observable. It is apparent that both ϕ v 0 ;x (u) and ϕ v 0 ;x,y (u) are gauge invariant, thanks to the condition that v 0 ∈ L * . The one-point observable ϕ v 0 ;x is appropriate for testing the breaking of translational symmetry (i.e., of symmetry 1 in Sect.2.1) in the presence of Dirichlet boundary conditions: in fact, it is peaked at u(x) = 0 mod L (in particular, it is not invariant under translations u(x) → u(x) + τ ), and it has zero average under translations of u(x). Similarly, the corresponding two-point observable ϕ v 0 ;x,y is appropriate for testing the existence of positional long range-order. We define the expectations We are interested in taking the thermodynamic limit Λ → L that, for boxes Λ = Λ (N ) like in (2.2), simply indicates the limit N → ∞.
There are positive constants C, β 0 , r 0 , which do not depend on x, y and β such that, if β > β 0 and |x − y| > r 0 , Eqs. (3.4) establish the existence of translational symmetry breaking and long range positional order in the three-dimensional setting: in particular, the first equation implies that, for β large enough, the weak limit as Λ → L of the Gibbs state E β,Λ breaks the translational symmetry u(x) → u(x) + τ . Conversely, at small enough β, it is expected that the limiting Gibbs state is invariant under such translational symmetry, that c β,Λ (v 0 ; x) decays exponentially to zero in the distance dist(x, Λ c ) as Λ → L, and that lim inf Λ→L c β,Λ (v 0 ; x, y) decays exponentially to zero as |x − y| → ∞, because of the screening phenomenon [9], but this remains to be proved for the Ariza-Ortiz model. The limiting value lim inf Λ→L c β,Λ (v 0 ; x), which passes from being positive at large β to being (at least conjecturally) identically zero at low β, has the interpretation of order parameter for positional order.
The reason why we write lim inf Λ→L rather than lim Λ→L in (3.4) is that a priori we do not know whether the limit exists: our system is of Coulomb-type and the standard theory of the existence of the thermodynamic limit does not apply directly. There are several results in the literature about the existence of the thermodynamic limit of Coulomb systems in 3D, but they do not apply literally to our case, see e.g. [16,32] and the review [9] and references therein. It is likely that they could be adapted to our context as well, but this is beyond the scope of our paper. Theorem 3.1 is a consequence of Theorem 5.1 which is stated with proof in Section 5. Theorem 5.1 and its proof provide a more detailed estimate than (3.4): in particular, they show that both c β,Λ (v 0 ; x) and c β,Λ (v 0 ; x, y) factor exactly into the product of two contributions, one associated with a Gaussian average (the 'spin wave contribution') and one associated with an effective theory for the dislocation cores. The first contribution is explicit, and asymptotically equal, as Λ → L and |x − y| → ∞, to e −C 0 /(2β) resp. e −C 0 /β for the one-point resp. two-point observable, for an explicit constant C 0 . The second is bounded via cluster expansion and the use of Jensen's inequality, following the same strategy of [19,29], and leads to an exponentially small correction to the spin-wave contribution. Note that the assumption d = 3 plays a key role both in the computation of the spin wave contribution (Section 5.1) and in the estimate of the correction due to the dislocation cores (Section 5.2). Theorem 3.1 is analogous to [6,Theorem 3] that, however, refers to long-range orientational order in a mesoscopic model for a solid with dislocations. An important difference between our setting and the one in [6] concerns the modeling part. While our model, even though simplified, has a direct microscopic interpretation, theirs involves an auxiliary set of currents, whose microscopic interpretation is not immediate. It is likely that the model in [6] could be obtained starting from a more fundamental atomistic one, via a suitable coarse graining procedure. It would be very interesting to substantiate this expectation by rigorous results. From a technical point of view, the tensorial structure of the Ariza-Ortiz Hamiltonian on the FCC lattice introduces some extra difficulties, compared to [6], in the reduction to an effective model of dislocations and in the treatment thereof, which we solve thanks to the tools of discrete exterior calculus, reviewed below. On the other hand, the general strategy of our proof is analogous to that in [6], in that both rely on the ideas of [19,29].
If the dimension is one or two then the existence of long-range positional order is prevented by the Mermin-Wagner theorem [34,35], see also [17,28,37,39]. However, the Mermin-Wagner theorem does not prevent the possibility of having orientational order in two dimensions; actually, spin wave theory suggest that orientational order should be present in two dimensions [34]. It is not a priori clear, not even heuristically or intuitively, whether the presence of dislocations, and in particular of grains, can destroy the prediction based on spin wave theory. Therefore, it would be extremely interesting to prove or disprove the existence of long-range orientational order in a concrete atomistic model for a two-dimensional elastic crystal with dislocation. Probably, the simplest such model is the analogue of the model studied in this paper, in a two-dimensional setting (e.g., in the case that the 3D FCC lattice is replaced by the 2D triangular lattice). We expect that the methods developed in [18] for the study of the Kosterlitz-Thouless transition in the 2D Villain rotator model may be adapted to such a case. We plan to come back to this problem in a future publication.

Energy scaling of grains
An important question which cuts to the core of the crystal problem is whether the Ariza-Ortiz model accounts for structures such as grain-boundaries. Grains only exist because of the additional symmetry of atomistic systems: (linearized) rotational symmetry, sometimes also referred to as 'objectivity'. More precisely, let S ∈ R 3×3 be a skew-symmetric matrix and G ⊂ L be the location of the grain, which we assume to be simply connected and bounded. We say that a pair (u, σ), with u : L → R 3 a displacement field and σ : E 1 → L a lattice-valued slip field (here E 1 is the set of nearest neighbor pairs of L and σ(x, y) is assumed to be odd under orientation flip (x, y) → (y, x)) supports a 'perfect grain' G with orientation S if it is gauge equivalent to a configuration (u , σ ) such that We say that (u , σ ) is gauge equivalent to (u, σ), if (u , σ ) = (u + v, σ + dv) for some lattice valued function v, see Sect.2.1. Note that (3.5)-(3.6) do not impose any constraint on the nearest neighbour bonds (x, y) such that x ∈ G and y ∈ G c , or viceversa. For later reference, we denote this set of bonds by E b 1 (G) ('b' for 'boundary'): and we let |E b 1 (G)| denote the number of elements of E b 1 (G) modulo orientation. It is not obvious from the outset whether a pair (u, σ) supporting a perfect grain G with orientation S can be chosen such that the associated energy is smaller than the volume of G. For example, the pair (u S , 0), with clearly supports a perfect grain G with orientation S, for any fixed τ ∈ R 3 . However, it can be easily checked that if |G| is large. The fact that the energy H AO (u S , 0) is larger than |G| is a consequence of the discontinuity across the boundary of G. However, one should not conclude from this that the 'optimal grain energy' scales more than extensively: on the contrary, it is remarkable that, for large grains, it scales proportionally to the size of the grain boundary E b 1 (G), as summarized in Theorem 3.2 below. Of course, before formulating our result, we first need to clarify what we mean by 'optimal grain energy'. Let us remark that there is no unique, well-established, notion of grain energy at a microscopic level. In fact, such a notion is well-defined at the mesoscopic level, in which case it is related to the distribution of dislocations at the boundary of the grain, and is known to scale proportionally to the grain boundary, see e.g. [40]. Here we propose a microscopic definition thereof, which we expect to reduce to the usual mesoscopic definition in an appropriate scaling limit. A rigorous connection of our microscopic definition with the continuum one is an interesting open problem, which goes beyond the purposes of this paper.
We let the optimal grain energy be defined as is the set of -minimizers of H AO (u S , σ), with u S as in (3.7), over the slip fields compatible with the grain G: more explicitly, M ( ) S (G) is the set of latticevalued slip fields σ such that there exists τ (see (3.7)) for which σ, τ realize the infimum of inf τ inf * σ H AO (u S , σ) within a precision , where the * on inf * σ indicates the constraint supp σ ⊆ E b 1 (G). We are ready to state the basic bound on E G (S), showing, as anticipated above, that the optimal grain energy scales like the grain boundary for large grains.
Theorem 3.2. For any skew symmetric matrix S ∈ R 3×3 and any bounded connected set G, It is very likely that our upper bound can be improved, i.e., it is not sharp. In fact, the construction of matching upper and lower bounds in the limit of a large grain constitutes an interesting mathematical problem.
Proof. In order to prove the theorem, we will construct a lattice valued slip field σ S sup- .7), with τ = 0. We first determine slip amplitudes ξ (l,n) ∈ R such that the matrix S can be decomposed into simple slip systems, i.e.
with the convention that b l ∈ L are the slip vectors and m n ∈ L * are the slip plane normals. The standard 12 slip systems of the FCC lattice are Once the slip amplitudes ξ (l,n) are fixed, we let and σ S (x, y) = 0 otherwise. Let us now compute H AO (u S , σ S ). We partition the set E 1 into three groups: is the set of bonds outside G. The partition of E 1 induces a decomposition of the energy: where we used the fact that σ S is zero on  Colored triangles indicate the support of dσ, with σ = σ S + dv S . The right panel shows the relaxed displacement field u σ which minimizes H AO (·, σ). Now, the difference in parentheses in the right side is between 0 and 1. Therefore, recalling that |x − y| = |b l | = 1 and that the sum over (l, n) runs over 12 different terms, we find as desired.
In order to visualize the 'optimal' location of the atoms within a grain, we remark that the pair (u S , σ S ) used in the proof of Theorem 3.2 is gauge equivalent to a configuration (u, σ) such that: (1) |u(x)| ≤ 6 for x ∈ G, and u(x) = 0 otherwise, (2) the support of σ is contained in E i 1 (G). In order to exhibit such an equivalent pair, we let u = u S + v S and σ = σ S + dv S , with (u S , σ S ) the same as those used in the proof of the theorem, and We visualise in left panel of Fig. 1 such a displacement field u in a two-dimensional setting where S = 1 and Neumann boundary conditions are used (see Sect.4.1 for a definition of Neumann boundary conditions). The colored triangles are the support of dσ. The minimizer of H AO (·, σ) is shown in the right panel of Fig. 1. The corresponding minimal energy, inf u H AO (u, σ), is the one that the system will reach after relaxation at fixed slip field σ. In the limit of large grain, it is supposed to provide a good approximation for the optimal grain energy E G (S) in (3.8). As it will be proved in the following sections, remarkably, the minimal energy inf u H AO (u, σ) only depends on the 'charge distribution' q = dσ, which, therefore, characterizes the grain from an energetic point of view.

Read-Shockley law
Theorem 3.2 shows that the optimal energy of a perfect grain scales like its boundary, but does not provide an explicit formula for the surface tension, that is, the proportionality constant in front of |E b 1 (G)|, in the limit of a large grain. Physically, there are explicit expectations for the surface tension, specifically in the limit of small rotation angles: according to the Read-Shockley formula [38], given a large grain, rotated by a small angle θ with respect to a reference crystalline background, its total energy is proportional to its boundary, with a proportionality constant γ(θ) of the form (1.1). An upper bound which is consistent with the logarithmic scaling can be found in [31].
In this section, we state two results about the exact, asymptotic, computation of the energy of a dislocation dipole and of two walls of dislocations with opposite charges, far away from each other. In particular, the energy of the two parallel walls of dislocations with opposite charges is expected to correspond to the optimal energy of a grain supported in the region between the two walls (the electrostatic analogue to keep in mind is a capacitor: dislocations correspond to the charges on the plates of the capacitor, and the intermediate region between the plates is where the elastic energy concentrates), in the sense of definition (3.8). The reader can convince herself/himself that the smaller the density of dislocations on the walls, the smaller the rotation angle of the grain, and that in the limit of small density of dislocations, the rotation angle goes to zero linearly with the density. Therefore, the computation of the energy of the 'dislocation capacitor' performed below provides information on the optimal energy of the corresponding grain. Our main result is that we recover the Read-Shockley law for a grain with such a simple, specific, geometrical shape.
The computations are reported in Sect.6. For simplicity, we perform the computations in two dimensions, but similar results, including the logarithmic dependence of the surface tension on the rotation angle, in the sense of (1.1), can be extended to three dimensions, by assuming that the distribution of dislocations under consideration is translationally invariant in the third coordinate direction; however, in three dimensions the computations become cumbersome and their key features would be hidden behind unimportant technical complications: therefore, we prefer to restrict to 2D and leave the tedious but straightforward extension to higher dimensions to the interested reader.
We denote by T the triangular lattice and, with some abuse of notation, we let its basis . Given a finite box Λ ⊂ T of side N (the 2D analogue of (2.2)), we let the 2D Ariza-Ortiz energy in Λ with Dirichlet boundary conditions be defined by the same formula (2.3); with some abuse of notation, we denote the 2D energy by the same symbol H AO (u, σ). Figure 2: Graphical illustration of the charge distribution q n dip , for n = 6. The shaded triangles, corresponding to faces f 0 and f n , indicate the support of q n dip . We also show in red the support of a slip field σ n dip such that dσ n dip = q n dip , see (6.3). The sites labelled j = 1, . . . , n on the bottom (resp. top) row have coordinates jb 1 (resp. jb We consider a dislocation dipole formed by a pair of opposite charges ±b 1 , separated by a distance n in direction b 1 , whose 'charge distribution' is: We also consider two parallel arrays of dislocations, formed by M dislocation dipoles as in (3.12), arranged one at a distance m √ 3 from the other in the direction orthogonal to b 1 , whose charge distribution is: In the limit M → ∞, the charge distribution q M,n,m grain tends to that of two infinite walls of dislocations, separated by a distance n, with charge density ∼ 1/m. As discussed in Appendix B, its energy is expected to coincide at dominant order with the optimal energy of a grain supported in the region between the walls, rotated by an angle θ ∼ 1/m, in the limit m → ∞.
be the energy of a dipole and the energy density of a grain boundary per unit length, in the thermodynamic limit. Then The proof of Theorem 3.3 is given in Sect.6. Eq.(3.17) is the desired Read-Shockley law for the energy of a grain boundary. Its remarkable feature is that it is asymptotically independent of the separation among the two arrays of charges it consists of. This is in sharp contrast with the 'capacitor law', i.e., with the formula for the energy of two parallel arrays of 'scalar' dipoles, i.e., of a similar arrangement of charges in the usual Coulomb lattice gas, which scales linearly in n at large separation n. For a technical comparison of the computations leading to the Read-Shockley and the capacitor laws, see Sect.6.2.1 below.
Note that the E grain (n, m) does not include a contribution from the dislocation cores. Of course, the inclusion of such a contribution, of the form W (q), see (3.2), can be done without any additional difficulty. Note that the extra energy from the dislocation cores would contribute O(1/m) to the right side of (3.17) and, therefore, would not modify the dominant asymptotics of the Read-Shockley law.

Exterior calculus
In this section, we review a few basic aspects of discrete exterior calculus, which is a fundamental tool used in the proof of the main results. In particular, an application of the Hodge decomposition to the Ariza-Ortiz model will allow us to decompose its energy in the sum of a 'spin wave' part plus a 'dislocation' part: such a decomposition is central to our analysis and will be used systematically in the following.

Cellular complex, discrete p-forms and discrete differential
The domain of the three-dimensional Ariza-Ortiz model is given by cells consisting of • oriented faces E 2 (polygons whose sides are consistently oriented edges) • oriented volumes E 3 (polyhedra whose faces are consistently oriented faces), which form a cellular complex, cf [24]. The orientation of a face f ∈ E 2 is defined by the direction of a reference vector, orthogonal to f ; the sides of f are said to be consistently oriented if their orientation satisfies the 'right-hand rule'. The orientation of a volume v ∈ E 3 is either 'outward' or 'inward': its faces are said to be consistently oriented if the directions of their reference vectors all point, correspondingly, in the outward or inward direction. The case that is of interest to us is where the vertices coincide with a Bravais lattice, a case that is commonly referred to as lattice cellular complex. We are specifically interested in the case that E 0 = L, with L the face centered cubic lattice, in which case we let, in particular, E 1 be the set of all ordered pairs of 'nearest neighbor' sites (those at smallest Euclidean distance), and E 2 the set of oriented triangular faces associated with the 3-cycles of nearest neighbor sites. A detailed description of the corresponding cellular complex is given in Appendix C.
The boundary operator ∂ p : E p → E p−1 with p > 0 returns the set of boundary cells with the appropriate orientation. By repeated applications of the boundary operator, any p-cell c with p > 0 is mapped to a set of vertices in E 0 , which we refer to as the 'set of vertices of c' and denote by V (c). We only require a small subset of cohomology theory and will use a minimalistic setup. In particular, the action of ∂ p is defined via explicit formulae in Appendix C, the reader is encouraged to confirm that it coincides with the standard definition [24, Sec 3].
The vector space C p is the set of p-forms, namely the set of functions u : E p → R 3 that are odd under orientation flip. The lattice-valued p-forms, that is, those that return values in L, will be denoted by C p L . We define for p = 0, 1, 2 the exterior derivative operators d p : C p → C p+1 . If p = 1, 2, they are given by the formula if p = 0 and e = (x, y) ∈ E 0 is an oriented edge then A straightforward calculation shows that d p+1 d p = 0, for p = 0, 1, see, e.g., [24, Lemma 2.1]. In some cases, it is useful to interpret d p+1 d p as being = 0 also for p = 2, in which case we let d 3 := 0. Whenever the notation is un-ambiguous, we will drop the label p from d p (i.e., if it is clear from the context that u is a p-form, then we will write du instead of d p u).
We are interested in the cellular complexes and the corresponding set of p-forms, obtained by taking finite portions Λ of L, with prescribed boundary conditions, namely Dirichlet, Neumann, or periodic. For simplicity, we restrict to cases in which such finite portions are parallelepipeds of size N , like in (2.2).
In the case of Neumann boundary conditions, we let Λ p be the subset of E p consisting of the p-cells c, whose set of vertices are contained in Λ; the p-forms of interest are those that depend only on the p-cells in Λ p . In the case of Dirichlet boundary conditions or periodic boundary conditions we maintain the same cellular complex as for L. For Dirichlet boundary conditions the relevant p-forms are those that assume non-zero values on cells whose vertices have non-empty intersection with Λ. For periodic boundary conditions the p-forms of interest are N -periodic in the directions b 1 , b 2 , b 3 . In all these cases, with some abuse of notation, we denote the cellular complex by (Λ 0 , Λ 1 , Λ 2 , Λ 3 ) and by C p the corresponding sets of p-forms.
For any given finite Λ as in (2.2) and all the three boundary conditions introduced above, the vector spaces C p are finite dimensional Hilbert spaces with canonical inner product 3 · , · . Thanks to the relation d p d p−1 = 0, for p = 1, 2, 3, one has that range d p−1 ⊂ null d p and we can define the cohomology groups with the conventions that: / denotes the standard quotient operator, null d 3 = C 3 , and H 0 = null d 0 . As usual, we say that: • if u ∈ C p has the property that du = 0 (i.e., if u ∈ null d p ), then u is closed; • if u ∈ C p has the property that u = dv for some v ∈ C p−1 (i.e., if u ∈ range d p−1 ), then u is exact.
In terms of these definitions, H p is the subspace of closed p-forms modulo the exact pforms (i.e., modulo the following equivalence relation for closed p-forms: , then any closed v is automatically exact. More generally, v ∈ C p is exact if and only if it is closed and additionally satisfies dim H p linear constraints. The cohomology groups associated with the box Λ with Dirichlet, Neumann and periodic boundary conditions are known, and are the following.

Periodic boundary conditions
(4.5) 3 We use the convention that, for where we recall that, for p > 0, Λ p is the set of oriented p-cells: therefore, the factor 1/2 in front of the sum is chosen so that every unoriented p-cell is effectively counted just once.
In order to prove these formulas, note that the box Λ is topologically equivalent to a 3-dimensional ball B 3 ⊂ R 3 . Therefore, the cohomology for Neumann boundary conditions corresponds to the de Rham cohomology of the ball, which is given by the Poincaré Lemma [7]. The cohomology with Dirichlet boundary conditions corresponds to the cohomology with compact support for the ball B 3 ⊂ R 3 ; the result then follows from Poincaré duality between de Rham cohomology and cohomology with compact support [7]. Finally, the cohomology with periodic boundary conditions is the de Rham cohomology of a three-dimensional torus T 3 ∼ S 1 × S 1 × S 1 ; the result is then an application of Künneth formula [7].
In the following, we will also need a quantitative version of the Poincaré Lemma for lattice valued 2-forms, in the form stated next.

Hodge decomposition
In this section, we obtain a representation of H AO (u, σ) in terms of dσ. Setting q = dσ we can choose a representative displacement-and slip-field (u σ , σ q ) ∈ C 0 × C 1 such that dσ q = q and du σ = σ − σ q . With such a choice Interestingly it is possible to choose (u σ , σ q ) so that the energy decomposes into a purely elastic part and a dislocation part cf Theorem 4.3 below. In general σ q is not L-valued, which means that a physical interpretation is not obvious. While the additive decomposition simplifies our analysis significantly, we believe that it is not central for the validity of our main results.
To establish the additive decomposition of the Ariza-Ortiz energy we employ the classic Hodge decomposition. The fundamental idea is to construct the relevant u q and σ q in terms of solutions of the Poisson equation.
The Laplace operator ∆ p : C p → C p is defined by where d * p−1 : C p → C p−1 is the adjoint of d p−1 , with respect to ·, · (for p = 0, (4.6) should be interpreted as ∆ 0 = d * 0 d 0 , i.e., d −1 = d * −1 := 0). Also in this case, as in Sect.4, we will drop the dependence upon p whenever the space, which ∆ p or d * p−1 act on, is clear from the context. Note that d and d * both commute with ∆ (that is, Proof. The invertibility of the Laplacian is an immediate consequence of the classical Hodge decomposition, which we prove next: In order to prove this decomposition, we first demonstrate that null d p and null d * p−1 are orthogonal. Let u, v ∈ C p be such that du = 0 and d * v = 0. Since H p = {0}, there exists w ∈ C p−1 such that u = dw. Therefore as desired. Next, we demonstrate that Suppose that: (i) ϕ, u = 0, ∀u ∈ null d p , that is, for all p-forms u such that u = dw for some w ∈ C p−1 ; (ii) ϕ, v = 0, ∀v ∈ null d * p−1 , in particular for all the p-forms such that v = d * z, for some z ∈ C p+1 . By using (i), ϕ, dw = d * ϕ, w = 0, ∀w ∈ C p−1 , that is, d * ϕ = 0. Moreover, by using (ii), ϕ, d * z = dϕ, z = 0, ∀z ∈ C p+1 , that is, dϕ = 0 ⇒ ϕ = dψ, for some ψ ∈ C p−1 . In conclusion, ϕ, ϕ = dψ, ϕ = ψ, d * ϕ = 0, as desired (in the last step we used that d * ϕ = 0). This concludes the proof of (4.7).
We are now in position of proving the invertibility of ∆. We first prove injectivity: assume that ∆u = 0, from which that is du = 0 and d * u = 0. In view of (4.7), this implies that u = 0 and, therefore, ∆ is injective. Next, we prove surjectivity: assume that u ∈ (range ∆) ⊥ , i.e. u, ∆v = 0 for all v. Then ∆u, v = 0 for all v, that is ∆u = 0, which implies u = 0, as we already saw.
The condition H p = {0} motivates the use of Dirichlet or Neumann boundary conditions, in which cases H 1 = H 2 = {0}, so that the Laplacian acting on 1and 2-forms are invertible.
With this notation the Ariza-Ortiz energy can be written as a functional where B ∈ Lin(C 1 , C 1 ) is the linear operator such that, for any v ∈ C 1 , (Bv)(e) = (v(e) · δe) δe, and, given e = (x, y), we denoted δe := y − x. We are now in a position to establish the decomposition of the Ariza-Ortiz energy into elastic and dislocation part. with Dirichlet or Neumann boundary conditions. Assume that q ∈ C 2 L satisfies dq = 0. For any σ ∈ C 1 L with the property dσ = q, the Ariza-Ortiz energy admits the additive decomposition

9)
where σ q is the minimizer of v → v, Bv on C 1 (rather than on C 1 L ) subject to the constraint that dv = q, and u σ is defined as (4.10) If we take Dirichlet boundary conditions, then the minimizer σ q is given by σ q = Gq with where A : C 0 → C 0 is the invertible operator A := d * Bd.
For later reference, we note that the adjoint of G can be explicitly written as As an immediate consequence from (4.9) one obtains the equation Proof. Equation (4.10) implies that du σ = dd * ∆ −1 (σ−σ q ). Moreover, thanks to the commutation relation d∆ = ∆d, and recalling that dσ = dσ q , we also find that Combining these two identities, we find (4.14) Furthermore, as σ q is the solution of a constrained minimization problem there exists a Lagrange parameter λ ∈ C 2 such that σ q satisfies the Euler-Lagrange equation Bσ q = d * λ. Hence, thanks to d * d * = 0, one obtains Therefore, which establishes (4.9).

Proof of Theorem 3.1
Our goal is to compute a lower bound on in the two cases that ϕ(u) = ϕ v 0 ;x (u) = cos(u(x) · v 0 ) and ϕ(u) = ϕ v 0 ;x,y (u) = cos((u(x) − u(y)) · v 0 ). These functions can be conveniently rewritten in terms of the functions g, h andg,h, defined as follows.
• d * h = g and d * h =g, or equivalently h, du = g, u and h , du = g, u .
To show that the equation d * h = g actually admits a solution we can consider a pairwise disjoint collection of oriented edges (e i ) i=1...n ⊂ E 1 that form an oriented connected path P x→xext from x to x ext , where x ext is some vertex outside Λ, which we fix once and for all at a distance 1 from Λ. With such a path we can define where in the second lineē is the edge with the same vertices as e but opposite orientation. Analogously, we leth = h v 0 ;x,y and note that, with this choice, d * h =g. In terms of these definitions, for any σ ∈ C 1 L , ϕ v 0 ;x,y (u) = cos g, u = cos h , du = cos h , du − σ , because both h, σ and h , σ are equal to 0 mod 2π. Thanks to the decomposition (4.9), in the cases of interest (5.1) can be rewritten as 3) where q = dσ, and similarly for E β,Λ (ϕ v 0 ;x,y ), with h replaced byh. Recalling (4.14), we can rewrite du − σ = d(u − u σ ) − σ q , with σ q = Gq, so that, renaming u − u σ ≡ u , 4) where C 2 * = {q ∈ C 2 L : dq = 0} is the set of closed, lattice-valued, 2-forms satisfying Dirichlet boundary conditions on Λ; again, E β,Λ (ϕ v 0 ;x,y ) admits an analogous representation, with h replaced byh. Note that the probability measure in the right side of (5.4) is factorized: it is the product of a Gaussian measure P sw β,Λ on u (the spin wave part of the measure) times a discrete measure P dis β,Λ on the dislocation cores q. This factorization property is due to the quadratic nature of the Ariza-Ortiz model, and makes our statistical mechanics version of the Ariza-Ortiz model reminiscent of the Villain model for classical rotators. Of course, the partition function inherits the same factorization property: Plugging this representation in (5.4) and in its analogue for E β,Λ (ϕ v 0 ;x,y ), and noting that P sw β,Λ and P dis β,Λ are even, we find where h = h v 0 ;x,xext andh = h v 0 ;x,y ; moreover, The rest of the section is devoted to the proofs of (5.7) and (5.8).

5.1
The spin wave contribution to the two-point function (proof of (5.7)) Recalling the definitions A = d * Bd and g = d * h= 1 x · v 0 , we find that the spin wave contribution to the one-point function is As proved in Appendix D, the thermodynamic limit of the right side can be explicitly written in Fourier space: where B = {ξ 1 m 1 + ξ 2 m 2 + ξ 3 m 3 : ξ i ∈ [0, 1)} is the Brillouin zone (recall that m 1 , m 2 , m 3 are the basis vectors of L * , see Appendix C), |B| is its volume, and In Appendix D we prove thatÂ(k) is singular iff k ∈ L * and that, if k is close to 0, where the positive constant c 0 can be chosen, e.g., to c 0 = (3 − √ 5)/4. We remark that this bound depends critically on the structure of the underlying lattice: changing FCC into cubic does not preserve the property thatÂ(k) behaves qualitatively like the Laplacian k 2 at low momenta. The inverse operator reads: so that the integral in the right side of (5.10) is convergent and we can rewrite Note that the fact that C 0 is finite crucially depends on the fact that we are in three (or, better, in more than two) dimensions: in two dimensions its analogue would be logarithmically divergent. This concludes the proof of the first of (5.7). The proof of the second of (5.7) is analogous: a repetition of the previous computation implies that which leads to the second of (5.7), thanks to the fact that for a suitable constant c 1 > 0, see Appendix D for details.

The dislocation contribution to the two-point function (proof of (5.8))
We now want to bound from below the dislocation contribution to the one-point function, namely and similarly for the two-point function. Recall that the probability weight e −βW (q) is of factorized form, e −βW (q) = f e −βw(q(f )) , where the product runs over the faces that have non zero intersection with Λ and w(q(f )) ≥ w 0 |q(f )| 2 for some positive w 0 . Note that this weight can be equivalently rewritten as The goal is to find a lower bound on Z dis β,Λ (h)/Z dis β,Λ (0), of lower order than the spin wave contribution. We now perform a sine-Gordon transformation: we introduce the Gaussian measure µ β (dφ) with covariance β (G * BG + w 0 1)> 0, so that e − β 2 q,(G * BG+w 0 1)q = µ β (dφ)e i q,φ , Remark 5.2. The presence of the w 0 1 term in the covariance β (G * BG + w 0 1) is crucial for making it positive definite, rather than just non-negative. In fact, G * BG has null directions. To see this, take σ ∈ C 1 L to be +b 4 resp. −b 4 on the edge (0, b 1 ) resp. (b 1 , 0), and zero otherwise, so that Bσ = 0 but q := dσ = 0. With these definitions, it is easy to check that (q, G * BGq) = (σ q , Bσ q ) = 0: in fact, σ q is the minimizer of (v, Bv) = Bv 2 among the v ∈ C 1 such that dv = q; taking v = σ, we have Bv = 0, so Bσ q = 0.
We now perform the sum over q in two steps: we first fix the support of q and then sum over the charge configurations compatible with that support: For short, we shall write K(X, φ+G * h) = K(X). Note that K(X 1 ∪X 2 ) = K(X 1 )K(X 2 ) if X 1 , X 2 ⊂ Λ 2 are disconnected, i.e., if there is no 3-cell whose boundary contains both an element of X 1 and an element of X 2 ; the key thing to observe is that in such a situation, if we let q = q 1 + q 2 with supp (q 1 ) = X 1 and supp (q 2 ) = X 2 , the constraint d(q 1 + q 2 ) = 0 'factorizes' in dq 1 = dq 2 = 0. Note that this factorization into locally neutral contributions is another point where the condition that we are in more than two dimensions enters crucially. Given X, we let X 1 , . . . , X n be its maximally connected components and note that The upper bound on λ implies that We use the factorization property of K to rewrite Z dis β,Λ (h) as where δ(X 1 , . . . , X n ) = 1≤i<j≤n δ(X i , X j ), and δ(X, Y ) = 1 if X and Y are disconnected, and = 0 otherwise. As well known, see, e.g., [ κ(X 1 ) · · · κ(X n )δ(X 1 , . . . , X n ) = = exp n≥1 X 1 ,...,Xn⊆Λ 2 κ(X 1 ) · · · κ(X n )ϕ(X 1 , . . . , X n ) , where ϕ(X 1 , . . . , X n ) is the Ursell function: if G n is the complete graph on the vertex set {1, . . . , n}, for n > 1, while ϕ(X 1 ) = 1, for n = 1. Note that κ(X 1 ) · · · κ(X n )ϕ(X 1 , . . . , X n ) is non-zero only if Y = X 1 ∪ · · · ∪ X n is connected. The sums in (5.19) are absolutely convergent, uniformly in Λ, provided that there exists a positive function a(X), independent of Λ, such that, for any fixed, connected, non-empty X * ⊆ Λ 2 , see [15,Theorem 5.4]. In our case, if β is sufficiently large, thanks to the upper bound on κ(X), eq.(5.18), we can choose a(X) = e −βw 0 /4 |X|. We now insert the definition of κ in (5.19) and rewrite it as (5.20) where in the right side X i := supp (q i ). Using the fact that λ(x) is exponentially small, as well as the fact that the Ursell function decays exponentially to zero at large distances, we get that, for β large enough, see Appendix E for a proof. Note also that z(β, q) is zero unless q has connected support.
We now apply Jensen's inequality, i.e., m β (dφ) exp (·) ≥ exp m β (dφ)(·) , and find (noting that m β (dφ) sin q, φ = 0 and | m β (dφ) cos q, φ | ≤ 1), We now need to manipulate q, G * h . Proposition 4.1 implies that there exists an L-valued 1-form n q such that dn q = q. Recall that 1. The support of n q is contained in B(q), the smallest parallelepiped containing the support of q, 2. The maximum of |n q | is bounded in terms of the 2-norm of q, as follows: n q ∞ ≤ c q 4 2 for some positive c.
Therefore, recalling the definition of G * , see (4.12), we get Now the first term in the right side is an integer multiple of 2π, and can be dropped for the purpose of computing the cosine. The last term, by using that d * commutes with ∆, which is zero, simply because d * BdA −1 = 1. We are left with the second term, which can be rewritten in terms of g = d * h; in conclusion: If we now plug this into (5.22) and bound 1 − cos x ≤ x 2 /2, we find Recalling that n q ∞ ≤ c q 4 2 and supp n q ⊆ B(q), we find , n q 2 ≤ c 2 q 8 2 e,e ∈B(q) (e) (e ).
Using this bound, we get: provided that β is large enough. Therefore, where in the last step we used Young's inequality and Finally, we note that , = g, A −1 g that, combined with (5.9), implies the desired estimate, the first of eq.(5.8). A step by step repetition of the argument with h and g replaced byh andg, respectively, leads to the second of (5.8).

Proof of Theorem 3.3
In this section we compute the energy of a dislocation dipole q n dip , see (3.12), asymptotically as n → ∞, and the energy of two parallel arrays of dipoles q M,n,m grain , see (3.13), asymptotically as M → ∞ first, then n → ∞, then m → ∞.
Note that, also in two dimensions, the Ariza-Ortiz energy satisfies the additive decomposition property (4.9), from which we get min{H AO (u, σ) : dσ = q} = 1 2 q, G * BGq , (6.1) see (4.13) (with some abuse of notation, in this section we denote by A, B, G the analogues for the triangular lattice of the operators A, B, G introduced in Sect.4.2 for the case of the FCC lattice).

6.1
The energy of a dipole (proof of (3.16)).
Let σ n dip be such that dσ n dip = q n dip . Then The optimal 0-form u satisfies Au = d * Bσ n dip and consequentially We choose else.
It is a simple exercise to check that dσ n dip = q n dip is satisfied. For a visualization of the support of the slip-field σ n dip see Fig. 2. We compute the energy in (6.2) by using this σ n dip . We start by computing Bσ n dip : which implies that the first term in the right side of (6.2) satisfies (recall that the number of red bonds equals 2n) In order to explicitly compute the second term in the right side of (6.2) in the thermodynamic limit, it is convenient to fix a convention for the Fourier transform: given functions u, σ, q on vertices, edges, faces of the infinite triangular lattice, respectively, we let is the (first) Brillouin zone, and |B| = 8π 2 / √ 3 its area. With this notation, the thermodynamic limit of the second term in the right side of (6.2) can be written as whereÂ(k) is the Fourier symbol of A. Note that, for any 0-form u, where Π l = b l ⊗ b l is the projector in direction b l . Therefore, passing to Fourier space, we getÂ and detÂ(k) = 0 ⇔ k = 0 mod L * . In order to compute d * Bσ n dip (k) in (6.4), we first note that Moreover, for any 1-form f , In conclusion, In the vicinity of the singularity, letting Π k = k ⊗ k/k 2 be the projector in direction k, ). Using these properties ofÂ −1 (k) we see that the function F (k) in (6.7) is even, uniformly bounded on B, and analytic in k away from k = 0. In the vicinity of the singularity, it behaves like In order to extract the dominant contributions from (6.6), we rewrite F (k) = F ((0, k 2 )) + [F (k) − F ((0, k 2 ))], with F ((0, k 2 )) = 2. The contribution from F ((0, k 2 )) reads, for any small > 0: where the remainder O(1) is uniformly bounded as n → ∞. By using (6.8), we can rewrite the contribution from [F (k) − F ((0, k 2 ))], for any small > 0, as where, again, the remainder O(1) is uniformly bounded as n → ∞. The dominant term in (6.9) can be computed explicitly, and gives Putting things together, we obtain (3.16), as desired.
6.2 The energy of a pair of infinite, parallel, grain boundaries (proof of (3.17)). Let with σ n dip defined in (6.3). By proceeding as in the previous subsection, we find that where F (k) is the same as in (6.7). If we now let M → ∞, where p j = p j (m) = 2πj m √ 3 . In order to compute this expression asymptotically, as n, m → ∞, it is convenient to rewrite F ((k 1 , p j )) as [F ((k 1 , p j )) − F ((0, p j ))] + F ((0, 2πj m √ 3 )), where F ((0, p j )) = 2 (in the case j = 0, this identity should be understood as lim k 1 →0 F ((k 1 , 0)) = 2). The contribution from F ((0, p j )) reads: while the one from [F ((k 1 , p j )) − F ((0, p j ))] reads: A computation shows that the difference in brackets is O(k 2 1 ) for k 1 close to 0 (possibly non-uniformly in j, m); correspondingly, if we let n → ∞, the term proportional to cos(nk 1 ) under the integral sign goes to zero as (log n)/n. Summarizing, The dominant contribution to the first term in the right side as m → ∞ comes from the region (k 1 , p j ) ∈ [− , ] 2 , the contribution from the complement being bounded uniformly in m (here is an arbitrary small, positive, constant). Moreover, by rewriting 1 − cos k 1 for k 1 small as k 2 1 2 (1 + O(k 2 1 )), and by expanding F (k) as in (6.8), we find, letting J = √ 3 m/(2π) , Finally, recalling that p j = 2πj m √ 3 , by summing over j and integrating over k 1 , we find: asymptotically as m → ∞. This is the desired 'Read-Shockley' law for the energy of a grain boundary.

Comparison of the Read-Shockley formula with the capacitor law.
As promised above, let us now make a technical comparison between the derivation of the Read-Shockley formula (6.13) and the analogous computation in the case that the operator B is replaced by the identity. In this case we lose the key feature of our discrete elasticity model, that is, invariance under linearized rotations. This is the physical reason why the scaling of the corresponding energies are completely different. More specifically, let u be a R 2 -valued function on T , σ a lattice-valued function on the nearest neighbor bonds of T , and q a lattice-valued function on the faces of T , with finite support and zero total charge, f q(f ) = 0. Consider the minimum energy defined by where the minimum over v is performed over R 2 -valued (rather than T -valued) functions on the nearest neighbor bonds of T . The minimizer σ q is characterized by the Euler-Lagrange equations dσ q = q and d * σ q = 0. Clearly σ q = d * ∆ −1 q satisfies those equations and is the unique minimizer. Hence i.e. the modified Ariza-Ortiz energy reduces to a lattice Coulomb interaction. To compute E for specific 2-forms q it is convenient to work with the Fourier representation of the Laplacian acting on 2-forms. A simple calculation shows that where we use the abbreviation q(x, Fig. 3 for a visualization of the two face types. We write the Fourier transform of 2-forms q by and obtain the Fourier symbol where Ω(k) = 1 + e −ik·b 2 + e ik·b 3 and we used the convention that each coefficient of the symbol is interpreted as a multiple of the 2-dimensional identity matrix (i.e. the symbol is actually a Hermitian 4 × 4-matrix). In conclusion, Let us now compute the energy of the dislocation dipole: If q = q n dip (cf. (3.12)) then the corresponding energy is: Note that, close to the singularity k = 0, we have |Ω(k)| 2 = 9 − 3 2 |k| 2 + O(k 3 ), so that, for any > 0, which is qualitatively the same as the energy of the dislocation dipole, (3.16). Since the energy of a single dipole is asymptotically the same at large distances, up to a multiplicative constant, both for this lattice Coulomb case and the standard case of the Ariza-Ortiz model, one may naively expect that the energy of two parallel arrays of dipoles is also qualitatively the same in the two models. However, this is not the case. If we consider a charge distribution which resembles two parallel capacitor plates where q = q M,n,m grain (cf. (3.13)), then The dominant contribution to the right side as n → ∞ at m fixed comes from the region (k 1 , p j ) ∈ [− /m, /m] 2 mod L * , for any small > 0, the contribution from the complement being bounded from above uniformly in n, as n → ∞ (this is an immediate consequence of the fact that the only zero of 9−|Ω(k)| 2 is in k = 0). On the other hand, the contribution from (k 1 , p j ) ∈ [− /m, /m] 2 mod L * grows linearly in n, as n → ∞, so that, noting that the 9 − |Ω(k 1 , 0) which is the usual electrostatic energy of an infinite capacitor. Note the linear behaviour of the energy in n, as n → ∞, to be compared with the asymptotic independence of the energy of a grain boundary in n in the Ariza-Ortiz model, see (3.17).

A Interpretation of the model
In this appendix we provide a heuristic interpretation of the Ariza-Ortiz model, which aims at clarifying its connection with more fundamental microscopic models for atomic crystals. For a more systematic discussion, from a different perspective, the reader is referred to the original paper [2] where the model has been introduced.
For simplicity, we restrict our discussion to two-dimensions. Assume that the particles interact via a classical, rotationally invariant, pair potential V , with a non-degenerate minimum at distance, say, 1. If such a minimum is deep and narrow, the potential energy of a particle configuration can be well approximated by a sum over nearest neighbors: where the sum runs over ordered pairs of nearest neighbor sites of the graph DT (z), the Delaunay triangulation of the particle configuration z (i.e., the dual of the Voronoi diagram of z); in (A.1), z(ξ) and z(η) indicate the coordinates in R 2 of the vertices of the graph DT (z) labelled ξ and η, respectively. Under the assumption that the pair potential V has a deep, narrow, minimum, located at 1, we expect that the low-energy particle configurations are such that the nearest neighbor pairs involved in the summation in (A.1) have distance |z(ξ) − z(η)| close to 1. See Fig. 4 for an example. Figure 4: A particle configuration z and its Delaunay triangulation DT (z). Notice that, in this example, the graph DT (z) is not a global deformation of the regular triangular lattice: there are sites with 7 neighbours ((0,0) and (3,1)) next to sites with 5 neighbours ((0,1) and (3,0)). These kinds of defects correspond to dislocations, which play a key role in our analysis.
The inconvenient feature of (A.1) is that the sum runs over the nearest neighbor sites of a graph, whose structure depends upon the configuration itself. A more convenient way of expressing the same energy is to reduce to a fixed reference graph, after appropriate redefinition of the nearest neighbor bonds. In our 2D setting, the natural reference graph is the regular triangular lattice of unit mesh, denoted T . Given a particle configuration 4 z, we establish a one-to-one correspondence of the vertices of DT (z) with those of T (call it the 'vertex correspondence') and a one-to-one correspondence of the ordered nearest neighbor pairs of DT (z) with those of T (call it the 'bond correspondence'), in such a way that (A.1) is re-expressed as a sum over ordered pairs of nearest neighbor sites of T (note that the bond correspondence we introduce is not the one induced by the vertex correspondence, see below for its definition and an explicit example): x,y∈T : x∼y V (|z(ϕ(x, y)) − z(ψ(x, y))|).
Here (ϕ(x, y), ψ(x, y)) is the image of an ordered nearest neighbor pair in DT (z) under the aforementioned bond correspondence 5 . Letting where u and σ are the displacement and slip fields, respectively, eq.(A.2) can be further rewritten as The slip-field σ is a key feature of our model, and it provides a direct way of measuring the change of the nearest neighbors, as well as the occurence of dislocations. For example, if the 'charge' q = i σ(x i , x i+1 ) is different from zero for some closed path x i ∈ T , then the path encloses a dislocation defect. The Burger's vectors are closely related to the charges q associated with elementary circuits (the boundaries of the triangular faces of T ), but also depend on the orientation and the location of the path; their specific definition is unimportant for our purposes and, therefore, we skip it.
With these conventions, the slip field on positively oriented edges is The Ariza-Ortiz model is obtained from eq.(A.5) under a couple additional approximations. Letting x − y ≡ 0 and u(ϕ(x, y)) − u(ψ(x, y)) + σ(x, y) ≡ δ 0 , we rewrite V (| 0 + δ 0 |) by expanding it in Taylor series around the minimum: recalling that | 0 | = 1 and assuming δ 0 to be small, we find: Neglecting the remainder, dropping an additive constant and rescaling the resulting energy, we obtain: This quadratic approximation corresponds to the standard small-strain assumption in continuum mechanics. Finally, we replace u(ϕ(x, y)) − u(ψ(x, y)) ≈ u(x) − u(y), (A.8) 6 We use the convention that the 'positively oriented' nearest neighbor pairs of T are those of the form ; of course, the nearest neighbor pairs of the form (x, x − b l ) are called negatively oriented. It is sufficient to define the functions ϕ, ψ and the slip field on positively oriented n.n. pairs: if (x, y) is negatively oriented, we let (ϕ(x, y), ψ(x, y)) = (ψ(y, x), ϕ(y, x)), so that σ(x, y) = −σ(y, x). which corresponds to the 'linearized plasticity' approximation in continuum mechanics. After these replacements, we obtain E(z) ≈ 1 4 x,y∈T : x∼y which is the Ariza-Ortiz model. All the approximations involved in the previous scheme are uncontrolled, and their validity should be (at least) checked a posteriori, by showing that the 'typical', low energy, configurations of the Ariza-Ortiz Hamiltonian are close (in a sense to be defined) to those of the original, realistic, Hamiltonian. This remains to be done: in fact, proving (even at heuristic level) the correctness of these approximations is a major challenge in the field and goes beyond the purposes of our paper.

B On the energy of a grain supported on an infinite vertical strip
In this appendix we discuss the connection between the energy E grain (n, m) defined in (3.15) and the optimal grain energy E G (S) defined in (3.8), in the case of a grain G supported in a vertical strip of width n, 'rotated' by an 'angle' θ ∼ 1/m. As anticipated in Sect.3.3, for simplicity, we restrict our attention to two-dimensions. We recall that T is the infinite 2D triangular lattice with basis vectors b 1 = 1 0 , b 2 = −1/2 √ 3/2 , and we let b 3 = −b 1 − b 2 . We also let m 1 = 4π √ 3 1 be a basis of the dual lattice T * , such that b i · m j = 2πδ i,j , for i, j = 1, 2; moreover, we define We consider a grain whose support is an infinite vertical strip of width n: Figure 7: A portion of the grain: the dark gray area represents the grain, and the circled sites are those belonging to its left boundary. The light gray area is the left 'boundary layer', and the red bonds (the 'boundary bonds') are those connecting sites in the grain with sites in its complement: those are bonds on which the slip field σ in the minimization problem Eq.(B.4) may be non zero; correspondingly, the light gray faces are those where the 'charge' dσ may be non zero.
Note that the complement of the grain is disconnected and consists of two semi-infinite components, the left one, denoted G c L , and the right one, denoted G c R . We let for some θ ∈ R (the normalization factor √ 3 is chosen for later convenience, e.g., for having a prefactor θ 2π , rather than θ 2π √ 3 , in eq.(B.5)) and τ, τ R to be fixed. We are interested in estimating the optimal grain energy per unit vertical length that, in analogy with the definition (3.8) for the finite-grain case, is defined as follows: where we recall that Λ = Λ (N ) is the 2D analogue of (2.2), whose vertical height is √ 3 2 N (which explains the normalization factor 1/( √ 3 2 N ) in the right side of (B.3)), and M ( ) S (G) is the set of lattice-valued slip fields σ such that there exist τ, τ R for which σ, τ, τ R realize the infimum of inf τ,τ R inf * σ H AO (u S , σ) within a precision , where the * on inf * σ indicates the constraint that the support of σ is over the bonds connecting the grain with its complement, such as the red bonds of Fig.7.
Let us focus on the minimization problem defining the set M ( ) Note that by definition the only non-zero terms in the sum in the right side are those associated with boundary bonds (x, y), with x ∈ G and y ∈ G (such as the red bonds of Fig.7): for such bonds u S (x) − u S (y) = Sx + τ , if x belong to the left boundary, and u S (x) − u S (y) = Sx + τ − τ R , if x belong to the right boundary. The goal is to find a minimizer σ(x, y) for (x, y) a boundary bond. For this purpose, it is convenient to decompose S into simple slip systems. Recall that any 2 × 2 skew-symmetric matrix A can be decomposed into simple slip systems, i.e., A = 3 l=1 ξ l b l ⊗ m n(l) , for suitable coefficients ξ l , where m n l ∈ T * are the slip plane normals, namely: m n(1) = m 2 , m n(2) = m 1 , and m n(3) = m 3 . In particular, a simple computation shows that Setting temporarily τ = 0, Eq.(B.5) suggests the following choice for the slip-field minimizer for (x, y) a boundary bond with x ∈ G and y ∈ G: where x := x + 1 2 denotes the 'nearest integer function'. On the left boundary (similar considerations hold for the right one), x equals n 2 (b 2 − b 3 ) or n 2 (b 2 − b 3 ) + b 2 for some n 2 ∈ Z, depending on whether x is in G a or G b . Note that, for x = n 2 (b 2 − b 3 ), we have: 1 2π x · m 2 = 2n 2 , 1 2π x · m 1 = n 2 , and 1 2π x · m 3 = −n 2 . Therefore, on the left boundary, if up to O(θ) corrections, which are present if x has the form n 2 (b 2 − b 3 ) + b 2 (and are small, for θ small). The computation leading to (B.7) neglected the presence of τ in the definition of u S and was not based on an exact minimization of the sum in the right side of (B.4). However, the patient reader can check that, by performing an exact minimization, one can choose τ of order θ and σ can be chosen as in (B.7), up to a bounded O(θ) correction (details left to the patient reader). Similarly, the exact minimization along the right boundary leads to a slip field equal to the opposite of (B.7), up to a bounded O(θ) correction.
In conclusion, neglecting O(θ) fluctuation terms in the boundary slip field, which are not expected to contribute to the optimal grain energy per unit vertical length at the dominant order in the limit of small θ and large n, the optimal slip field for the minimization problem (B.4) equals ∓b 1 ( 2θn 2 + θn 2 ) on the boundary bonds of the grain with vertical coordinate n 2 (the minus and plus signs are for the left end right boundaries, respectively). Notice that the charge distribution dσ of such a slip field consists of isolated charges equal to ±b 1 or ±2b 1 on suitable faces of the left and right boundaries (with opposite signs on the two boundaries), vertically separated in average by a distance 1 2θ . The average density of such boundary charge distribution equals ±b 1 /(3θ), the same as the charge distribution (3.13), provided we identify m with 1/(3θ).

C The cellular complex of the FCC lattice
We recall that the three-dimensional FCC lattice L is the Bravais lattice with basis vectors b 1 , b 2 , b 3 , as in (2.1). We also define the dual lattice L * as the Bravais lattice with basis vectors m 1 , m 2 , m 3 , defined as Note that b i · m j = 2πδ i,j , with i, j = 1, 2, 3. For later reference, we also let m 4 = In terms of these definitions, the cellular complex associated with the FCC lattice is defined in terms of the following cells: 1. The vertices x ∈ E 0 are the vertices of L, of the form x = n 1 b 1 + n 2 b 2 + n 3 b 3 .
2. The edges e ∈ E 1 are the ordered pairs of nearest neighbour vertices of L, namely pairs (x, x ) with x −x = ±b l , l = 1, . . . , 6: here b 1 , b 2 , b 3 are the same as (2.1), and we recall that The action of the boundary operator on E 1 is defined by: ∂(x 1 , x 2 ) = {x 1 , x 2 }, for any (x 1 , x 2 ) ∈ E 1 . Note that, in the notation of Sect.4.1, ∂e = V (e), ∀e ∈ E 1 , where V (e) is the set of vertices of e. The primitive unit cell of the FCC lattice is shown in the right panel. It can be dissected into a regular octahedron and two regular tetrahedra, which are the 3-cells of our cellular complex. The figure shows that there are two inequivalent type of tetrahedra: red and green. We shall call r-tetrahedra (resp. g-tetrahedra) those that can be translated into the red (resp. green) tetrahedron.

D On the operator A and its inverse
In this appendix we discuss and prove a few basic properties of the operator d * 0 Bd 0 , both in the case that it acts on the 0-forms associated with the infinite FCC lattice L, and in the case that it acts on those associated with a finite box Λ= Λ (N ) ⊂ L, of the form (2.2), with Dirichlet boundary conditions. In order to avoid confusion between the two cases, in this appendix (contrary to the rest of the paper) we denote by C 0 , resp. C 0 Λ , the set of 0-forms associated with the infinite lattice, resp. with the box Λ with Dirichlet boundary conditions. Correspondingly, we denote by A, resp. A Λ , the operator d * 0 Bd 0 acting on C 0 , resp. C 0 Λ . Note that A Λ can be rewritten as with k l = 1 2π k · b l and e j the j-th standard Euclidean basis vector. We have:Ã Λ u k,j (x) = 2 3 l=1 α l Π l ij u k,i (x), where α l := 1 − cos(2πk l ), which is positive for k ∈ Λ * D . Of course, 2 3 l=1 α l Π l ≥ 2 min{α 1 , α 2 , α 3 } 3 l=1 Π l . By using the explicit form of Π l , we get 2 3 l=1 Π l =   2 1 1 1 2 1 1 1 2   , whose smallest eigenvalue is 1, that is, 2 3 l=1 Π l ≥ 1.
D.2 Proof of (5.10) In order to prove (5.10), we derive upper and lower bounds on g, A −1 Λ g , that is, the argument of the limit in the left side of (5.10), in the notation of this appendix. For the reader's convenience, we recall that g = g v 0 ;x = 1 x v 0 , where x ∈L and v 0 ∈ L * . With no loss of generality (since we are interested in the thermodynamic limit Λ→L), we assume that x ∈ Λ. The important feature of g to be used in the following is that it is compactly supported, with support contained in Λ. Note that − g, A −1 Λ g = min We recall that the minimum in the right side is over the compactly supported 0-forms u : L → R 3 , whose support is contained in Λ. In order to get a lower bound, we write the quadratic function u, Au − 2 u, g in Fourier space, by using the convention u(z) = withĝ(k) andÂ(k) defined as in (5.11)-(5.12). As anticipated in Sect.5.1,Â −1 (k) is singular only at k = 0, close to which it behaves like ∼ k −2 , see below for a proof: therefore, the right side of (D.3) is finite for any compactly supported g.
In order to get an upper bound, we use the test function u * (z) := χ Λ (z)u ∞ (z), where χ Λ (z) := min{1, 4 dist(z, Λ c )/N } (recall that N is the side of the box Λ, see (2.2)) and du * (e) · δe 2 − 2 u ∞ , g , (D.5) where we recall that, for an edge e = (x, y), δe = y − x, and in the last identity we used the fact that, thanks to the definition of χ Λ , for N sufficiently large the support of g = 1 x v 0 is contained in supp1(χ Λ = 1) ⊂Λ, so that in particular u * = u ∞ on the support of g. Moreover, letting u(z e ) := 1 2 (u(x) + u(y)) for an edge e = (x, y), Now, for Λ = Λ (N ) large enough (see (2.2)), the right side of (D.9) vanishes as N → ∞, which concludes the proof of (5.10).
This completes the proof thatÂ(k) is invertible iff k = 0 mod L * .
Then, we rewrite it as: Note that ∂ k jÂ −1 (k) = −Â −1 (k)·∂ k jÂ (k)·Â −1 (k), with ∂ k jÂ (k) = 2 6 l=1 Π l (b l ) j sin(k· b l ). Recalling thatÂ(k) is even, it is singular iff k = 0, and, for k close to zero, it can be bounded from above and below by (const.)k 2 , we find that ∂ k jÂ −1 (k) is odd, it is singular iff k = 0 and, close to the singularity, it can be bounded from above by (const.)|k| −3 . Therefore, the right side of (D.15) can be rewritten as and, in order to bound it from above, we multiply the integrand by 1 = χ(k) + (1 − χ(k)), where χ(k) is a positive, monotone, C ∞ radial function, equal to 1 for |k| ≤ and equal to 0 for |k| ≥ 2 . Now, the term associated with (1 − χ(k)) is the Fourier transform of a C ∞ function and, therefore, it decays faster than any power in real space. The term associated with χ(k) can be bounded as follows: Putting things together, we obtain that, if |x| ≥ −1 , then (Â −1 (k)v 0 ) l e −ik·(x−y) ≤ (const.) log |x|.
Summing over j from 1 to 3, we get the desired estimate, (D.14).
E.2 Proof of (5.26) Plugging (5.21) in the left side of (5.26), and using the fact that | supp (q)| ≤ q 1 and q 2 ≤ q 1 , we find We now weaken the constraint that B(q) e, e into q 1 ≥ dist(e, e ) and find that, for β large enough, as desired.