N Rougerie

S Serfaty

Higher Dimensional Coulomb Gases and Renormalized Energy Functionals

We consider a classical system of n charged particles in an external confining potential, in any dimension d ≥ 2. The particles interact via pairwise repulsive Coulomb forces and the coupling parameter is of order n -1 (mean-field scaling). By a suitable splitting of the Hamiltonian, we extract the next to leading order term in the ground state energy, beyond the mean-field limit. We show that this next order term, which characterizes the fluctuations of the system, is governed by a new "renormalized energy" functional providing a way to compute the total Coulomb energy of a jellium (i.e. an infinite set of point charges screened by a uniform neutralizing background), in any dimension. The renormalization that cuts out the infinite part of the energy is achieved by smearing out the point charges at a small scale, as in Onsager's lemma. We obtain consequences for the statistical mechanics of the Coulomb gas: next to leading order asymptotic expansion of the free energy or partition function, characterizations of the Gibbs measures, estimates on the local charge fluctuations and factorization estimates for reduced densities. This extends results of Sandier and Serfaty to dimension higher than two by an alternative approach.

Introduction

We study the equilibrium properties of a classical Coulomb gas or "one-component plasma": a system of n classical charged particles living in the full space of dimension d ≥ 2, interacting via Coulomb forces and confined by an external electrostatic potential V . We will be interested in the mean-field regime where the number n of particles is large and the pair-interaction strength (coupling parameter) scales as the inverse of n. We study the ground states of the system as well its statistical mechanics when temperature is added. Denoting x 1 , . . . , x n the positions of the particles, the total energy at rest of such a system is given by the Hamiltonian

H n (x 1 , . . . , x n ) = i =j w(x i -x j ) + n n i=1 V (x i ) (1.1) where    w(x) = 1 |x| d-2 if d ≥ 3 w(x) = -log |x| if d = 2 (1.2)
is a multiple of the Coulomb potential in dimensions d ≥ 2, i.e. we have

-∆w = c d δ 0 (1.3) with c 2 = 2π, c d = (d -2)|S d-1 | when d ≥ 3 (1.4)
and δ 0 is the Dirac mass at the origin. The one-body potential V : R d → R is a continuous function, growing at infinity (confining potential). More precisely, we assume

     lim |x|→∞ V (x) = +∞ if d ≥ 3 lim |x|→∞ V (x) 2 -log |x| = +∞ if d = 2.
(1.5)

Note the factor n in front of the one-body term (second term) in (1.1) that puts us in a meanfield scaling where the one-body energy and the two-body energy (first term) are of the same order of magnitude. This choice is equivalent to demanding that the pair-interaction strength be of order n -1 . One can always reduce to this situation in the particular case where the trapping potential V has some homogeneity, which is particularly important in applications.

We will not treat at all the case of one-dimensional Coulomb gases (where the interaction kernel w is |x|), since this case has been shown to be essentially completely solvable a long time ago [START_REF] Aizenman | Structure of Gibbs States of one Dimensional Coulomb Systems[END_REF][START_REF] Lenard | Exact statistical mechanics of a one-dimensional system with Coulomb forces[END_REF][START_REF] Lenard | Exact statistical mechanics of a one-dimensional system with Coulomb forces III: Statistics of the electric field[END_REF][START_REF] Brascamp | Some Inequalities for Gaussian Measures and the long-range order of the one-dimensional plasma[END_REF][START_REF] Kunz | The One-Dimensional Classical Electron Gas[END_REF]. Classical Coulomb systems are fundamental systems of statistical mechanics, since they can be seen as a toy model for matter, containing the truly long-range nature of electrostatic interactions. Studies in this direction include [SM, LO, JLM, PS], see [Ser] for a review. Another motivation is that, as was pointed out by Wigner [Wi2] and exploited by Dyson [Dys], two-dimensional Coulomb systems are directly related to Gaussian random matrices, more precisely the Ginibre ensemble, and such random matrix models have also received much attention for their own sake. A similar connection exists between "log-gases" in dimension 1 and the GUE and GOE ensembles of random matrices, as well as more indirectly to orthogonal polynomial ensembles. For more details on these aspects we refer to [For], and for an introduction to the random matrix aspect to the texts [AGZ, Meh, Dei]. A recent trend in random matrix theory is the study of universality with respect to the entries' statistics, i.e. the fact that results for Gaussian entries carry over to the general case, see e.g. [TV, ESY].

We are interested in equilibrium properties of the system in the regime n → ∞, that is on the large particle number asymptotics of the ground state and the Gibbs state at given temperature. In the former case we consider configurations (x 1 , . . . , x n ) that minimize the total energy (1.1). We will denote E n := min

R dn H n (1.6)
the ground state energy. It is well-known (we give references below) that to leading order

E n = n 2 E[µ 0 ](1 + o(1)) (1.7)
in the limit n → ∞ where

E[µ] = R d ×R d w(x -y) dµ(x) dµ(y) + R d V (x) dµ(x) (1.8)
is the mean-field energy functional defined for Radon measures µ, and µ 0 (the equilibrium measure) is the minimizer of E amongst probability measures on R d . In this paper we quantify precisely the validity of (1.7) and characterize the next to leading order correction. We also study the consequences of these asymptotics on minimizing and thermal configurations. By the latter we mean the Gibbs state at inverse temperature β, i.e. the probability law P n,β (x 1 , . . . , x n ) = 1 Z β n e -β 2 Hn(x 1 ,...,xn) dx 1 . . . dx n (1.9)

where Z β n is a normalization constant, and we are again interested in obtaining next order expansions of the partition function Z β n as well as consequences for the distributions of the points according to the temperature. This program has been carried out in [START_REF] Sandier | 2D Coulomb gases and the renormalized energy[END_REF] in dimension d = 2 and here we extend it to arbitrary higher dimension -in particular the more physical case d = 3, and provide at the same time a simpler approach to recover (most of) the results of dimension 2. In [START_REF] Sandier | 2D Coulomb gases and the renormalized energy[END_REF] it was shown that the next order corrections are related to a "renormalized energy" denoted W -so named in reference to the procedure used in its definition and related functionals used in Ginzburg-Landau theory [START_REF] Bethuel | Progress in Nonlinear Differential Equations and their Applications 13[END_REF][START_REF] Sandier | Vortices in the Magnetic Ginzburg-Landau Model[END_REF], but the derivation of this object was restricted to dimension 2 due to an obstruction that we still do not know how to overcome (more precisely the derivation relies on a "ball construction method" which crucially uses the conformal invariance of the Coulomb kernel in two dimensions). Here we again connect the problem to a slightly different "renormalized energy," this time denoted W (it is the same in good cases, but different in general) and the approach to its definition and derivation are not at all restricted by the dimension. They rely on smearing out point charges and Onsager's lemma [Ons], a celebrated tool that has been much used in the proof of the stability of matter (see [LO] and [START_REF] Lieb | The stability of matter in quantum mechanics[END_REF]Chapter 6]).

We choose to use an electrostatic/statistical mechanics vocabulary that is more fit to general dimensions but the reader should keep in mind the various applications of the Coulomb gas, especially in two dimensions: Fekete points in polynomial interpolation [SaTo], Gaussian random matrices which correspond to β = 2 [For], vortex systems in classical and quantum fluids [START_REF] Caglioti | A Special Class of Stationary Flows for Two-Dimensional Euler Equations: A Statistical Mechanics Description[END_REF][START_REF] Sandier | From the Ginzburg-Landau Model to Vortex Lattice Problems[END_REF][START_REF] Correggi | Energy and Vorticity in Fast Rotating Bose-Einstein Condensates[END_REF][START_REF] Correggi | The transition to a Giant Vortex Phase in a Fast Rotating Bose-Einstein Condensate[END_REF], fractional quantum Hall physics [START_REF] Girvin | Introduction to the fractional quantum Hall effect[END_REF][START_REF] Rougerie | Quantum Hall states of bosons in rotating anharmonic traps[END_REF][START_REF] Rougerie | Quantum Hall phases and plasma analogy in rotating trapped Bose gases[END_REF] ... The easiest way to think of the limit (1.7) is as a continuum limit: the one-body potential nV confines the large number n of particles in a bounded region of space, so that the mean distance between points goes to zero and the empirical measure

µ n := 1 n n i=1 δ x i (1.10)
of a minimizing configuration converges to a density ρ. The functional (1.8) is nothing but the continuum energy corresponding to (1.1): the first term is the classical Coulomb interaction energy of the charge distribution µ and the second the potential energy in the potential V .

If the interaction potential w was regular at the origin we could write

H n (x 1 , . . . , x n ) = n 2 R d V (x) dµ n (x) + R d ×R d w(x -y) dµ n (x) dµ n (y) -nw(0) = n 2 E[µ n ] 1 + O(n -1 )
and (1.7) would easily follow from a simple compactness argument. In the case where w has a singularity at the origin, a regularization procedure is needed but this mean-field limit result still holds true, meaning that for minimizers of H n , the empirical measure µ n converges to µ 0 , the minimizer of (1.8). This is standard and can be found in a variety of sources: e.g. [START_REF] Saff | Logarithmic Potentials with External Fields[END_REF]Chap. 1] for the Coulomb kernel in dimension 2, [CGZ] for a more general setting including possibly non-coulombian kernels, or [Ser] for a simple general treatment. This leading order result is often complemented by a much stronger large deviations principle in the case with temperature: a large deviations principle with speed n 2 and good rate function βE holds; see [PH, BZ, Har] for the two-dimensional Coulomb case (with β = 2), which can be adapted to any finite temperature and any dimension [CGZ, Ser]. This is also of interest in the more elaborate settings of complex manifolds, cf. e.g. [Ber, BBN] and references therein.

Another way to think of the mean-field limit, less immediate in the present context but more suited for generalizations in statistical and quantum mechanics, is as follows. In reality, particles are indistinguishable, and the configuration of the system should thus be described by a probability measure µ(x) = µ(x 1 , . . . , x n ), which is symmetric under particle exchange: µ(x 1 , . . . , x n ) = µ(x σ(1) , . . . , x σ(n) ) for any permutation σ.

(1.11)

An optimal (i.e. most likely) configuration µ n is found by minimizing the n-body energy functional

I n [µ] := R dn
H n (x)µ(dx) (1.12) amongst symmetric probability measures µ ∈ P s (R dn ) (probability measures on R dn satisfying (1.11)). It is immediate to see that µ n must be a convex superposition of measures of the form δ (x 1 ,...,xn) with (x 1 , . . . , x n ) minimizing H n (in other words it has to be a symmetrization of some δ x for a minimizing configuration x). The infimum of the functional (1.12) of course coincides with inf

µ∈Ps(R dn ) R dn H n (x)µ(dx) = E n
and a way to understand the asymptotic formula (1.7) is to think of the minimizing µ n as being almost factorized

µ n (x 1 , . . . , x n ) ≈ ρ ⊗n (x 1 , . . . , x n ) = n j=1 ρ(x i ) (1.13)
with a regular probability measure ρ ∈ P(R d ). Plugging this ansatz into (1.12) we indeed obtain

I n [ρ ⊗n ] = n 2 E[ρ] 1 + O(n -1 )
and the optimal choice is ρ = µ 0 . The mean-field limit can thus also be understood as one where correlations amongst the particles of the system vanish in the limit n → ∞, which is the meaning of the factorized ansatz.

In this paper we shall pursue both the "uncorellated limit" and the "continuum limit" points of view beyond leading order considerations. That is, we shall quantify to which precision (and in which sense) the empirical measure (2.35) of a minimizing configuration can be approximated by µ 0 and the n-body ground state factorizes in the form (1.13) with ρ = µ 0 . Previous results related to the "uncorellated limit" point of view may be found in [START_REF] Caglioti | A Special Class of Stationary Flows for Two-Dimensional Euler Equations: A Statistical Mechanics Description[END_REF][START_REF] Kiessling | On the Equilibrium Statistical Mechanics of Isothermal Classical Self-Gravitating Matter[END_REF][START_REF] Kiessling | Statistical mechanics of classical particles with logarithmic interactions[END_REF][START_REF] Kiessling | A Note on the Eigenvalue Density of Random Matrices[END_REF][START_REF] Messer | Statistical mechanics of the isothermal Lane-Emden equation[END_REF]. Another way of viewing this is that we are looking at characterizing the "fluctuations" of the distribution of points around its limit measure, in other words the behavior of n(µ n -µ 0 ) where µ n is the empirical measure (1.10). In the probability literature, such questions are now understood in dimension 2 for the particular determinantal case β = 2 [START_REF] Ameur | Fluctuations of eigenvalues of random normal matrices[END_REF][START_REF] Ameur | Random normal matrices and Ward identities[END_REF], and in dimension 1 with the logarithmic interaction [START_REF] Valkó | Continuum limits of random matrices and the Brownian carousel[END_REF][START_REF] Bourgade | Universality of general β-ensembles[END_REF][START_REF] Bourgade | Bulk Universality of General β-ensembles with non-convex potential[END_REF]. Our results are less precise (we do not exhibit exact local statistics of spacings), but they are valid for any β, any V , and any dimension d ≥ 2 (we could also treat the log gas in dimension 1, borrowing ideas from [START_REF] Sandier | 1D Log Gases and the Renormalized Energy: Crystallization at Vanishing Temperature[END_REF] to complete those we use here). We are not aware of any previous results giving any information on such fluctuations for ground states or thermal states in dimension ≥ 3.

Most of our results follow from almost exact splitting formulae for the Hamiltonian (1.1) that reveal the corrections beyond leading order in (1.7), in the spirit of [START_REF] Sandier | 2D Coulomb gases and the renormalized energy[END_REF]. Let us first explain what physics governs these corrections. As already mentioned, points minimizing H n tend to be densely packed in a bounded region of space (the support of µ 0 , that we shall denote Σ) in the limit n → ∞. Their distribution (i.e. the empirical measure) has to follow µ 0 on the macroscopic scale but this requirement still leaves a lot of freedom on the configuration at the microscopic scale, that is on lengthscales of order n -1/d (the mean inter-particle distance). A natural idea is thus to blow-up at scale n -1/d in order to consider configurations where points are typically separated by distances of order unity, and investigate which microscopic configuration is favored. On such length scales, the equilibrium measure µ 0 varies slowly so the points will want to follow a constant density given locally by the value of µ 0 . Since the problem is electrostatic in nature it is intuitive that the correct way to measure the distance between the configuration of points and the local value of the equilibrium measure should use the Coulomb energy. This leads to the idea that the local energy around a blow-up origin should be the electrostatic energy of what is often called a jellium in physics: an infinite collection of interacting particles in a constant neutralizing background of opposite charge, a model originally introduced in [START_REF] Wigner | On the interaction of electrons in metals[END_REF]. At this microscopic scale the pair-interactions will no longer be of mean-field type, their strength will be of order 1. The splitting formula will allow to separate exactly the Coulomb energy of this jellium as the next to leading order term, except that what will come out is more precisely some average of all the energies of the jellium configurations obtained after blow-up around all possible origins.

Of course it is a delicate matter to define the energy of the general infinite jellium in a mathematically rigorous way: one has to take into account the pair-interaction energy of infinitely many charges, without assuming any local charge neutrality, and the overall energy may be finite only via screening effects between the charges and the neutralizing background that are difficult to quantify. This has been done for the first time in 2D in [START_REF] Sandier | From the Ginzburg-Landau Model to Vortex Lattice Problems[END_REF][START_REF] Sandier | 2D Coulomb gases and the renormalized energy[END_REF], the energy functional for the jellium being the renormalized energy W alluded to above. As already mentioned, one of the main contributions of the present work is to present an alternate definition W that generalizes better to higher dimensions. The precise definition will be given later, but we can already state our asymptotic formula for the ground state energy (minimum of H n ), where α d denotes the minimum of W for a jellium of density 1 in dimension d:

E n =        n 2 E[µ 0 ] + n 2-2/d c d α d µ 2-2/d 0 (x)dx + o n 2-2/d if d ≥ 3 n 2 E[µ 0 ] - n 2 log n + n α 2 2π - 1 2 µ 0 (x) log µ 0 (x) dx + o (n) if d = 2.
(1.14)

This formula encodes the double scale nature of the charge distribution: the first term is the familiar mean-field energy and is due to the points following the macroscopic distribution µ 0 -we assume that the probability µ 0 has a density that we also denote µ 0 by abuse of notation.

The next order correction which happens to lie at the order n 2-2/d governs the configurations at the microscopic scale and the crystallization, by selecting configurations which minimize W (on average, with respect to the blow-up centers).

The factors involving µ 0 in the correction come from scaling, and from the fact that the points locally see a neutralizing background whose charge density is given by the value of µ 0 . The scaling properties of the renormalized energy imply that the minimal energy of a jellium with neutralizing density µ 0 (x) is α d µ 2-2/d 0 (x) (respectively µ 0 (x)α 2 -µ 0 (x) log µ 0 (x) when d = 2). Integrating this energy density on the support of µ 0 leads to the formula.

The interpretation of the correction is thus that around (almost) any point x in the support of µ 0 there are approximately µ 0 (x) points per unit volume, distributed so as to minimize a jellium energy with background density µ 0 (x). Due to the properties of the jellium, this implies that, up to a µ 0 (x)-dependent rescaling, the local distribution of particles are the same around any blow-up origin x in the support of µ 0 . This can be interpreted as a result of universality with respect to the potential V in (1.1), in connection with recent works on the 1D log gas [START_REF] Bourgade | Universality of general β-ensembles[END_REF][START_REF] Bourgade | Bulk Universality of General β-ensembles with non-convex potential[END_REF][START_REF] Sandier | 1D Log Gases and the Renormalized Energy: Crystallization at Vanishing Temperature[END_REF].

We remark that even the weaker result that the correction in (1.14) is exactly of order n2-2/d (with an extra term of order n log n in 2D due to the scaling properties of the log) did not seem to have been previously noticed (except in [START_REF] Sandier | 2D Coulomb gases and the renormalized energy[END_REF] where the formula (1.14) is derived in the case d = 2).

The next natural question is of course that of the nature of the minimizers of the renormalized energy W. It is widely believed that the minimizing configuration (at temperature 0) consists of points distributed on a regular lattice (Wigner crystal). A proof of this is out of our present reach: crystallization problems have up to now been solved only for specific short range interaction potentials (see [The, BPT, HR, Sut, Rad] and references therein) that do not cover Coulomb forces, or 1D systems [START_REF] Brascamp | Some Inequalities for Gaussian Measures and the long-range order of the one-dimensional plasma[END_REF][START_REF] Kunz | The One-Dimensional Classical Electron Gas[END_REF][START_REF] Sandier | 1D Log Gases and the Renormalized Energy: Crystallization at Vanishing Temperature[END_REF]. In [START_REF] Sandier | From the Ginzburg-Landau Model to Vortex Lattice Problems[END_REF] it was shown however that in dimension 2, if the minimizer is a lattice, then it has to be the triangular one, called the Abrikosov lattice in the context of superconductivity. In dimension 3 and higher, the question is wide open, not only to prove that minimizers are crystalline, but even to identify the optimal lattice configurations. The FCC (face centered cubic) lattice, and maybe also the BCC (body centered cubic) lattice, seem like natural candidates, and the optimisation of W is related to the computation of what physicists and chemists call their Madelung constants.

This question is in fact of number-theoretic nature: in [START_REF] Sandier | 2D Coulomb gases and the renormalized energy[END_REF] it is shown that the question in dimension 2 reduces to minimizing the Epstein Zeta function p∈Λ |p| -s with s > 2 among lattices Λ, a question which was in turn already solved in dimension 2 in the 60's (see [START_REF] Cassels | On a problem of Rankin about the Epstein zeta-function[END_REF][START_REF] Rankin | A minimum problem for the Epstein zeta function[END_REF][START_REF] Ennola | A remark about the Epstein zeta function[END_REF][START_REF] Ennola | On a problem about the Epstein zeta-function[END_REF][START_REF] Diananda | Notes on two lemmas concerning the Epstein zeta-function[END_REF] and also [Mon] and references therein), but the same question is open in dimension d ≥ 3 -except for d = 8 and d = 24 -and it is only conjectured that the FCC is a local minimizer, see [SaSt] and references therein. The connection with the minimization of W among lattices and that of the Epstein Zeta functions is not even rigorously clear in dimension d ≥ 3. For more details, we refer to Section 3.1 where we examine and compute W in the class of periodic configurations and discuss this question.

Modulo the conjecture that the minimizer of the renormalized energy is a perfect lattice, (1.14) (and the consequences for minimizers that we state below) can be interpreted as a crystallization result for the Coulomb gas at zero temperature.

Concerning the distribution of charges at positive temperature, one should expect a transition from a crystal at low temperature to a liquid at large enough temperature. While we do not have any conclusive proof of this fact, our results on the Gibbs state (1.9) strongly suggest that the transition should happen in the regime β ∝ n 2/d-1 . In particular we prove that (see Theorem 3 below):

1. If β ≫ n 2/d-1 (low temperature regime) in the limit n → ∞, then the free-energy (linked to the partition function Z β n ) is to leading order given by the mean-field energy, with the correction expressed in terms of the renormalized energy as in (1.14). In other words, the free energy and the ground state energy agree up to corrections of smaller order than the contribution of the jellium energy. We interpret this as a first indication of crystallization in this regime. Note that for the 1D log gas in the corresponding regime, a full proof of crystallization has been provided in [START_REF] Sandier | 1D Log Gases and the Renormalized Energy: Crystallization at Vanishing Temperature[END_REF]. correction to the free energy is no longer given by the renormalized energy, but rather by an entropy term. We take this as a weak indication that the Gibbs state is no longer crystalline, but a more detailed analysis would be required.

Note the dependence on d of the critical order of magnitude: it is of order 1 only in 2D. Interestingly, in the main applications where the Gibbs measure of the 2D Coulomb gas arises (Gaussian random matrices and quantum Hall phases), the inverse temperature β is a number independent of n, i.e. one is exactly at the transition regime.

In the next section we proceed to state our results rigorously. Apart from what has already been mentioned above they include:

• Results on the ground state configurations that reveal the two-scale structure of the minimizers hinted at by (1.14), see Theorem 2.

• A large deviations -type result in the low temperature regime (Theorem 4), that shows that in the limit β ≫ n 2/d-1 the Gibbs measure charges only configurations whose renormalized energy converges to the minimum. When β ∝ n 2/d-1 only configurations with a certain upper bound on W are likely.

• Estimates on the charge fluctuations in the Gibbs measure, in the low temperature regime again (Theorem 5). These are derived by exploiting coercivity properties of the renormalized energy.

• Precise factorization estimates for the ground state and Gibbs measure, giving a rigorous meaning to (1.13), that we obtain as corollaries of our estimates on charge fluctuations, see Corollary 6.

These results, let us note, should be seen as coming in two groups. The first one reduces the crystallisation question for the trapped Coulomb gas to the same question for the jellium, in a suitable weak sense. Roughly speaking, we show that configurations for the trapped Coulomb gas look almost everywhere locally crystalline at small enough temperature, in an average sense that we make precise in Sections 2.2 and 2.3. Our second set of results aims at quantifying the deviations from mean-field theory in a stronger sense than that we use in Sections 2.2 and 2.3. Here our results (stated in Section 2.4), are not optimal, for example we can only show that particles are localized within a precision O(n -1/d+2 ), whereas crystallization happens at scale n -1/d . The estimates we obtain are however, as far as we know, the best of their kind in d ≥ 3. Concerning this second group of results we finally mention that, although we focus on the low temperature regime, our methods (or close variants) also yield explicit estimates on the Gibbs measure for any temperature. They can be used if the need for robust and explicit information on the behavior of the Gibbs measure arises, as e.g. in [START_REF] Rougerie | Quantum Hall states of bosons in rotating anharmonic traps[END_REF][START_REF] Rougerie | Quantum Hall phases and plasma analogy in rotating trapped Bose gases[END_REF] where we studied quantum Hall states with related techniques.

Acknowledgments: S.S. was supported by a EURYI award and would like to thank the hospitality of the Forschungsinstitut für Mathematik at the ETH Zürich, where part of this work was completed. N.R. thanks Gian Michele Graf and Martin Fraas for their hospitality at the Institute for Theoretical Physics of ETH Zürich, and acknowledges the financial support of the ANR (project Mathostaq ANR-13-JS01-0005-01).

Statement of main results

We let the mean-field energy functional be as in (1.8). The minimization of E among P(R d ), the space of probability measures on R d , is a standard problem in potential theory (see [Fro] or [SaTo] for d = 2). The uniqueness of a minimizer, called the equilibrium measure is obvious by strict convexity of E and its existence can be proven using the continuity of V and the assumption (1.5) -note that it only depends on the data of V and the dimension. We recall it is denoted by µ 0 , assumed to have a density also denoted µ 0 (x) by abuse of notation, and its support is denoted by Σ and usually called the droplet. We will need the following additional assumptions on µ 0 :

∂Σ is C 1 (2.1) µ 0 is C 1 on its support Σ and 0 < m ≤ µ 0 (x) ≤ m < ∞ there. (2.2)
It is easy to check that these are satisfied for example if V is quadratic, in which case µ 0 is the characteristic function of a ball. More generally, if V ∈ C 2 then µ 0 (x) = ∆V (x) on its support. The regularity assumption is purely technical: essentially it makes the construction needed for the upper bound easier. The upper and lower bounds to µ 0 on its support ensure that the "jellium energy with background density µ 0 (x)", which will give the local energy fluctuation around x ∈ Σ, makes sense . We note that µ 0 is also related to the solution of an obstacle problem (see the beginning of Section 3), and if V ∈ C 2 and ∆V > 0 then the droplet Σ coincides with the coincidence set of the obstacle problem, for which a regularity theory is known [Caf].

The renormalized energy

This section is devoted to the precise definition of the renormalized energy. It is defined via the electric field E generated by the full charge system: a (typically infinite) distribution of point charges in space in a constant neutralizing background. Note first that the classical Coulomb interaction of two charge distributions (bounded Radon measures) f and g,

D(f, g) := R d ×R d w(x -y) df (x) dg(y) (2.3) is linked to the (electrostatic) potentials h f = w * f , h g = w * g that they generate via the formula 1 D(f, g) = R d f h g = R d gh f = 1 c d R d ∇h f • ∇h g (2.4)
where we used the fact that by definition of w,

-∆h f = c d f, -∆h g = c d g.
The electric field generated by the distribution f is given by ∇h f , and its square norm thus gives a constant times the electrostatic energy density of the charge distribution f :

D(f, f ) = 1 c d R d |∇h f | 2 .
The electric field generated by a jellium is of the form described in the following definition.

1 Here we assume that all quantities are well-defined.

Definition 2.1 (Admissible electric fields).

Let m > 0. Let E be a vector field in R d . We say that E belongs to the class

A m if E = ∇h with -∆h = c d p∈Λ N p δ p -m in R d (2.5)
for some discrete set Λ ⊂ R d , and N p integers in N * .

Note that in [START_REF] Sandier | From the Ginzburg-Landau Model to Vortex Lattice Problems[END_REF] a different convention was used, the electric field being rotated by π/2 at each point to represent a superconducting current. This made the analogy with 2D Ginzburg-Landau theory more transparent but does not generalize easily to higher dimensions. In the present definition h corresponds to the electrostatic potential generated by the jellium and E to its electric field, while the constant m represents the mean number of particles per unit volume, or the density of the neutralizing background. An important difficulty is that the electrostatic energy D(δ p , δ p ) of a point charge density δ p , where δ p denotes the Dirac mass at p, is infinite, or in other words, that the electric field generated by point charges fails to be in L 2 loc . This is where the need for a "renormalization" of this infinite contribution comes from.

To remedy this, we replace point charges by smeared-out charges, as in Onsager's lemma: We pick some arbitrary fixed radial nonnegative function ρ, supported in B(0, 1) and with integral 1. For any point p and η > 0 we introduce the smeared charge

δ (η) p = 1 η d ρ x η * δ p . (2.6)
Even though the value of W will depend in general on the precise choice of ρ, the results in the paper will not depend on it, implying in particular that the value of the minimum of W will not either (thus we do not try to optimize over the possible choices). A simple example is to take ρ = 1 |B(0,1)| 1 B(0,1) , in which case

δ (η) p = 1 |B(0, η)| 1 B(0,η) .
We also define

κ d := c d D(δ (1) 0 , δ (1) 0 ) for d ≥ 3, κ 2 = c 2 , γ 2 = c 2 D(δ (1) 0 , δ (1)
0 ) for d = 2.

(2.7)

The numbers κ d , γ 2 , depend only on the choice of the function ρ and on the dimension. This nonsymmetric definition is due to the fact that the logarithm behaves differently from power functions under rescaling, and is made to ease the formulas below. Newton's theorem [START_REF] Lieb | Analysis[END_REF]Theorem 9.7] asserts that the Coulomb potentials generated by the smeared charge δ (η) p and the point charge δ p coincide outside of B(p, η). A consequence of this is that there exists a radial function f η solution to

-∆f η = c d δ (η) 0 -δ 0 in R d f η ≡ 0 in R d \B(0, η). (2.8)
and it is easy to define the field E η generated by a jellium with smeared charges starting from the field of the jellium with (singular) point charges, using f η :

Definition 2.2 (Smeared electric fields).

For any vector field E = ∇h satisfying

-div E = c d p∈Λ N p δ p -m (2.9)
in some subset U of R d , with Λ ⊂ U a discrete set of points, we let

E η := ∇h + p∈Λ N p ∇f η (x -p) h η = h + p∈Λ N p f η (x -p).
We have

-div E η = -∆h η = c d p∈Λ N p δ (η) p -m (2.10)
and denoting by Φ η the map E → E η , we note that Φ η realizes a bijection from the set of vector fields satisfying (2.9) and which are gradients, to those satisfying (2.10) and which are also gradients.

Note that the above definition in principle depends implicitly on the set U , whose choice will be clear from the context in the sequel (most of the time we will take U = R d ).

For any fixed η > 0 one may then define the electrostatic energy per unit volume of the infinite jelium with smeared charges as lim sup

R→∞ - K R |E η | 2 := lim sup R→∞ |K R | -1 K R |E η | 2 (2.11)
where E η is as in the above definition and K R denotes the cube [-R, R] d . Note that the quantity -K R |E η |2 may not have a limit (this does occur for somewhat pathological configurations). This motivates the use of the lim sup instead, which gives the strongest possible control on the configuration (see e.g. Lemma 3.1 below). In addition, for the configurations we obtain as asymptotic limits, the ergodic theorem which we use to exhibit these quantities will always ensure the existence of a true limit. This energy is now well-defined for η > 0 and blows up as η → 0, since it includes the self-energy of each smeared charge in the collection, absent in the original energy (i.e. in the Hamiltonian (1.1)). One should thus renormalize (2.11) by removing the self-energy of each smeared charge before taking the limit η → 0. We will see that the leading order energy of a smeared charge is κ d w(η), and this is the quantity that should be removed for each point. But in order for the charges to efficiently screen the neutralizing background, configurations will need to have the same charge density as the neutralizing background (i.e. m points per unit volume). We will prove in Lemma 3.1 that this holds. We are then led to the definition Definition 2.3 (The renormalized energy).

For any E ∈ A m , we define

W η (E) = lim sup R→∞ - K R |E η | 2 -m(κ d w(η) + γ 2 1 d=2 ) (2.12)
and the renormalized energy is given by 2

W(E) = lim inf η→0 W η (E) = lim inf η→0 lim sup R→∞ - K R |E η | 2 -m(κ d w(η) + γ 2 1 d=2 ) . It is easy to see that if E ∈ A m , then E ′ := m 1/d-1 E(m -1/d .) belongs to A 1 and W η (E) = m 2-2/d W ηm 1/d (E ′ ) W(E) = m 2-2/d W(E ′ ) if d ≥ 3 W η (E) = m W ηm 1/d (E ′ ) -κ 2 2 log m W(E) = m W(E ′ ) -κ 2 2 log m if d = 2,
(2.13) thus the same scaling formulae hold for inf Am W. One may thus reduce to the study of W(E) on A 1 , for which we have the following result:

Theorem 1 (Minimization of the renormalized energy).

The infimum

α d := inf E∈A 1 W(E) (2.14)
is achieved and is finite. Moreover, there exists a sequence (E n ) n∈N of periodic vector fields (with diverging period in the limit n → ∞) in A 1 such that

W(E n ) → α d as n → ∞. (2.15)
We should stress at this point that the definition of the jellium (renormalized) energy we use is essential for our approach to the study of equilibrium states of (1.1). In particular it is crucial that we are allowed to define the energy of an infinite system via a local density (square norm of the electric field). For a different definition of the jellium energy, the existence of the thermodynamic limit was previously proved in [LN], using ideas from [START_REF] Lieb | Existence of thermodynamics for real matter with Coulomb forces[END_REF][START_REF] Lieb | The Constitution of Matter: Existence of Thermodynamics for Systems Composed of Electrons and Nuclei[END_REF]. This approach does not transpose easily in our context however, and our proof that α d is finite follows a different route, using in particular an unpublished result of Lieb [Lie2]. The existence of a minimizer is in fact a consequence of our main results below (see the beginning of Section 6). A direct proof is also provided in Appendix A for convenience of the reader.

In [START_REF] Sandier | From the Ginzburg-Landau Model to Vortex Lattice Problems[END_REF][START_REF] Sandier | 2D Coulomb gases and the renormalized energy[END_REF] a different but related strategy was used for the definition of W (in dimension 2). Instead of smearing charges out, the "renormalization" was implemented by cutting-off the electric field E in a ball of radius η around each charge, as in [BBH]. This leads to the following definition, which we may present in arbitrary dimension (the normalizing of constants has been slightly modified in order to better fit with the current setting): Definition 2.4 (Alternative definition [START_REF] Sandier | From the Ginzburg-Landau Model to Vortex Lattice Problems[END_REF][START_REF] Sandier | 2D Coulomb gases and the renormalized energy[END_REF]). We let A m denote the subclass of A m for which all the points are simple (i.e. N p = 1 for all p ∈ Λ.) For each E ∈ A m , we define

W (E) = lim sup R→∞ 1 |K R | W (E, χ K R ) (2.16)
where, for any E satisfying a relation of the type (2.9), and any nonnegative continuous function χ, we denote

W (E, χ) = lim η→0 R d \∪ p∈Λ B(p,η) χ|E| 2 -c d w(η) p∈Λ χ(p),
(2.17)

and {χ K R } R>0 denotes a family of cutoff functions, equal to 1 in K W R-1 , vanishing outside K R and of universally bounded gradient.

In addition to the way the renormalization is performed, between the two definitions the order of the limits η → 0 and R → ∞ is reversed. It is important to notice that for a given discrete configuration, the minimal distance between points is bounded below on each compact set, hence the balls B(p, η) in K R appearing in (2.17) are disjoint as soon as η is small enough, for each fixed R. On the contrary, the smeared out charges δ

