Third-Order Phase Transition: Random Matrices and Screened Coulomb Gas with Hard Walls

Consider the free energy of a d-dimensional gas in canonical equilibrium under pairwise repulsive interaction and global confinement, in presence of a volume constraint. When the volume of the gas is forced away from its typical value, the system undergoes a phase transition of the third order separating two phases (pulled and pushed). We prove this result (i) for the eigenvalues of one-cut, off-critical random matrices (log-gas in dimension \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d=1$$\end{document}d=1) with hard walls; (ii) in arbitrary dimension \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d\ge 1$$\end{document}d≥1 for a gas with Yukawa interaction (aka screened Coulomb gas) in a generic confining potential. The latter class includes systems with Coulomb (long range) and delta (zero range) repulsion as limiting cases. In both cases, we obtain an exact formula for the free energy of the constrained gas which explicitly exhibits a jump in the third derivative, and we identify the ‘electrostatic pressure’ as the order parameter of the transition. Part of these results were announced in Cunden et al. (J Phys A 51:35LT01, 2018).


Introduction and Statement of Results
Phase transitions-points in the parameter space which are singularities in the free energygenerically occur in the study of ensembles of random matrices, as the parameters in the joint probability distribution of the eigenvalues are varied [9]. The aim of this paper is to characterise the pulled-to-pushed phase transition (defined later) in random matrices with hard walls and, more generally, in systems with repulsive interaction in arbitrary dimensions.
B Fabio Deelan Cunden fabio.cunden@ucd.ie 1 Given an interaction kernel Φ : R d → (−∞, +∞] and a potential V : R d → R, we define the energy associated to a gas of N particles at position x i ∈ R d (i = 1, . . . , N ) as We regard −∇Φ(x − y) as the force that a particle at x exerts on a particle at y, and V (x) as a global coercive potential energy, V (x) → +∞ as |x| → ∞. The typical interactions we have in mind are repulsive at all distances, i.e. −∇Φ(x) · x ≥ 0. It is natural to assume that the repulsion is isotropic Φ(x) = ϕ(|x|) and that the potential is radial V (x) = v(|x|).
The normalisation of the energy is done in such a way that both terms (the sum over pairs and the sum of one-body terms) are of same order O(N 2 ) for large N . Indeed, in terms of the 'granular' normalised particle density, the energy (1.1) reads The minimisers ρ N of the discrete energy should achieve the most stable balance between the repulsive effect of the interaction term and the global confinement. Finding global and constrained minimisers of the discrete energy E N is a question of major interest in the theory of optimal point configurations.
For a large class of interaction kernels, the sequence of minimisers ρ N converges toward some non-discrete measure ρ when N → ∞. It is therefore convenient to reframe the optimisation problem in terms of a field functional E : P(R d ) → (−∞, +∞] defined on the set of probability measures ρ ∈ P(R d ) by (1.4) where ρ(x) represents the normalised density of particles around the position x ∈ R d . Mean field energy functionals of the form (1.4) and their minimisers have received attention in the study of asymptotics of the partition functions of interacting particle systems. Consider the positional partition function of a particle system defined by the energy (1.1) at inverse temperature β > 0 (1.5) For several particle systems including eigenvalues of random matrices, Coulomb and Riesz gases, the leading term in the asymptotics of the free energy is the minimum of the mean field energy functional [7,44,56] − 1 β N 2 log Z N → min as N → ∞. (1.6) Remark 1 We stress that physically the change of picture {x i } → ρ(x) is accompanied by an entropic contribution of the gas (the 'number' of microstates {x i } of the gas that contribute to a given macroscopic density profile ρ(x)). However, the scaling in N of the energy is such that the entropic contribution is always sub-leading in the large N limit, and the free energy is dominated by the internal energy component. To see this, we remark that the Boltzmann factor corresponding to the energy (1.1) can be written as (1.7) Therefore, the mean-field limit N → ∞ is simultaneously a limit of large number of particles and a zero-temperature limit, with an energy O(N ) extensive in the number of particles as in standard statistical mechanics. (We learned this argument from a paper by Kiessling and Spohn [40].) One then expects that in the large-N limit, the free energy of N particles at inverse mean-field temperature β M F = β N approaches the minimum energy since at zero temperature the entropy is absent. This also explains the rescaling in (1.6), as β M F N = β N 2 .
Having clarified the physical meaning of the mean-field functional (1.4), we now turn to the problem addressed in this paper.
Notational remark Throughout the paper P(B) denotes the set of probability measures whose support lies in B ⊂ R d . The Euclidean ball of radius R centred at 0 is denoted by B R = {x ∈ R d : |x| = (x 2 1 + · · · + x 2 d ) 1/2 ≤ R}. If F is a formula, then 1 F is the indicator of the set defined by the formula F. We also use the notation a ∧ b = min{a, b}.

