Theory of Non-Interacting Fermions and Bosons in the Canonical Ensemble

We present a self-contained theory for the exact calculation of particle number counting statistics of non-interacting indistinguishable particles in the canonical ensemble. This general framework introduces the concept of auxiliary partition functions, and represents a unification of previous distinct approaches with many known results appearing as direct consequences of the developed mathematical structure. In addition, we introduce a general decomposition of the correlations between occupation numbers in terms of the occupation numbers of individual energy levels, that is valid for both non-degenerate and degenerate spectra. To demonstrate the applicability of the theory in the presence of degeneracy, we compute energy level correlations up to fourth order in a bosonic ring in the presence of a magnetic field.


I. INTRODUCTION
In quantum statistical physics, the analysis of a fixed number N of indistinguishable particles is difficult, even in the non-interacting limit. In such a canonical ensemble, the constraint on the total number of particles gives rise to correlations between the occupation probability and the corresponding occupation numbers of the available energy levels. As a result, the canonical treatment of non-interacting fermions and bosons is in general avoided, and it is often completely absent in introductory texts addressing quantum statistical physics [1][2][3]. The standard approach is to relax the fixed N constraint and instead describe the system using the grand canonical ensemble. While in most cases this approximation can be trusted in the large−N limit, it can explicitly fail in the low-temperature regime [4,5]. The situation is even worse if the low-temperature system under study is mesoscopic, containing a relatively small number of particles.
In this paper, we present a unified framework for the calculation of physical observables in the canonical ensemble for N fermions or bosons that are applicable to general non-interacting Hamiltonians that may contain degenerate energy levels. We present an analysis of the mathematical structure of bosonic and fermionic canonical partition functions in the non-interacting limit that leads to a set of recursion relations for exactly calculating the energy level occupation probabilities and average occupation numbers. Using argument from linear algebra, we show how higher-order correlations between the occupation numbers can be factorized, allowing them to be obtained from the knowledge of the occupation numbers of the corresponding energy levels, a canonical generalization of Wick's theorem [82]. The key observation yielding simplification of calculations in the canonical ensemble is that occupation probabilities and occupation numbers can be expressed via auxiliary partition functions (APFs) -canonical fermionic or bosonic partition functions that correspond to a set of energy levels that are obtained from the full spectrum of the targeted system by making a subset of the energy levels degenerate (increasing its degeneracy if its already degenerate) or alternatively by excluding it from the spectrum. Results obtained via auxiliary partition functions are validated by demonstrating that a number of previously reported formulas can be naturally recovered in a straightforward fashion within this framework.
In a non-interacting system, observables such as the average energy and magnetization can be obtained solely from the knowledge of average occupation numbers, however, the calculation of the corresponding statistical fluctuations in such quantities, i.e., specific heat and mag-arXiv:2007.15661v2 [cond-mat.stat-mech] 8 Aug 2020 netic susceptibility, requires knowledge of the fluctuations in occupation numbers and (in the canonical ensemble) correlations between them. Therefore, the factorization of correlations between occupation numbers in terms of the average occupation numbers of individual energy levels provides a simplified approach to calculate quantities such as specific heat and magnetic susceptibility of a given system. Here, we highlight a method to compute such correlations by considering a degenerate system of a finite periodic chain of non-interacting bosons that is influenced by an external uniform magnetic field.
The auxiliary partition function approach presented here provides a set of tools that can be used to analyze experimental data in low-density atomic gases where the number of particles is fixed. The application of the physically relevant canonical ensemble can eliminate errors introduced via the grand canonical approximation (especially in the inferred temperature) and lead to a more accurate interpretation of experimental results, including an improved diagnosis of the role of weak interaction effects. It is hoped that the relative simplicity of the mathematical approach presented in this paper may encourage the inclusion of the interesting topic of the canonical treatment of Fermi and Bose gases in collegelevel textbooks.
In Sec.II we present the general formalism of the theory, where we introduce APFs and write occupation probability distributions, occupation numbers, and their correlations in terms of the APFs. In the same section, we also show how some of the previously known results could be recovered in a straightforward fashion. In Sec.III, considering fermionic and bosonic systems, we derive the decomposition of level higher-order occupation numbers correlations into individual energy levels occupation numbers for non-degenerate and degenerate energy spectra, alike. To illustrate the applicability of the theory, In Sec.IV, we consider a bosonic ring in the presence of a magnetic field. We conclude in Sec.V.

II. NON-INTERACTING INDISTINGUISHABLE PARTICLES IN THE CANONICAL ENSEMBLE
As the thermodynamic properties of a system of noninteracting particles are governed by the single-particle spectrum and the underlying particle statistics, we begin by considering the general one-particle spectrum i , with i ∈ S = {1, 2, , . . . , M }. For an unbounded spectrum M → ∞. In the canonical ensemble defined by fixing the total particle number N , the canonical partition function Z N ≡ Z N (S) for N indistinguishable particles is defined by: where are the Boltzmann factors at inverse temperature β = 1/k B T and the components {n i } of the vector n| N = (n 1 , . . . , n M )| N are the occupation numbers for the corresponding energy levels satisfying i∈S n i = N . The summation in Eq. (1) runs over all the possible occupation vectors n| N which, in addition to conserving the total number of particles N , obeys occupation limits for each of the energy levels i : Consider the spectrum defined by S as the union of disjoint subsets (subspectra) S (1) and S (2) = S \ S (1) , i.e., S = S (1) ∪ S (2) . Under this decomposition the Boltzmann factors in S can be factorized as: and N − k particles in S (1) and S (2) , respectively. Summing X(n| N ) over all possible n (1) | k and n (2) | N −k gives: where we have introduced the APFs: For Eqs. (4) and (5)  are the maximum numbers of particles allowed in S (1) and S (2) , respectively. Additionally, any vector n| N , under these constraints can be decomposed into two allowed vectors n (1) | k and n (2) | N −k and vice versa. Thus the sum in Eq. (1) can be similarly decomposed as: max ) yielding the full partition function Employing the convention that Z N (S) = 0 whenever N is negative, or when it exceeds the maximum number of particles set by S, the limits in the above summation can be simplified to N k=0 . The above notation can be made more explicit by specifying the subset of levels that are not included in the partition function and thus for S (1) = {j 1 , j 2 , . . . , j } containing levels, Eq. (6) is equivalent to To obtain a physical interpretation of Eq. (8), recall that in the canonical ensemble, the likelihood of the Nparticle system being in a microstate defined by the occupation vector n| N is given by the ratio X(n| N )/Z N . Accordingly, for the subset of energy levels with indices {j 1 , j 2 , . . . , j }, the joint probability distribution of the corresponding occupation numbers P nj 1 ,nj 2 ,...,nj can be obtained by performing the summation n (2) It follows that the probability P k ({j 1 , j 2 , . . . , j }) of finding k particles in {j 1 , j 2 , . . . , j } and, of course, N − k particles in S \ {j 1 , j 2 , . . . , j }, can be obtained by applying the summation n (1) | k : where the normalization of P k ({j 1 , j 2 , . . . , j }) is guaranteed by Eq. (8).
So far the analysis of the partition functions has been completely general and we have not specified what type of particles are being described. Now, let's be more specific and consider a system that solely consists of either fermions or bosons.

A. An Inverted Analogy Between Fermionic and Bosonic Statistics
Having developed an intuition for the general structure of the canonical partition function under a bipartition into sub-spectra, we now observe how this can provide insights into the relationship between fermionic (n max i = 1) and bosonic (n max i → ∞) statistics of N non-interacting particles. To distinguish the two cases we introduce a new subscript on the partition function (F for fermions and B for bosons).
If the set S (1) represent a single energy level with an index j 1 = j then using our convention we have: for fermions and for bosons. Substituting into Eq. (8) we immediately find: The relations in Eq. (13) and (14) formally describe the procedure for generating the canonical partition function of the N particle system after introducing an energy level j to the preexisting spectrum S \ {j}.
Examining the structure of Eq. (13) suggests a simple matrix form: , . . . ) and the matrix A n,m = δ n,m +e −β j δ n,m+1 is bidiagonal and can be inverted such that: Comparing this expression with Eq. (14), we observe an identical structure apart from exchanging the factor e −β j with −e −β j . Thus we can obtain the inversion of Eq. (14) by replacing e −β j with −e −β j in Eq. (13), i.e., The relations Eq. (15) and Eq. (16) exemplify the elimination of an energy level as they represent the inverse of Eq. (13) and Eq. (14) respectively. If we absorb the negative signs in Eqs. (15) by shifting the energy j by ±iπ/β, we can write Z \{j} F,N = N k=0 e −β j k Z F,N −k , where j = j ± iπ/β. As the general bipartition into sub-spectra introduced in Eq. (6) holds for energy levels with mixed statistics, and doesn't require real entries for the i , we can build the bosonic partition functions Z B,N ({j 1 , j 2 , . . . , j }) using the shifted energies { j1 , j2 , . . . , j } and then combine it with the Z F,N to generate a mixed-form: The effect of shifting the single particle spectrum by a constant ω on the canonical partition function Z N is captured by a rescaling factor e −βωN and the resulting N particle partition function of the shifted spectrum is e −βωN Z N . Using Z B,k ({j 1 , j 2 , . . . , j }) = e ±iπk Z B,k ({j 1 , j 2 , . . . , j }) then yields The last two equations can be seen as a generalization of Eqs. (15) and (16). However, they also recover the symmetry between fermionic and bosonic statistics. To this end, we show how to obtain the partition function of a given spectrum using the APFs of two complementary subsets of the spectrum through Eq. (8). Also, Eqs. (18) and (19) show that this relation can be inverted, i.e., we can calculate the APFs of a subset of energy level using the APFs of its complement with the opposite statistics and the partition functions of the full spectrum.

B. Energy Level Occupations and Correlations
For fermions, the Pauli exclusion principle restricts the number of particles occupying an energy level j to n j = 0 or 1. Equivalently, n p j = n j for any p > 0, simplifying the calculation of energy level occupation numbers and the correlations between them, including higher moments. More specifically, the aver- (20) using the probability: defined in Eq. (9). The only term that survives in Eq. (20) has n j1 = n j2 = · · · = n j = 1 giving: Following the same procedure, Eq. (20) can be generalized to describe both correlations and anti-correlations between energy levels: where γ jr = 1, 0. For the latter with γ jr = 0, the occupation numbers in Eq. (20) have been replaced with their complements, 1 − n jr . Thus, we find that fermionic level occupations and correlations can be directly written in terms of APFs without resorting to the usual definition n j1 n j2 . . . n j F,N = 1 When considering a single level ( = 1), the occupation probability of the j th fermionic level immediately follows with the associated probability For bosons, the occupation numbers n j B,N can be calculated from the corresponding occupation probability distribution P nj which are obtained from Eq. (9): and for the -point correlations: However, unlike the fermionic case, such an approach requires performing the unrestricted summation n({j1,j2,...,j }) . An alternative method which avoids this difficulty can be developed by exploiting the inverted analogy between fermionic and bosonic statistics introduced in Sec.II A. In the fermionic case, the occupation number of an energy level j is proportional to the APF Z \{j} F,N −1 (Eq. (24)) which corresponds to the actual spectrum of the system missing the energy level j . This suggests a route forward for bosons via the analogous inversion of doubly including the energy level j instead of removing it, i.e., we construct an APF where this level is twofold-degenerate. We denote the corresponding N -boson APF by Z ∪{j} B,N and distinguish the two levels using the dressed indices j (0) and j (1) such that the resulting combined spectrum has level indices {1, . . . j − 1, Returning to the general definition of the canonical partition function in Eq. (1), we can write Z ∪{j} B,N = n | N X(n | N ), where the occupation vectors n | N have one extra component: n j in n| N is replaced by n j (0) and n j (1) . The modified Boltzmann factors are: and thus their value is dependent only on the total occupancy of the j th level, n j (0) + n j (1) . As a result, X(n| N ) = X(n | N ) for any occupation vectors n| N and n | N with all i = j components equal as well as [n| N ] j = n j = n j (0) + n j (1) . For fixed n| N , the number of vectors n | N that satisfies the previous conditions is equal to n j + 1, or, the number of ways in which n j bosons can occupy two energy levels. The APF can then be written in terms of the original occupation vector n| N by inserting a frequency factor to account for the extra level j: = Z B,N n j + 1 B,N and thus we can write: which can be substituted into Eq. (30) to arrive at which is in the same form as Eq. (24) for fermions.
To generalize this expression to -level correlations with > 1 we examine the numerator of Eq. (31), recalling that we have added an extra copy of energy level j to the partition function for N − 1 particles such that it now appears m j +1 = 2 times in the associated spectrum: Here we obtain the second line using the same trick as in Eq.
is the total number of particles occupying the level j as it appears in the Boltzmann factor X(n| N −1 ) = e −β jñj i =j e −β ini . The complicated looking binomial is the multiset coefficient that counts the number of waysñ j bosons can be distributed amongst the m j + 1 levels with energy j . Finally, the last line is obtained by using the fact that Eq. (32) can be immediately extended to the case where we add copies of not 1 but levels {j 1 , . . . , j }: where the conditions n jr ≥ 1 in the occupancy vector can be neglected as any n jr = 0 terms do not contribute to the sum due to the multiplicative string. Finally, we can write the desired result: An immediate extension of Eq. (34) will turn out to be useful, which introduces an APF with higher-order degeneracy. We consider m jr extra copies of the r th level jr with r ∈ {1, . . . } and find: r=1 n jr − q jr + m jr m jr where the q jr ≤ m jr allow for the added freedom of choice of how many particles are associated with each of the degenerate levels and allow us to write the left-hand side in terms of the desired occupations n j = [n| N ] j . Note: to simplify notation we only include a superscript on the levels in the bosonic APF and only if we are adding more than one extra copy per original level. Eqs. (23) and (35) are the major results of this section, and demonstrate that for both fermions and bosons,level correlations can be written in terms of APFs for N − particles with energy levels removed (added) for fermions (bosons).

C. Recovering Known Results Via the APF Theory
In this section, we illustrate the utility of our auxiliary expressions in simplifying the derivation of known recursion relations that govern the canonical partition functions and occupation numbers for fermions and bosons.
Beginning with fermionic statistics, if we use the definition of n j F,N in Eq. (24) to substitute for Z \{j} F,N and Z \{j} F,N −1 in Eq. (13), we obtain the well-known recursion relation for occupation numbers [57,59] n j F,N can be written explicitly in term of the partition functions by substituting for Z \{j} F,N −1 in Eq. (24) using Eq. (15) as: and using the canonical condition i n i F,N = N we find the original 1993 result of Borrmann and Franke [58]: where C k = j e −β j k . Bosonic statistics can be treated in analogy to the fermionic case. Assigning the roles played by Eqs. (24), (13) and (15) to Eqs. (31), (16) and (14)[83], respectively, we obtain the known bosonic equivalents [57][58][59]: Due to the Pauli exclusion principle for fermionic statistics, the occupation number n j F,N of an energy level j gives us direct access to the occupation probability distribution P F,nj of the level. Despite the absence of any such simplification in the bosonic case, the occupation probability distribution P B,nj of a single energy level can be related to the corresponding partition functions as which is obtained by using Eq. (26) to substitute for Z \{j} B,N −nj in Eq. (16), after replacing N with N − n j [61]. In summary, within the unified framework of APFs, it is straightforward to obtain most of the well-known general relations in the fermionic and bosonic canonical ensemble that were previously derived using a host of different methods. This highlights the utility of this approach as a unifying framework when studying N indistinguishable non-interacting particles.

D. General Expressions for Probabilities and Correlations
De facto, the APFs can also be used to generalize previous results, and in a form that is highly symmetric with respect to particle statistics. To accentuate this, let us introduce the notation: Then, the recursive relations for energy level correlations can be obtained using Eqs. (13) and (16) for fermions and bosons, respectively. The number of initial values of correlations needed is equal to the number of the involved points, e.g., for the two-level correlations we find: (44) Further, for the set of levels S = {j 1 , . . . , j }, Eqs. (21) and (27), define the joint probability distributions of the occupation numbers in terms of the APFs Z \S F,k and Z \S B,k , which can be re-expressed using Eq. (18) and Eq. (19) as 45) where, E tot = r=1 jr n jr and n tot = r=1 n jr . To avoid confusion we note again the convention in use that Z −ζ,k (S ) = 0 whenever k exceeds the maximum number of particles set by S for fixed particle statistics. Now, as the level occupation correlations of fermions and bosons are represented by the APFs Z \S F,k and Z ∪S B,k in Eqs. (22) and (34), we can write where we note that the k-particle bosonic partition function for the levels S appears in both fermionic and bosonic correlations.
Alternatively, the APFs Z and thus such correlations for fermions and bosons can be computed directly from Eq. (46) as

Expectation values of higher moments of degenerate levels
Let us now focus on the case of bosons and revisit Eq. (35) considering a single energy level j : Similarly, using Eq. (6), we can write the APF in terms of the system partition function and the bosonic APF Z B,N −q ({j (1) , . . . , j (m j ) }) of m j degenerate levels. As a result, we obtain where q ≤ m j . Note that if we set both of m j = and q = , then comparing with Eq. (48) we see that for s ∈ {0, . . . , − 1}. This demonstrates that the correlations between degenerate levels can be expressed in terms of moments of the occupation number of any of the degenerate levels. This also helps to simplify the calculation of such moments, for example, we can write which can be simplified using Eq. (50) as In the same fashion, we can write and thus

III. DECOMPOSITION OF LEVEL CORRELATIONS INTO OCCUPATION NUMBERS
In the previous section, we illustrated that the joint probability distributions of the occupation numbers and corresponding level correlations can be represented by auxiliary partition functions. The APF is distinguished from the actual partition function of the N -particle system through either the inclusion or exclusion of a set of levels from/to the complete spectrum under study. The resulting complexity of performing an actual calculation thus depends on the size of the modified set as demonstrated by, e.g. Eq. (46).
It is known that the resulting complexity can be reduced by relating higher-order correlations between nondegenerate levels to the related level-occupation numbers [63,64], representing an approach similar to Wicks theorem which only holds in the grand canonical ensemble. In this section, we directly obtain many known results using the APF method and, more importantly, generalize them to deal with degenerate energy levels.
Before we introduce the general and systematic approach to this problem, let us return to Eqs. (23) and (35) and consider two specific examples of increasing difficulty.
A. Examples

Two-level correlations
Consider the expectation value of two bosonic energy levels j 1 and j 2 where j1 = j2 : n j1 n j2 B,N . Employing Eq. (35) with m j1 = m j2 = 1, q j1 = 0 and q j2 = 1, we find and upon exchanging the values of q j1 and q j2 , we have Next, eliminating the APF from the two equations yields our final result: Similarly, if the levels are fermionic, we use Eq. (23) with (γ j1 , γ j2 ) = (0, 1) and (1, 0) to find: These known results [63,64] are thus obtainable within the APF approach with a few lines of algebra by generating a set of independent equations. We now extend this idea to three energy levels.

Three-level correlations
The previous example for bosons is modified by adding a third level j 3 with m j3 = 1 and j3 that is different than both of j1 and j2 . Setting (q j1 , q j2 , q j3 ) = (1, 0, 1) and (1, 1, 0) in Eq. (35) gives the two equations respectively. Solving for n j1 n j2 n j3 B,N leads to n j1 n j2 n j3 B,N = − e β j 2 n j1 n j2 B,N − e β j 3 n j1 n j3 B,N e β j 2 − e β j 3 .
(60) which can be further broken down into single-level occupation numbers by application of Eq. (56).
A slightly modified approach can be used if two of the energy levels are degenerate, j2 = j3 = j1 . As above, we denote degenerate levels via superscript and we relabel j 2 = j , say (1, 1, 0), with (0, 1, 1), which gives We now turn to the general decomposition of -level correlations into functions of the occupation numbers of the energy levels for both the non-degenerate and degenerate spectra.

B. A Systematic Approach
In the following we show that Eqs. (6), (18) and (19) can be directly employed to systematically relate higherorder correlations to levels occupation numbers. We begin by considering the set of levels S = {j 1 , j 2 , . . . , j } and relate the corresponding -point correlations to the r-point correlations for any nonempty subset of levels Imposing fermionic level statistics on the spectrum S, we can relate the correlations n j1 n j2 . . . n j F,N to n i1 n i2 . . . n ir F,N by relating their corresponding APFs, i.e., Z We note the upper limit in the previous summation is − r, where in general it should be k max = min( − r, N − r) (Eq. (6)). This is because the fermionic APF Z \Sr F,k (S ) cannot describe more than − r particles, as this is the number of levels that are in S \ S r . Further, for N < the correlations n j1 n j2 . . . n j N = 0 for general particle statistics. As a result, we only consider -points correlations with N ≥ . Moreover, we can obtain Z \Sr F,k (S ) by removing the contribution of the levels S r from the APF of S , using Eq. (18), as If we isolate the last term in the summation in Eq. (62), i.e., Z After changing the order of the summations and shifting the indexes k → k − r + 1 and m → m − r, we obtain the unwieldy expression Next, we substitute for Z where Y F,0 (S ) = n j1 n j2 . . . n j F,N are independent of S r . Therefore, for each choice of the subset S r we can write the linear nonhomogeneous equation (66) in the variables Y F,m with coefficients: A 0 = 1, A 0<m<r = 0 and A r≤m≤ −1 (S r ) = (−1) r−1 e −β iν ∈Sr iν Z B,m−r (S r ). The homogeneity of the linear equation is violated by the term b F (S r ) = n i1 n i2 . . . n ir F,N .
With this formulation, we observe that for any of the 2 − 2 choices of S r , we can write a linear equation in the same variables Y F,m , where Y F,0 = n j1 n j2 . . . n j F,N is the -point correlation while the remaining −1 variables are auxiliary, and depend symmetrically on the levels in S . Also, the r-point correlation n i1 n i2 . . . n ir F,N of the levels in S r plays the role of the nonhomogeneous term in the linear equation and the coefficients A m of the equation can be determined by the bosonic APFs of S r . An analogous expression can be obtained for bosonic statistics, with the same coefficients A m where, in this case, the variables are with Y B,0 (S ) = (−1) −1 n j1 n j2 . . . n j B,N and b B (S r ) = (−1) r−1 n i1 n i2 . . . n ir B,N (see Appendix A for a complete derivation).

Non-degenerate levels
Consider the set S specifying a set of distinct energy levels and choose S r=1 such that it contains only one of the levels in S . We can use Eq. (66) or (69) to construct a set of linear equations each corresponding to one level j s ∈ S with energy js . For fermions the equations are where the coefficients A m were obtained from the singlelevel bosonic APF Z B,m ({j s }) = e −mβ js in Eq. (12). Therefore, using the set of the independent linear equations in variables defined by Eq. (71), we can solve for Y F,0 (S ) = n j1 n j2 . . . n j F,N as This result was recently obtained by Giraud, Grabsch and Texier [64], using the properties of the Schur func-tions and it can be simplified using Vandermonde determinants: where Thus the fermionic -level correlation can be simplified as: . (75) This expression was also recently derived using an elegant second quantization scheme [63]. An equivalent procedure can be performed for bosons, again using the set S r=1 to yield linearly independent equations such that n j1 n j2 . . . n j B,N = (−1) −1 Y B,0 (S ) can also be expressed in terms of determinants leading to [64] n j1 n j2 . . . n j B,N = (−1) We now study the most general possible -level correlation function defined by the set S which could include both degenerate and non-degenerate levels. Consider the subset of level indices S mi = {i (0) , . . . , i (mi−1) } ⊂ S which contains m i > 1 degenerate levels (for m i = 1, we reproduce the non-degenerate analysis discussed above).
As a result, the corresponding m i equations, out of the total set of linear equations defined by Eqs. (71) and (76) for fermions and bosons respectively are identical, as they are distinguished from each other only via the energies of the involved levels and their occupation numbers. Moreover, S , could contain multiple subsets of degenerate energy levels, further complicating the problem. Eqs. (75) and (77) can not be applied in this case, as is apparent from their vanishing denominators whenever j k = jr .
To resolve the complication introduced by degenerate subsets, we generalize the procedure in Sec.III A 2 to treat the case where a three-level correlation contained a subset of 2 degenerate levels. More explicitly, we relate the -points correlations not only to the occupation numbers of the degenerate subsets, but to all of the m i distinct r-points correlations between the degenerate levels with 2 ≤ r ≤ m i . This is useful as we have already introduced Eq. (48), which simplifies the calculation of the correlations between degenerate energy levels for fermionic and bosonic statistics. In addition, it will result in the generation of m i new independent equations that could be used to calculate -points correlations.
Accordingly, for any choice of S r = {i (0) , . . . , i (r−1) } ⊂ S mi ⊂ S , the corresponding coefficients in the constructed linear equations are A 0 = 1, A 0<m<r = 0 and where i is the energy of all degenerate levels in S mi and we have substituted for Z B,m−r ({i (0) , . . . , i (r) }) = m−1 r−1 e −(m−r)β i , a bosonic partition function of r degenerate levels.
To illustrate how this works in practice, consider the 4-point correlation of the levels {j 1 (0) , j 1 (1) , j 1 (2) , j 2 }, labeling distinct energies 1 and 2 . The resulting new set of equations can be solved for both fermions and bosons to give:

IV. APPLICATIONS
To illustrate the applicability of our results for degenerate non-interacting systems of particles with fixed number, and highlight the practical usage of Eq. (80) we consider a one-dimensional tight-binding chain of N spinfull bosons hopping over L lattice sites. Code, scripts, and data used to produce all figures for the bosonic chain results can be found online [84]. The bosonic chain is subject to a static external magnetic field B applied along the z-axis. Spin-S bosons are described by the Hamilto- whereâ † α,σ andâ α,σ are creation and annihilation operators for a boson at site α with σ ∈ {−S, . . . , 0, . . . , S} satisfying [a α,σ , a † α ,σ ] = δ α,α δ σ,σ and t measures the hopping amplitude. S z σ,σ = σδ σ,σ are the matrix elements of the diagonal z-projection of the spin-S representation of the spin operatorŜ. Here, h = gµ B B, where g is the corresponding spin-S g-factor and µ B is the Bohr magneton. We employ periodic boundary conditions, such thatâ L+1,σ =â 1,σ , and to avoid having an unbalanced non-degenerate excited state, we fix the parity of L to be odd.
The tight-binding Hamiltonian in Eq. (81) can be diagonalized wheren j,σ counts the number of bosons with energy and j runs over the finite set S = {− L−1 2 , . . . , 0, . . . , L−1 2 } when L is odd such that |S| = L. An examination of the single-particle spectrum shows that each energy level, except the ground state 0,σ = −2t − hσ, is 2-fold degenerate, where −j,σ = j,σ ; a result of the right-left symmetry of the chain. Turning off the magnetic field and fixing S > 0, gives rise to an extra degeneracy factor of (2S + 1) that affects all levels.
For all numerical results presented in this section, we fix L = 1001, N = 1000 and measure the inverse temperature β = 1 kBT in units of 1/t, where T is the absolute temperature and k B is the Boltzmann constant.

A. Spinless bosons (S = 0)
We begin with the study of spinless bosons, where the model is insensitive to the applied magnetic field and we can drop the subscript σ without loss of generality. Using the single-particle spectrum defined in Eq. (83) with σ = 0 and h = 0 in combination with Eq. (27), we calculate the joint probability distribution P B,n0,n1 of the occupation numbers of the ground state and the first excited state, where, we choose the level j = 1 out of the two degenerate levels j = ±1. Note that we do not bother to use the superscript notation to distinguish degenerate level indices (1 ≡ 1 (0) , −1 ≡ 1 (1) ) here as there are no ambiguities due to the sign of the index j.
The calculation proceeds by obtaining the APFs Z \{0,1} B,k for 0 ≤ k ≤ N using the recursion relation Eq. (41), where the factors C k are calculated using the spectrum S \ {0, 1}. The resulting distribution is where Z B,N can be found by enforcing normalization. We expect Eq. (84) to exhibit interesting features at low temperature where the particles are mostly occupying the ground state with some fluctuations amongst the low lying energy levels. To obtain an estimate of low in this context, we choose a value of the inverse temperature β such that the ground state has a macroscopic occupation corresponding to 50% of the particles. We compare the ratio of the Boltzmann factors of having all particles in the ground state with that of having N/2 particles in the first excited state and the rest in the ground state. Setting the ratio of these factors e −β( 0− 1)N/2 to ∼ 0.1, suggests β ∼ 100/t. The results are illustrated in the left panel of Fig. 1, where the relative broadness of the distribution can be attributed to the degeneracy of the first exited level j = ±1.
If we now calculate the three-level joint probability distribution P B,n0,n1 ,n−1 and consider the fixed slice with n −1 = 0, as presented in the right panel of Fig. 1, we see that the distribution becomes significantly sharper, as blocking the level j = −1 makes the resulting nonnormalized conditional distribution more sensitive to the conservation of the total number of particles.
We now turn to the calculation of the two-level connected correlation function for our bosonic system C(n i , n j ) = n i n j B,N − n i B,N n j B,N .
The first step is to obtain the system partition function Z B,k , recursively, using Eq. (41) starting from Z B,0 up to Z B,N . The occupation numbers n j B,N can then be easily calculated using Eq. (40). All that remains is to calculate the two-points correlations using Eq. (56) for the non-degenerate levels. For the correlations between degenerate levels ( n −j n j B,N ), we use Eq. (48), with ζ = +1.
The results for C(n i , n j ) for all levels i and j at inverse temperature β = 1/t are shown as a heat-map in the lower panel of Fig. 2. Here, we chose temperature β = 1/t in order to distribute the correlations amongst higher energy levels. The upper panel of Fig. 2 shows C(n i , n j ) as a function of n i for fixed 0 ≤ j ≤ 500 corresponding to horizontal cuts through the lower panel were the former is also consistent with surrounding data, as expected from Eq. (51), where n −j n j B,N = nj 2 B,N and the symmetry C(n i , n j ) = C(n −i , n j ) due to the degeneracy.

B. Spin-1 bosons (S = 1)
To illustrate the utility of auxiliary partitions functions in studying correlations in a highly degenerate spectrum, we consider the case of spin-1 bosons. In the absence of a magnetic field (h = 0), each level picks up a degeneracy factor of 2S + 1 = 3, such that the ground state is threefold degenerate and all of the excitation levels are six-fold degenerate.
Degeneracy effects are apparent in the two-level connected correlations C(n i,σ , n j,σ ) at β = 1/t, which we calculate for various values of h as shown in Fig.3. The top-left panel of the figure presents C(n i,σ , n j,σ ) for h = 0, where the choice of σ and σ matters only in the presence of a magnetic field. A comparison with the heat-map of Fig. 2, shows an overall broadening and a reduction of one order of magnitude in the maximum of C(n i,σ , n j,σ ) for S = 1, as compared to S = 0 case.
Removing the energy-spin degeneracy by applying a strong magnetic field of h = 5t, results in a splitting of the spectrum into three bands, each with bandwidth 4t and separated from each other via an energy bandgap of t. In this case, we first focus on correlations between the levels in the lower energy band (σ = 1, bottom-left panel of Fig. 3), and see a partial recovery of the spinless bosons case (Fig. 2). Correlations that involve higher energy levels (σ = 0 and σ = −1) are orders of magnitude weaker, at the considered temperatures, as shown in the right panels of Fig .3.
Finally, we turn to correlations between a set of levels that is partially degenerate, an interesting feature of spin-1 bosons. We consider the four-level (disconnected) correlations n i,1 n j,1 n j,0 n j,−1 B,N between the levels i,1 , j,1 , j,0 and j,−1 , where, in the absence of a magnetic field, the last three levels are degenerate for any j. We employ the bosonic version of Eq. (80) with ζ = 1 with results shown in Fig. 4. According to Eq. (35), the results that we obtain, in this

V. DISCUSSION
In summary, we have presented a statistical theory of non-interacting identical quantum particles in the canonical ensemble, providing a unified framework that symmetrically captures both fermionic and bosonic statistics. Table I includes a listing of our most important results for fermions and bosons. We achieve this by: (1) Representing correlations (Eqs. (23) and (35)) and joint probability distributions ( (21) and (27)) via auxiliary partition functions. (2) Deriving general relations between the canonical partition function of a given spectrum and that of the auxiliary partition function describing a spectral subset, as captured by Eqs. (8), (18) and (19).
These key equations can be manipulated to simplify the derivation of the known recursive relations for partition functions in the canonical ensemble and lead immediately to generalizations, and more importantly, provide useful formulas for calculating the correlations between degenerate energy levels and for calculating higher moments of the occupation numbers distribution. Also, Eqs. (23) and (35) can be used to reduce the complexity order of the desired correlations, or, to relate them to the occupation numbers of the involved levels and correlations between entirely degenerate levels (see Eq. (48)). Moreover, the ability to manipulate the way an auxiliary partition function is built out of other ones, allows us to construct a systematic approach towards the decomposition of many-energy level correlations in terms of individual level occupancies. This reflects the additional constraints between energy levels due to fixed N even in the absence of interactions that are not present in a grand canonical description. Thus, we present an approach to working in the canonical ensemble that includes a generalization of Wick's theorem, where we obtain previous results for non-degenerate levels [63,64] and extend them to the case of a degenerate spectrum.
Interestingly, despite the substantial difference between fermionic and bosonic statistics, the resulting formulas show evident similarity. If we compare Eqs. (36), (37) and (38). with Eqs. (39), (40) and (41), respectively, we see that the differences between the fermionic and the bosonic formulas can be captured by simple ±1 factors. In view of the current theory, such similarity is associated with the interplay between fermionic and bosonic auxiliary partition functions. The inverted symmetry between the two distinct statistics is apparent via a comparison of Eqs. (13) and (14) with Eqs. (15) and (16) reflective of the fact that adding a fermionic energy level to the partition function is similar to excluding a bosonic one and vice-versa.
The presented formulas for combining and resolving auxiliary partition functions allows for their construction via different routes which we have utilized to obtain exact expressions for the decomposition of correlations in terms Here nj is the occupation of the j th level, γj r = 0, 1 and 0 ≤ qj r ≤ mj r ∈ Z. All computations rely on the introduction of auxiliary partition functions that describe a modified spectra or subset of levels connected to S through the removal of levels or the addition of degeneracy.
of single-level occupation numbers. These different forms may also have value in overcoming the known numerical instabilities of the recursive formula for the fermionic partition function due to influence alternating signs [63,85]. In addition, the simplicity of the presented theory suggests a possible generalization to cover different energylevels occupation-constraints beyond the fermionic and bosonic ones. We envision the results presented herein could have applications in the computation of entanglement entropy in the presence of super-selection rules, as well as in modelling cold atom experiments. In the context of quantum information, the spectrum of the reduced density matrix corresponding to a mode bipartition of a state of conserved number N of itinerant particles on a lattice can be associated with that of a fictional entanglement Hamiltonian. For non-interacting particles, the entanglement entropy can be obtained via the so-called correlation matrix method [86][87][88][89] which requires the evaluation of the canonical partition function of the resulting non-interacting entanglement Hamiltonian. For trapped ultra-cold atoms at low densities where N is fixed and interactions can be neglected, the analysis of experimental results in the physically correct canonical ensemble provides improved thermometry, especially for the case of fermions.
Finally, the ability to directly study level statistics in the canonical ensemble for bosons and fermions may have pedagogical value in the teaching of statistical mechanics, where the more physical concept of a fixed number of particles is quickly jettisoned and replaced with a grand canonical reservoir for the sake of simplifying derivations.