Xavier Blanc

Claude Le Bris

Pierre Louis

P.-L Lions

From Molecular Models to Continuum Mechanics

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

It is commonly admitted in the continuum mechanics literature that the stored energy of an elastic crystal is of the form:

E(u) = E(∇u(x))dx, (1)
where ⊂ R 3 is the reference configuration of the solid, and u is the deformation to which the solid is subjected. Moreover, the stored-energy density E is assumed to reflect the microscopic symmetries of the crystal [START_REF] Ashcroft | Solid State Physics[END_REF][START_REF] Ball | Convexity conditions and existence theorems in nonlinear elasticity[END_REF][START_REF] Chipot | Equilibrium configurations of crystals[END_REF][START_REF] Fonseca | Variational methods for elastic crystals[END_REF][START_REF] Zanzotto | On the material symmetry group of elastic crystals and the Born rule[END_REF], in addition to the standard frame invariance [START_REF] Ciarlet | Mathematical theory of elasticity[END_REF][START_REF] Love | A treatise on the mathematical theory of elasticity[END_REF]. In other words, denoting by G the invariance group of the underlying crystalline structure, we have:

∀M ∈ M 3×3 + , ∀Q ∈ G, ∀R ∈ SO(3), E(RMQ) = E(M), (2)
where M 3×3 + denotes the set of three-by-three matrices having positive determinant, on which the energy density E is supposed to be defined. This invariance property is at the origin of many important properties of crystalline solids, as well as mathematical difficulties in the use of energies of the form [START_REF] Abramowitz | Handbook of mathematical functions[END_REF]. See for instance [START_REF] Chipot | Equilibrium configurations of crystals[END_REF][START_REF] Ericksen | Weak martensitic transformations in Bravais lattices[END_REF][START_REF] Fonseca | Variational methods for elastic crystals[END_REF][START_REF] Üller | Variational models for microstructure and phase transitions[END_REF][START_REF] Zanzotto | On the material symmetry group of elastic crystals and the Born rule[END_REF]. However, an exact expression, or even an approximation for E is rarely available [START_REF] Weiner | Statistical mechanics of elasticity[END_REF]. Closely linked with this problem is the question of the associated functional space: To what space A should u belong in order to properly define the minimization problem I A = inf E(u), u ∈ A ? The definition of A should include, in addition to regularity properties, suitable boundary conditions.

In order to answer these questions, a rigorous derivation of ( 1) and ( 2) is needed, starting from the atomic level, since it is the scale at which the crystalline lattice is present. A standard approach to this kind of problem is -convergence [START_REF] Maso | An introduction to -convergence[END_REF], which assumes implicitly that at an intermediate level, the microscopic (or mesoscopic) deformation, which is not necessarily equal to the macroscopic deformation u, should locally minimize the energy. Roughly speaking, this approach consists in setting a microscopic version of the minimization problem I A in terms of the microscopic deformation u ε , depending on the interatomic distance ε. If we let this distance go to zero, this yields a (weak) limit deformation u = u 0 together with an energy, possibly of the form [START_REF] Abramowitz | Handbook of mathematical functions[END_REF]. However, it is far from clear what should be the microscopic counterpart of I A : what regularity should the deformation exhibit at this level, and how should the boundary conditions (which are by nature macroscopic) be translated microscopically? This kind of strategy is used in [START_REF] Braides | Limit of discrete systems with long-range interactions[END_REF][START_REF] Pagano | A simple model for phase transitions: from the discrete to the continuum problem[END_REF][START_REF] Paroni | Homogenization of polycrystalline aggregates[END_REF] in the case of a two-body finite-range energy. In [START_REF] Friesecke | A scheme for the passage from atomic to continuum theory for thin films, nanotubes and nanorods[END_REF], the case of a more complex microscopic model is considered.

Another (and seemingly more naive) approach is the following: the whole point is in fact to link the macroscopic deformation u appearing in [START_REF] Abramowitz | Handbook of mathematical functions[END_REF] with the deformation truely experienced by the atoms of the solid, that is, a sort of microscopic deformation. Since this link is far from being clear physically, let us assume the simplest link, that is, equality. This is the approach used in [START_REF] Bardenhagen | On higher order gradient continuum theories in 1-D nonlinear elasticity. Derivation from and comparison to the corresponding discrete models[END_REF][START_REF] Bardenhagen | Derivation of higher order gradient continuum theories in 2,3-D non-linear elasticity for periodic lattice models[END_REF], and also the one we adopt here. Let us emphasize that we are aware of the lack of physical justification of this assumption [START_REF] Friesecke | Validity and failure of the Cauchy-Born hypothesis in a 2D mass-spring lattice[END_REF]. With a view to tackling more realistic cases in the future, the present work seems to us a necessary and useful preliminary mathematical step, before we can relax the present simplifying assumptions. Let us give the example of a two-body interaction, in order to fix the ideas: the energy E({X i }) of N identical atoms of positions X i is then given by

E({X i }) = 1 i<j N W (X i -X j ), ( 3
)
where W is the interaction potential. Let us now assume that a solid is defined by a domain , and a lattice , and that the interatomic distance is equal to ε > 0, which in the end will tend to zero. Assuming that the solid experiences a deformation u, which is a C ∞ diffeomorphism defined on , then the positions of the atoms of the deformed solid are X i = u(X 0 i ), {X 0 i } 1 i N = ε ∩ , the set of points {X 0 i } being the reference state. Consequently, the energy per atom of the deformed configuration is: Next, considering the potential W , its characteristic length δ should be comparable with ε. In other words, considering the reference configuration, which minimizes the energy (3), its atomic spacing should be proportional to δ. (This is clear since if {X 0 i } is a ground state for W , {δX 0 i } is a ground state for W ( • δ ).) We will however also study other cases apart from the ε = δ case, that is, cases when ε δ and δ ε, respectively. The total energy of the deformed system is thus equal to

E(u(ε ∩ )) = 1 2N i | =j ∈ ∩ 1 ε W (u(εi) -u(εj )). (4)
E ε,δ (u) = 1 2N i | =j ∈ ∩ 1 ε W 0 u(εi) -u(εj ) δ , ( 5
)
where N = #( ∩ 1 ε ) is the total number of atoms, and W 0 is the rescaled potential, and does not depend on ε nor on δ. The energy E appearing in (1) should then be the limit, as ε and δ go to zero, of E ε,δ :

E(u) = lim ε,δ→0 E ε,δ (u). ( 6
)
Despite the seemingly crude assumption we make in this strategy (namely, that the microscopic deformation is equal to the macroscopic one), there are some advantages to be found in this approach.

First of all, the above limiting process [START_REF] Benguria | The Thomas-Fermi-von Weizsäcker theory of atoms and molecules[END_REF] is clearly not linked with the fact that the energy originates from a two-body interaction. Indeed, we will see in Section 3 below that more complex microscopic models, taking (at least partially) the quantum nature of the electrons into account, give rise to the same kind of convergence results. This allows us to hope that similar results will hold for more complex microscopic theories.

Another point is that since it is a mere limiting process, involving noconvergence properties, we may very well go further: considering the limit as a zero-order term in the development of the energy with respect to ε, it is possible to compute the next orders. This is what is done in Theorem 3.

In addition, let us point out that for numerical purposes, this approach gives a way of computing the sum (4) in the limit ε → 0, which is the exact expression of the energy, u being the microscopic deformation. Of course, its link with the macroscopic one is still unknown, but independently of this link, the above limit (6) may be seen as a good approximation in the case where ε is small, hence N large, circumventing any direct computation of (4), which indeed may be out of reach of computer facilites for very large N [START_REF] Cleri | Atomistic simulations of materials fracture and the link between atomic and continuum length scales[END_REF][START_REF] Ortiz | Quasicontinuum analysis of defects in solids[END_REF].

Finally, let us point out that our approach is intimately linked with the concept of a thermodynamic limit [START_REF] Lieb | The Thomas-Fermi theory of atoms, molecules and solids[END_REF][START_REF] Catto | Mathematical Theory of thermodynamic limits : Thomas-Fermi type models[END_REF], as explained in Section 3, which is a -limit process, although normally not presented as such.

The article is organized as follows: Section 2 gives the computation of the limit (6), together with the corresponding development, up to order two with respect to ε (Theorem 3). Section 3 then studies the limit [START_REF] Benguria | The Thomas-Fermi-von Weizsäcker theory of atoms and molecules[END_REF] in the case of some quantum models, namely Thomas-Fermi type theories. Let us emphasize that the use of these simplified models is only for mathematical purposes: the same problem with more sophisticated models should be addressed, but it seems for now beyond our reach, technically. We give a few possible extensions of the present work in Section 4, while the last section is devoted to the derivation of a few simple properties of the computed elastic energies. We hope to come back to these questions in the near future.

Most of the results detailed here have been announced in [START_REF] Blanc | Convergence de modèles moléculaires vers des modèles de la mécanique des milieux continus[END_REF].

The simplest case: two-body potentials

We present here the homogenization scheme described in Section 1 in the simplest case, that is, when the microscopic energy is defined by a two-body potential.

In other words, we study the limit [START_REF] Benguria | The Thomas-Fermi-von Weizsäcker theory of atoms and molecules[END_REF], with the energy E ε,δ being given by [START_REF] Bardenhagen | Derivation of higher order gradient continuum theories in 2,3-D non-linear elasticity for periodic lattice models[END_REF], that is:

E ε,δ (u) = 1 2N i | =j ∈ ∩ 1 ε W 0 u(εi) -u(εj ) δ .
Assuming that the reference configuration is an equilibrium state, i.e., that an infinite system with minimal energy is periodic, it is physically reasonable to assume that δ and ε have comparable size, i.e., that in the limit process ε → 0, we should take δ = ε. We will nevertheless study the other possible cases, that is, when ε δ or δ ε. Note also that the fact that the equilibrium configuration is periodic is not proven, so far as we know, except in one dimension [START_REF] Blanc | Periodicity of the infinite-volume ground state of a onedimensional quantum model[END_REF][START_REF] Gardner | The infinite-volume ground state of the Lennard-Jones potential[END_REF][START_REF] Ventevogel | On the configuration of a one-dimensional system of interacting particles with minimum potential energy per particle[END_REF][START_REF] Nijboer | On the configuration of systems of interacting particles with minimum potential energy per particle[END_REF][START_REF] Nijboer | On the configuration of systems of interacting particles with minimum potential energy per particle[END_REF][START_REF] Radin | Classical ground states in one dimension[END_REF][START_REF] Radin | Periodicity of classical ground states[END_REF], or in very simple cases in two dimensions [START_REF] Radin | Ground states for soft disks[END_REF].

Zero-order term

Hereafter, we denote by ε a sequence of positive real numbers converging to zero, and δ = δ(ε) also going to zero. The property ε δ means that lim ε→0 ε δ(ε) = 0, and δ ε means that lim ε→0 δ(ε) ε = 0, these limits being taken along the corresponding sequence. [START_REF] Bardenhagen | Derivation of higher order gradient continuum theories in 2,3-D non-linear elasticity for periodic lattice models[END_REF]. Then the following statements hold:

(i) If ε = δ, then lim ε→0 E ε,δ (u) = 1 2| | j ∈ \{0} W 0 (∇u(x)j )dx. ( 7
) (ii) If ε δ, and if W 0 ∈ L 1 (R 3 ), then lim ε→0 ε δ d E ε,δ (u) = 1 2| | R 3 W 0 dx | det(∇u(x))| . ( 8
) (iii) If δ ε, and if for some p ∈ R, lim |x|→∞ |x| p W 0 (x) = a, then lim ε→0 ε δ p E ε,δ (u) = 1 2| | j ∈ \{0} a |∇u(x)j | p dx. ( 9
)
Proof. We first prove (i). Note that this result seems to be well known, but since we found no rigorous proof of it in the literature, we provide the reader with one here. Consider [START_REF] Bardenhagen | Derivation of higher order gradient continuum theories in 2,3-D non-linear elasticity for periodic lattice models[END_REF], with ε = δ, that is:

E ε,δ (u) = 1 2N i | =j ∈ ∩ 1 ε W 0 u(εi) -u(εj ) ε .
We first split this sum into two sums, using a cut-off radius A > 0, which in the end will go to infinity:

E ε,δ (u) = 1 2N |i-j | A W 0 u(εi) -u(εj ) ε + 1 2N |i-j |>A W 0 u(εi) -u(εj ) ε , ( 10
)
both sums being reduced to couples

(i, j ) in 1 ε ∩ 2 such that i | = j . The deformation u being a C ∞ diffeomorphism, we have α|i -j | |u(εi)-u(εj )| ε
β|ij | for some constants α and β depending only on u. Consequently, the second sum in [START_REF] Bohr | Almost periodic functions[END_REF] is easily seen to vanish as A goes to infinity, uniformly with respect to ε. Turning to the first term, we use a Taylor expansion to write:

u(εi) -u(εj ) ε -∇u(εi)(j -i) CA 2 ε,
where C depends only on u. Next, using the fact that W 0 is Lipschitz on the exterior of any ball, and that αd( ) α|ij | |u(εi)-u(εj )| ε for any i and j appearing in the sum, where d( ) is the minimal distance between two points of the lattice , we may write:

W 0 u(εi) -u(εj ) ε -W 0 ∇u(εi)(j -i) C d(α) A 2 ε.
Using this inequality in [START_REF] Bohr | Almost periodic functions[END_REF], we thus have:

E ε,δ (u) - 1 2N i∈ ∩ 1 ε 0<|i-j | A W 0 ∇u(εi)(j -i) C d(α) A 2+d ε + R(A), ( 11
)
where R(A) corresponds to an upper bound of the second term of [START_REF] Bohr | Almost periodic functions[END_REF], and may thus be chosen equal to

R(A) = 1 2N |i-j |>A C |i-j | d+a C
A a , where C are various constants independent of ε. Next, the sum in (11) may be written as

1 2N i∈ ∩ 1 ε 0<|j | A W 0 ∇u(εi)j ,
where the sum over j is restricted to j ∈ (as in [START_REF] Braides | Limit of discrete systems with long-range interactions[END_REF]). Then, using the regularity of u, we may write |∇u(εi)j | α|j |, where α does not depend on ε and i. Consequently, the sum over j is, up to a term going to zero as A goes to infinity, uniformly with respect to ε, the sum over the whole lattice:

E ε,δ (u) - 1 2N i∈ ∩ 1 ε j ∈ \{0} W 0 ∇u(εi)j C d(α) A 2+d ε + C A a .
The last step is to point out that, since N = #ε ∩ , the sum over i is a Riemann sum, converging to the desired integral since the function

x → j | =0 W 0 (∇u(x)j ) is Lipschitz on . Choosing A = ε -d 2 -1
, and then letting ε go to zero, we prove (i). Let us now prove (ii). For the sake of simplicity, we restrict ourselves to the case of a finite-range potential, the generalization to the present case being only a technical matter, similar to the cut-off trick in the proof of (i). We then have, for i

, j ∈ ∩ 1 ε , |u(εi) -u(εj )| δ α ε δ |i -j |
for some constant α depending only on u. Consequently, given i ∈ ∩ 1 ε , the number of terms in the sum over j contributing to the energy is of order δ d ε d , in view of the inequality |ij | C δ ε . Therefore, we easily show, using the fact that W 0 and u are Lipschitz, that

ε 2d i | =j ∈ ∩ 1 ε W 0 u(εi) -u(εj ) δ - W 0 u(x) -u(y) δ dxdy C ε δ + ε δ d ,
where C depends only on W 0 , u and . Now, changing variables in the integral, we have

W 0 u(x) -u(y) δ dxdy = u( ) u( ) W 0 ξ -η δ | det(∇u(u -1 (ξ ))) det(∇u(u -1 (η)))| dξ dη.
Next, we use the fact that as δ goes to infinity, 1

δ d R d W 0 W 0 ( • δ ) converges to a
Dirac mass at zero, so that

lim δ→0 1 δ d u( ) u( ) W 0 ξ -η δ | det(∇u(u -1 (ξ ))) det(∇u(u -1 (η)))| dξ dη = R d W 0 u( ) dξ | det(∇u(u -1 (ξ )))| 2 = R d W 0 dx | det(∇u(x))| .
Putting all this together, we thus have:

lim ε δ→0 Nε 2d δ d E ε,δ (u) = 1 2 R d W 0 dx | det(∇u(x))| .
Observing that, at leading order, N ≈ | | |Q( )|ε d , where |Q( )| is the volume of the unit cell of , this proves (ii).

We next prove (iii): fixing a positive parameter α (which in the end will go to zero), we know that for |x| large enough,

|W 0 (x)-a |x| p | α |x| p .
Using the regularity of y, we thus have, for ε δ large enough,

W 0 u(εi) -u(εj ) δ - aδ p |u(εi) -u(εj )| p Cαδ p |i -j | p ,
where C depends only on u. Hence, we deduce easily that

ε δ p E ε,δ (u) - 1 2N i | =j aε p |u(εi) -u(εj )| p Cα 1 2N i | =j 1 |i -j | p Cα j ∈ \{0} 1 |j | p . ( 12
)
We next apply (i) with W 0 (x) = 1 |x| p to conclude that the second term of the righthand side of (12) converges to the desired formula. It only remains to point out that [START_REF] Catto | Mathematical Theory of thermodynamic limits : Thomas-Fermi type models[END_REF] is valid for any α > 0, as long as ε δ is chosen large enough, which concludes the proof.

Remark 2. The regularity assumption we have made on u could be slightly relaxed. Indeed, u being a C 1 diffeomorphism is clearly sufficient to make the above proof available. However, and although this assumption is reasonable (see [START_REF] Ciarlet | Mathematical theory of elasticity[END_REF]), the Taylor expansion we use here could not be carried out if only lower regularity was assumed. If for example we allow jumps in the gradient, it might be possible to obtain a concentration of energy on the jump set (see [START_REF] Braides | Limit of discrete systems with long-range interactions[END_REF]).

Higher order terms

We have studied in Subsection 2.1 the limit of the energy as ε and δ go to zero, that is, the zero-order term of a development of this energy in powers of ε. We give now a derivation of higher order terms of this development, limiting ourselves to order two, although the computations could be carried out at any order.

Before stating this result in Theorem 3 below, we need a definition. Let be a piecewise C 1 open bounded set, and let x ∈ ∂ . Denote by (x) the tangent plane at x, and define the non-negative measure µ (x), on R + by

µ (x), = lim ε→0 ε d-1 k∈A ε ( (x)) N ε (k)δ k , ( 13
) with A ε ( (x)) = {d(i, 1 ε (x)), i ∈ ∩ 1 ε ( (x) -∩ )}, the number is N ε (k) = #{i ∈ ∩ 1 ε ( (x) -∩ ), d(i, 1 ε (x))
= k}, and the set (x) -= {z + ν(x)t, t ∈ (-∞, 0), z ∈ (x)}, where ν(x) is the outer normal of at x. Theorem 3. Let W 0 be a function defined on R d \ {0}, such that for all x | = 0, W 0 (x) = W 0 (-x), which is smooth on the exterior of the ball B R , for any R > 0, and such that there exists C 0, R 0 > 0 and a > 2 satisfying

∀k ∈ N, ∀x ∈ (B R 0 ) c , |D k W 0 (x)| C |x| d+k+a . ( 14
)
Let be a piece-wise C 1 open bounded subset of R d , and let u be a C ∞ diffeomorphism defined on , with values in R d . Assume in addition that there exists a sequence ε n > 0 converging to zero such that, for all integers n, #(

ε n ∩ ) = | | ε d n
, and that |Q( )| the volume of the unit cell of the lattice is equal to one. Consider E δ,ε the energy defined by [START_REF] Bardenhagen | Derivation of higher order gradient continuum theories in 2,3-D non-linear elasticity for periodic lattice models[END_REF]. Then, restricting ε to the sequence ε n the following statements hold:

(i) If δ = ε, then E ε,δ (u) = 1 2| | j ∈ \{0} W 0 (∇u(x)(j )) dx - ε 2| | ∂ ∞ 0 j ∈ , j ν k W 0 (∇u(x)j )dµ (x), (k) dσ (x) - ε 2 24| | j ∈ \{0} D 2 W 0 (∇u(x)j )(D 2 u(x)(j, j ),D 2 u(x)(j, j )) dx + ε 2 F 1 (u) + o(ε 2 ), ( 15
)
where (x) is the tangent plane of ∂ at x, the non-negative measure dµ , (k) is defined by [START_REF] Catto | Limite thermodynamique pour des modèles de type Thomas-Fermi[END_REF], and F 1 contains only boundary terms;

(ii) If ε δ, and if W 0 ∈ L 1 (R 3 ), then ε δ d E ε,δ (u) = R d W 0 2| | dx | det(∇u(x))| - ε | | ∂ ∞ 0 ∞ k W 0 (∇u(x)z)dµ (x), (z) dydσ (x) - ε 2 24| | R d D 2 W 0 (∇u(x)y)(D 2 u(x)(y, y), D 2 u(x)(y, y)) dydx + ε 2 F 2 (u) + o(ε 2 ), ( 16
)
the measure dµ (x), being defined by [START_REF] Catto | Limite thermodynamique pour des modèles de type Thomas-Fermi[END_REF], and the term F 2 (u) containing only boundary terms.

Let us point out that, although formulas [START_REF] Chipot | Equilibrium configurations of crystals[END_REF] and ( 16) seem rather complicated, in the special case where is the unit cube of R d , and the lattice is equal to Z d , the plane (x) is exactly the face of the cube to which x belongs, and the measure dµ (x),Z d is easily computed to be p 1 δ p , the sum involving only p ∈ N. Consequently, in this case, the integrand of the second term of (15) reads:

k 1 j ∈Z d , jn k W 0 (∇u(x)j ),
which may be interpreted, when ∇u is the identity matrix, as the surface energy of a crystal. In a more general setting, the measure dµ (x), may be seen as the average number of points in and in the half-space containing which are at distance k from (x).

Proof. We only give the proof of (i), the proof of (ii) following the same line of arguments. For the sake of simplicity, we assume that the constant a of [START_REF] Catto | On the thermodynamic limit for Hartree-Fock type models[END_REF] satisfies a > 8. This is only a technical assumption, which allows us to simplify the argument below. We will next indicate how to deal with the general case (a > 2).

We use a cut-off radius A > 0, as in the proof of Theorem 1, writing

E ε,δ (u) = 1 2N i∈ ∩ 1 ε 0<|i-j |<A W 0 u(εi) -u(εj ) ε + i∈ ∩ 1 ε |i-j | A W 0 u(εi) -u(εj ) ε , (17)
where all the sums over j concern only j ∈ . Using [START_REF] Catto | On the thermodynamic limit for Hartree-Fock type models[END_REF], the second sum is easily bounded by a term of the form C A a , so that, taking A = ε -β , with 2 a < β < 1 4 , which is possible since a > 8, this sum is of strictly lower order than ε 2 . Next, we use a Taylor expansion of the expression W 0 u(εi)-u(εj ) ε

:

W 0 u(εi) -u(εj ) ε = W 0 ∇u(εi)(j -i) + ε 2 ∇W 0 ∇u(εi)(j -i) × D 2 u(εi)(j -i, j -i) + ε 2 1 8 D 2 W 0 ∇u(εi)(j -i) × D 2 u(εi) j -i, j -i , D 2 u(εi) j -i, j -i + 1 6 ∇u(εi)(j -i) D 3 u(εi) j -i, j -i, j -i + ε 3 72 D 2 W 0 ∇u(εi)(j -i) × D 2 u(εi)(j -i, j -i), D 3 u(εi)(j -i, j -i, j -i) + O ε 3 A 4 |j -i| d+a+1 , ( 18
)
where the term O

ε 3 A 4
|j -i| d+a+1 involves constants depending only on W 0 , u and . Since the term depending on j may be summed up over uniformly with respect to i, this will lead, when summed up with respect to i and j , to a quantity of order O(ε 3 A 4 ) = o(ε 2 ). We next study the term of order three: we need to show that, summed with respect to i and j , it remains of lower order than ε 2 . For this purpose, we estimate it as follows, setting j = ji:

0<|j |<A D 2 W 0 ∇u(εi)(j ) D 2 u(εi)(j , j ), D 3 u(εi)(j , j , j ) 0<|j |<A C |j | d+a+2 |j | 5 = 0<|j |<A C |j | d+a-3 C A a-3 . (19)
This shows that the third-order term in ( 18) is of order ε 3 A a-3 , which is negligible compared to ε 2 , thanks to the definition of A. We may thus write the energy as:

E ε,δ (u) = 1 2N i∈ ∩ 1 ε 0<|i-j |<A W 0 ∇u(εi)(j -i) + ε 2 ∇W 0 ∇u(εi)(j -i) D 2 u(εi) j -i, j -i + ε 2 6 ∇W 0 ∇u(εi)(j -i) D 3 u(εi) j -i, j -i, j -i + ε 2 8 D 2 W 0 ∇u(εi)(j -i) × D 2 u(εi) j -i, j -i , D 2 u(εi) j -i, j -i + o(ε 2 ). ( 20
)
We now notice that, since the missing terms sum up to a lower order term, we may sum up over all j ∈ \{i} in the terms of order ε 2 . We may in fact do the same thing for the term of order ε. Indeed, the missing terms may be estimated as follows:

ε 2 |j | A ∇W 0 ∇u(εi)j D 2 u(εi) j, j |j | A Cε|j | 2 |j | d+1+a = O ε A a-1 = o(ε 2 ), ( 21
)
from the definition of A. This estimate is valid for i far enough from the boundary of 1 ε , since otherwise the sum is truncated not at A, but at d(i, 1 ε c ). On the other hand, these terms are easily seen, through the same kind of estimate, to be boundary terms of higher order.

We leave the first term as it stands for now, dealing with it afterwards. The other terms are easily seen to be Riemann sums, converging, up to higher order boundary terms, to the corresponding integrals. We thus have:

E ε,δ (u) = 1 2N i∈ ∩ 1 ε 0<|i-j |<A W 0 ∇u(εi)(j -i) + ε 4| | j ∈ \{0} ∇W 0 ∇u(x)j D 2 u(x)(j, j ) dx + ε 2 1 6| | j ∈ \{0} ∇W 0 (∇u(x)j )D 3 u(x)(j, j, j )dx + 1 8| | j ∈ \{0} D 2 W 0 (∇u(x)j )(D 2 u(x)(j, j ), D 2 u(x)(j, j ))dx + ε 2 F 1 (u) + o(ε 2 ), ( 22
)
where F 1 (u) is a boundary term. Integrating by parts the first part of the term of order ε 2 , we easily see that it is exactly the corresponding term of [START_REF] Chipot | Equilibrium configurations of crystals[END_REF]. Turning to the term of order ε, we see that its integrand is an exact derivative with respect to x, so that it is equal to

j ∈ \{0} ∇W 0 ∇u(x)j D 2 u(x)(j, j ) dx = ∂ j ∈ \{0} W 0 (∇u(x)j )(j n)dσ (x). ( 23
)
Using the fact that W 0 (x) = W 0 (-x), together with the fact that the set \ {0} is symmetric with respect to 0, we see that this term cancels. We now deal with the term of order zero: it will naturally give the zero-order term of ( 15), and also approximating terms of order 1 and 2 a priori. The point here is to show that these terms give the boundary term of [START_REF] Chipot | Equilibrium configurations of crystals[END_REF]. Denoting by E 0 ε,δ (u) the zero-order term of ( 22), we write:

E 0 ε,δ (u) = 1 2N i∈ ∩ 1 ε j ∈ \{0} W 0 (∇u(εi)j )dx - 1 2N i∈ ∩ 1 ε j | ∈ 1 ε W 0 (∇u(εi)(j -i)) - 1 2N i∈ 1 ε j ∈ 1 ε ∩B A (i) c W 0 (∇u(εi)(j -i)). ( 24
)
Here, all the sums are restricted to points belonging to . The third term is easily bounded as follows:

1 2N i∈ ∩ 1 ε j ∈ 1 ε ∩B A (i) c W 0 (∇u(εi)(j -i)) C |j |>A 1 |j | a+d C A a = o(ε 2 ).
The two remaining terms will respectively give the terms of order zero and one in [START_REF] Chipot | Equilibrium configurations of crystals[END_REF]. In order to prove this claim for the first term, we denote by f the function

f (x) = j | =0 W 0 (∇u(x)j
), and write it as:

ε 3 2| | i∈ ∩ 1 ε f (εi) = 1 2| | f (x)dx + 1 2| | i∈ ∩ 1 ε εi+εQ( ) (f (εi) -f (x))dx, + i∈ ∩ 1 ε (εi+εQ( ))∩ c f - i∈ ∩ 1 ε c (εi+εQ( ))∩ f, ( 25
)
where Q( ) is the primitive unit cell of . We next make a Taylor expansion of the second term, writing f (εi

)-f (x) = ∇f (εi)(i-x)-D 2 f (εi)(x-i, x-i)+O(ε 3 ).
The first term of this expansion is linear with respect to xi, and thus cancels when integrated over i + Q( ). We thus have

ε 3 2| | i∈ ∩ 1 ε f (εi) = 1 | | f (x)dx - ε 5 2| | i∈ ∩ 1 ε Q( ) D 2 f (εi)(y, y)dy + O(ε 3 ).
We recognize here again in the second term a Riemann sum, so that it may be replaced, up to terms which are negligible before ε 2 , by the quantity

ε 2 Q( ) D 2 f (x)(y, y)dydx.
This term being an integral of a derivative with respect to x, it will give only a boundary term of order ε 2 . We now turn to the last terms of [START_REF] Gardner | The infinite-volume ground state of the Lennard-Jones potential[END_REF], which are boundary terms, since they clearly involve only terms such that d(i, 1 ε ∂ ) < 2ε. Next, we notice that each term of this difference is a Riemann sum of the boundary, and thus may be replaced, up to boundary terms of order ε 2 , by the corresponding integrals, which are equal.

We finally deal with the remaining term, that is, the second term of [START_REF] Friesecke | Validity and failure of the Cauchy-Born hypothesis in a 2D mass-spring lattice[END_REF]. Since the treatment of this term is rather delicate, we provide a proof in the case when is a polyhedron. The extension to the general case follows the same argument with some straightforward technical adaptations that we omit. The first point is that, since j is restricted to belong to 1 ε c , using once again that ∇u -1 is bounded uniformly on , we may reduce this sum to terms which are not far from the boundary:

- 1 2N i∈ ∩ 1 ε j | ∈ 1 ε W 0 ∇u(εi)(j -i) = - 1 2N i∈ ∩∂ ε j | ∈ 1 ε W 0 ∇u(εi)(j -i) + o(ε 2 ),
where

∂ ε = {y ∈ 1 ε , d(y, ∂( 1 ε )) < 1 √ ε }.
This expression shows that this is going to be a boundary term. Next, we separate the boundary of into P different faces, denoted by p , with 1 p P . We now consider only one face, putting all of them together in the end:

E p (u) = - 1 2N i∈ ∩ ε p j +i∈ ∩ 1 ε c W 0 ∇u(εi)j ,
where ε p is the set of points in 1 ε which are at a distance smaller than 1 √ ε from 1 ε p . Now, considering the set to which j belongs, it may be described, up to a term of order ε 2 , by the constraints j ∈ and (j + iπ 1 ε p (i))n 0, where n is the outer normal of the face p , and π 1 ε p is the orthogonal projection on the hyper-plane 1 ε p . Now, since d(i,

1 ε p ) = (π 1 ε p (i) -i)n, this property also reads jn d(i, 1 ε p ). Hence, E p (u) = - ε d 2| | i∈ ∩ ε p j ∈ , j n d(i, 1 ε p ) W 0 ∇u(εi)j + O(ε 2 ).
We next approximate, up to order ε, the sum over j by the same one where εi is replaced by π p (εi). This may be done using a Taylor expansion, and the correcting term is shown to be of order ε 2 when summed up with respect to i using the fact that a > 2. The term we want to treat now is thus:

E p (u) = - ε d 2| | i∈ ∩ ε p j ∈ , j n d(i, 1 ε p ) W 0 ∇u(π p (εi))j + O(ε 2 ).
Assume for a while that ∇u is constant and equal to F on p . We now define the set

A p ε ( p ) = {d(i, 1 ε p ), i ∈ ∩ ε p }, and N p ε (k) = #{i ∈ ε p ∩ , d(i, 1 ε
p ) = k}, and using this notation, we write the above energy as

E p (u) = - ε d 2| | k∈A p ε ( p ) N p ε (k) j ∈ , j n k W 0 Fj + O(ε 2 ).
It only remains to point out that the number [START_REF] Catto | Limite thermodynamique pour des modèles de type Thomas-Fermi[END_REF], up to a correcting term of order ε 2-d . Here | p | is the (d -1)dimensional measure of the face p . We therefore conclude, in this special case, that

N p ε (k) is equal to the number | p |N ε (k) appearing in
E p (u) = - ε| p | 2| | ∞ 0 j ∈ , j n k W 0 (∇u(x)j ) dµ p , (k),
which matches exactly the second term of [START_REF] Chipot | Equilibrium configurations of crystals[END_REF]. Here, x is any point of p . In order to finish the proof, we only need to point out that if ∇u is not a constant on p , a similar but more tedious analysis leads to the second term of (15) through a Riemann sum over p . Next, using an approximation of by a polyhedron, we conclude the proof in the general case.

We now explain how the above argument may be adapted to the case where a 8. Note that the points where we have used the assumption a > 8 concern only the bulk term, and not the surface terms. In the bulk term, the point is the choice of β, i.e., of A, the cut-off radius used in [START_REF] Cleri | Atomistic simulations of materials fracture and the link between atomic and continuum length scales[END_REF]. We need β > 2 a so that the remainder term in (17) be of order strictly lower than ε 2 . If a 8, this implies that the remainder of ( 18) is no more of order o(ε 2 ). Hence, we need to expand at an order higher than 2, say q. And the remainder, which is of order O( ε q A q+1 |j -i| d+a+1 ), leads, when summed up over j , to a term of order O(ε q-1-βq ). Hence, if β is chosen strictly lower than 1, which is possible since a > 2, this term is of order strictly lower than ε 2 for q large enough. Then we only have to deal with the additional terms of the development, which we do using their exact expression, as was done in [START_REF] Maso | An introduction to -convergence[END_REF] for the term of order three.

Note that, as is clear in the above proof, the development could be carried on at any order. However, even the surface term of order ε is a little cumbersome to establish, and its expression seems to be difficult to use except in special cases. This indicates that higher order surface terms are likely to be hard to compute and to use. On the contrary, if we only look for bulk terms, the development is within reach. (In such a case, the assumption ( 14) would involve a > q, q being the order at which we want to develop.) Another point is that the term of order ε is a boundary term, even if the assumption W 0 (x) = W 0 (-x) is dropped. However, in this case, the right-hand side of [START_REF] Friesecke | A scheme for the passage from atomic to continuum theory for thin films, nanotubes and nanorods[END_REF] would appear in addition to the term of [START_REF] Chipot | Equilibrium configurations of crystals[END_REF].

Finally, let us point out that this order-one term shows a minus sign because (and it is apparent in the above proof) it is the correction of the bulk term of order zero, to which we have added terms so as to have the whole sum over when x is close to the boundary. It may be therefore interpreted as the opposite of the surface energy of the crystal. However, this interpretation is subjected to the assumption that the atoms are distributed on the lattice even near the boundary. This assumption is highly questionable, since it is known that relaxation effects near the boundary may change or even destroy the periodicity locally. These effects might be of same order as the surface energy we have derived.

Thomas-Fermi type models

We deal in this Section with Thomas-Fermi type models, where the quantum nature of the electrons is (partially) taken into account, whereas the nuclei are supposed to be classical particles.

We first present the Thomas-Fermi-von Weizsäcker (TFW) model, define it for molecules, and then explain how it translates into TFW theory for infinite periodic (solid-state) systems through the thermodynamic limit process [START_REF] Catto | Mathematical Theory of thermodynamic limits : Thomas-Fermi type models[END_REF]. We also show the link between this thermodynamic limit process and the present homogenization scheme.

Presentation of the models

We briefly present in this subsection the Thomas-Fermi-von Weizsäcker (TFW) model, recalling that everything we are going to do for this model is clearly adaptable to Thomas-Fermi (TF), and even to Thomas-Fermi-Dirac-von Weizsäcker model. We refer to [START_REF] Lieb | Thomas-Fermi and related theories of atoms and molecules[END_REF] for details on these models.

The TFW energy is defined, for a set of M nuclei of positions {X i } 1 i M and charges {Z i } 1 i M , and a set of electrons defined by their total density ρ 0 (such

that √ ρ ∈ H 1 (R 3 )), by E TFW ({X i , Z i }, ρ) = h2 m R 3 |∇ √ ρ| 2 + R 3 ρ 5/3 + e 2 4πε 0 - M i=1 R 3 Z i ρ(x) |x -X i | dx + 1 2 R 3 ρ(x)ρ(y) |x -y| dxdy + 1 2 i | =j Z i Z j |X i -X j | , ( 26
)
where h is Planck's constant, e the elementary charge, m the mass of an electron and ε 0 the dielectric permittivity constant. Note that the TFW model is usually stated in a unit system where the coefficients h 2 m and e 2 4πε 0 of ( 26) are both equal to 1, but here we need to scale the characteristic length of the model with respect to the atomic spacing, as was done in Section 2 for two-body energies. Let us also point out that we have skipped here dimensionless constants which should appear in front of the first two terms of (26) (see [START_REF] Parr | Density-Functional Theory of Atoms and Molecules[END_REF] for the details), since they are mathematically irrelevant here.

The integer N being the total number of electrons, the density ρ is subjected to the constraint ρ = N . When the electrons are in their ground state, they minimize the above energy, i.e., they are a solution of the minimization problem

E TFW ({X i , Z i }) = inf E TFW ({X i , Z i }, ρ), ρ 0, √ ρ ∈ H 1 (R 3 ), R 3 ρ = N . ( 27
)
Here, we assume that we are dealing with a set of N identical atoms, and fix the nuclear charge to Z i = 1, although none of these assumptions are limitations. The important assumption is that the nuclei are periodically distributed, that is, as in Section 2, {X i } = ε ∩ , where is the periodic lattice on which the nuclei are distributed, ε the inter-atomic distance, and the Lipschitz open set defining the solid we are studying.

We next point out that the characteristic length of the system is easily shown by a dimensional analysis to be δ = h2 4πε 0 e 2 m . Hence, setting E 0 = h 2 m , [START_REF] Lieb | Thomas-Fermi and related theories of atoms and molecules[END_REF] translates into

E TFW δ ({X i }, ρ) =E 0 R 3 |∇ √ ρ| 2 + R 3 ρ 5/3 + 1 δ - M i=1 R 3 ρ(x) |x -X i | dx + 1 2 R 3 R 3 ρ(x)ρ(y) |x -y| dxdy + 1 2 i | =j 1 |X i -X j | , (28)
and ( 27) similarly. Assuming now that the set of nuclei {X i } = ε ∩ is deformed by a C ∞ diffeomorphism u, we thus define

E TFW ε,δ (u) = 1 N inf E TFW δ (u(ε ∩ ), ρ), ρ 0, √ ρ ∈ H 1 (R 3 ), R 3 ρ = N , ( 29
)
where N is the total number of nuclei (equal to the number of electrons), that is, 29) is the equivalent of ( 5) in the present case of TFW theory.

N = #(ε ∩ ). Equation (
We recall that given a set of any finite nuclei {X i } 1 i N , the minimization problem ( 27) has a unique solution ρ. We refer to [START_REF] Benguria | The Thomas-Fermi-von Weizsäcker theory of atoms and molecules[END_REF] for the proof of this result.

The function E TFW ε,δ (u) is thus well defined by (29).

Infinite periodic systems

The problem of the thermodynamic limit of the above (TFW) model is closely linked with the problem we are dealing with here. Indeed, setting ρ(x) = ε 3 ρ(εx) and changing variables in [START_REF] Lieb | The Hartree-Fock theory for Coulomb systems[END_REF], we get

E TFW δ (u(ε ∩ ), ρ) = E 0 1 ε 2 R 3 |∇ ρ| 2 + R 3 ρ5/3 + 1 δε - i∈ ∩ 1 ε R 3 ρ(x) |x -u(εi) ε | dx + 1 2 R 3 R 3 ρ(x) ρ(y) |x -y| dxdy + 1 2 i | =j ∈ ∩ 1 ε 1 | u(εi)-u(εj ) ε | . ( 30
)
In the special case where the function u is linear, u(εi) ε simplifies to u(i), so that the problem we address here is exactly the problem of finding the limit of the ground-state TFW energy as the set of nuclei fills in the lattice u( ), which is the thermodynamic limit problem for the TFW model, as dealt with in [START_REF] Catto | Mathematical Theory of thermodynamic limits : Thomas-Fermi type models[END_REF][START_REF] Catto | Limite thermodynamique pour des modèles de type Thomas-Fermi[END_REF], given the fact that the constants involved in (30) scale properly, namely 1 ε 2 ∼ 1 δε , i.e., δ ∼ ε, which is the equivalent of (i) of Theorem 1. We thus recall here the results of [START_REF] Catto | Mathematical Theory of thermodynamic limits : Thomas-Fermi type models[END_REF]: it is possible to define the (renormalized) TFW energy of an infinite periodic system in which the nuclei are distributed on a lattice and the electronic density ρ is such that ρ ∈ H 1 per ( ), the set of functions in H 1 loc which are -periodic, by:

E TFW ( ρ, ) = Q( ) |∇ ρ| 2 + Q( ) ρ5/3 - Q( ) ρG + 1 2 Q( ) Q( ) ρ(x)G (x -y) ρ(y)dxdy, (31)
where Q( ) is the primitive unit cell (or any unit cell) of the lattice , and G may be seen as the -periodic version of the Coulomb potential, and is defined by the following:

-G = 4π δ 0 -1 |Q( )| in Q( ), G is -periodic, lim |x|→0 G (x) - 1 |x| = 0. ( 32
)
We then have the following theorem (see [START_REF] Catto | Mathematical Theory of thermodynamic limits : Thomas-Fermi type models[END_REF]):

Theorem 4 (Thermodynamic limit of the TFW model, [START_REF] Catto | Mathematical Theory of thermodynamic limits : Thomas-Fermi type models[END_REF]).

Let be an open Lipschitz bounded subset of R 3 . Denote by E TFW ( ∩ 1 ε , ρ) the rescaled TFW energy, that is, [START_REF] Love | A treatise on the mathematical theory of elasticity[END_REF] with E 0 = ε 2 , δ = ε and u = Id:

E TFW ∩ 1 ε , ρ = R 3 |∇ ρ| 2 + R 3 ρ5/3 - i∈ ∩ 1 ε R 3 ρ(x) |x -i| dx + 1 2 R 3 R 3 ρ(x) ρ(y) |x -y| dxdy + 1 2 i | =j ∈ ∩ 1 ε 1 |i -j | . ( 33
)
And define

E TFW ∩ 1 ε = 1 N inf E ε ( ∩ 1 ε , ρ), ρ 0, ρ ∈ H 1 (R 3 ), R 3 ρ = N , with N = #( ∩ 1 ε ). Then, it follows that lim ε→0 E TFW ∩ 1 ε = E TFW ( ),
where E TFW ( ) is defined by the following minimization problem:

E TFW ( ) = inf E TFW ( , ρ), ρ 0, ρ ∈ H 1 per ( ), Q( ) ρ = 1 , ( 34
)
the periodic energy E TFW ( , ρ) being defined by [START_REF] Man | On the constitutive equations of some weakly-textured materials[END_REF].

In addition, the solution ρε of the minimization problem E TFW ( ∩ 1 ε ) converges to the solution of E TFW ( ), uniformly on any set of the form

( 1 ε -γ ε ) , where 1 γ ε 1 ε .
The proof of Theorem 4 as presented in Chapter 5 of [START_REF] Catto | Mathematical Theory of thermodynamic limits : Thomas-Fermi type models[END_REF] is based on the convergence and uniqueness of the solution of the Euler-Lagrange equation of minimization problems E TFW ( ∩ 1 ε ) and E TFW ( ). The first one converges in some sense to the second, so that their solutions also converge. This is why we now write down the Euler-Lagrange equation of the minimization problem defining E TFW ε,δ (u) [START_REF] Lions | Solutions of Hartree-Fock equations for Coulomb systems[END_REF]. Setting v = ρ, it reads:

-v + 5 3 v 7/3 + v 2 * 1 |x| - i∈ ∩ 1 ε 1 |x -i| + θ v = 0,
where θ is the Lagrange multiplier associated with the constraint ρ = N . Setting

φ = -v 2 * 1 |x| + i∈ ∩ 1 ε 1
|x-i|θ, we thus have:

-v + 5 3 v 7/3 -φv = 0, -φ = 4π i∈ ∩ 1 ε δ i -v 2 . ( 35
)
In the periodic case, the sum over ∩ 1 ε is replaced by a sum over the whole of . Note that as ε goes to 0, the system [START_REF] Ortiz | Quasicontinuum analysis of defects in solids[END_REF] converges to the periodic one. Our aim in the following subsection is to adapt the proof of Theorem 4 to the case when u is not a linear transformation, but a general C ∞ diffeomorphism.

Convergence theorem

We give in this subsection the equivalent of Theorem 1, together with its proof: Theorem 5. Let be a Lipschitz open set of R 3 , and let u be a C ∞ diffeomeorphism defined on , and a periodic lattice. Suppose that the volume of the primitive unit cell of is normalized, i.e |Q( )| = 1. Consider the energy E ε,δ (u) defined by [START_REF] Lieb | The Hartree-Fock theory for Coulomb systems[END_REF], [START_REF] Lions | Solutions of Hartree-Fock equations for Coulomb systems[END_REF]. Then, the following statements hold:

(i) If ε = δ, and if E 0 = ε 2 ,
then there is a convergence:

lim ε→0 E ε,δ (u) = 1 | | E TFW (∇u(x) )dx, ( 36
)
where E TFW is the rescaled TFW energy defined by [START_REF] Nijboer | On the configuration of systems of interacting particles with minimum potential energy per particle[END_REF], [START_REF] Man | On the constitutive equations of some weakly-textured materials[END_REF], [START_REF] Üller | Variational models for microstructure and phase transitions[END_REF]. (ii) If ε δ, and if E 0 = ε 2 , then there is a convergence:

lim ε→0 E ε,δ (u) = 1 | | dx | det(∇u(x))| 2/3 . ( 37
)
(iii) If δ ε, and if E 0 = δ 2 , then E ε,δ (u) converges to a constant independent of u. This constant may be identified as the rescaled atomic TFW energy [START_REF] Lieb | The Thomas-Fermi theory of atoms, molecules and solids[END_REF].

Proof. The ingredients of the present proof are mainly present in [START_REF] Catto | Mathematical Theory of thermodynamic limits : Thomas-Fermi type models[END_REF]. We nevertheless provide them for the sake of completeness. We start with the proof of (i). We denote by ρ the solution of the minimization problem ( 28), ( 29) defining E ε,δ (u), and by ρ the rescaled electronic density, that is, ρ(x) = ε 3 ρ(εx). The proof will be carried out in three steps:

Step 1. We have the convergence