Formulation of the Problem
Consider a particle systems with Boltzmann factor as in (1.7) and assume that, for large N , the partition function behaves like (1.6). Under general hypotheses, the balance between mutual repulsion and external confinement allows for the existence of a compactly supported global minimiser of the energy functional E in (1.4). In radially symmetric systems (Φ(x) = ϕ(|x|) and V (x) = v(|x|)), the minimiser ρ R is supported on a ball (1. 8) We say that the gas is at equilibrium in a ball of radius R with density ρ R . Suppose that we want to compute the probability that the gas is contained in a volume B R with R = R , i.e. (1.9) The denominator is nothing but the partition function of the gas Z N (∞) = Z N , while the integral in the numerator is the partition function Z N (R) of the same gas constrained to stay in the ball B R . In view of the asymptotics (1.6), for large N where ρ R is the equilibrium measure of the gas confined in B R  Fig. 1 The pulled-to-pushed transition for a log-gas in dimension d = 1 in a quadratic potential (GUE) (Note that ρ R = ρ R for all R ≥ R .) We conclude that the probability (1.9) decays as where the large deviation function F(R) is (1.13) Clearly, F(R) ≥ 0 and is non-increasing. The physical interpretation of F(R) and ρ R is clear: F(R) is the excess free energy of the gas, constrained within B R , with respect to the situation where it occupies the unperturbed volume B R ; the measure ρ R describes the equilibrium density of the constrained gas in the limit of large number of particles N . The general picture is as follows (see Fig. 1): (i) In the unconstrained problem (B R = R d ) the global minimiser ρ R is supported on the ball B R ; (ii) If R > R , the constraint in (1.13) is immaterial (B R contains B R ), and hence the equilibrium measure is ρ R = ρ R and F(R) = 0. This is the so-called pulled phase, borrowing a terminology suggested in [48]; (iii) If R < R the system is in a pushed phase, the constraint is effective, and the equilibrium energy of the system increases E [ρ R ] ≥ E [ρ R ]. (iv) At R = R the gas undergoes a phase transition and the free energy F(R) displays a non-analytic behaviour. Typically at microscopic scales one expects a crossover function separating the pushed and pulled phases.
The goal of this work is to investigate the properties of the excess free energy at the critical point R for certain systems with pairwise repulsive interactions. Explicitly solvable models related to random matrices suggest that in the vicinity of the critical point 15) implying that the transition between the pushed and pulled phases of the gas is third-order.
Here we demonstrate that (1.15) is generically true for a large class of systems with repulsive interactions. More precisely, we prove (1.15): (i) For the log-gas Φ(x) = − log |x| in dimension d = 1 (eigenvalues of one-cut, offcritical random matrices with hard walls); a (1.16) in arbitrary dimension d ≥ 1, including its limiting cases m → 0 (Coulomb gas) and a → 0 (Thomas-Fermi gas). In (1.16), K ν denotes the modified Bessel function of the second kind.
The precise assumptions on the confining potential V (x) are presented together with the statements of Theorems 2 and 5 below.

Hermitian Random Matrices
Many random matrix theory (RMT) phenomena have been first discovered for invariant measures on the space of N × N complex Hermitian matrices M of the form where V is a scalar function referred to as the potential of the matrix model and β = 2. Expectations of conjugation-invariant random variables with respect to these measures can be reduced, via the Weyl denominator formula, to an integration against the joint density of the eigenvalues x 1 , . . . , x N , which has the form This energy is of the form (1.1) in dimension d = 1 and with Φ(x) = − log |x|. Hence, Z N =´R e −β E N dx 1 · · · dx N can be interpreted as the partition function of a log-gas (or a 2d Coulomb gas) of N repelling particles on the line in a global confining potential V . This physical interpretation also suggests that we may consider generic values of the Dyson index β > 0, interpreted as inverse temperature. In the large-N limit the eigenvalue empirical measure ρ N = 1 N i δ x i weakly converges to a deterministic density. This limit is the equilibrium measure (the minimiser) of the functional  [18,25], 20) where F β (t) is known as the β-Tracy-Widom distribution [59] and can be expressed in terms of the Hastings-McLeod solution of the Painlevé II equation. The macroscopic (atypical) fluctuations of max |x i | are instead described by a large deviation function. More precisely, for all β > 0 the following limit exists For the log-gas in d = 1, the asymptotics of the partition function has been rigorously established in several works. When R ≥ R , the equilibrium density is supported on the ball of finite radius R , so that ρ R = ρ R (the volume constraint is ineffective). Hence, the large deviation function is the excess free energy (1.14) of the gas of eigenvalues forced to stay between two hard walls at ±R. It is clear that F(R) = 0 for R ≥ R (pulled phase), while F(R) ≥ 0 for R < R (pushed phase) when the log-gas gets pushed by the hard walls at ±R. The calculation of F(R) for the GUE and its β > 0 extensions was performed in detail by Dean and Majumdar [19,20] who found explicit expressions for the density and for the excess free energy (1.24) A closer inspection of the latter formula provides a thermodynamical characterisation of the pulled-to-pushed transition. Indeed, we see that as R → R . Therefore, the third derivative of the free energy of the log-gas at the critical point R = √ 2 is discontinuous. Similar phase transitions of the pulled-to-pushed type have been observed in several physics models related to random matrices [13,47], including large-N gauge theories [3,36,57,64], longest increasing subsequences of random permutations [39], quantum transport fluctuations in mesoscopic conductors [12,33,34,62,63], non-intersecting Brownian motions [28,58], entanglement measures in a bipartite system [17,26,27,49], random tilings [10,11], random landscapes [29], and the tail analysis in the KPZ problem [41]. (See also the recent popular science articles [5,65]).
An explanation of the critical exponent '3' has been put forward by Majumdar and Schehr [47] (see also [2]) based on a standard extreme value statistics criterion and a matching argument of the large deviation function behaviour in the vicinity of the critical value R and the left tail of the Tracy-Widom distribution [52] The criterion predicts that if the equilibrium density of a log-gas in the pulled phase vanishes as ρ R (x) ∼ R 2 − x 2 at the edges-the so-called off-critical case-then the pulled-topushed phase transition is of the third order. This conjectural relation between the particular behaviour of the gas density and the arising non-analyticities in the free energies has been verified in several examples, even though each particular case (i.e. each matrix ensemble defined by a potential V ) requires working out explicitly the model-dependent F(R) to compute the critical exponent.
In this paper we derive a general explicit formula for the free energy F(R) of a log-gas in dimension d = 1 in presence of hard walls, and we prove the universality of the third-order phase transition for one-cut, off-critical matrix models. This proves the prediction arising from the extreme value statistics criterion formulated in [47].
While here we examine problems with radial symmetry in both the potential and the hard walls as they constitute a paradigmatic framework and allow for a systematic treatment, the electrostatic interpretation we present in Sect. 1.4.2 below enjoys a wider range of application. In Remark 5, indeed, we will show how the formulae and conclusions concerning the symmetric case carry over without extra efforts to non-symmetric potentials or random matrices with a single hard-wall as well.
The constrained minimisation problem for the log-gas in d = 1 is usually solved by using complex-analytic methods or Tricomi's formula, which are well-known to those working in potential theory and random matrices. Nevertheless, we found a particularly convenient (and perhaps not so well-known) method based on the decomposition into Chebyshev polynomials which is well-suited to this class of problems. At the heart of the method is the following pointwise multipole expansion of the two-dimensional Coulomb interaction (proven in Appendix B) where the T n 's are the Chebyshev polynomials of the first kind. They are defined by the orthogonality relation and they form a complete basis of L 2 ([−1, 1]) (with respect to the arcsine measure). The above identity was recently used and discussed in [30,32] (the authors refer to some unpublished lecture notes by U. Haagerup). We consider potentials V (x) satisfying the following assumptions.
We remark that strictly convex and super-logarithmic V (x)'s are in the class of one-cut, off-critical potentials.
We are going to present and discuss two theorems for the log-gas that will be proven in Sect. 3.