(η)
p in the definition of E η may overlap. In fact, we can prove (see Section 3.1) at least in dimension 2, that if the distances between the points is bounded below by some uniform constant (not depending on R) -we will call such points "well separated", in particular they must all be simple -, then the order of the limits η → 0 and R → ∞ can be reversed and W and W coincide (in addition the value of W then does not depend on the choice of ρ by Newton's theorem). An easy example is the case of a configuration of points which is periodic. But if the points are not well separated, then in general W and W do not coincide (with typically W ≤ W ). An easy counter example of this is the case of a configuration of well-separated points except one multiple point. Then computing the limit η → 0 in W immediately yields +∞. On the contrary, the effect of the multiple point gets completely dissolved when taking first the limit R → ∞ in the definition of W. At least an immediate consequence of (1.14) (or Theorem 2 below), by comparison with the result of [START_REF] Sandier | 2D Coulomb gases and the renormalized energy[END_REF], is that we know that min

A 1 W = min A 1 W in dimension 2.
One of the advantages of W is that it is more precise: in 2D it can be derived as the complete Γ-limit of H n at next order, while W is not (it is too "low"), see below. It also seems more amenable to the possibility of charges of opposite signs (for example it is derived in [START_REF] Sandier | From the Ginzburg-Landau Model to Vortex Lattice Problems[END_REF] as the limit of the vortex interactions in the Ginzburg-Landau energy where vortices can a priori be positive or negative). The main advantage of W is that it is much easier to bound it from below: W η is bounded below from its very definition, while proving that W is turned out to be much more delicate. In [START_REF] Sandier | From the Ginzburg-Landau Model to Vortex Lattice Problems[END_REF] it was proven that even though the energy density that defines W (E, χ) is unbounded below, it can be shown to be very close to one that is. This was a crucial point in the proofs, since most lower bound techniques (typically Fatou's lemma) require energy densities that are bounded below. The proof in [START_REF] Sandier | From the Ginzburg-Landau Model to Vortex Lattice Problems[END_REF] relied on the sophisticated techniques of "ball construction methods", which originated with Jerrard [Jer] and Sandier [San] in Ginzburg-Landau theory (see [START_REF] Sandier | Vortices in the Magnetic Ginzburg-Landau Model[END_REF]Chapter 4] for a presentation) and which are very much two-dimensional, since they exploit the conformal invariance of the Laplacian in dimension 2. It is not clear at all how to find a replacement for this "ball construction method" in dimension d ≥ 3, and thus not clear how to prove that the local energy density associated to W is bounded below then (although we do not pursue this, it should however be possible to show that W itself is bounded below with the same ideas we use here). Instead the approach via W works by avoiding this issue and replacing it with the use of Onsager's lemma, making it technically much simpler, at the price of a different definition and a less precise energy. However we saw that W and W have same minima (at least in dimension 2) so that we essentially reduce to the same limit minimization problems.

Main results: ground state

Our results on the ground state put on a rigorous ground the informal interpretation of the two-scale structure of minimizing configurations we have been alluding to in the introduction. To describe the behavior of minimizers at the microscopic scale we follow the same approach as in [START_REF] Sandier | 2D Coulomb gases and the renormalized energy[END_REF] and perform a blow-up: For a given (x 1 , . . . , x n ), we let x ′ i = n 1/d x i and

h ′ n (x ′ ) = w * n i=1 δ x ′ i -µ 0 (n -1/d x ′ ) , (2.18)
Note here that the associated electric field

∇h ′ n is in L p loc (R d , R d ) if and only if p < d d-1
, in view of the singularity in ∇w around each point. One of the delicate parts of the analysis (and of the statement of the results) is to give a precise averaged formulation, with respect to all possible blow-up centers in Σ, and in this way to give a rigorous meaning to the vague sentence "around almost any point x in the support of µ 0 there are approximately µ 0 (x) points per unit volume, minimizing W with background density µ 0 (x)".

Our formulation of the result uses the following notion, as in [START_REF] Sandier | From the Ginzburg-Landau Model to Vortex Lattice Problems[END_REF][START_REF] Sandier | 2D Coulomb gases and the renormalized energy[END_REF], which allows to embed (R d ) n into the set of probabilities on X = Σ × L p loc (R d , R d ), for some 1 < p < d d-1 . For any n and x = (x 1 , . . . , x n ) ∈ (R d ) n we let i n (x) = P νn , where ν n = n i=1 δ x i and P νn is the push-forward of the normalized Lebesgue measure on Σ by

x → x, ∇h ′ n (n 1/d x + •) .
Explicitly:

i n (x) = P νn = - Σ δ (x,∇h ′ n (n 1/d x+•)) dx. (2.19)
This way i n (x) is an element of P(X), the set of probability measures on X = Σ×L p loc (R d , R d ) (couples of blown-up centers, blown-up electric fields) that measures the probability of having a given blown-up electric field around a given blow-up point in Σ. As suggested by the above discussion, the natural object we should look at is really i n (x) and its limits up to extraction, P .

Due to the fact that the renormalized jellium functional describing the small-scale physics is invariant under translations of the electric field, and by definition of i n , we should of course expect the objects we have just introduced to have a certain translation invariance, formalized as follows:

Definition 2.5 (T λ(x) -invariance). We say a probability measure

P on X is T λ(x) -invariant if P is invariant by (x, E) → (x, E(λ(x) + •)), for any λ(x) of class C 1 from Σ to R d .
Note that T λ(x) -invariant implies translation-invariant (simply take λ ≡ 1).

Definition 2.6 (Admissible configurations). We say that P ∈ P(X) is admissible if its first marginal is the normalized Lebesgue measure on Σ, if it holds for P -a.e. (x, E) that E ∈ A µ 0 (x) , and if P is T λ(x) -invariant.

Assumption 2.2 ensures here that µ 0 (x) > 0 for any x ∈ Σ, so that the class A µ 0 (x) makes sense.

Our main result on the next order behavior of minimizers of H n is that min

x n 2/d-2 H n (x) -n 2 E[µ 0 ] + n 2 log n 1 d=2 → min P admissible W(P )
where we define

W(P ) := |Σ| c d X W(E) dP (x, E), (2.20)
if P is admissible, and +∞ otherwise. The function W is precisely what we meant by "the average with respect to blow-up centers" of W.

We note that, in view of the scaling relation (2.13), we may guess that min

P admissible W = ξ d :=        1 c d min A 1 W R d µ 2-2/d 0 if d ≥ 3 1 2π min A 1 W - 1 2 R 2 µ 0 log µ 0 if d = 2.
(2.21)

Note that the minima on the right-hand side exist thanks to Theorem 1. Also, the left-hand side is clearly larger than the right-hand side in view of the definition of W, that of admissible, and (2.13). That there is actually equality is a consequence of the following theorem on the behavior of ground state configurations (i.e. minimizers of (1.1)) :

Theorem 2 (Microscopic behavior of ground state configurations).

Let (x 1 , . . . , x n ) ∈ (R d ) n minimize H n . Let h ′ n be associated via (2.18) and P νn = i n (x 1 , . . . , x n ) be defined as in (2.19).
Up to extraction of a subsequence we have P νn ⇀ P in the sense of probability measures on X, where P is admissible and,

if d ≥ 3, lim n→∞ n 2/d-2 H n (x 1 , . . . , x n ) -n 2 E[µ 0 ] = W(P ) = ξ d = 1 c d min A 1 W µ 2-2/d 0 (x) dx, (2.22) respectively for d = 2, lim n→∞ n -1 H n (x 1 , . . . , x n ) -n 2 E[µ 0 ] + n 2 log n = W(P ) = ξ 2 = 1 c 2 min A 1 W - 1 2 µ 0 (x) log µ 0 (x) dx, (2.23)
P is a minimizer of W and E minimizes W over A µ 0 (x) for P -a.e. (x, E).

Note that the configuration of points itself depends on n, i.e. we have slightly abused notation by writing (x 1 , . . . , x n ) instead of (x n 1 , . . . , x n n ). The statement above implies (1.14) and is a rigorous formulation of our vague sentence in the introduction about the double scale nature of the charge distribution. Note that since P is always admissible, the result includes the fact that the blown-up field around x is in the class A µ 0 (x) , i.e. has a local density µ 0 (x).

If the crystallization conjecture is correct, and if the only minimizers of W are periodic configurations, this means that after blow-up around (almost) any point in Σ, one should see a crystalline configuration of points, the one minimizing W, packed at the scale corresponding to µ 0 (x). Note that in [RNS] a stronger result is proved in dimension 2: exploiting completely the minimality of the configuration, it is shown that this holds after blow up around any point in Σ (not too close to ∂Σ however), and that the renormalized energy density as well as the number of points are equidistributed (modulo the varying density µ 0 (x)). The same results are likely to be proveable in dimension d ≥ 3 combining the present approach and the method of [RNS].

The proof of Theorem 2 relies on the following steps:

• First we show the result as a more general lower bound, that is, we show that for an arbitrary configuration (x 1 , . . . , x n ), it holds that

lim inf n→∞ n 2/d-2 H n (x 1 , . . . , x n ) -n 2 E[µ 0 ] + n 2 log n 1 d=2 ≥ W(P ) -o(η)
with P = lim n→∞ P νn . This relies on the splitting formula alluded to before, based on Onsager's lemma, and the general ergodic framework introduced in [SS3, SS4] and suggested by Varadhan, which allows to bound first from below by W η instead of W, for fixed η.

• Taking the lim inf η→0 requires showing that W η is bounded from below by a constant independent of η (and then using Fatou's lemma). This is accomplished in a crucial step, which consists in proving that minimizers of W η have points that are well-separated: theirdistances are bounded below by a constant depending only on µ 0 L ∞ . This relies on an unpublished result of Lieb [Lie2], which can be found in dimension d = 2 in [START_REF] Nodari | Renormalized Energy Equidistribution and Local Charge Balance in 2D Coulomb Systems[END_REF]Theorem 4] and which we readapt here to our setting. This allows to complete the proof of the lower bound.

• We prove the corresponding upper bound inequality only at the level of the minimal energy. The reason is that the proof relies on an explicit construction, based on the result of the previous step: we take a minimizer of W and we "screen" it, as done in [START_REF] Sandier | From the Ginzburg-Landau Model to Vortex Lattice Problems[END_REF][START_REF] Sandier | 2D Coulomb gases and the renormalized energy[END_REF][START_REF] Sandier | 1D Log Gases and the Renormalized Energy: Crystallization at Vanishing Temperature[END_REF]. This means we truncate it over some large box and periodize it, in order to be able to copy-paste it after proper rescaling in order to create an optimal configuration of points. This screening uses crucially the preliminary result that the points are well separated. At this point, we see that arbitrary configurations cannot be screened: for example if a configuration has a multiple point, then truncating and periodizing it would immediately result in an infinite W, and obviously to an infinite H n ! We note also that if the upper bound inequality was holding for any general admissible P , this would mean that W is the complete Γ-limit (at next order) of H n . But in dimension 2 at least, it has been proven in [START_REF] Sandier | 2D Coulomb gases and the renormalized energy[END_REF] that the corresponding Γ-limit is W.

By uniqueness of the Γ-limit, this would mean that W dP = W dP for all admissible probabilities P . But we know that (at least in dimension 2) W and W are not equal, so there remains to know whether their average over admissible probabilities can coincide.

• Combining the lower bound and the matching upper bound immediately yields the result for minimizers.

Main results: finite temperature case

We now turn to the case of positive temperature and study the Gibbs measure (1.9). We shall need the common assumption that there exists

β 1 > 0 such that e -β 1 V (x)/2 dx < ∞ when d ≥ 3 e -β 1 ( V (x) 2 -log |x|) dx < ∞ when d = 2.
(2.24)

Note that this is only a slight strengthening of the assumption (1.5). By definition, the Gibbs measure minimizes the n-body free energy functional

F n,β [µ] := R dn µ(x)H n (x)dx + 2 β R dn µ(x) log(µ(x))dx (2.25)
over probability measures µ ∈ P(R dn ). It is indeed an easy and fairly standard calculation to show that the infimum is attained3 at P n,β (cf. (1.9)) and that we have

F n,β := inf µ∈P(R dn ) F n,β [µ] = F n,β [P n,β ] = - 2 β log Z β n . (2.26)
where Z β n is the partition function normalizing P n,β . As announced in the introduction, we will mainly focus on the low temperature regime β n 2/d-1 where we can get more precise results. We expect that crystallization occurs when β ≫ n 2/d-1 and that β ∝ n 2/d-1 should correspond to the liquid-crystal transition regime. Since they are the main basis of this conjecture, let us first state precisely our estimates on the n-body free energy defined above. In view of (2.26), estimating log Z β n and F n,β are equivalent and we shall work with the later. We need to introduce the mean field free energy functional

F[µ] = E[µ] + 2 nβ µ log µ (2.27)
with minimizer, among probabilities, µ β . It naturally arises when taking µ ⊗n as a trial state in (2.25). Note that for βn ≫ 1, this functional is just a perturbation of (1.8) and µ β agrees with µ 0 to leading order (results of this kind are presented in [RSY2, Section 3]). For convenience we will restrict the high temperature regime in Theorem 3 to β ≥ Cn -1 in order that the last term (entropy) in the mean-field free energy functional stays bounded above by a constant4 . The results below do not depend on C.

Theorem 3 (Free energy/partition function estimates).

The following estimates hold.

1. (Low temperature regime). Let β := lim sup n→+∞ βn 1-2/d , and assume 0 < β ≤ +∞.

There exists C β > 0 depending only on V and d, with

lim β→∞ C β = 0 and C ∞ = 0 such that lim sup n→∞ n 2/d-2 F n,β -n 2 E[µ 0 ] -n 2-2/d ξ d ≤ C β ,
(2.28)

respectively if d = 2 and β ≥ c(log n) -1 for n large enough with c > 0, lim sup n→∞ n -1 F n,β -n 2 E[µ 0 ] + n 2 log n -nξ 2 ≤ C β , (2.29)
where ξ d is as in (2.21).

(High temperature regime

). If d ≥ 3 and c 1 n -1 ≤ β ≤ c 2 n 2/d-1 for some c 1 , c 2 > 0, we have for n large enough, F n,β -n 2 E[µ β ] - 2n β R d µ β log µ β ≤ Cn 2-2/d , (2.30) respectively if d = 2 and c 1 n -1 ≤ β ≤ c 2 (log n) -1 for some c 1 , c 2 > 0 F n,β -n 2 E[µ β ] - 2n β R d µ β log µ β ≤ Cn log n, (2.31)
where C depends only on V and d.

The leading order contribution to (2.28)-(2.31) has been recently derived, along with the corresponding large deviation principle, in [CGZ] under the assumption β ≫ n -1 log n in our units. In this regime one may use either E[µ β ] or E[µ 0 ] for leading order considerations, but when β ∝ n -1 it is necessary to use the former.

In the regime β ≫ n 2/d-1 (2.28) says that the free energy agrees with the ground state energy (compare with (1.14)) up to a negligible remainder, and this is the only case where we prove the existence of a thermodynamic limit. We conjecture that this corresponds to a transition to a crystalline state. On the other hand, (2.30) shows that for β ≪ n 2/d-1 the effect of the entropy at the macroscopic scale prevails over the jellium energy at the microscopic scale. This alone is no indication that there is no crystallization in this regime: One may perfectly well imagine that a term similar to the subleading term in (1.14) appears at a further level of approximation (one should at least replace µ 0 by µ β but that would encode no significantly different physics at the microscopic scale). We believe that this is not the case and that the entropy should also appear at the microscopic scale, but this goes beyond what we are presently able to prove. Note that when d ≥ 3, the fixed β regime is deep in the low-temperature regime, and so we do not expect any particular transition to happen at this order of inverse temperature.

The case d = 2, in which we re-obtain the result of [START_REF] Sandier | 2D Coulomb gases and the renormalized energy[END_REF] is a little bit more subtle due to the particular nature of the Coulomb kernel, which is at the origin of the n log n term in the expansion. As in higher dimensions, comparing (1.14) and (2.29) we see that free energy and ground state energy agree to subleading order when β ≫ 1, which we conjecture to be the crystallization regime. The estimate (2.31) shows that the entropy is the subleading contribution only when β ≪ (log n) -1 however. It is not clear from Theorem 3 what exactly happens when (log n) -1 β 1. We expect entropy terms at both macroscopic and microscopic levels to enter. For the regime β ∝ 1 (the most studied regime, e.g. [BZ]), one can see the further discussion in [START_REF] Sandier | 2D Coulomb gases and the renormalized energy[END_REF] and Theorem 4 below.

Our next result exposes the consequences of (2.28)-(2.29) for the Gibbs measure itself. Roughly speaking, we prove that it charges only configurations whose renormalized energy W is below a certain threshold, in the sense that other configurations have exponentially small probability in the limit n → ∞. This is a large deviations type upper bound at speed n 2-2/d , but a complete large deviations principle is missing. The threshold of renormalized energy vanishes in the limit β ≫ n 2/d-1 , showing that the Gibbs measure charges only configurations that minimize the renormalized energy at the microscopic scale, which is a proof of crystallization, modulo the question of proving that minimizers of W are really crystalline. Our precise statement is again complicated by the double scale nature of the problem and the fact that the renormalized energy takes electric fields rather than charge configurations as argument. Its phrasing uses the same framework as Theorem 2 in terms of the limiting probability measures P , which should now be seen as random, due to the temperature. In fact we also prove the existence of a limiting "electric field process", as n → ∞, i.e. a limiting probability on the P 's, which is the limit of the push-forwards of P n,β by i n .

We will consider the limit of the probability that the system is in a state (x 1 , . . . , x n ) ∈ A n for a given sequence of sets A n ⊂ (R d ) n in configuration space. Associated to this sequence we introduce

A ∞ = {P ∈ P(X) : ∃x n ∈ A n , P νn ⇀ P up to a subsequence} (2.32)
where P νn is as in (2.19) with x = x n and the convergence is weakly as measures. A ∞ can be thought of as the limsup in the sense of sets of the sequence of sets i n (A n ).

Theorem 4 (Microscopic behavior of thermal states in the low temperature regime).

For any n > 0 let A n ⊂ (R d ) n and A ∞ be as above. Let β > 0 be as in Theorem 3, ξ d as in (2.21). There exists C β depending only on V and d such that C β = 0 for β = ∞ and lim sup

n→∞ log P n,β (A n ) n 2-2/d ≤ - β 2 inf P ∈A∞ W(P ) -ξ d -C β .
(2.33)

Moreover, for fixed β > 0, letting P n,β denote the push-forward of P n,β by i n (defined in (2.19)), { P n,β } n is tight and converges as n → ∞, up to a subsequence, to a probability measure on P(X) which is concentrated on admissible probabilities satisfying W(P ) ≤ ξ d +C β .

Note that the error term C β becomes negligible when β ≫ n 2/d-1 . This means as announced that the Gibbs measure concentrates on minimizers of W in that regime. When β = βn 2/d-1 , we have instead a threshold phenomenon: the Gibbs measure concentrates on configurations whose W is below the minimum plus C β . For the 2D case, a partial converse to (2.33) is proved in [START_REF] Sandier | 2D Coulomb gases and the renormalized energy[END_REF], establishing a kind of large deviation principle at speed n in the regime β ≫ 1, but with W instead of W. A full LDP should involve an additional entropy term and is still missing.

Estimates on deviations from mean-field theory

Our methods also yield some quantitative estimates on the fluctuations from the equilibrium measure. We again focus on the low temperature regime for concreteness but if need arises our methods could work just as well in the opposite regime, with less hope of optimality however. The following estimates are insufficient for a crystallisation result, as noted before, but they set bounds on the possible deviations from mean-field theory in a stronger sense than what has been considered so far, namely by estimating charge deviations and correlation functions.

We will estimate the following quantity, defined for any configuration (x 1 , . . . ,

x n ) ∈ R dn D(x, R) := ν n (B(x, R)) -n B(x,R) µ 0 , (2.34)
i.e. the deviation of the charge contained in a ball B(x, R) with respect to the prediction drawn from the equilibrium measure. Here

ν n = n i=1 δ x i . (2.35)
Note that since µ 0 is a fixed bounded function, the second term in (2.34) is of order nR d , and typically the difference of the two terms is much smaller than this, since the distribution of points at least approximately follows the equilibrium measure.

For pedagogical purposes we only state two typical results:

1. the probability of a large deviation of the number of points in a ball B(x, R) (i.e. a deviation of order nR d ) is exponentially small as soon as R n -1/(d+2) .

2. the probability of a charge deviation of order n 1-1/d in a macroscopic ball B(x, R) with R = O(1) is exponentially small.

We deduce these results from sharper but more complicated estimates similar to [SS4, Eq.

(1.47) and (1.49)], see Section 7.3. Actually these contain a whole continuum of results of the type above: the larger the scale at which one considers the deviations, the smaller they are allowed to be. If need be, any scale R satisfying n -1/(d+2) R ≤ C can be considered, and deviations on this scale are exponentially small, with a rate depending on the scale. We believe that the two results mentioned above are a sufficient illustration of our methods and we proceed to state them rigorously, along with an estimate on the discrepancy between ν n and the equilibrium measure nµ 0 in weak Sobolev norms. Below W -1,q (Ω) denotes the dual of the Sobolev space W 1,q ′ 0 (Ω) with 1/q + 1/q ′ = 1, in particular W -1,1 is the dual of Lipschitz functions.

Theorem 5 (Charge fluctuations). Assume there is a constant C > 0 such that β ≥ Cn 2/d-1 . Then the following holds for any x ∈ R d .

1. (Large fluctuations on microscopic scales). Let R n be a sequence of radii satisfying

R n ≥ C R n -1 d+2 for some constant C R .
Then, for any λ > 0 we have, for n large enough,

P n,β |D(x, R n )| ≥ λnR d n ≤ Ce -Cβn 2-2/d (C R λ 2 -C) , (2.36)
for some C depending only on dimension.

2. (Small fluctuations at the macroscopic scale). Let R > 0 be a fixed radius. There is a constant C depending only on dimension such that for any λ > 0, for n large enough,

P n,β |D(x, R)| ≥ λn 1-1/d ≤ Ce -Cβn 2-2/d (min(λ 2 R 2-d ,λ 4 R 2-2d )-C) .
(2.37)

3. (Control in weak Sobolev norms). Let R > 0 be some fixed radius. Let 1 ≤ q < d d-1 and define

t q,d = 2 - 1 d - 1 q < 1 tq,d = 3 - 1 d - 2 dq > 0.
There is a constant C R > 0 such that the following holds for n large enough, and any λ large enough,

P n,β ν n -nµ 0 W -1,q (B R ) ≥ λn t q,d ≤ Ce -βC R λ 2 n tq,d . (2.38)
As regards (2.38), note that only estimates in W -1,q norm with q < d d-1 make sense, since a Dirac mass is in W -1,q if and only if q < d d-1 (hence the same for ν n ). In view of the values of the parameters t q,d and tq,d defined above, our results give meaningful estimates for any such norm: a large deviation in a ball of fixed radius would be a deviation of order n.

Since t q,d < 1, equation (2.38) above implies that such deviations are exponentially unlikely in W -1,q for any q < d d-1 , in particular in W -1,1 for any space dimension. The proof of Item 3 is based on a control in W -1,2 of νn -µ 0 , where νn is the regularization of ν n with charges smeared-out on a scale n -1/d . This is another instance where the method of smearing out charges makes the proof siginificantly easier than in [START_REF] Sandier | 2D Coulomb gases and the renormalized energy[END_REF], replacing the use of a "displaced" energy density [START_REF] Sandier | Improved lower bounds for Ginzburg-Landau energies via mass displacement[END_REF] and Lorentz space estimates [SeTi].

Finally we state some consequences for the marginals (reduced densities) of the probability (1.9). Let us denote

P (k) n,β (x 1 , . . . , x k ) = x ′ ∈R d(n-k) P n,β (x 1 , . . . , x k , x ′ ) dx ′ .
(2.39)

Remark that since P n,β is symmetric w.r.t. exchange of variables, it does not matter over which n -k particles we integrate to define

P (k)
n,β (particles are indistinguishable). The value P (k) n,β (x 1 , . . . , x k ) is interpreted as the probability density for having one particle at x 1 , one particle at x 2 , . . ., and one particle at x k .

Corollary 6 (Marginals of the Gibbs measure in the low temperature regime).

Let R > 0 be some fixed radius, 1 ≤ q < d d-1 . Under the same assumptions as Theorem 5, there exists a constant C > 0 depending only on the dimension such that the following holds:

• (Estimate on the one-particle reduced density).

P (1) n,β -µ 0 W -1,q (B R ) ≤ Cn 1-1/d-1/q = o n (1).
(2.40)

• (Estimate on k-particle reduced densities). Let k ≥ 2 and ϕ : R dk → R be a smooth function with compact support, symmetric w.r.t. particle exchange. Then we have

R dk P (k) n,β -µ ⊗k 0 ϕ ≤ C kn 1-1/d-1/q + k 2 n -1 sup x 1 ∈R d . . . sup x k-1 ∈R d ∇ϕ(x 1 , . . . , x k-1 , . ) L p (R d ) , (2.41)
where 1/p = 1 -1/q.

As our other results, Corollary 6 concerns the low temperature regime. One may use the same technique to estimate the discrepancy between the n-body problem and mean-field theory in other regimes. When βn becomes small however (large temperature), in particular when β ∼ n -1 so that entropy and energy terms in (2.27) are of the same order of magnitude, a different method can give slightly better estimates. Since this regime has been considered for related models in several works [START_REF] Messer | Statistical mechanics of the isothermal Lane-Emden equation[END_REF][START_REF] Kiessling | On the Equilibrium Statistical Mechanics of Isothermal Classical Self-Gravitating Matter[END_REF][START_REF] Kiessling | Statistical mechanics of classical particles with logarithmic interactions[END_REF][START_REF] Caglioti | A Special Class of Stationary Flows for Two-Dimensional Euler Equations: A Statistical Mechanics Description[END_REF], it is worth mentioning that quantitative estimates in the spirit of (2.40)-(2.41) can be obtained in total variation norm. We refer to Remark 7.3 in Section 7.1 for details on this approach, which is in a slightly different spirit from what we have presented so far.

Organization of the paper

In Section 3, we smear out the charges at scale η for fixed η and we use Onsager's lemma to obtain a sharp "splitting formula" for the Hamiltonian, in which the leading order and next order contributions decouple. We also obtain a control on the charge fluctuations and on the electric field in terms of the next order term in the Hamiltonian.

In Section 4 we start taking the limit η → 0 in the estimates of the previous section, and prove the lower bound on the energy via the ergodic framework of Varadhan presented in [START_REF] Sandier | From the Ginzburg-Landau Model to Vortex Lattice Problems[END_REF]. This assumes the lower bound on W η .

In Section 5, we use the result of [Lie2] to reduce to points that are well-separated and deduce a lower bound on W η independent on η. This requires a "screening result", which will also be used in Section 6, where we prove an upper bound for the minimal energy by the construction of a precise test-configuration. This finishes the proof of Theorem 2.

In Section 7, we apply all the previous results to the case with temperature, and deduce Theorems 3, 4 and 5.

Appendices A and B respectively contain a direct proof of the existence of a minimizer for the renormalized energy functional and a discussion of the relation of the two versions of the renormalized energy functional.

Splitting formulae and control on fluctuations

In this section, we start to exploit the idea of smearing out the charges and Onsager's lemma in a way similar to [START_REF] Rougerie | Quantum Hall phases and plasma analogy in rotating trapped Bose gases[END_REF]. We also explore easy corollaries that can be obtained for fixed η. Since the smearing our procedure turns out to have no effect on the energy of configurations with well-separated points we start by discussing these particular configurations in the following subsection.

Preliminaries and well-separated configuraions

We start with a lemma that shows that if E ∈ A m and W(E) < ∞ then the density of points is indeed equal to that of the neutralizing background, i.e. m. From now on, for any E ∈ A m we denote by ν the corresponding measure of singular charges, i.e.

p∈Λ N p δ p .

Lemma 3.1 (Density of points in finite-energy configurations).

Let E ∈ A m be such that W η (E) < ∞ for some η ≤ 1, and let ν = -div E + m. Then we have lim R→∞ ν(K R ) |K R | = m. 5
Proof. First we show that

ν(K R-2 ) ≤ m|K R |+ CR d-1 2 E η L 2 (K R ) ν(K R+1 ) ≥ m|K R-1 |-CR d-1 2 E η L 2 (K R ) . (3.1)
To prove this, first by a mean value argument, we find t

∈ [R -1, R] such that ∂Kt |E η | 2 ≤ K R |E η | 2 . (3.2)
Let us next integrate (2.10) over K t and use Stokes's theorem to find

Kt p∈Λ N p δ (η) p -m|K t | = - ∂Kt E η • ν, (3.3)
where ν denotes the outer unit normal. Using the Cauchy-Schwarz inequality and (3.2), we deduce that Kt p∈Λ

N p δ (η) p -m|K t | ≤ CR d-1 2 E η L 2 (K R ) . (3.4)
Since η ≤ 1, by definition of ν and since the δ

(η) p are supported in B(p, η), we have ν(K R-2 ) ≤ Kt p∈Λ N p δ (η) p ≤ ν(K R+1 ) in view of the definition of ν = N p δ p . The claim (3.1) follows. Since W η (E) < +∞ we have K R |E η | 2 ≤ C η R d for any R > 1.
Inserting this into (3.1), dividing by |K R | and letting R → ∞, we easily get the result.

We next turn to configurations with well-separated points. We will need the following scaling relation, which can be obtained from (2.3) and (2.7) by a change of variables:

D(δ (ℓ) 0 , δ (ℓ) 0 ) = κ d c d w(η) if d ≥ 3 D(δ (ℓ) 0 , δ (ℓ) 0 ) = w(η) + γ 2 c 2 = κ 2 c 2 w(η) + γ 2 c 2 if d = 2. (3.5)
Lemma 3.2 (The energy of well-separated configurations).

Assume that E = ∇h satisfies

-∆h = c d p∈Λ δ p -a(x) (3.6)
in some subset U ⊂ R d , for some a(x) ∈ L ∞ (U ), and Λ a discrete subset of U , and min min

p =p ′ ∈Λ |p -p ′ |, min p∈Λ dist(p, ∂U ) ≥ η 0 > 0. (3.7)
Then, we have

U |E η | 2 -#(Λ ∩ U )(κ d w(η) + γ 2 1 d=2 ) = W (E, 1 U ) + #(Λ ∩ U )o η (1) a L ∞ (U ) , (3.8)
where o η (1) → 0 as η → 0 is a function that depends only on the dimension. Moreover,

W (E, 1 U ) ≥ -C#(Λ ∩ U ), (3.9)
where C > 0 depends only on the dimension, γ 2 (hence the choice of smearing function ρ), a L ∞ and η 0 .

Proof. We recall that by definition of E η (cf. Definition 2.2) we have

E η = E + p∈Λ ∇f η (x -p).
Since the B(p, η 0 ) are disjoint and included in U , and f η is identically 0 outside of B(0, η) we may write for any η < η 0 , and any 0 < α < η,

U \∪ p∈Λ B(p,α) |E η | 2 = U \∪ p∈Λ B(p,α) |E| 2 + #(Λ ∩ U ) B(0,η)\B(0,α) |∇f η | 2 + 2 p∈Λ B(p,η)\B(p,α) ∇f η (x -p) • E. (3.10)
First we note that, using Green's formula, and ν denoting the outwards pointing unit normal to ∂B(0, α) we have

B(0,η)\B(0,α) |∇f η | 2 = - ∂B(0,α) f η ∇f η • ν + c d B(0,η)\B(0,α) f η δ (η) 0 .
By Green's formula again and the definition of f η we have

∂B(0,α) ∇f η • ν = -c d B(0,α) δ (η) 0 + c d = c d + o α (1)
as α → 0, and combining with the fact that f η = w * δ (η) 0 -w (see its definition (2.8)) we find

B(0,η)\B(0,α) |∇f η | 2 = -c d f η (α) + c d R d w * δ (η) 0 δ (η) 0 -c d R d wδ (η) 0 + o α (1). ( 3
|∇f η | 2 = -c d w * δ (η) 0 (α) + c d w(α) + κ d w(η) -c d w * δ (η) 0 (0) + o α (1) = -2c d w * δ (η) 0 (0) + c d w(α) + c d D(δ (η) 0 , δ (η) 0 ) + o α (1), (3.12)
in view of the fact that for fixed η, w * δ (η) 0 is continuous at 0. On the other hand, using Green's theorem and (3.6) we have

B(p,η)\B(p,α) ∇f η (x -p) • E = -c d B(p,η)\B(p,α) f η (x -p)a(x) dx -f η (α) ∂B(p,α) E • ν.
First we note that

B(p,η)\B(p,α) f η (x -p)a(x) dx ≤ c d a L ∞ B(0,η) |f η |(x)dx ≤ a L ∞ o η (1),
where o η (1) depends only on ρ and d. To see this, just notice that we have |f η | ≤ |w| and the Coulomb kernel w is integrable near the origin. Secondly, by Green's theorem again we have

- ∂B(p,α) E • ν = c d + O( a L ∞ α d ).
Inserting these two facts we deduce

B(p,η)\B(p,α) ∇f η (x -p) • E = c d f η (α) + O( a L ∞ α d w(α)) + a L ∞ o η (1).
Combining this and (3.12), (3.10), (3.5), and again f η (α) = w * δ (η) 0 (0) -w(α) + o α (1), we find

U \∪ p∈Λ B(p,α) |E η | 2 = U \∪ p∈Λ B(p,α) |E| 2 + #(Λ ∩ U )(κ d w(η) + γ 2 1 d=2 -c d w(α) + o α (1)) + a L ∞ o η (1)#(Λ ∩ U ) + O( a L ∞ α d f η (α))#(Λ ∩ U ).
Letting α → 0, (3.8) follows by the definition (2.17).

The proof of (3.9) is a toy version of "ball construction" lower bounds in Ginzburg-Landau theory, made much simpler by the separation of the points. From (3.6) and the Cauchy-Schwarz inequality, we have, for any p ∈ Λ,

B(p, 1 2 η 0 )\B(p,η) |E| 2 ≥ η 0 /2 η 1 |S d-1 |t d-1 ∂Bt E • ν 2 dt ≥ c 2 d η 0 /2 η 1 |S d-1 |t d-1 (1 -a L ∞ |B d |t d ) 2 dt ≥ c d (w(η) -w(η 0 /2)) -C (3.13)
where |B d | is the volume of the unit ball in dimension d, and we have used the definition of c d , and C depends only on a L ∞ and d. We may then absorb c d w(η 0 /2) into a constant C > 0 depending only on a L ∞ , η 0 and d. Since the B(p, 1 2 η 0 ) are disjoint and included in U , we may add these lower bounds, and obtain the result.

Splitting lower bounds

We start by discussing the problem of minimization of E defined in (1.8). Direct variations of the form (1 -t)µ 0 + tν for ν ∈ P(R d ) and t ∈ [0, 1] yield that the unique minimizer of E (note that E is strictly convex), denoted µ 0 , solves

h µ 0 + V 2 = c := 1 2 (E[µ 0 ] + D(µ 0 , µ 0 )) q.e. in Supp µ 0 h µ 0 + V 2 ≥ c := 1 2 (E[µ 0 ] + D(µ 0 , µ 0 )) q.e.
(3.14)

where q.e. means "outside of a set of capacity 0", and h µ 0 = w * µ 0 is the potential generated by µ 0 . More precisely, the variations first yield that (3.14) holds for some constant c on the right-hand side of both relations. Then, integrating the first relation against µ 0 yields that

c = R d (h µ 0 + V 2 ) dµ 0 = D(µ 0 , µ 0 ) + V 2 dµ 0
which identifies c as the right-hand side in (3.14). It turns out that these relations can also be shown to characterize uniquely µ 0 (by convexity). For more details, as well as a proof of the existence of µ 0 if the assumptions (1.5) are verified, one can see [Fro], [START_REF] Saff | Logarithmic Potentials with External Fields[END_REF]Chap. 1] in dimension 2, or [Ser]. It turns out, although this is rarely emphasized in the literature, that the solution µ 0 is also related to an obstacle problem, in the following sense: the potential h µ 0 generated by µ 0 can be shown (for a proof, one can refer to [ASZ], the proof is presented in dimension 2 and for V quadratic but carries over to dimension d ≥ 3 and general V with no change) to solve the following variational inequality:

∀u ∈ K, R d ∇h µ 0 • ∇(u -h µ 0 ) ≥ 0 (3.15)
where

K = {u ∈ H 1 loc (R d ), u -h µ 0 has bounded support and u ≥ - V 2 + c, q.e.}.
This happens to characterize a classical obstacle problem with obstacle ϕ := -V 2 + c (for general background on the obstacle problem, one can see [KiSt]). We then denote

ζ = h µ 0 + V 2 - 1 2 (E[µ 0 ] + D(µ 0 , µ 0 )) , (3.16)
and note that in view of (3.14), ζ ≥ 0 and ζ = 0 in Σ. Because h µ 0 is a solution to the obstacle problem with obstacle ϕ, the support of µ 0 , that we denote Σ, is contained in the so-called "coincidence set" where h µ 0 = ϕ i.e. the set {ζ = 0}. Equality between these sets happens for example V ∈ C 2 and ∆V does not vanish on the coincidence set. We observe that in all cases the set {ζ = 0} is bounded since h µ 0 behaves like w at infinity and (1.5) holds.

We next turn to recalling our blow-up procedure:

x ′ = n 1/d x, x ′ i = n 1/d x i .
For a configuration of points (x 1 , . . . , x n ) we let here and in the sequel, as in (2.18),

h n = w * n i=1 δ x i -nµ 0 h ′ n (x ′ ) = n 2/d-1 h n (x) = w * n i=1 δ x ′ i -µ 0 (n -1/d x ′ ) , (3.17)
and for ℓ and η related by ℓ = n -1/d η, we let

h n,ℓ := w * n i=1 δ (ℓ) x i -nµ 0 h ′ n,η (x ′ ) = n 2/d-1 h n,ℓ (x) = w * n i=1 δ (η) x ′ i -µ 0 (n -1/d x ′ ) .
(3.18) We also denote

ν ′ n = n i=1 δ x ′ i . (3.19)
Our first splitting result is the equivalent of [SS4, Lemma 2.1].

Lemma 3.3 (Exact splitting).

For any n ≥ 1 and any (x 1 , . . . , x n ) ∈ (R d ) n , letting h ′ n be as in (3.17), we have

H n (x 1 , . . . , x n ) = n 2 E[µ 0 ] + 2n n i=1 ζ(x i ) - 1 2 n log n + 1 c d W (∇h ′ n , 1 R 2 ) if d = 2 = n 2 E[µ 0 ] + 2n n i=1 ζ(x i ) + n 1-2/d c d W (∇h ′ n , 1 R d ) if d ≥ 3. (3.20)
Proof. Exactly as in [SS4, Lemma 2.1], we can show that

H n (x 1 , . . . , x n ) = n 2 E[µ 0 ] + 2n n i=1 ζ(x i ) + 1 c d W (∇h n , 1 R d ).
(3.21)

A change of variables yields that

R d \∪ n i=1 B(x i ,η) |∇h n | 2 = n 1-2/d R d \∪ n i=1 B(x ′ i ,ηn 1/d ) |∇h ′ n | 2
and subtracting off nc d w(η) from both sides, writing w(η) = w(ηn 1/2 ) + 1 2 log n in dimension 2 and w(η) = w(ηn 1/2 )n d/2-1 in dimension d ≥ 3, and letting η → 0, we are led to

W (∇h n , 1 R d ) = W (∇h ′ n , 1 R d ) -c d n 2 log n if d = 2 W (∇h n , 1 R d ) = n 1-2/d W (∇h ′ n , 1 R d ) if d = 3. (3.22)
Inserting this into (3.21) yields the result.

In [START_REF] Sandier | 2D Coulomb gases and the renormalized energy[END_REF] W was then bounded below locally. As already mentioned, we do not know how to do this in dimension d ≥ 3 (due to the lack of an efficient "ball construction" method) so we resort to a second type of splitting, using the smearing out of charges.

A crucial ingredient in our approach is that the electrostatic energy of a configuration of positive smeared charges is always a lower bound for the energy of the corresponding configuration of point charges, with equality if the smeared charges do not overlap. This is the so-called Onsager lemma, which one can find for example in [START_REF] Lieb | The stability of matter in quantum mechanics[END_REF]Chapter 6]. We reproduce it here for sake of completeness.

Lemma 3.4 (Onsager's lemma).

For any nonnegative distribution µ such that R d µ = n, any n, any x 1 , . . . , x n ∈ R d , and any ℓ > 0, we have

i =j w(x i -x j ) ≥ D µ - i δ (ℓ) x i , µ - n i=1 δ (ℓ) x i -D(µ, µ) + 2 n i=1 D(µ, δ (ℓ) x i ) -nD(δ (ℓ) 0 , δ (ℓ) 0 ) (3.23) with equality if B(x i , ℓ) ∩ B(x j , ℓ) = ∅, for all i = j.
Proof. It is based on Newton's theorem, see [START_REF] Lieb | Analysis[END_REF]Theorem 9.7], which can easily be adapted to any dimension, and asserts in particular that for any ℓ, w * δ

(ℓ) 0 ≤ w = w * δ 0 pointwise. It follows from this that i =j w(x i -x j ) ≥ i =j D(δ (ℓ) x i , δ (ℓ) x j ) (3.24)
with equality if min i =j |x i -x j | ≥ 2ℓ. Indeed, by Newton's theorem we have

w * δ (ℓ) x i δ (ℓ) x j ≤ (w * δ x i ) δ (ℓ) x j = w * δ (ℓ) x j δ x i ≤ w * δ x j δ x i .
We then write

D µ - n i=1 δ (ℓ) x i , µ - n i=1 δ (ℓ) x j = D(µ, µ) -2 n i=1 D(µ, δ (ℓ) x i ) + i =j D(δ (ℓ) x i , δ (ℓ) x j ) + n i=1 D(δ (ℓ) x i , δ (ℓ) x i )
and from this relation and (3.24), the lemma easily follows.

Smearing-out charges comes with a cost, that we quantify in the following lemma.

Lemma 3.5 (The cost of smearing charges out).

For any µ ∈ L ∞ (R d ) and any point x, we have

D µ, δ x -δ (ℓ) x ≤ Cℓ 2 µ L ∞ , (3.25)
where C depends only on the choice of ρ in (2.6) and the dimension.

Proof. Without loss of generality, we may assume x = 0, and write

D(µ, δ 0 -δ (ℓ) 0 ) = (w -w * δ (ℓ) 0 )(x)µ(x) dx.
By Newton's theorem, the function w -w * δ

(ℓ)
0 is nonnegative and supported in B(0, ℓ). In dimension 3, we may just write that in particular it is smaller than w1 B(0,ℓ) , and so we may write D(µ, δ 0 -δ

(ℓ) 0 ) ≤ µ L ∞ B(0,ℓ) dx |x| d-2 ≤ Cℓ 2 µ L ∞ .
In dimension 2, we write

R 2 |w -w * δ (ℓ) 0 | = B(0,ℓ) -log |x| + 1 ℓ 2 log |x -y|ρ( y ℓ ) dy dx = ℓ 2 B(0,1) -log |x ′ | + log |x ′ -y ′ |ρ(y ′ ) dy ′ dx = Cℓ 2
where we have used the changes of variables x = ℓx ′ , y = ℓy ′ , and we conclude in the same way.

From these two lemmas we deduce Lemma 3.6 (Splitting lower bound).

For any n, any x 1 , . . . , x n ∈ R d , letting h n,ℓ be as in (3.18), we have

H n (x 1 , . . . , x n ) ≥ n 2 E[µ 0 ]+ 1 c d R d |∇h n,ℓ | 2 -n(κ d w(ℓ) + γ 2 1 d=2 ) +2n n i=1 ζ(x i )-Cn 2-2/d η 2 ,
(3.26) where C depends only on the dimension. Moreover, equality holds if min i =j |x i -x j | ≥ 2ℓ.

Note that compared to Lemma 3.3 it is an inequality, and not an equality, and it has an error term, however it achieves the same goal. Indeed, as we will see, points in minimizing configurations are well-separated so that there is equality in (3.26) and no information is lost in the end.

Proof. We proceed as in [RSY2, Proof of Thm 3.2]. First, applying Onsager's lemma above with ℓ and µ = nµ 0 , and using (3.5), we find

i =j w(x i -x j ) ≥ D nµ 0 - n i=1 δ (ℓ) x i , nµ 0 - n i=1 δ (ℓ) x i -n 2 D(µ 0 , µ 0 ) + 2n n i=1 D(µ 0 , δ (ℓ) x i ) -nD(δ (ℓ) 0 , δ (ℓ) 0 ), (3.27)
with equality if the B(x i , ℓ) are disjoint. We then use Lemma 3.5, the fact that µ 0 is a fixed L ∞ function and (3.16) to obtain

2n n i=1 D(µ 0 , δ (ℓ) x i ) = 2n n i=1 D(µ 0 , δ x i ) + O n 2 ℓ 2 = 2n n i=1 h µ 0 (x i ) + O n 2 ℓ 2 = n 2 (E[µ 0 ] + D(µ 0 , µ 0 )) + 2n n i=1 ζ(x i ) -n n i=1 V (x i ) + O n 2 ℓ 2 . (3.28)
Combining (3.27) and (3.28), and observing that in view of (2.4) and (3.18), we have

D nµ 0 - n i=1 δ (ℓ) x i , nµ 0 - n i=1 δ (ℓ) x i = 1 c d R d |∇h n,ℓ | 2
we thus obtain

H n (x 1 , . . . , x n ) = i =j w(x i -x j ) + n n i=1 V (x i ) ≥ n 2 E[µ 0 ] + 1 c d R d |∇h n,ℓ | 2 + 2n n i=1 ζ(x i ) -nD(δ (ℓ) 0 , δ (ℓ) 0 ) + O n 2 ℓ 2 .
Since ℓ = n -1/d η and (3.5) holds, we obtain the result.

Combining this lower bound with our blow-up and noting that by a change of variables we have

R d |∇h n,ℓ | 2 = n 1-2/d R d |∇h ′ n,η | 2 . (3.29)
we obtain the following rephrasing:

H n (x 1 , . . . , x n ) -n 2 E[µ 0 ] + n 2 log n 1 d=2 ≥ n 2-2/d J n (x 1 , . . . , x n ) -Cη 2 + 2n n i=1 ζ(x i )
(3.30) where we have written

J n (x 1 , . . . , x n ) := 1 c d 1 n R d |∇h ′ n,η | 2 -(κ d w(η) + γ 2 1 d=2 ) .
Remark 3.7. Taking for example η = 1, and using that ζ ≥ 0, it immediately follows from the above that

H n (x 1 , . . . , x n ) -n 2 E[µ 0 ] + n 2 log n 1 d=2 ≥ -Cn 2-2/d (3.31)
where the constant depends only on the dimension. This provides a very simple proof of this fact.

Controlling fluctuations

The key ingredient in the proof of Theorems 5 is the fact that the energy density (square of the local L 2 norm) of the electric field generated by the smeared charges provides a control of the deviations we are interested in. We start by formalizing this idea in two lemmas which, thanks to the smearing out method, provide significantly simpler alternatives to the estimates of [START_REF] Sandier | 2D Coulomb gases and the renormalized energy[END_REF][START_REF] Serfaty | Lorentz Space Estimates for the Coulombian Renormalized Energy[END_REF]. Note that the following considerations do not require that µ 0 be the equilibrium measure.

We start with a control on the fluctuations of the number of points (i.e. total charge) in a given ball:

D(x ′ , R) = ν ′ n (B(x ′ , R)) - B(x ′ ,R) µ ′ 0 . (3.32)
where µ ′ 0 = µ 0 (n -1/d •). Note that the quantity (2.34) used in Theorem 5 is

D(x, R) = D(x ′ , Rn 1/d ) (3.33)
where x ′ is the blown-up of x. For this reason we do not feel the need for a new notation and we will use the above for fluctuations either at the microscopic or macroscopic scale.

Lemma 3.8 (Controlling charge fluctuations).

For any x 1 , . . . , x n and h ′ n,η given by (3.17), for any 0 < η < 1, R > 2 and x ′ ∈ R d , we have

R d |∇h ′ n,η | 2 ≥ C D(x ′ , R) 2 R d-2 min 1, D(x ′ , R) R d , (3.34)
where C is a constant depending only on d.

Proof. In the proof, C will denote a constant depending only on d that may change from line to line. We distinguish two cases, according to whether D := D(x ′ , R) > 0 or D ≤ 0, and start with the former. We first claim that for all t such that

R + η ≤ t ≤ T := (R + η) d + D 2C 1/d (3.35)
for some well-chosen constant C, it holds that

- ∂B(x ′ ,t) ∇h ′ n,η • ν ≥ c d 2 D (3.36)
where ν is the outwards pointing normal to ∂B(x ′ , t). Indeed, by Green's formula

- ∂B(x ′ ,t) ∇h ′ n,η • ν = B(x ′ ,t) -∆h ′ n,η = c d B(x ′ ,t) n i=1 δ (η) x ′ i -µ ′ 0 ≥ c d D(x ′ , R) -c d B(x ′ ,t)\B(x ′ ,R) µ ′ 0 ≥ c d D(x ′ , R) -C t d -R d ≥ c d 2 D
if t satisfies (3.35). We have used the positivity of the smeared charges and the fact that µ ′ 0 is bounded in L ∞ (2.2). Then, integrating in spherical coordinates and using the Cauchy-Schwarz inequality as in (3.13), we find

R d |∇h ′ n,η | 2 ≥ B(x ′ ,T )\B(x ′ ,R+η) |∇h ′ n,η | 2 ≥ T t=R+η 1 t d-1 |S d-1 | ∂B(x ′ ,t) ∇h ′ n,η • ν 2 dt ≥ CD 2 T t=R+η 1 t d-1 |S d-1 | = CD 2 log T R+η if d = 2 1 (R+η) d-2 -1 T d-2 if d ≥ 3.
There only remains to note that in view of (3.35) and the fact that η < R/2, we have

log T R + η = log 1 + D 2C(R + η) d 1/d ≥ C ′ min(1, D R d ) 1 (R + η) d-2 - 1 T d-2 = 1 (R + η) d-2 1 -1 + D 2C(R + η) d d d-2 ≥ C ′ R d-2 min(1, D R d )
to conclude the proof in the case D > 0.

If D ≤ 0 the computation is the same, except that we set

T = (R -η) d - D 2C 1/d
, and use that for any

T ≤ t ≤ R -η - ∂B(x ′ ,t) ∇h ′ n,η • ν = B(x ′ ,t) -∆h ′ n,η ≤ c d D -C(R d -t d ) ≤ c d 2 D so that ∂B(x ′ ,t) ∇h ′ n,η • ν 2 ≥ CD 2
again. We then integrate the energy density on B(x ′ , R-η)\B(x ′ , T ) and argue as before.

The method of smearing out charges also provides a very convenient way to control the L q norms of the electric field. Lemma 3.9 (Controlling electric field fluctuations). For any x 1 , . . . , x n and h ′ n,η , h ′ n given by (3.17), for any 1 ≤ q < d d-1 , any 1 > η > 0 and any R > 0, we have

∇h ′ n L q (K R ) ≤ |K R | 1/q-1/2 ∇h ′ n,η L 2 (K R ) + C q,η ν ′ n (K R+η ) (3.37) ≤ C(R d q -d 2 + R d-1 2 ) ∇h ′ n,η L 2 (K R+2 ) + C q,η µ 0 L ∞ R d , (3.38)
where C q,η depends only on q and η and satisfies C q,η → 0 when η → 0 at fixed q, and C depends only on d.

Proof. Recalling the definition (2.8) we have

∇h ′ n = ∇h ′ n,η - n i=1 ∇f η (x -x i )
and thus

∇h ′ n L q (K R ) ≤ ∇h ′ n,η L q (K R ) + ν ′ n (K R+η ) ∇f η L q (R d )
where we used that if x ∈ K R and η < 1, then

f η (x -x i ) = 0 if x i ∈ (K R+η ) c
. A simple application of Hölder's inequality then yields

∇h ′ n,η L q (K R ) ≤ |K R | 1/q-1/2 ∇h ′ n,η L 2 (K R )
and concludes the proof of the first inequality, with C q,η := ∇f η L q (R d ) . Bounding then ν ′ n (K R+η ) exactly as we did for (3.1), we deduce

∇h ′ n L q (K R ) ≤ |K R | 1/q-1/2 ∇h ′ n,η L 2 (K R ) + C q,η µ 0 L ∞ |K R | + CR d-1 2 ∇h ′ n,η L 2 (K R+2 ) ,
and (3.38) follows.

4 Lower bound to the ground state energy

In this section, we will start to take the limits n → ∞ and η → 0, and we provide the lower bound part of Theorem 2 by proving the following:

Proposition 4.1 (Lower bound to the ground state energy).

Let (x 1 , . . . , x n ) ∈ (R d ) n , let h ′ n be associated via (3.17), and let P νn = i n (x 1 , . . . , x n ) be defined as in (2.19). Then {P νn } n is compact for the weak topology of probability measures on X and any P which is the limit of an extracted subsequence of P νn is admissible and satisfies

lim inf n→∞ n 2/d-2 H n (x 1 , . . . , x n ) -n 2 E[µ 0 ] + n 2 log n 1 d=2 ≥ W(P ). (4.1)
Lower bounds corresponding to (2.22) and (2.23) follow from the above and (2.21).

The proof of Proposition 4.1 will occupy all this section and a large part of the next. More precisely, this section is devoted to the connection between the original problem and the renormalized energy at fixed η, using the general strategy for providing lower bounds for two-scale energies developed in [START_REF] Sandier | From the Ginzburg-Landau Model to Vortex Lattice Problems[END_REF][START_REF] Sandier | 2D Coulomb gases and the renormalized energy[END_REF]. For a lower bound, it is then sufficient to use a lower semi-continuity argment (Fatou's lemma) to pass to the limit η → 0. In this last step we will need the fact that W η is bounded below independently of η. The proof of this is postponed to Section 5.

The local energy

The next step is to study the lim inf when n → 0 of 1 n R d |∇h ′ n,η | 2 and make the two scale structure apparent. This is a general fact, simple application of Fubini's theorem. Let us pick a smooth cut-off function χ with integral 1, support in B(0, 1) and which equals 1 in B(0, 1/2). We use it like a smooth partition of unity, writing

R d |∇h ′ n,η | 2 = R d R d χ(x -y)dx ∇h ′ n,η 2 dy = R d R d χ(x) ∇h ′ n,η (x + y) 2 dydx ≥ n 1/d Σ R d χ(x) ∇h ′ n,η (x + y) 2 dydx = n|Σ|- Σ R d χ(x) ∇h ′ n,η yn 1/d + x 2 dxdy (4.2)
where we simply dropped a part of the integral that we guess will be irrelevant and changed variables. The last line should be interpreted as the average over blow-up centers y of the local electrostatic energy around the center. Indeed, because of the cut-off χ the integral in x is limited to a bounded region, which means we are integrating the square of the original (non blown up) electric field on a region of size n -1/d around y. The last line is in a form that we can treat using the abstract framework of [START_REF] Sandier | 2D Coulomb gases and the renormalized energy[END_REF]Theorem 7]. We will interpret the integral in x as a local energy functional, and deduce from a Γ-convergence and compactness result (with limit depending on y) on this functional a lower bound to (4.2). Note that this is also reminiscent of ideas of Graf-Schenker [GS] when they average over simplices.

The following result shows the coercivity of the "local" energy

χ ∇h ′ n,η yn 1/d + • 2 ,
which is needed to apply the framework of [START_REF] Sandier | From the Ginzburg-Landau Model to Vortex Lattice Problems[END_REF][START_REF] Sandier | 2D Coulomb gases and the renormalized energy[END_REF]. In particular we need to prove compactness of the electric fields and that their limits are in the admissible classes of Definition 2.1.

Lemma 4.2 (Weak compactness of local electric fields).

We pick a sequence of blow-up centers y n → y ∈ R d . Assume that for every R > 1 and for some η ∈ (0, 1), we have where Λ is a discrete set and N p ∈ N * .

sup n B R |∇h ′ n,η (n 1/d y n + •)| 2 ≤ C η,R . (4.3) Then {ν ′ n (n 1/d y n + •)} n is
In addition, there exists

E ∈ L p loc (R d , R d ), p < d d-1 , E η ∈ L 2 loc (R d , R d ), with E η = Φ η (E), such that, up to extraction of a subsequence, ∇h ′ n (n 1/d y n + •) ⇀ E weakly in L p loc (R d , R d ) for p < d d -1 , as n → ∞, (4.4) and ∇h ′ n,η (n 1/d y n + •) ⇀ E η weakly in L 2 loc (R d , R d ) as n → ∞. (4.5)
Moreover E = ∇h, and if y / ∈ ∂Σ, we have

-∆h(x) = c d (ν(x) -µ 0 (y)) (4.6) hence E ∈ A µ 0 (y) .
Proof. Given R > 1, by a mean-value argument, there exists t ∈ (R -1, R) such that for all n,

∂Bt |∇h ′ n,η (n 1/d y n + •)| 2 ≤ C η,R ,
from which it follows that (with a different constant)

∂Bt ∇h ′ n,η (n 1/d y n + •) • ν ≤ C η,R .
But, by Green's formula and (3.17) the left-hand side is equal to

Bt δ (η) x ′ i (n 1/d y n + x ′ ) - Bt µ 0 (y n + n -1/d x ′ ) dx ′ . By (2.2) it follows that, letting ν ′ n := ν ′ n (n 1/d y n + •), we have ν ′ n (B R-1 ) ≤ C d µ 0 L ∞ R d + C η,R .
This establishes that {ν ′ n } is locally bounded independently of n. In view of the form of ν ′ n , its limit can only be of the form ν = p∈Λ N p δ p , where N p are positive integers, and Λ is a discrete set.

Up to a further extraction we also have (4.5) by ( 4.3) and weak compactness in L 2 loc . The compactness and convergence (4.4) follows from Lemma 3.9. The weak local convergences of both ν ′ n and ∇h ′ n (n 1/d y n + •) together with the continuity of µ 0 away from ∂Σ (cf. (2.2)), imply after passing to the limit in (3.17) that E must be a gradient and that (4.6) holds. Finally E η = Φ η (E) because Φ η commutes with the weak convergence in L p loc (R d , R d ) for the ∇h ′ n (n 1/d y n + •) described above. Indeed by definition

Φ η (∇h ′ n (n 1/d y n + •)) = ∇h ′ n,η (n 1/d y n + •) = ∇h ′ n (n 1/d y n + •) + p∈Λn ∇f η (. -p)
where Λ n is the set of points associated with ν ′ n (n 1/d y n + •). Since by assumption all these points have limits, one may check that the sum in the right-hand side converges to p∈Λ ∇f η (. -p), at least weakly in L p loc . Using in addition the convergences (4.5) and (4.4) we deduce

E η = E + p∈Λ ∇f η (. -p),
i.e. E η = Φ η (E) as desired.

4.2 Large n limit: proof of Proposition 4.1

We now have the tools to pass to the double-scale limit at fixed η. The subsequent limit η → 0 will use the following result, whose proof is postponed to Section 5.

Proposition 4.3. W η is bounded below on A 1 by a constant depending only on the dimension.

Proof of Proposition 4.1. As announced, we first fix η > 0 and let n → +∞. We start from (4.2) and apply the framework of Theorem 7 in [START_REF] Sandier | 2D Coulomb gases and the renormalized energy[END_REF]. In the notation of that paper we let G = Σ and X = Σ × L 2 loc (R d , R d ), and take ε = n -1/d . For λ ∈ R d we let θ λ denote both the translation x → x + λ and the action

θ λ E = E • θ λ .
Accordingly the action T n λ is defined for λ ∈ R d by

T n λ (x, E) = x + λn -1/d , E • θ λ .
We then define a functional over X starting from the local electrostatic energy (the yintegral in (4.2))

f n (x, Y) :=    R d χ(y)|Y| 2 (y) dy if Y = ∇h ′ n,η (n 1/d x + •), +∞ otherwise.
(4.7)

Then (4.2) is naturally associated with an average of f n over blow-up centers in Σ:

F n (Y) := - Σ f n (x, θ xn 1/d (Y)) dx. (4.8) Indeed if F n (Y) = +∞, we have F n (Y) ≤ 1 |Σ|n R d |∇h ′ n,η | 2 . (4.9)
Theorem 7 in [START_REF] Sandier | 2D Coulomb gases and the renormalized energy[END_REF] was precisely designed to deduce results at the macroscopic scale (on F n ) from input at the lower scale (on f n ). We now check that its assumptions are satisfied, i.e. that coercivity and Γ-liminf properties follow from

∀R, lim sup n→∞ B R f n (T n λ (x n , Y n )) dλ < +∞, x n ∈ Σ.
(4.10)

But using the definitions above this condition is equivalent to

∀R and ∀n ≥ n 0 , Y n = ∇h ′ n,η (n 1/d x n + •) and lim sup n→∞ χ * 1 B R |Y n | 2 < +∞ (4.11)
where n 0 is a large enough number. This implies that the assumption (4.3) of Lemma 4.2 is satisfied. Up to extraction of a subsequence we may also assume, since x n ∈ Σ which is compact, that x n → x * , and we may apply Lemma 4.2. Thus we have

Y n ⇀ Y * weakly in L 2 loc (R d , R d )
, with all the other results of that lemma, which proves that the compactness assumption of [SS4, Theorem 7] is satisfied. Moreover, this weak convergence implies that lim inf

n→∞ f n (x n , Y n ) ≥ R d χ(y)|Y * | 2 = f (x * , Y * ) (4.12)
by lower semi-continuity, where

f (x * , Y * ) :=        R d χ(y)|Y * | 2 (y) dy if x * ∈ Σ\∂Σ and Y * = Φ η (E) for some E ∈ A µ 0 (x * ) , 0 if x * ∈ ∂Σ +∞ otherwise. (4.13)
This is the Γ-lim inf assumption of [SS4, Theorem 7] and we may thus apply this theorem to pass to the limit n → ∞. Let

P n,η := - Σ δ (x,∇h ′ n,η (n 1/d x+•)) dx
or in other words the push-forward of the normalized Lebesgue measure on Σ by

x → x, ∇h ′ n,η (n 1/d x + •) .
Theorem 7 of [START_REF] Sandier | 2D Coulomb gases and the renormalized energy[END_REF] yields that

• P n,η converges to a Borel probability measure P η on X

• P η is T λ(x) -invariant and its marginal with respect to

x is 1 |Σ| dx |Σ • P η -a.e. (x, E) is of the form lim n→∞ (x n , θ n 1/d xn Y n ),
and moreover

lim inf n→∞ F n (Y n ) ≥ f (x, Y)dP η (x, Y) = lim R→∞ - K R f (x, θ λ Y) dλ dP η (x, Y) (4.14)
where f is defined in (4.13). It is also a part of the result (consequence of the ergodic theorem) that the limit R → ∞ in the above exists. From (4.13), and since ∂Σ is of Lebesgue measure 0 (by assumption (2.1)) we also see that for P η -a.e. (x, Y) it must be that

Y = Φ η (E), E ∈ A µ 0 (x) , and f (x, Y) = R d χ(y)|Y| 2 (y) dy. Combining (4.9) with (4.14) we get lim inf n→∞ 1 n R d |∇h ′ n,η | 2 ≥ |Σ| lim R→∞ - K R χ * 1 K R |Y| 2 dP η (x, Y). (4.15) Since χ * 1 K R ≥ 1 K R-1 we can replace this by lim inf n→∞ 1 n R d |∇h ′ n,η | 2 ≥ |Σ| lim R→∞ - K R |Y| 2 dP η (x, Y). (4.16)
Now the above is true for all η > 0. From the convergence P n,η → P η , and since the bijection Φ η commutes with the convergence in (4.4) (see the proof of Lemma 4.2), we may check that P νn (as defined in (2.19)) converges weakly to P (independent of η), the push-forward of P η by Φ -1 η , for any η. By the above results on P η , and Lemma 4.2, we deduce that P is admissible. Moreover, by definition of the push-forward, (4.16) means that lim inf

n→∞ 1 n R d |∇h ′ n,η | 2 -(κ d w(η) + γ 2 1 d=2 ) ≥ |Σ| lim R→∞ - K R |Φ η (E)| 2 -(κ d w(η) + γ 2 1 d=2 )µ 0 (x) dP (x, E),
where we have used the fact that µ 0 = 1 and the first marginal of P is the normalized Lebesgue measure on Σ. Combining with (3.30) and the definition of W η , we obtain

lim inf n→∞ n 2/d-2 H n (x 1 , . . . , x n ) -n 2 E[µ 0 ] + n 2 log n 1 d=2 ≥ |Σ| c d W η (E) dP (x, E) -Cη 2 .
(4.17)

To conclude, there remains to let η → 0 in the above. Indeed, for P -a.e. (x, E), E ∈ A µ 0 (x) . Since µ 0 is bounded and independent of η, the scaling property (2.13), (2.2) and the uniform lower bound on A 1 of Proposition 4.3 thus imply that W η (E) is bounded below independently of η for P -a.e. (x, E). We can then simply use Fatou's lemma, to deduce that

lim inf n→∞ n 2/d-2 H n (x 1 , . . . , x n ) -n 2 E[µ 0 ] + n 2 log n 1 d=2 ≥ |Σ| c d lim inf η→0 W η (E) dP (x, E),
which is (4.1), by definition of W and W.

Remark 4.4. Note that the positive term 2n n i=1 ζ(x i ) has been discarded in (3.30) and never used in the lower bound, so it can be reintroduced as an additive term in the right-hand side of (4.1), which will be useful when studying the Gibbs measure in Section 7.

Screening and lower bound for the smeared jellium energy

In this section we prove Proposition 4.3, whose statement we recall: Proposition 5.1 (Lower bound on the smeared jellium functional). W η is bounded below on A 1 by a constant depending only on the dimension.

In view of the definition (2.12), the method will consist in bounding from below the minimal energy over a sequence of cubes of size R → ∞. For that we prove

• that by minimality we can reduce to configurations with simple and well-separated points (at least away from the boundary of the cube), i.e. at distances bounded below independently of η and R. This is done in Section 5.1 via an argument adapted from a Theorem of E. Lieb [Lie2] for the case η = 0 (whose statement and proof can be found in [RNS]).

• that we can "screen" efficiently such a configuration, i.e. modify it close to the boundary of the cubes to make the normal component of the electric field vanish on the boundary of the cube and make the points well-separated all the way to the boundary of the cube, at a negligible energy cost. (The vanishing normal component will in particular impose the total number of points in the cube). This is carried out in Section 5.2.

• Once this is done, we will be able to immediately bound from below the minimal energy on these large cubes via Lemma 3.2.

In all this section K R is a hyperrectangle whose side-lengths are in [2R, 3R]. For the actual proof of Proposition 4.3 we need consider only hypercubes but we will need more general shapes later in the paper.

Separation of points

Let us consider the following variational problem:

F η,R := inf K R |∇h| 2 , -∆h = c d p∈Λ N p δ (η) p -1 in K R , for some discrete set Λ ⊂ K R .
(5.1) Obviously F η,R ≥ 0 and the infimum exists. We do not address the question of whether there exists a minimizer. Note that the points p may in principle depend on η.

Proposition 5.2 (Points "minimizing" F η,R are well-separated). Let Λ be a discrete subset of K R and h satisfy

-∆h = c d p∈Λ N p δ (η) p -1 := c d (µ h -1) in K R .
(5.2)

Denote Λ R = Λ ∩ K R-1 .
There exists three positive constants η 0 , r 0 , C such that if η < η 0 , R is large enough and one of the following conditions does not hold: (5.4) then there exists Λ a discrete subset of K R and an associated potential h satisfying

∀p ∈ Λ R , N p = 1 (5.3) ∀p ∈ Λ R , dist(p, Λ R \ {p}) ≥ r 0 ,
-∆ h = c d p∈ Λ N p δ (η) p -1 in K R (5.5) such that K R |∇ h| 2 ≤ K R |∇h| 2 -C.
What this proposition says is that if two points in the configuration are too close to one another (independently of R and η), one can always find another configuration with strictly less energy. In other words, in a minimizing sequence for F η,R one can always assume that (5.3)-( 5.4) hold for some r 0 > 0 (depending only on d). Note that the failure of (5.3) may be seen as an extreme case of the failure of (5.4). Without loss of generality we may thus prove the result when (5.4) fails and N p = 1 for all p.

We start with a lemma that gives explicitly the variation of K R |∇h| 2 when one point in the configuration is moved and the variation of the potential is chosen to satisfy Dirichlet boundary conditions. We shall denote G R the constant c d times the Green-Dirichlet function of the hyperrectangle K R :

-∆ x G R (x, y) = c d δ y (x) if x ∈ K R and y ∈ K R G R (x, y) = 0 if x ∈ ∂K R or y ∈ ∂K R .
(5.6)

Note that G R is a symmetric function of x and y: G R (x, y) = G R (y, x). Associated to the Green-Dirichlet function is the quadratic form

D R (µ, ν) = K R ×K R µ(dx)G R (x, y)µ(dy) (5.7)
defined for measures µ, ν, that plays a role analogue to (2.3). We will consider specific variations of F η,R that can be expressed in terms of G R :

Lemma 5.3 (Variations of F η,R ).
Let h, Λ and µ h be as in Proposition 5.2. Let x ∈ Λ and y ∈ K R and define

Λ = (Λ \ {x}) ∪ {y}
with an associated h = h + h where the variation h is defined as the unique solution to

-∆ h = c d (δ (η) y -δ (η) x ) in K R h = 0 on ∂K R .
(5.8)

We have

K R |∇ h| 2 = K R |∇h| 2 + c d D R δ (η) y , δ (η) y -c d D R δ (η) x , δ (η) x + 2c d D R µ h -δ (η) x -1, δ (η) y -δ (η)
x .

(5.9)

Proof. Let us first note that h is an admissible trial state for F η,R since we immediately check that it satisfies (5.5). Expanding the square, integrating by parts the cross-term and using the Dirichlet boundary condition for h in addition to Equation (5.2) we obtain

K R |∇ h| 2 = K R |∇h| 2 + K R |∇ h| 2 + 2 K R ∇ h • ∇h = K R |∇h| 2 + K R |∇ h| 2 -2 K R h∆h = K R |∇h| 2 + K R |∇ h| 2 + 2c d K R h(µ h -1).
We now use the Green representation of h to write

h(z) = K R G R (z, z ′ )δ (η) y (z ′ )dz ′ - K R G R (z, z ′ )δ (η)
x (z ′ )dz ′ := hy (z) -hx (z).

We then have (using the Dirichlet boundary condition again to integrate the cross term by parts)

K R |∇ h| 2 = K R |∇ hx | 2 + K R |∇ hy | 2 -2 K R ∇ hx • ∇ hy = c d D R (δ (η) x , δ (η) x ) + c d D R (δ (η) y , δ (η) y ) -2c d D R (δ (η) x , δ (η) y )
and

K R h(µ h -1) = D R (δ (η) y -δ (η) x , µ h -1) = D R (δ (η) y -δ (η) x , µ h -δ (η) x -1) + D R (δ (η) y -δ (η) x , δ (η) x ).
Putting everything together yields the desired result.

The next step is to remark that for large R we should have G R (x, y) ≃ w(x -y). It is then natural to expect that the difference between the self-interactions

D R δ (η) y , δ (η) y -D R δ (η) x , δ (η)
x should be small, at least if the points x and y do not approach the boundary of the domain and thus do not see the Dirichlet boundary condition. On the other hand we can rewrite

D R (δ (η) y -δ (η) x , µ h -δ (η) x -1) = δ (η) y h x -δ (η) x h x ≃ h x (y) -h x (x) (5.10)
where

h x := K R G R (., y) µ h -δ (η)
x -1 (y)dy.

(5.11)

Were we dealing with the case of point charges η = 0, the variation would thus reduce to h x (y) -h x (x). Lemma 5.3 is then a rephrasing in our context where charges are smeared out of the well-known fact that the condition for optimality is that the charge at x lies at the minimum of the potential h x generated by all the other charges and the neutralizing background.

We thus carry on with estimates of the Green-Dirichlet function that will allow us to make this intuition rigorous in the case of smeared out charges. We will use the well-known fact that G R (x, y) may be written as the sum of the Green function of the whole space w(x -y) and a regular part. Important for us is the fact that the regular part is uniformly bounded on K R with a bound independent of R, such that we will be able to focus on the singular part w(x -y) when estimating the variation in the right-hand side of (5.9).

Actually, when d = 2, things are a little bit more subtle due to the fact that the Green function does not decay at infinity. Our estimate (5.13) is thus different in this case, but still sufficient for our proof.

Lemma 5.4 (Estimates on the Green function).

• If d ≥ 3, there exists a constant C G , independent of R, such that |G R (x, y) -w(x -y)| ≤ C G (5.12) for all x, y ∈ K R with min(dist(x, ∂K R ), dist(y, ∂K R )) ≥ 1. • If d = 2, there exists a constant C G , independent of R, such that G R (x, y) -w(x -y) + w(x -y ′ ) ≤ C G (5.13)
where y ′ is the reflection of y with respect to ∂K R . In particular, if dist(y, ∂K R ) ≥ 1 and x, x ′ ∈ B(y, c) for some constant c, then

G R (x, y) -G R (x ′ , y) -w(x -y) -w(x ′ -y) ≤ C c (5.14)
for some C c depending only on c.

Proof. The proof of (5.12) is a well-known argument: Since G R (x, y) = G R (y, x) we may restrict to the case where dist(y, ∂K R ) ≥ 1 and consider G R (x, y) and w(x -y) as functions of x only. Now, by definition G R (., y)-w(.-y) is harmonic in K R and hence (by the maximum principle) reaches its maximum and minimum on the boundary of K R . The Green-Dirichlet function is zero there and in view of the assumption that dist(y, ∂K R ) ≥ 1 we have

min x∈∂K R (-w(x -y)) ≥ -C independently of R and max x∈∂K R (-w(x -y)) ≤ -w(CR) → 0
when R → ∞. We used the asumption d ≥ 3 in the last equation.

To prove (5.13), we fix y ∈ K R and note that the function

f (x, y) := G R (x, y) -w(x -y) + w(x -y ′ ) is harmonic as a function of x ∈ K R since ∆ x G R (x, y) = ∆ x w(x-y), ∆ x w(x-y ′ ) = -c d δ y ′ (x)
and y ′ / ∈ K R . We thus know that f (x, y) reaches its maximum and minimum on the boundary ∂K R :

inf x∈∂K R log |x -y| |x -y ′ | ≤ f (x, y) ≤ sup x∈∂K R log |x -y| |x -y ′ | , (5.15)
where we used the Dirichlet boundary condition for G R (., y) and the definition of w when d = 2. Next we observe that, y ′ being the reflection of y with respect to ∂K R it holds that

|x -y| ≤ |x -y ′ |
for any x ∈ ∂K R and thus f (x, y) ≤ 0 for any x, y ∈ K R . For a lower bound we note that |x -y ′ | ≤ |x -y| + |y -y ′ | with equality when x, y and y ′ are aligned, in which case

|x -y| |x -y ′ | = 1 - |y -y ′ | |x -y ′ | = 1 - |y -y ′ | L + |y -y ′ |/2
where L is one of the side-lengths of the hyperrectangle K R . By assumption 2R ≤ L ≤ 3R, which implies that |y -y ′ | ≤ 3R and we easily conclude that inf x∈∂K R log |x-y| |x-y ′ | is bounded below independently of R. Plugging in (5.15) concludes the proof of (5.13). The estimate (5.14) follows immediately by noting that, if dist(y

, ∂K R ) ≥ 1 then |∇w(. -y ′ )| is bounded above independently of R in K R .
We can now give the main idea of the proof of Proposition 5.2, following [Lie2] for singular charges. Suppose there are two points, say 0 and x in Λ, very close to one another6 . Then one can decrease the energy by sending x to some point y with dist(0, y) > dist(0, x) and changing the potential as described in Lemma 5.3. As we will prove below, it is sufficient to consider only the contribution of the charge at 0 and of the background at some fixed distance thereof to bound from above the variation of the potential h x defined in (5.10). Then, if the distance dist(0, x) is really small to begin with, Lemma 5.4 tells us that we can approximate the variation (5.10) in the potential h x using the Coulomb potential w(x -y). This observation leads us to the definition

V η = R d w(. -z) δ (η) 0 -1 B(0,r 3 ) (z)dz (5.16)
of the potential generated (via the Green function w of the full space) by a unit charge at 0 smeared out at scale η and a disc of constant opposite density of charge, where we choose r 3 < 1 so that |B(0, r 3 )| < 1. We now observe that V η is radial and decreasing. It is also obviously independent of R and gets closer and closer to V 0 , independent of η, when η is small. Since V 0 (r) → +∞ when r → 0, one can make V η (y) -V η (x) arbitrarily negative if x → 0 while y stays bounded at a fixed distance away from 0. This is the content of the following Lemma 5.5 (The jellium potential close to a point charge). For any M > 0 there exists min( 1 2 , r 3 ) > r 2 > r 1 > 0 such that for any small enough η and any r < r 1 , we have

V η (r 2 -η) -V η (r + η) ≤ -M.
(5.17)

Proof. By Newton's theorem, we have V η (r) = V 0 (r) for r ≥ η, thus since V η is decreasing, we have V η (r) ≥ V 0 (η) for r < η. Thus, since we will always choose η < r 1 < r 2 , we can study only V 0 , which is readily computed by taking advantage of the radiality of the charge distribution:

V 0 (r) = (1 -|B(0, r)|) w(r) + |S d-1 | B(0,r 3 )\B(0,r) w(t)t d-1 dt if 0 ≤ r ≤ r 3 V 0 (r) = (1 -|B(0, r 3 )|) w(r) if r ≥ r 3 ,
by Newton's theorem again. Simple considerations show that V 0 (r) → +∞ when r → 0 and that V 0 is decreasing. The result follows.

We can now proceed to the Proof of Proposition 5.2.

As announced, we work assuming N p = 1 and (5.4) is violated. Without loss of generality (we may always translate K R without changing the energy) we assume that 0 ∈ Λ R and that there exists some point x ∈ Λ R close to 0. We will prove that if dist(0, x) < r 1 , where r 1 is as in Lemma 5.5 one can decrease the energy by moving x to a certain point y ∈ ∂B(0, r 2 ).

Since we assume dist(0, ∂K R ) ≥ 1, and since r 1 < 1 2 , we have B(0, r 1 ) ⊂ K R-1/2 . We pick some y, to be chosen later, and let Λ be the corresponding modified configuration, with an associated potential h as in Lemma 5.3. The variation of the energy is given by (5.9), which we estimate using (5.12) (respectively (5.14) when d = 2): for the first two terms we may write

D R δ (η) y , δ (η) y -D R δ (η) x , δ (η) x ≤ C δ (η) x + δ (η) y ≤ C
since the contribution of the singular part w(x -y) to the Green function G R (x, y) yields the same contribution for both terms. We thus have

K R |∇ h| 2 - K R |∇h| 2 ≤ 2c d D R µ h -δ (η) x -1, δ (η) y -δ (η) x + C (5.18)
where C is independent of R and η. We may then focus on proving that the first term of the right-hand side can be made arbitrarily negative. For that purpose, we start from (5.10), (5.11) and, following [Lie2], we decompose

h x = V η,R + h rem V η,R = K R G R (., z) δ (η) 0 -1 B(0,r 3 ) (z)dz h rem = K R G R (., z)   p∈Λ\{x,0}) δ (η) p -1 + 1 B(0,r 3 )   (z)dz. (5.19)
We further define the η-averaged function h rem a (z) := h rem * δ (η) 0 , and note that since -∆h rem ≥ 0 in B(0, r 3 ), and r 2 < r 3 , we have -∆h rem a ≥ 0 in B(0, r 2 ), at least if η is small enough. By the maximum principle h rem a reaches its minimum on B(0, r 2 ) at the boundary and we now choose y ∈ ∂B(0, r 2 ) to be (one of) its minimum point(s). Then, in view of (5.10), we have

D R µ h -δ (η) x -1, δ (η) y -δ (η) x = h rem a (y) -h rem a (x) + δ (η) y V η,R -δ (η) x V η,R ≤ δ (η) y V η,R -δ (η) x V η,R ≤ δ (η) y V η -δ (η) x V η + C ′ (5.20)
where on the second line we have used (5.12) (respectively (5.14) when d = 2) by replacing G R (x, y) by w(x -y) in the terms involving V η,R . Assuming now that we have chosen M = 2(C + c d C ′ ) in Lemma 5.5, where C is the constant in (5.18) and C ′ that in (5.20), we have

V η (r 2 -η) -V η (r 1 + η) ≤ -2C -2C ′ .
Using that V η is radial and decreasing, we deduce that, since we assumed |x| ≤ r 1 and |y| = r 2 , we have

δ (η) y V η -δ (η) x V η ≤ -2C -2C ′ .
Combining this with (5.20) and (5.18), we have thus obtained

K R |∇ h| 2 -K R |∇h| 2 ≤ -C
for some C > 0. This proves the proposition, with r 0 = r 1 .

Screening

Starting with a configuration of points modified according to Proposition 5.2 we are able to modify it further and "screen it" as announced.

Proposition 5.6 (Screening a configuration of points).

There exists η 0 > 0 such that the following holds for all η < η 0 . Let K R be a hyperrectangle such that |K R | is an integer and the sidelengths of K R are in [2R, 3R]. Let E η = ∇h η be a gradient vector field with h η satisfying (5.2) with (5.3)-(5.4) satisfied and

K R |E η | 2 ≤ CR d .
(5.21)

Then there exists Λ a configuration of points and Ē an associated gradient vector field (both possibly also depending on η) defined in K R and satisfying

-div Ē = c d p∈ Λ δ p -1 in K R (5.22) Ē • ν = 0 on ∂K R (5.23) such that for any p ∈ Λ min dist(p, Λ \ {p}), dist(p, ∂K R ) ≥ r 0 10 (5.24)
with r 0 as in (5.4), and

K R | Ēη | 2 ≤ K R |E η | 2 + o(R d ) (5.25)
as R → ∞, where Ēη = Φ η Ē and the o(R d ) depends only on η and the constants in (5.3), (5.4), (5.21).

Property (5.23) is crucial in order to be able to paste together configurations defined in separate hyperrectangles. It also implies (integrating (5.22) over K R and using Green's formula) that the number of points of the modified configuration is exactly equal to the volume of the domain:

# Λ = |K R |, (5.26)
which is also important for the proof of Proposition 4.3.

The proof proceeds by modifying the configuration of points in the vicinity of ∂K R :

1. We first select by a mean-value argument a "good boundary" ∂K t on which we control the L 2 norm of E η (i.e. the energy density), and at a distance 1 ≪ L ≪ R from the boundary of the original hyperrectangle.

2. We will not move the points whose associated smeared charges intersect ∂K t . Instead, we isolate them in small cubes and leave unchanged all the points lying in Γ= the union of K t with these small cubes. A new configuration of points is built in K R \Γ. It is built by splitting this set into hyperrectangles and using on each some test vectorfields (obtained by solving explicit elliptic PDEs) whose energies is evaluated by elliptic estimates similarly to [START_REF] Sandier | From the Ginzburg-Landau Model to Vortex Lattice Problems[END_REF]Sec. 4.3]. The test vector-fields are then pasted together.

3. It is more convenient to first "straighten" the boundary of Γ. We pick some α such that Γ ⊂ K t+α and put points in K t+α \ Γ so that neither the energy nor the flux at the boundary are increased too much. To obtain a control on those quantities, it is important that no smeared charge intersects the boundary of Γ to begin with, whence the second step.

We start with the preliminary lemmas containing the elliptic estimates which provide the elementary bricks of our construction.

Lemma 5.7 (Adding a point without flux creation). Let K be a hyperrectangle of center 0 and sidelengths in [l 1 , l 2 ], with l 2 ≥ l 1 > 0, and let m = 1/|K|. The mean zero solution to

-∆u = c d (δ 0 -m) in K ∇u • ν = 0 on ∂K satisfies lim η→0 K |∇u η | 2 -κ d w(η) ≤ C
where C depends only on d, l 1 , l 2 and m.

Proof. It is a very simple computation, similar for example to arguments employed in the proof of Lemma 3.2.

Lemma 5.8 (Correcting fluxes on hyperrectangles).

Let K be a hyperrectangle with sidelengths in [l/2, 3l/2]. Let g ∈ L 2 (∂K) be a function which is 0 except on one face of K. Let m be defined by

m = 1 -c -1 d |K| -1 ∂K g (5.27)
and assume |m -1| < 1/2. Then the mean zero solution to (5.29) where C depends only on d.

-∆u = c d (m -1) in K ∇u • ν = g on ∂K (5.28) satisfies K |∇u| 2 ≤ Cl ∂K g 2 ,
Proof. We may write u = u 1 + u 2 where

-∆u 1 = c d (m -1) in K ∇u 1 • ν = ḡ on ∂K
with ḡ a constant on the face of ∂K where g is nonzero, and zero otherwise (integrating the equation we find ḡ = c d (m -1)l) and

-∆u 2 = 0 in K ∇u 2 • ν = g -ḡ on ∂K.
First we note that by (5.27) and the Cauchy-Schwarz inequality we have

c d |m -1| ≤ 1 |K| ∂K |g| ≤ Cl -d-1 2 ∂K g 2 1/2 .
(5.30)

Next we observe that u 1 can be solved explicitly: if the face where g is nonzero is included in {x 1 = 0} (after translation and rotation), then u 1 (x) = m-1 2 (x 1 -l) 2 and we may compute

K |∇u 1 | 2 ≤ C(m -1) 2 l d-1 l 3 = C(m -1) 2 l d+2 .
(5.31)

Secondly, by scaling we check that

K |∇u 2 | 2 ≤ Cl ∇u 2 • ν 2 L 2 (∂K) ≤ Cl g 2 L 2 (∂K) + (m -1) 2 l 2 l d-1
(5.32)

Combining (5.31), (5.32), and inserting (5.30), we obtain the result.

Proof of Proposition 5.6. We set two lengths l and L such that 1 ≪ l ≪ L ≪ R as R → ∞, and to be determined later. In all the sequel we denote by C 1 a generic constant which is a multiple (by a universal factor) of the constant C of (5.21).

-Step 1: We find a good boundary. First we claim that there exists

t ∈ [R -2L, R -L] such that ∂Kt |E η | 2 ≤ C 1 R d /L (5.33)
and

K t+1 \K t-1 |E η | 2 ≤ C 1 R d /L.
(5.34) Indeed, since (5.21) holds, we may split K R-L \K R-2L into L/2 "annular" regions of width 2. On one of them (5.34) must hold, and by a mean-value argument we may also shift t so that (5.33) holds.

Since we assume (5.3)-(5.4) and since L ≥ 1, the boundary ∂K t intersects balls B(p, η) with p ∈ Λ which are all at distance ≥ 1 2 r 0 > 0 from all other such balls (for that it suffices to take η 0 < r 0 /4). Without l.o.g. we assume r 0 < 1. Let Λ 0 denote the set of p's for which B(p, η) ∩ ∂K t = ∅. For each p ∈ Λ 0 , let K s (p) denote the hypercube of sidelength s centered at p. By a mean-value argument similar to the above, and using (5.34) we can find s ∈ [r 0 /16, r 0 /8] such that

p∈Λ 0 ∂Ks(p) |E η | 2 ≤ C 1 R d /L. (5.35) We now denote Γ = K t ∪ (∪ p∈Λ 0 K s (p))
. By construction ∂Γ does not intersect any B(p, η) and we have

∂Γ |E η | 2 ≤ C 1 R d /L.
(5.36)

By perturbing this construction a little we may define the resulting set Γ such that |Γ| ∈ N and (5.36) still holds.

-Step 2: definition of E in K t+α . First we define E to be Φ -1 η (E η ) in Γ, that is we set

E = E η - p∈Λ∩Γ f η (x -p) in Γ
where Λ is the set of points corresponding to E η as in (5.2) and f η is defined as in (2.8). We note that this definition is unambiguous and coincides with the restriction to Γ of Φ -1 η (E η ) since ∂Γ does not intersect any B(p, η).

We next wish to extend E to K t+α \Γ. We take α > 1 to be such that |K t+α \Γ| is an integer. Since |Γ| ∈ N by construction and |K R | ∈ N by construction we thus have

|K t+α | ∈ N and |K R \ K t+α | ∈ N.
(5.37)

Then K t+α \Γ can be split into a disjoint union of hyperrectangles of volume 1, whose sidelengths are all bounded below by r 0 /8, and bounded above by 1, i.e. by universal constants.

In each of these hyperrectangles we apply Lemma 5.7 with m = 1. Since the normal derivatives of the functions obtained that way are zero on the boundary, their gradients can be glued together into a global vector field, and no divergence will be created at the boundaries. More precisely this means that there exists a vector field X defined in K t+α \Γ and satisfying

-div X = c d p∈Λ 1 δ p -1 in K t+α \Γ
for some set Λ 1 equal to the union of the centers of these hyperrectangles, and

K t+α \Γ |X η | 2 ≤ C η R d-1 (5.38)
since the energy cost given by Lemma 5.7 is proportional to the volume concerned. Next we claim that we can find a vector field Y satisfying div Y = 0 in K t+α \Γ and g defined over (5.39) and (5.40) To see this, we may split K t+α \Γ into hyperrectangles K i which all have one face included in ∂Γ and one face included in ∂K t+α , and whose sidelengths are all bounded above by 2 and below by r 0 /8. Then we solve in each

∂K t+α such that Y • ν = g on ∂K t+α , while Y • ν = E η • ν on ∂Γ and ∂K t+α g 2 ≤ C ∂Γ |E η | 2 ≤ C 1 R d /L,
K t+α \Γ |Y | 2 ≤ C 1 R d /L.
K i        ∆u i = 0 in K i ∇u i • ν = E η • ν on ∂K i ∩ ∂Γ ∇u i • ν = g i on ∂K i ∩ ∂K t+α ∇u i • ν = 0 else on K i
where g i is the unique constant that makes the equation solvable. It is straightforward to check that ∂K i ∩∂K t+α g 2 i ≤ C ∂K i ∩∂Γ |E η | 2 (again by Green's formula applied to the equation). Defining g on ∂K t+α as g = g i on ∂K t+α ∩ ∂K i , we then have (5.39). Since u i is harmonic, we also have the estimate

K i |∇u i | 2 ≤ Cl ∇u i • ν 2 L 2 (∂K i ) ≤ C ∂K i ∩∂Γ |E η | 2 .
Defining now the vector field Y to be ∇u i on each K i , we see that Y satisfies the desired properties and that (5.40) holds. We finally set E = X + Y in K t+α \Γ. It satisfies -div E = c d p∈Λ 1 δ p -1 and, since no divergence is created at the interface ∂Γ, the vector field E now defined in the whole on K t+α still satisfies

-div E = c d p∈Λ 1 ∪(Γ∩Λ) δ p -1 in K t+α E • ν = g on ∂K t+α .
(5.41)

Moreover, in view of (5.38) and (5.40), we have

K t+α | E η | 2 ≤ Γ |E η | 2 + C η R d-1 + CR d /L.
(5.42)

Note that by construction the distances between the points in Λ 1 ∪ (Γ ∩ Λ) are all bounded below by r 0 /8. We now discard the notations used for this step, except for the conclusions (5.41)-(5.42) and (5.39).

-Step 3: correcting the flux. We construct a new domain 7 Γ ′ on the boundary of which we make the flux vanish. This uses Lemma 5.8 in an essential way to construct a configuration in Γ ′ \ K t+α . First we split ∂K t+α into O((R/l) d-1 ) hyperrectangles (of dimension d-1) I i of sidelength ∈ [l/2, 3l/2]. For each of them we consider a hyperrectangle (of dimension d) included in K R \K t+α which has one side equal to I i . By perturbing the sizes of the sides, we may have a hyperrectangle with aspect ratios in [1/2, 3/2] (so all sides have sizes in [l/2, 3l/2]). This forms a disjoint collection K i (we use the same notation as in the previous step, even though it does not correspond to the same rectangles). We let g i be the restriction of g on ∂K t+α to I i , and let m i be defined by

c d (m i -1)|K i | = - ∂K i g i .
(5.43) Let us check that |m i -1| < 1 2 . Using the Cauchy-Schwarz inequality and (5.39) we have

c d |m i -1| ≤ l -d-1 2 ∂K i g 2 i 1/2 ≤ l -d-1 2 ∂K t+α g 2 1/2 ≤ Cl -1/2-d/2 R d/2 L 1/2 .
(5.44)

It is clear that we may choose (5.45) such that this is o(1) as R → ∞, which we do from now on. Since m i ∼ 1 and

1 ≪ l ≪ L ≪ R
|K i | ∼ Cl d ≫ 1
it is also clear that we may perturb once more the sizes of the sides of K i to ensure in addition that m i |K i | ∈ N, with the previous conditions preserved. We may then apply Lemma 5.8 over

K i , this yields a u i satisfying      -∆u i = c d (m i -1) in K i ∇u i • ν = g i on I i ∇u i • ν = 0 on ∂K i \ I i and K i |∇u i | 2 ≤ Cl I i g 2 ,
where C depends only on d. We next split each K i into hyperrectangles of sides ∈ [1/2, 3/2] and of volume m -1 i , which is possible since by construction m i |K i | ∈ N. On each of these hyperrectangles we apply Lemma 5.7. Pasting together the ∇u's given by that lemma yields an X i defined over K i , such that -div

X i = c d ( k δ x k -m i ) in K i X i • ν = 0 on ∂K i . Also, denoting (X i ) η = Φ η (X i ), it follows from the fact that |m i -1| < 1 2 that K i |(X i ) η | 2 ≤ C η l d .
(5.46)

We note that by construction, all of these x k 's are at distance at least 1/2 of ∂K i hence, since η ≤ 1/2, the balls B(x k , η) do not intersect ∂K i . Defining now Y i = X i + ∇u i in K i we have obtained a solution of

-div Y i = c d ( k δ x k -1) in K i Y i • ν = g i on ∂K i and K i |(Y i ) η | 2 ≤ 2 K i |(∇f i ) η | 2 + 2 K i |∇u i | 2 ≤ C η l d + Cl ∂K i g 2 i .
(5.47)

The new domain Γ ′ is K t+α ∪ ∪ i K i and we may paste the

Y i constructed above to obtain a Y satisfying    -div Y = c d ( k δ x k -1) in Γ ′ Y i • ν = g on ∂K t+α Y i • ν = 0 on ∂Γ ′ ∂K t+α -
Step 4: completing the construction. There remains to complete the construction in

K R \Γ ′ = K R \(K t+α ∪ ∪ i K i ).
To this end, note that (5.43) implies

c d i m i |K i | -c d i |K i | =
|K R \(K t+α ∪ ∪ i K i )| ∈ N.
We may thus tile K R \(K t+α ∪ ∪ i K i ) by hyperrecangles volume 1 ands sidelengths ∈ [1/2, 3/2] on which we apply Lemma 5.7 with m = 1. We paste together the gradients of the functions obtained this way. It gives a contribution to the energy proportional to the volume i.e. of CR d-1 L. We also paste this with the Y i 's of the preceding step. It gives a final vector field

E defined in K R \K t+α , satisfying      -div E = c d ( k δ x k -1) in K R \K t+α E • ν = g on ∂K t+α E • ν = 0 on ∂K R .
Here all the points x k are at distance ≥ 1 2 > η from ∂K t+α . The energy of E can now be controlled as follows:

K R \K t+α | E η | 2 ≤ CR d-1 L + C η R d-1 l 1-d l d + l i ∂K i g 2 i = C η R d-1 L + l ∂K t+α g 2 ≤ C η (R d-1 L + lR d /L). (5.48)
Combining with (5.41), we have thus obtained an E on all K R satisfying (5.22)-(5.23). Using (5.42), (5.48) and the conditions l ≪ L ≪ R we have

K R | E η | 2 ≤ K R |E η | 2 + C η (R d-1 L + lR d /L) ≤ K R |E η | 2 + o(R d ),
-Step 5: making E a gradient. E satisfies (5.22) but is not necessarily a gradient. We claim that we may add to it a vector-field X to make it a gradient, while still having (5.22)-(5.23), and decreasing the energy. It is standard by Hodge decomposition that we may find X satisfying div X = 0 in K R and X • ν = 0 on ∂K R , and such that Ē := E + X is a gradient. Then Ē still satisfies (5.22)-(5.23), and we note that

K R | E η | 2 = K R | Ēη | 2 + K R |X | 2 -2 K R Ēη • X .
Since Ē is a gradient and thus Ēη too, and since div X = 0 and X • ν = 0 on ∂K R , we find that the last integral vanishes, and so

K R | Ēη | 2 ≤ K R | E η | 2 .
This last operation has thus not increased the total energy and the vector-field Ē satisfies all the desired conclusions.

An immediate consequence of the screening process is the following Proposition 5.9 (Periodic minimizing sequences).

Given η > 0 small enough, for any R large enough and any hyperrectangle K R of sidelengths in [2R, 3R] such that |K R | ∈ N, there exists an Ē satisfying (5.22), (5.23), (5.24) and lim sup

R→∞ - K R | Ēη | 2 -(κ d w(η) + γ 2 1 d=2 ) ≤ inf A 1 W η (5.49)
Also, min A 1 W η admits a minimizing sequence made of periodic vector fields, and so does min A 1 W.

The construction of periodic minimizing sequences proves the corresponding claim in Theorem 1. The existence of a minimizer is a direct consequence of Theorem 2 as discussed in Section 6. We also give a self-contained proof in Appendix A.

Proof. It suffices to take K R = [-R, R] d and ∇h η which approximates F η,R and apply Propositions 5.2 and 5.6. This yields the desired Ē satisfying (5.22), (5.23), (5.24), (5.26) and lim sup

R→∞ - K R | Ēη | 2 -(κ d w(η) + γ 2 1 d=2 ) ≤ lim sup R→∞ F η,R |K R | -(κ d w(η) + γ 2 1 d=2 )
and then (5.49) is an obvious consequence of the definitions of W η and F η,R . In addition, Ē is the gradient of a function h R over K R , with ∇h R • ν = 0 on ∂K R . We may then reflect h R with respect to ∂K R to make it into a function over [-R, 3R] d which can then be extended periodically to the whole R d . Its gradient then defines a periodic vector field Y R over R d , which belongs to A 1 . By periodicity of Y R and definition of W η , we obviously have

W η (Y R ) = -K R | Ēη | 2 -(κ d w(η) + γ 2 1 d=2
). In view of (5.49), this implies that we have lim sup R→∞ W η (Y R ) ≤ inf A 1 W η . We thus conclude that inf A 1 W η admits a minimizing sequence made of periodic vector fields. Using a diagonal argument to deal with the η → 0 limit, inf A 1 W also does.

Proof of Proposition 5.1. In view of Proposition 5.9, to bound from below W η , it suffices to bound from below -K R | Ēη | 2 with Ē satisfying (5.22), (5.23), (5.24) and (5.26). But then, we are in a situation where Lemmas 3.2 applies. With (5.26), the combination of these results easily yields lim sup

R→∞ - K R | Ēη | 2 -(κ d w(η) + γ 2 1 d=2 ) ≥ -C
where C depends only on d and on r 0 , which itself depends also only on d. In view of (5.49) this concludes the proof.

6 Upper bound to the ground state energy

In this section we complete the proof of Theorem 2 by proving the following Proposition 6.1 (Energy upper bound).

For any ε > 0 there exists r 1 > 0 and for any n a set A n ⊂ (R d ) n such that

|A n | ≥ n! π(r 1 ) d /n n (6.1)
and for any (y 1 , . . . , y n ) ∈ A n we have

lim sup n→∞ n 2/d-2 H n (y 1 , . . . , y n ) -n 2 E[µ 0 ] ≤ ξ d + ε if d ≥ 3 (6.2) lim sup n→∞ n -1 H n (y 1 , . . . , y n ) -n 2 E[µ 0 ] + n 2 log n ≤ ξ 2 + ε if d = 2. (6.3)
To deduce Theorem 2 it suffices to let ε → 0 in (6.2) (respectively (6.3)). This yields the upper bound corresponding to (2.22) (respectively (2.23)). On the other hand, we have a lower bound for the ground state energy, given by Proposition 4.1. Comparing the two yields the result of Theorem 2, where it also comes as a by-product that P minimizes W, which can only happen if E minimizes W over A µ 0 (x) for P -a.e. (x, E) ∈ X. Note that this also proves the existence result in Theorem 1: for P obtained in Theorem 2, we have P -a.e. (x, E) minimizes W in A µ 0 (x) , hence this implies the existence of a minimizer for W in some A m , hence in all A m for all m by (2.13). The difficult part is of course the existence of the probability measure P , which employs the ergodic framework of [START_REF] Sandier | 2D Coulomb gases and the renormalized energy[END_REF]Theorem 7]. In Appendix A we provide a more direct proof where this method is applied directly to the renormalized energy. Remark 6.2. In view of Remark 4.4, comparing upper and lower bounds also yields that for minimizers (or almost minimizers) of

H n we have n i=1 ζ(x i ) = o(n 1-2/d ) as n → ∞.
Since ζ is expected to typically grow quadratically away from Σ, this provides a control on how far from Σ the points can be. In fact, arguing as in [RNS] one can show that for minimizers, there are no points outside Σ.

The fact that the upper bound holds, up to a small error ε, for any configuration in A n , a set which has a reasonably large volume in configuration space, will be crucial when studying the fluctuations of the Gibbs measure at finite temperature. It would be very interesting to know more precisely the behavior of W around a minimizing configuration.

The proof relies on an explicit construction, using the test charge configurations and electric fields of Corollary 5.9. Roughly speaking, we fix η > 0 small and split Σ into hyperrectangles of sidelengths of order Rn -1/d for some R that will ultimately tend to ∞. Equivalently By standard elliptic regularity we have

∇u K L 2 (K) ≤ C R m K -µ ′ 0 L ∞ (K) ≤ C R n -1/d . (6.8)
Indeed, µ 0 is assumed to be C 1 on its support (this is where we use this part of (2.2)) hence after scaling ∇µ ′ 0 L ∞ (Σ ′ ) ≤ Cn -1/d . Next we denote σ m the rescaling to scale m:

σ m E = m 1-1/d E(m 1/d .)
and define E K to be

E K = ∇u K + σ m K E if x ∈ K (6.9)
where E is provided by Proposition 5.9 applied to m

1/d K η over the hyperrectangle m 1/d K K (and suitably translated), hence satisfies m 1/d K K | E m 1/d K η | 2 -|m 1/d K K|(κ d w(m 1/d K η) + γ 2 1 d=2 ) ≤ |m 1/d K K|(min A 1 W m 1/d K η + o R (1)
). (6.10)

The vector field E K then satisfies

K | (E K ) η | 2 ≤ K σ m K E η 2 + ∇u K 2 L 2 (K) + 2 ∇u K L 2 (K) (σ m K E) η L 2 (K) . (6.11)
But by change of scales, we have

K σ m K E η 2 = m 1-2/d K m 1/d K K E m 1/d K η
2 so (6.10) gives that this is bounded above by

m 1-2/d K m K |K|(κ d w(ηm 1/d K ) + γ 2 1 d=2 ) + m K |K| min A 1 W m 1/d K η + o R (1)) .
Indeed since µ 0 is bounded above and below on its support (see (2.2)), m K also is. Combining with the above, (6.8) and (6.11), and using the definition of w, (2.13) and the Cauchy-Schwarz inequality we are led to

K |(E K ) η | 2 ≤ m K |K|κ d w(η) + |K| min Am K W η + o R (1) + C R n -2/d + C R,η n -1/d (6.12) if d = 3; respectively, if d = 2, K |(E K ) η | 2 ≤ m K |K| κ 2 w(η) + γ 2 - 1 2 κ 2 log m K +|K| min Am K W η + κ 2 m K 2 log m K + o R (1) = m K |K|(κ 2 w(η) + γ 2 ) + |K| min Am K W η + o R (1) + C R n -1 + C R,η n -1/2 (6.13)
The electric field in Σ ′ int is then set to be E int = K∈Kn E K . We can extend it by 0 outside of Σ ′ int , and it then satisfies

-div E int = c d p∈Λ int δ p -µ ′ 0 in R d (6.14)
for some discrete set Λ int . Indeed, no divergence is created at the interfaces between the hyperrectangles since the normal components of ∇u K and E are zero. In view of (5.24) all the points in Λ int are of simple multiplicity and at a distance > r 0 10 µ 0 L ∞ from all the others. Note that integrating (6.14) we have

#Λ int = Σ ′ int µ ′ 0 ,
which is an integer.

-Step 2: we define the configuration near the boundary. Since ∂Σ ′ ∈ C 1 , the set

Σ ′ bound := Σ ′ \Σ ′ int is a strip near ∂Σ ′ of volume ≤ Cn d-1 d and width ≥ cR by Lemma 6.3. Since Σ ′ µ ′ 0 = n, Σ ′ bound µ ′ 0 is also an integer. We just need to place Σ ′ bound µ ′ 0 points in Σ ′
bound , all separated by distances bounded below by some constant r 0 > 0 independent of n, η, and R (up to changing r 0 if necessary). Proceeding as in [START_REF] Sandier | 2D Coulomb gases and the renormalized energy[END_REF]Section 7.3,Step 4], using the fact that ∂Σ ′ is C 1 , we may split Σ ′ bound into regions C i such that C i µ ′ 0 ∈ N and C i is a set with piecewise C 1 boundary, containing a ball of radius C 1 R and contained in a ball B i of radius C 2 R, where C 1 , C 2 > 0 are universal. We then place C i µ ′ 0 points in C i , in such a way that their distances (and their distance to ∂C i ) remain bounded below by r 0 > 0, and call Λ i the resulting set of points. We then define v i to solve (6.15) and extend ∇v i by 0 outside B i , this way, we have globally

