Inclusive spectra and Bose-Einstein correlations in small thermal quantum systems

The spectra and correlation of identical particles emitted from small locally equilibrium sources are considered. The size of the system is defined by negative part of parabolic falling chemical potential. The analytical solution of the problem is found for the case of inclusive measurements. It is shown that in the case when the size of the system is comparable to the thermal wave-length of the particles, the spectra and correlation functions are far from the quasi-classical approximation expecting for large systems, and observed femtoscopy scales (interferometry radii) will be essentially smaller than the Gaussian radii of the source. If the maximal value of chemical potential is closing to the critical one, specific for the system, one can consider a possibility of the Bose-Einstein condensation. Then the reduction of the intercept of the correlation function for inclusive measurements takes place. The results can be used for searching of femtoscopy homogeneity lengths in proton-proton collisions at the LHC energies.


I. INTRODUCTION
During a few last years an intensive femtoscopy study of proton-proton collisions at the LHC are provided by the CMS [1], ATLAS [2] and ALICE [3] Collaborations. Some interesting results, such as saturation of the femtoscopy scales for increasing particle multiplicities, peculiarities of the intercept behavior for correlation function, and also anti-correlations of identical pions were observed. A decrease of the interferometry radii with increase of pair transverse momentum in p + p collisions was found if specific selection of events (e.g., according to sphericity criteria [3]) is not performed. One of the interpretation of the radii behavior is the hydrodynamization of the systems created in very high energy p + p events with large multiplicities [1]. In this way a successive description of the femto-data on p + p collisions at √ s = 7 TeV in the Hydro-Kinetic Model (HKM) has been reached in Refs. [4].
This is one of the points that allows CMS Collaboration to interpret the obtained results [1] for √ s = 13 TeV as a consequence of hydrodynamic expansion of the thermal systems formed in p + p collisions at the LHC energies.
One of the peculiarities of the correlation femtoscopy for p+p collisions is the smallness of the strongly interacting system created in these processes: its typical size is about 1 fm that is comparable with mean wave-lengths of emitted particles. As it was considered in Ref. [5], the standard method of independent sources [6] is violated because of uncertainty principle: one cannot consider the emission of the particles from different parts of a small system as independent, if the wave packets from close emitting centers are essentially overlapping (overlap integral is not small). As the result, the visible interferometry scale is reduced as compare to geometrical system's size and correlation function is suppressed: its intercept decreases [5]. It is worth noting, that the approach to the problem of correlation femtoscopy for small system, developed in Ref. [5], that bring a good description of the 7 TeV p + p data [4], deasl, however, with events having small and fixed multiplicity, and does not use the hypotheses of thermalization.
In this note we propose the new results for correlation femtoscopy in analytically solved model of small thermal quantum systems. These results could be applied for inclusive correlation measurements of the homogeneity lengths [7,8] in p + p collisions at large mean multiplicities in the way similar to what is used in Ref. [4].

