Quantum local-equilibrium state with fixed multiplicity constraint and Bose-Einstein momentum correlations

The one- and two-boson momentum spectra are derived in the quantum local-equilibrium canonical ensemble of noninteracting bosons with a fixed particle number constraint. We define the canonical ensemble as a subensemble of events associated with the grand-canonical ensemble. Applying simple hydro-inspired parametrization with parameter values that correspond roughly to the values at the system's breakup in $p+p$ collisions at the LHC energies, we compare our findings with the treatment which is based on the grand-canonical ensembles where mean particle numbers coincide with fixed particle numbers in the canonical ensembles. We observe a significantly greater sensitivity of the two-particle momentum correlation functions to fixed multiplicity constraint compared to one-particle momentum spectra. The results of our analysis may be useful for interpretation of multiplicity-dependent measurements of $p+p$ collision events.


I. INTRODUCTION
Inasmuch as mean particle multiplicities in relativistic heavy ion collisions are large, the whole set of collision events at a fixed energy of nuclear collisions is typically divided into subsets with fixed charged-particle multiplicities.Corresponding multiplicity classes are associated with collision centralities and, thereby, with the initial system's geometry which is primarily characterized by the overall shape of the interaction region.This makes it possible to study the multiplicity dependence of various observables measured at the same energy of collisions.In particular, the fixed particle multiplicity technique has been utilized for analysis of the Bose-Einstein momentum correlations of identical particles.These correlations are typically represented in terms of the interferometry radii.They are the result of the Gaussian fit of the correlation function defined as a ratio of the two-particle spectra to the product of the single-particle ones.These radii reflect the space-time structure and dynamical evolution of the systems created in nuclear collisions (for review of the correlation femtoscopy method see e.g.Ref. [1]).One notable feature of these measurements is that the effective system's volume, when extracted from the Gaussian interferometry radii, appears to scale nearly linearly with charged particle multiplicity (see e.g.Ref. ([2])).This observation is in agreement with the hydrodynamical picture of nuclear collisions.
Recently, because of the start of LHC experiments, the fixed particle multiplicity technique has been utilized for analysis of the Bose-Einstein momentum correlations of identical particles in proton-proton collisions at a fixed energy of collisions.It was observed, in particular, that measured in these collisions interferometry correlation radius parameters do not increase with multiplicity at high charged-particle multiplicities [3,4].While an explanation of this effect is still absent, it is suggestive to assume that the saturation effect in the multiplicity dependence of the interferometry correlation radius parameters takes place once the maximal overlap of colliding nucleons is achieved in most central collisions.Indeed, the color glass condensate effective theory predicts that once maximal overlap is achieved higher multiplicities can only be reached by certain color charge fluctuations, which do not increase the initial size of the system [5].Then, one can speculate that an individual system created in a high-multiplicity p + p collision can be regarded as an element of a quantum-statistical ensemble of systems with various numbers of particles produced under the same initial-state geometry.
In a quantum-statistical framework, observables are the expectation values of the corresponding quantum operators with respect to a suitable statistical operator.For example, successful applicability of almost perfect relativistic hydrodynamics for the description of a particle production in relativistic heavy ion collisions (for a recent review see, e.g., Ref. [6]) indicates that actual state of a system created in collisions with the same centrality can be approximated by a local-equilibrium statistical operator ρ leq , Tr[ρ leq ] = 1, which is obtained by maximizing the von Neumann entropy, S = −Tr[ρ ln ρ], with constrained mean values of energy-momentum and conserved charge densities on a given three-dimensional hypersurface (see, e.g., Ref. [7]).It is noteworthy that high-multiplicity proton-proton collisions exhibit collective behavior similar to that observed in relativistic nuclear collisions.It indicates that a hydrodynamic description of matter formed in these collisions might also be possible [6].
Application of fixed high-multiplicity constraint to p+p collision events means then selecting some subensemble of events with the same initial-state geometry to which the considered system belongs.To assign a quantum statistical state to a subensemble of events with fixed multiplicity, one can utilize the projection operator P N , which automatically invokes such a constraint.The aim of this work (see also Ref. [8]) is to clarify how imposed particle number constraint affects the one-particle spectra and two-boson momentum correlations in a quantum-field local-equilibrium state.It is worth noting that for fairly high particle numbers a canonical ground-state Bose-Einstein condensation can occur.Such a condensation could, in principle, lead to noticeable effects in particle momentum spectra and correlations at fixed multiplicities.This issue is, however, beyond the scope of this paper. 1

