Fluctuations of observables for free fermions in a harmonic trap at finite temperature

We study a system of 1D noninteracting spinless fermions in a confining trap at finite temperature. We first derive a useful and general relation for the fluctuations of the occupation numbers valid for arbitrary confining trap, as well as for both canonical and grand canonical ensembles. Using this relation, we obtain compact expressions, in the case of the harmonic trap, for the variance of certain observables of the form of sums of a function of the fermions' positions, $\mathcal{L}=\sum_n h(x_n)$. Such observables are also called linear statistics of the positions. As anticipated, we demonstrate explicitly that these fluctuations do depend on the ensemble in the thermodynamic limit, as opposed to averaged quantities, which are ensemble independent. We have applied our general formalism to compute the fluctuations of the number of fermions $\mathcal{N}_+$ on the positive axis at finite temperature. Our analytical results are compared to numerical simulations. We discuss the universality of the results with respect to the nature of the confinement.


Introduction
The recent experimental progresses in cold atoms [1][2][3], which have made accessible new types of observables, have stimulated a renewed interest in the study of fermionic systems over the past few years. Although some emphasis has been put on many body physics, many physical aspects of the problem are captured by a simple non-interacting picture. Moreover, there are practical ways to reach the non-interacting regime experimentally [1,2,4]. Here, we will restrict ourselves to this case, and consider a system of N fermions in a one dimensional trap described by the HamiltonianĤ For specific choices of the confining potential V (x), the positions of the fermions at zero temperature can be mapped onto the eigenvalues of random matrices. This relation is based on the fact that the ground state of this system takes the form of a Slater determinant: where ψ k (x) is the one-body eigenfunction of energy ε k , with k ∈ N. For example, in the case of a harmonic confinement V (x) = 1 2 mω 2 x 2 , the joint distribution of the positions, given by the modulus square of the many-body wave function, reads: where Z N is a normalisation constant. It is exactly the distribution of the eigenvalues of matrices in the Gaussian unitary ensemble (GUE) [5,6]. Similarly, in the case of an infinite square well, the distribution of the positions can be mapped onto the Jacobi unitary ensemble (JUE) [7,8]. This connection to random matrices has allowed to study many properties of the ground state, like the density, the correlations, the number fluctuations and entanglement entropy [9][10][11][12]. It is interesting to remark that in the cold atom literature [13][14][15], some global results were derived using the local density approximation (LDA), i.e. the Thomas-Fermi approximation, without realizing the connection to RMT. A remarkable recent achievement is the development of Fermi quantum microscopes [16][17][18], allowing the direct measurement of the fermions' positions in a confining trap. Motivated by this context, several theoretical studies have focused on different observables counting the fermions in a given spatial domain (see e.g. Refs. [9,10,12,19]). One such observable is the number N + of fermions on the positive axis: Thanks to the mapping to the random matrix theory (RMT), the number N + of fermions in the domain x > 0 at zero temperature is precisely the number of positive eigenvalues of GUE random matrices. In the RMT literature, this number is known as the index, and its statistical properties have been studied for various RMT ensembles, including the GUE [20,21] and the Cauchy ensemble [22]. For GUE, the mean value of N + is trivially N/2, but the variance is nontrivial, given by where γ the Euler-Mascheroni constant. In the context of fermions, this finite value characterizes the quantum fluctuations. The logarithmic behaviour can be related to the anticorrelation of fermions discussed below, see Eq. (29). However in most experiments, the measurements are done at low but finite temperature. The zero temperature results obtained from RMT are therefore not sufficient to address these finite temperature properties. Our goal here will be to characterise the effect of thermal fluctuations on the variance of N + and to generalise the result (5) to finite temperature. Statistical physics provides several tools to analyse thermal fluctuations. When quantum correlations are dominant, like for the problem we aim to study, it is well known that the most efficient approach is supplied by the grand canonical ensemble in which the temperature T and the chemical potential µ are fixed, while the energy and the number of fermions fluctuate. Many-body quantum eigenstates are conveniently labelled by a set of occupation numbers {n k }, where n k = 1 if the individual eigenstate ψ k (x) is occupied by one fermion and n k = 0 otherwise. The grand canonical weight is: where is the grand canonical partition function. β = 1/(k B T ) is the inverse temperature and ϕ = e βµ the fugacity. In atomic traps, the number of atoms is fixed, which is best described by the microcanonical or canonical ensembles. Furthermore, due to the evaporative cooling techniques, the number of atoms is only moderately large, N ∼ 10 4 to 10 7 , and the equivalence between statistical physics ensembles is questionable. This has recently motivated several works where basic quantities, such as occupation numbers, energy, specific heat, etc, have been analysed in the canonical and microcanonical ensembles (see for instance Refs. [23][24][25][26][27][28]). In the canonical ensemble, the number N of fermions is fixed (instead of the chemical potential) and the Gibbs weight is where is the canonical partition function. The constraint on the number of particles included in the distribution however makes the calculations much more difficult in practice. For large N , deviations from the thermodynamic limit, and thus differences between predictions from the ensembles, are supposed to be small. However, it is worth stressing that the equivalence of ensembles is valid only for thermodynamic quantities (averaged observables), and not for their fluctuations [29,30], which are our main interest here. In the present article, we will introduce a general formalism allowing to study the fluctuations of a wide class of observables of the form known as linear statistics of the positions of the fermions, where h is any given function, not necessarily linear. For example, the potential energy in a harmonic trap corresponds to h(x) = x 2 , whereas the index N + , under consideration here, corresponds to h(x) = Θ(x). Most theoretical studies of fluctuations at finite temperature were performed in the grand canonical ensemble [2,11,[31][32][33], with the exception of Ref. [28], where the specific form n x 2 n has allowed an exact calculation at all temperatures in both ensembles (see also Ref. [34]). Our aim here is to introduce a general framework to analyse the fluctuations of arbitrary linear statistics within both the grand canonical and canonical ensembles.