Theorem 1
In dimension d = 1, let Φ(x) = − log |x|, and V (x) be a potential satisfying Assumption 1. Then (i) there exists a unique probability measure that is solution of the constrained minimisation problem (1.11) for the energy functional (1.4), and it takes the form Then, the equilibrium measure (1.29) is uniquely determined as The critical radius R is the smallest positive solution of the equation In other words, Eq. (1.29) shows that in the pulled phase the equilibrium density ρ R (x) = ρ R (x) is supported on a single interval on the real line (one-cut property), is strictly positive in the bulk, and vanishes as a square root at the edges ±R (off-critical case). Moreover, since Q(R ) > 0 one gets a fact that will be used later on.
On the other hand in the pushed phase, the density is strictly positive in its support and has an integrable singularity at the hard walls ±R (cf. with the GUE case (1.23)). See Fig. 1.
A corollary of formulae (1.29)-(1.31) is the first main result of this paper.

Theorem 2
With the assumptions of Theorem 1, the excess free energy (1.14) of the log-gas is Moreover, F(R) displays the following non-analytic behaviour at R :

Remark 2
In the statement of the results, the family of strictly convex potentials is not the largest class of potentials where the above result holds true. What is really required is that in the pulled phase the associated equilibrium measure is supported on a single interval, and vanishes as a square root at the endpoints. (An explicit characterisation of these conditions is quite complicated.) More precisely, Theorem 1 is true when V (x) is a one-cut potential. 'One-cut' means that the equilibrium measure is supported on a single bounded interval. Theorem 2 is valid under the hypotheses that V (x) is one-cut and off-critical potential. 'Off- This is certainly true for strictly convex potentials. For 'critical' potentials the pulled-topushed transition is weaker than third-order (see Eq. (3.18) in the proof). These potentials are however 'exceptional' in the 'one-cut' class. The problem for 'multi-cut' matrix models remains open.

Example 1
We reconsider the free energy of the GUE (the log-gas on the real line in a quadratic potential V (x) = x 2 /2). There are only two nonzero Chebyshev coefficients (1.30) From (1.32) we read that the critical radius R is the positive solution of R 2 /2 = 1, i.e., R = √ 2. From (1.31), in the pushed phase P R (x) = (2 + R 2 − 2x 2 )/2, so that P r (r ) = (2 − r 2 )/2. Applying the general formula (1.34), one easily computes the large deviation function as for R ≤ √ 2, and zero otherwise. This coincides with the known result (1.24).
Theorem 2 and the previous discussion might lead to conclude that the universality of the third-order phase transition is inextricably related to the presence of a Tracy-Widom distribution separating the pushed and pulled phases. A hint that this is not the case comes from the study of extreme statistics of non-Hermitian matrices whose eigenvalues have density in the complex plane.

Non-Hermitian Random Matrices
Consider a log-gas in dimension d = 2 (the 'true' 2d Coulomb gas in the plane) The gas density for large N converges to the equilibrium measure of the energy functional (1.18) extended to measures supported on the complex plane. When V (x) = |x| 2 /2 and β = 2, with the identification C R 2 , the model corresponds to the density of eigenvalues of the complex Ginibre (GinUE) ensemble of random matrices, a non-Hermitian relative of the GUE. The equilibrium measure in this case is uniform in the unit disk (circular law) with R = 1, but the typical fluctuations of the extreme statistics max |x i | are not in the Tracy-Widom universality class. More precisely, setting γ N = log N − 2 log log N − log 2π, Rider [55] proved that where the limit is the Gumbel distribution G(t) = exp(− exp(−t)). This result is universal [6] for the log-gas in the plane at inverse temperature β = 2. The atypical fluctuations are described by a large deviation function F GinUE (R) which can be computed by solving the constrained variational problem of a log-gas in the plane. For the Ginibre ensemble, this was computed by Cunden et al. [14] who found 1 Of course, the explicit form of the large deviation function is specific to the model (for instance F GUE (R) = F GinUE (R)). The surprising fact is that even in dimension d = 2, the pushed-to-pulled transition of the log-gas is of the third order We remark that, in this case, the order of the phase transition is not predicted by the classical 'matching argument'. In fact, for the log-gas in d = 2, the matching between the typical fluctuations (Gumbel) and the large deviations is more subtle due to the presence of an intermediate regime, as found recently in [42,43]. This suggests that the critical exponent '3' is shared by systems with repulsive interaction whose microscopic statistics belongs to different universality classes. 2 For which interactions can Theorem 2 be extended?

