Number statistics for β-ensembles of random matrices: applications to trapped fermions at zero temperature Ricardo Marino Department of Physics of Complex Systems, Weizmann Institute of Science, 76100 Rehovot, Israel. Satya N. Majumdar, Gr´egory Schehr LPTMS, CNRS, Univ. Paris-Sud, Universit´e Paris-Saclay, 91405 Orsay, France. Pierpaolo Vivo 6 King’s College London, Department of Mathematics, Strand, London WC2R 2LS, United Kingdom. 1 (Dated: January 14, 2016) 0 2 Let Pβ(V)(NI) be the probability that a N ×N β-ensemble of random matrices with confining potential V(x) has N eigenvalues inside an interval I = [a,b] of the real line. We introduce a n I general formalism, based on the Coulomb gas technique and the resolvent method, to compute a J analyticallyPβ(V)(NI)forlargeN. WeshowthatthisprobabilityscalesforlargeN asPβ(V)(NI)≈ (cid:16) (cid:17) 3 exp −βN2ψ(V)(N /N) ,whereβ istheDysonindexoftheensemble. Theratefunctionψ(V)(k ), I I 1 independentofβ,iscomputedintermsofsingleintegralsthatcanbeeasilyevaluatednumerically. The general formalism is then applied to the classical β-Gaussian (I = [−L,L]), β-Wishart (I = ] h [1,L])andβ-Cauchy(I =[−L,L])ensembles. Expandingtheratefunctionarounditsminimum,we c findthatgenericallythenumbervarianceVar(NI)exhibitsanon-monotonicbehaviorasafunction e of the size of the interval, with a maximum that can be precisely characterized. These analytical m results,corroboratedbynumericalsimulations,providethefullcountingstatisticsofmanysystems whererandommatrixmodelsapply. Inparticular,wepresentresultsforthefullcountingstatistics - t of zero temperature one-dimensional spinless fermions in a harmonic trap. a t s . t a m - d n o c [ 1 v 8 7 1 3 0 . 1 0 6 1 : v i X r a 2 CONTENTS I. Introduction 2 II. Setting and summary of results 5 III. β-Gaussian ensemble 9 A. The Coulomb gas 9 B. Resolvent method 11 C. Calculation of the rate function ψ(G)(k ) 13 I D. Number variance 15 1. Extended bulk regime 16 2. Edge regime 18 3. Tail regime 19 E. Comparison with numerics 21 IV. β-Wishart ensemble 21 A. Calculation of the rate function ψ(W)(k ) 22 I B. Number variance 23 1. Extended bulk regime 23 2. Edge regime 23 3. Tail regime 25 C. Comparison with numerics 26 V. β-Cauchy ensemble 26 A. Analysis for [ L,L] 28 − B. Number variance 29 C. Comparison with numerics 30 VI. Conclusions 30 A. Asymptotic analysis of Gaussian extended bulk regime 31 References 33 I. INTRODUCTION Therehasbeenanintenseactivityinthefieldofcoldatomsinthelasttwodecades[1,2]. Sincethefirstrealizations ofBose-Einsteincondensation[3–5],newexperimentaldevelopmentsinvolvingthecoolingandconfinementofparticles inopticaltrapsledtoanewchapterinmany-bodyquantumsystems,wheretheparticlestatisticsisthemaininterest rather than individual atoms. After the achievement of cooling trapped fermions to the point where Fermi statistics becomes dominant [6], experiments were able to explore the remarkable properties of Fermi gases [7]. Non-interacting fermions exhibit non-trivial quantum effects, arising from Pauli exclusion principle. While confined bosons may collapse to the lowest level of the trap, fermions are forced to dilute and behave as a strongly correlated system, regardless of their original interaction. This induces a Fermi motion of the particles, which is a purely quantum phenomenon: it downplays to some extent the importance of individual interactions among the fermions, which can be, in many cases, neglected or treated as a small perturbation [2]. This turns the ideal Fermi gas into the natural playground to explore the properties of more general Fermi gases. Consideraone-dimensionalFermigasofN particlesconfinedbyaharmonicpotentialV (x)= 1mω2x2, wherethe Q 2 subscript Q stands for the quantum potential. For simplicity, we set m=ω =(cid:126)=1. The many-body Hamiltonian is then given by H =(cid:80)N (cid:2) 1∂2 +V (x )(cid:3). We want to calculate the many-body ground state wave function Ψ ((cid:126)x), i=1 −2 xi Q i 0 where(cid:126)x x ,x , ,x arethepositionsoftheparticlesonaline. Thiswavefunctioncanbeconstructedfromthe 1 2 N ≡{ ··· } single-particleeigenfunctionsofaharmonicoscillatorφ (x) e−x2/2H (x),(withenergyeigenvalues(cid:15) =(n+1/2)). n n n Here H (x) are Hermite polynomials of degree n. For exam∝ple, H (x) = 1, H (x) = 2x, H (x) = 2x2 4 etc. One n 0 1 2 − then writes Ψ ((cid:126)x) as a Slater determinant Ψ ((cid:126)x) = det[φ (x )]/√N!, with 0 i N 1 and 1 j N. By 0 0 i j ≤ ≤ − ≤ ≤ construction, this wave function vanishes whenever x = x for i = j, thus satisfying the Pauli exclusion principle. i j (cid:54) 3 This corresponds to filling each single-particle energy level from n = 0 to n = N 1 with a single fermion. The − highestoccupiedenergylevel,N 1,isthustheFermienergy. Thegroundstateenergyassociatedtothismany-body wave function is E = (cid:80)N−1(n−+1/2) = N2/2. Amazingly, this Slater determinant Ψ ((cid:126)x) = det[φ (x )]/√N! = 0 n=0 0 i j √1N!e−(cid:80)Ni=1x2i/2det[Hi(xj)] can be evaluated explicitly, by noting that det[Hi(xj)] ∝ (cid:81)i<j(xi−xj), which is just a Vandermonde determinant. Thus the squared many-body ground state wave function can be written as [8] Ψ0((cid:126)x)2 = 1 e−(cid:80)Ni=1x2i (cid:89)(xk xj)2, (1) | | Z − N j<k where Z is the normalization constant. This quantity Ψ ((cid:126)x)2 quantifies the quantum fluctuations in this system N 0 | | at zero temperature. It can indeed be interpreted as the probability density function of this fermionic system at zero temperature. NotethattheVandermondeterm(cid:81) (x x )2 inEq. (1)makesthevariablesx ’sstronglycorrelated. i<j i− j i Here the physical origin of this strong correlation is due to the Pauli exclusion principle satisfied by the fermions. Precisely the same joint distribution also appears in the Gaussian Unitary Ensemble (GUE) of Random Matrix Theory (RMT). Consider an N N random Hermitian matrix X with complex entries drawn from the distribution × Pr(X) e−Tr(X2) which is invariant under a unitary transformation. For each realization of X, one has N real ∝ eigenvalues x ,x , ,x . The joint distribution of the eigenvalues x ,x , ,x can be computed explicitly by 1 2 N 1 2 N ··· ··· making a change of variables from the entries to the eigenvalues and eigenvectors [9]. The eigenvectors decouple and thejointdistributionofeigenvalueshasexactlythesameformasinEq. (1),wheretheVandermondetermcomesfrom the Jacobian of this change of variables. Thus N fermions in a harmonic trap at zero temperature provide a physical realization of the GUE eigenvalues [8], the origin of the Vandermonde term in the two problems being however very different. This remarkable connection has allowed recently to explore various observables in the ground state of the trapped fermions system such as counting statistics and entanglement entropy [8, 10–16]. InRMTvariousensembleshavebeenstudiedextensively[9,17]. Inparticular,forrotationallyinvariantensembles, the joint probability distribution function (jpdf) of eigenvalues can be written explicitly as P((cid:126)x)= 1 e−βN(cid:80)Ni=1V(xi)(cid:89) xk xj β, (2) Z | − | N,β j<k where β > 0 is the Dyson index of the ensemble and V(x) is a potential growing suitably fast at infinity to ensure that the jpdf is normalizable. Usually, for rotationally invariant ensembles, the Dyson index is quantized: β =1,2,4 respectively for real symmetric, complex Hermitian and quaternionic self dual matrices. However there are other matrix ensembles where β >0 can be arbitrary [18]. In the case V(x)=x2/2 one recovers the β-Gaussian ensembles. Theotherwellstudiedexamplesareβ-Wishart(V(x)=x/2 αlnx)andβ-Cauchy(V(x)=((N 1)/2+1/β)ln(1+ x2)/N), described in detail in section II. The fermions in a−harmonic well at T = 0 thus corres−pond to the special case of the RMT with β = 2 and V(x) = x2/2 (GUE). Note that in Eq. (2) we have rescaled the potential V(x) by a factor N. This convention ensures that the eigenvalues are of order (1) (see the discussion later). With this O convention, which we will use in the rest of the paper, the ground state squared wave function for fermions in Eq. (1) reads |Ψ0((cid:126)x)|2 = Z˜1 e−N (cid:80)Ni=1x2i (cid:89)(xk−xj)2, (3) N j<k where Z˜ is the new normalization constant (in this particular case, this expression (3) is obtained from Eq. (1) N just by a rescaling of the fermions’ positions x √Nx ). In the context of step fluctuations in vicinal surfaces of a i i → crystal, such a fermionic representation of GUE was already noticed (see [19] for a review). For the steps in presence of a hard wall, a similar fermionic representation was found that corresponds to the Wishart case of RMT with β =2 [20]. One of the most natural observables in both RMT as well as in many body quantum physics is the global density N 1 (cid:88) ρ (x)= δ(x x ), (4) N i N − i=1 where x ’s are the positions of eigenvalues or fermions on the line. Note that ρ (x) is normalized to unity: i N (cid:82) { } ρ (x)dx = 1 and thus measures the fraction of eigenvalues in the interval [x,x+dx] in any given sample. The N average density of eigenvalues ρ (x) , where stands for an average over the jpdf in Eq. (2), converges to the N (cid:104) (cid:105) (cid:104)···(cid:105) celebrated Wigner semicircle law [9, 21] in the large N limit ρ (x) N(cid:29)1 ρ (x), ρ (x)= 1(cid:112)2 x2 . (5) N sc sc (cid:104) (cid:105)−−−→ π − 4 N 2/3 N 1 N 2/3 − − − edge edge ) x ( sc bulk ρ √2 √2 − x Figure1. AveragedensityofscaledGaussianeigenvalueswitharepresentationofbulkandedgeregimes. Thebulkisthescale of typical fluctuations of eigenvalues, whose width is of order 1/N. The edge is the scale of typical fluctuations of the largest eigenvalue. Therefore the eigenvalues are spread, on an average, over a finite interval [ √2,+√2]. The typical spacing between − eigenvalues, near the origin, is thus of order (1/N). This regime near the center is known as the bulk, as explained O in Fig. 1. In contrast, as one approaches the edges √2 the density becomes smaller, indicating that the eigenvalues ± aresparserneartheedges. Indeed, itiswellknownthatthetypicalspacingbetweeneigenvaluesneartheedgesscales as (N−2/3) for large N. This regime is known as the edge regime (see Fig. 1). O An important observable to characterize quantum fluctuations of cold fermions is their number statistics, or full counting statistics. This is the number N of fermions inside an interval at zero temperature. For the ideal bosonic I I case at T =0, all particles occupythe ground state, centered at the minimum of thetrap. They are uncorrelated and their number statistics can be easily determined. In the fermionic case, where Pauli principle applies, the statistics of the random variable N is instead particularly interesting, as it reflects the non trivial quantum correlations between I thefermionpositions. Asasimpleexample, considerasymmetricinterval =[ L,L]forthefermionsinaharmonic I − trap. In the context of RMT, this is just the number of eigenvalues in [ L,L] in the GUE. This number N is a [−L,L] − random variable whose mean can be easily computed, for large N, from the average density given in Eq. (5) (cid:90) +L N (cid:18) (cid:112) (cid:18) L (cid:19)(cid:19) N =N ρ (x) dx L 2 L2+2sin−1 , 0 L √2. (6) [−L,L] N (cid:104) (cid:105) −L (cid:104) (cid:105) ≈ π − √2 ≤ ≤ For L √2, N saturates to its maximum value, i.e., N N. [−L,L] [−L,L] ≥ (cid:104) (cid:105) (cid:104) (cid:105)≈ However, the random variable N fluctuates around this mean value from sample to sample. It is therefore [−L,L] naturaltostudythehighercumulantsofthisrandomvariable, forinstancethe varianceVar(N ), andeventually [−L,L] the full distribution of N . In RMT, the statistics of N in Gaussian ensembles (β =1,2,4) is actually well [−L,L] [−L,L] known, but only in the bulk limit when L (1/N). In this case, the variance was computed by Dyson in [22] ∼O 2 Var(N ) ln(NL)+C +o(1), L (1/N), (7) [−L,L] ≈ βπ2 β ∼O where C is a β-dependent constant that was computed explicitly by Dyson and Mehta [23]. For instance, for β =2 β 1 C = (1+γ+ln2)=0.230036... , (8) 2 π2 where γ =0.577215... is the Euler constant. In this bulk regime, even the full distribution of N was computed [−L,L] inRefs. [24]and[25]anditwasshowntobeasimpleGaussian. However,beyondthebulkregime,i.e.,L (1/N), (cid:29)O the behavior of the variance Var(N ) (and also the full distribution of N ) was not studied in the RMT [−L,L] [−L,L] literature. In the context of fermions trapped in a confining potential, this number variance has found a recent renewed interest, as it was found to be closely related to the entanglement entropy at T =0 of the subsystem [ L,L] − with the rest of the system [10–16]. In that context, Var(N ) was studied numerically as a function of L for [−L,L] various confining potentials and a rather striking non-monotonous dependence on L was found [12, 13] (see Fig. 2). In a recent Letter [8], we studied Var(N ) analytically for fermions in a harmonic potential. Exploiting the [−L,L] mapping to the GUE eigenvalues and using a Coulomb gas technique, we were able to compute Var(N ) for any [−L,L] 5 L and large N. As a function of increasing L, we found the following behavior for the variance [8]  1 ln(cid:0)NL(2 L2)3/2(cid:1), N−1 <L<√2 N−2/3 π2 − − Var(N[−L,L])≈eVx2p(s[),2Nφ(L)], LL>=√√22++N√s2−N2/−3,2/3 (9) − where the functions V (s) and φ(L) were computed explicitly and are given in Eqs. (81) and (86) respectively. In 2 the bulk limit, using L (1/N) in the first line of Eq. (9), we recover the results of Dyson in Eq. (7) for β = 2. ∼ O In Fig. 2 these theoretical results (9) are compared to numerical simulations and one finds a very good agreement. In addition, the full distribution of N was also computed in Ref. [8] using the Coulomb gas method. The large [−L,L] deviation function associated with the distribution of N was found to have a non-trivial logarithmic singularity [−L,L] where the distribution has its peak. This singularity was attributed to a phase transition in the associated Coulomb gas problem [8]. 1.2 Theory NumericsforN=5000 1.0 0.8 ) I N ( 0.6 0.7 r a 0.6 V 0.5 0.4 s) 0.4 ( 20.3 V 0.2 0.2 0.1 0.0 −6 −5 −4 −3 −s2 −1 0 1 2 0.0 10−3 10−2 10−1 √2 L Figure 2. Results for the number variance of GUE eigenvalues when I =[−L,L]. The theoretical result is in equation (9). The goal of this paper is twofold. First, we present the details of the Coulomb gas calculations for β-Gaussian ensembles which were just announced in our Letter [8]. This is a generalization of the results for the β = 2 case mentioned in Eq. (9). Secondly, we show that our results can be generalized to compute the number variance of an interval for two other classical random matrix ensembles: β-Wishart and β-Cauchy ensembles. The rest of the I paper is organized as follows. We begin by recalling a few definitions of random matrix theory and summarizing our main results in section II. We then describe the details of the Coulomb gas method and the derivation of the results for the three ensembles: β-Gaussian in section III, β-Wishart in section IV and β-Cauchy in section V. We present our conclusions in section VI and Appendix A contains details of several calculations performed in section III. II. SETTING AND SUMMARY OF RESULTS We consider ensembles of random matrices for which the joint distribution of eigenvalues is given in Eq. (2) where V(x) denotes the matrix potential. In this paper we focus on three classical random matrix ensembles, β-Gaussian (V(x)=x2/2), β-Wishart (V(x)=x/2 αlnx) and β-Cauchy (V(x)=((N 1)/2+1/β)ln(1+x2)/N), described − − in details below. Our quantity of interest is the number N of eigenvalues inside an interval . This is a random variable which can I (cid:80) I be expressed as N = 1 (x ), where 1 (x ) is the indicator function = 1 when x and zero otherwise. Note I i I i I i i ∈ I thatfortheGaussiancaseifwechoosetheinterval =[0,+ ),thenN isjusttheindex,i.e.,thenumberofpositive I I ∞ eigenvalues. In this special case, the distribution of N was computed using Coulomb gas method in Refs. [26, 27]. I The same method was later used to compute the distribution of the number of Wishart eigenvalues in the interval 6 =[l,+ ) [28]. In all these cases, the interval was unbounded. Here, we show that this Coulomb gas method can I ∞ I also be used to compute the distribution of N where the interval is bounded, i.e., =[a,b]. For such an interval, I I I the probability density of N is given by I (cid:90) N (cid:32) N (cid:33) (V)(N )= (cid:89)dx P((cid:126)x)δ N (cid:88)1 (x ) , (10) Pβ I i I − I l i=1 l=1 where P(x ,...,x ) is the jpdf of the eigenvalues (2). The superscript (V) in Eq. (10) just refers to the potential 1 N V(x). The average N is easily obtained in the large N limit by integrating the average density over : I (cid:104) (cid:105) I (cid:90) N N ρ(x)dx, (11) I (cid:104) (cid:105)≈ I whereρ(x)=lim ρ (x) isthelimitingaveragedensityofeigenvalues,whichdependsexplicitlyonthepotential N→∞ N (cid:104) (cid:105) V(x). Whilethecomputationofρ(x)isrelativelystraightforwardandiswellknownforthethreeensemblesmentioned above, it is much harder to compute the variance, Var(N )= N2 N 2, and the full distribution (V)(N ). We I (cid:104) I(cid:105)−(cid:104) I(cid:105) Pβ I show that for large N, and as long as the interval is contained within the support of ρ(x), the distribution of N I I admits a large deviation form (V)(N =k N) e−βN2ψ(V)(kI), (12) Pβ I I ≈ where k =N /N denotes the fraction of eigenvalues in the interval and ψ(V)(k ) is the rate function independent I I I I of both β and N that can be derived by a general method described below. This rate function ψ(V)(k ) reaches its minimum at k = k = N /N (cid:82) ρ(x)dx [see Eq. (11)]. In standard I I I (cid:104) I(cid:105) ≈ I large deviation setting, usually the rate function exhibits a smooth quadratic behavior around this minimum. In our case, we find that for all three potentials the rate function, near its minimum, behaves as ψ(V)(k ) (k I I k )2/[ ln(k k /a )] where a is a non-universal constant that depends on the interval and on the p∝otentia−l. I I I I I − | − | I Theratefunctionthusindeedhasaquadraticbehaviorclosetok =k ,exceptthatitisactuallymodulatedbyaweak I I logarithmic singularity. This singularity occurs when we investigate the fluctuations on a scale N N (N). I I −(cid:104) (cid:105) ∼ O Writingln(k k /a )=ln N N ln(Na ),itisclearthatif N N N onecanreplaceln(k k ) I I I I I I I I I I | − | | −(cid:104) (cid:105)|− | −(cid:104) (cid:105)|(cid:28) | − | by ln(Na ). Thus on this much smaller scale, it follows from Eq. (12) that the fluctuations just become Gaussian I − (cid:34) (cid:0) (cid:1)2(cid:35) N Nk (V) I I (N =k N) exp − for N Nk , (13) Pβ I I ≈ − 2(Var(N )) I → I I with the variance growing, for large N, as Var(N ) ln(Na ). (14) I I ∝ Indeed, thisquadraticbehavior(modulatedbylogarithmicsingularity)oftheratefunctionnearitsminimumalready appearedinthespecialcasesofsemi-infiniteintervals =[0,+ )(fortheGaussiancase)[26,27]andfor =[l,+ ) I ∞ I ∞ (fortheWishartcase)[28]mentionedearlier. Herewefindthatthislogarithmicsingularityismoregeneralandholds also for arbitrary bounded intervals . I In summary, there are two scales associated with the fluctuations of N . There is a shorter scale, N N I I I (cid:112) −(cid:104) (cid:105) ∼ ( ln(Na )), which describes the typical fluctuations of N around its average and the distribution of these typical I I O fluctuations is Gaussian with a variance ln(Na ) for large N. However, the atypically large fluctuations where I ∝ N N (N) are not described by the Gaussian distribution, rather by the large deviation form in Eq. (12). I I The−lo(cid:104)gari(cid:105)th∼mOicsingularityofψ(V)(k )aroundk ensuresthematchingofthesetwoscales. Inthispaperwecompute I I theratefunctionψ(V)(k )forthethreeensemblesmentionedabove. Eventhoughtheratefunctionψ(V)(k )depends I I on the choice of the ensemble, we will show that its logarithmic singularity near its minimum is universal. As a consequence, the variance of the typical Gaussian fluctuations around the mean grows universally as lnN for large N. Below we provide a summary of the results for each of these three ensembles separately. β-Gaussian ensemble: V(x) = x2/2. This ensemble includes real symmetric (β = 1), complex Hermitian (β = 2) or quaternion self-dual (β = 4) matrices whose entries are Gaussian variables such as the probability density of the matrix X is given by Pr(X) e−βNTrX2/2. The average spectral density for large N and for any β >0 converges to the celebrated Wigner’s semi∝-circle law ρ (x), described in equation (5). The rate function ψ(G)(k ) in Eq. (12) is sc I given explicitly in Eq. (60) below. 7 N 1 − ) ) x x ( ( sc sc ρ ρ √2 LL √2 √2 L L √2 − −x − − x (i)Bulkregime. (ii)Extendedbulkregime. N 2/3 N 2/3 − − edge edge ) ) x x ( ( sc sc ρ ρ L L − √2 √2 L √2 √2 L − x − − x (iii)Edgeregime. (iv)Tailregime. Figure 3. Regimes of behavior of the number variance for the Gaussian ensemble. The solid blue line represents the interval I =[−L,L]. Our results for the number variance of the Gaussian ensemble are given by  2 ln(cid:0)NL(2 L2)3/2(cid:1), N−1 <L<√2 N−2/3 βπ2 − − Var(N[−L,L])≈eVxβp(s[),βNφ(L)], LL>=√√22++N√s2−N2/−3,2/3 (15) − and were announced in the Letter [8]. We were able to calculate V (s) for β =2, equation (81), and φ(L) is given by β equation(86). Aplotofequation(15)forβ =2,whichisequivalenttothestatisticsoftheharmonicallyconfinedideal Fermi gas when distances are scaled by √N and (cid:126) = m = ω = 1, is presented in figure 2 compared with numerical results. As L increases, the number variance for the interval [ L,L] changes behavior, and we identify four regimes, − highlightedinfigure3: (i)abulk regime, (ii)anextended bulk regime, (iii)anedge regimeand(iv)atail regime. The variance has a non-monotonic behavior, reaches its maximum at a non-trivial value (L(cid:63) = 1/√2) and drops sharply near the edge of the semicircle (L = √2). From the first line of equation (15), setting L = s/N, where s is of order (1), andtakingthelargeN limitwegetVar(N ) (2/βπ2)lns, thusrecoveringtheresultofDysoninEq. (7). [−L,L] O ∼ This is also evident in figure 2 where we see a linear growth on a logarithmic scale in the early regime. β-Wishart ensemble: In the Wishart ensemble [29] we consider a product matrix W =X†X where X is a M N × rectangular random matrix with Gaussian entries. W is thus an N N square random matrix whose entries are distributed according to Pr(W) e−βNTr(X†X), where β = 1,2 or 4 r×espectively when the entries are real, complex ∝ and quaternionic. The eigenvalues x ’s of W are all non-negative real numbers whose joint distribution is known to i beoftheforminEq. (2)withV(x)=x/2 αlnx, whereα isaconstantgivenbyα=β(1+M N)/2N 1/N [30]. − − − Without any loss of generality, we consider the case where M N. The case N M can be easily analyzed from ≥ ≥ the result M N (see for instance the discussion in [31]). We consider the limit N and M with the ratio ≥ →∞ →∞ N/M =c 1 fixed. In this limit, the average spectral density is given by the Marˇcenko and Pastur law [32] (see Fig. ≤ 4(a) for the c=1 case) 1 (cid:112) ρ (x)= (x x )(x x), (16) mp − + 2πx − − 8 where x =(1 1/√c)2. Here we focus, for convenience, on the c=1 case where the density takes the form ± ± (cid:114) 1 4 x ρ (x)= − . (17) mp 2π x In this case, the eigenvalues are supported over the interval [0,4] with their center of mass located at x=1. TheobservableofinterestisagainthenumberN ofeigenvaluesbelongingtoanyinterval =[a,b]whereb>a 0. I I ≥ For the semi-infinite interval [a,+ ), the statistics of N was studied in Ref. [28] using the Coulomb gas method. I ∞ Here we study the statistics of N for any finite interval using a similar Coulomb gas method. For convenience, we I provide explicit results in the case where =[x,x+l]=[1,1+l], where l is the size of the interval. I 1.2 Theory NumericsforN=5000 1.0 ) 0.8 ) I x N (p ( 0.95 ρm Var 0.6 00..8950 ) 0.80 s ( 0.75 2 V 0.70 0.4 0.65 0.60 0.55 −5 −4 −3 −s2 −1 0 1 2 0 1 x 1+l 4 0.120−3 10−2 l10−1 100 (a) (b) Figure4. (a)AveragedensityofeigenvaluesfortheWishartensemble,whichistheMarˇcenko-Pasturdistribution,forthec=1 case and (b) Results for the Wishart number variance when I =[1,1+l]. The theoretical result is in equation (18). In this case, the rate function Ψ(W)(k ) associated with the distribution of N is given in equation (104) and I [1,1+l] the number variance is given by  2 ln(cid:0)Nl(3 l)3/4(cid:1), N−1 <l<3 N−2/3, βπ2 − − Var(N[1,1+l])≈β11π2 llnnNN ++T21V(βl()s,), ll=>33,+ (N/4s)2/3, (18) βπ2 β where T (l) e−βNΦ(1+l), where Φ(w) is given in equation (117) and V (s) is given for β =2 in equation (81), as in β β ≈ the Gaussian case. The number variance is compared to numerical results in figure 4(b). As in the Gaussian case, we notice that beyond l=3, where the right edge of the chosen interval coincides with the edge of the Marˇcenko-Pastur distribution at x=4 in Eq. (17), the variance again drops sharply as a function of l and approaches a constant value 1 lnN when l , in agreement with the calculations in Ref. [28] for the semi-infinite interval [1,+ ). As in the βπ2 →∞ ∞ Gaussian case, the variance also has a non-monotonic behavior with a maximum at l(cid:63) =12/7. β-Cauchy ensemble: Another well studied ensemble in RMT is the so-called β-Cauchy ensemble. In this case an N N square matrix X whose entries are distributed via Pr(X) det(1 + X2)−β(N−1)/2−1, where 1 N N × ∝ is the N N identity matrix and β = 1,2 or 4 respectively for real, complex or quaternionic entries. The × eigenvalues of X are real and are distributed over the full real axis with a joint density given in Eq. (2) with V(x) = ((N 1)/(2N)+1/(βN))ln(1+x2). The average density of eigenvalues is well known [17] and is given, for − any N, by 1 1 ρ (x) ρ (x)= . (19) (cid:104) N (cid:105)≡ ca π1+x2 Here again, we are interested in the statistics of the number of eigenvalues N for an interval . For the semi-infinite I I interval[0,+ ), thestatisticsofN wasstudiedforlargeN usingtheCoulombgasmethod[33]. Here, weconsidera I ∞ 9 finitesymmetricinterval[ L,L]. Theratefunctionψ(C)(k )associatedtothedistributionofN isgivenbyequation I I − (135). For the variance, we find (cid:18) (cid:19) 2 NL Var(N ) ln , N−1 <L<N. (20) [−L,L] ≈ βπ2 1+L2 The variance rises to its maximum value at L(cid:63) =1, and falls to zero as L increases further. We note that the average position of the largest eigenvalue, L N, behaves, for the purpose of the behavior of the number variance, as an ∼ effective edge for the Cauchy ensemble, even though its average density (19) has no edge. Results are shown in figure 5(b). 1.0 NumericsforN=500 Theory 0.8 ) 0.6 ) I x N ( a ( c r ρ a V 0.4 0.2 0.0 L 0 L 10−2 10−1 100 101 102 103 104 − x L (a) (b) Figure 5. (a) Average density of Cauchy ensemble, which is the Cauchy distribution and (b) Results for the Cauchy number variance when I =[−L,L]. The theoretical result is equation (20). Westartbyconsideringtheβ-Gaussiancaseindetail. TheformalismbasedontheCoulombgastechniquecombined with the resolvent method is presented for this case first and later applied to the other ensembles. III. β-GAUSSIAN ENSEMBLE For Gaussian ensembles, the joint distribution is given by Eq. (2) with V(x) = x2/2 and Dyson index β = 1,2,4 depending on whether the entries are real, complex or quaternionic. However, the same joint distribution with arbitrary β > 0 also appears in a class of tri-diagonal matrices found by Dumitriu and Edelman [18]. Therefore, we will consider a general β > 0. Below we first develop the Coulomb gas method for β-Gaussian ensembles (valid for all β >0). The same method will be applicable to the other two cases. A. The Coulomb gas Weconsidertheinterval =[a,b]forV(x)=x2/2. OurgoalistoderivetheprobabilitydensityofN ,thenumber I I of eigenvalues falling inside . It reads by definition I Pβ(G)(NI =kIN)= Z1 (cid:90) (cid:89)N dxke−βN(cid:80)Nk=1 x22k (cid:89)|xi−xj|βδ(cid:32)kIN −(cid:88)N 1I(xl)(cid:33) N,β k=1 i>j l=1 Z (k ) N,β I = , (21) Z N,β where 1 (x) is the indicator function for the interval , defined as 1 when x and zero otherwise. I I ∈I 10 Usingtheexponentialrepresentationofthedelta, itcanbewrittenastheratiooftwocanonicalpartitionfunctions (cid:90) ∞ (cid:89)N (cid:90) i∞ dµ (cid:90) ∞ (cid:89)N Z (k )= dx e−βE[x,µ,NI], and Z = dx e−βE[x,NI] , (22) N,β I k N,β k 2π −∞k=1 −i∞ −∞k=1 with energy functions (cid:88)N x2 (cid:88) (cid:32)(cid:88)N (cid:33) E[x,µ,N ]=N i ln x x +µ 1 (x ) k N , (23) I i j I l I 2 − | − | − i=1 i>j l=1 (cid:88)N x2 (cid:88) E[x,N ]=N i ln x x . (24) I i j 2 − | − | i=1 i>j Equation (24) can be interpreted as the energy of a Coulomb gas of N charged particles on a line, each subjected to an external harmonic potential V(x) = x2/2 and repelling each other via a logarithmic interaction. Similarly, Eq. (23) represents another Coulomb gas which is subjected to an additional constraint, the µ-dependent term, that prescribes the fraction of particles (i.e., eigenvalues) inside the box . I We introduce now the two-point density 1 (cid:88) ρ (x,x(cid:48))= δ(x x )δ(x(cid:48) x ) . (25) 2 i j N(N 1) − − − i(cid:54)=j This function and the one-point density (4) ρ (x) (hereafter denoted simply by ρ(x)) allow us to write the energy in N integral form, using the identities (cid:90) (cid:88) f(x )=N dxρ(x)f(x) , (26) i i (cid:90)(cid:90) (cid:88) N(N 1) f(x x )= − dxdx(cid:48)ρ (x,x(cid:48))f(x x(cid:48) ) . (27) i j 2 | − | 2 | − | i<j For the energy E[x,µ,N ] we obtain I (cid:90) +∞ x2 N(N 1)(cid:90)(cid:90) +∞ E[x,µ,N ]=N2 ρ(x) dx − dxdx(cid:48)ρ (x,x(cid:48))ln x x(cid:48) I 2 2 − 2 | − | -∞ −∞ (cid:32) (cid:90) b (cid:33) (cid:18) (cid:90) +∞ (cid:19) +µ N ρ(x)dx N +η N ρ(x)dx N , (28) I − − a -∞ where a supplementary Lagrange multiplier η was introduced to enforce the normalization condition of the density ρ(x). In the large-N limit, we may replace the two-point density ρ (x,x(cid:48)) by the product of one-point densities of x 2 and x(cid:48): ρ (x,x(cid:48))= N ρ(x)ρ(x(cid:48)) 1 ρ(x)δ(x x(cid:48)) ρ(x)ρ(x(cid:48)). Rescaling the Lagrange multipliers with N yields 2 N−1 − N−1 − ∼ eventually E[x,µ,N ]N(cid:29)1N2S(G)[ρ,µ], (29) I ∼ where (cid:90) +∞ x2 1(cid:90)(cid:90) +∞ S(G)[ρ,µ]= ρ(x) dx dxdx(cid:48)ρ(x)ρ(x(cid:48))ln x x(cid:48) 2 − 2 | − | -∞ −∞ (cid:32)(cid:90) b (cid:33) (cid:18)(cid:90) +∞ (cid:19) +µ ρ(x)dx k +η ρ(x)dx 1 (30) I − − a -∞ is called the action.