Summary of the main results
The fluctuations of the linear statistics L = n h(x n ) have two origins at finite temperature. For a given many-body quantum state | {n k } , labelled by the occupation numbers, the positions x n 's fluctuate due to quantum fluctuations. In addition, the occupation numbers n k 's themselves fluctuate at finite temperature -this is the thermal fluctuations. To characterize these thermal fluctuations, we have found a general relation for the correlator of occupation numbers: where · · · c,g denotes the thermal average (see also [35]). We stress that this relation (11) is a very general one, independent of the confining trap. Moreover it is valid both in the canonical (c) and grand canonical ensemble (g). Note that in the grand canonical case, where the mean occupation numbers are given respectively by the Bose-Einstein and Fermi-Dirac distributions, this relation leads to: n k n l g = n k g n l g . This is the well-known independence of energy levels. This relation (11) plays a crucial role for the derivation of our subsequent results.
We have studied a 1D system of N fermions in a harmonic trap V (x) = 1 2 mω 2 x 2 , where N is either fixed (canonical ensemble) or fluctuating (grand canonical ensemble). We are interested in the fluctuations of the number N + of particles in the domain x > 0. We can already identify two temperature scales, which will play a role below: -a quantum scale T Q = ω/k B : when T ∼ T Q , the small thermal energy allows only a few excitations near the Fermi level and the discreteness of the spectrum matters. We refer to this case as the quantum regime; -the Fermi temperature T F = N ω/k B : when T ∼ T F , the system is dominated by large thermal fluctuations. All energy levels contribute and the spectrum can be considered as continuous. Thus we call it the thermal regime. The quantum regime covers the transition between the regime where the spectrum should be considered as discrete (T T Q ) and the regime where it can be considered as continuous (T T Q ), while the thermal regime describes the crossover between the regime dominated by quantum fluctuations (T T F ) and the classical regime (T T F ) [36], cf. Fig.1. Note that, while in the canonical ensemble N is fixed, in the grand canonical ensemble it is a fluctuating random variable. Then, in the grand canonical ensemble, the Fermi temperature is defined as T F = N g ω/k B . Henceforth, we will denote the variance of the total number of particles by Var g (N ) (the properties of this variance, in particular its T -dependence, are discussed in Appendix A; it is plotted in Fig. 4 as a function of temperature). We denote Var(N + ) the variance of N + which is computed in these two different regimes and the different statistical ensembles.
In the quantum regime T ∼ T Q , we have obtained: where the zero temperature variance Var(N + )| T =0 is given by Eq. (5), and we introduced In the grand canonical case, the variance receives a contribution proportional to the fluctuations of the total number of particles, Var g (N ). This observation is also valid in the thermal regime T ∼ T F , where our result reads: where In the thermal regime, we have not included the contribution of the quantum fluctuations at zero temperature, as it is subdominant, of order ln N . We checked that these two regimes smoothly match together for T Q T T F . Plots of these different expressions for the variance are shown in Fig. 2.
We stress an important difference between the two scaling functions F Q and F T : while the latter is not universal, i.e. depends on the precise form of the confining potential, the function F Q is universal.

Outline of the paper
The paper is organised as follows: Section 2 introduces the linear statistics for fermions. A general formula for the variance of linear statistics is obtained. In Section 3 we derive universal relations on the occupation numbers, which are useful to derive our results for the variance. We check our formulae in Section 4 by recovering the results of Ref. [28] on the potential energy. The case of the index variance for fermions in a harmonic well is discussed in Section 5. We end the paper with some concluding remarks. Some technical details are relegated to Appendices.

Observables of the form of linear statistics
Consider a system of spinless and non-interacting fermions in a confining potential V (x). We denote ψ k (x) the one-particle wave-function associated to the energy ε k , with k ∈ N. These energy levels are non-degenerate in one dimension. In the absence of interaction, a many-body quantum state can be conveniently labelled by the set of occupation numbers {n k }, where n k = 0 or 1 for fermions. The total number of fermions is then N = k n k fermions. The associated many-body wave function takes the form of a Slater determinant combining the one-particle wave functions of the occupied levels: where {k i } i=1,...,N is the list of occupied levels: n k = 1 if k ∈ {k i }.
We will consider the situation where the fermions are subjected to thermal fluctuations. In the following, we describe these fluctuations either in the canonical or the grand canonical ensemble. Since the description of our system of fermions involves two different sets of random variables, the positions and the occupation numbers, we are led to define two types of averaging: -the "quantum average", denoted · · · | {n k } , which corresponds to averaging over the positions of the fermions. This averaging is defined for a given quantum state | {n k } , with fixed number N = k n k of particles, as: for any function F of p N positions. In the following we will often omit the subscript | {n k } for simplicity. -the "thermal average", which corresponds to averaging over the occupation numbers.
This averaging depends on the ensemble under consideration. We denote it · · · c in the canonical ensemble and · · · g in the grand canonical ensemble. For any function G of the occupation numbers, it is defined as where P c,g ({n k }) represents the canonical or grand canonical measures given by Eqs. (8) and (6) respectively. The quantum and thermal average will involve first averaging over the positions of the fermions {x n }, Eq. (17), and then over the occupation numbers, Eq. (18). At zero temperature, the system is frozen in its ground state Ψ 0 and only the quantum average remains: Such integrals can be evaluated by making use of the determinantal structure of Ψ 0 , as we will see in section 2.1. When going to finite temperature, the excited states contribute to the thermal averaging: The integral over the positions can still be computed by using the determinantal structure, but the summation over the quantum states makes the problem much more challenging.
In the following, we will study a wide class of observables which take the form of linear statistics L of positions of the fermions, Eq. (10). Our aim is to study the variance In order to compute this variance, we first need to evaluate the quantum averages L and L 2 in a given quantum state | {n k } . This is the object of the next section. The thermal averaging will be discussed in Section 3.