Beyond Random Matrices: Yukawa Interaction in Arbitrary Dimension
The logarithmic interaction Φ(x − y) = − log |x − y| is the Coulomb potential in dimension d = 2, namely the Green's function of the Laplacian on the plane Based on this observation, in [15] we put forward the idea of considering energy functionals whose interaction potential Φ is the Green's function of some differential operator D:

Coulomb and Thomas-Fermi Gas
When D = −Δ in R d , the system corresponds to a d-dimensional Coulomb gas (a system with long range interaction) The constrained variational problem for a d-dimensional Coulomb gas (in R d ) can be solved using macroscopic electrostatic considerations. Another solvable model is the gas with delta potential corresponding to D = I (the identity operator) in R d . In this case, the repulsive kernel is proportional to a delta function (zero range interaction) 49) and the energy is an energy functional in the Thomas-Fermi class.
with v increasing and strictly convex. In the case of Coulomb interaction, assume also the growing condition lim |x|→∞ Then, the constrained equilibrium measure (1.11) of the energy functional (1.

4) for Coulomb and Thomas-Fermi gases is unique and is given by
(1.51) The critical radius R is the unique positive solution of the equation

52)
while c(R) and μ(R) are given by (1.53) In the pulled phase, the equilibrium density of the Coulomb gas is supported on the ball of radius R and there is no accumulation of charge on the surface (c(R) = 0 for R ≥ R ). In the pushed phase, the equilibrium density in the bulk does not change, while an excess charge (c(R) > 0 for R < R ) accumulates on the surface.
For the Thomas-Fermi gas, R -the edge in the pulled phase-is determined by the condition that the gas density vanishes on the surface, i.e. μ(R ) = v(R ) for R ≥ R . In the pushed phase R < R , the chemical potential increases to keep the normalisation of ρ R , but there is no accumulation of charge on the surface, i.e. singular components in the equilibrium measure (otherwise the energy would diverge).
A direct calculation yields the free energy [15,16] (1.54) From the exact formulae above, one can check that F(R) has a jump in the third-derivative at R = R . Therefore, the critical exponent '3' is shared by systems with long-range (Coulomb) and zero-range (delta) interaction. This suggests that the third-order phase transition is even more universal than originally expected.