-∆v i = c d p∈Λ i δ p -µ ′ 0 1 C i in B i ∇v i • ν = 0 on ∂B i ,
-div (∇v i ) = c d p∈Λ i δ p -µ ′ 0 1 C i in R d .
We then set

E bound := i ∇v i .
We can also check that, arguing as in [SS4, Section 7.3], the energy of E bound is bounded by a constant depending on the number of points involved, times the volume of the boundary strip, that is

i B i (E bound ) η 2 ≤ C R,η n 1-1/d . (6.16)
-Step 3: we define the configuration globally and evaluate the energy. We set

E = E bound + E int in Σ ′
and extend it by 0 outside Σ ′ , and we let (6.17) We also have #Λ = Σ ′ µ ′ 0 = n and we can define thus our test configuration as

Λ = Λ int ∪ ∪ i Λ i . Then E satisfies -div E = c d p∈Λ δ p -µ ′ 0 in R d .
x = {x i = n -1/d x ′ i } n i=1 where Λ = (x ′ 1 , . . . , x ′ n ). (6.18)
Note that it depends on R and η. There remains to bound H n (x 1 , . . . , x n ) from above. Since all the points in Λ are separated by distances > r 0 fixed, the points in the collection (x 1 , . . . , x n ) are separated by distances > r 0 n -1/d , where r 0 is independent of η. We may then apply Lemma 3.6 with ℓ = ηn -1/d , and if η is small enough, we have |x i -x j | ≥ 2ℓ. We are then in the case of equality in that lemma: (6.19) where we used that all the points are in Σ where the function ζ vanishes, and by definition of J n , and letting h ′ n,η be as in (3.18),