Quantum averages and determinantal structure
Let us first consider a quantum state | {n k } , which contains N = k n k particles. It is convenient to use the fact that the positions of the fermions, in that given state, is a determinantal point process [11,33]. This means that the modulus square of the many-body wave function can be rewritten as a determinant where we have introduced the kernel Since the wave functions are orthogonal, this kernel verifies the reproducibility property where we used that n 2 k = n k since n k = 0 or 1 for fermions. The direct consequence of this property is that the n-points correlation function also takes a simple determinantal form: In particular, the one-point function is given by When averaged over n k 's in the grand canonical ensemble, using Eq. (6), this is precisely the mean density of fermions where ρ(x) = n δ(x − x n ) and · · · g is the usual quantum and statistical average. The two-point correlation function reads: When averaged over n k 's in the grand canonical ensemble, this is related to the familiar relation for the density-density correlation function in the Fermi gas where the minus sign is related to the effective repulsion between fermions due to the Pauli principle. When the density can be considered constant, equal toρ, the zero temperature kernel is the famous sine-kernel Using these properties we can easily express the quantum average of the linear statistics (10): where we used the symmetry under the exchange of particles. Therefore, only the one-point function, Eq. (26), is needed for this computation: Using the expression of the kernel, Eq. (23), one can express this average in terms of the one-particle wave functions: We can similarly compute the mean square: where we have introduced the matrix elements Using these definitions, we can rewrite (33) as: This derivation shows that one only needs to compute the matrix elements A k,l and B k to perform the quantum averages involved in the variance of a linear statistics L.

A general formula for the variance
Using the results of the quantum average, Eqs. (34,36), one can straightforwardly take the thermal average in the canonical or grand canonical ensemble. Combining these two relations, we obtain a general expression for the variance of a linear statistics: where Cov c,g (n k , n l ) = n k n l c,g − n k c,g n l c,g . The difference between the two ensembles clearly appears on this relation. Indeed, the covariance of occupation numbers is zero in the grand canonical ensemble: Cov g (n k , n l ) = 0 for k = l. We will see that this term gives an additional contribution of the same order as the first two.
In order to use this formula for the variance, one first needs to compute the coefficients A k,l and B k , which will be analysed in Appendix C. Then, the sums over the levels need to be evaluated. An important piece of the calculation is a relation on the occupation numbers which will be derived in the next section.

Occupation numbers
In this section, we derive general relations involving the occupation numbers, valid both in the canonical and grand canonical ensembles. Although this article focuses on fermionic systems, we will also consider the bosonic case for the sake of generality.

Grand canonical ensemble
Let us first consider the case of the grand canonical ensemble. In this ensemble, the mean occupation numbers are the well-known Bose-Einstein and Fermi-Dirac distributions: for fermions (40) where µ is the chemical potential, which controls the mean number of particles The variance of the occupation numbers is also known: with the upper sign for bosons and the lower sign for fermions. The simplicity and the success of the grand canonical ensemble relies on the independence of individual energy levels, i.e. the absence of correlations between occupation numbers:

Canonical ensemble
In the canonical ensemble, the total number of particles is fixed to N . This constraint induces correlations between the occupation numbers associated to different levels, which are in general very difficult to handle. The mean canonical occupation numbers are obtained by averaging with the canonical measure (8): The general strategy of statistical physics is to introduce generating functions, i.e. consider the sum where we recognised the grand canonical measure, Eq. (6). We can then deduce n k c by a contour integral which selects the ϕ N term in Eq. (45): where the integrals run over a closed contour winding once around the origin in the counterclockwise direction. The partition functions Z c and Z g being related by we can rewrite the occupation numbers as A similar relation clearly holds for any average quantity, for instance We will first derive general expressions for (48,49). In a second step we will analyse the large N limit by a saddle point method. The two results will be useful in the next sections.

General representation for the covariance
To evaluate the numerator in Eq. (48), we only need to determine the coefficient of the term ϕ N in the expansion of Z g (ϕ)n k g as a power series in ϕ. This series can be obtained by using Eqs. (40,47). Then, isolating the term proportional to ϕ N yields: where the upper sign is for bosons, and the lower sign for fermions. This relation was known in the literature, see e.g. Refs. [26,37,38]. Similarly, we can obtain the expression for the product of occupation numbers: This last relation can be found in the case of bosons in Ref. [26]. It can be simplified by introducing s = p + q: Separating this expression into two sums, we recognise the expressions of n k c and n l c , Eq. (50), thus, We derived this relation in the canonical ensemble, but one can easily check that it also holds in the grand-canonical one, where it becomes n k n l g = n k g n l g . Therefore, this is a very general relation, valid for any system of non interacting bosons or fermions, either in the canonical or grand canonical ensemble. To the best of our knowledge, this relation (53) involving the occupation numbers is new.
We have extended this analysis to the p-point correlation functions in Ref. [35].

