Spectral order statistics of Gaussian random matrices: large deviations for trapped fermions and associated phase transitions

# Spectral order statistics of Gaussian random matrices: large deviations for trapped fermions and associated phase transitions

## Abstract

We compute the full order statistics of a one-dimensional gas of fermions in a harmonic trap at zero temperature, including its large deviation tails. The problem amounts to computing the probability distribution of the th smallest eigenvalue of a large dimensional Gaussian random matrix. We find that this probability behaves for large as , where is the Dyson index of the ensemble. The rate function , computed explicitly as a function of in terms of the intensive label , has a quadratic behavior modulated by a weak logarithmic singularity at its minimum. This is shown to be related to phase transitions in the associated Coulomb gas problem. The connection with statistics of extreme eigenvalues of random matrices is also elucidated.

\DeclareBoldMathCommand\boldlangle\DeclareBoldMathCommand\boldrangle

Introduction - Recent spectacular progresses in the fabrication of devices for the manipulation of cold atoms [[1], [2]] have given a formidable impulse towards the theoretical understanding of the behavior of quantum many-body systems using tools and techniques borrowed from classical statistical physics.
Perhaps the simplest conceivable setting is a system of fermions confined by optical laser traps into a limited region of space [[2], [3], [4], [5], [6], [7], [8]]. In this experimentally rather common setup, there is a well-known connection between the ground state many-body wavefunction and the statistics of eigenvalues of a certain class of random matrices.
More precisely, in presence of a harmonic potential , the single particle eigenfunctions are given by , where are Hermite polynomials and . For a non-interacting Fermi gas of particles in a 1D harmonic trap (in the same spin state) the spatial ground-state manybody wavefunction is the Slater determinant of the first eigenstates. The explicit evaluation of this determinant yields

 |Ψ0(x)|2=CNe−mωℏ∑Ni=1x2i∏j

where are the positions of the particles on the line and is the normalization constant such that . As firstly noticed in [[9]], (1) can be interpreted as the joint probability density (jpd) of the real eigenvalues of a matrix belonging to the Gaussian Unitary Ensemble (GUE) [[10]]. The general jpd for Gaussian matrices is indeed written as

 Pβ(λ)=1ZN,βe−βN2∑Ni=1λ2i∏j

with the Dyson index of the ensemble ( for GUE) and a normalization constant. The jpd of real eigenvalues (2) can be written as

 Pβ(λ)=e−βE(λ)ZN,β ,E(λ)=N2N∑i=1λ2i−12∑j≠kln|λj−λk| . (3)

Written in this form, (2) can be interpreted as the canonical Boltzmann-Gibbs weight of an associated thermodynamical system: a gas of charged particles confined to a line, and in equilibrium at inverse temperature under competing interactions, the quadratic confining potential and a logarithmic (2D-Coulomb) repulsive interaction. Then, through Random Matrix Theory (RMT), observables of a (quantum) non-interacting 1D Fermi gas at zero temperature are mapped to thermodynamical properties of a (classical) 2D-Coulomb gas confined on the line at temperature . In the Fermi gas picture, the logarithmic repulsion of the Coulomb gas is essentially a manifestation of the exchange interaction between fermions.
In this Letter, we employ this mapping to study the order statistics of trapped fermions, i.e. the distribution of the th leftmost particle of a Fermi gas on the line. In the RMT language, this amounts to computing the full statistics of the th smallest eigenvalue, including its large deviation (LD) tails, of large Gaussian matrices (for all ). Results on the order statistics abound for i.i.d. random variables [[23]], but are much scarcer for correlated variables (see [[24], [25], [26], [27], [28]] for most interesting recent developments). Our approach provides a unified and transparent framework to probe interesting features of these distributions, such as non-Gaussian tails and a non-analytic behavior at the peak, which is a direct consequence of phase transitions of the underlying Coulomb gas system. We find that the spectral order statistics are governed by a rate function that is a function of two variables, which we are able to compute explicitly, with the “intensive” label of the th eigenvalue. Quite remarkably, by taking different limits, this single function is able to recover at once i.) the law of small fluctuations of the th smallest eigenvalue [[29], [30]] in the bulk, ii.) the distribution of the number of eigenvalues lying to the left of a barrier at (), and iii.) known results on left and right LD for the extreme eigenvalues [[32], [33], [34], [35]].
Setting and summary of results - We consider -Gaussian ensembles ( being the Dyson index) of random matrices (see [[10], [31]])) whose eigenvalues have law (2). Using the 2D-Coulomb gas picture, it is well known [[10]] that, for large , the average density of eigenvalues (normalized to unity) of (2) approaches the celebrated Wigner’s semicircle law on the single support , . We now arrange the eigenvalues in increasing order and we denote them by . It is easy to notice that the cumulative distribution of the th smallest eigenvalue is related in a simple way to the tail-cumulative distribution of (the number of eigenvalues smaller than ),

 P[λ(k)k] . (4)