H n (x 1 , . . . , x n ) -n 2 E[µ 0 ] + n 2 log n 1 d=2 ≤ n 2-2/d J n (x 1 , . . . , x n ) + Cη 2 ,
J n (x 1 , . . . , x n ) := 1 c d 1 n R d |∇h ′ n,η | 2 -(κ d w(η) + γ 2 1 d=2 ) ≤ 1 c d 1 n R d |E η | 2 -(κ d w(η) + γ 2 1 d=2 ) . (6.20)
The last inequality is a consequence of (6.17):

R d |E η | 2 = R d |∇h ′ n,η | 2 + R d |E η -∇h ′ n,η | 2 + 2 R d ∇h ′ n,η • E η -∇h ′ n,η = R d |∇h ′ n,η | 2 + R d |E η -∇h ′ n,η | 2 -2 R d h ′ n,η div E η -∇h ′ n,η = R d |∇h ′ n,η | 2 + R d |E η -∇h ′ n,η | 2
where the integration by parts is justified by the decay at infinity 9 of E η and ∇h ′ n,η , and we have div (E η -∇h ′ n,η ) = 0 by definition (3.17). Next we recall that by construction of E we may write

R d |E η | 2 ≤ K∈Kn K |(E K ) η | 2 + C R n 1-1/d (6.21)
where we used (6.16) to estimate the contribution of the boundary terms, and summing the bounds (6.12)-(6.13) over K and inserting into (6.21), we obtain (6.22) 9 Remark that the right-hand side of (3.18) always has zero total charge and is compactly supported.