Saddle point estimate for large N
Despite having obtained general relations (50,51), their large N analysis remains a challenge. For this purpose, we perform a saddle point analysis of Eqs. (48,49). For all the integrals, the saddle point ϕ is given by the condition which fixes the fugacity (or the chemical potential) such that the grand canonical mean number of particles in the trap is equal to the fixed canonical number N . Then, the integral (48) yields: This indicates that the statistics of the occupation numbers are the same at leading order in N in the canonical and grand-canonical ensembles. This is not surprising since we expect both ensembles to be equivalent for averaged quantities in the thermodynamic limit. The equivalence of thermodynamic results however only holds for averages and not for variances or covariances, as it is well-known [29]. Performing the same saddle point analysis with Eq. (49) yields Since we will need the covariance of occupation numbers, we push the saddle point approximation of Eqs. (48,49) to the next order, see Appendix B, in order to get the O(N −1 ) corrections. After many simplifications, we obtain the compact expression: Cov c (n k , n l ) = − Var g (n k )Var g (n l ) p Var g (n p ) where Var g (n k ) is given by Eq. (42). Note that since Var g (n k ) = Var c (n k ) + O(N −1 ), we can express this covariance in terms of the variance of occupation numbers in any ensemble. We chose to express it in terms of the grand canonical one since the expressions are simpler in this case. In the canonical ensemble, N = k n k is fixed, which implies the sum rule It is clear that Eq. (57) is consistent with this sum rule, up to O(N 0 ). However relation (57) is valid as long as the covariance is a O(N −1 ) correction to the saddle point approximation. It is the case when the denominator, which corresponds to the variance of the total number of particles in the grand-canonical ensemble, is of order N : p Var g (n p ) = Var g (N ) = O(N ). As discussed in Appendix A, this is verified only in the thermal regime T ∼ T F . Therefore, this relation for the covariance can only be used this regime. In the quantum regime T ∼ T Q we will instead rely on the exact relation (53).

A symmetry relation for fermion occupation numbers
Let us now consider fermions in a harmonic trap V (x) = 1 2 mω 2 x 2 , which will be the main focus of the paper. In this case, the spectrum is linear, ε n = (n + 1 2 ) ω, n ∈ N.
In the grand canonical ensemble, the Fermi-Dirac distribution (40) has a the well-known symmetry around the chemical potential µ: which is usually interpreted as the particle-hole symmetry. In the case of a discrete spectrum, setting the chemical potential in the middle of a gap, µ = N f ω, where N f is an integer, the relation reduces to: Due to the linearity of the spectrum, we have N g = N f up to exponentially small correction In the canonical ensemble the occupation numbers (50), are expressed in terms of the canonical partition function. For the harmonic oscillator, it can be computed analytically, and the result can be found in a few textbooks [30,39]: Using this result, the expression of the occupation numbers (50) greatly simplifies near the Fermi level in the large N limit: with corrections exponentially small with N . Using this expression, it is straightforward to show that the canonical occupation numbers also exhibit the particle-hole symmetry around the Fermi level: which is exactly the same as the grand canonical relation (60). These symmetries (60,63), along with relations (53,57) on the occupation numbers will be essential for our study of linear statistics for fermions in a harmonic trap. We will first check our approach by recovering the recent results of Ref. [28] on the potential (or kinetic) energy of the fermions. Then, we will apply our method to the observable introduced in Section 1: the number N + of fermions on the positive axis, given by Eq. (4).

A first check: potential energy of fermions in a harmonic trap
In the harmonic trap, the one-particle wave functions are expressed in terms of Hermite polynomials: with energies ε n = (n + 1/2) ω, n ∈ N. We denote {x n } the positions of the fermions.
Recently, the distribution of the potential energy E p , or equivalently the kinetic energy, was obtained in Ref. [28] for any temperature T or fixed number N of fermions (canonical ensemble). However their method is restricted to the study of the potential energy which is a specific linear statistics (10), with h(x) = x 2 . They have obtained two different scaling functions describing the quantum and thermal regimes: and Var c (I) where the two scaling functions are where Li s (z) = ∞ k=0 z k /k s is the polylogarithm function. One can check that these two expressions (66,67) smoothly match in the intermediate regime T Q T T F = N T Q , as the two scaling function present the limiting behaviours V q (z) z for z → ∞ and V th (z) z for z → 0. In this section we will check using our more general approach that we recover these results.
In this particular case of linear statistics with h(x) = x 2 , the matrix elements (35) can be computed exactly:

Quantum regime
In this regime, the small temperature T ∼ T Q allows only a few excitations above the Fermi level ε N . Therefore, we expect the main contribution to come from the proximity of this level. At leading order in N , for fixed k and l, the matrix elements become: Using these expressions in Eq. (39), the variance of I becomes, at leading order: where we extended the summation to −∞ instead of −N , since the corrections are exponentially small. Using then relation (53), this becomes This last sum cannot be separated into two sums because they would both diverge. Using the symmetry of the mean occupation numbers around the Fermi level, see Eq. (63), we get Under this form, the sum can be separated into two sums. This yields Finally, using again Eq. (63) gives Therefore, we recover the result of Ref. [28] in the quantum regime, Eq. (66):

Thermal regime
We now consider the regime where the temperature is of the order of the Fermi temperature T ∼ T F . We fix y = βN ω = T F /T and let N → ∞. In this case, we can use the results of section 3.2.2, for instance where the chemical potential µ is fixed by k n k g = N . This sum over k can be replaced by an integral over x = k/N since the discreteness of the spectrum plays no role in this regime: Imposing that this sum is the total number of particles yields: 1 Therefore, we can rewrite the occupation numbers (80) as where we introduced the notation for the Fermi-Dirac distribution in terms of convenient variables. In most cases, the approximation (83) is sufficient. However, the corrections are essential when the fluctuations of the occupation numbers contribute. For instance, the covariances in the last term of Eq. (39) must be evaluated using Eq. (57). Replacing now the sums over k by integrals over x = k/N , and keeping only the leading order of the matrix elements (70,71) yields: Evaluating these integrals, we finally obtain Hence we recover the result of Ref. [28] in the thermal regime, Eq. (67). This verification validates our approach to compute the variance of linear statistics.

Discussion
It is interesting to comment on the physical content of these results, and in particular compare the fluctuations of the potential energy with the fluctuations of the total energy E = E c + E p , which, in the canonical ensemble, can be related to the heat capacity studied in detail in [30] by Var c (E) = k B T 2 C V . From (61), we get 1. A more precise treatement of the steepest descent equation k n k g = N , with the help of the Euler- Bp p! f (p−1) (0) with Bp the Bernoulli number (B1 = −1/2, B2 = 1/6, etc), shows that the zero temperature limit of the chemical potential is precisely in the middle of the energy gap above the last occupied level: µ = N ω + O(N −1 ) = (εN + εN−1)/2.
In the low temperature regime, the exponential suppression of the fluctuations can be related to the existence of a gap ω in the excitation spectrum. The result in the intermediate regime T Q T T F has been rewritten in terms of the Fermi temperature to make clear that it corresponds to the classical result Var c (E) N (k B T ) 2 (given by the equipartition theorem) multiplied by the small factor T /T F 1. This well-known suppression factor originates from the Pauli principle, which restricts thermal fluctuations to take place in a small window of width k B T around the Fermi level.
The fluctuations of the potential energy are given by Eq. (66) and (67): . This observation has interesting consequences for the correlations of kinetic and potential energies. We express the variance of the total energy E = E c + E p , and use the fact that E c and E p have the same statistical properties for harmonic confinement: In the classical regime T T F , we have obtained Var c (E p ) = Var c (E c ) (1/2)Var c (E), which is related to the well-known fact that the kinetic and potential energy are uncorrelated: Cov c (E c , E p ) 0. In the regime T T F , we have obtained that Var c (E p ) Var c (E), implying that potential and kinetic energies are anti-correlated so that the fluctuations can be related as δE p −δE c .

Index variance for fermions in a harmonic trap
We now apply the general considerations of sections 2 and 3 to the study of the index N + , corresponding to the number of fermions on the positive axis. It is given by Eq. (4). This quantity is a linear statistics: it is of the form (10), with h(x) = Θ(x). Therefore, we can use the results of section 2. In particular, the variance of N + is given by Eq. (39), with and This last relation is exact, due to the symmetry of the potential. Using this result, we can rewrite Eq. (39) as: The first term gives the mean number of particles. Moreover the second term has a simple structure thanks to the fact that the matrix elements B k and A k,k are equal and independent of k, in Eq. (92). Note that this property in Eq. (92) is specific to the choice of the observable considered here, namely the index N + . Writing N = k n k , the second term can be simply identified as the variance of the total number of particles where obviously Var c (N ) = 0 by definition, while Var g (N ) is finite. Thus: Before studying into detail the variance of N + at any temperature, we will first discuss the limit of high temperature in which the fermions behave as classical particles.

High temperature limit: the Maxwell-Boltzmann regime
We start by considering the simplest limiting case, the limit of high temperature T T F . In this case the fermions can be considered as classical particles as the thermal fluctuations dominate. Therefore, their positions {x n } are independent, and they follow the Maxwell-Boltzmann distribution: The probability that a particle is in the domain x > 0 is p + = 1 2 . We now need to distinguish the statistical ensembles: -In the canonical ensemble, the number N of particles in the trap is fixed. Therefore, the mean number of particles with position x n > 0 is: To compute the variance, we also need the square of N + : from which we deduce: From these results, we can deduce the variance: -Let us now consider the grand canonical ensemble in which the number N of fermions fluctuates. In this case, the expressions are simply obtained by averaging Eqs. (97,99) over N : From which we deduce: where we have introduced the variance of the total number of particles Var g (N ) = The properties of this variance and its temperature dependence are discussed in Appendix A. Let us remark that even though the relation between canonical and grand canonical variances (103) has been derived in the classical Maxwell-Boltzmann regime it actually turns out to be much more general, as we discuss in Subsection 5.4 below.
In the limit of high temperature T T F , the index variance thus reads: Var g (N + ) N g 2 grand canonical. (104)

Canonical ensemble
We have derived in the previous sections a general expression for Var(N + ), Eq. (95), which involves the matrix elements A k,l and occupation numbers. The coefficients A k,l are computed in Appendix C, and we studied the occupation numbers in section 3. In this section we will combine these results to derive a more explicit expression for the index variance of fermions in a harmonic trap Var c (N + ), in the canonical ensemble. In this case, since the total number N of fermions is fixed, Var c (N ) = 0 and Eq. (95) reduces to: We will compute this variance first in the quantum regime T ∼ T Q , then in the thermal regime T ∼ T F .

