Bose-Einstein momentum correlations at fixed multiplicities: Lessons from an exactly solvable thermal model for $pp$ collisions at the LHC

Two-particle momentum correlations of $N$ identical bosons are studied in the quantum canonical ensemble. We define the latter as a properly selected subensemble of events associated with the grand canonical ensemble which is characterized by a constant temperature and a harmonic-trap chemical potential. The merits of this toy model are that it can be solved exactly, and that it demonstrates some interesting features revealed recently in small systems created in $p+p$ collisions at the LHC. We find that partial coherence can be observed in particle emission from completely thermal ensembles of events if instead of inclusive measurements one studies the two-boson distribution functions related to the events with particle numbers selected in some fixed multiplicity bins. The corresponding coherence effects increase with the multiplicity.


I. INTRODUCTION
Femtoscopic study results on the two-particle momentum correlations (see, e.g., Ref. [1]) in p + p collisions at the CERN Large Hadron Collider (LHC) have been presented recently by the ALICE [2], ATLAS [3], CMS [4], and LHCb [5] Collaborations. It was found that the femtoscopic radii measured by the ATLAS and CMS Collaborations decrease with the increasing momentum of a pair. It can be interpreted in the hydrodynamical approach as the decrease of "homogeneity lengths " [6] (sizes of the effective emission region) due to generation by the collective flow x − p correlations. Also, one found that the "correlation strength" parameter λ is essentially less than unity. This is at variance with the expected behavior for emission from thermalized systems [1].
Another very interesting observation is the saturation of the multiplicity dependence of the interferometry correlation radius parameters for very high charged-particle multiplicity.
Such an effect was observed recently by the ATLAS [3] and CMS [4] Collaborations. Then, while there is some evidence that hydrodynamics can be successfully applied to describe particle momentum spectra in high-multiplicity p + p collisions (for recent review see, e.g., Ref. [7]), it is still unclear whether the reported results on Bose-Einstein momentum correlations can be attributed to hydrodynamic evolution like in A + A collisions.
In our opinion, observed peculiarities of Bose-Einstein momentum correlations in highmultiplicity p+p collisions do not indicate inapplicability of hydrodynamics but can be partly associated with quantum coherence effects in small systems, when the effective system size is comparable with typical wavelength of the thermal bosons. Recall that the effective geometrical size is associated with the length of homogeneity in the system [6].
Recently, a detail analysis of inclusive spectra and Bose-Einstein correlations in small thermal quantum systems was done for the analytically solved model in Ref. [8]. It is shown that if one deals (even locally) with a grand canonical ensemble, a nontrivial coherence parameter appears in inclusive two-boson spectra only in the case of coherent condensate formation. Without the latter, no coherence-induced suppression of the inclusive correlation function is possible because of the thermal Wick's theorem.
As for nonthermal or quasithermal emission with fixed particle multiplicity, the traditional pair-correlation function is distorted for events with high phase-space density, in particular, suppression of the Bose-Einstein correlations arises. The special algorithms for symmetriza-tion of multiboson N-particle states with independent particle emissions, and subsequent calculations of one-and two-particle spectra were developed in Refs. [9][10][11][12][13][14]. The situation, when particle radiation from different source points are not independent because the wave packets of emitted bosons are overlapping, was considered in Ref. [15].
Coming back to the thermal sources, in Ref. [16] the coherence effects in Bose-Einstein correlation functions in thermal systems are studied in subensembles of events with fixed multiplicities. The analytical calculations were done in the corresponding canonical ensemble. It was found that the correlation functions are suppressed in a finite system in a large volume and low particle number density approximation. In the present paper, we study the two-boson momentum correlations in small systems with high particle number densities at the moment when the system breaks up. Such almost sudden freeze-out can happen due to very fast expansion (when the homogeneity lengths are around 1 fm) of the matter formed in high-multiplicity p + p collisions at the LHC. To make the problem tractable we utilize a model of the finite system with smooth edges to avoid strong boundary effects. Keeping in mind the collective expansion inherent to systems created in particle and nucleus collisions, one can associate the corresponding system's scale-parameter with the homogeneity length.
The particle momentum spectra at a sharp freeze-out are formed according to Ref. [17], which is a reasonable approximation for p + p collisions. In order to keep things as simple as possible, we consider nonrelativistic ideal gas of bosons at fixed temperature trapped by means of a harmonic chemical potential. Such an exactly solvable toy model of inhomogeneous and finite-sized systems is mathematically identical to an ideal bosonic gas trapped by a harmonic potential. Then we apply the fixed particle number constraint to the corresponding grand-canonical statistical operator and discuss the influence of such constraints on one-particle momentum spectra and two-boson momentum correlations.