Yukawa Gas in Generic Dimension
The ubiquity of this transition calls for a comprehensive theoretical framework, which should be valid irrespective of spatial dimension d and the details of the confining potential V , and for the widest class of repulsive interactions Φ. Fix two positive numbers a, m > 0, and define Φ = Φ d as the solution of Note that this kernel naturally interpolates between the Coulomb electrostatic potential in free space (long-range, for a = 1 and m = 0), and the delta-like interaction (zero-range, for a = 0 and m = 1); intermediate values a, m > 0 correspond to the Yukawa (or screened Coulomb) potential.
with v increasing and strictly convex.
The condition that v(r ) is increasing implies the confinement of the gas whenever m > 0. Indeed, using the asymptotic expansion of K ν (z) for large argument, the following limit holds [50,Eq. 10.25 as long as there is screening (m > 0). By a routine argument, this implies the existence and uniqueness of the minimiser of E in P(R d ). The minimiser ρ R is compactly supported In particular, the support is simply connected. The solution of the constrained equilibrium problem for the Yukawa interaction (announced in [16]) is stated below. 3 In the case of Coulomb interaction, m = 0, assume also the growing condition lim |x|→∞ Then, the constrained equilibrium measure (1.11) of the energy functional (1.4) for the Yukawa gas is unique and is given by (1.60) The critical value R is the smallest positive solution of c(R ) = 0, i.e. the solution of In particular: (i) in the pulled phase the equilibrium measure is absolutely continuous with respect to the Lebesgue measure on R d ; (ii) when the gas is 'pushed', the density in the bulk increases by a constant and a singular component builds up on the surface of the ball B R . Theorem 4 implies the second main result in this paper : the universality of the jump in the third derivative of excess free energy of a Yukawa gas with constrained volume. This universality extends to the limit cases m → 0 (Coulomb gas) and a → 0 (Thomas-Fermi gas) solved in [15,16], respectively. See Remark 3 below.

Remark 3
One can check that, in the limit cases of Coulomb (m = 0) and delta interactions (a → 0), the expression (1.57) for the equilibrium measure simplifies as in (1.51). Indeed, Eq. (1.57) for m = 0 and a = 1 yields for m > 0. Therefore, in the Thomas-Fermi gas there is no condensation, c(R) → 0 as a → 0. When a → 0 and m = 1, from equation (1.60) we recover the value of the chemical potential in the Thomas-Fermi gas (1.53), and Moreover, Eq. (1.61) for the critical radius R reduces to its corresponding expression (1.52) for the Thomas-Fermi gas.
Theorem 5 Under the assumptions of Theorem 4, the excess free energy (1.14) for a Yukawa gas is given by This formula implies that

Remark 4
The formula of the free energy F(R) clearly matches (1.54) for m = 0 and a = 1.
The limit a → 0 can be obtained as above by using the limit (1.63). Indeed from the expression (1.59) of c(R) one gets for m = 1 and R ≤ R and the critical radius is R = 1. The excess free energy can by easily computed using (1.65) for R ≤ 1, and zero otherwise, in agreement with (1.44).

Order Parameter of the Transition
The universal formulae (1.34) and (1.65) for the free energy F(R) are remarkably simple. It feels natural to ask whether it is possible to derive them from a simpler physical argument. Moreover, it would be desirable to express the free energy in terms of a quantity that captures the non-analytic behaviour in the vicinity of the transition. This quantity, traditionally called order parameter of the transition, must be zero in one phase and nonzero in the other phase.
In the following, we outline a heuristic argument that reproduces formulae (1.34) and (1.65), and identifies the 'electrostatic pressure' as the order parameter of the transition.

Electrostatic Pressure: Screened Coulomb Interaction
For the pushed-to-pulled phase transition we can identify the order parameter as follows. Note that the derivative of F(R) is essentially the variation of the free energy with respect to the volume of the system. The volume of the ball is vol( where P = ∂ F ∂vol is the pressure of the gas confined in B R . The increase in free energy of the constrained gas should match the work W R →R done in a compression of the gas from the initial volume vol i = vol(B R ) to the final volume vol f = vol(B R ), the system being in equilibrium with density ρ r at each intermediate stage R ≤ r ≤ R . In formulae, where P = p(r ) is the pressure on the gas confined in B r . In other words, p(r )Ω d r d−1 dr is the elementary work done on the surface of the ball of radius r being compressed from r + dr to r . We now show that the 'electrostatic' pressure for a Yukawa gas is quadratic in the excess charge on the surface in the pushed phase, and zero in the pulled phase.
How to compute the pressure exerted by the surface's field on itself? The argument that follows is similar to the one used to evaluate the 'electrostatic pressure' on a layer of charges, e.g. the surface of a charged conductor. This is a problem that some textbooks in electrostatics occasionally mention (see the classical books of Jackson [38, Section 3.13] and Purcell [51, Section 1.14]), but which is rarely discussed due to the difficulty of making the argument rigorous for conductors of generic shapes. To stress the analogy and for the lack of a better terminology, in the following we will keep using the expressions 'electrostatic' pressure, force, field, Gauss's law, etc. even though the interaction we are considering is Yukawa (screened Coulomb).
When the system is confined in a ball B r at density ρ r , the pressure is given by the normal force per unit area, where ΔA = r d−1 ΔΩ is a small area on the sphere of radius r . An intuitive guess for the force ΔF n experienced by the surface element ΔA during the compression is  Fig. 2 The electrostatic pressure on the spere. Left: consider removing a small disk of area ΔA from the surface of the sphere. The electric field experienced by a point in ΔA on the sphere corresponds to the field on the hole generated by the sphere with the small disk removed (in the limit where ΔA → 0). Right: the electric field on the hole can be computed from the total electric field generated by the sphere and the field generated by the disk using the superposition principle. Note that in the bulk the electric field must be zero at equilibrium. Therefore E hole must be equal in magnitude to E disk and directed outward This is almost right, but it contains a serious flaw. Indeed, note that the electric field across ΔA includes contributions from the total amount of charge on the sphere-thus including the charge that is being acted upon by the field! We must therefore be 'over-counting', as the charge on ΔA cannot act upon itself. To fix the over-counting, we may imagine to open a small hole corresponding to ΔA and compute the electrostatic field produced by the charge distribution ρ r minus the hole. A glance at Fig. 2 may be helpful. The electrostatic field in the hole is the field experienced by the charges in ΔA. Hence, the force acting on ΔA is ΔF n = (charge in ΔA) × (electrostatic field in the hole). (1.76) When the gas density is ρ r , the amount of charge in ΔA is clearly given by (1.77) as the charge in the bulk is a continuous distribution and does not contribute. There is a clever argument to compute the electric field in the hole: we know that the field in the hole plus the field in the disk gives the field produced by the sphere just inside the sphere.
From basic considerations, it is clear that the field E bulk is zero inside the ball B r (a consequence of electrostatic equilibrium) and perpendicular to the surface immediately outside the ball (a necessary condition for electrostatic equilibrium). Therefore E hole = −E disk .
What is E disk ? Denote by Φ r (x) the potential generated at x ∈ R d by the charge in ΔA. The electrostatic field produced by this tiny amount of charge is −∇Φ r (x). Fortunately, we only need to know this in the vicinity of the disk E disk = −∇Φ r (x) with |x| = r , where the disk can be approximated by a planar surface with uniform density c(r ) ΔA Ω d r d−1 , see Fig. 2. Consider a small cylinder C(h, ΔA) of height 2h and base ΔA cutting across the surface of the ball of radius r as in Fig. 3. To compute the field on the surface, one should integrate the screened Poisson equation (1.78) If we integrate (1.78) over the small cylinder (see Fig. 3 (1.79) By symmetry, the electrostatic field −∇Φ r is perpendicular to the bases, alongx = x/r , and, using the divergence theorem, (1.80) Letting h → 0, the volume integral on the left hand side vanishes and we get immediately inside the disk.
(1.81) Therefore, we conclude that As a byproduct, we see instead that E surface = c(r ) a 2 r d x; this is the familiar statement that, at equilibrium, the electrostatic field generated by a charged conductor immediately outside is perpendicular to its surface and proportional to the charge density.

Electrostatic Pressure: Random Matrices
The argument outlined above can be repeated almost verbatim for the log-gas on the line (eigenvalues of random matrices). However there is a twist in the computation. Again, by conservation of energy, the increase in free energy must match the work W R →R done in a compression of the gas from the initial volume (length) vol i = vol(B R ) = 2R to the final volume vol f = vol(B R ) = 2R, with the system in equilibrium with density ρ r at each intermediate stage R ≤ r ≤ R . In formulae, where 2 p(r )dr is the elementary work done in an infinitesimal compression (the factor 2 comes from axial symmetry).
When the system is confined in a ball B r at density ρ r , the pressure is given by the normal force per unit length. The force F n (x) at point x is equal to the charge ρ r (x)dx in the infinitesimal segment dx around x times the electric field. To proceed in the computation it is convenient to use complex coordinates (recall that − log |x| is the Coulomb interaction in dimension d = 2).
The electric field generated by ρ r (y)dy at z ∈ C \ [−r , r ] is the Stieltjes transform Note that ρ r (z) = 0 when z / ∈ [−r , r ]. Using Gauss's theorem (Plemelj formula), when z approaches the real axis, the field generated by ρ r is when |x| < r , where ffl denotes Cauchy's principal value. Note that, at equilibrium, the real part of G(x ± i ) cancels with the field −V (x) generated by the external potential V (x) (the net tangential field must be zero). Therefore, the electric field experienced by a point in the vicinity of x is (1.89) The force (= electric field× charge) per unit length is In the vicinity of x it becomes (1.91) Note, in particular, that the tangential force experienced by a point x inside the conductor is zero (as it should be at equilibrium). At the edge x = r (similar considerations for x = −r ) the situation is more delicate. By symmetry, the total field E edge generated by ρ r (x) at x = r must be directed along the x-axis, but must be zero for x < r − . Therefore, repeating the argument of the previous section, the field experienced by the 'hole' at the edge x = r is half the field generated by ρ r (x) at r + , i.e. E hole = (1/2)E edge .
To compute the pressure, we look at the edge x = r , we sum the forces on a small circular contour of radius centred at r , and then we take the limit → 0 (see Fig. 4) (1.92) where the factor 1/2 comes from the fact that E hole = (1/2)E edge . Remembering that the force is eventually zero on the semicircular part of the contour in the bulk, the integral is given by πi times the residue at z = r : Inserting this formula in (1.85) we indeed recover (1.34).

Remark 5
These electrostatic considerations indicate a route to compute large deviation functions for extreme eigenvalues of random matrices more general that those fulfilling Assumption 1. For instance, one may ask whether it is possible to obtain an electrostatic formula for the large deviation of the top eigenvalue x max (i.e. the rightmost particle) of a random matrix from a β-ensemble. While this question does not fit into the symmetric setting considered so far, it can be nevertheless easily answered within the electrostatic framework developed above. For a one-cut matrix model with typical top eigenvalue equal to b , the rate function function F(b) in the large deviation decay where the electrostatic 'pressure' experienced by the gas is now given by is the density of the log-gas (with support in [a, b]) constrained to stay on the left of the hard wall at x = b. One can easily check the validity of this formula for a few cases already considered in previous works. For instance, the constrained density of the GUE ensemble with one hard wall is [19,20] In this case Computing the residue (1.100) The rate function for b ≤ 4 is

Remark 6
The calculation of the work done in a compression of a log-gas in dimension d = 1 bears a strong resemblance to a way to calculate the work (energy) per unit fracture length for a crack propagating in a continuous medium. In linear elastic fracture mechanics, Cherepanov [8] and Rice [54] have independently developed a line integral called the Jintegral that is contour independent. The usefulness of this integral comes about when the contour encloses the crack-tip region, as this is where the most intense (actually divergent) fields are found (c.f. the edge of the cut in the log-gas). Evaluating the J -integral then gives the variation of elastic energy. An analogous integral in electrostatics has been discovered later [31]. It is likely that the computations outlined above can be recast in the language of linear elastic mechanics/electrostatics. The link between 'eigenvalues of random matrices' and the 'theory of fractures' suggested here will be explored in future works.

Outline of the Paper
The rest of the article is organised as follows. In Sect. 2 we recall some general variational arguments for the solution of the constrained equilibrium problem. Section 3 contains the proof of Theorems 1 and 2 for the log-gas. In Sect. 4 we present the proof of Theorem 4 (equilibrium problem for Yukawa interaction) and Theorem 5 (universality of the jump in the third derivative of the free energy). Finally, the Appendices A and B contain the proof of some technical lemmas.

Variational Approach to the Constrained Equilibrium Problem
We resort to a variational argument to derive necessary and sufficient conditions for ρ R to be the minimiser of the energy functional E over P(B R ). (These arguments are not new at all. They appear in many different forms and specialisations in the literature, see, e.g., [4,21,46].) Denote by ρ R ∈ P(B R ) a local equilibrium and let Here σ is a (small) perturbation of zero mass. Of course, the perturbation σ must be nonnegative on (supp ρ R ) c , the complement of supp ρ R .
The functional E is quadratic in ρ, hence where E 1 and E 2 are the first and second variations, respectively. They are given explicitly by 3) A sufficient condition for ρ R to be the global minimiser in P(B R ) is that E 1 [ρ R , σ ] ≥ 0 and E 2 [σ, σ ] > 0 for all perturbations σ . Consider the first variation (2.3) and suppose that σ varies among the perturbations whose support lies in supp ρ R . Because σ is arbitrary and zero mass, for the first variation E 1 to vanish, it must be true that for some constant μ(R). Consider now perturbations σ with support in B R . Remembering that σ ≥ 0 in (supp ρ R ) c , we see that a sufficient condition for E 1 ≥ 0 is that The conditions (2.5)-(2.6) are known as Euler-Lagrange (E-L) conditions and can be summarised by saying that if ρ R is a minimiser of E in P(B R ) then there exists a constant (2.7) The constant μ(R) is called chemical potential. In general, the E-L conditions provide the saddle-point(s) of the energy functional. To show that ρ R is actually the minimiser it remains to check that the E is strictly convex, i.e. that the second variation E 2 > 0. Denote by σ the Fourier transform of σ . Then, For the class of interactions considered in this paper, Φ(k) > 0 (see Appendix A), and this ensures that E 2 [σ, σ ] > 0. Therefore, for all R > 0, the solution of the E-L conditions is the unique minimiser of the energy functional E in P(B R ).

Proof of Theorems 1 and 2
Part i) of Theorem 1 is a classical result in potential theory, see [21]. There are two cases: -if the walls are not active (pulled phase) the density is supported on supp and the density is given by Tricomi's formula [60] ρ and is given by (3.2) with the replacement R → R.
The new content of the theorem is part ii). To prove it, we expand the potential V and the regular part of the density into Chebyschev polynomials where a n (R) A priori, the above expansions are in L 2 ([−1, 1]). In fact, V ∈ C 3 implies that c n (R) = O(n −3 ) so that the series n≥0 c n (R)T n (u) and its derivative are pointwise convergent almost everywhere to V and V , respectively. We will see in the course of the proof that the absolute convergence of n nc n (R) implies the pointwise convergence of n a n (R)T n (u), too. Note also that c n = 0 if n is odd (the potential V (x) is symmetric by Assumption 1). To proceed we use the following identity.
Lemma 1 Let n ≥ 0 be an even integer. Then, uT n (u) = nT 0 (u) + nT n (u) + 2n (T 2 (u) + T 4 (u) + · · · + T n−2 (u)) . (3.5) We first express the equation (3.1) for the critical radius R in terms of the c n 's. After the change of variable x = R u, (3.1) becomes where we used Lemma 1 and the orthogonality relation (1.28). Note that nc n (R ) = O(n −2 ) and hence the series is absolutely convergent. This proves that R is the solution of (1.32). The Chebyshev polynomials satisfy the following electrostatic formula (for a proof see Appendix B).
From the E-L equation, an application of the Chebyshev electrostatic formula (3.6) gives This equation and the normalisation of ρ R (x) imply that a 0 (R) = 1, 1 n a n (R) + c n (R) = 0 ∀n ≥ 1. (3.7) In particular, we have an explicit formula for the chemical potential The sequence nc n (R) is O(n −2 ) and hence the series is pointwise convergent almost everywhere. This concludes the proof of Theorem 1.
To prove Theorem 2 we begin by computing F(R): (3.10) (The first identity follows from the definition of F(R) as energy difference (1.18), and the E-L condition; then we expanded in Chebyshev polynomials and used their orthogonality relation; the last equality follows from (3.8) and (3.7).) Therefore We want to prove that the above expression is equal to −P R (R) 2 /(2R). First, notice that (3.12) Comparing with (3.11), the identity to show to complete the proof is (3.14) Therefore we find This concludes the proof of the integral formula (1.34).
We proceed to prove that F(R) has a jump in the third derivative. First, note that F(R) is identically zero for R ≥ R , while F(R) ≥ 0 for R ≤ R . From (1.34) and the fact that P R (R ) = 0, one sees that On the other hand, Indeed, by Assumption 2 on the potential V (x), it is easy to check that P R (R ) < 0 strictly. .