Quantum regime
This regime can be reached for finite β = 1/k B T by letting the number N of fermions become large. Since we already know the variance of N + at zero temperature, Eq. (5), we will focus on the difference between the variance at temperature T and the one at zero temperature T = 0: At T = 0, Eq. (105) is still valid, but with fixed occupation numbers. Only the N lowest energy levels are occupied, thus Therefore, we have from Eq. (105): Since the main differences between the occupation numbers at zero and finite temperature are visible near the Fermi level N − 1, we shift the indices in the sums to start the summation from the Fermi level: First, for large N , we can let the summations go to infinity, as the corrections are exponentially small with N . Then, since the coefficients A k,l given by Eq. (168) are non zero only if k and l have different parity, we get: We can then use Eq. (53) to evaluate the thermal averages of products of occupation numbers. Introducing k − l = 2n − 1 in the first sum, k + l = 2n − 2 in the second one, and making use of Eq. (63) many cancellations occur, yielding a compact expression: involving a universal function F Q , as argued below in Section 5.5. This is our final result for the variance in the quantum regime for the canonical ensemble. Remarkably, this formula for fermions involves Bose-Einstein factors, like in the mean total energy [30] or its variance (87). This observation has a simple origin: the system of fermions have particle-hole excitations of bosonic nature. This well-known fact is as the heart of the bosonization technique for 1D Fermi liquids (see [40] for a general reference and [23] for a discussion of bosonization in the presence of a harmonic well).
From our result (111), we can extract the asymptotic behaviours of the temperature dependent part of the variance: The low temperature behaviour can be simply associated with the existence of a gap ω in the excitation spectrum.

Thermal regime
We now consider to the thermal regime, where the temperature is of the order of the Fermi temperature T ∼ T F . We fix y = βN ω = T F /T and let N → ∞. In this case, we proceed as in Subsection 4.2 where we analysed the potential energy. The study of the index is however more simple thanks to the simplification mentioned in the beginning of Section 5, which makes the second term in parenthesis in Eq. (93) vanish in the canonical ensemble. As a result, in this case it is sufficient to replace the occupation numbers by the rescaled Fermi-Dirac distribution: n k c f y (k/N ) + O(N −1 ) where f y is given by Eq. (84): Let us first rewrite the double sum as: Since A k,l is non zero only if k and l have different parity, see Eq. (168), the sum over p involves only odd integers p = 2n − 1. Replacing the summation over k by an integral over x = k/N gives: Using this result in Eq. (113), along with yields where Var g (N ) is the variance of the total number of particles in the grand canonical ensemble, where N g must be replaced by N . This non trivial relation between Var c (N + ) and Var g (N ) relies on the specific properties of the matrix elements A k,l , thus on the nature of the observable N + . Using the expression of this variance from Appendix A, we obtain the final expression for the variance of the index in this regime: where the subleading T = 0 contribution has been omitted. This is our final result in the thermal regime for the canonical ensemble. From this general expression of the variance, we can extract its asymptotic behaviours as function of the temperature: First note that the high temperature limit T T F matches with the Maxwell-Boltzmann case, Eq. (100), as it should. In addition, the low temperature limit in this thermal regime, T T F , smoothly matches the high temperature limit from the quantum regime, T T Q , Eq. (112). This indicates that there is no intermediate regime of temperature between these two. The low temperature result can be simply understood as the classical result, N/4, reduced by the factor T /T F characteristic of a degenerate Fermi gas [30], as already mentioned for the potential energy.

Grand canonical ensemble
In the previous section we derived expressions for Var(N + ) in the canonical ensemble in both quantum and thermal regimes. We now perform a similar computation in the grand canonical ensemble. In this case, the chemical potential µ is fixed, while the number of fermions in the trap fluctuates. In order to easily compare the results between the two ensembles, we will use the mean number of particles N g as a parameter instead of the chemical potential. The two are related by Eq. (41). The mean occupation numbers n k g are given by the Fermi-Dirac distribution (40). Occupations are uncorrelated between different energy levels. This allows to rewrite the general expression (95) as where the variance of the total number of particles Var g (N ) is studied in Appendix A. As before, we will first discuss the quantum regime and then the thermal one.

Quantum regime
Again, we fix β = 1/(k B T ) and let N g → ∞. As before, we focus on the difference Using Eq. (120), we can express this as We evaluated the same double sums in section 5.2.1, using only relations (11) and (63) which hold in both ensembles. Therefore, our previous derivation is still valid, and we have: We obtain the final expression for the variance of the particle number: where F Q (ξ) is given by Eq. (111). The variance in the grand canonical ensemble is thus obtained from the canonical one by adding a term proportional to the variance of the total number of particles. Using the limiting behaviours of Var g (N ) given in Appendix A along with Eq. (112), we can straightforwardly deduce the asymptotic behaviours: We again obtain an essential singularity at zero temperature, but different from the canonical case (e −T Q /T vs e −T Q /2T ). This is due to the fact that in the grand-canonical case, the leading contribution comes from the term proportional to the total number of particles Var g (N ).