II. STATEMENT OF THE PROBLEM AND BASIC EQUATIONS
Let us consider a locally equilibrium system on some hypersurface σ with the time-like normal vector n µ (x), which is described by the statistical operatorρ: where S max (σ) is a maximum of entropy on σ under conditions fixed by the local distributions of energy, momentum and charges densities (see [9] for detail). For simplicity, we consider a real free scalar field in a (d + 1) dimensional space-time with the energy-momentum tensorT µν (x) = ∂ µφ ∂ νφ − 1 2 g µν ∂ ρφ ∂ ρφ − m 2φ2 and the current of particle number densitŷ J µ (x) = −iφ + ←→ ∂ µφ− (x), whereφ ± (x) are positive and negative frequency parts of the field which are defined as followŝ Then the statistical operator (1) takes form [9], [11], [12]: where β(x) = 1 T (x) and µ(x) are Lagrange multipliers, corresponding to the inverse temperature and the chemical potential respectively, and Z is the normalization factor making T r[ρ] = 1.
The creation and annihilation operators obey the commutation relations: To examine quantum effects in small systems we consider an exact-solved model without internal flows on the a hypersurface σ with a uniform temperature distribution T (x) = T in the moment of time t = 0. Corresponding to σ normal vector is n µ = (1, 0) so dσ µ = n µ d d x.
Using Eqs. (3) and (2) we obtain: In the non-relativistic limit energy and chemical potential can be decomposed as p 0 = m+ p 2 2m , µ(x) = m + µ 0 + µ ′ (x) (restrictions for chemical potential value will be discussed later). For the simplicity we take chemical potential in Gaussian form µ . It is easy to see that terms which contain m in the non-relativistic limit of (5) are reduced.
Following the Gaudin's idea [13], modified for the case of local equilibrium systems [10], we introduce new operators which depend on the dimensionless parameter α : The operatorÂ here is defined in the next way: Using new operators (6), an inclusive spectrum can be calculated [16]: Explicit dependence of the operator a + p (α) can be obtained from the next equation, which follows from definition (6): Substituting here the expression (7) and taking into account the commutation relations (4) we get: Here it is useful to represent coordinates x i in the form of derivative of the exponent with respect to momenta: Then integrating by parts two times over k ′ 2 allows to integrate over x i : It is our basic equation that allow to find solutions for inclusive thermal mean values a + k 1 a k 2 that define single-and double-particle spectra in the locally equilibrium systems with parabolic falling chemical potential.

III. ANALYTIC SOLUTION OF THE PROBLEM
Since the density matrixρ (5) acting on any state does not change its particle number, the solution of Eq. (13) can be expressed as an integral over all creation operators. Moreover, due to its linearity, the general solution can be written as where C n ( k, k ′ ) are solutions of oscillator-like equation: Since a + k (α = 0) = a + k , C n ( k, k ′ ) satisfy additional condition: From Eqs. (15) and (16) it follows that C n can be factorized: where δ i,j is Kronecker delta. Withal Eq. (15) allows separation of variables . So in terms of the variable k i it is Schrödinger equation for a harmonic oscillator. Its solution is represented by the Hermitian functions while Eq. (16) is a completeness of the orthonormal basis Eqs. (16)- (19) yield consist of d components n = {n 1 , n 2 , ..., n d } each of them runs from 0 to infinity, and Altogether the next notations are used in the paper: Now we are ready to write the solution (14) precisely: which allows us to calculate the inclusive spectrum (8) here we introduce a kernel M( k 1 , k 2 ) Eq. (23) with the commutation relations (4) leads to integral equation with a separable kernel with respect to the spatial components of momenta: The solution of this equation can be founded in the form 1 : with the recurrent equation on K (s) The kernel M i (k 1 , k) can be calculated from the definition (24) using Mehler's formula [14], [15]: One can verify that the solution of the Eq.(27) takes the form: Equations (26), (27), and (24) lead to the expression of the inclusive spectra, which in variables p = k 1 + k 2 2 and q = k 1 − k 2 takes the form: The total number of particles in system can be obtained after integration over momenta A necessary condition for convergence of the series is which gives a restriction for the maximum value of the chemical potential µ 0 : The corresponding Wigner function can be obtained in the following way: One can expect that at some kind of thermodynamic limit when the thermal wavelength of the emitting bosons much smaller than the homogeneity length -source size in our case, a quasi-classical limit for the Wigner function (36) should be reached. Naive expectation is that such a function takes the form of the Bose-Einstein distribution with the corresponding coordinate dependent chemical potential. To demonstrate this, one has to consider the thermal (Compton) wavelengths of boson quanta Λ T to be much smaller then the size of For the simplicity we investigate the isotropic case (R 1 = ... = R d = R). A linear approximation to the hyperbolic functions in (36) is applied for s less than some value s 0 , say, such that s 0 βω ≈ 1 2 . Another criterion of possibility to get the quasi-classical limit is non-positiveness of the chemical potential, µ 0 < 0. Then one can get from (36) Extending of the first sum up to infinity (and subtracting the added terms), we obtain a quasi-classical approximation with corrections that vanish in the thermodynamic limit when βµ 0 = const < 0 and βω = Λ T R → 0:

IV. FEMTOSCOPY ANALYSIS
To investigate correlations in our model we have to calculate the two-particle inclusive spectra [16] on the freeze-out hypersurface: The thermal average of four operators is typically reduced to a sum of products of twooperators averages by means of the Wick's theorem. It was generalized from globally to locally equilibrium thermal systems in [5], where the problem was reduced to the integral Aiming to show an importance of quantum effects in small systems for the femtoscopy analysis we compare the correlation functions in quantum and quasi-classical approaches.  The femtoscopic analysis of high-energy proton-proton collisions accompanied by relatively large multiplicities has to be carried out more carefully, since the particle density is high (µ → µ max ), and many bosons occupy the lowest energy level. In Appendix C the average number of particles on this level is found and has the form (C9): which also follows from general Eq. (32), when the thermal wavelengths of quanta exceed the geometrical size of the system Λ T R, so that one can use a linear approximation for sinh(nβω). In fact, when this condition is satisfied, most of particles occupy the lowest energy level. Note that their total number is strongly depends on the chemical potential µ 0 .
It is worth noting that the operator average for a lowest level (C8) (see Appendix C) does not contain momenta difference q, so that the correlation function (41) in one-level approximation is constant, C(q) = 2. That leads to the important consequence, which is shown in Fig. 2, specifically a broadening of the correlation function with increase of the chemical potential (so more and more particles occupy one level). It means that interferometry radii obtained from the Gaussian fit of the real correlation function can be noticeably smaller than these formally related to the isotropic Gaussian source, f G ∼ exp(−x 2 /2R 2 ) with naive CF for independent boson emission : C(q) = 1 + exp (−R 2 q 2 ). Now we approach to very important point: when occupation numbers at the ground state become dominate, it can lead to the significant overlap between wave packets of bosons, that leads to coherence in the system [5], [20], [21]. The strongly overlapping bosons can hardly be considered as fully independently emitted, and even rather small interaction between them can bring correlations of the phases of the wave packets [5]. Thereby small systems with high multiplicities (R 1 √ mT and µ 0 close to µ max ) have to be described by a density matrix of partially coherent thermal state 2 . According to the general idea described in [17] (see also [18]) creation (annihilation) operators at the post-thermal freeze-out stage split into a quantum part, associated with q-numbers (operators) b p and c-numbers: Since an area of our interests is small systems, it is reasonable to consider a fast freezeout scenario (d coh (p, t f ) = d p ) with simultaneous decay of boson field into free particles.
In the first approximation the lowest energy state (see details in Appendix C) on the post 2 Note that the grand canonical ensemble covers events with any number of particles, so the coherence develops only in some part of them. We estimate the appearance frequency for high-multiplicity events freeze-out stage is occupied by the coherent particles, while the quantum part corresponds to excitation in thermal medium: The Wick's theorem is applicable only for q-operators b p (one needs to substitute a p /a + p → b p /b + p in Appendix A). In order to decompose a product of the operators a p (43), an additional term, related to the coherent part of the operators has to be subtracted: It is easy to see that in this approach the intercept is less then two, and determined by the condensate contribution to the inclusive spectrum. The latter is controlled by the chemical potential and the ratio of the thermal wave length to the size of the system Λ T R = βω. accounting for coherence are calculated with the same parameters as for the Fig. 2 are shown in Fig. 3. Note, that the Λ T R ratio is fixed. However as was mentioned before, for smaller systems the contribution from the ground state is larger, so the intercept of the correlation functions goes down, and vice verse for larger systems.
V. CONCLUSIONS In this paper, we have studied the Bose-Einstein correlations in small locally-equilibrium systems in a simple model having exact analytic solution. We have considered a free scalar field on the freeze-out hypersurface with the uniform temperature. It is shown, that in systems, comparable in size with thermal wavelength of emitted bosons, quantum corrections to the two-particle correlation functions of identical particles are substantial. Qualitatively, interferometry radii of the considered systems is smaller than those formally related to the Gaussian source with the same radii as geometrical sizes of the system. This difference increases in systems with higher multiplicities.
In the case of strong overlapping of the wave packets in the ground state in the most part of the events -the overlapping that happens because of thermal wave-length of the quanta are larger or similar comparing with geometric size of the system and/or due to chemical potential becomes close to its maximal value -the coherent Bose-Einstein condensate can appear. It leads to reduction of the intercept of the inclusive correlation function.
Note that the effect of reduction of the femto-scales and suppression of the correlation functions comparing with naive picture of independent boson emission from the Gaussian source of the same effective size, was found in non-thermal model in Ref. [5]. The exact analytic solutions found in this note for the model of small thermal source are planning to be applied for analysis of femtoscopic phenomena in p + p collisions at the LHC.