II. LOCAL-EQUILIBRIUM STATISTICAL OPERATOR
As a starting point, we consider the quasiequilibrium state (see, e.g., Ref. [7]) of a real relativistic scalar field.This state is represented by the statistical operator ρ q (σ) as (we use 1 For such an analysis, the ground-state of the local-equilibrium statistical operator should be specified, and canonical Bose condensation in the corresponding ground state should be taken into account.For simple nonrelativistic quantum-field models, it was done in Ref. [9], in which the relations of the ground-state Bose-Einstein condensation at a fixed particle number constraint to the particle momentum spectra and correlations were discussed. the convention g µν = diag(+1, −1, −1, −1)) where σ is a three-dimensional spacelike hypersurface with a timelike normal vector n µ (x); is the inverse temperature, and u µ is the 4-velocity) on the hypersurface σ, adjusted such as to satisfy the actual mean values of energy and momentum density at this hypersurface; Z q (σ) is the normalization factor making Tr[ρ q (σ)] = 1; and T µν (x) is a scalar-field energy-momentum tensor.For simplicity, we disregard field self-interactions and consider a noninteracting scalar quantum field model.Then, the T µν (x) reads where the Lagrangian density is Here and The quantization prescription means that a † (p) and a(p) are creation and annihilation operators, respectively, which satisfy the following canonical commutation relations: Before proceeding further, let us digress for a moment and consider the simple case of the covariant global-equilibrium state, where the β µ does not depend on spacetime coordinates across the infinite three-dimensional hypersurface.Then the statistical operator reads where P µ = t d 3 rT µ0 (x) is 4-momentum of the field defined at t = const hypersurface.
Note that a(p, 0) = a(p).Expression (13) implies that a(p, α) satisfies equation Taking into account that this yields then The solution of this equation is Our next step is to combine the cyclic invariance of the trace, Tr[ρ eq (σ)a † (p 1 )a(p 2 )], and Eqs. ( 13) and (17).Using the cyclic invariance of the trace and Eq. ( 13), we obtain Taking into account Eqs. ( 7) and ( 17), the r.h.s. of the above equation can be rewritten as Substituting this into Eq.( 18) we finally have Utilization of the Fourier transformation of Eq. ( 20) with respect to ∆p = p 2 − p 1 immediately results in the ideal gas Bose distribution function (12).Now, going back to the quasiequilibrium statistical operator (1), (2), we suppose that β(x) and u µ (x) are slowly varying functions across the three-dimensional hypersurface σ.This makes it possible to apply a local thermal equilibrium approximation (see, e.g., Refs.[7,12]) of the statistical operator (1), (2).The local thermal equilibrium is an approximate concept which is usually associated with the possibility of defining a fluid cell, i.e., with the existence of a scale at which the system appears to be at homogeneous equilibrium.Therefore, this scale should be much smaller than the distance over which the β µ (x) = β(x)u µ (x) varies essentially.On the other hand, this scale has to be assumed large enough from a microscopic point of view, meaning that the typical microscopic correlation lengths are much smaller than the size of a cell.
To avoid additional complications and formulate the idea more concretely, we restrict ourselves to the case when the timelike normal vector n µ (x) of the hypersurface σ coincides with the 4-velocity field u µ (x), Then, we replace the integral in Eq. ( 2) by the sum as where and the integral in the above equation is taken over the homogeneity region of the β ν (x) around some point x µ s .The homogeneity region is defined as a region of the three-dimensional hypersurface σ where β ν (x) does not vary in a noticeable way.It is instructive to rewrite β ν (x s )P ν (σ s ) in the comoving coordinate system where ũµ (x s ) = (1, 0).Then, and ts = const.The key assumption underlying the local-equilibrium approximation is that characteristic size, L, of the corresponding volume element is large enough, i.e., L ≫ 1/m.
This assumption has important consequences.In particular, by using Eqs.( 3), (4), and (5), one can show that contributions of aa and a † a † terms to the P 0 (σ s ) can be neglected.
In a sense, this provides the local thermal equilibrium in the region s around x µ s . 2 The corresponding local-equilibrium statistical operator is where P µ reg (σ s ) = P µ (σ s ) − 0|P µ (σ s )|0 .By using Eqs.( 3), (4), and (5), we get Going to the local rest frame for a cell, we can write Eq. (29) in the following form: 2 Then, in particular, an ideal fluid approximation with a corresponding form of the energy-momentum tensor is approximately valid; see, e.g., Ref. [7].For quasiequilibrium states characterized by strong β µ (x) gradients, corrections to local thermal equilibrium approximation and, thereby, to ideal fluid approximation are sizeable and need to be taken into account.
Equation (29) makes possible to rewrite the operator s β ν (σ s )P ν reg (σ s ) as where In this section, we calculate one-particle and two-particle momentum spectra in the grandcanonical ensemble, which is described by the local-equilibrium statistical operator.For this aim, it is convenient to compute fist a † (p 1 )a(p 2 ) , where ... = Tr[...ρ leq (σ)].It can be done by adapting the Gaudin's method to our problem.We start by defining a(p, α), where s β ν (σ s )P ν reg (σ s ) is defined by Eqs.(31), (32), and (33).Applying the operator identity we can write the result as Taking into account the nonoverlapping of different cells and neglecting the surface effect on the boundaries of neighboring cells, we get Substituting Eq. (33) into Eq.( 37) and going to the local rest frame for each cell, we can perform an approximate integration over the momenta assuming that L ≫ 1/m, where L is the characteristic length scale of a cell.The result written in the laboratory coordinate system is Taking into account that main contribution in the integral over k is given by k ≈ p, it is where we introduced notation Note that the first term in Eq. (39) may be written as Substituting (41) into (39), we obtain It is convenient to introduce notation and rewrite Eq. ( 42) in the form From Eq. ( 43), we have that and that3 We can now employ the cyclic invariance of the trace to get expression for a † (p 1 )a(p 2 ) .
By using Eq.(34) and the cyclic invariance of the trace, one can write Here, a(p 2 , 1) is given by Eq. (44).By using Eq.(44), we can write Then, Eq. (47) becomes This equation can be solved by iteration.The result is Taking into account (46) we get where G * n is given by Eq. ( 43).Substituting G * n into Eq.( 51), we have Our next step is to replace sums over cells with integral over the hypersurface σ.This leads to where p µ = (p µ 1 + p µ 2 )/2.One-particle momentum spectra then read where f leq (x, p) is the grand-canonical distribution function, which has the familiar form of the local-equilibrium distribution function of the relativistic ideal gas of bosons, It is worth noting that our derivation can be readily extended to the local-equilibrium grand-canonical ensemble with nonzero constant chemical potential, µ, associated with mean number of particles.Then, Evidently, our derivation is rather heuristic and nonrigorous.But, in our opinion, it is instructive and adds some insights into the consistency of the approximations needed to associate quasiequilibrium statistical operator ρ q with the local-equilibrium ideal Bose gas distribution f leq (x, p).
Proceeding in the same way as above, one can readily derive an expression for the twoparticle momentum spectra, We start by using the cyclic invariance of the trace.This leads to The above equation is solved by iteration.We obtain Using Eq. ( 51), we can write the result as where a † a are given by Eq. (54).Equation ( 60) is the particular case of the thermal Wick's theorem [13].
Results of this section, in particular Eq.(53), will be used in the next section to evaluate particle momentum spectra and correlations at a fixed particle number constraint.