Indeed, the th eigenvalue is smaller than if and only if there are at least eigenvalues to the left of . Using (2), the probability density of is

 P[Nx=cN]=∫RNdλPβ(λ)δ(cN−N∑i=1Θ(x−λi)) . (5)

The full statistics of was computed for large in [[17], [18]], while for a symmetric interval around the origin , the full distribution of (including large deviation tails) and its variance in all regimes of have been computed in [[11]] (see [[13], [14], [15], [16], [19], [20], [21]] for related results). These results are in agreement for with numerical and analytical estimates provided in recent literature on 1D Fermi systems [[22]].
Computing the large behavior of (5) and using (4), we establish the following LD estimates

 P[Nx=cN]≈e−βN2ψ(c,x);P[λ(k)=x]≈e−βN2ψ(k/N,x) , (6)

where stands for logarithmic equivalence. The central result of our paper is the calculation of the large deviation (or rate) function (with and ), given explicitly in Eq. (15) and plotted in Fig. 1. Interestingly, the same function (viewed in turn as a function of or of , the other variable acting as a parameter) governs both the LD behavior of the cumulative distribution of and of the th eigenvalue . We will show below that, as an additional bonus, it also recovers in certain limits the right and left LD laws for the extreme eigenvalues of Gaussian matrices [[32], [33], [34]].

The rate function is convex, it satisfies and has a minimum (zero) on the red line in Fig. 1, identified by the equation

 c⋆(x)=∫x−∞dzρsc(z)=1−θx−sinθx2π, (7)

with for , while (resp. ) if (resp. ). Thus e.g. the distribution of is peaked around which is precisely its mean value for large , while the distribution of is peaked around , where and its generalized inverse are given by (7).
Asymptotics - Expanding the rate function (15) around its minima in the two directions, we get access to the law of typical fluctuations of and . We find

 Missing or unrecognized delimiter for \right (8)

where we have denoted and is a constant independent of and . The fact that is not simply harmonic in and close to its minima but contains also a logarithmic contribution has important consequences on the variance of the two random variables, as it is customary in this type of problems [[11], [17], [18]]. Indeed, inserting this behavior back into (6) [[36]], we find that the extensive variables and in the bulk have Gaussian fluctuations

 P[Nx=cN]≃e−(Nc−Nc⋆(x))22Δ1;Δ1=1βπ2lnρ3⋆N (9) P[yk=y]≃e−(y−√Nx⋆(c))22Δ2;Δ2=1βπ2lnNρ2⋆N, (10)

up to corrections for large , in agreement with earlier results in both cases [[29], [30], [17], [18], [37]]. These estimates are respectively valid for such that , i.e. for “deeply” within the semicircle edges, and , i.e. . However, for larger fluctuations, the Gaussian form is no longer valid, and the distribution has non-Gaussian LD tails described by the rate function (15). Furthermore, the limit (or symmetrically ) of (15) at fixed (corresponding to the largest or smallest eigenvalue [[38], [39]]) is interesting. We find

 limc→1−ψ(c,x)1−c =Ψ+(x),x>√2 (11) limc→1−ψ(c,x) =Ψ−(x),x<√2 , (12)

where ([[34]], r.h.s. of Eq. 13 with ) and ([[33]], Eq. 59) are respectively the right and left LD functions for the largest eigenvalue of Gaussian matrices. The different scalings in for the two branches of the largest eigenvalue ( on the left and on the right) can be recovered using a simple argument (see [[11]]). We find particularly interesting (as already discovered for the number statistics problem in a bounded interval [[11]]) that a rate function obtained with a Coulomb gas calculation is able to recover not only the “Coulomb gas” branch , but also the “energetic” (or instantonic) branch , which was originally calculated using an entirely different approach [[34]] (see also [[40]] for a subsequent elaboration in the context of quantum gravity and string theory). In next section, we sketch the derivation of the rate function (15).

Derivation - In order to evaluate the multiple integral (5) for large , we resort to the “constrained” Coulomb gas method (introduced in [[32], [33]] and subsequently used in many different problems [[41], [42], [43], [44], [45], [46], [47], [48], [20], [21], [49]]), where a multiple integral over an eigenvalue jpd (like (5)) is converted into the partition function of the associated 2D-Coulomb gas with density , in equilibrium at inverse temperature . Skipping details [[36]], the multiple integral (5) can be written as a functional integral over , , where

 E[ρ]=12∫dλρ(λ)λ2−12∬dλdλ′ρ(λ)ρ(λ′)ln|λ−λ′| +A1(∫dλρ(λ)Θ(x−λ)−c)+A2(∫dλρ(λ)−1) (13)

