Universality of the third-order phase transition in the constrained Coulomb gas

The free energy at zero temperature of Coulomb gas systems in generic dimension is considered as a function of a volume constraint. The transition between the 'pulled' and the 'pushed' phases is characterised as a third-order phase transition, in all dimensions and for a rather large class of isotropic potentials. This suggests that the critical behaviour of the free energy at the 'pulled-to-pushed' transition may be universal, i.e., to some extent independent of the dimension and the details of the pairwise interaction.


Introduction
In this paper, we are interested in some aspects of the mean field approximation in the statistical mechanics of systems with long-range interactions. For these systems, it is possible to characterise the ground state in the thermodynamic limit and compute the (free) energy at zero temperature. In the simultaneous limit of large number of particles and zero temperature, the system typically concentrates in a bounded region Γ ⊂ R n and the free energy per particle attains a finite value. In presence of volume constraints-i.e. if the system is forced within a specified region of space-the ground state and the free energy may be altered. Since phase transitions are associated with the breakdown of the analyticity of thermodynamic potentials, it is legitimate to ask whether the dependence of the free energy on a volume constraint is analytic or not. If not, is it possible to establish what the regularity of the free energy-i.e., the order of the phase transition-is?
Suppose, for instance, that one constrains the system to be completely contained in a certain region, say a ball B R ⊂ R n of radius R > 0. If R is large enough so that Γ ⊂ B R , the constraint is ineffective (pulled phase) and the free energy does not change.
On the other hand, if R is so small that Γ ⊂ B R , then the equilibrium configuration of the system changes (pushed phase) and the free energy increases. It is clear that there exists a critical radius R that separates the two phases; for R > R the system is in the pulled phase; for R < R the system is in a pushed phase. In order to go beyond this qualitative picture, sketched in Fig. 1, it would be desirable to assign a thermodynamic meaning to this 'pulled-to-pushed' phase transition and compute its order.
Examples of this transition abound in the physics literature of particles interacting via logarithmic potential (log-gas) in dimension n = 1 and n = 2. These systems are related to the eigenvalue statistics of unitarily invariant random matrix models [33] (see Section 2). In those examples, a thermodynamic potential typically turns out to have a discontinuity in the third derivative at the critical value R . The conclusion is that the 'pulled-to-pushed' transition is a third-order phase transition (in the sense of Ehrenfest).
Perhaps, the most popular third-order transitions are the so called Gross-Witten-Wadia [22,40] and Douglas-Kazakov [14] large-N phase transitions. In both cases, the authors realised that even the simplest integrals involving the unitary group were nonanalytic in the coupling constant. In particular, at the critical value of the coupling constant separating the weak coupling and the strong coupling phases, the free energy has a jump in the third derivative. These models are naturally mapped onto the statistical physics of log-gases, and the strong-to-weak coupling transition of lattice gauge theories can be translated as a pulled-to-pushed transition for classical particles. Several other third-order phase transitions have been discovered over the last decade in many physics problems related at various levels to random matrices, such as the distribution of the maximum height of non-intersecting Brownian excursions [19,38], conductance and shot noise in chaotic cavities [7,39,10], Rényi entanglement entropy of bipartite random pure states [17,35,13], wireless telecommunication [24], complexity Figure 1. Illustration of the pulled-to-pushed transition for a constrained gas. At equilibrium, the gas is concentrated in a region Γ (here the disk within the red dotted line). Right: For R > R the constraint (solid blue line) is ineffective, and the free energy is constant. Left: When R < R , the free energy increases as the gas gets more and more pushed within a much narrower region than it would normally occupy at equilibrium. At the critical point R = R the free energy is not analytic.
of spin glass landscapes [20], and discrete log-gases related to random tilings [6] (see the review [31] and references therein for more examples and a detailed discussion).
Based on the available evidence so far, the natural question that arises is: to what extent is the third-order pulled-to-pushed phase transition universal ? In other words, what are the ingredients-dimension, type of interactions, properties of the confining potential, etc.-that are necessary to induce this type of transition? In a recent work [8], three of us suggested that the order of the phase transition should be related to the regularity of the equilibrium measure as a function of the external constraint. In this work, we elaborate further on this idea. The prominent role played by log-gases in R or R 2 in the literature on the topic might lead to believe that the appearance of such weak transitions must be inextricably linked to this specific type of pairwise interaction and/or to low dimensions. Yet, we show here that this is not the case.
More precisely, consider a system of N classical charged particles in R n . The particles interact via the d-dimensional Coulomb potential Φ d (x), the free space Green's function for −∆ in dimension d (see Eq. (14) below). Notice that the Coulomb interaction is repulsive and long-range. To ensure stability in the absence of constraints, the particles are subject to an isotropic confining potential V (x) which may be generated, e.g., by a fixed neutralising background of opposite total charge −N (this is related to the celebrated jellium model of Wigner). Here we consider the natural case n = d, i.e. d-dimensional Coulomb gases in R d . These systems have recently received much attention (see, for instance, the works of Chafaï, Gozlan and Zitt [4], Rougerie and Serfaty [37] and Leblé and Serfaty [26]).
For these Coulomb gas models, we establish that, if the confining potential V (x) is radially symmetric, a constraint on the volume available to the system induces a thirdorder phase transition, irrespective of the dimension and under mild assumptions on V (x). Furthermore, a general formula for the excess free energy in all dimensions is also derived (see Eq. (31) below). This third-order transition was established by an explicit computation in a previous work [9] for the jellium model in d = 2 (see also [2,23]). In this paper, we make the most of the isotropy assumption and the basic properties of the d-dimensional Coulomb interaction in R d to reduce the many-body problem in dimension d to an integrable one-dimensional system.
The paper is organised as follows. In the next Section, we review two examples of third-order phase transitions in random matrix theory. Then in Section 3 we define the model and its thermodynamic limit, and we present and discuss the main result. In Section 4 we derive the main formula on the large deviation functions. As an illustrative example, in Section 5 we present in details explicit formulae for the quadratic case, related to random matrices and the jellium model. The final Section is a summary with a pointer to open problems.

