Correlations of occupation numbers in the canonical ensemble and application to BEC in a 1D harmonic trap

# Correlations of occupation numbers in the canonical ensemble and application to BEC in a 1D harmonic trap

Olivier Giraud LPTMS, CNRS, Univ. Paris-Sud, Université Paris-Saclay, F-91405 Orsay, France    Aurélien Grabsch LPTMS, CNRS, Univ. Paris-Sud, Université Paris-Saclay, F-91405 Orsay, France    Christophe Texier LPTMS, CNRS, Univ. Paris-Sud, Université Paris-Saclay, F-91405 Orsay, France
February 7, 2018
###### Abstract

We study statistical properties of non-interacting identical bosons or fermions in the canonical ensemble. We derive several general representations for the -point correlation function of occupation numbers . We demonstrate that it can be expressed as a ratio of two determinants involving the (canonical) mean occupations , …, , which can themselves be conveniently expressed in terms of the -body partition functions (with ). We draw some connection with the theory of symmetric functions, and obtain an expression of the correlation function in terms of Schur functions. Our findings are illustrated by revisiting the problem of Bose-Einstein condensation in a 1D harmonic trap, for which we get analytical results. We get the moments of the occupation numbers and the correlation between ground state and excited state occupancies. In the temperature regime dominated by quantum correlations, the distribution of the ground state occupancy is shown to be a truncated Gumbel law. The Gumbel law, describing extreme value statistics, is obtained when the temperature is much smaller than the Bose-Einstein temperature.

05.30.-d

## I Introduction

The theory of noninteracting identical quantum particles is a fundamental block of the basic education in statistical physics Hua63 ; PatBea11 ; TexRou17book . In the standard approach, calculations are performed in the grand canonical ensemble, as it provides the clearest and most efficient tools to relate single-particle and thermodynamic properties. One can then use the equivalence between statistical physics ensembles in the thermodynamic limit in order to get the thermodynamic observables as a function of the relevant parameters. One should however keep in mind that the correspondence between ensembles in the thermodynamic limit only holds for averages of observables, and not for their fluctuations PatBea11 ; TexRou17book . Recently, remarkable progress in atomic physics of ultracold atoms has raised many questions concerning many-body effects in those systems BloDalZwe08 . Simpler questions related to quantum correlations in noninteracting gases have also been put forward, as experiments deal with extremely diluted gases, for which the noninteracting limit is often a good starting point. Due to cooling techniques by evaporation, a trapped ultracold atomic gas only contains a moderately large number of atoms (few thousands to few millions) ; this has led to re-examine the differences between the various statistical physics ensembles BroDevLem96 ; SchMed96 ; GroHol96 ; WeiWil97 ; WilWei97 ; GajRza97 ; HolKalKir98 ; ChaMekZam99 ; Pra00 , as the microcanonical or canonical ensembles are more relevant in this case.

### i.1 Occupation numbers

Not only thermodynamic properties and global observables are of interest, but, motivated by the remarkable achievement of the “atomic microscope” CheNicOkaGerRamBakWasLomZwi15 ; HalHudKelCotPeaBruKuh15 ; ParHubMazChiSetWooBlaGre15 , also local observables have been studied recently (for a review see Ref. DeaLeDMajSch16 ). A basic ingredient of such studies is the knowledge of the number of particles in each individual eigenstate . The grand canonical mean occupation is given by the Bose-Einstein or Fermi-Dirac distribution

 ¯¯¯¯¯¯nλg=11/(xλφ)∓1,xλ=e−βελ, (1)

where is the fugacity, are the individual energy levels and is the inverse temperature. is the grand canonical average. In this formula, the upper () and lower () signs stand for bosons and fermions, respectively. For a fixed number of particles, the canonical mean occupation numbers can be expressed in terms of the -body canonical partition functions with . The canonical partition function for bosons or fermions can be obtained by recursion from the formula Lan61 ; For71 ; BorFra93 ; SchSch02 ; MulFer03

 ZN(β)=1NN∑k=1(±1)k−1Z1(kβ)ZN−k(β), (2)

