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 kth smallest eigenvalue λ (k) of a large dimensional Gaussian random matrix. We find that this probability behaves for large N as, where β is the Dyson index of the ensemble. The rate function ψ(c, x), computed explicitly as a function of x in terms of the intensive label c = k/N , 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.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 (T = 0) many-body wavefunction and the statistics of eigenvalues of a certain class of random matrices. More precisely, in presence of a harmonic potential U (x) =