Thermal regime
We now fix y = N g β ω and let N g → ∞. As in the quantum regime, the last term in Eq. (120) was already computed in section 5.2.2, and is given by Eq. (115). Therefore, we straightforwardly obtain: where Var c (N + ) is here given by Eq. (118) with N substituted by N g . In the r.h.s., we have neglected the subleading T = 0 contribution, Eq. (5). In this regime, the variance Var g (N ) can be computed explicitly and is given by Eq. (151): Therefore, we can rewrite the index variance as In this thermal regime, the index variance takes twice its canonical value in the grand canonical case. We can thus straightforwardly obtain its asymptotic behaviours as a function of T : Again, we check that the low temperature limit T T F smoothly matches the limit T T Q in the quantum regime, Eq. (125). We also recover the Maxwell-Boltzmann limit, Eq. (104) in the high temperature limit, as expected, whereas the quantum regime again shows the usual reduction factor T /T F .

Relation between the canonical and the grand canonical variances
We stress here that the relation between the variances in the canonical and grand canonical ensembles is completely general: compare (100) Although this argument is not quite precise, in the case of the energy E of the system, it gives which is a well-known relation in statistical mechanics [29,30]. Applied to the case of N + , Eq. (130) becomes: which is clearly the relation obtained several times above (103,124,126).

Universality
All our discussion so far has focused on the example of the harmonic trap as in this case the one body wave-functions are known, which makes the computations more explicit. Some of our results can however be extended to any type of confining potential, assuming a single minimum for simplicity.
Let us first discuss the zero temperature result, Var(N + )| T =0 , given by Eq. (5) for the harmonic trap. A similar result has been obtained for a system of fermions confined in an infinite square well. In this case, the variance of the number of particles in the right half part of the trap can be found in Eq. (40) of Ref. [41]: The variations of N + correspond to particles crossing the origin, thus Var(N + )| T =0 measures the quantum fluctuations through the origin. The leading log-terms in (5) and (133) coincide, and thus is a bulk property (which can be related to (30)). However,the subleading constants are different. We thus conclude that they are sensitive to the precise form of the potential. Therefore, we expect only the leading log-term to be universal. At finite temperature, the fluctuations are controlled by the two scaling functions F Q and F T , for the quantum and thermal regimes respectively. These functions are determined by the occupation numbers n k c,g , which depend on the spectrum {ε n } and thus on the potential, and the matrix elements A k,l . These latter depend only on the values of the wave functions ψ n at the center of the trap, see Eq. (165). These can thus be derived by semiclassical methods, such as the WKB approximation, for any kind of potential and exhibit universal properties. In particular, one can show that the expression (168) of these matrix elements for large quantum numbers is universal.
In the quantum regime, the derivation of Sections 5.2.1 and 5.3.1 is based on this universal expression of the matrix elements A k,l and on two properties of the occupation numbers, Eqs. (53) and (63). The first relation is universal, as discussed in Section 3. The second property (63), which implies the symmetry of the occupation numbers around the Fermi level, holds only for a linear spectrum. However, any regular spectrum {ε n } can be linearised near the Fermi level and thus (63) is also universal in the vicinity of the Fermi level. The temperature scale T Q is then defined from the gap at the Fermi level: Therefore, the derivation of Sections 5.2.1 and 5.3.1 can be extended to any confining potential, and the scaling function F Q is thus universal, and given by expression (13).
In the thermal regime, controlled by the scale all the spectrum contributes to the variance Var(N + ). Therefore, we do not expect the function F T to be universal. This can be shown explicitly from the relation (117) which states that this scaling function is proportional to the variance of the total particle number in the grand canonical ensemble: F T (T F /T ) = Var g (N )/4. We do not know this function explicitly, however we show in Appendix A that its limiting behaviours are only controlled by the exponent governing the one body density of states ρ(ε) ∝ ε α−1 : which thus depends explicitly on α.

Numerical simulations
We now compare our results with numerical simulations. Generating numerically realisations of the positions of the fermions is quite difficult. But in the grand canonical ensemble, we can use the fact that the positions of the fermions form a determinantal point process [11,33,42], with kernel: where ψ k are the one-particle wave functions of the harmonic oscillator (64). We reproduce in Appendix D an algorithm described in Refs. [43,44] which allows to sample such point processes. Therefore, we performed simulations only in the grand canonical case.  We generate realisations of the positions of fermions and compute the corresponding value of N + for each of them. From the set of data, we compute numerically Var(N + ). We studied the quantum regime, with N g = 100. We did not investigate the thermal regime because the computational efforts increase with the temperature, as more energy levels contribute. This makes it difficult to obtain enough statistics in the thermal regime where T = O(N g ).
Our result in the grand canonical ensemble (124), combined with the previously known zero temperature expression (5) read: where F Q (β ω) is given by Eq. (111) and c is a constant, see Eq. (5). We used this expression to compare the numerical data to our analytical result. They show an excellent agreement, see Fig. 3.