where is the canonical partition function for one particle, and the upper () and lower () signs stand for bosons and fermions, respectively. The canonical mean occupation number is then given by

 ¯¯¯¯¯¯nλ=N∑k=1(±1)k−1ZN−kZNxkλ (3)

with WeiWil97 ; BorHarMulHil99 ; ChaMekZam99 . The canonical average is simply denoted . This expression is expected to coincide with (1) in the thermodynamic limit, provided that the fugacity in Eq. (1) is chosen in such a way that the condition is fulfilled. The power of the grand canonical ensemble lies in the independence of individual energy level properties (which leads to many useful additivity properties for the thermodynamic observables), i.e. the absence of correlations between occupation numbers

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯nλ1⋯nλpg=¯¯¯¯¯¯¯¯nλ1g×⋯×¯¯¯¯¯¯¯¯nλpg. (4)

However, in the canonical ensemble, the constraint on the total number of particles implies non-trivial correlations between occupation numbers, . In the present paper, we focus on the canonical ensemble and study these correlations.

### i.2 Main results

Our main results are complementary expressions for the correlation function. In order to lighten notations, we consider the correlations of the first levels ; this does not imply any restriction on the generality of the discussion, as levels all play an equivalent role. The first expression we obtained is in terms of the canonical partition functions and of the , and reads

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯n1⋯np =(±1)pZN (5) ×N∑k1,…,kp=1k1+⋯+kp⩽N(±x1)k1⋯(±xp)kpZN−k1−⋯−kp.

This generalizes (3). The second expression is a representation in terms of two determinants,

 (6)

where the denominator is a Vandermonde determinant. A remarkable observation is that the -point correlation function can be expressed only in terms of the mean values . It can also be expressed in terms of -point functions with , see Eq. (36) below. Equation (6) turns out to be valid also in the grand canonical ensemble. A particular instance of this general relation for ,

 ¯¯¯¯¯¯¯¯¯¯¯n1n2=∓eβε1¯¯¯¯¯n1−eβε2¯¯¯¯¯n2eβε1−eβε2, (7)

was recently used for the study of fluctuations of certain observables for trapped noninteracting one-dimensional fermions GraMajSchTex17 .

In the case of bosons, we could express the general correlation function for arbitrary integer powers as

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯nr11⋯nrpp=1ZN (8) ×N∑k1,…,kp=1 k1+⋯+kp⩽Nak1(r1)xk11⋯akp(rp)xkppZN−k1−⋯−kp,

where

 ak(r)=kr−(k−1)r. (9)

Equation (8) generalizes (5).

We illustrate these formulae by considering the study of Bose-Einstein condensation in a one-dimensional harmonic trap. We obtain several simple analytical results for the occupation numbers of single-particle levels. In particular, we show that, in the temperature regime , where is the trap frequency, the individual ground state occupancy is distributed according to a truncated Gumbel law :

 Proba{n0=N−Tln(T/ω)ω+Tωζ}∝θH(ζ−lnz)eζ−eζ, (10)

where and is the Heaviside step function.

### i.3 Outline

The outline of the article is as follows : in Section II, we recall the connection between our problem and the theory of symmetric functions, which will allow us to introduce some useful tools. Our main results expressing correlations between occupation numbers, Eqs. (5) and (6), are derived in Section III. Finally, in Section IV, we illustrate our results on the problem of Bose-Einstein condensation in a one-dimensional harmonic trap.

## Ii Symmetric functions

The connection between the problem of identical particles and the theory of symmetric functions has been discussed in Refs. For71 ; Bal01 . In Ref. SchSch02 , Schmidt and Schnack have pointed out that the relation (2) for fermions, which is attributed to Landsberg Lan61 in many articles, is nothing else but the well-known Newton identity. In this section, we introduce some useful notation. As a simple illustration, we will recover the relations (2). For a recent reference on the mathematical theory of symmetric functions, see the monograph Mac95 .