is the action (the integral version of the energy function ) and and are Lagrange multipliers constraining a fraction of eigenvalues to the left of , and enforcing normalization of the density to . The functional integral will then be dominated by the equilibrium density of eigenvalues (Coulomb gas particles) minimizing the action (13).
Using the resolvent method [[36]], we find that the equilibrium (saddle-point) density has the general form

 ρ⋆c,x(λ)=1π√(λ+−λ)(λ0−λ)(λ−λ−)x−λ (14)

for such that the radicand is positive. The edge points (depending parametrically on ) are the three roots of the polynomial , and the function is determined by the solution of the constraint . It turns out that, depending on the values of and , five different phases of the gas of eigenvalues are possible (see Fig. 2). For a generic position of the barrier , whenever we recover the semicircle law (phase ). For (resp. ) and (resp. ) we have a single support solution (phases and , respectively). For all the other values the solution has double support (phases and ) and generalizes the density in [[17], [18]] to the case . One can see that the double-support phases and are separated by the unperturbed phase (the semicircle law) that lies on the curve . This tiny separation in the plane between the double support phases is ultimately responsible for the weak non-analytic behavior of the rate function close to its minimum in both the and directions (see (8)).
Once the saddle-point density is known, we precisely find the law (6) with (where comes from the normalization factor in (2) and needs to be subtracted), where

 ψ(c,x)=∫c⋆(x)cdc′(λ20(c′)−x22+∫xλ0(c′)dzS(z,c′)), (15)

and . The rate function in (15) can be written in closed form in terms of elliptic integrals [[36]]. We have checked that for we recover the rate function for the so called index problem [[17], [18]], and all the results have been checked numerically with excellent agreement [[36]].

Conclusions - In summary, we provided a complete characterization of the spectral order statistics of large dimensional Gaussian random matrices. We showed that the problem is amenable to a Coulomb gas treatment through a simple relation between the distribution of the th smallest eigenvalue and the number of eigenvalues smaller than a threshold . Details on the derivations, numerical checks and the outlook for future research will be provided elsewhere [[36]]. In the future, it will be interesting to study the order statistics for other ensembles and the crowding effects close to a specific eigenvalue in the bulk (see [[50]] for the first eigenvalue of GUE), as well as to investigate finite corrections. In the Fermi gas picture (), our results provide the full statistics of particle number on a semi-infinite line (extending recent results [[11], [17], [18]]) and single-particle fluctuations in the bulk of a system of 1D fermions in a harmonic trap, including their LD tails. We expect the results presented here to apply also to bosonic systems in presence of very strong repulsive interactions between bosons [[51]]. The link between 1D Fermi gases and RMT makes it possible to speculate about the possibility of “simulating” RMT eigenvalues in laboratory by manipulating Fermi gas systems, and the experimental verification of LD results is an intriguing challenge whose realization is very much called for.

Acknowledgments - IPC acknowledges hospitality by the LPTMS (Univ. Paris Sud) where this work was initiated. He also acknowledges multiple discussions with Pierpaolo Vivo and Fabio Cunden as well as the work done on the manuscript.

### References