Proof of Theorems and 5
A systematic analysis of the equilibrium problem for the screened Coulomb interaction does not seem to have appeared in the existing literature. Here we solve the problem (Theorem 4). In this Section we put forward a sensible ansatz for ρ R depending on two parameters (chemical potential μ and surface charge c), that we then prove to be the minimiser by imposing the E-L conditions that fix μ = μ(R) and c = c(R). The only thing that needs to be checked is the positivity of the candidate solution (Remark 7). In solving the problem, we will make the most of its spherical symmetry. A key technical ingredient in the derivation will be a shell integration lemma (Lemma 3), the analogue of Newton's formula (4.20) for a Yukawa interaction in generic dimension d ≥ 1.

General Form of the Constrained Minimisers
The usual strategy in these minimisation problems is to look for a candidate solution of the E-L conditions (2.7). The condition E 2 > 0 guarantees that the saddle-point is the minimiser in P(B R ).
In absence of volume constraints, i.e. B R = R d , the global minimiser is supported on a ball of radius R . For R > R , the density of the gas does not feel the hard walls and ρ R = ρ R . We anticipate here that for R < R , supp ρ R = B R (thus explaining the name 'pushed phase') so that the second E-L condition (2.7) is immaterial. A first idea to find ρ R is to use the property DΦ d (x) = Ω d δ(x). Applying D to both sides of the first E-L condition we formally get dρ R (x) However the above ansatz is, in general, incorrect. In particular, the equality (4.1) is not true if ρ R contains singular components.
One can prove the absence, at equilibrium, of condensation of particles in the bulk (absence of δ-components). To see this, one compares the energy of a density containing a δ-function in the bulk to one where the δ-function has been replaced by a narrow, symmetric mollification δ (see [4, Section 3.2.1] for details). This argument fails on the boundary where it is not possible to consider symmetric mollifications δ contained in the support.
As a matter of fact, condensation of particles, though not possible in the bulk, can and do occur on the boundary of the support! For Coulomb gases (D = −Δ) this statement is the well-known fact that, at electrostatic equilibrium, any excess of charge must be distributed on the surface of a conductor. In a rotationally symmetric problem, any accumulation of charge must be uniformly distributed on the 'surface', i.e. the boundary of the support of ρ R . Using the same argument, one can argue that in the unconstrained problem-pulled phasethe accumulation of charge on the surface is not possible (otherwise it would be possible to mollify the δ-components on the surface, lowering the energy). Therefore, a more appropriate ansatz is for some constants μ(R) and c(R) to be determined. In the following we will show that, for all R, there exists a unique choice of μ(R) and c(R) such that (4.2) is a saddle-point.