NUMBER CONSTRAINT
We begin with a brief overview of the properties of the grand-canonical ensemble of noninteracting nonrelativistic quantum-field bosons at fixed temperature, T , trapped by a harmonic chemical potential. For such a quantum field the Hamiltonian is given by where the operators Ψ † (r) and Ψ(r) are the creation and annihilation operators, respectively.
They fulfill the commutation relations and The Fourier transformed operators are defined as They satisfy the following canonical commutation relations: and The grand-canonical ensemble of such a system can be represented by the thermal statistical operator ρ, where Z is the grand-canonical partition function, andρ = e −β H , where β = 1/T is inverse temperature. The chemical potential, µ(r), reads whereμ = const. The expectation value of an operator O can be expressed as It is well known that H is not diagonal in momentum (plane-wave) representation but can be diagonalized in the oscillator representation. Decomposing Ψ(r) and Ψ † (r) in terms of the harmonic oscillator eigenfunctions we get where the creation, α † (n, k, l), and annihilation, α(n, k, l), operators satisfy the commutation and [α(n, k, l), α(n ′ , k ′ , l ′ )] = [α † (n, k, l), α † (n ′ , k ′ , l ′ )] = 0.
Functions φ n (x), φ k (y), φ l (z) are the harmonic oscillator eigenfunctions satisfying corresponding equations, e.g., The normalized solution of Eq. (17) reads where H n (x/b x ) is the Hermite polynomial, and Eigenfunctions (18) are complete, and orthonormal Then, from Eq. (14) it immediately follows that In such a basis the H reads Equation (24) of the particle number operator j α † (j)α(j), and the identity which express the completeness and normalization of this basis, one can insert Eq. (24) into Eq. (10) and writeρ in the harmonic oscillator basis, We denote here Then, using an elementary operator algebra and Eq. (27) one can see that Using trace invariance under the cyclic permutation of an operator, we get The Kronecker delta in the above equation, δ j 1 j 2 , is From Eq. (30) we have which is a familiar Bose-Einstein distribution. It follows then that In a similar way, one can get Then, taking into account Eq. (32) we have which is nothing but the particular case of the thermal Wick's theorem. Then, utilizing Eq. Now, let us apply the fixed particle number constraint to the grand-canonical statistical operator (8) to define canonical statistical operator ρ N . For this aim, one can utilize the projection operator P N , which automatically invokes the corresponding constraint. Using Eqs. (14), (22) and (25) one can see that It is worth noting that such a projection is accompanied by the proper normalization in order to insure the probability interpretation of the ensemble obtained in result of this projection.
Then, using (37) we assert that the canonical statistical operator is whereρ and Z N is the corresponding canonical partition function, It follows from Eq. (39) that The vacuum state, N = 0, yields Z 0 = 0|0 = 1. Let us denoteρ N associated withμ = 0 asρ 0 N . Then one can readily see that ρ N = e β µNρ0 N and Therefore [see Eq. (38)] e βμN is factored out and ρ N does not depend onμ: The expectation value of an operator O is defined as It follows from Eqs. (39) and (44) that To evaluate the expectation values of operators α † (j 1 )α(j 2 ) and α † ( with the canonical statistical operator ρ N , one can adopt the procedure which was used above to calculate expectation values with the grand-canonical statistical operator ρ. It can be done in a similar way as it was done, e.g., in Ref. [16]. For the reader's convenience, below we adjust the derivation from Ref. [16] for our model. A starting point is the relation which follows from Eq. (39) and commutation relations (15) and (16). Then one can exploit invariance under cyclic permutation and get the iteration equation With the starting value α † (j 1 )α(j 2 ) 0 = 0 one can get from the above equation that It follows from the definition of ρ N [see Eqs. (38) and (39) Utilizing relation (46) we have Then the same procedure leads to One can show by induction that Eq. (51) can be written as Then, taking into account that α † (j 1 )α(j 2 ) 0 = 0 and Eq. (48), we get The above expressions explicitly demonstrate deviations from the Wick's theorem in the canonical ensemble for a system of noninteracting bosons.
Canonical partition functions in Eqs. (48) and (53) can be calculated by means of the recursive formula of the canonical partition function for a system of N noninteracting bosons as given in Ref. [19] (an elementary derivation of it can be seen in Ref. [16]): where Z 0 0 = 0|0 = 1 and n = 1, ..., N. As a final comment we would like to point out that there is an essential difference between states defined by the grand-canonical statistical operator, ρ [see Eqs.  (40)]. While the former is a mixture of all N-particle states including vacuum state with N = 0, the latter is a mixture of states with N fixed to some value. In a sense, the quantum canonical state, ρ N , can be interpreted as a state which is not completely chaotic but has some quantum coherent properties. In what follows we demonstrate that such a coherence is enhanced in the case of the Bose-Einstein condensation, when the number of particles in the ground state, N 0 , is of the order of the total number of particles, N, 2 and discuss possible relations of our results to two-boson momentum correlations measured in p + p collisions at the LHC.

MULTIPLICITIES
In this section we relate the model with physical observables in relativistic particle and nucleus collisions. To keep things as simple as possible, below we assume that ω x = ω y = ω z ≡ ω. Note that the mean particle number, N , defined by the grand canonical ensemble, as well as the particle number, N, in the canonical ensemble are the same for Ψ particles and α quasiparticles because unitary transformation (14) does not mix creation and annihilation operators.
First, let us estimate spatial size of the system at fixed multiplicities. It is defined as Ψ † (r)Ψ(r) N is the mean particle number density in the canonical ensemble, and dxdydz Ψ † (r)Ψ(r) N = N. From Eqs. (14) and (48) we get where the eigenfunctions are defined by Eq. (18), see Eq. (20). Then, utilizing integral representation of the Hermite function (see e.g. Ref. [21]), one can perform summations over n, k, l in Eq. (56). A lengthy but straightforward calculation results in .
Utilizing identity (tanh A) −1 − (sinh A) −1 = tanh(A/2), we have from Eq. (59) that mean particle number density in the canonical ensemble reads Substituting the above expression in Eq. (55) we readily find To relate parameters of the model with physically meaningful parameters in relativistic particle and nucleus collisions, it is convenient to introduce parameter R such as then mω 2 2 = 1 2βR 2 ; see Eq. (12). In what follows we treat R as free parameter instead of ω. As we will see below, R can be approximately associated with the spatial size of the system, where Λ T is the thermal wavelength, which we defined as We now turn to the two-particle momentum correlation functions. Two-particle momentum correlation function is defined as ratio of two-particle momentum spectrum to one-particle ones and can be written in canonical ensemble at fixed multiplicities as Here k = (p 1 + p 2 )/2, q = p 2 − p 1 , and 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: C exp (k, q) → 1 for |q| → ∞.
Expressions in the denominator of Eq. (66) can be written immediately using Fourier transform of Ψ † (r 1 )Ψ(r 2 ) N ; see Eq. (59). We thus have where we introduced shorthand notation Utilizing Eq. (53) and the same technique which was used to derive Ψ † (r 1 )Ψ(r 2 ) N , we get after somewhat lengthy but straightforward calculations, where we introduced notation Φ 2 (k, q, βωs) = b 3 (2π sinh(βωs)) 3/2 exp −k 2 b 2 tanh( To estimate normalization constant G N in Eq. (72), one needs to utilize the limit |q| → ∞ at fixed k in the corresponding expression. One can readily see that when |q| → ∞ at fixed . It follows then that proper normalization is reached (73)

IV. RESULTS AND DISCUSSION
In this section, we calculate one-particle momentum spectra and two-particle Bose-Einstein momentum correlations in the model. For specificity, we assume that m is equal to pion mass and we take the set of parameters corresponding roughly to the values at the system's breakup in p + p collisions at the LHC energies: The temperature T is set to 150 MeV, and for R we use 1.5 and 3 fm. The thermal wavelength Λ T = 1/ √ mT ≈ 1.36 fm.
We varied N in the range 1, ..., 20. Our aim here is to investigate how particle momentum spectra and correlations in the canonical ensemble with the fixed particle number constraint differ from the ones in the corresponding grand-canonical ensemble.
We start with calculations of the one-particle momentum spectra in the canonical ensemble, n N (p) ≡ Ψ † (p)Ψ(p) N ; see Eq. (67). We compare these calculations with the ones performed in the corresponding grand-canonical ensembles whereμ were found numerically to guarantee proper values of N , such as N = N. One-particle momentum spectra in the grand-canonical ensembles are calculated utilizing Eq. (67) after substitution N s=1 The results are plotted in Fig. 1 as a function of the particle momentum for several different values of the radius parameter R and particle number N.  by the two-Gaussian expression where λ 1 > 0 and λ 2 > 0. Here 1 − λ 1 (k, N)e −q 2 R 2 1 (k,N ) is associated with the first term in Eq. (72), and λ 2 (k, N)e −q 2 R 2 2 (k,N ) with the second one. The results of fittings are plotted in Fig. 2. It is evident that C N (k, q) is rather well fitted by Eq. (74). This suggests that much of the non-Gaussian deviations observed in Fig. 2 arises from such a two-scale structure of the correlation function. If the fitting procedure is restricted to the correlation peak region, then one observes from Fig. 2 that the correlation function is well fitted by the one-Gaussian where λ is equal to the intercept of the correlation function, C N (k, 0). corresponds to the mean spatial size of the system in the varied range of N. It means that the decrease of R at fixed N results in an increase of the mean particle number density, ∝ N/R 3 .
It is well known that utilizing the above expression for approximate evaluation of the canonical partition function, in the leading order of the saddle-point approximation one obtains whereμ σ is solution of the equation d dμ (−μβN + ln Z(μ)) = 0. For an ideal gas it means that µ σ is such that N = N. Equation (77) becomes exact for N → ∞. Then, using Eqs. (41), (42), and (45), one can expect that for finite but large N we get N 0 /N ≈ N 0 / N where N 0 / N is the condensate fraction in the grand-canonical ensemble with N = N. Let us compare our results for the canonical condensate fraction, N 0 /N, with the grand-canonical condensate fraction for a finite mean number of particles in a three-dimensional harmonic potential [23], calculated in the approximation βω ≪ 1. Here ζ(x) is the Riemann zeta function, ζ(2) ≈ 1.645 and ζ(3) ≈ 1.202. In the above expression we have approximated e β(μ−(3/2)ω) ≈ 1.
Identifying N with the actual particle number N, and βω with Λ T /R, see Eq. (63), we compare N 0 /N with N 0 / N in Fig. 5 for R = 1.5 fm and R = 3 fm.
One observes that the approximate grand-canonical formula shows a rather good agreement with the exact canonical results even for not very large values of N. Loosely speaking, the canonical condensate fraction of the large system becomes noticeable (say, about 1/2) when the mean interparticle distance, (N/R 3 ) −1/3 , becomes smaller than the correlation length, for an ideal gas the latter coincides with the thermal wavelength, Λ T .
While the quantitatively accurate description of the canonical condensate fraction within the grand-canonical approximation is manifest, it is not the case for fluctuations. It is well known that fluctuations in the ground state differ in the canonical and grand-canonical ensembles, and that for the latter the condensate fluctuations are very large; see, e.g., Ref.
[24] and references therein. In the canonical ensemble with a fixed number of particles such large fluctuations are impossible, and therefore an increase of the ground-state fraction N 0 /N increases "coherence" of the state. The latter distinguishes the ideal gas Bose-Einstein condensation in the canonical ensemble from the ideal gas Bose-Einstein condensation in the grand-canonical ensemble. It is well known that the intercept of the two-boson momentum correlation function for a maximally mixed (chaotic) state is equal to 2, and that the one for a pure state is equal to 1; see, e.g., Ref. [1]. Therefore, an increase of the ground-state fraction, N 0 /N, results in a decreasing of the λ parameter.
Finally, in Fig. 6 we plot the R HBT as a function of the pair momenta, k, for different R and N. One observes a consistent trend: by increasing k the interferometry radii, R HBT , become independent of N.

V. CONCLUSIONS
Usually one does not care so much about quantum coherence in the canonical ensemble at fixed multiplicities. 3 However, utilizing the simple analytically solvable model, we demonstrated that the formulas derived in the fixed-N canonical ensemble for a small inhomogeneous thermal system are not always accurately approximated by the grand-canonical ones with N = N. Namely, we noticed that while the one-particle momentum spectra can be well approximated by the corresponding grand-canonical ensemble expressions, it is not the case for the two-boson momentum correlations. Interestingly, we observed that the most significant deviations arise if the particle number density in the canonical ensemble can increase with N. In the considered simple model it implies that interferometry radii are independent on N at moderately high pair momenta. Then for fairly high N the particle number density exceeds some limit value leading to the noticeable Bose-Einstein condensation in the corresponding ground state of the fixed-N canonical ensemble state. Such a condensation strengthens the coherence properties of the canonical ensemble state, and results in the decreasing of the intercept of the two-boson momentum correlation function when N increases. This may explain the observed phenomenon of partial quantum coher-ence in high-multiplicity p + p collisions events in fixed multiplicity bins at the LHC energies [3,4]. It would be very interesting to revisit the results of experimental studies in view of our findings.
The main lesson from this study is that the canonical and grand-canonical ensembles can yield different results for two-boson momentum correlations of particles emitted by small inhomogeneous systems. The results of our analysis can be useful to elucidate the influence on the shape of the measured correlation function of both factors: an experimental selection of events with fixed multiplicity and the effects of thermalization and flow. Therefore, determination of the extent to which our results can be generalized for a realistic model of heavy ion and particle collisions could be of great interest.