c d nJ n (x 1 , . . . , x n ) ≤ (κ d w(η) + γ 2 1 d=2 ) K∈Kn m K |K| -n + K∈Kn |K| min Am K W η + o R (1) + C R n 1-1/d .

Low temperature regime

The result is proved by finding upper and lower bounds to F n,β , which are themselves direct consequences of the upper and lower bounds we have obtained for H n . We recall that F n,β is linked to the partition function via (2.26), so that

F n,β = - 2 β log R dn e -β 2 Hn . ( 7
) ∈ R dn H n (x 1 , . . . , x n ) ≥ n 2 E[µ 0 ] - n 2 log n 1 d=2 + n 2-2/d ξ d + 2n n i=1 ζ(x i ) + o n (n 2-2/d )
where o n (1) → 0 when n → ∞, and ξ d is defined in (2.21). Consequently,

e -β 2 Hn(x 1 ,...,xn) ≤ exp - n 2 β 2 E[µ 0 ] + βn 4 log n 1 d=2 - n 2-2/d β 2 (ξ d + o n (1)) n i=1 exp (-nβζ(x i ))
and we deduce from (7.1), separating variables when integrating n i=1 exp (-nβζ(

x i )) over (R d ) n , that F n,β ≥ n 2 E[µ 0 ] - n 2 log n 1 d=2 + n 2-2/d (ξ d + o n (1)) - 2n β log R d exp(-nβζ(x)) dx .
On the other hand, with the assumption (2.24), the dominated convergence theorem yields that