lim ε→0 ρ • + 1 ε u(εi) -ρ ∇u(εi) L ∞ (∇u(εi)Q( )) = 0, (38)
which is uniform with respect to i ∈ ∩

( 1 ε -1 √ ε ) .
Here, for any periodic lattice , the density ρ is the corresponding TFW electronic density, i.e., the solution of ( 34), [START_REF] Man | On the constitutive equations of some weakly-textured materials[END_REF], [START_REF] Üller | Variational models for microstructure and phase transitions[END_REF].

In order to show [START_REF] Paroni | Constitutive equations of elastic polycrystalline materials[END_REF], we argue by contradiction (this proof is an adaptation of that of Theorem 5.9 in [START_REF] Catto | Mathematical Theory of thermodynamic limits : Thomas-Fermi type models[END_REF]). Assuming that it does not hold, we deduce the existence of some

x ε ∈ ∇u(εi)Q( ) such that v(x ε + u(εi) ε ) -v ∇u(εi) (x ε ) α > 0
for some α independent of ε, where v = ρ, and v ∇u(εi) = √ ρ ∇u(εi) . Using the regularity of u, we point out that the sequence x ε is compact, so that, up to a subsequence, we may assume that it converges to some

x 0 ∈ 0<ε 1 i∈ ∩ 1 ε ∇u(εi)Q( ).
We recall now that v satisfies the Euler-Lagrange equation of the minimization problem ( 28), [START_REF] Lions | Solutions of Hartree-Fock equations for Coulomb systems[END_REF], that is, the equivalent of ( 35):

-v + 5 3 v 7/3 -φv = 0, -φ = 4π j ∈ ∩ 1 ε δ u(εj ) ε -v 2 , ( 39
)
where the effective potential φ may also be defined by

φ = j ∈ ∩ 1 ε 1 |x -u(εj ) ε | -v 2 * 1 |x| -θ,
the constant θ being the Lagrange multiplier associated with the constraint v 2 = N in [START_REF] Lions | Solutions of Hartree-Fock equations for Coulomb systems[END_REF]. Using elliptic regularity results, it is then possible, using the method of [12] (Propositions 3.8 and 3.12), to show that there exists a constant C independent of ε such that

v L ∞ (R 3 ) + φ L p unif (R 3 ) + φ L ∞ ( 1 ε c ) C, ( 40
)
for any p < 3. Inserting this information in the first equation of [START_REF] Parr | Density-Functional Theory of Atoms and Molecules[END_REF], this also implies that v is uniformly continuous on R 3 , uniformly with respect to ε. Hence, for a sufficiently small ε, we have

v x 0 + u(εi) ε -v ∇u(εi) (x 0 ) α 2 . ( 41
)
Now, the bounds we have on v and φ are valid for v(• + u(εi) ε ) and φ(• + u(εi) ε ), and thus allow us to assume, using elliptic regularity and Rellich's theorem, that they converge, up to a subsequence once more, to some v and φ in L ∞ loc (R 3 ) and L p loc (R 3 ), respectively. We also point out that εi belongs to , which is compact, so that we may assume that it converges to some y ∈ . We thus have:

v • + u(εi) ε -→ v in L ∞ loc (R 3 ), φ • + u(εi) ε -→ φ in L p loc (R 3 ), εi -→ y ∈ .
We may then pass to the limit in the system satisfied by v

•+ u(εi) ε and φ •+ u(εi) ε , that is, -v • + u(εi) ε + - 5 3 v • + u(εi) ε 7/3 -φv • + u(εi) ε = 0, -φ • + u(εi) ε = j ∈ ∩ 1 ε δ u(εj ) ε -u(εi) ε -v • + u(εi) ε 2 .
Note here that the measure

j ∈ ∩ 1 ε δ u(εj ) ε -u(εi) ε converges in D (R 3 ) to the measure j ∈ δ ∇u(y)j , because lim ε→0 d(i, 1 ε c ) = +∞.
Therefore, we have:

-v + 5 3 v 7/3 -φ -v = 0, -φ = 4π j ∈ δ ∇u(y)j -v 2 .
This system has a unique solution (v ∇u(y) , φ ∇u(y) ) in L ∞ (R 3 ) × L 1 unif (R 3 ), according to Theorem 6.5 of [START_REF] Catto | Mathematical Theory of thermodynamic limits : Thomas-Fermi type models[END_REF], and therefore we should have v = v ∇u(y) . We now reach a contradiction with (41) if we can pass to the limit in [START_REF] Radin | Classical ground states in one dimension[END_REF], that is, if v M is a continuous function of the matrix M. This result is easily shown by repeating the same argument as above. This completes the proof of (38).

Step 2. We have the convergence

lim ε→0 φ • + u(εi) ε -φ ∇u(εi) L ∞ (∇u(εi)Q( )) = 0, (42)
which is uniform with respect to i ∈

( 1 ε -1 √ ε ) .
Here, for any lattice , φ is the effective potential associated with the density ρ , solution of (34), [START_REF] Man | On the constitutive equations of some weakly-textured materials[END_REF], [START_REF] Üller | Variational models for microstructure and phase transitions[END_REF]. In other words, v = √ ρ and φ are the unique solutions in L ∞ (R 3 ) and L 1 unif (R 3 ) respectively, of the system

-v + 5 3 v 7/3 -φv = 0, -φ = 4π j ∈ -δ j -v 2 .
We skip the proof of [START_REF] Radin | Ground states for soft disks[END_REF], since it is an easy adaptation of the one of [START_REF] Paroni | Constitutive equations of elastic polycrystalline materials[END_REF].

Step 3: Convergence of the energy. We are now in position to show the desired convergence result. In order to do so, we split the expression of the energy into different terms, which we treat separately. We first show:

E 1 ( ρ) := 1 N R 3 ρ5/3 -→ 1 | | ∇u(x)Q( ) ρ 5/3 ∇u(x) (z)dz dx, ( 43
)
as ε → 0. We separate this integral into a sum of integrals over domains which will, in the end, converge to the unit cell ∇u(x)Q( ) = Q(∇u(x) ) :

E 1 ( ρ) := 1 N j ∈ ∩ 1 ε 1 ε u(εj +εQ( )) ρ5/3 + 1 N j ∈ ∩ 1 ε c 1 ε u(εj +εQ( ))
ρ5/3 .

We begin by dealing with the first sum. We have

1 ε u(εj +εQ( )) ρ5/3 = 1 ε (u(εj +εQ( ))-u(εj )) ρ5/3 (z + 1 ε u(εj ))dz = ∇u(εj )Q( ) ρ 5/3 ∇u(εj ) (z)dz + o(1),
uniformly with respect to j ∈ ( 1 ε -1 √ ε ) according to [START_REF] Paroni | Constitutive equations of elastic polycrystalline materials[END_REF] and the fact that ∇u(εj

)Q( ) \ 1 ε u(εj + εQ( )) -u(εj ) + 1 ε u(εj + εQ( )) -u(εj ) \ ∇u(εj )Q( ) converges to 0 uniformly with respect to j ∈ 1 ε -1 √ ε
. Hence, using the fact that ρ is bounded independently of ε to show that boundary terms

j ∈ 1 ε \ ( 1 ε -1 √ ε
) make a negligible contribution, we deduce that

E 1 ( ρ) = 1 N j ∈ ∩ 1 ε ∇u(εj )Q( ) ρ 5/3 ∇u(εj ) + 1 N j ∈ ∩ 1 ε c 1 ε u(εj +εQ( )) ρ5/3 + o(1).
The first sum may be identified as the Riemann sum converging to the desired integral. Hence, in order to conclude the proof of ( 43), we only need to show that the second sum converges to 0. In order to do so, we point out that the same argument as above allows us to show that

1 = 1 N R 3 ρ = 1 N j ∈ ∩ 1 ε ∇u(εj )Q( ) ρ ∇u(εj ) + 1 N j ∈ ∩ 1 ε c ∇u(εj )Q( ) ρ + o(1) = 1 + 1 ε c ρ + o(1), so that lim ε→0 1 ε c ρ = 0. ( 44
)
This, together with the fact that ρ is uniformly bounded in R 3 , allows us to show that 1 ε c ρ5/3 = o(1), concluding the proof of [START_REF] Radin | Periodicity of classical ground states[END_REF]. We now turn to the second term of the energy, and show the following convergence:

E 2 ( ρ) := 1 N R 3 |∇ ρ| 2 -→ 1 | | ∇u(x)Q( ) |∇ √ ρ ∇u(x) | 2 dx. ( 45
)
We proceed here in the same way as for [START_REF] Radin | Periodicity of classical ground states[END_REF], showing first that the exterior contribution is negligible. In order to do so, we recall that, using the same notation as in

Step 1, we have

-v + 5 3 v 7/3 -φv = 0,
so that, multiplying this equation by v and integrating over 1 ε c , we have

1 N 1 ε c (-v)v = 1 N 1 ε c φ ρ - 1 N 1 ε c 5 3 ρ5/3 = o(1),
according to ( 40) and [START_REF]Quasicrystals and geometry[END_REF]. Next, we use elliptic regularity to show that

1 N ∂( 1 ε ) v ∂v ∂n = o(1), thereby proving that lim ε→0 1 N 1 ε c |∇v| 2 = 0.
Now, using elliptic regularity here again, it can easily be seen that ∇v is bounded in L 2 unif (R 3 ), so that

E 2 ( ρ) = 1 N j ∈ ∩( 1 ε -1 √ ε ) 1 ε u(εj +εQ( )) |∇v| 2 + o(1).
Following the proof of ( 43), what we need here in order to conclude the proof is the following convergence result:

∇v -∇v ∇u(εi ) L 2 (∇u(εi)Q( )) -→ 0 as ε → 0, uniformly with respect to i ∈ ( 1 ε -1 √ ε )
. This result is easily proved using ( 38) and ( 42) and the first equation of [START_REF] Parr | Density-Functional Theory of Atoms and Molecules[END_REF]. This allows to conclude the proof of [START_REF] Stave | The structure of Ni N and P d N clusters[END_REF].

We finally deal with the electrostatic terms of the energy, which is more intricate than ( 43) and [START_REF] Stave | The structure of Ni N and P d N clusters[END_REF]. We recall that the electrostatic energy reads:

1 N - i∈ ∩ 1 ε R 3 ρ(x) |x -u(εi) ε | dx + 1 2 R 3 R 3 ρ(x) ρ(y) |x -y| dxdy + 1 2 i | =j ∈ ∩ 1 ε ε |u(εi) -u(εj )| = - 1 2N R 3 φ ρ + 1 2N i∈ ∩ 1 ε lim x→ u(εi) ε φ(x) - 1 |x -u(εi) ε | . ( 46
)
Here we have used the definition of φ, that is,

φ = i∈ ∩ 1 ε 1 |x-u(εi) ε | -ρ * 1 |x| -θ.
We may deal separately with each term of the right-hand side of [START_REF] Ventevogel | On the configuration of a one-dimensional system of interacting particles with minimum potential energy per particle[END_REF]. The first term is easily treated using the same arguments as for [START_REF] Radin | Periodicity of classical ground states[END_REF], based on ( 38), ( 42) and ( 40), showing: [START_REF] Wales | What can calculations employing empirical potentials teach us about bare transition metal clusters?[END_REF] as ε → 0. We now deal with the last term of [START_REF] Ventevogel | On the configuration of a one-dimensional system of interacting particles with minimum potential energy per particle[END_REF], and show the following:

E 3 ( ρ) := - 1 2N R 3 φ ρ -→ - 1 2| | ∇u(x)Q( ) φ ∇u(x) ρ ∇u(x) dx
E 4 ( ρ) := 1 2N i∈ ∩ 1 ε lim y→ u(εi) ε φ(y) - 1 |y -u(εi) ε | -→ 1 2| | lim y→0 φ ∇u(x) (y) - 1 |y| dx (48)
as ε → 0. Here again, we copy the proof of ( 43), separating into terms i ∈ ( 1 ε -

1 √ ε ) and those belonging to 1 ε \ ( 1 ε -1 √ ε )
. For the first part, the use of ( 42) allows us to conclude. In order to show that the remaining terms are negigible, i.e., that

1 N i∈ ∩( 1 ε \ 1 √ ε ) lim y→0 φ(y + u(εi) ε ) - 1 |y| = o(1), ( 49
)
we use the second equation of [START_REF] Parr | Density-Functional Theory of Atoms and Molecules[END_REF], and (40), getting

-φ(x + u(εi) ε ) - 1 |x| C
on some ball B α , with α > 0 and C independent of i and ε, provided ε is sufficiently small. Hence, using the fact that φ is bounded in L 1 unif (R 3 ) and the mean value inequality, this implies that

φ(x + u(εi) ε ) - 1 |x| C ∀x ∈ B α .
This estimate allows us to conclude that (49) holds, ending the proof of ( 48). Now we only have to collect ( 43), ( 45), ( 47) and ( 48), and point out that, for all x ∈ ,

- ∇u(x)Q( ) G ∇u(x)Q( ) ρ ∇u(x) + 1 2 ∇u(x)Q( ) ∇u(x)Q( ) ρ ∇u(x) (z)G ∇u(x) (z -y)ρ ∇u(x) (y)dydz = - 1 2 ∇u(x)Q( ) φ ∇u(x) ρ ∇u(x) + 1 2 lim y→0 φ ∇u(x) (y) - 1 |y| ,
to conclude the proof of (i).

We now prove (ii). Here again, we first rescale the density ρ, setting ρ(x) = ε 3 ρ(εx). The energy may then be written as:

E ε,δ (u) = 1 N R 3 |∇ ρ| 2 + R 3 ρ5/3 + ε δ - j ∈ ∩ 1 ε R 3 ρ(x) |x -u(εj ) ε | + 1 2 R 3 R 3 ρ(x) ρ(y) |x -y| dxdy + 1 2 i | =j ∈ ∩ 1 ε ε |u(εi) -u(εj )| .
Next, we follow exactly the steps of the proof of (i), except that here, ε and δ appear as parameters in the system satisfied by v = ρ and φ

= j ∈ ∩ 1 ε 1 |x-u(εj ) ε | -v 2 * 1 |x| + δ ε θ.
We give the main ingredients, skipping their proofs since they involve the same arguments as above:

v L ∞ (R 3 ) + ε δ φ L p unif (R 3 ) + φ L ∞ ( 1 ε c ) C,
for any p < 3, and

lim ε δ→0 ρ • + u(εi) ε - ε |u(εi + εQ( ))| L ∞ (∇u(εi)Q( )) = 0, lim ε δ→0 ε δ φ • + u(εi) ε - ε 2/3 |u(εi + εQ( ))| 2/3 L 2 (∇u(εi)Q( )) = 0, uniformly with respect to i ∈ ∩ ( 1 ε -1 √ ε )
. With these convergence results, we can then easily show that all terms of the energy converge to 0 except the term ρ5/3 , which converges to the desired quantity.

Let us now prove (iii). Here again, it is possible to study the convergence of the density and deduce from it the convergence of the energy. But we will provide an alternate and more direct proof using the fact that we are dealing with a minimization problem. Once again we rescale ρ, but here the rescaling parameter will be δ instead of ε. We thus set ρ(x) = δ 3 ρ(δx). Then, the energy reads:

E ε,δ (u) = 1 N R 3 |∇ ρ| 2 + R 3 ρ5/3 - j ∈ ∩ 1 ε R 3 ρ(x) |x -u(εj ) δ | + 1 2 R 3 R 3 ρ(x) ρ(y) |x -y| dxdy + 1 2 i | =j ∈ ∩ 1 ε δ |u(εi) -u(εj )| = 1 N E TFW u(ε ∩ ) δ ,
the term E TFW u(ε ∩ ) δ denoting the solution of the minimization problem

E TFW u(ε ∩ ) δ = inf E TFW ( u(ε ∩ ) δ , ρ), ρ 0, ρ ∈ H 1 (R 3 ), R 3 ρ = N , ( 50
)
the energy E TFW u(ε ∩ ) δ , ρ being the rescaled TFW energy [START_REF] Nijboer | On the configuration of systems of interacting particles with minimum potential energy per particle[END_REF]. Let now ρ 0 be a non-negative, radial, compactly supported, smooth function of total mass 1. Using the function η(x) = j ∈ u(ε ∩ ) δ ρ 0 (xj) as a test function for the minimization problem (50), we have:

E ε,δ (u) 1 N E TFW u(ε ∩ ) δ , η .
We now point out that, for ε δ large enough, the terms of the sum defining η have disjoint supports, so that we have

R 3 |∇ √ η| 2 + R 3 η 5/3 = N R 3 |∇ √ ρ 0 | 2 + N R 3 ρ 5/3 0 .
Next, in order to compute the electrostatic terms, we out that since ρ 0 is radially symmetric, has compact support and total mass 1, the function ρ 0 * 1 |x| -1 |x| has its support included in the support of ρ 0 . Therefore, we have, if ε δ is large enough,

- 1 2 j ∈ ∩ 1 ε R 3 η(x) |x -u(εj ) δ | dx + 1 2 R 3 R 3 η(x)η(y) |x -y| dxdy = - 1 2 R 3 η(x) j ∈ ∩ 1 ε 1 |x -u(εj ) δ | -ρ 0 * 1 |x| x - u(εj ) δ dx = N 2 R 3 ρ 0 (x) - 1 |x| + ρ 0 * 1 |x| dx.
Similarly,

- 1 2 j ∈ ∩ 1 ε R 3 η(x) |x -u(εj ) δ | dx + 1 2 i | =j ∈ ∩ 1 ε δ |u(εi) -u(εj )| = - N 2 R 3 ρ 0 (x) |x| dx.
Collecting these computations, we thus have:

E ε,δ (u) R 3 |∇ √ ρ 0 | 2 + R 3 ρ 5/3 0 - R 3 ρ 0 (x) |x| dx + 1 2 R 3 R 3 ρ 0 (x)ρ 0 (y) |x -y| dxdy (51)
for any smooth radially symmetric ρ 0 0 having compact support. We now define the TFW atomic energy:

I TFW at (λ) = inf E TFW at (ρ), ρ 0, √ ρ ∈ H 1 (R 3 ), R 3 ρ = λ , ( 52)
with

E TFW at (ρ) = R 3 |∇ √ ρ| 2 + R 3 ρ 5/3 - R 3 ρ(x) |x| dx + 1 2 R 3 R 3 ρ(x)ρ(y) |x -y| dxdy. ( 53
)
This minimization problem has been studied in [START_REF] Benguria | The Thomas-Fermi-von Weizsäcker theory of atoms and molecules[END_REF], where it is shown to have a unique minimizer for all λ λ c , for some λ c > 1, which is radially symmetric. Hence, in (51), we may use a sequence of functions ρ 0 converging to the unique minimizer of (52), getting lim sup

δ ε→0 E ε,δ (u) I TFW at (1). ( 54
)
It only remains for us to prove the reverse inequality. In order to do so, we come back to the function ρ which achieves the minimum defining E ε,δ (u), and the corresponding rescaled density ρ = δ ρ(δx). Setting v = ρ, we write down the Euler-Lagrange equation satisfied by v, that is,

-v + 5 3 v 7/3 -φv = 0, -φ = 4π j ∈ ∩ 1 ε δ u(εj ) δ -v 2 .
Using this system of equations, it is possible to adapt the method introduced in [START_REF] Catto | Mathematical Theory of thermodynamic limits : Thomas-Fermi type models[END_REF] (Propositions 3.8, 3.10 and 3.12) in order to have

v L ∞ (R 3 ) + φ L p unif (R 3 ) C,
for any p < 3, with C independent of ε and δ. Using elliptic regularity and Rellich's theorem, we may thus assume that for all j ∈ ∩ 1 ε , ρ(• + u(εj ) δ ) converges weakly in H 1 loc (R 3 ) and strongly in L 2 loc (R 3 ) to some v j = √ ρ j . Fixing a radius R > 0, we have, for ε δ large enough:

lim inf δ ε→0 B √ ε δ ( u(εj ) δ ) |∇ ρ| 2 lim inf δ ε→0 B R ( u(εj ) δ ) |∇ ρ| 2 B R |∇ √ ρ j | 2 ,
and similarly for ρ 5/3 . Hence, letting R go to infinity, we have lim inf

δ ε→0 1 N R 3 |∇ ρ| 2 + R 3 ρ 5/3 lim inf ε→0 1 N j ∈ ∩ 1 ε R 3 |∇ √ ρ j | 2 + R 3 ρ 5/3 j .
Next, we use the bounds we have on φ and ρ, finding that φ -

1 |x-u(εi) δ | is bounded in L ∞ B √ ε δ , so that φ x + u(εj ) δ -1 |x| converges in L ∞ loc (R 3 ) to some φ j -1 |x| ∈ L ∞ (R 3
) such that φ j satisfiesφ j = δ 0ρ j . We now point out that the method used in [START_REF] Blanc | Geometry optimization for crystals in Thomas-Fermi type theories of solids[END_REF] (Propositions 4.2 and 4.3) is easily adapted to the present case, and yields:

max φ(x), |v(x)| 4/3 C |x -u(εi) δ | 2 + j ∈ ∩ 1 ε C |x -u(εj ) δ | 4 , for all x such that |x -u(εi) δ | inf j ∈ ∩ 1 ε \{i} |x -u(εj ) δ |.
This allows us to pass to the limit in the termρφ, getting lim inf

δ ε→0 1 N - R 3 ρφ = lim inf δ ε→0 1 N j ∈ ∩ 1 ε - B √ δ ε ( u(εj ) δ ) ρφ = lim inf ε→0 1 N j ∈ ∩ 1 ε - R 3 ρ j φ j . .
Finally, it is easy to adapt the proof of Step 2 of (i) in order to show that the convergence of φ(x

+ u(εj ) δ ) -1 |x| to φ j -1 |x| is uniform with respect to j ∈ ∩ 1 ε -1 √ ε , getting lim 1 N j ∈ ∩ 1 ε lim x→0 φ x + u(εj ) δ - 1 |x| = lim 1 N j ∈ ∩ 1 ε lim x→0 φ j (x) - 1 |x| .
Hence, pointing out that φ j = 1 |x| -ρ j * 1 |x| +θ j , for some constant θ j , we deduce that

-1 2 R 3 ρ j φ j + 1 2 lim x→0 (φ j -1 |x| ) =-ρ(x) |x| dx + 1 2 ρ(x)ρ(y)
|x-y| dxdy. Gathering all these results, we have lim inf

δ ε→0 E ε,δ (u) lim inf ε→0 1 N j ∈ ∩ 1 ε E TFW at (ρ j ) lim inf ε→0 1 N j ∈ ∩ 1 ε I TFW at (λ j ), with λ j = ρ j . Now, the function I TFW at (λ) is convex non-increasing with respect to λ ∈ R + . Since in addition 1 N λ j 1, we conclude that lim inf δ ε→0 E ε,δ (u) I TFW at (1). ( 55
)
This concludes the proof of (iii).

Remark 6.

Let us point out that the method used in the proof of (iii), based on variational methods, could be used in the proof of (i) and (ii) as well. However, although it seems more natural, it would be considerably more delicate.

Possible extensions

We indicate in this section some direct consequences of the above results, and natural extensions of them. There are mainly two types of extensions: those concerning the microscopic model (the electronic problem), and those concerning the geometry of the atoms.

Changing the model describing the electrons

In Section 3, we have dealt with Thomas-Fermi type models. These theories are very crude compared to the Schrödinger equation, which they are supposed to approximate. (The approximation of two-body potentials, used in Section 2, is even worse.)

However, the methods used in this section are intimately linked with the thermodynamic limit problem [START_REF] Catto | Mathematical Theory of thermodynamic limits : Thomas-Fermi type models[END_REF][START_REF] Lieb | The Thomas-Fermi theory of atoms, molecules and solids[END_REF]. More precisely, it seems that the existence of a thermodynamic limit allows us to pass to the limit ε = δ → 0, at least in the two cases of two-body potentials and Thomas-Fermi type models. Hence, formulas [START_REF] Blanc | Geometry optimization for crystals in Thomas-Fermi type theories of solids[END_REF] and ( 36) are likely to be adaptable to other models, such as for instance Hartree-Fock models [START_REF] Lieb | The Hartree-Fock theory for Coulomb systems[END_REF][START_REF] Lions | Solutions of Hartree-Fock equations for Coulomb systems[END_REF]. In this case, although the thermodynamic limit has not yet been fully justified, it is possible (see [START_REF] Catto | On the thermodynamic limit for Hartree-Fock type models[END_REF]) to derive it formally. Therefore, the elastic energy

E(u) = 1 | | E HF (∇u(x) )dx,
where E HF ( ) is the Hartree-Fock energy of the lattice , defined in [START_REF] Catto | On the thermodynamic limit for Hartree-Fock type models[END_REF], is a good candidate for the elastic Hartree-Fock energy. The same remarks hold for any quantum model.

In the case of the true Schrödinger equation, the difficulty is, as far as we know, that there is no derivation, even formally, of any thermodynamic limit of the model. More precisely, it is clearly possible to derive an energy functional, but the associated variational space is not so easy to guess [START_REF] Coleman | Reduced density matrices[END_REF].

Let us make a final remark about these quantum models: in the case of the Thomas-Fermi model, that is, when we forget the term |∇ √ ρ| 2 in the energy [START_REF] Lieb | Thomas-Fermi and related theories of atoms and molecules[END_REF], we have:

E TF ({X i }, ρ) = h2 m R 3 ρ 5/3 + e 2 4πε 0 - M j =1 R 3 ρ(x) |x -X j | dx + 1 2 R 3 R 3 ρ(x)ρ(y) |x -y| dxdy + 1 2 i | =j 1 |X i -X j | , ( 56
)
where we have used the same notation as in [START_REF] Lieb | Thomas-Fermi and related theories of atoms and molecules[END_REF], setting the nuclear charges Z j to 1. Let us for a while forget the physical constant, and assume that h2 m = 1 and e 2 4πε 0 = 1 in the system of units we use. A scaling argument (see [START_REF] Parr | Density-Functional Theory of Atoms and Molecules[END_REF]) easily shows that the power 5 3 appearing in (56) is in fact equal to d+2 d , where d is the dimension of the space (three, in our case). Therefore, if we consider the corresponding twodimensional model, the power we should use is 2, so that in this case (56) reads:

E TF ({X i }, ρ) = R 2 ρ 2 + M j =1 R 2 ρ(x) log(|x -X j |)dx - 1 2 R 2 R 2 ρ(x) log(|x -y|)ρ(y)dxdy - 1 2 i | =j log(|X i -X j |). (57)
Note that, dealing with a two-dimensional model, we have replaced the threedimensional Coulomb potential 1 |x| by the two-dimensional one, namelylog(|x|). The energy being quadratic with respect to ρ, the Euler-Lagrange equation of the problem

E TF ({X i }) = inf E TF ({X i }, ρ), ρ 0, ρ ∈ L 1 (R 2 ) ∩ L 2 (R 2 ), log(2 + |x|)ρ ∈ L 1 (R 2 ), R 2 ρ = N (58)
is linear with respect to ρ. Assuming neutrality, that is, N = M, it reads 2ρ

+ N j =1 log(|x -X i |) -ρ * log(|x|) + θ = 0, ( 59
)
where θ is the Lagrange multiplier associated with the mass constraint. Now, taking the Laplacian of this equation, we have

-ρ + 1 2 ρ = 2π N j =1 δ X j , ( 60
) with ρ ∈ L 1 (R 2
). This equation is easily solved, using the Yukawa potential in dimension two, that is the solution W of -W + 1 2 W = 2πδ 0 going to zero at infinity. The potential W is nothing else, in fact, than K 0 ( |x| √ 2

), where K 0 is the Bessel function of the second kind as defined in [START_REF] Abramowitz | Handbook of mathematical functions[END_REF]. Hence, we have the equality

ρ(x) = N j =1 W (x -X j ). ( 61
)
We now go back to the expression of the energy (56), and using (59) and (61), we have

E TF ({X i }, ρ) = R 2 ρ 2 + 1 2 R 2 ρ(x) N i=1 log |x -X i | -ρ * log |x| dx + 1 2 N j =1 R 2 ρ(x) log |x -X j |dx - i | =j log |X i -X j | = R 2 ρ 2 + 1 2 R 2 ρ(-2ρ -θ) + 1 2 N j =1 2ρ + θ + log(| • -X j |) (X j ) = N j =1 ρ + 1 2 log(| • -X j |) (X j ) = i | =j W (X i -X j ) + N lim x→0 W (x) + 1 2 log(|x|) .
Hence, up to an additive constant, the energy (58) may be expressed in terms of the two-body potential W : the two-dimensional TF model may be recast into a twobody model, and therefore enters the scope of Section 2. Incidentally, although the above computations seem rather basic, we have not found them in the literature. Note that, according to the proof of (ii) of Theorem 5, the link between the power p appearing in [START_REF] Paroni | Homogenization of polycrystalline aggregates[END_REF] and the power q appearing in the microscopic energy ( 26) is p = q -1, which ensures that in the present case formulas ( 8) and (37) become equivalent.

Changing the microscopic geometry of the atoms

We now make a few remarks about the microscopic arrangement of the atoms: so far, we have assumed that they are periodically distributed. Since this assumption is not always physically satisfactory, the same problem should be addressed in some other cases.

Our first point is a direct improvement of the preceding sections, and is concerned with polycrystalline materials. In this type of solid, we have a mix of different lattices 1 , 2 , . . . K , with volume ratios a 1 , a 2 , . . . , a K . The characteristic length of this mixing is far larger than the atomic one, and far smaller than the macroscopic one. We refer the interested reader to [START_REF] Man | On the constitutive equations of some weakly-textured materials[END_REF][START_REF] Paroni | Homogenization of polycrystalline aggregates[END_REF][START_REF] Paroni | Constitutive equations of elastic polycrystalline materials[END_REF] and the references therein. Therefore, introducing an intermediate scale γ such that ε γ 1 (we deal here only with the case ε = δ), we need to introduce a tiling of size γ of the set , for instance the unit cells of the lattice γ Z 3 , setting Q j = ∩ (γ Q + γj ) for j ∈ Z 3 , and Q being the unit cube. Then, separating each Q j into K sets Q 1 j , Q 2 j , . . . , Q K j of volume ratios a 1 , a 2 , . . . , a K respectively and assuming that in the set Q k j , the atoms are distributed on the set k , we thus may use exactly the same computations as those of Theorem 1(i) or Theorem 5(i), getting as elastic energy:

E(u) = 1 | | K k=1 a k E micro (∇u(x) k )dx, ( 62
)
where the energy functional E micro denotes the corresponding rescaled microscopic energy. We can also allow the coefficients a k to depend on x without any change.

Another possible improvement consists in replacing the measure K k=1 a k δ k implicitly used in (62) by any probability measure µ defined in the space of lattices L 3 (R 3 ). In this case, the Krein-Milmann theorem allows us to approximate µ by a sum of Dirac masses as above. The integer K becoming a parameter depending on ε, and going to infinity as ε → 0, with the condition ε γ

1 K
1, it is here again possible to adapt our method, finding: Let us point out that the result does not depend on the tiling we choose, as long as it is sufficiently regular (one can for instance replace the lattice Z 3 by any other lattice in the above argument).

E(u) = 1 | | L 3 (R 3 ) E micro (∇u(x) )dµ x ( )dx. ( 63
The point here is in fact that everything we did in the preceding sections extends to this situation, except that we have created interfaces between the grains of each lattice (the cubes Q k j in the present case). Now, the condition ε γ implies that the bulk in each Q j is far more important energetically than these interfaces. However, this remark allows us to predict that this presence of different phases should trigger additive terms in the higher order expansion (Theorem 3), and in particular a bulk term of order one accounting for these grain boundaries. Worst of all, this term will probably exhibit a dependence on the tiling used. Hence, the expansion at higher order should involve a tiling which is consistent with the physics of grain boundary contacts. All this is clearly beyond our reach today.

It should be noted that in the proof of Theorem 1 (and implicitly in that Theorem 5), the important feature about periodicity is translation invariance. Indeed, the crucial point, appart from the Taylor expansion, is the fact that for any lattice , we have:

∀i ∈ , j ∈ \{i} W 0 (j -i) = j ∈ \{0} W 0 (j ),
as long as W 0 decays fast enough at infinity. We cannot break this property without breaking translation invariance. Conversely, if a sort of translation invariance holds, for instance in the case of almost periodic systems (see [START_REF] Bohr | Almost periodic functions[END_REF][START_REF]Quasicrystals and geometry[END_REF]), it is possible to adapt our argument. Indeed, in this case, we would find the same kind of result:

E(u) = 1 | | E W 0 (∇u(x) )dx,
where E W 0 (∇u(x) ) denotes the two-body energy of the almost periodic set ∇u(x) , that is,

E W 0 ( ) = lim R→∞ 1 #(B R ∩ ∇u(x) ) i | =j ∈B R ∩∇u(x) W 0 (i -j),
which exists because , hence ∇u(x) , is almost periodic. Note that this quantity does not depend on the center of the ball B R , even in the case when this center depends on R. This property of existence and uniqueness of an average energy seems to be the crucial one in order to use the method of Theorem 1. Of course, even some sets which are not almost periodic enjoy this property (think for instance of the set Z 3 \ {0}), and the question of the characterization of such sets seems to be open.

Convexity and related properties of the homogenized energy

We make here a few remarks about the homegenized energies we have obtained in Section 2 and Section 3. These are only basic remarks, and we hope to come back to them in a more general setting in the near future. Starting with the zeroorder terms (7), ( 8), ( 36) and (37), we then study higher order terms obtained in Theorem 3.

Zero-order term

We start with the zero-order terms ( 7), ( 8), [START_REF] Pagano | A simple model for phase transitions: from the discrete to the continuum problem[END_REF] and [START_REF] Paroni | Homogenization of polycrystalline aggregates[END_REF], first pointing out that each of them exhibits the invariance (2), and in particular,

∀M ∈ M 3×3 + , ∀Q ∈ GL 3 (Z), E(MQ) = E(M). ( 64
)
Here M 3×3 + denotes the set of 3-by-3 matrices having positive determinant, and GL 3 (Z) the set of matrices having integer entries, positive determinant, and whose inverses have integer entries (this set is also equal to the set of matrices having integer entries and determinant equal to ±1). Equation (64) simply expresses the fact that the lattice defining the microscopic structure of the solid is invariant under a change of basis. Note that this remains valid in the high-density limit because

| det Q| = 1.
A first point is that this invariance prevents any growth at infinity of the energy, therefore excluding the use of standard tools of the calculus of variations. This, however, is not such a big limitation compared to some remarked on in [START_REF] Chipot | Equilibrium configurations of crystals[END_REF] and [START_REF] Fonseca | Variational methods for elastic crystals[END_REF], where it is stated that the invariance (64) prevents any quasiconvexity property. We provide here a simple proof of this fact in the case of a radially symmetric two-body interaction, referring to [START_REF] Fonseca | The lower quasiconvex envelope of the stored energy function for an elastic crystal[END_REF] and [START_REF] Fonseca | Variational methods for elastic crystals[END_REF] for a more general result. Proposition 7. Let W 0 : R + → R ∪ {+∞} satisfy the following:

(a) W 0 is of class C 1 on [R, +∞) for any R > 0; (b) ∃a > 0 / ∀t ∈ [1, +∞), |W 0 (t)| C t 3+a and |W 0 (t)| C t 4+a
, for some constant C 0; and (c) there exists some t 0 > 0 such that W 0 (t 0 ) | = 0, and W 0 is monotone on [t 0 , +∞).

Define the energy

E(M) = j ∈ \{0} W 0 (|Mj |)
for some fixed lattice , and for any matrix M ∈ M 3×3 + . Then E is not rank-one convex, thus not quasiconvex.

Let us point out that the conditions imposed on W 0 are fairly general, and include in particular almost all two-body potentials currently used in solid-state physics.

Proof. We assume for the sake of simplicity that = Z 3 . For any z ∈ R, λ, µ 0, we define

A(λ, µ, z) =   λ z 0 0 µ 0 0 0 µ   ,
and define the function ϕ(λ, µ, z) = E(A(λ, µ, z)). Assuming for the sake of contradiction that E is rank-one convex, ϕ is convex with respect to z. Now, we can easily show, using (64), that ϕ(λ, µ, λn) = ϕ(λ, µ, 0) for any integer n. Hence, ϕ must be independent of z. Differentiating it with respect to z, we thus find:

0 = j | =0 W 0 (λj 1 + zj 2 ) 2 + µ 2 j 2 2 + µ 2 j 2 3 (λj 1 + zj 2 )j 2 (λj 1 + zj 2 ) 2 + µ 2 j 2 2 + µ 2 j 2 3 .
We now let λ go to infinity, getting:

0 = j ∈Z 2 \{0} W 0 (µ 2 + z 2 )j 2 1 + µ 2 j 2 2 zj 2 1 (µ 2 + z 2 )j 2 1 + µ 2 j 2 2 .
Dividing this equality by z and then taking z to zero, implies that 0

= j ∈Z 2 \{0} W 0 (µ|j |) j 2 2 |j | .
Now, choosing µ = t 0 , we see that all terms of this sum have the same sign, while some of them are equal to W 0 (t 0 ) | = 0, which is contradictory.

However, in the high-density case, that is,

E(u) = 1 | | dx | det(∇u(x))| p , ( 65
)
with p = 1 or p = 2 3 , the function t → 1 t p being convex in either case with respect to t on R + , the energy E is polyconvex [START_REF] Ball | Convexity conditions and existence theorems in nonlinear elasticity[END_REF], hence quasiconvex. The energy being highly degenerate, it is necessary to add confining terms in order to have some equilibrium state. This kind of problem is dealt with in [START_REF] Podio-Guidugli | Cavitation and phase transition of hyperelastic fluids[END_REF].

Note that in the proof of Proposition 7, we have essentially used the fact that a convex function cannot be periodic, unless it is a constant. It is a direct consequence of [START_REF] Fonseca | The lower quasiconvex envelope of the stored energy function for an elastic crystal[END_REF] that a rank-one convex function satisfying (64) must be of the form (65): Theorem 8 (Fonseca, [START_REF] Fonseca | The lower quasiconvex envelope of the stored energy function for an elastic crystal[END_REF]). Let E be a function defined on M 3×3 + , the set of threeby-three matrices having positive determinant.

(i) Assume that E is bounded below and satisfies (64). Then, E is rank-one convex if and only if there exists a convex function g : (0, +∞) → R such that

E(F ) = g(det F ) ∀F ∈ M 3×3 + .
(ii) If E satisfies lim det F →0 + E(F ) = +∞ and (64), then its lower quasiconvex envelope QE is equal to its rank-one-convex envelope, and there exists a convex function g : (0, +∞) → R such that

QE(F ) = g(det F ) ∀F ∈ M 3×3 + .
Therefore, the application of this theorem to the TFW case allows us to conclude that a result similar to that of Proposition 7 holds in this case. Indeed, should the corresponding energy be quasiconvex, it would only depend on the determinant of the gradient deformation, according to (i) above. Hence, using

M =    λ 0 0 0 1 √ λ 0 0 0 1 √ λ   
as a gradient deformation, the corresponding energy should be independent of λ. This is in contradiction with the fact that as λ goes to zero, the energy goes to infinity (because a "great amount" of nuclei get closer and closer in the process; see [START_REF] Blanc | Geometry optimization for crystals in Thomas-Fermi type theories of solids[END_REF], Proposition 5.2 for a rigourous proof).

The boundary term

We now turn to the term of order one in [START_REF] Chipot | Equilibrium configurations of crystals[END_REF] and [START_REF] Ciarlet | Mathematical theory of elasticity[END_REF]. Here again, we give only a simple example of the role this term may play.

Consider the case ε = δ, that is [START_REF] Chipot | Equilibrium configurations of crystals[END_REF], discarding all terms of order strictly higher than 1. Assume in addition that = Q is the unit cube, and that = Z 3 . Then, the elastic energy reads:

E(u) = 1 2 j | =0 W 0 (∇u(x)j )dx - ε 2 ∂ k 1 j •n(x) k W 0 (∇u(x)j )dσ (x),
the sums over j being restricted to j ∈ Z 3 . The consequence of the presence of this first-order term is the breaking of property (64). Indeed, let us consider the example

u(x) =   x 1 + px 2 x 2 x 3   ,
where p is an integer. Then, we have Observe that the last two sums do not depend on p. We are now going to assume that the potential W 0 is radially symmetric, and satisfies

E(u) = 1
W 0 (x) = W 0 (|x|) < 0 ∀|x| > 1.
Hence, we have

E(u) = A 0 - ε 2 k 1 |j 1 | k W 0 (j 1 + pj 2 ) 2 + j 2 2 + j 2 3 A 0 -ε p k=1 j 3 ∈Z W 0 (p -p) 2 + 1 + j 2 3 A 0 -εp m∈Z W 0 1 + m 2 = A 0 + εBp,
where A 0 does not depend on p, and B > 0 is independent of ε and p. Therefore, as p goes to infinity, the energy grows like +εp, going to infinity.

The second-order term

We now study the second-order term. Let us start with the one-dimensional case (where ε = δ, that is, (i) of Theorem 3), in which the energy reads, if we neglect boundary terms and assume that = (0, 1):

E(u) = 1 2 1 0 j | =0 W 0 (u (x)j )dx - ε 2 24 1 0 j | =0 j 4 W 0 (u (x)j )(u (x)) 2 dx, ( 66
)
where the sums are over Z \ {0}, since we assume for the sake of simplicity that the lattice is equal to Z. Hence, we have an energy of the form

E(u) = 1 0 E 0 (u (x))dx + ε 2 h(u (x))(u (x)) 2 dx,
with h(y) = -1 24 j | =0 j 4 W 0 (yj ), and E 0 is the standard zero-order energy. Thus, if W 0 satisfies the inequality ∀y ∈ R, j | =0 j 4 W 0 (yj ) < 0, then the energy (66) exhibits a convexification term of order 2. The influence of this term on the energy has been studied in detail in [START_REF] Bardenhagen | On higher order gradient continuum theories in 1-D nonlinear elasticity. Derivation from and comparison to the corresponding discrete models[END_REF] for the case of the pure displacement problem, that is, I = inf E(u), u(0) = 0, u(1) = 1 + , for some ∈ R, corresponding to imposing a displacement of length at the right end of the solid, the other end remaining still. It is shown in [START_REF] Bardenhagen | On higher order gradient continuum theories in 1-D nonlinear elasticity. Derivation from and comparison to the corresponding discrete models[END_REF] that when ε = 0, I exhibits discontinuous critical points, and has no absolute minimizer in a classical sense. When the second-order term is added, and if h > 0, then bifurcation phenomena (with only smooth critical points) occur.

Let us now point out that there exist potentials for which the corresponding quantity h is indeed positive. This is the case for the Morse potential (well suited for a wide range of materials [START_REF] Stave | The structure of Ni N and P d N clusters[END_REF][START_REF] Wales | What can calculations employing empirical potentials teach us about bare transition metal clusters?[END_REF]). Indeed, we have the following lemma: Lemma 9. Consider the Morse potential, that is, W 0 (x) = e -2(x-r 0 ) -2e -(x-r 0 ) , and assume that the charactestic length r 0 satisfies r 0 < 5 log(5) -8 log(2) ≈ 2.502. Then it follows that ∀y > 0, j | =0 j 4 W 0 (yj ) < 0.

(67)

Proof. We have W 0 (x) = 4e -2(x-r 0 ) -2e -(x-r 0 ) . Hence, since we have j 1 j 4 e -αj = e -α 1 + 11e -α + 11e -2α + e -3α

(1e -α ) 5 , j | =0 j 4 W 0 (yj ) = 4e r 0 j 1 e r 0 2j 4 e -2yjj 4 e -yj = 4e r 0 2e r 0 e -2y 1 + 11e -2y + 11e -4y + e -6y

(1e -2y ) 5

e -y 1 + 11e -y + 11e -2y + e -3y

(1e -y ) 5 = 4e r 0 e -y (1e -y ) 5 × 2e r 0 e -y 1 + 11e -2y + 11e -4y + e -6y

(1 + e -y ) 5

-1 -11e -y -11e -2ye -3y .

Setting t = e -y and t 0 = e -r 0 , our aim is to show that ϕ(t) := t (1 + 11t 2 + 11t 4 + t 6 ) (1 + t) 5 (1 + 11t + 11t 2 + t 3 ) < t 0 ∀t ∈ (0, 1)

Since t ∈ (0, 1), t 2 < t, so that ϕ(t) < t (1+t) 5 . Studying the variations of t → t (1+t) 5 on [0, 1], we find that it has a maximum at 1 4 , where it is equal to 4 4 5 5 . Hence, if t 0 > 4 4 5 5 , condition (67) is satisfied. This inequality is equivalent to e -r 0 > 4 4 5 5 , that is, r 0 < 5 ln 5 + 8 ln 2.

On the other hand, it can be shown that the Lennard-Jones potential W 0 (x) = 1 x 12 -1 x 6 cannot satisfy (67). A similar analysis could be carried out in the high-density case, that is, (ii) of Theorem 3, the second-order term enjoying the same kind of property.

A similar study should be possible in the three-dimensional case, but calculations are a lot more involved, and it is not clear whether an ellipticity property can be derived in this case. However, the high-density case is more tractable: looking at the second-order term of ( 16), it is possible to change variables in the integral with respect to y, and obtain: This expression is certainly easier to study. For instance, if we only look for radially symmetric deformations, then assuming that W 0 is radially symmetric, the total energy is bounded from below by a norm of the second derivative of v, hence by a norm of the second derivative of u.

E 2 (u) : = - ε 2 24| | R 3
The question remains: In the general three-dimensional setting, and both for ( 15) and ( 16), is it possible to find some situations where the second-order bulk term exhibits ellipticity properties, as was assumed for instance in [START_REF] Bardenhagen | Derivation of higher order gradient continuum theories in 2,3-D non-linear elasticity for periodic lattice models[END_REF]?

Fig. 1 .

 1 Fig. 1. The reference configuration and the deformed one

Theorem 1 .

 1 Consider a periodic lattice such that its periodic cell Q( ) is of unit volume: |Q( )| = 1. Let W 0 be a function defined on R d \ {0}, which is Lipschitz on the exterior of any ball B R , with R > 0, and such that there exists C 0 and a > 0 such that |W 0 (x)| C |x| d+a . Let be a bounded Lipschitz open subset of R d , and let u be a C ∞ diffeomorphism defined on , with values in R d . Consider E ε,δ the energy defined by

) ε γ Fig. 2 .

 γ2 Fig. 2. The limiting process allowing the derivation of (62), with K = 2

D 2 W 3 D 2 1 D

 2321 0 (∇u(x)y)(D 2 u(x)(y, y), D 2 u(x)(y, y))dydx = -ε 2 24| | R W 0 (z) M(x, z), M(x, z) dz | det(∇u(x))| dx = -ε 2 24| | u( ) R 3 D 2 W 0 (z) N(ξ, z), N(ξ, z) dzdξ, where v = u -1 , M(x, z) = D 2 u(x)(∇u(x) -1 z, ∇u(x) -1 z), N(ξ, z) = ∇v(ξ ) det ∇v(ξ ) -2 v(ξ )(z, z).