### ii.1 Families of symmetric polynomials

A function in variables is said to be symmetric if for any permutation of the indices. We now introduce three useful families of symmetric polynomials.

The elementary symmetric polynomials are defined as

 eN(x1,…,xM)=∑λ1<λ2<⋯<λNxλ1xλ2⋯xλN, (11)

with and for . For example, and . Their generating function is given by

 ΞF(φ)=M∑N=0eNφN=M∏λ=1(1+φxλ). (12)

The complete homogeneous symmetric polynomials are defined as

 hN(x1,…,xM)=∑λ1⩽λ2⩽⋯⩽λNxλ1xλ2⋯xλN, (13)

with . For example, . Their generating function is given by

 ΞB(φ)=∞∑N=0hNφN=M∏λ=1(1−φxλ)−1. (14)

Note that, contrary to the sum in the definition of , the sum here extends to infinity, as for .

The power sum polynomials are defined as

 (15)

Their generating function is given by

 P(φ)=∞∑k=1pkφk−1=M∑λ=1xλ1−φxλ (16)

(for convenience the definition is slightly different from those of the previous generating functions as we did not introduce a ).

##### Correspondence with the problem of identical particles.—

Let us consider particles in energy levels (possibly infinite). Setting as in Eq. (1), one readily sees that the canonical partition function for bosons coincides with , while the canonical partition function for fermions coincides with . Moreover, one obviously has , and , so that the single-particle canonical partition function at inverse temperature , which is or , coincides with . One can thus establish the following dictionary between the mathematician’s and the physicist’s notations :

 xλ ⟶e−βελ eN ⟶ZFN(β) hN ⟶ZBN(β) pk ⟶ZB1(kβ)=ZF1(kβ).

The generating functions and coincide with the grand canonical partition functions for bosons and fermions, respectively.

### ii.2 Newton identity

There exist various identities relating the generating functions , and . Expanding these identities in powers of provides some relations between the three families of symmetric polynomials defined above. For instance, the duality relation

 ΞB(φ)ΞF(−φ)=1, (17)

readily obtained from (12) and (14), allows to express the ’s in terms of the ’s, or conversely.

Another identity that can be easily obtained from the expressions of the previous subsection is

 P(φ)=ddφlnΞB(φ)=−ddφlnΞF(−φ), (18)

that is,

 −ddφΞF(−φ)=P(φ)ΞF(−φ) (19)

and

 ddφΞB(φ)=P(φ)ΞB(φ). (20)

Expanding explicitly (19) in powers of yields

 M∑N=1N(−φ)N−1eN=∞∑k=1pkφk−1M∑j=0(−φ)jej. (21)

Identification of terms in in the r.h.s. gives the relation

 eN=1NN∑k=1(−1)k−1pkeN−k. (22)

This is precisely Eq. (2) for fermions, according to the above dictionary. The relation (22) was derived by Isaac Newton in his book, published in 1666 (see Ref. Whi67 , p. 519). It is known as the Newton identity Mac95 . Several instances of the relation, for and , were obtained earlier by Albert Girard in 1629 (see Ref. Gir1629 page F2, i.e. page  of the manuscript, where the elementary polynomial is called “-th meslé”).

Similarly, expanding (20) in powers of gives

 ∞∑N=1NφN−1hN=∞∑k=1pkφk−1∞∑j=0φjhj, (23)

 hN=1NN∑k=1pkhN−k. (24)

This corresponds to Eq. (2) for bosons.

The theory of symmetric functions also provides a determinantal representation of elementary symmetric polynomials in terms of power sum polynomials (p. 28 of Mac95 ), as

 eN=1N!∣∣ ∣ ∣ ∣ ∣ ∣∣p110⋯0p2p120⋮⋮⋱⋱⋮pN−1pN−2N−1pNpN−1pN−2⋯p1∣∣ ∣ ∣ ∣ ∣ ∣∣. (25)

This provides an explicit expression of the -body partition function in terms of the one particle partition function . The homogeneous polynomials can also be expressed by the r.h.s. of Eq. (25), but with the determinant replaced by a permanent For71 . Alternatively, they can be expressed in terms of a determinant Mac95 , as

 hN=1N!∣∣ ∣ ∣ ∣ ∣ ∣∣p1−10⋯0p2p1−20⋮⋮⋱⋱⋮pN−1pN−2−N+1pNpN−1pN−2⋯p1∣∣ ∣ ∣ ∣ ∣ ∣∣, (26)

which provides an expression of in terms of .

## Iii Correlation functions

The tools developed in the previous section allow us to easily derive Eqs. (5) and (6), as we now show.

### iii.1 Canonical and grand canonical ensembles

We consider identical particles in (possibly infinite) energy levels. The occupation number of the individual eigenstate of energy is denoted by . Thus we have for fermions and for bosons. A basis of Fock space is given by the many-body states , which are fully specified by the knowledge of all occupation numbers. We can express the canonical (Gibbs) distribution at inverse temperature as

 PcN({nλ})=∏λxnλλZNδN,∑λnλ,xλ=e−βελ, (27)

which gives the probability of occupying the quantum state . Here is the -body canonical partition function, and the Kronecker symbol constrains the number of particles to be . On the other hand, the grand canonical distribution is controlled by the fugacity , and reads

 Pg({nλ};φ)=∏λ(xλφ)nλΞ(φ) (28)

with the grand canonical partition function given by Eqs. (12) or (14). For bosons, the convergence of the series is ensured by the condition , where is the individual ground state ().

Using these distributions, one can relate the grand canonical average and the canonical average for particles (the superscript will only be introduced if needed). Indeed, if is any function of the occupation numbers, then, from (27)–(28) one has

 Ξ(φ)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯A({nλ})g=∞∑N=0ZN¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯A({nλ})(N)φN. (29)

### iii.2 p-point correlation functions

#### iii.2.1 Proof of Eqs. (5) and (6)

We now apply (29) to . In the grand canonical ensemble, the occupation numbers are independent (see Eq. (4)), and they are given by Eq. (1). We thus get from (29)

 ∞∑N=pZN¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯n1⋯np(N)φN = Ξ(φ)[1/(x1φ)∓1]⋯[1/(xpφ)∓1]. (30)

The expansion of (30) and the identification of each power of directly gives Eq. (5). Consider for example . We have explicitly

 ∞∑N=1ZN¯¯¯¯¯n1(N)φN=x1φ∞∑k=0(±x1φ)k∞∑m=0Zmφm. (31)

Identification of the terms on both sides gives , which is Eq. (3).

We now introduce the determinant

 ∞∑N=1ZN∣∣ ∣ ∣ ∣ ∣ ∣∣¯¯¯¯¯n1(N)x1x21⋯xp−11¯¯¯¯¯n2(N)x2x22⋯xp−12⋮⋮⋮⋱⋮¯¯¯¯¯np(N)xpx2p⋯xp−1p∣∣ ∣ ∣ ∣ ∣ ∣∣φN=∣∣ ∣ ∣ ∣ ∣ ∣∣∑NZN¯¯¯¯¯n1(N)φNx1x21⋯xp−11∑NZN¯¯¯¯¯n2(N)φNx2x22⋯xp−12⋮⋮⋮⋱⋮∑NZN¯¯¯¯¯np(N)φNxpx2p⋯xp−1p∣∣ ∣ ∣ ∣ ∣ ∣∣. (32)

Inserting (31) in the right-hand side of this expression, we get

 ∣∣ ∣ ∣ ∣ ∣ ∣ ∣∣x1φΞ(φ)1∓x1φx1x21⋯xp−11x2φΞ(φ)1∓x2φx2x22⋯xp−12⋮⋮⋮⋱⋮xpφΞ(φ)1∓xpφxpx2p⋯xp−1p∣∣ ∣ ∣ ∣ ∣ ∣ ∣∣=x1⋯xpφΞ(φ)(1∓x1φ)⋯(1∓xpφ)∣∣ ∣ ∣ ∣ ∣ ∣∣1(1∓x1φ)x1(1∓x1φ)⋯xp−11(1∓x1φ)1(1∓x2φ)x2(1∓x2φ)⋯xp−12(1∓x2φ)⋮⋮⋮⋱⋮1(1∓xpφ)xp(1∓xpφ)⋯xp−1p(1∓xpφ)∣∣ ∣ ∣ ∣ ∣ ∣∣. (33)

From linear combinations of columns in the determinant we then obtain

 x1⋯xpφΞ(φ)(1∓x1φ)⋯(1∓xpφ)(∓φ)p−1∣∣ ∣ ∣ ∣ ∣ ∣∣1x1x21⋯xp−111x2x22⋯xp−12⋮⋮⋮⋱⋮1xpx2p⋯xp−1p∣∣ ∣ ∣ ∣ ∣ ∣∣=(∓1)p−1∣∣ ∣ ∣ ∣ ∣ ∣∣1x1x21⋯xp−111x2x22⋯xp−12⋮⋮⋮⋱⋮1xpx2p⋯xp−1p∣∣ ∣ ∣ ∣ ∣ ∣∣∞∑N=pZN¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯n1⋯np(N)φN, (34)

where we have used (30) for the last equality. Identifying terms in in the left-hand side of Eq. (32) and in the right-hand side of Eq. (34) demonstrates Eq. (6). In order to prove (6) in the grand canonical case, it suffices to replace by 1 in the left-hand sides of (33) and (34) and use the absence of correlation between occupation numbers Eq. (4).

#### iii.2.2 Generalization

The same technique allows us to generalize Eq. (6) straightforwardly : the -point correlation function can be expressed in terms of -point correlation functions for any . For instance for we have

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯n1n2⋯np=∓∣∣∣¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯n1n3…npx1¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯n2n3…npx2∣∣∣∣∣∣[]cc1x11x2∣∣∣, (35)

and other such relations obtained by picking different indices among the pairs of indices. More generally, as a function of -point correlation functions we have

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯n1n2⋯np=(∓)q−1∣∣ ∣ ∣ ∣ ∣ ∣∣¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯n1nq+1…npx1⋯xq−11¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯n2nq+1…npx2⋯xq−12⋮⋮⋱⋮¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯nqnq+1…npxq⋯xq−1q∣∣ ∣ ∣ ∣ ∣ ∣∣∣∣ ∣ ∣ ∣ ∣ ∣∣1x1⋯xq−111x2⋯xq−12⋮⋮⋱⋮1xq⋯xq−1q∣∣ ∣ ∣ ∣ ∣ ∣∣ (36)

and other such relations obtained by picking any indices among the possibilities. The steps are exactly the same as in (32)–(34), replacing the correlators in the determinant by their expression (30). For we obviously recover (6). Again, substituting by 1 in the proof shows that Eq. (36) is valid also for the grand canonical ensemble.

#### iii.2.3 Relation with Schur polynomials

The ratio of determinants in Eq. (6) allows us to identify another interesting connection with the theory of symmetric functions, namely with a fourth family of symmetric functions, the so-called Schur polynomials Mac95 . The definition of these polynomials is given in Appendix A; they are expressed as a ratio of two determinants.

The Schur polynomials naturally appear if one replaces the in the numerator of Eq. (6) by their expression (3). By linear expansion of the determinant with respect to its first column we directly obtain

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯n1⋯np=x1⋯xpZNN∑k=p(±1)p+kZN−ksϖk(x1,⋯,xp), (37)

where are the Schur polynomials for the partition , cf. Appendix A. We recall that is identified with either the elementary symmetric polynomial , or the homogeneous symmetric polynomial , where is the dimension of the one-particle Hilbert space. Introducing the Schur functions has allowed us to reduce the -fold sum in Eq. (5) to a single sum in (37).

### iii.3 Correlation functions with higher powers (bosons) : proof of Eq. (8)

In the case of bosons, we can also consider higher moments of the occupation numbers (for fermions we have of course ). This question has attracted a lot of attention for the characterisation of the number of condensed bosons in a BEC Pol96 ; GroHol96 ; GajRza97 ; WilWei97 ; HolKalKir98 . In the grand canonical ensemble, the integer moments of each occupation number can be obtained simply from the individual grand partition function for individual eigenstate as

 ¯¯¯¯¯¯nrλg=1ξλ(φddφ)rξλ=∞∑k=1ak(r)(xλφ)k, (38)

with and . Applying (29) to , and making use of the independence of the occupation numbers in the grand canonical ensemble as in (4), we readily obtain (8). For instance for the second moment we have

 ¯¯¯¯¯¯n2λ=N∑k=1(2k−1)ZN−kZNxkλ. (39)

This representation will be of practical use in the following section.

## Iv Condensation of bosons in a 1D harmonic trap

As a simple illustration of our results, we consider bosons in a 1D harmonic trap with frequency . The problem has been studied within the grand canonical BroDevLem96 ; KetDru96 ; Mul97 ; PetSmi02 ; PetGanShl04 , the canonical BroDevLem96 ; WilWei97 ; WeiWil97 ; BorHarMulHil99 and the microcanonical Tem49 ; GroHol96 ; GajRza97 ; WeiWil97 ensembles 111 Several simple thermodynamic properties (energy, heat capacity) are discussed in Ref. TexRou17book (see also Ref. GraMajSchTex17 ). . In particular some limiting behaviors of the ground state occupancy for , recalled below, were obtained in several of these references. The probability distribution of the occupation number of the -th level was obtained in Ref. WilWei97  :

 Pk,N(n)=¯¯¯¯¯¯¯¯¯¯δnk,n(N)=xkn[1−xk+xN−n+k]ZN−nZN, (40)

where . It is however not straightforward to extract simple information, like moments or cumulants, or to analyze the large- asymptotics of this distribution. Here we show that the canonical formulae obtained in the previous sections lead to simple analytical results appropriate to discuss the large limit.

### iv.1 Thermodynamic properties

Up to a shift in energy, the one-body spectrum is for (we set ). The -body partition function TodKubSai92 ; TexRou17book

 ZN=N∏n=1(1−e−nβω)−1 (41)

corresponds to independent bosonic modes with frequencies for .

The problem involves three characteristic temperature scales. (i) The lowest scale, , separates the regime where the spectrum should be considered discrete () from the one where it can be described as a continuous spectrum (). (ii) The scale separates the quantum regime , where the upper modes are frozen in their ground state (see Eq. (41)), from the classical regime where all the modes can be described as classical oscillators, in which case we recover the Maxwell-Boltzmann partition function corresponding to neglecting the effect of quantum correlations. (iii) The third temperature scale is the Bose-Einstein temperature

 TB=NωlnN, (42)

below which a macroscopic fraction of bosons accumulates in the individual ground state. It can be obtained from the analysis of the canonical chemical potential or the fugacity . Introducing (incorrectly) this expression in the grand canonical expression of the ground state occupancy, , Eq. (1), shows that for i.e. .

In the following, we will not describe the effect of the discrete nature of the spectrum, and will always consider the limit (the condition will be implicit in the rest of the paper). In particular, we can treat the sum over the spectrum in as an integral, so that (41) yields . The latter expression can be reformulated in terms of the polylogarithm function (see §25 of DLMF ). The free energy is then given by

 FN=−1βlnZN≃Li2(e−Nβω)−Li2(1)β2ω. (43)

### iv.2 Mean occupation numbers

#### iv.2.1 Classical regime: T≫T∗

In the classical regime , i.e. , one can use the asymptotic behavior for . From Eq. (43) we recover the classical (Maxwell-Boltzmann) result , which coincides with the expression given above since . The mean occupation number, given by Eq. (3), is dominated by the first term, so that for the -th level it is given by . Since the canonical fugacity behaves as , we get , which coincides with the well-known grand canonical behavior.

#### iv.2.2 Quantum regime: T≪T∗

We now turn to the more interesting regime where (and of course ). The sum in (3) can be replaced by an integral. Using (43) for the expression of , the mean occupation number can be reexpressed as

 ¯¯¯¯¯¯nk≃N∫10dy (44) ×exp{−Nβεky+Li2(e−Nβω)−Li2(e−Nβω(1−y))βω},

where . While the exact form (3) is only tractable for small in pratice, the integral representation (44) has the advantage that it allows to study the occupation without restriction on . One must however keep in mind that (44) only holds in the intermediate regime .

The integral expression (44) can be further simplified using the behavior of the polylogarithm function in the vicinity of 0, for . Indeed, in the regime where , one can replace by . Moreover, when one can also replace by , since in the vicinity of , where this approximation breaks down, the integrand becomes proportional to . For the same reason one can extend the integral to infinity. Equation (44) thus reduces to

 ¯¯¯¯¯¯nk≃N∫∞0dyexp{−Nβεky+e−Nβω−e−Nβω(1−y)βω}. (45)

Introducing the parameter

 z=e−Nβωβω=TωN−TB/T (46)

and making the change of variables , Eq. (45) yields a representation in terms of the incomplete Gamma function,

 ¯¯¯¯¯¯nk≃1βωzkezΓ(−k,z). (47)

In the inset of Fig. 1 we compare (47) with the exact expression (3) : the difference is below on the interval for relatively small and we can see that the temperature range over which the difference remains small grows as increases.

#### iv.2.3 Ground state

Setting in Eq. (47), the incomplete Gamma function reduces to the exponential integral . The regime corresponds to the limiting behavior for , where is the Euler-Mascheroni constant. Hence

 ¯¯¯¯¯n0N≃1−Tln(eγT/ω)Nωfor T≪TB. (48)

This behavior was already obtained in Ref. HolKalKir98 by a different approach (see also Appendix B, where we recall how the limiting behavior can be obtained within a grand canonical treatmeant PetGanShl04 ). In Fig. 1 we compare the approximate expression (48) with the exact sum (3). For large enough the behavior (48) is indistinguishable from the exact result up to .

#### iv.2.4 Excited states

For the excited states, because of the factor , the integral (44) is dominated by the neighbourhood of the lower boundary. In this case we can use for . This is equivalent to replacing (47) by the approximation , which corresponds to the interpolation between the two limiting behaviors for and for (the agreement with becomes excellent at large ).

As a result we get the approximate form

 ¯¯¯¯¯¯nk≃Tεk+Te−Nω/T (49)

for . For we get the linear behavior . However, unlike (48), Eq. (49) also describes the crossover from this linear behavior to the decaying behavior above , as illustrated in Fig. 1. Equation (49) shows that the mean occupation reaches its maximum for , with

 ¯¯¯¯¯¯nk∣∣max≃N/kln(N/k). (50)

The presence of the logarithm in the denominator shows that only the ground state has a macroscopic occupation below , as expected when BEC occurs.

For the highest excited states, such that , the continuous approximation of the sum (3), leading to the integral (44), fails, and the occupation decays exponentially as . This is the expected classical behavior (for ), as weakly occupied levels can be considered in the classical regime.

### iv.3 Variance of the ground state occupation

We now study the fluctuations around the mean occupation number. We restrict ourselves to the ground state, which has attracted some attention in higher dimension Pol96 ; HolKalKir98 . The exact expression for is given by (39). In the most interesting regime, , we get an integral representation similar to (44) : the coefficient in the sum (39) translates into a factor in the integral (44), so that

 ¯¯¯¯¯n20≃2N2∫10dyye[Li2(e−Nβω)−Li2(e−Nβω(1−y))]/(βω). (51)

Performing the same approximations as above and using the asymptotics as , we obtain

 Var(n0