IV. QUANTUM LOCAL-EQUILIBRIUM CANONICAL ENSEMBLE WITH FIXED PARTICLE NUMBER CONSTRAINT
We begin this section by defining the local-equilibrium canonical ensemble with a fixed particle number constraint as a subensemble of the corresponding grand-canonical ensemble.
For this aim, we apply the constraint to the statistical operator given by Eq. ( 28).It implies utilization of the projection operator P N , which automatically invokes the corresponding constraint.Then, the local equilibrium statistical operator with the constraint, and we define ... N = Tr[ρ N ....].To evaluate two-boson momentum spectra at a fixed multiplicity, a † (p 1 )a † (p 2 )a(p 1 )a(p 2 ) N , we will follow the same strategy as in the previous section.We begin with evaluation of a † (p 1 )a † (p 2 ) N .This can be done by using its invariance under cyclic permutations.One gets Utilizing elementary operator algebra, one can prove that We also have where a(p 2 , 1) and ρ(σ) are defined by Eqs. ( 34) and (28), respectively.Therefore the r.h.s.
of Eq. (66) can be rewritten as Next, using Eqs.( 44) and (67), we obtain a(p 2 , 1)P N = P N −1 a(p 2 , 1).Therefore, Furthermore, accounting for Eq. ( 44) one can see that and that Substituting Eqs. ( 71) and (72) into Eq.( 70) and then into the r.h.s. of Eq. (66), we finally obtain the iteration relation, which yields Now, let us derive an expression for the two-particle momentum spectra at a fixed particle number constraint, Using the cyclic invariance of the trace, we obtain One can prove by induction that Combining Eqs. ( 74) and (77), we finally obtain It is immediately apparent from the above expression that the computation of a † (p 1 )a † (p 2 )a(p 1 )a(p 2 ) N involves summations over n and s, and these summations fail to factorize.Then, the question may arise as to whether this expression is invariant with respect to permutation of particles, i.e., with respect to permutation p 1 ↔ p 2 .To address this question, let us note that sums N −1 n=1 N −n s=1 can be rewritten as N −1 s=1 N −s n=1 .This means that Eq. ( 78) is invariant with respect to permutation s ↔ n and, therefore, is invariant with respect to permutation p 1 ↔ p 2 .
To evaluate Eqs. ( 74) and (78), we need explicit expressions for G * n (p 2 , p 1 , σ) and the partition functions Z n (σ).The former has been evaluated in the previous section; see Eq. (53).As for the latter, it can be evaluated as follows.First, note that the definition of ρ N (σ) means that Then, accounting for Eq. ( 74), we get the recursive formula where Z 0 = 1 by definition.
It is now a simple matter to write explicit expressions for the one-and two-particle momentum spectra.First, using Eqs.( 53) and (80), we get the recurrence relation that can be easily implemented numerically, Substituting Eq. ( 53) into Eq.( 74), we get Consequently, the one-particle momentum spectra at a fixed multiplicity constraint take the form where f leq N (x, p) is the local-equilibrium canonical distribution function at fixed N, It is worth noting that constant chemical potential of the grand-canonical ensemble, whose subensemble is the canonical fixed-N ensemble, does not influence on particle momentum spectra and correlations calculated at fixed multiplicity.It follows from the recurrence and therefore e βµn is factored out from expressions for particle momentum spectra and correlations.
Comparing Eq. ( 84) with Eq. ( 56), one can conclude that the selection of a fixed-N subensemble of the corresponding local-equilibrium grand-canonical ensemble results in nontrivial modifications of distribution functions.In particular, the one-particle distribution function (84) demonstrates multiplicity-dependent deviations in spacetime and momentum dependencies from the familiar local-equilibrium Bose ideal gas distribution function; see Eq. ( 56).In the next section, we compare particle momentum spectra and correlations calculated in the local-equilibrium grand-canonical and canonical ensembles for some simple but reliable for p + p collisions model.