Remark 7 In spherical coordinates
where we split the contribution on the surface c(R) and the density f (r ) in the bulk The equilibrium measure ρ R is continuous in the bulk |x| < R, and contains, in general, a singular component on the surface |x| = R. From the explicit formulae (1.59)-(1.60) for c(R) and μ(R), it is not yet obvious to see that ρ R is a positive measure, nor that the critical radius R is nonzero. We show here that the equilibrium measure is both positive in the bulk and on the surface.
Consider first the pushed phase R < R . Starting from the surface, note that c(R) is continuous for R ≥ 0 and differentiable at least twice for R > 0 and R = R . Moreover, and is nondecreasing (for R > 0). Since v(r ) and v (r )r d−1 are both strictly increasing, we conclude that c(0) = 1, and c(R) is strictly decreasing. Therefore, there exists a positive radius R > 0 such that c(R ) = 0. The pulled phase is characterised by the absence of condensation of charge on the surface: c(R) = 0 for R ≥ R . Hence the singular component on the surface is nonnegative for all R ≥ 0. We now analyse the bulk. This amounts to showing that f (r ) in (4.4) is nonnegative. First, note that a 2 (r d−1 v (r )) is positive by assumption. Since v(r ) is strictly increasing, it is enough to show that μ(R) − v(R) ≥ 0 to imply that μ(R) − v(r ) ≥ 0 for all r ≤ R. In the pushed phase (4.6) The ratio ϕ d (R) ϕ d (R) is negative while c(R) and v (R) are both positive. Therefore, f (r ) ≥ 0 for r ≤ R ∧ R .
Before embarking in the calculations leading to the explicit formulae for μ(R) and c(R), we derive the universal formula of the free energy (Theorem 5) assuming Theorem 4.

Proof of Theorem 5
First, note that for R ≥ R , For R < R , we compute E [ρ R ] by inserting the explicit form of the minimiser ρ R in the functional. This calculation is made possible by using the explicit formulae (1.58)-(1.60) for ρ R and, crucially, a shell theorem for Yukawa interaction (see Lemma 3 below). Eventually, we find for R < R : Combining (4.7) and (4.8) we obtain the claimed formula (1.65).
To prove the jump in the third derivative, we remark again that c(R) is continuous for R ≥ 0 and differentiable at least twice for R > 0 and R = R . Moreover, R > 0, and c(R ) = 0. From the explicit formula (1.65): On the other hand, since c(R) is strictly decreasing in the pushed phase.

Chemical Potential and Excess Charge
The chemical potential μ(R) and the excess charge c(R) are fixed by the normalisation of ρ R and the E-L conditions. Assume R ≤ R (note that at this stage R is not known). We show here that in the pushed phase the chemical potential and the excess charge are solution of the following linear system therefore the solution of the linear system (4.13) is unique. The normalisation condition iŝ The E-L condition in the support of ρ R iŝ At this stage, we need an analogue of the electrostatic shell theorem to perform the angular integration in (4.16). For Yukawa interaction the potential outside a uniformly 'charged' sphere is the same as that generated by a point charge at the centre of the sphere with a dressed charge (dependent on the radius of the sphere); inside the spherical shell the potential is not constant. The precise statement is a 'shell theorem' for screened Coulomb interaction. where we set z = |x| ≤ R.
Differentiating with respect to z, (4.22) Using integration by parts and the properties of ψ d and ϕ d we find Therefore, Eq. (4.22) reads which can be simplified using the following lemma. Using the above lemma, Eq. (4.22) simplifies as The above equality must be true for all z. Therefore, the E-L equation can be eventually written as c(R) (4.28) The above relation between μ(R) and c(R) is the first equation in (4.13). We finally get the explicit formulae Here the gas is confined in a quadratic potential v(r ) = r 2 /2, and a = m = 1/2. Bottom: energy density (on a convenient scale) corresponding to the equilibrium measures. The energy density is constant, and equal to the chemical potential μ, in the support of the equilibrium measure . (4.30) At this stage, what would remain to do is just checking the positivity of ρ R . We did this in Remark 7. This concludes the proof of the theorem (Fig. 5).
using the standard trigonometric addition formula. After simplifications The last integral can be evaluated with a partial fraction expansion and we find (for t > 1 and −1 < z < 1) ∂ ∂t n≥1 J n (t) z n = − 1 t Comparing with the right hand-side of (B.10), we conclude that n≥1 J n (t) z n = log 1 − z t + k(z). (B.17) To fix the constant k(z), we may appeal to the large-t behaviour from the right-hand side of (B.15) lim Therefore, k(z) = 0 and the proof of (B.7) is complete.