1. I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
2. S. Giorgini, L. P. Pitaevski, and S. Stringari, Rev. Mod. Phys. 80, 1215 (2008).
3. P. Calabrese, M. Mintchev, and E. Vicari, Phys. Rev. Lett. 107, 020601 (2011); J. Stat. Mech. P09028 (2011).
4. E. Vicari, Phys. Rev. A 85, 062104 (2012).
5. M. Campostrini and E. Vicari, Phys. Rev. A 82, 063636 (2010).
6. A. Angelone, M. Campostrini, and E. Vicari, Phys. Rev. A 89, 023635 (2014).
7. V. Eisler, Phys. Rev. Lett. 111, 080402 (2013).
8. V. Eisler and Z. Racz, Phys. Rev. Lett. 110, 060602 (2013).
9. E. Brézin, C. Itzykson, G. Parisi, and J. B. Zuber, Comm. Math. Phys. 59, 35 (1978).
10. M. L. Mehta, Random Matrices (Academic Press, Boston, 1991).
11. R. Marino, S. N. Majumdar, G. Schehr, and P. Vivo, Phys. Rev. Lett. 112, 254101 (2014).
12. N. S. Witte, P. J. Forrester, and C. M. Cosgrove, Nonlinearity 13, 1439 (2000).
13. O. Costin and J. L. Lebowitz, Phys. Rev. Lett. 75, 69 (1995).
14. M. M. Fogler and B. I. Shklovskii, Phys. Rev. Lett. 74, 3312 (1995).
15. A. Soshnikov, J. Stat. Phys. 100, 491 (2000).
16. P. J. Forrester and J. L. Lebowitz, Preprint [arXiv:1311.7126] (2013).
17. S. N. Majumdar, C. Nadal, A. Scardicchio, and P. Vivo, Phys. Rev. Lett. 103, 220603 (2009).
18. S. N. Majumdar, C. Nadal, A. Scardicchio, and P. Vivo, Phys. Rev. E 83, 041105 (2011).
19. N. S. Witte and P. J. Forrester, Random Matrices: Theory Appl. 01, 1250010 (2012).
20. R. Marino, S. N. Majumdar, G. Schehr, and P. Vivo, J. Phys. A: Math. Theor. 47, 055001 (2014).
21. S. N. Majumdar and P. Vivo, Phys. Rev. Lett. 108, 200601 (2012).
22. V. Eisler and I. Peschel, J. Stat. Mech. P04005 (2014).
23. H. A. David and H. N. Nagaraja, Order Statistics (Wiley, New Jersey, 2003), 3rd ed.
24. E. Brunet and B. Derrida, Europhys. Lett. 87, 60010 (2009); J. Stat. Phys. 143, 420 (2011).
25. N. R. Moloney, K. Ozogány, and Z. Rácz, Phys. Rev. E 84, 061101 (2011).
26. G. Schehr and S. N. Majumdar, Phys. Rev. Lett. 108, 040601 (2012).
27. K. Ramola, S. N. Majumdar, and G. Schehr, Phys. Rev. Lett. 112, 210602 (2014).
28. G. Schehr and S. N. Majumdar, First-Passage Phenomena and Their Applications, Eds. R. Metzler, G. Oshanin, S. Redner (World Scientific, 2013), Preprint [arXiv:1305.0639] (2013).
29. J. Gustavsson, Ann. Inst. H. Poincare Probab. Stat. 41(2), 151 (2005).
30. S. O’Rourke, J. Stat. Phys. 138, 1045 (2010).
31. I. Dumitriu and A. Edelman, J. Math. Phys. 43, 5830 (2002).
32. D. S. Dean and S. N. Majumdar, Phys. Rev. Lett. 97, 160201 (2006).
33. D. S. Dean and S. N. Majumdar, Phys. Rev. E 77, 041108 (2008).
34. S. N. Majumdar and M. Vergassola, Phys. Rev. Lett. 102, 060601 (2009).
35. G. Borot, B. Eynard, S. N. Majumdar, and C. Nadal, J. Stat. Mech. P11024 (2011).
36. Details will be published elsewhere.
37. O. Bohigas and M. P. Pato, J. Phys. A: Math. Theor. 43, 365001 (2010).
38. C. A. Tracy and H. Widom, Commun. Math. Phys. 159, 151 (1994).
39. C. A. Tracy and H. Widom, Commun. Math. Phys. 177, 727 (1996).
40. M. R. Atkin and S. Zohren, Journal of High Energy Physics 04, 118 (2014).
41. P. Vivo, S. N. Majumdar, and O. Bohigas, J. Phys. A: Math. Theor. 40, 4317 (2007).
42. P. Vivo, S. N. Majumdar, and O. Bohigas, Phys. Rev. Lett. 101, 216809 (2008).
43. P. Facchi, U. Marzolino, G. Parisi, S. Pascazio, and A. Scardicchio, Phys. Rev. Lett. 101, 050502 (2008).
44. E. Katzav and I. Pérez Castillo, Phys. Rev. E 82, 040104(R) (2010).
45. C. Nadal, S. N. Majumdar, and M. Vergassola, Phys. Rev. Lett. 104, 110501 (2010).
46. C. Nadal, S. N. Majumdar, and M. Vergassola, J. Stat. Phys. 142, 403 (2011).
47. Y. V. Fyodorov and C. Nadal, Phys. Rev. Lett. 109, 167203 (2012).
48. S. N. Majumdar, G. Schehr, P. Vivo, and D. Villamaina, J. Phys. A: Math. Theor. 46 022001 (2013).
49. For a recent review see S. N. Majumdar and G. Schehr, J. Stat. Mech. P01012 (2014).
50. A. Perret and G. Schehr, J. Stat. Phys. 10.1007/s10955-014-1044-5, pp. 1-34 (2014).
51. E. H. Lieb and W. Liniger, Phys. Rev. 130, 1605 (1963); M. Girardeau, J. Math. Phys. (NY) 1, 516 (1960); Phys. Rev. 139, B500 (1965).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters