Large deviations and phase transitions in spectral linear statistics of Gaussian random matrices

We evaluate, in the large-N limit, the complete probability distribution P(A,m) of the values A of the sum ∑i=1N|λi|m , where λ i ( i=1,2,…,N ) are the eigenvalues of a Gaussian random matrix, and m is a positive real number. Combining the Coulomb gas method with numerical simulations using a matrix variant of the Wang–Landau algorithm, we found that, in the limit of N→∞ , the rate function of P(A,m) exhibits phase transitions of different characters. The phase diagram of the system on the (A, m) plane is surprisingly rich, as it includes three regions: (i) a region with a single-interval support of the optimal spectrum of eigenvalues, (ii) a region emerging for m < 2 where the optimal spectrum splits into two separate intervals, and (iii) a region emerging for m > 2 where the maximum or minimum eigenvalue ‘evaporates’ from the rest of eigenvalues and dominates the statistics of A. The phase transition between regions (i) and (iii) is of second order. Analytical arguments and numerical simulations strongly suggest that the phase transition between regions (i) and (ii) is of (in general) fractional order p=1+1/|m−1| , where 0<m<2 . The transition becomes of infinite order in the special case of m = 1, where we provide a more complete analytical and numerical description. Remarkably, the transition between regions (i) and (ii) for m⩽1 and the transition between regions (i) and (iii) for m > 2 occur at the ground state of the Coulomb gas which corresponds to the Wigner’s semicircular distribution.

Understanding statistical properties of sums of large numbers of strongly correlated random variables is a classical problem in probability theory with important applications in natural sciences and beyond.Important examples of strongly correlated random variables are provided by the real eigenvalues λ i of the classical β-Gaussian, β-Wishart and β-Jacobi ensembles of random matrices [1].Random matrix theory has numerous applications ranging from nuclear physics to condensed matter physics, quantum information and quantum chromodynamics, deep learning and finance, to name just a few [2,3].Here a natural question is about spectral linear statistics, that is the probability distribution of the values A of sums of the type N i f (λ i ), where f (. . . ) is a given function, and N ≫ 1.Not surprisingly, spectral linear statistics of large random matrices have attracted much attention in the last 20 years, see Refs.[4][5][6][7] and references therein.Another important class of systems for which large deviations of linear statistics is of interest, and can be studied by similar methods, is provided by a family of non-interacting and interacting spinless fermions in confining potentials.Their close connection to random matrices is well documented [8,9].
In this work we consider β-Gaussian ensembles of large random matrices and study the probability distribution P(A, m) of the values A taken by the sum of absolute values of the eigenvalues raised to an arbitrary real power m > 0: The linear statistics (1) is a particular case of the general linear statistics A typical behavior of the fluctuating quantity A from Eq. ( 2) -as described by the mean Ā and the variance var(A) -has been studied in many papers [10][11][12][13][14][15][16][17][18][19][20][21][22] starting from the classical work of Dyson and Mehta [10].Here we focus on large deviations of A, defined by Eq. ( 1), which correspond to the large-A tails of P(A, m).They have received much less attention, except for the special case of m = 2, where P(A, 2) -a gamma-distribution -is known exactly, see e.g.Ref. [23] 1 .Our approach to this problem is based on a standard Coulomb gas technique, see e.g.Ref. [1].It exploits the fact that the joint probability distribution of the eigenvalues of the (unconstrained) Gaussian β-ensemble, P β (λ 1 , λ 2 , . . ., λ N ) ∝ can be interpreted as a measure of a 2D Coulomb gas of N particles confined on the line by a one-dimensional harmonic potential: where β plays the role of inverse temperature, and the energy of the Coulomb gas is Rescaling the eigenvalues, λ = √ N x, we obtain 1 A similar in spirit problem appears in the case of a Coulomb gas of N particles confined by the potential V (x) = x − ln x at x > 0.
Large deviations of the linear statistics N i x −1 i of such a gas were studied in Ref. [24].
up to a constant shift.We now define the linear statistics in the rescaled variables as Our goal is to derive the probability density of A which, by definition, can be written as where is the partition function of the Coulomb gas, unconstrained by Eq. (7).Using the exponential representation of the delta function, we can rewrite Eq. ( 8), up to an overall constant factor, as where is the constrained energy.Exploiting the large parameter N → ∞, we can reformulate the problem in a continuum language by introducing the density of rescaled eigenvalues which can be treated as a continuous field.The constrained energy can be expressed as a functional while the multiple integrals over x i turn into functional integrals over the (normalized to unity) density ρ(x): The O(N ) terms in this expression describe contributions from entropy and self-interactions [25,26].In the large-N limit the energy terms O(N 2 ) are dominant, so that we can neglect the O(N ) terms.Using again the exponential representation of delta functions, we obtain where the modified energy, now includes both Λ and µ.At large N the functional integrals in Eq. ( 15) can be evaluated by the saddle-point method.This brings about a variational problem, where one should minimize the modified energy E[ρ(x), Λ, µ] over all possible density configurations ρ(x), and ultimately express parameters Λ and µ (which play the role of Lagrange multipliers) through A by using the continuum version of Eq. ( 7), and the condition that total mass of the gas, M , is equal to unity: Once the solution of the minimization problem, ρ(x) = ρ A (x), Λ = Λ A and µ = µ A , is found, we can evaluate P(A, m), up to a pre-exponential factor, as with a rate function Here ρ W (x) is the density profile unconstrained by A: the famous Wigner's semicircular distribution [27] ρ while µ 0 is the corresponding unconstrained value of µ [1].It is evident from Eq. ( 20) that the rate function Ψ(A, m) vanishes when ρ A (x) = ρ W (x). In this case A is equal to its (unconstrained) expected value A = Ā ≡ Ām .We perform the minimization of the modified energy functional (16) and calculate the collective rate function Ψ(A, m) in Sec.II.Importantly, the Wigner's semicircular distribution (21) lives on a compact support, |x| < √ 2. As we shall see, the conditioned distributions ρ A (x) also have compact support.Somewhat surprisingly, however, we find that for 0 < m < 2 a single-connected compact support of ρ A (x) breaks up into two separate compact supports once A exceeds an m-dependent critical value.Such a breakup corresponds to the appearance of a macroscopic gap in the conditioned spectrum of the matrix.In the supercritical region the gap grows with A. A qualitatively similar, but quantitatively different, scenario occurs in the Ginibre ensemble (non-Hermitian random matrices with Gaussian entries) and the associated 2D one-component plasma, where the disc-annulus phase transition of the optimal charge density was observed [28,29] 2 .We discuss the emergence of the gap and a related phase transition in the limit of N → ∞ in Sections II B and II C.This phase transition is related to the nature of a singularity of the optimal charge density at x = 0.The order p of this phase transition is m-dependent, and we present analytical and numerical arguments which show that it is equal to p = 1 + 1/|m − 1|.In particular, the transition becomes of infinite order for m = 1.
For 0 < m < 1 the critical value of A for the breakup coincides with the expected value of the linear statistics A = Ā.In other words, the transition occurs in this case at the ground state of the Coulomb gas, which corresponds to the Wigner's semicircular distribution (21).
Crucially, the collective large-deviation scaling, described by Eq. ( 19), does not hold for all values of A and m.As we show in Sec.III, for m > 2 the collective scaling can give way to a single-particle scaling of the form Here, as well as in Eq. ( 19) for the collective rate function, we have neglected subleading terms -those with smaller powers of N -in the exponent.The single-particle scaling (22) holds, in the limit of N → ∞, at A > Ā and fixed A − Ā = O(1).Noticeable in Eq. ( 22) is the factor N 1+ 2 m instead of N 2 which appears in Eq. (19).At m > 2 the probability ( 22) is much higher than the collective one (19).Similarly to the breakup transition, this dramatic change of scaling behavior occurs at the ground state of the Coulomb gas.However, in contrast to the breakup, here it has the character of a secondorder phase transition for all m > 2. The mechanism behind the single-particle scaling (22) is the "evaporation" phenomenon, where either the maximum, or the minimum eigenvalue "escapes" from the rest of the eigenvalues, and it dominates the linear statistics, while the rest of the eigenvalues stay inside the Wigner's semicircle ρ W (x). Similar evaporation phenomena have been encountered in problems of statistics of the maximum or minimum eigenvalue of random matrices [29,32,33].To our knowledge, the present work is the first one to have observed the evaporation scenario for a linear statistics of Gaussian random matrices 3 .
By its very nature, the present work deals with rare events that are virtually impossible to sample by means of conventional Monte-Carlo simulations.Therefore, we verified our main theoretical predictions by using the multicanonical sampling of GOE matrices (β = 1) [34], based on the Wang-Landau algorithm [35].
We shall present the calculation of the single-particle rate function Ψ s (A, m) in Sec.III.Our main results are briefly summarized and discussed in Sec.V. Figure 1 shows a phase diagram of this system on the (A, m) plane.The curved line corresponds to the ground state of the system A = Ā, when the Coulomb gas density is described by the Wigner's semicircular distribution (21).

II. COLLECTIVE SCALING OF RATE FUNCTION
The minimization of the modified energy (16) brings about three conditions [7].One of them, yields a linear integral equation for ρ(x) which describes the mechanical equilibrium of charges in an effective potential This equation holds in the region(s) of space where ρ(x) > 0. The other two conditions are The first of them imposes the condition (17) and yields a relationship between Λ and A, whereas the second one imposes the normalization condition (18).It is impossible to satisfy Eq. (24) if ρ(x) > 0 on whole line |x| < ∞, because otherwise the l.h.s. and r.h.s.would have different asymptotics at |x| → ∞.(In the particular case of Λ = 0, which corresponds to the Wigner's semicircle, this fact is well known, see e.g.Ref. [1].)Therefore, the solution of Eq. ( 24) must have finite support.Once we have realized it, we should add one more condition to Eq. ( 23), which follows from variation of the a priori unknown positions of the boundaries of the compact support: ρ = 0 at support boundaries .(26) Equation (24), subject to the normalization condition (18), is equivalent to the equation subject to the same normalization condition.Equation ( 27) can be obtained by differentiating Eq. ( 24) with respect to x.It describes the force balance of the charge at point x: the l.h.s is the net force from the other charges, while the r.h.s. is the force exerted by the effective potential V eff (x, Λ).
It will be useful for the following to consider a more general family of solutions ρ(x) of Eq. ( 24) [or, equivalently, of Eq. ( 27)] which are not subject to the normalization condition (18) and can have an arbitrary mass M .These solutions can be parameterized by Λ and µ, and they have an important scale-invariance property.Suppose that ρ(x) is such a solution.Then the corresponding mass, moment, and energy of the Coulomb gas are where the last equality is a consequence of Eq. ( 24), and I m is defined in Eq. (17).Introducing the rescaling we can see that ρ solves Eq. ( 24) with rescaled values of Λ and µ: The rescaled mass, moment and energy are Now let us return to Eq. ( 24) subject to the normalization condition (18).Before dealing with arbitrary A, we note that the Wigner's distribution (21) follows from these equations in the absence of any constraint on A, that is for Λ = 0.In order to see it, one should assume that ρ(x) has a single-connected support x ∈ (− √ B, √ B) [where B is to be found alongside with ρ(x)] and solve the Carleman's equation (24).For √ B ̸ = 2 the solution is [36] ρ where the effective potential is The first integral in Eq. ( 34), and similar integrals in the following, are understood as Cauchy principal value integrals.The integrations can be performed explicitly, and we obtain a two-parameter solution: By virtue of the boundary condition (26), the numerator in Eq. ( 35) must vanish at x = ± √ B. This condition eliminates one of the parameters, say µ in favor of B, and we arrive at a one-parameter solution: The remaining parameter B is determined from the normalization condition (18).This gives B = 2 leading to the Wigner's semicircular distribution (21).The expected value Ā = Ā(m) of the linear statistics, corresponding to the unconstrained density ρ 0 (x), is where Γ(. . . ) is the gamma function.The dependence of Ā on m is shown in Fig. 2.

A. No spectral gap
Let us consider the full Eq.( 24), which is conditioned on A via the Lagrange multiplier Λ, and again assume that the solution ρ A (x) of Eq. ( 24) is defined on a single-connected compact support |x| < √ B. For √ B ̸ = 2 the solution [36] can be written as where ℜ[ 2 F1 1, m+1 2 ; m+4 2 ; B x 2 ] is the real part of the regularized hypergeometric function 2 F1 (a, b, c; z).As in the case of unconstrained density, we can get rid of the extra parameters by using the constraints ( 17), ( 18) and (26).Equation (18) yields After some algebra, we can eliminate µ.Then Eq. ( 17) yields Finally, the condition ρ(|x| As this equation for B is implicit, it is convenient to represent the solution in the parametric form with B as the single parameter: The rest of parameters A, Λ, and µ can be expressed in terms of B with the help of Eqs. ( 39) and (40): Several examples of the optimal density for different values of B and m are plotted in Fig. 3. To remind the reader, (i) B = 2 corresponds to the semicircular Wigner distribution ( 21), (ii) for any m we have A(B = 2) = Ā as described by Eq. ( 37), and (iii) Λ(B = 2) = 0.For m = 2 one obtains a rescaled Wigner distribution . This is an example of the scale-invariance property (31 The rescaled effective potential Eq. ( 24) is V eff = (1/2)x 2 leading to the Wigner semicircle.A < Ā corresponds to Λ < 0 and to "squeezed" distributions.Indeed, the effective potential V eff (x, Λ) confines the gas stronger as Λ decreases, and vice versa.42).The left panel corresponds to B = 1, for which A < Ā, where Ā is given by Eq. (37).The shown examples are: m = 4 so that A ≃ 0.16 and Ā = 0.5 (black); m = 2 so that A ≃ 0.25 and Ā = 0.5 (blue), and m = 5/4 so that A ≃ 0.33 and Ā ≃ 0.56 (red).The right panel corresponds to B = 2.5, for which A > Ā.The shown examples are m = 6 so that A ≃ 0.68 and Ā = 0.5 (black), m = 2 so that A ≃ 0.63 and Ā = 0.5 (blue), and m = 5/4 so that A ≃ 0.68 and Ā ≃ 0.56 (red).
Noticeable in Fig. 3 are a positive (at A < Ā) or negative (at A > Ā) peak of the optimal density (42) at x = 0.The asymptotic behavior of ρ A (x; B) in the vicinity of x = 0 is the following: That is, the single-interval solution (42) develops a power-law singularity (for 0 < m < 1) or a logarithmic singularity (for m = 1) at the origin.The character of the singularities (46) is such that the solution for 0 < m ≤ 1 becomes negative -and therefore inadmissible -for any Λ > 0. For 1 < m < 2 there is no singularity at x = 0, but the leading term in the expansion (46) becomes negative when B exceeds a critical value B cr > 0, so the single-interval solution breaks down here as well.For all m the critical value B cr > 0 is the minimum value of B at which the solution vanishes at x = 0. Using Eq. ( 43), we can also find the corresponding critical value of A, so that the single-interval solution breaks down at A > A cr .The critical values of B and A are for 1 < m < 2; (47) As we found, at B > B cr and 0 < m < 2 the single-interval solution (42) breaks down, and the true optimal density of the Coulomb gas splits into two domains, separated by a gap around x = 0. Remarkably, for 0 < m ≤ 1, this transition occurs already in the ground state of the system which corresponds to the Wigner semicircle and to A = Ā, see Fig. 4. 1.

4.
m Rescaled critical value of the linear statistics Acr/ Ā, see Eq. ( 47), for the transition from a single-to a double-interval optimal density, as a function of m.
For m > 2 the single-interval solution (42) continues to be formally valid for all B < B * = 2m m−2 4 .According to Eq. ( 43), the corresponding critical value of A is However, at N → ∞, the true minimum of the action here, already at Λ > 0, that is A > Ā, is reached via oneparticle evaporation.We will analyze these regimes, and the corresponding phase transitions, in Sections II B and III, respectively.In the remainder of this Section we will complete the solution of the problem for B ≤ B * by calculating the rate function Ψ(A) which determines the probability distribution (19).Before doing that, we plot in Fig. 5  Our results for A = A(B) and Λ = Λ(B) in Eqs. ( 43) and (44), respectively, make it possible to calculate the rate function Ψ 1 (B) in a parametric form (we use the subscript 1 to indicate single-interval support) by employing the "shortcut relation" dΨ 1 /dA = Λ, see e.g.Ref. [6].We can write dΨ 1 /dA = (dΨ 1 /dB)(dB/dA) and obtain the equation Using Eq. (43) for A(B) and integrating Eq. ( 49) from B = 2 to B , we obtain The same result (50) can be also obtained from Eq. ( 20) by evaluating the total energy of the Coulomb gas in Eq. ( 16), or by substituting the second moment, into the r.h.s. of Eq. (30).In both cases one should subtract the energy, corresponding to the unconstrained Wigner distribution, E[ρ W (x)] = 1 8 (3 + ln 4) ≃ 0.54828 . . . .Equation (50) alongside with A = A(B) from Eq. ( 43) give Ψ 1 = Ψ 1 (A) in the parametric form.The plot of Ψ 1 (A) for a particular value of m = 4 is shown in Fig. 6.As to be expected on the physical grounds, the minimum of the rate function is at A = Ā.A Taylor expansion of Ψ 1 (A) in the vicinity of A = Ā yields a quadratic asymptotic Alongside with Eq. ( 19), this asymptotic describes typical, Gaussian fluctuations of A around the mean with the variance N −2 σ 2 , where For even integer values of m this variance coincides with the one obtained in Ref. [4], where typical fluctuations of the quantity For m = 2, the collective rate function, determined by Eqs. ( 43) and (50), can be easily expressed in terms of A: When multiplied by β = 2, this rate function coincides with the one describing large-deviation statistics of the potential energy, (1/2) N i=1 x 2 i , of N ≫ 1 non-interacting fermions in a harmonic trap [23].The latter rate function arises from the large-N asymptotic of a gamma distribution [23].

B. Spectral gap
To determine the optimal density and the rate function in the regime with a gap, that is at A > A cr , we have to solve the double-interval integral equation In the rescaled variables Eq. (55) becomes where k < |x| < 1 and 0 , and we dropped the tildes.The solution, according to Ref. [37], can be written in the following form: Here the sum of integrals over the two symmetric intervals is denoted for brevity by ), E(k), K(k) are the complete elliptic integrals, and k ′ = √ 1 − k 2 .A simplification comes from the fact that f (x), and as a result ρ(x), are even functions.Using this property, we can rewrite Eq. (58) as where We were able to perform the integration in the r.h.s. of Eq. ( 59) analytically (with the help of "Mathematica") only for m = 1.This calculation, presented in the next subsection, leads to the prediction of infinite-order phase transition at A = A cr in this special case.For other 0 < m < 2 we present a non-rigorous perturbative argument, which is based on the asymptotic behavior (46) of the single-interval solution ρ A (x; B) in the vicinity of x = 0.It predicts a phase transition at A = A cr with an m-dependent and, in general, fractional order.Where possible, we will support these predictions by numerical simulations.

m = 1: general
For m = 1 the integrals in Eq. ( 59) can be evaluated analytically.Noticing that and taking into account that R(−t) = R(t), we can rewrite Eq. (59) as follows: The integrals in the r.h.s of Eq. ( 62) can be evaluated analytically, and they are presented, together with the resulting expression for µ, in Appendix A. Combining Eqs. ( 62), (A1) and (A2), we obtain the following four-parameter solution for ρ(x; c, λ, k, B) (for x ≥ 0): The parameters c, λ, k, and B are determined from the normalization condition (18), the vanishing of the density at the boundaries of support (26), and the condition (17) on A: The first two conditions in Eq. ( 64) allow us to eliminate the parameters c and λ, which yields a two-parameter distribution ρ(x; k, B), some examples of which are shown in Fig. 7.The last two conditions in Eq. ( 64) are hard to implement analytically.Therefore, we calculated the rate function by numerically evaluating the integrals: varying the parameter k, determining the parameter B from the normalization condition, and then calculating A and the rate function (20) for each set of parameters: Figure 8 shows the full rate function Ψ(A) for m = 1.It includes two branches, corresponding to the single-interval solution (50) for A ≤ Ā and the double-interval solution for A > Ā, that we have just presented.Also shown are results of the Wang-Landau simulations, and excellent agreement is observed.

Large A asymptotics
In the general case of 0 < m < 2 one can obtain the large-A asymptotic of the rate function without the need for explicit evaluation of the integrals in the r.h.s. of Eq. ( 59).This is because when k approaches 1 (and thus A ≫ Ā), the double-interval solution ρ(x; k, B) splits into two "droplets" located far apart, as it is clearly seen from Fig. 7 for m = 1.Remarkably, this feature makes it possible to evaluate the rate function's leading term at large A without the need to know the exact form of the solution ρ(x) of Eq. (55).
Consider a solution ρ(x) of Eq. (55) which satisfies the normalization condition (18) and the boundary conditions (26) Taking into account the strong inequality √ b ≫ √ B − √ b and the normalization condition, we obtain the following estimates: In the leading order, we obtain A ≃ b m/2 .Further, we notice that the energy of the two-droplet configuration is dominated by the potential energy x 2 ρ(x)dx ∼ b, while the droplets' interaction energy ln |x−y|ρ(x)ρ(y)dxdy ∼ ln b when the integration involves different droplets.The contribution to the energy of the Coulomb gas coming from the self-interaction term ln |x − y|ρ(x)ρ(y)dxdy, i.e., when the coordinates x and y belong to the same droplet, is of order ∼ O(1); it does not depend on the distance between the droplets.The resulting estimate of the rate function's leading term is where we have used the relation A ≃ b m/2 , following from Eq. ( 67).In Fig. 9 we compare the asymptotic (68) with results of the Wang-Landau simulations of the linear statistics of 20 × 20 GOE matrices for different values of m, and observe a very good agreement.We reiterate that the asymptotic (68) was obtained without using the solution for the density ρ(x).In the case of m = 1, however, one can go further and obtain a large-A asymptotic of the optimal density itself, which turns out to be quite instructive.To this end, it is convenient to consider the force balance equation ( 27): For definiteness, let us consider x ∈ √ b, √ B .Then the first and the second terms in the l.h.s. of Eq. ( 69) describe the inter-droplet and intra-droplet interactions, respectively.We can establish the lower and upper bounds on the forces ρ(y)dy x−y acting between droplets by replacing x − y under the integrals by 2 , respectively.We have or, in the limit That is, the contribution from the inter-droplet interaction is small compared to that of the intra-droplet interactions.Therefore, the interaction between the droplets is insignificant.This leads us to a simple integral equation We immediately notice that the solution of this equation is just a shifted and rescaled semicircular Wigner distribution.Normalizing the total mass to unity, we obtain In Fig. 7 we compared the full double-interval solution ρ(x; k, B) (63) with the shifted and rescaled Wigner semicircle (A + 1)ρ[x(A + 1), A] from Eq. ( 73).As one can see, the simple semicircle approximation works surprisingly well even for a relatively small value A ≃ 2.98.

General
The change of topology of the support of the optimal density ρ(x) implies a phase transition at A = A cr .In the absence of a closed form solution for the rate function Ψ 2 (A), it is still possible to predict the order of the transition by using some quantitative but non-rigorous argument that we will present shortly.When possible, we will test these predictions numerically.
Our argument makes use of three types of solutions of Eq. ( 24) for ρ(x) which depend on m and are parametrized by Λ and µ.For all the solutions we assume that ρ(x) vanish at the support boundaries.In solutions of one type we will abandon the non-negativity condition.In the solutions of another type we will allow the gas to have mass M different from 1.These can be obtained from the solutions with M = 1 via the rescaling transformation (33).
The argument starts with the single-interval solutions, ρ 1 (x), with M = 1 but without the non-negativeness condition.They are described by Eqs. ( 41) and ( 43)-( 45), regardless of the sign of A − A cr .According to Eq. ( 46), such a solution can be presented in the vicinity of x = 0 as Here C(m) is a function of m only, and the sign of C(m) coincides with the sign of 1 − m.The parameter ρ 0 depends on m, and Λ.Further arguments depend on whether 0 < m < 1 so that A cr = Ā, or 1 < m < 2 so that A cr > Ā, see Eq. (47).

0 < m < 1
Here the second term in Eq. ( 74) is the leading one at x → 0, and it determines the sign of ρ(x = 0).Also, here we have Λ cr = 0, B cr = 2, and A cr = Ā.Consider now a solution ρ 1 (x) with 0 < Λ ≪ 1.For this solution the density ρ 1 (x) is negative on the interval of the length around x = 0.This interval gives a negative contribution, δM 1 , to the mass of this solution, M 1 = 1: This implies that the contribution of the rest of the gas (where ρ(x) ≥ 0) to the total mass M 1 = 1 is Now consider the solution ρ a (x) with the same µ = µ a = µ 1 and Λ = Λ 1 as for the solution ρ 1 (x), but with the non-negativeness condition restored.As we are keeping the same Λ and µ, this solution will have a larger mass M a > M 1 = 1.We do not know exact two-interval solutions of Eq. ( 24) in a closed form.Here we make an assumption that the gap of this solution has a length of the same order of magnitude as in Eq. ( 75) and, most importantly, the excess of mass of this solution, M 2 − 1, can be estimated as in Eq. ( 77): (The coefficients O(1) in Eqs. ( 77) and ( 78) are quite possibly different, but this does not affect our predictions below.)Finally, we consider the two-interval solution ρ 2 (x) of Eq. ( 24), which obeys the non-negativeness condition and has the same Λ = Λ 2 = Λ a = Λ 1 .It has, however, a different µ = µ 2 ̸ = µ a = µ 1 which is tuned so that M 2 = 1.By virtue of Eq. ( 78) it is natural to assume that, in order to bring the total mass M 2 to its correct value 1, one should change µ by a quantity O(1) Λ 1/(1−m) : As a result, the relevant functionals A and Ψ of the solution ρ 2 (x), get corrected in the same way as M 2 : and Importantly, A 1 − Ā and Ψ 1 are analytic functions of Λ, which Taylor series start with the terms ∝ Λ and ∝ Λ 2 , respectively.Using these properties and eliminating Λ from Eqs. ( 80) and (81), we obtain Here we keep only the leading term in Ψ(A) − Ψ 1 (A) coming from the leading term of Ψ 1 as function of Λ.The latter one is: Since Ψ 2 (A) obeys one additional constraint -non-negativeness -the second term in Eq. ( 82) must be positive.It immediately follows from Eq. ( 82) that the order p of the phase transition at 0 < m < 1 is equal to In this case the sign of ρ(x = 0) for the single-interval solution is determined by the first term, ρ 0 , in Eq. (74).This term vanishes at Λ = Λ cr , corresponding to B = B cr , see Eq. (47).As before, we start from the single-interval solution ρ 1 (x) (again, without the non-negativeness condition) in a small vicinity of B = B cr .According to Eq. (74), this solution can be presented in the vicinity of x = 0 as where δΛ = Λ − Λ cr , and both (∂ Λ ρ 0 ) Λ=Λcr and C(m) are negative.The parameter µ of this solution is equal to µ 1 (Λ), determined by Eq. (45).For δΛ > 0 the solution ρ 1 (x) becomes negative on an interval of length ∆ which can be evaluated, with the help of Eq. ( 84), as This interval gives a negative contribution, δM 1 , to the total mass M 1 = 1.This contribution can be estimated as This implies that the rest of the gas gives the following contribution to the total mass M 1 = 1: Now we consider the second solution, ρ a (x), of Eq. ( 24) with the same Λ and µ = µ a = µ 1 , but obeying the nonnegativeness condition.This is a two-interval solution with a gap of the size O∆.As previously, our main assumption is that the total mass M a of this solution can be estimated as in (87), maybe up to a coefficient O(1).
Finally, we consider the third solution, ρ 2 (x): a two-interval solution which obeys the non-negativeness condition, which has the same Λ = Λ 2 = Λ a = Λ 1 , but a different µ = µ 2 ̸ = µ a = µ 1 chosen so as to obey the normalization condition M 2 = 1.Again, it is natural to assume that, in order to make the total mass equal to 1, one should change the parameter µ by a quantity O(1) Λ 1+ 1 m−1 : As a result, the functionals A and Ψ are corrected in the same way as M : and Now, A 1 (Λ) and Ψ 1 (Λ) are analytic functions However, in contrast to the case 0 < m < 1, the Taylor series of A 1 (Λ) and Ψ 1 (Λ) at the point Λ = Λ cr both start from the terms A cr + O(1) δΛ and Ψ cr + O(1) δΛ, respectively.Using these properties and eliminating Λ from Eqs. ( 89) and (90), we obtain: where we have kept only the leading term in Ψ(A) − Ψ 1 (A).Again, the second term in Eq. ( 91) should be positive.The resulting order p of the phase transition in this case is equal to Combining Eqs. ( 83) and (92), we obtain A plot of the order of this (in general, fractional) phase transition as a function of m is shown in Fig. 10.In particular, Eq. (93) indicates that, at m = 1, the order of the transition is infinite.In this special case, which is formally beyond the validity of our arguments in subsections II C 2 and II C 3, the single-interval solution has a logarithmic singularity at x = 0, see Eq. ( 46).Note that similar logarithmic singularities of the optimal density profile have already been associated with infinite-order phase transitions in large deviations of some different (truncated) linear statistics of the Gaussian [30] and Laguerre [31] ensembles of random matrices.Albeit quantitative, our arguments leading to Eq. ( 93) are far from rigorous.Therefore, we tested our predictions numerically.For the special case of m = 1 we used numerical evaluations of the analytical expressions from Sec. II B 1.For m ̸ = 1 we employed the Wang-Landau simulations.

III. SINGLE-PARTICLE SCALING OF RATE FUNCTION
As we have seen, for m < 2 large values of A are achieved when the optimal charge density of the ensemble splits into two separate "droplets", and the distance between them grows with A. Such a scenario, however, is impossible for m > 2 due to a qualitative change in the form of the effective potential V eff (x, Λ), see Eq. (24).Indeed, for Λ > 0 and m > 2, the double-well form of the effective potential disappears.Instead, the effective potential now has a local minimum at x = 0, and it is unbounded from below at |x| → ∞.As a result, for sufficiently large Λ > 0 Eq. ( 24) does not have physically reasonable solutions.
In this case a different type of optimal solution emerges via the evaporation scenario [38].In this scenario, one point charge "evaporates" from the rest of charges and is located at an a priori unknown position x = X to be found, while the rest of charges can still be described as a continuous "sea": The evaporation changes in the modified energy of the Coulomb gas ( 16): see Eq. ( 50).As one can see, at A− Ā = O(1) > 0 and m > 2 the probability (107) is higher than the probability (108).Overall, for m > 2, N → ∞, and arbitrary A − Ā = O(1) we obtain which summarizes our leading-order predictions for m > 2. In the limit of N → ∞ the rescaled rate function Ψ(A) = β −1 N −2 ln P(A, m) vanishes for A > Ā, signaling a second-order phase transition at the ground state A = Ā.How do these leading-order predictions compare with the results of our GOE matrix simulations?The latter are shown, for m = 4 and a moderately large value N = 60, in the left panel of Fig. 13.The dots represent the simulation results for the rescaled distribution N −2 ln [P(A, m)].The blue lines show the collective rate function Ψ 1 (A) from Eq. (50), and the red line shows the single-particle rate function N −1 Ψ s {X = N (A − Ā) 1/m } from Eq. (105).As one can see, the single-particle rate function describes the right tail of the numerical data remarkably well.However, it deviates from the simulated data for A closer to Ā.In its turn, the collective rate function Ψ 1 continues to be accurate, in a limited range of A, for A > Ā where, according to Eq. ( 109), it should not apply.These differences between the asymptotic N → ∞ theory and simulations result from insufficiently large N , and they decrease with an increase of N .The latter fact is evident on the right panel of Fig. 13, which shows the results of matrix simulation for different N .As one can see, the range of values of A, not covered by either the collective rate function Ψ 1 , or the single-particle rate function Ψ s , shrinks with an increase of N5 .Unfortunately, it seems impractical to reach the asymptotic regime of Eq. (109) in the matrix simulations.IV. m ≫ 1: FROM SECOND-TO THIRD-ORDER TRANSITION We have just seen that, at m > 2 the rate function of spectral linear statistics (1) in the limit of N → ∞ exhibits a phase transition, described by evaporation scenario.This transition occurs at the ground state and, for any finite m > 2, it is always of second order.A natural question arises on what happens in the limit of m → ∞ when the statistics (1) is dominated by the maximum value of |λ i | m , i = 1, 2, . . ., N .The statistics of the largest eigenvalue λ max itself (without the absolute value) has been studied in several works, see Ref. [32] and references therein.For this and other related statistics the evaporation scenario is also at work.In the limit of N → ∞, the evaporation scenario predicts a vanishing largest-eigenvalue rate function β −1 N −2 ln P(λ max ) for λ max > √ 2. However, the associated phase transition, governed by the collective rate function for λ max < √ 2, is of third order, see e.g.Ref. [32].To clarify this issue, here we consider in some detail the limit of m ≫ 1 of the single-interval rate function, described by Eqs. ( 43) and (50).The limit of m → ∞ can be taken in different ways.Let us introduce ∆A ≡ A − Ā.First we consider the limit of m → ∞ at fixed ∆A/ Ā = O(1).To this end we use Eqs.(37) and (43) and present A/ Ā in the following form: It is clear from Eq. ( 110) that, in order to have ∆A/ Ā = O(1) at m → ∞, we must demand that 2 − B → 0 so that u ≡ m(2 − B) = O(1).In the limit of m → ∞ at fixed u, we obtain where we have used the identity Under the same scaling u ≡ m(2 − B) = O(1), Eq. ( 50) becomes Equations ( 111) and (112) determine in implicit form the m ≫ 1 asymptotic of the rate function Ψ 1 as a function of ∆A/ Ā: The function f (. . . ) is defined by Eqs. ( 111) and (112).At small arguments it can be expanded as where we have used the small-u expansion of Eq. ( 111 We recall that B is the square of the radius of support of the optimal density distribution for a given A < Ā.On the other hand, using Eq.(43), we obtain That is, in the limit of m → ∞, the statistics of A 2/m is equivalent to the statistics of B itself.It is not surprising, therefore, that the rate function (116) can be also obtained, in a straightforward way, from the known rate function [26] of the joint probability distribution of the minimum and the maximum eigenvalues λ min and λ max of the same class of matrices.One should only set there, for the rescaled eigenvalues, x min = −x max = − √ B, where 0 < x max < √ 2. For small 2 − B > 0 Eq.(116) yields As one can see, the transition does become of third order in the limit of m → ∞, but for the order parameter A 2/m rather than A6 .

V. SUMMARY AND DISCUSSION
As we have shown here, the large-deviation behavior of the probability distribution P(A, m) at N → ∞ is quite intricate, and the phase diagram of the system on the (A, m) plane, depicted in Fig. 1, is unexpectedly rich.This is especially true in the region 0 < m < 2 of the phase diagram, where we observed a change in topology of the optimal density of the Coulomb gas at A = A cr accompanied by a phase transition of a collective nature.This phase transition occurs at a critical value of A that we determined.The transition has an m-dependent order p, such that 2 < p ≤ ∞, see Eq. (93) and Fig. 10.
In the future work one can try to improve our non-rigorous analytical arguments leading to the prediction of p = p(m) in the regime of 0 < m < 2.An obvious obstacle on the way is the difficulty in the evaluation, for arbitrary 0 < m < 2, of the integrals that appear in the expression (58) for the exact solution of the two-interval Carleman's equation (57).One can try to circumvent this difficulty by developing a perturbative method of evaluating these integrals which employs from the start the small parameter k ≪ 1.
According to our non-rigorous arguments, supported by numerics, the transition order at 0 < m < 2 is determined by the type of singularity at x = 0, exhibited by the single-interval optimal density profile.It is reasonable to assume that a phase transition of the same order can occur for other choices of linear statistics, confining potential and inter-particle interaction if those lead to the same type of density singularity.It would be very interesting therefore to look for similar phase-transition phenomena in large deviations of linear spectral statistics in additional systems which are describable in terms of an ensemble of particles confined by an external potential.Immediate examples are provided by non-interacting and interacting spinless fermions in confining potentials [8,9], as well as interacting classical gases, not necessarily related to random matrices, such as those studied in Refs.[40,41].These examples may include systems where a rigorous determination of the transition order is more straightforward.
For m > 2, we found a different kind of phase transition, which is described by evaporation scenario.The transition occurs at the ground state of the Coulomb gas and, for any finite m > 2, it is always of second order.If one keeps A as the order parameter, the transition remains of second order even in the limit of m → ∞.However, when switching, e.g., to the order parameter A 2/m , the transition becomes of third order in this limit.This happens because, as m → ∞, the order parameter A 2/m coincides with B, see Eq. (117).A third-order phase transition of the same nature occurs in the statistics of the maximum (or minimum) eigenvalue, see Ref. [32] and references therein.
To conclude, we believe that the unexpected richness of the phase diagram and of the phenomenology of the spectral linear statistics (1), uncovered in this work, justifies further study of these statistics.In particular, it poses interesting questions about the possibility of similar behaviors in other systems of interacting particles.

I.
Introduction and formalism II.Collective scaling of rate function A. No spectral gap B. Spectral gap 1. m = 1: general 2. Large A asymptotics C. Order of the splitting phase transition 1.General 2. 0 < m < 1 3. 1 < m < 2 4. Numerical results: m = 1 5. Numerical results: 0 < m < 1 III.Single-particle scaling of rate function IV. m ≫ 1: from second-to third-order transition V. Summary and discussion I. INTRODUCTION AND FORMALISM

FIG. 1 .
FIG.1.Phase diagram of the large-deviation statistics (1) on the (A, m) plane.The curved line corresponds to the ground state of the system A = Ā, when the Coulomb gas density is described by the Wigner's semicircular distribution(21).

FIG. 2 .
FIG. 2. Expected value Ā of the linear statistics as a function of m.
where B > b > 0, ρ(x) is normalized to unity, and the parameters B and b are a priori unknown.Let us rescale x/ √ B = x, denote b/B = k > 0, and introduce a rescaled density ρ

FIG. 8 .
FIG. 8.The full rate function Ψ(A) for m = 1.The left branch (blue) and the right branch (red) correspond to the single-and double-interval solutions, respectively.The black circles show the results of Wang-Landau simulations of the linear statistics of 20 × 20 GOE matrices.The blue dashed curve shows analytic continuation of a single interval-rate function to the region of A > Ā.It deviates from the red curve very slowly, suggesting that the phase transition at A = Ā is of infinite order.The green curve in the right inset is the large A asymptotic (68).
. Let the support of ρ(x) have the form of two droplets at (− √ We can establish simple lower and upper bounds on the linear statistics A of this configuration by replacing |x| m under the integrals by b m/2 and B m/2 = b m/2 1 +

FIG. 10 .
FIG.10.The order of phase transition between the gapless and gapped spectra as a function of m.

FIG. 13 .
FIG. 13.Left panel: the simulated rate function Ψ(A) = N −2 ln P(A, m) (dots) for m = 4 and N = 60, the collective rate function Ψ1(A) (blue), and the single-particle rate function N −1 Ψs(A, N ), corresponding to the evaporation scenario (red).The inset shows a blowup of the vicinity of A = Ā for m = 4. Right panel: the right tail of the rate function N −2 ln (P(A, m)).The dots represent the results of GOE matrix simulations for N = 30 (purple), N = 60 (black), and N = 120 (orange).The solid lines represent the single-particle rate functions (105) for the corresponding values of N .Note that, for m = 4, the critical value A * from Eq. (48) is equal to 1.