VI. ACKNOWLEDGEMENT
The research was carried out within the scope of the EUREA: European Ultra Relativistic Energies Agreement (European Research Network "Heavy ions at ultrarelativistic energies").
N-particle "distribution functions" f N ( p 1 , p 2 ) = a + p 1 a p 2 N . To establish this relation we expand the grand canonical ensemble in a set of N-particle canonical ensembles by means of projection operators P (N ) Inserting the following completeness relation note that T r(ρP (N) ) T r(ρ) = p(N) is a probability that system consist of N particles, and T r(ρa + p 1 a p 2 P (N) ) T r(ρP (N) ) = f N ( p 1 , p 2 ) is a distribution function in canonical ensemble. Then If p 1 = p 2 = p, after integration over the momentum p we get obvious relation Similarly to the calculations in the main part of the article we get an integral equation for T r ρa + p 1 a p 2 P (N ) . For this aim we use the following permutation relation Strait calculations (see Eqs. (23)-(25)) lead to For the case N = 1 (f 0 ( k 1 , k 2 ) = 0) we can see that or in terms of the number of particles in the system which can be used as a recurrent equation for p(N). We can also derive Eq.(26) summing up by the N Eq. (B9): Appendix C: Ground state spectrum In order to obtain a condensate contribution to spectrum and correlations we select bosons in the ground state (zero level of harmonic oscillator) from the grand canonical ensemble (5) by means of a projection on this state operator P (0) a + k 1 a k 2 cond = a + k 1 a k 2 P (0) = a k 2 P (0) a + k 1 (α = 1) , 1 N! |n 1 = m, ..., n N = m n 1 = m, ..., n N = m| (C2) The general idea in further calculations is to get an integral equation similar to (25) but with a different kernel M cond ( k 1 , k 2 ). For this we need to determine permutation relation between the creation operator a + k and the projection operator P (n) . Explicit calculations gives Due to the orthonormality of Hermitian functions (18) the recurrent relations (27) simplifies to the following: which together with Eq. (18) yields: N 0 = dp a + p a p cond = 1 e β(µmax−µ 0 ) − 1 (C9) As soon as microscopical description of the lowest energy state is relevant only when it is occupied by large number of particles N 0 , it is reasonable to find distribution of this number in the grand canonical ensemble (5). For this purpose we can use Eq. (B10), but with the kernels (C6), which is the same as analytical continuation of the low temperature limit. Assume that the overlap between the various wave-packets became sufficiently enough when the number of bosons on the same energy level exceeds some critical value N cr . The appearance frequency P coh of events with such property can be obtained by the next sum: which approaches unity when N 0 ≫ 1.