Examples of third-order phase transitions for log-gases
Third-order phase transitions have been observed in one-dimensional and twodimensional systems with logarithmic repulsion, i.e. for eigenvalues of unitarily invariant matrix models. A rather extensive discussion for Hermitian models is contained in [31]. In this section we discuss two paradigmatic examples of random matrices [15,21,33]. (ii) M is a N × N complex matrix whose N 2 entries are independent complex standard Gaussian. The eigenvalues of M are generically complex. This ensemble is called Ginibre Unitary Ensemble (GinUE).
The key fact of these ensembles is that the joint probability density of the eigenvalues (x 1 , . . . , x N ) of M is explicitly known [33]. After a suitable rescaling on the variances of the entries M ij , one finds that where (Of course, the value of the normalization constant Z N is different in the two ensembles. With a slight abuse of notation, we will not distinguish between a measure and its density, and we identify C R 2 .) As remarked by Dyson [15], the joint density (1) can be seen as the canonical measure of a system of particles interacting via logarithmic repulsion in a quadratic external potential at inverse temperature β = 2. In the GUE ensemble, the particles are constrained on the real line (n = 1); in the GinUE ensemble the particles live on the plane (n = 2). With this picture in mind, Z N in both cases is the partition function of a d = 2 Coulomb gas (− log |x| is the electrostatic potential in dimension two): where n = 1 if M ∈ GUE, and n = 2 if M ∈ GinUE. Already at this early stage, a moment of reflection brings up a rather curious feature of these particle systems: for the GUE ensemble, there is an evident mismatch between the physical dimension where the particles live (n = 1) and the Coulomb potential the particles feel (d = 2). We will come back to this simple observation in the Conclusion section. To avoid confusion, we will always denote as log-gas a system of particles in R n (for generic n ≥ 1) repelling with a logarithmic pairwise interaction.
For the GUE and GinUE, the eigenvalues empirical distribution N = 1 N N k=1 δ x k converges to the celebrated 'Wigner law' and 'circular law', respectively: Here 1 A denotes the indicator function of the set A, having the value 1 for all x ∈ A and 0 for all x / ∈ A. The important feature is that, in both cases, the log-gas density concentrates on a bounded region (a ball of radius R ). It is possible to show that, as N → ∞ (and without any scaling), Pr{max |x i | ≤ R} converges to a step function: 0 if R < R , and 1 if R > R . The large-N scaling of these distributions approaching the step function is where the symbol ≈ stands for equivalence at logarithmic scales. The exponential decay is decorated with a 'large deviation function' F (R) which is nothing but the excess free energy of the log-gas constrained to stay in the ball B R = {|x| ≤ R}: The computation of the large deviation functions has been carried out explicitly by Dean and Majumdar [11] for the GUE, and by Cunden, Mezzadri and Vivo [9] for the GinUE ‡. The final result is The form of the large deviation function is specific to the model (F GUE (R) = F GinUE (R)). Nevertheless, in both cases the pulled-to-pushed transition is a third-order phase transition, i.e., the excess free energies are non-analytic with a discontinuous third derivative exactly at the critical point R = R : There is a long list of matrix models exhibiting this third-order singularity. It is therefore tempting to suspect that non-universal large deviation functions of generic statistical models share the same universal critical exponent in presence of volume constraints. However, it also seems that logarithmic interactions must play a prominent, if not essential, role in this respect. Our findings below will instead show that a logarithmic repulsion is not essential to obtaining a third-order phase transition, and the source of this universality must be sought elsewhere.