V. PARTICLE MOMENTUM SPECTRA AND CORRELATIONS: COMPARI-SON OF THE ENSEMBLES
It is instructive to compare our findings with the treatment which is based on the grandcanonical ensembles (GCE) where chemical potential, µ = const < m, is taken such that mean particle number, N , is equal to particle number, N, in the canonical ensembles (CE) with a fixed multiplicity.Such an approach is often used for the sake of calculational convenience.Our simulations are performed for a simple hydro-inspired [14] local-equilibrium model of the longitudinally boost-invariant expanding system.In this model, the longitudinal direction (Z axis) coincides with the beam direction, and the 4-velocity is given by5 where τ = √ t 2 − z 2 is the proper time.We assume that a local-equilibrium state is defined at a hypersurface with constant energy density in the comoving coordinate system.Then, β(x) is constant on the corresponding hypersurface, and such a three-dimensional hypersurface σ is defined by a constant τ .It is convenient to parametrize t and z at this hypersurface as where η is the longitudinal spatial rapidity, tanh η = v L , and v L = z/t is the longitudinal velocity.This implies that where r T = (r x , r y ) are the transverse Cartesian coordinates.
This picture of an ultrarelativistic collision is, of course, not valid for large values of the spatial rapidity and for large transverse distances.We assume that the system has a finite transverse size encoded in the limits of integration over r T : 0 < r T < R T .As for the longitudinal direction, the finiteness of the system is provided by limits of integration over spatial rapidity η: The on-mass-shell particle 4-momentum p µ can be expressed through the momentum rapidity y, tanh y = p z /p 0 ; transverse momentum p T ; and transverse mass m T = p 2 T + m 2 , Then, For specificity and in order to compare the ensembles at the extreme small-system limits, we utilize for numerical calculations the set of parameters corresponding roughly to the values at the system's breakup in p + p collisions at the LHC energies.We take the particle's mass as of a charged pion, m = 139.57MeV, and the temperature T = 150 MeV (then the inverse temperature β = 1/T = 1/150 MeV −1 ).For τ , we use 1.5 fm/c.To account for finiteness of the system we assume that R T = 2 fm and η f = 2.
One-particle momentum spectra in the canonical ensemble with fixed multiplicity constraint, p 0 d 3 N d 3 p , are calculated utilizing Eqs. ( 81), (83), and (84): Employing the longitudinally boost invariant parametrization, we get where Z n are defined by the recurrence relation One-particle momentum spectra in the grand-canonical ensemble with N = N are calculated utilizing Eqs. ( 55) and (56) after substituting Then, For the considered model it implies that We now turn to the two-particle momentum correlations.The two-particle momentum correlation function at fixed multiplicities is defined as ratio of two-particle momentum spectrum to one-particle ones and in the canonical ensemble with fixed particle number constraint can be evaluated as where a † (p 1 )a † (p 2 )a(p 1 )a(p 2 ) N and a † (p 1 )a(p 1 ) N are defined in Eqs. ( 53), ( 78), (83), and (84).Here, G N is the normalization constant.The latter is needed to normalize the theoretical correlation function in accordance with normalization that is applied by experimentalists: It is convenient to evaluate the correlation function in terms of the relative momentum q = p 2 −p 1 and the pair momentum k = (p 1 +p 2 )/2.The correlation function takes a particular simple form for pairs with vanishing longitudinal pair momentum k z = (p 1z + p 2z )/2 = 0 and with k T = p 1T = p 2T , where k T is the pair momentum projected onto the transverse plane.Then, momentum rapidities of the particles in pairs are y 1 = −y 2 .Explicitly, the longitudinal projection (q T = 0) of the correlation function, C N (k T , q L ) (the subscript L ="long" indicates the longitudinal direction), is given by where (2) Here Ψ(n, m T , y 1 , q L ) = σ τ dηm T cosh η cosh y 1 e iq L τ sinh η e nβm T cosh η cosh y It is of interest to estimate the significance of the differences between the correlation functions calculated in the canonical and grand-canonical ensembles.We define the correlation function in the grand-canonical ensemble with N = N as where p µ = (p µ 1 + p µ 2 )/2.Then, To compare the ensembles, we begin by calculating N as function of µ/m.The results are presented in Fig.Then, to calculate particle momentum spectra and correlations in the canonical ensembles, we need to evaluate Z N for various N.The results are plotted in Fig. 2. Now, we are ready to compare spectra and correlations calculated in the grand-canonical and canonical ensembles.First, we compare particle number rapidity densities, dN/dy, for N = N.As illustrated by Fig. 3, the grand-canonical particle number rapidity densities are virtually indistinguishable from their canonical counterparts.
The transverse particle momentum spectra are compared in Fig. 4. Aside from the low transverse momenta region of the spectra with N = N = 11, where the grand-  canonical spectrum is above the canonical one due to the Bose-Einstein enhancement (µ is approximately equal to m; see Fig. 1), we see no significant differences.It distinguishes two-boson correlation functions in the canonical ensembles from the ones in the corresponding grand-canonical ensembles where the correlation functions (not shown here) approach to 1 from above and the intercepts are equal to 2.
Notwithstanding the essential non-Gaussianity of the canonical correlation functions, if the fitting procedure is restricted to the correlation peak region, then the correlation function is well fitted by the Gaussian expression It is instructive to compare canonical radius parameters extracted according to this expression with the ones calculated in the grand-canonical ensembles for N = N, For definiteness, for both ensembles we apply the fitting procedures in the q L range 0 < q L < q max L , where q max L is such that e −q 2 L R 2 long = 0.4.Our results are depicted in Figs.7 and 8.One can see that the canonical radius parameters slightly decrease with N, and the same trend, i.e., a decrease with N = N, is also observed for the grand-canonical radius parameters that are slightly smaller than the canonical ones.This decrease with N = N can be interpreted as increasing deviations from the Boltzmann approximation. Figure 8 shows R long as a function on k T = |k T | for several different values of N. One can see that R long in both ensembles is much smaller than the actual longitudinal size of the system (∼ τ sinh η f ) and decreases when k T increases.Such a smallness of the correlation radius parameters and a decline with increasing pair momentum are typical for locally equilibrated expanding systems [1].In Fig. 8 we plot for comparison the approximate analytical formula for R long , 2βm T [15].The latter approximate equality is obtained by means of the asymptotic (large argument) expansion of the Macdonald functions.
In the limit, βm T ≫ 1, this reduces to the formula R long ≈ τ 1 βm T [16] (see also [17]).All of the two figures reveal a consistent trend: if radius parameters are fitted in the region of the correlation peak, then deviations of the canonical radius parameters from their grandcanonical counterparts are rather small.It is also the case for N = N = 11 and small k T , because the effects of the Bose-Einstein enhancement (µ is approximately equal to m at N = 11) are nearly canceled out in the ratio (104).See the text for details.

VI. CONCLUSIONS
In this paper, we derived analytical expressions for one-and two-particle momentum spectra of a noninteracting relativistic boson field in the canonical ensemble described by the local-equilibrium statistical operator with a fixed particle number constraint.To see the effect of this constraint, we considered a corresponding grand-canonical state and compared the one-particle spectra and two-particle Bose-Einstein correlation functions.The correspondence was fixed by the condition that particle numbers, N, in the canonical states and mean particle numbers, N , in the grand-canonical states are the same.Then, applying hydrodynamically motivated parametrization and parameter values that correspond roughly to the values at the system's breakup in p + p collisions at the LHC energies, we compare our results with the grand-canonical ensemble where artificial chemical potential, µ = const < m, is taken such that N = N.We have found that, calculated in both ensembles, one-particle momentum spectra are rather close to each other except for low transverse momenta region of the spectra with N = N = 11, where the grand-canonical spectrum is above the canonical one due to the Bose-Einstein enhancement (µ is approximately equal to m).Then, we compared the two-particle Bose-Einstein momentum correlations.We demonstrated that there are small quantitative but qualitative differences between the correlation radius parameters in both ensembles if they are fitted in the region of the correlation peak: the canonical radius parameters are slightly larger than the grand-canonical ones.Furthermore, we showed that, in contrast to the predictions of the grand-canonical ensemble, the intercepts of the canonical correlation functions are not equal to 2 and depend on particle multiplicities and momenta and that the canonical correlation functions can be less than unity in some intermediate region of relative momentum of particles.Such features should be taken into account when theoretical models are compared with the multiplicity-dependent measurements of the Bose-Einstein momentum correlations.As a final comment, we wish to note that the apparent independence of correlation radius parameters on the particle number densities in high-multiplicity p + p collisions at a fixed energy of the LHC [3,4] still remains unexplained, inviting further studies.