R d exp(-nβζ(x)) dx → |{ζ = 0}| = C as nβ → +∞ (7.2)
where C is a constant. Note that the assumptions on β in Theorem 3, Item 1, ensure that βn → ∞ as n → ∞. The lower bound corresponding to (2.28) -(2.29) then follows.

Upper bound. The key tool for the upper bound is Proposition 6.1. For any ε we have an r 1 (depending on ε) and a set A n as described therein and we may write

F n,β ≤ - 2 β log An e -β 2 Hn .
This corresponds to taking as a trial state for the free energy functional F n,β the probability measure An) , where 1 An is the characteristic function of the set A n . Using (6.2)-( 6.3) we deduce

1 An P n,β P n,β (
F n,β ≤ - 2 β log |A n | + n 2 E[µ 0 ] - n 2 log n 1 d=2 + n 2-2/d (ξ d + ε) .
Using the estimate on |A n | and Stirling's formula we have

log |A n | ≥ n log πr d 1 e -C
and thus

F n,β ≤ n 2 E[µ 0 ] - n 2 log n 1 d=2 + n 2-2/d (ξ d + ε) - 2n β log πr d 1 e + C β . (7.3)
The upper bound corresponding to (2.28)-(2.29) follows, with C ε depending only on log r 1 and thus bounded when ε is bounded away from 0.

High temperature regime

In this regime we will use the mean-field density at positive temperature µ β . We recall some of its properties in a lemma. The proof uses arguments already provided in [RSY2, Section 3.2] and is omitted.

Lemma 7.1 (The mean-field density at positive temperature).

For any β > 0 the functional (2.27) admits a unique minimizer µ β among probability measures. It satisfies the bounds

0 < µ β ≤ C on R d (7.4)
where C is some constant depending only on the dimension and the potential V . Moreover we have the variational equation

2h µ β + V + 2 nβ log(µ β ) = F[µ β ] + D(µ β , µ β ) on R d (7.5) where h µ β = w * µ β . (7.6)
An upper bound to the n-body free energy is obtained by taking the trial state µ ⊗n β in (2.25). Independently of the dimension this yields

F n,β ≤ F n,β [µ ⊗n β ] = n 2 F[µ β ] -nD(µ β , µ β ) ≤ n 2 F[µ β ] -Cn (7.7)
where the second term is due to the fact that there are n(n -1)/2 and not n 2 /2 pairs of particules. The upper bound in (2.30), (2.31) then follows from the definition (2.27).

For the lower bound we use a variant of Lemma 3.6.

Lemma 7.2 (Alternative splitting lower bound).

For any x 1 , . . . , x n ∈ R d , we have

H n (x 1 , . . . , x n ) ≥ n 2 F[µ β ] - 2 β i=1 n log µ β (x i ) - n 2 log n 1 d=2 -Cn 2-2/d , (7.8)
where C only depends on the dimension.

Proof. We only sketch the proof since it follows exactly that of Lemma 3.6. Applying Onsager's Lemma 3.4 with ℓ = n -1/d and µ = nµ β , we find

i =j w(x i -x j ) ≥ D nµ β - n i=1 δ (ℓ) x i , nµ β - n i=1 δ (ℓ) x i -n 2 D(µ β , µ β ) + 2n n i=1 D(µ β , δ (ℓ) x i ) -nD(δ (ℓ) 0 , δ (ℓ) 0 ) ≥ -n 2 D(µ β , µ β ) + 2n n i=1 D(µ β , δ (ℓ) x i ) -n κ d c d w(ℓ) + γ 2 c 2 1 d=2 ≥ -n 2 D(µ β , µ β ) + 2n n i=1 h µ β (x i ) -n κ d c d w(ℓ) + γ 2 c 2 1 d=2 -Cn 2 ℓ 2 ≥ n 2 F[µ β ] -n n i=1 V (x i ) + 2 βn n i=1 log(µ β (x i )) -n κ d c d w(n -1/d ) -Cn 2-2/d .
But the result of Proposition 4.1 and Remark 4.4 give us that for any

x n ∈ A n , lim inf n→∞ n 2/d-2 H n (x n ) -n 2 E(µ 0 ) + n 2 log n 1 d=2 ≥ |Σ| c d W(E) dP (x, E)
where P = lim n P νn and P νn is defined as in (2.19). Since P νn ∈ i n (A n ) we have P ∈ A ∞ , where A ∞ is as in the statement of the theorem, by definition, and thus lim inf

n→∞ n 2/d-2 H n (x n ) -n 2 E(µ 0 ) + n 2 log n 1 d=2 ≥ |Σ| c d inf P ∈A∞ W(E) dP (x, E).
Inserting this into (7.12) and using (2.28)-(2.29), we find, for every ε > 0 lim sup

n→∞ log P n,β (A n ) n 2-2/d ≤ - β 2 |Σ| c d inf P ∈A∞ W(E)dP (x, E) -(ξ d + ε) -C ε lim n→∞ n 2/d-1 β + lim sup n→∞ 1 n 2-2/d log An exp -β n i=1 nζ(x i ) dx 1 . . . dx n
and there only remains to treat the term involving ζ as before using (7.2) to arrive at the desired result (2.33).

We next turn to the tightness of P n,β . Starting from (7.12) and inserting the upper bound on F n,β implied by (2.28) and the lower bound (3.30), we obtain that for any event A n , for any η ≤ 1,

P n,β (A n ) ≤ An exp - β 2 n 2-2/d 1 c d n R d |∇h ′ n,η | 2 -(κ d w(η) + γ 2 1 d=2 ) -Cη 2 n i=1 exp(-nβζ(x i ))dx 1 . . . dx n . (7.13) Let now A n,M := x ∈ (R d ) n , ∀η < 1 2 , 1 n R d |∇h ′ n,η | 2 -(κ d w(η) + γ 2 1 d=2 ) ≤ M .
Inserting into (7.13) and using (7.2), we obtain

P n,β ((A n,M ) c ) ≤ exp Cn + β 2 (Cn 2-2/d - M c d n 2-2/d ) .
Using that β ≥ cn 2/d-1 , it follows that we can find a large enough M (independent of n and β) for which P n,β ((A n,M ) c ) → 0 as n → ∞. To prove the tightness of P n,β , in view of [START_REF] Sandier | 2D Coulomb gases and the renormalized energy[END_REF]Lemma 6.1], it then suffices to show that if P n ∈ i n (A n,M ) then P n has a convergent subsequence. But this has been precisely established in Section 4.2. The fact that the limits of P n,β are concentrated on admissible P 's satisfying W(P ) ≤ ξ d + C β is an easy consequence of what precedes.

Since q ≥ 1, we have 1 d -2 qd + 1 ≥ 0 and thus, for n large enough, inserting into (7.13) we deduce that P n,β ∇h n L q (B R ) ≥ λn t q,d ≤ C exp -C R λ 2 βn tq,d (7.15)

where t q,d and tq,d are as in the statement of the theorem. This ends the proof since, in view of the definition (3.17), -∆h n = c d (ν n -nµ 0 ), and thus

∇h n L q (B R ) = c d ν n -nµ 0 W -1,q (B R ) .

Estimates on reduced densities: proof of Corollary 6

We first prove a simple lemma that formalizes in our setting the idea that a control on the fluctuations of the empirical measure in a symmetric probability measure implies a control of the marginals of that measure. With this additional ingredient in hand, Corollary 6 becomes a consequence of Theorem 5, Item 3. We have used a similar idea in [START_REF] Rougerie | Quantum Hall phases and plasma analogy in rotating trapped Bose gases[END_REF]Lemma 3.6].

Lemma 7.4 (Control of the empirical measure implies control of the marginals). Let µ n be a symmetric probability measure over (R d ) n with reduced densities µ (k) . Let