Definition of the model
We consider here a d-dimensional Coulomb gas in R d . The canonical distribution of N charges in dimension d is with the energy and the partition function Here Φ d (x) is the Coulomb electrostatic potential in dimension d, i.e., the solution of the distributional equation where Ω d = 2π d/2 /Γ(d/2) is the surface area of the unit sphere S d−1 in R d , and Γ is the Euler gamma function. One gets where | · | denotes the Euclidean norm. To ensure the finiteness of the partition function for sufficiently large N , we will require the integrability and the growth conditionŝ

Thermodynamic limit and main result
It is known [30] that entropy plays no role at leading order in N in the asymptotics of the partition function. This can be understood easily by recasting (11) as which shows, at least formally, that the limit N → ∞ is a simultaneous thermodynamic and zero-temperature limit in the mean-field regime. (This readily explains the familiar rescaling in β > 0 of the large deviation functions in random matrix theory.) The energy of a configuration (x 1 , . . . , x N ) can be written in terms of the empirical distribution, as where the mean field energy functional is defined as Under the assumptions (16) on V (x), it is known (see, for instance [4]) that the functional E d has weak- * compact level sets (and is thus lower semicontinuous) and that E d is strictly convex where it is finite. The proof of the convexity is quite standard in d = 1 and d = 2. For d ≥ 3, the proof is based on the fact that the external potential part is linear while the interaction potential can be written as conic combination of Gaussian kernels Compactness and strict convexity imply that E d has a unique minimiser eq , called equilibrium measure, among the probability measures on R d . It turns out that, in the large-N limit, the empirical distribution N converges to the deterministic equilibrium measure eq [4,26,30,37].
In the following, we consider the case of a radial external potential V (x) = v(|x|) satisfying the following hypotheses.
Assumptions A-1: The external radial potential v(r) (r ≥ 0) is of class C 3 . Moreover, we assume that v(r) and r d−1 v (r) are both strictly increasing (the last is true, e.g., if v(r) is strictly convex).
Note that the problem is radially symmetric and this symmetry is inherited by eq . By our assumptions on v(r), it follows that eq is supported on a ball. Indeed, introducing spherical coordinates x = rω ∈ R d , with r ≥ 0 and ω ∈ S d−1 , by an application of Gauss's law, we get and thus The equilibrium measure is concentrated on a ball of radius R , which is given by the unique positive solution of a direct consequence of normalization´d eq = 1.
Following the same line of reasoning of Section 2, the probability that the Coulomb gas is contained in a ball of radius R converges to a step function On the other hand, similarly to (5), one expects where the large deviation function F d (R) decorating the asymptotics (27) is the excess free energy of the constrained d-dimensional Coulomb gas where Z N,β is the free partition function (13), while is the partition function of the Coulomb gas constrained to stay inside the ball of radius R From (26) and (27), we infer that F d (R) is identically zero for R ≥ R and strictly positive when R < R . Hence, the excess free energy is non-analytic at the critical value R = R that separates the pulled and the pushed phases. What is the order of this transition? Our main result is the existence of a third-order phase transition with respect to the parameter R. The precise statement is as follows.
Main result. Let the radial potential V (x) = v(|x|) satisfy the coercivity conditions (16) and the Assumptions A-1. Then, (i) the excess free energy (28) is given by the explicit formula: (ii) the system undergoes a third-order phase transition at R = R , i.e., F d (R) ∈ C 2 (0, ∞), and at the critical value the third derivative is discontinuous: and lim If we accept the claim (i), we can prove the third-order phase transition (ii) as follows. F d (R) given in (31) is manifestly continuous. Note that F d (R) and all its derivatives F (k) d (R) are identically zero for R > R . It remains to study the behaviour of F (k) d (R) as R ↑ R . From (31) we have that, for R < R , We also need the first three derivatives of the Coulomb potential Therefore, from (34), using (25) and (37) we find For the second derivative, using again (25) and (37), we get We analyse now the third derivative. Proceeding as before we find If we write explicitly the condition that Hence we get v (R ) > (1 − d)/R d and using (37) we conclude that This concludes the proof of the third-order transition It remains to prove the claim (i) on the explicit form of the large deviation function.
Remarks. The assumptions on the potential v(r) are not technical. Of course, we require that v ∈ C 3 in order to make sense of the derivatives of F d (R) for R < R [see (34)- (36)]. The condition that v(r) and r d−1 v (r) are strictly increasing implies that the support of equilibrium measure (pulled phase) is simply connected (a ball). It is known that the 'transitions' when the pulled phase is supported on a non-simply connected region display, in general, different exponents (see, e.g., multi-cut Hermitian matrix models). The exponent 3 in (43) is therefore 'universal'. [The non-universal constant in front of (43) is the right hand-side of (40).] For d = 2, the explicit computation of F d (R) in the case of quadratic potential (Ginibre ensemble in random matrix theory) has been done in [9]. In that paper, the authors posed a question on the universality of the third-order transition for normal matrix models. The answer is 'Yes'.