Conclusion
In this paper, we have introduced a general method to study certain many body observables of the form of sums of one body observables of the positions of fermions in a confining trap.
We have applied our method to the study of the number N + of fermions on the positive axis, in the case of a harmonic well. We have obtained explicit expressions for the variance of these observables, both in the quantum regime T ∼ T Q = ω/k B and the thermal regime T ∼ T F = N ω/k B in terms of two scaling functions, one universal and the other not. We have shown that these expressions smoothly match. We note that the fluctuations of linear statistics were recently studied by Johansson and Lambert in the grand canonical ensemble [45]. They considered the case where each of the fermions positions x i 's are scaled like N −δ while the temperature scales like T ∼ N α . The results presented in this paper thus corresponds to δ = 0 (which in their terminology corresponds to "macroscopic scale"). In addition, here, we have studied the crossover between the quantum regime (α = 0 in the notation of Ref. [45]) and the thermal regime (α = 1).
We have emphasised the difference between the statistical ensemble by computing the variance of N + both in the canonical and grand-canonical ensembles. This difference is due to the fact that we have considered the fluctuations of a global observable, N + = n Θ(x n ). On the other, local quantities, such as the two point density-density correlation functions, are ensemble independent [11].
The computation of these fluctuations in the microcanonical ensemble is still an open question. Indeed, our derivation made extensive use of relation (11) which no longer holds in the microcanonical case. We have studied the first moments of the distribution of the observable N + at finite temperature. Computing the full distribution of N + would be an interesting but much more challenging question.
A Variance of the total particle number in the grand canonical ensemble for any confining potential As discussed in Section 3, in the grand-canonical ensemble the mean occupation numbers n k g are given by the Fermi-Dirac distribution, Eq. (40). In addition, occupations are uncorrelated: n k n l g = n k g n l g . This allows to easily evaluate the first moments of the fluctuating total number of particles N = k n k : This first relation links the mean number of particles N g to the chemical potential µ. The variance reads: Cov g (n k , n l ) Therefore, we obtain: We can use these expressions to study the limiting behaviours of this variance for any confining potential of the form V (x) ∝ |x| p . The spectrum can be determined from a WKB approximation where x t is the turning point V (x t ) = ε n . This condition gives the scaling We can check that α = 1 for the harmonic potential while α = 1/2 corresponds to the infinite square well. In the continuum limit, this spectrum gives a density of states of the form where δ is an energy scale and A is a dimensionless parameter. For a trap containing on average N g particles, the Fermi energy is given by ε F = (αN g /A) 1/α δ. Therefore, we define the Fermi temperature as The temperature scale T Q can be defined from the gap at the Fermi level (134), which gives In the quantum regime, the spectrum can be linearised near the Fermi level: The variance Var g (N ) is thus universal in this regime. In the low temperature limit T T Q , Eq. (139) imposes that the chemical potential is fixed to the middle of the gap µ ε F + k B T Q /2. Using this value in Eq. (141), we can study the low temperature limit of the variance of N . The leading contribution comes from the two levels ε N g −1 and ε N g which are the closest to the chemical potential µ. This gives Var g (N ) 2 e −T Q /2T , T T Q .
In the thermal regime T ∼ T F , the sums can be replaced by integrals over the energy ε. Eq. (139) becomes: where we recall the Li α (z) = ∞ k=1 z k /k α . This last relation is more conveniently expressed in terms of dimensionless variables: −y −α Γ(α + 1) Li α (−e yμ ) = 1 , y = This relation fixes the rescaled chemical potentialμ in terms of the rescaled inverse temperature y. A similar computation for the variance (141) gives: These two relations (150) and (151) allow to plot this variance for different confining potentials, corresponding to different values of α, see Fig. 4, right. In the high-temperature limit, the variance reaches the classical limit The other limiting case is: which shows once again the reduction factor T /T F , as usual in the degenerate Fermi gas [30]. See

B Saddle point estimate
Let us consider integrals of the form: where C is any contour in the complex plane, g and φ are any given smooth functions. We want to estimate this integral in the limit of large N by using a saddle point method. The saddle point z is given by φ (z ) = 0 .
Let us make the change of variable z = z + t/ √ N and deform the contour C such that t is real. We can expand φ and g near z : Using these expansions, the integral I(N ) becomes: Computing the Gaussian integrals yields: We used this expression in section 3.2.2 to obtain the O(N −1 ) corrections to the covariance n k n l c − n k c n l c .

C Matrix elements
In this section we compute the coefficients A k,l associated to the index variance, Eq. (91), which are obtained from the quantum average. Since the index N + is a linear statistics (10) for h(x) = h(x) 2 = Θ(x), we can straightforwardly obtain the "diagonal terms" determinants involving the kernel K, see Eq. (25). We consider such a process with a kernel where 0 λ k 1 and A general method to generate numerically realisations of this process was introduced in Ref. [43]. We reproduce here a similar algorithm described in [44]. This algorithm was also used in the physics literature, see e.g. Refs. [47,48]. It relies on the following theorem: introduce a set of Bernoulli random variables n k = 0 or 1, with mean value n k = λ k . Then, the determinantal point process with kernel has the same statistics as the original process with kernel (169). In terms of fermions, this means that picking a realisation of the positions of the particles in the grand canonical ensemble is equivalent to first picking a quantum state {n k } from the Gibbs measure (6) and then generating a realisation of the positions from that quantum state. Using this property, one can generate realisations of the determinantal point process using the following procedure: 1. Generate the index M of the highest occupied level, using that 2. Generate the occupation numbers for k < M , from the measure Proba(n k = 1) = λ k .
Set n M = 1 and n p = 0 for p > M . Note that this realisation will contain N = k n k points. Denote also {k n } n=1,...,N the indices of the occupied levels (n k i = 1) and v(x) = (ψ k 1 (x), . . . , ψ k N (x)) T . 3. Pick the first point X 1 from the distribution and introduce e 1 = v(X 1 )/|| v(X 1 )||. 4. Knowing the positions {X 1 , . . . , X n } of the first n points and the set of unit vectors ( e 1 , . . . , e n ) generate the position X n+1 of the next point from the distribution Construct e n+1 = w n+1 /|| w n+1 ||, where ( e j * · v(X n+1 )) e j .