FIG. 1 :
FIG.1:The µ/m dependence on N .See the text for details.

FIG. 2 :
FIG.2:The Z N dependence on N .See the text for details.

5 FIG. 4 :
FIG. 4: The transverse momentum spectra calculated in the canonical and grand-canonical ensembles with different N = N .

Figures 5 and 6
Figures 5 and 6 display two-boson momentum correlation functions C N (k T , q L ) calculated in the canonical ensembles as a function of the momentum difference.From these figures, it is evident that the intercepts of the canonical correlation functions, C N (k T , 0), are not equal to 2 and that the canonical correlation functions approach to 1 from below when |q L | → ∞.

FIG. 5 :FIG. 6 :
FIG. 5: The canonical correlation functions for N = 5 and several different values of k T .

FIG. 7 :
FIG. 7: The R long dependence on N = N in the canonical and grand-canonical ensembles for several different values of k T .See the text for details.

FIG. 8 :
FIG.8:The R long dependence on m T in the canonical and grand-canonical ensembles for several different values of N = N and the R long calculated from the approximate analytic expression.
N .It can be realized by means of the limit |q L | → ∞ at fixed k T in the corresponding expression.One can readily see that proper normalization is reached if 1.One observes from this figure that µ is about m when N is near 11.Because we do not aim to calculate here the Bose-Einstein condensation in the grandcanonical and canonical ensembles, in what follows, we do not consider canonical ensembles with N larger than 11.