1 ≤ q < d d -1
and p = q q -1 . Recall the definition (3.17) of h n as a function of x = (x 1 , . . . , x n ). where we used (3.17), the assumption that ϕ has compact support to justify the integration by parts and Hölder's inequality. This proves Item 1 since only the term ∇h n L q in the right-hand side of (7.18) depends on x.

We now turn to Item 2, for which a little bit more algebra is required. We first note that, using the symmetry under particle exchange, (7.19) This is a classical combinatorial trick that seems to go back to [Gru] and can be found e.g. in [START_REF] Hauray | On Kac's chaos and related problems[END_REF]Proof of Lemma 8] or [START_REF] Rougerie | Théorèmes de de Finetti[END_REF]Proof of Theorem 2.2]. In fact we are here using the de Finetti-type theorem of [DF] in disguise. Then, in view of the condition on q we have p > d and thus the Sobolev embedding W 1,p ֒→ L ∞ , so that, in view of (7.19), Proof of Corollary 6. From (7.15) it is easy to deduce that for some large enough λ x∈R dn ∇h n L q dP n,β (x) = { ∇hn L q ≤λn t q,d } ∇h n L q dP n,β (x) + { ∇hn L q ≥λn t q,d } ∇h n L q dP n,β (x) ≤ λn t q,d (1 + o(1)) .

R dk µ (k) ϕ -
Indeed, the second term is negligible because of (7.15) and for the first one we only use the fact that P n,β is a probability. The results of Corollary 6 then follow from Lemma 7.4.

A Existence of a minimizer for W: direct proof

We start by proving that there exists a minimizer for W η . Picking a minimizing sequence E n of admissible electric fields, a standard concentration argument yields a sequence R n of radii such that where P η is a probability measure on L 2 loc (R d , R d ). Recalling (A.1) we thus have inf

A 1 W η = lim n→∞ W η (E n ) ≥ lim R→∞ 1 |K R | χ * 1 K R |E| 2 dP η (E) -(κ d w(η) + γ 2 1 d=2 ). (A.3)
We denote P the push-forward of P η by Φ -1 η , i.e. dP (Φ η (E)) = dP η (E). Recalling that P η is defined using [SS3, Theorem 3], we can argue as in the proof of Proposition 4.1 to prove that P does not depend on η. We then rewrite (A.3) as inf

A 1 W η ≥ W η (E)dP (E) (A.4)
and since W η is bounded below independently of η by Proposition 5.1 we can use Fatou's lemma to pass to the lim inf as η → 0:

inf

A 1 W = lim inf η→0 inf A 1 W η ≥ lim inf η→0 W η (E) dP (E) = W(E)dP (E);
from which it follows that P is concentrated on the set of minimizers of W, which has to be nonempty.

B Comparison of W and W and the periodic case

In this appendix we discuss the relation between the two version of the renormalized energy functionals. In the case of periodic configurations, both coïncide with a simplified expression dependind only on the points. Lemma 3.2 gives a convenient way to bound from below the energy of well-separated charge configurations. We can also use it to compare W and W when d = 2: Proposition B.1 (W and W coincide in 2D for well-separated points). Assume d = 2, and let E ∈ A 1 be such that W(E) < +∞ and the associated set of points satisfies min p =p ′ ∈Λ |p -p ′ | ≥ η 0 > 0 for some η 0 > 0. Then W(E) = W (E).

Proof. Without loss of generality, we may assume that η 0 < 1 2 . By well-separation of the points, for any R > 1 we may find a set U R such that K R ⊂ U R ⊂ K R+1 , and (3.7) is verified in U R with η 0 /2. We may then apply Lemmas 3.2 in U R and obtain

U R |E η | 2 -#(Λ ∩ U R )(κ d w(η) + γ 2 1 d=2 ) = W (E, 1 U R ) + o η (1)#(Λ ∩ U R ). (B.1)
with C depending only on η 0 and d. Next, let χ K R be as in Definition 2.4. We have

W (E, χ K R ) -W (E, 1 U R ) = W (E, χ K R 1 U R \U R-1 ) -W (E, 1 U R \U R-1 ). (B.2)
Note that by construction, the B(p, η 0 /2) do not intersect ∂U R and ∂U R-1 . We may next write, by definition of W , for any 0 < r ≤ η 0 /2, Similarly as above, we have This result should be compared with [START_REF] Sandier | From the Ginzburg-Landau Model to Vortex Lattice Problems[END_REF]Lemma 1.3].

W (E, 1 U R \U R-1 ) = (U R \U R-
W (E, 1 U R \U R-1 ) -W (E, χ K R 1 U R \U R-1 ) ≤ (U R \U R-1 )\∪pB(p,η 0 /2) (1 -χ K R )|E| 2 + p∈Λ∩(U R \U R-1 ) (1 -χ K R (p)) lim
Proof. In view of Lemma 3.2 and the periodicity of h, we easily check that W(∇h) = 1 |T| W (∇h, 1 T ). This is a renormalized energy computation à la [BBH]: first, using Green's formula and the equation satisfied by h, we compute

(c d R(x -a i ) + f (x -a i ))∇h • ν
where f is a C 1 function equal to 0 at 0. We then use that, by Green's theorem, Inserting into (B.8), in view of the definition of W (∇h, 1 T ), we get the result.

If there is a multiple point in a periodic configuration, it is easy to see that both W and W are +∞ for this configuration, e.g. as a limit case of the above result. Both ways of computing the renormalized energy are thus perfectly equivalent for any periodic configuration.

Configurations that form a simple lattice correspond to this situation but with only one point in the torus, in that case W is thus equal to c 2 d R = lim x→0 c 2 d G(x) -c d w(x) . In addition, we may compute explicitly the Green's function of the torus, using Fourier series. Using the normalization of the Fourier transform f (y) = where Λ denotes the lattice.

The series that appears here is an Eisenstein series E Λ (x) = p∈Λ * \{0} e 2iπp•x 4π 2 |p| 2 . In dimension 2, using the "first Kronecker limit formula" (see also [START_REF] Sandier | From the Ginzburg-Landau Model to Vortex Lattice Problems[END_REF] for an analytic proof), one can show that this series is related to the Epstein Zeta function of the lattice Λ * defined by Thus minimizing W over lattices reduces to minimizing the Zeta function over lattices Λ of volume 1, and this question was solved by Cassels, Rankin, Ennola, Diananda, in the 60's (see a self-contained proof in [Mon]) in dimensions 2, where the unique minimizer is the triangular lattice (the minimizer is also identified in dimensions 8 and 24). Minimizing the Zeta function over lattices remains an open question in dimension ≥ 3, even though the face-centered cubic (FCC) lattice is conjectured to be a local minimizer (cf. [SaSt]). In dimension d ≥ 3 however, (B.10) is not proven and it is not clear how to give it a meaning since the series involved are all divergent (one would at least need to use the meromorphic extension of the Zeta function), see e.g. [Lan, Sie].

  .11) We next observe that c d R d w * δ

  locally bounded and up to extraction converges weakly as n → ∞, in the sense of measures, to ν = p∈Λ N p δ p

  first equation of (5.41) yields ∂K t+α g ∈ c d N because by construction |K t+α | ∈ N. Since we have also ensured m i |K i | ∈ N we deduce that i |K i | ∈ N. Recalling (5.37) we deduce that

1.

  For any ϕ ∈ C ∞ c (R d ) we have R d (µ (1) -µ 0 )ϕ ≤ 1 c d n ∇ϕ L p x∈R dn∇h n L q µ(x)dx. (7.16)2. For any ϕ ∈ C ∞ c (R dk ) symmetric in the sense that for any permutation σ ϕ(x 1 , . . . , x k ) = ϕ(x σ(1) , . . . , x σ(k) )we haveR dk (µ (k) -µ ⊗k 0 )ϕ ≤ k c d n x∈R dn ∇h n L q µ(x)dx + C k 2 n sup x 1 ∈R d . . . sup x k-1 ∈R d ∇ϕ(x 1 , . . . , x k-1 , . ) L p (R d ) . (7.17)Proof. Using the symmetry of µ we write

δ

  x i -nµ 0 ϕ = 1 c d R d ∇ϕ • ∇h n ≤ 1 c d ∇ϕ L p ∇h n L q (7.18)

RRδ

  1 , . . . , x k-1 , . ) L p (R d ) (7.20) follows from the above.We then use this and the triangular inequality to obtainR dk (µ (k) -µ ⊗k 0 )ϕ ≤ n -k dk 1≤i 1 ,...,i k ≤n δ x i 1 ⊗ . . . ⊗ δ x i k -(nµ 0 ) ⊗k ϕ x k-1 ∈R d ∇ϕ(x 1 , . . . , x k-1 , . ) L p (R d ) dk 1≤i 1 ,...,i k ≤n δ x i 1 ⊗ . . . ⊗ δ x i k -(nµ 0 ) ⊗k ϕ .On the other hand, recalling that nµ 0 TV = n i=1 δ x i TV = n and using the symmetry of ϕ we have x i -nµ 0 ϕ(x 1 , . . . , x k-1 , . )≤ kn k-1 c d sup x 1 ∈R d . . . sup x k-1 ∈R d ∇ϕ(x 1 , . . . , x k-1 , . ) L p (R d ) ∇h n L q (R d ) (7.23)by definition of the total variation norm and an estimate exactly similar to (7.18). Combining (7.21), (7.22) and (7.23) we obtain the desired result.

  compactly supported χ and introducef n (E) := χ(y)|E| 2 if ∃ λ ∈ K Rn : E = θ λ (E n,η ) +∞ otherwisewhere θ λ is the action of the translation by λ. θ λ (E n,η ))dλ and we may apply the framework of[START_REF] Sandier | 2D Coulomb gases and the renormalized energy[END_REF] Theorem 7], or even the simpler version in[START_REF] Sandier | From the Ginzburg-Landau Model to Vortex Lattice Problems[END_REF] Theorem 3]. We use the latter. Arguing as in the proof of Lemma 4.2, Assumptions 1 and 2 about the functional f n simply follow from compactness in L 2 loc and lower semi-continuity. Indeed, ifE n ⇀ E in L 2 loc (R d , R d ) we have lim inf n→∞ f n (E n ) ≥ f (E)wheref (E) := χ(y)|E| 2 if ∃ E ′ : E = Φ η (E ′ ) +∞ otherwise. Applying [SS3, Theorem 3] we deduce lim n→∞ -K Rn |E n,η | 2 = lim inf n→∞ -K Rn f n (θ λ (E n,η ))dλ ≥ lim R→∞ -K R f (θ λ E)dλ dP η (E).

  1 )\∪pB(p,r)|E| 2 + p∈Λ∩(U R \U R-1 ) lim η→0 B(p,r)\B(p,η) |E| 2 -c d w(η) ,and using (3.13) (applied to r instead of η 0 /2) we havelim η→0 B(p,r)\B(p,η) |E| 2 -c d w(η) ≥ -c d w(r) -Cfor each p (where C depends only on d). It follows that we may write(U R \U R-1 )\∪pB(p,r) |E| 2 + p∈Λ∩(U R \U R-1 ) lim η→0 B(p,r)\B(p,η) |E| 2 -c d w(η) ≤ W (E, 1 U R \U R-1 ) + (c d w(r) + C)#(Λ ∩ (U R \U R-1 )). (B.3)

  η→0 B(p,η 0 /2)\B(p,η) |E| 2 -c d w(η) + p∈Λ∩(U R \U R-1 ) lim η→0 B(p,η 0 /2)\B(p,η) |χ K R -χ K R (p)||E| 2 (B.4)In view of (B.3) applied with r = η 0 /2 we can bound the first two terms on the right-hand side by W (E,1 U R \U R-1 ) + C#(Λ ∩ (U R \U R-1 )),where C depends only on d and η 0 . We turn to the last term. Let us set φ(r) = (U R \U R-1 )\∪pB(p,r) |E| 2 .

  outer unit normal. The second term tends to 0 as η → 0 since h has mean zero. For the first term we note that h = c d N i=1 G(x -a i ), and that G(x) = c -1 d w(x) + R(x) with R a C 1 function, and insert this to find∂B(a i ,η) h∇h • ν = w(η) + c d j =i G(a i -a j ) ∂B(a i ,η) ∇h • ν + ∂B(a i ,η)

  ∂B(ai ,η) ∇h • ν = B(a i ,η) ∆h = -c d + o η (1), and that |∇h|(x) ≤ C|∇w|(x -a i ) + C ≤ C|x -a i | 1-d to conclude that lim η→0 ∂B(a i ,η) -h∇h • ν -c d w(η) = c 2 d j =iG(a i -a j ) + c 2 d R(0).

  is the dual lattice of Λ i.e. the set of q's such that p • q ∈ Z for every p ∈ Λ. By Fourier inversion formula we obtain the expression of G in Fourier seriesG(x) = p∈Λ * \{0} e 2iπp•x 4π 2 |p| 2 .

E

  Λ 1 (x) -E Λ 2 (x) = C lim s>0,s→0 ζ Λ * 1 (s) -ζ Λ * 2 (s). (B.10)

  .1) Lower bound. Proposition 4.1 and Remark 4.4 yield that for any (x 1 , . . . , x n

If β ≪ n 2/d-1 (high temperature regime) in the limit n → ∞, the next-to leading order

As in[START_REF] Sandier | From the Ginzburg-Landau Model to Vortex Lattice Problems[END_REF] we could define the renormalized energy with averages on more general nondegenerate (Vitali) shapes, such as balls, etc, and then prove that the minimum of W does not depend on the shapes.

Remark that the symmetry constraint (1.11) is automatically satisfied even if do not impose it in the minimization.

The opposite case is in fact somewhat easier, see[START_REF] Rougerie | Quantum Hall phases and plasma analogy in rotating trapped Bose gases[END_REF] Section 3].

We could easily prove the same result with averages on balls or other reasonable shapes.

Translating the whole system we may always assume that 0 ∈ Λ

We have Γ ⊂ Kt+α ⊂ Γ ′ ⊂ KR. Of these sets, only Kt+α and KR are hyperrectangles.

we split Σ ′ , the blown-up of Σ at scale n 1/d , into hyperrectangles of size R. In each hyperrectangle K we put a configuration of points constructed via Corollary 5.9, properly scaled so that the local electric field is in A m K where m K = -K µ ′ 0 is the mean value of µ ′ 0 in the hyperrectangle K. Thanks to the screening property (5.23) the electric field is then defined globally in the domain by just gluing together the fields defined in each hyperrectangle. Since all the points in the test configuration so constructed are well-separated, we are in the case of equality of Lemma 3.6 so the next to leading order in the energy upper bound will exactly be given by the electrostatic energy of charges smeared on a length ηn -1/d around the positions of our configuration, up to remainder terms that will ultimately become negligible. The conclusion in the limit R → ∞ (tiling on a scale much larger than the interparticle distance) and η → 0 (smearing charges on a scale much smaller than the interparticle distance) then follows by splitting this energy into the contribution of each hyperrectangle, using the scaling property (2.13) and a Riemann sum argument to recover the µ 0 -dependent factors in the right-hand sides of (6.2), (6.3). Of course there is a boundary layer near ∂Σ ′ that we cannot tile properly, and we will have to complete the construction there in a way that does not cost too much energy.

Proof. -Step 1: we define the configuration in the interior of Σ ′ and estimate its energy. We need to tile the interior of Σ ′ into hyperrectangles K ∈ K n and put in each a number of points equal to the charge of the background µ 0 . This requires that K µ ′ 0 be an integer 8 for any K. The following lemma provides such a tiling. It corresponds to [START_REF] Sandier | 2D Coulomb gases and the renormalized energy[END_REF]Lemma 7.5] (with q = 1) where the proof is provided in dimension 2. The adaptation to higher dimension is immediate and left to the reader. Lemma 6.3 (Tiling the interior of Σ ′ ). There exists a constant C 0 > 0 such that, given any R > 1, there exists for any n ∈ N * a collection K n of closed hyperrectangles in Σ ′ with disjoint interiors, whose sidelengths are between 2R and 2R + C 0 /R, and which are such that

Note that the sizes of the hyperrectangles are controlled thanks to the positive lower bound assumption on the density µ ′ 0 . We apply this lemma, which yields a collection K n . For each K ∈ K n we denote

We need to correct for the difference between m K and µ ′ 0 by setting u K to be the solution to

First we note that from (6.7)

In view of the regularity of µ 0 (which implies m

) and the properties of our tiling K n , with (6.7) again, we can then also recognize a Riemann sum to see that

using the continuity of m → min Am W η which can be checked from (2.12), (2.13) and the continuity of η → Φ η . The proof is concluded by dividing (6.22) by c d n, using (6.19), and taking successively n → ∞, R large enough and then η small enough (and changing the configuration of points accordingly). Using that min Am W η → min Am W along a sequence η → 0, by definition, we can make

Finally, inserting (2.13) gives the upper bound result at the configuration (x 1 , . . . , x n ). There remains to perturb the points and prove the statement about the volume of the set A n . We note that (still for the R large enough and η small enough that (6.23) holds), the construction of (x ′ 1 , . . . , x ′ n ) has been made (at the blown up scale) by gluing together configurations over an order n/R d of hyperrectangles K ∈ K n (or cells near the boundary C i ) of volumes O(R d ), each containing a number of points O(R d ) which are at distances bounded below by r 0 > 0. If each of these points is moved by r 1 < r 0 /2 small enough, by continuity of the energy in each box 10 , this induces an error in |E η | 2 in each hyperrectangle which can be made < ε 2CR d if r 1 is chosen small enough (depending on ε, R, η), for example in each estimate (6.12). Choosing C large enough depending on R, and summing these errors over all the hyperrectangles leads to an error 1 2 nε, which, once divided by n, yields an error term ε 2 in (6.20). Continuing on leads to an error ε instead of ε/2 in (6.23). We have thus shown that we have the desired estimate in the set A n consisting of those (y 1 , . . . , y n ) such that sup i (y

), it is raised to the power n because there are n points in the configuration, and multiplied by n! because permuting y 1 , . . . , y n does not change the energy.

7 Applications to the partition function and large deviations 7.1 Estimates on the partition function: proof of Theorem 3

In this section we prove Theorem 3.

10 which follows from the techniques of Section 5.

The second inequality is obtained by dropping a positive term and using (3.5), the third follows from Lemma 3.5 and (7.4) as in (3.28) and the fourth from the variational equation (7.5) and computing w(n -1/d ) explicitly. Plugging this inequality in the expression of the Hamiltonian, we get the result.

Using the results of Lemma 7.2 and the expression (2.25) we deduce that for any probability measure µ

The integral term in this equation is (minus) the entropy of µ relative to µ ⊗n β and is thus positive (see e.g. [RSY2, Lemma 3.1]) for any probability measure µ ∈ P(R dn ). Dropping this term, combining this lower bound with the upper bound (7.7), and noting that n 2-2/d ≫ n for n large and d ≥ 2, we conclude the proof of Theorem 3, Item 2.

Remark 7.3 (Total variation estimates for reduced densities). Instead of simply dropping the relative entropy in (7.9) for our final estimate of the n-body free-energy we may combine our upper and lower bounds to control this term. Interesting estimates then follow from the coercivity properties of the relative entropy. Indeed, using subadditivity of entropy (see e.g. [Kie2, Proposition 2])

where µ is a symmetric probability on R dn , µ (k) is its k-th marginal, ⌊ . ⌋ stands for the integer part and n[k] is n modulo k. On the other hand

By positivity of relative entropies the contribution of the difference of the second terms in equations (7.10) and (7.11) is non negative. One can then use the Csiszár-Kullback-Pinsker inequality (see e.g. [RSY2, Lemma 3.1] and [BV] for a proof ) to bound from below the difference of the first terms and obtain

where . TV stands for the total variation norm. A control on the n-body relative entropy thus provides estimates in the spirit of Corollary 6, but it turns out (combining our free energy upper and lower bounds) that those are meaningful only in the high temperature regime. See [START_REF] Rougerie | Quantum Hall phases and plasma analogy in rotating trapped Bose gases[END_REF] (in particular Remark 3.4 and Section 3.4) where this method is used in such a regime.

The Gibbs measure at low temperature: proof of Theorem 4

Let A n be any event, i.e. subset of (R d ) n . Using (2.26) we may write

. . dx n . (7.12)

Charge fluctuations: proof of Theorem 5

We start from (7.13) applied with some arbitrary 0 < η < 1 and insert the result of Lemma 3.8. We deduce that for any event A n ,

where we recall that D(x ′ , R) depends on (x 1 , . . . , x n ) ∈ A n via the "empirical measure" ν ′ n . Then, using (7.2), we find Proof of Item 1. We pick a sequence of microscopic balls of radii R n . According to (3.33) and using (7.14), we have d+2) . If λ ≤ 1/n then the inequality is trivially true for n large enough anyway. It follows that (2.36) holds.

Proof of Item 2. The argument is similar as above: Fixing some radius R, we have

Proof of Item 3. We start again from (7.13). Using Lemma 3.9 in the ball of radius Rn 1/d , the fact that the total mass of ν ′ n is n, and then a change of variables, we have

In particular, for λ ≥ 2C, we have

Since χ K R is Lipschitz, we may write

Using (B.3) to bound φ(r), and using the integrability of w(r) = -log r near 0 (this is the one point where we use that d = 2), we deduce that the third term in (B.4) can also be bounded by

Combining with (B.1) and (B.2) we are led to

On the other hand, since W(E) < ∞, Lemma 3.1 applies and gives that lim R→∞

Remark B.2. The fact that Lemma 3.2 holds in any dimension d ≥ 2 but Proposition B.1 only for d = 2 seems to indicate that the cutoff procedure with χ K R is probably not adapted for dimension d ≥ 3 and the value of W defined in (2.16)-(2.17) may depend on the choice of cutoff χ K R , contrarily to what happens for d = 2 as proven in [START_REF] Sandier | From the Ginzburg-Landau Model to Vortex Lattice Problems[END_REF]. This has no consequence for us however as we will never need the result of Proposition B.1.

We now turn to consequences of the previous result: when configurations are periodic with simple points, then these are automatically well-separated, and we can deduce an explicit expression for W in terms of the points only.