Free energy for radial potentials
In order to compute the large deviation function (28), one needs to extract the leading order in the asymptotics of the partition function of the constrained Coulomb gas (29). Recently, the large-N asymptotics of the free energies of Coulomb gases in dimension d with quite generic potential has been investigated in great detail [4,26,30,37,41]. It has been proved that, for large N , the free energy is dominated by the minimum of the energy functional E d in (20). Then, the following saddle-point approximation holds true where R is the equilibrium measure of the constrained Coulomb gas, i.e., the minimiser of the mean field functional over the set of probability measures supported on the ball of radius R Therefore, the excess free energy (28) is given by Thus the technical problem is i) to find the minimiser R , and ii) compute the minimum energy E d [ R ].

Finding the equilibrium measure
To solve the problem, we exploit the radial symmetry of the system and the fundamental properties of the Coulomb potential. First, we define the total energy density of the constrained gas at Observe that ε d (x, R) is the electrostatic potential at x jointly created by the distribution of charge R and the external potential V . The unknown measure R satisfies the condition of electrostatic equilibrium [4] for some constant C R (independent of x). The 'Euler-Lagrange equations' (48) can be understood, in electrostatic terms, as follows. The fact that ε d (x, R) is constant inside the support of R indicates that the electric field (−∇ε d (x, R)) is zero inside this support. The fact that ε d (x, R) takes greater values outside this support implies that the charges are confined within the support, and that taking any small amount of charge outside the support would increase the total electrostatic energy of the system. In the unconstrained problem, the equilibrium measure eq is supported on the ball of radius R [see (24)- (25)]. Hence, as long as R ≥ R , the constraint is ineffective and the functional is minimised by R = eq . Let us consider now the case R < R , when the Coulomb gas gets pushed by the volume constraint. The minimiser can be identified by imposing two basic conditions valid at electrostatic equilibrium: the Gauss's law must be satisfied and any excess of charge gets localised at the surface of a conductor. Hence, the volume constraint does not change the charge distribution in the interior of the ball and we have Of course´| x|<R d R (x) =´| x|<R d eq (x) < 1. The missing charge is localized on the surface and, by symmetry, is uniformly distributed on the sphere of radius R. Using spherical coordinates, the radial and angular distributions factorise and eventually we find valid when R < R . In the pushed phase, the Coulomb gas equilibrium measure has a singular component (non-fluid phase). At R = R , using (25) one recovers eq .
The correctness of (50) can be ascertained by checking the equilibrium conditions (48). Here we need the following electrostatic integration formula (The identity is a one-line proof for d = 1, a classical formula in complex analysis for d = 2 and a familiar fact for electrostatics in dimension d = 3. In fact, it is possible to prove (51) for generic d ≥ 1 by using the same textbook argument valid in R 3 based on Gauss's law. See [27,Thm. 9.7].) Using (51), we easily verify that (50) satisfies the condition of constant energy density ε d (x, R) on the support of R . For |x| ≤ R, we have where we applied (in order) the electrostatic formula (51), and integration by parts.

Computation of the large deviation function
The next step is to compute the minimum energy E d [ R ]. The simplest way is by observing that Since, from (52), ε d (x, R) = ϕ d (R) + v(R) for |x| ≤ R, and´B R d R (x) = 1, using (50) we find that, for R < R , Thus the large deviation function is given by Note again that R = eq for R ≥ R , and hence F d (R) = 0. Therefore F d (R) can be written as (31) as claimed. This concludes the proof of part (i) of the main result.

Quadratic confinement: random matrices and the jellium model
In this section, we specialise the general result to the case of a quadratic external potential As already discussed in Section 2, in dimension d = 2, the log-gas in a quadratic potential is statistically equivalent to the eigenvalues of the complex Ginibre ensemble (GinUE). In fact, the Coulomb gas in a quadratic potential is also related to the jellium and has a long history in the tradition of classical and quantum statistical physics in dimension d = 1, 2 and 3 (see [1,3,18,23,25,28,42]). The jellium model was invented by Wigner [42] to investigate qualitative properties of electrons in metals and plasmas. Wigner's approximation is quite simple: instead of studying electrons hosted in an ion lattice, one considers an electron gas embedded in a positive continuum and uniform gel (a 'positive jelly'). It is a simple exercise to show that a uniformly charged ball generates a quadratic electric potential inside the ball. For this model in generic dimension d, it is easy to see that, for large N , the normalised charge density at equilibrium converges to the uniform measure on the unit ball (the critical radius is R = 1) This is the unique distribution that neutralises the uniform background, and is the d-dimensional counterpart of the circular law of random matrix theory (d = 2). In presence of a volume constraint, the minimiser of the energy E d among the probability measures supported on the ball of radius R is given by (50) with v(r) = r 2 /2. Explicitly we have Using the general result (31), we eventually find the unified formula for the large deviation function in generic dimension (see Fig. 2) (for d = 2 we take F 2 (R) = lim d→2 F d (R)). For instance, we have for R ≤ 1. For d = 2, Eq. (61) agrees with the result in [9]. For all d ≥ 1:

Conclusions
From a wide perspective, Coulomb gases are, perhaps, the simplest example of statistical mechanical systems for which it is possible to assign a thermodynamic meaning to and to characterise in full the 'pulled-to-pushed' transition. The key elements in the proof of our result are the basic properties of the Coulomb interaction and the isotropy of the models considered. As a consequence, the equilibrium measure of those systems can be written rather explicitly, even in presence of hard constraints. The difficulty in trying to generalise our findings to non-isotropic Coulomb gases is that explicit formulae for the equilibrium measure are not readily available. We conclude with some open problems and a couple of further remarks, inspired by our discussion in Section 2: (i) The excess free energy of the pushed phase provides the left large deviation tail of Pr{max |x i | ≤ R}, when R ≤ R . A standard argument [9,29,32,34] allows to obtain the right large deviation tail as energy cost in pulling one particle from its equilibrium position inside the support of the equilibrium measure and relocating it in a generic position at distance R > R ('pulled' phase). Skipping the details of the computations, one finds where F d (R) is given in (31), and While the functions F d (R) and H d (R) govern the rate of atypical fluctuations away from the unconstrained situation, the question of typical fluctuations is also interesting. One indeed expects the existence of scaling constants a N , b N (possibly dependent on d) and a cumulative function G d (t) (the latter independent of N ) such that In other words, the scaling function G d (t) describes the typical fluctuations of max |x i | in a neighbourhood of the critical point R . In a certain sense, G d is the crossover function between the large deviation functions F d and H d that describe fluctuations of order O(1). It would be interesting to compute G d (t) explicitly for various d and address the question of its (microscopic) universality. For d = 2 (log-gas in the plane), it is known that G d=2 (t) is the Gumbel distribution [5,36] (again under the Assumptions A-1 on V (x)). To the best of our knowledge, the scaling function G d (t) for d = 2 has not been investigated. [Note that for the log-gas on the line (e.g. the GUE) the scaling function is a squared Tracy-Widom [12,16].] (ii) In the GUE model, the equilibrium density (the semicircle law) is continuous and vanishes as a square root (R 2 − x 2 ) 1/2 at the edges (3) as long as R > R ; in the 'pushed' phase R < R , the density diverges as (R 2 − x 2 ) −1/2 at the pushing walls ±R. For the GinUE ensemble, the equilibrium density (the circular law) is uniform in the unit circle (4), and acquires a singular component at the pushing constraining walls. Clearly, the phenomenology of constrained densities for loggases in dimension n = 1 and n = 2 is quite different (see [9,11] for details). Nevertheless the critical exponent is the same. Can a general relation be found between the behaviour of the equilibrium density at the edges and the critical exponent? An intriguing 'scaling argument' was put forward to justify the critical exponent for eigenvalues of Hermitian random matrices (see the review by Majumdar and Schehr [31]). The argument, based on a matching between 'typical' and 'large' deviations, predicts that if the equilibrium density (without constraints) vanishes as (R − x) γ at the edge, then the excess free energy at the critical point behaves as F (R) (R − R) 2(1+γ) . For off-critical Hermitian random matrices, when γ = 1/2, one indeed recovers the exponent 3. Unfortunately, this argument seems unable to explain the occurrence of a thirdorder phase transition in either the GinUE ensemble, or the d-dimensional Coulomb gases considered in this paper (where the equilibrium density is discontinuous at the edge). Hence, the quest for a more viable connection between these seemingly unrelated objects is still open.
(iii) We noticed in Section 2 that for unitary invariant random matrices a mismatch occurs between the dimension of the Coulomb potential (d = 2) and the physical dimension the system lives in (n = 1). Yet, the critical exponent is still 3 (although this has only been proved for a quadratic confining potential). It would be interesting to understand if this feature persists in higher dimensions: for a gas in R n interacting through the d-dimensional Coulomb repulsion (the solution Φ d of the d-dimensional distributional equation (14)), is the critical exponent always 3?
The underlined checkmarks correspond to the log-gases living in dimension n = 1 and n = 2 (e.g. the GUE and GinUE, respectively). In the present work, we have covered the diagonal situation n = d (checkmarks). The above discussion suggests that for n < d the singularity of the pushed density at the constraining walls should become milder and milder as one moves away from the diagonal n = d, but whether the transition is still third-order remains an interesting open problem.
More generally, the main moral to be drawn from the above is that the analytic properties of the excess free energy of constrained particle systems with repulsive interaction may be universal, i.e., to some extent independent of the detailed behaviour of the constrained equilibrium density. One could investigate this point by considering particle systems with interaction kernel G(x) which is the fundamental solution DG(x) = δ(x) of some differential operator D (in this paper D = (−1/Ω d )∆ is the Laplacian). This however is another story which will have to await a future work. It may not be too much to hope that the method and the ideas developed here will find useful applications in other statistical models.