# Bose-Einstein correlations in high energy multiple particle production processes

## 1 INTRODUCTION

In this paper we will consider Bose-Einstein correlations in high-energy particle production processes, i.e. the correlations among identical bosons in the final state, which follow from Bose-Einstein statistics. When hundreds of identical bosons are being produced in a single scattering act, as happens e.g. in heavy ion collisions at high energy, such correlations can lead to spectacular phenomena. They are also, most probably, the best way of getting information about the space-time structure of the region, where the final state particles are produced. Let us begin with a very simple example.

Consider the elastic scattering of two alpha particles with initial momenta equal in magnitude, opposite and parallel to a horizontal axis, say the $`x`$-axis. Suppose that the detectors register the final state particles if and only if the scattering is at $`90^{}`$ and one of the final particles goes up and hits the upper detector ($`U`$), while the other goes down and hits the lower detector ($`L`$). There are two possibilities. Either particle $`1`$, say the particle coming from the left, hits detector $`U`$ and particle $`2`$ hits detector $`L`$, or particle $`1`$ hits detector $`L`$ and particle $`2`$ hits detector $`U`$. Let us denote the probability amplitudes for these two processes by $`A`$ and $`B`$ respectively. Since a rotation around the $`x`$-axis can convert these two processes into each other, $`|A|=|B|`$. Since the alpha particles are identical bosons and an exchange of the two final state particles converts the two processes into each other, $`A=B`$. The detection probability is $`|A+B|^2=4|A|^2`$. If the particles were distinguishable, the probability would be $`|A|^2+|B|^2=2|A|^2`$. Thus, the fact that the particles are indistinguishable increases the probability by a factor of two.

Let us make a few comments about this simple result.

* In the example the two amplitudes interfere constructively, because they are coherent. This is sometimes called first order interference. We will see in the following that the Bose-Einstein correlations of interest for us are due to the incoherence of the production process, and are a manifestation of the so called second order interference.
* The statement that the scattering probability for identical particles is twice the corresponding probability for distinguishable particles is not possible to check experimentally, because non identical alpha particles are not available. The best one can do is to compare the experimental result for the identical alpha particles with the calculation for the non identical ones. In the present example, where the calculation is simple and non controversial, this is not much of a problem, but in multiple production processes a calculation from first principles is not possible and the definition of the distribution for distinguishable particles, which should be modified by the Bose-Einstein correlations to yield the distribution which can be compared with experiment, is a difficulty.
* The final state can be represented by the density operator

$$\widehat{\rho }=\frac{1}{2}|U_1L_2+U_2L_1U_1L_2+U_2L_1|,$$
(1)   
where $`U_iL_k`$ is the state, where particle $`i`$ is registered by the detector $`U`$ and particle $`k`$ by detector $`L`$. Expanding the left-had-side one obtains four terms, if, however, the density operator is to be used only for calculating averages of operators symmetric with respect to exchanges of the identical particles, which is sufficient for all practical applications, one can use the simpler form

$$\widehat{\rho }=|U_1L_2U_1L_2|+|U_1L_2U_2L_1|.$$
(2)   
It is useful to rewrite this formula in the form

$$\widehat{\rho }=\underset{P}{}|U_1L_2U_{P1}L_{P2}|,$$
(3)   
where the summation is over all the permutations P of the indices $`1`$ and $`2`$ and $`Pi`$ denotes the index obtained from index $`i`$ under permutation $`P`$. For our simple example this formula is ridiculously complicated, for more difficult cases, however, its analogues are very convenient.

## 2 HBT contribution

An interesting application of the Bose-Einstein interference to find the sizes of the emitting objects was discovered in the fifties by two astronomers, R. Hanbury Brown and R.Q. Twiss . By studying the second order Bose- Einstein interference of photons, they were able to measure the radii of some stars. The idea may seem obvious today, but it was not so at the time it was put forward. In the seventies Hanbury Brown wrote  (quoted after )

> Now to a surprising number of people this idea seemed not only heretical but patently absurd and they told us so in person, by letter, in publications, and by actually doing experiments which claimed to show that we were wrong. At the most basic level they asked how, if photons are emitted at random in a thermal source can they appear in pairs at the two detectors. At a more sophisticated level the enraged physicist would brandish some sacred text, usually by Heitler, and point out that …our analysis was invalidated by the uncertainty relation …

Since the distances to many stars are known, their radii could be determined, if the opening angles between the light rays coming from the stars could be measured. These angles, however, are in most cases too small for a direct measurement. Hanbury Brown and Twiss suggested the following procedure. Consider two light rays coming from two points on the surface of the star – ray $`a`$ from point $`a`$ and ray $`b`$ from point $`b`$. The problem is to measure the angle $`\mathrm{\Theta }`$ between the two rays. Each of the rays falls on two photodetectors denoted $`1`$ and $`2`$. The distance between the photodetectors is $`d`$. Elementary trigonometry yields to first order in $`\mathrm{\Theta }`$ the relation

$$\mathrm{\Theta }=\frac{\mathrm{\Delta }_a\mathrm{\Delta }_b}{d\mathrm{sin}\alpha },$$
(4)   
where $`\mathrm{\Delta }_i`$ is the difference of distances between point $`i`$ on the star and the two photodetectors, while $`\alpha `$ is the angle between the line connecting the two photodetectors and the direction of the two rays. Thus the problem of measuring the opening angle $`\mathrm{\Theta }`$ reduces to the problem of measuring $`\mathrm{\Delta }_a\mathrm{\Delta }_b`$. Of course in practice, in order to find the radius of the star a suitable averaging over the possible emission points is necessary, this, however, is rather simple and we shall not discuss it any further.

The current generated in the photodetector is proportional to the intensity of the incident light. Thus for counter $`1`$ it is

$$i_{1u}=K_1\left[E_a\mathrm{sin}(\omega _at+\varphi _a)+E_b\mathrm{sin}(\omega _bt+\varphi _b)\right]^2,$$
(5)   
where $`\varphi _i,E_i,\omega _i`$ denote respectively the phase at the star surface, the amplitude and the frequency for ray $`i`$, and $`K_1`$ is a proportionality coefficient dependent on the working of the photodetector $`1`$. For simplicity, the polarization effects have been ignored and the time necessary to reach photodetector $`1`$ has been put equal $`t`$ for both point $`a`$ and point $`b`$. In the apparatus the current $`i_{1u}`$ is further filtered so that only frequencies from 1 Hz to 100 Hz survive. Thus finally, the current from the first photodetector is

$$i_1=K_1E_aE_b\mathrm{cos}[(\omega _a\omega _b)t+(\varphi _a\varphi _b)].$$
(6)   
This current is zero on the average and does not look particularly interesting. The analysis for the second photodetector is similar except that the time necessary to reach the detector for the ray from point $`i`$ on the star is increased by $`\mathrm{\Delta }_i/c`$ . One obtains

$$i_2=K_2E_aE_b\mathrm{cos}\left[(\omega _a\omega _b)t+\frac{\omega }{c}(\mathrm{\Delta }_a\mathrm{\Delta }_b)+(\varphi _a\varphi _b)\right],$$
(7)   
where $`\omega \omega _a\omega _b`$. This is again a rather uninteresting current, but the average of the product of the filtered currents from the two photodetectors

$$i_1i_2=K_1K_2E_aE_b\mathrm{cos}\left[\frac{\omega }{c}(\mathrm{\Delta }_b\mathrm{\Delta }_a)\right],$$
(8)   
which is measurable, yields $`\mathrm{\Delta }_a\mathrm{\Delta }_b`$ and consequently the necessary opening angle $`\mathrm{\Theta }`$.

Note that the result is obtained in spite of the fact that the presence of the random phases $`\varphi _i`$ means that light from $`a`$ is incoherent with respect to light from $`b`$. Because of these phases the product of two amplitudes, one for the ray $`a`$ and one for the ray $`b`$ averages to zero, The product of four amplitudes, two from $`a`$ and two from $`b`$, however, can survive. For this reason one calls this effect second order interference or intensity interferometry.

## 3 The GGLP contribution

The first application of intensity interferometry in particle physics was made by the Goldhabers Lee and Pais . Their problem was somewhat different from that that in the HBT case. The interfering particles were like sign pion produced at two points within the interaction region of a hadron - hadron collision. The interference of interest was not between the measurements at two points in space, but between momentum measurements. Assuming that at the production space-time point $`x_k`$ the pion wave function has phase $`\varphi _k`$ and that the momentum of the pion $`p_k`$ is well-defined, one expects at the registration point $`x`$ an amplitude proportional to $`\mathrm{exp}[ip_k(x_kx)+i\varphi _k]`$. The probability of finding the two pions produced at points $`x_1`$ and $`x_2`$ with momenta $`p_1`$ and $`p_2`$, after proper symmetrization of the wave function, is proportional to

$$\frac{1}{2}\left|e^{i(p_1x_1+p_2x_2)}+e^{i(p_1x_2+p_2x_1)}\right|^2=1+\mathrm{cos}[(p_1p_2)(x_1x_2)].$$
(9)   
We assume now that the production process is incoherent, so that the averaging over the times and positions $`x_1,x_2`$ should be made at the level of probabilities and not of amplitudes. Then the distribution of the difference in momenta should be approximately given by the formula

$$C(p_1p_2)=1+\mathrm{cos}[(p_1p_2)(x_1x_2)],$$
(10)   
where the averaging is over $`x_1`$ and $`x_2`$. Qualitatively, the result does not depend much on the actual prescription being used for the averaging. For $`p_1p_2`$ the argument of the cosine is close to zero and consequently $`C2`$. For large momentum differences, the argument of the cosine is a rapidly oscillating function of $`x_1x_2`$, which is strongly suppressed by the averaging process, and $`C1`$. If the weight function used for the averaging contains just one parameter with the dimension of length, let us denote it R, the width of region in $`p_1p_2`$, where $`C`$ is significantly bigger than one, must be of the order of $`R^1`$. There are many specific recipes how to perform the averaging. The results obtained for the correlations of momenta and for the sizes and shapes of the interaction regions are reasonable. For reviews see  and . In spite of this success many difficulties remain.

* Since the energy of a pion is determined by its momentum, one has data only on the three-dimensional distribution of the differences of spacial momenta. This is not enough to derive the four-dimensional distribution of the sources in space-time. Therefore, the results are strongly model dependent.
* The averaging over the square of the wave function corresponds to the assumption that the density matrix of the final pions in coordinate representation is diagonal. This in turn implies that the momentum distribution should be flat, which contradicts experiment. A closely related question is, how the pion can be initially localized at the production point and then represented by a plane wave corresponding to well-defined momentum.
* Information about the production region is, in practice, obtained only from pairs of pions with similar momenta. Consequently, what is being measured is not the whole interaction region, but the region, where the pions with similar momenta are produced. This, incidentally, explains the fact that the interaction regions usually come out roughly spherical, while one believes that the full interaction region is more string like.
* There are many corrections, which probably should be applied, but it is controversial how. Here belong the corrections for coulomb repulsion between the charged like sign pions, the corrections for the final state interactions due to strong coupling, the corrections due to resonance production, the corrections due to partial coherence of the source etc.

## 4 Density matrix approach

In order to obtain a more general formulation for the GGLP problem it is convenient to use the formalism of density matrices. This has been described by a number of people, here we are the closest to the formulation used by Bialas and Krzywicki . We introduce an auxiliary, unphysical process, where all the particles produced are distinguishable. We assume that for this process simple, intuitive ideas work. Then we correct for the Bose-Einstein correlations in order to obtain results comparable with experiment. This approach has its defects as discussed in the Introduction (see also , ), but for lack of a better idea it is widely used. As our starting point for the distinguishable particles we use an independent production model (cf , ,  and references contained there). In this model the multiplicity distribution for the particles is poissonian

$$P_N^{(0)}=\frac{\nu ^N}{N!}e^\nu $$
(11)   
and for each multiplicity the density matrix is a product of single particle density matrices

$$\rho _N^{(0)}(q,q^{})=\underset{i=1}{\overset{N}{}}\rho _1^{(0)}(q_i,q_i^{}).$$
(12)   
The momentum distribution is given, as usual, by the diagonal elements of the density matrix in the momentum representation

$$\mathrm{\Omega }_{0N}(q)=\rho _N^{(0)}(q,q).$$
(13)   
It is convenient to normalize it to unity

$$\sigma _N^{(0)}=𝑑q\mathrm{\Omega }_{0N}(q)=1.$$
(14)   
For identical particles the density matrix should be symmetrized as explained in the Introduction

$$\rho _N(q,q^{})=\underset{P}{}\rho _N^{(0)}(q,q_P^{}).$$
(15)   
The corresponding momentum distribution for a given multiplicity is

$$\mathrm{\Omega }_N=\rho _N(q,q).$$
(16)   
This, however, is no more normalized to unity, because

$$\sigma _N=𝑑q\mathrm{\Omega }_N(q)=1+\mathrm{}$$
(17)   
The first term in the last expression corresponds to the identity permutation, but there are $`(N1)!`$ further terms. This yields the multiplicity distribution

$$P_N=𝒩P_N^{(0)}\sigma _N,$$
(18)   
where $`𝒩`$ is an $`N`$-independent normalizing factor, which ensures that $`P_N=1`$.

## 5 Simple case: pure final state

In order to present simply the qualitative features of the result, let us consider first the case, when for each multiplicity the final state is pure. This is a grossly oversimplified model, but we will find that it contains some features of the much more realistic approach presented in the following section. For the pure state model

$$\widehat{\rho }_N^{(0)}=|\psi _N^{(0)}\psi _N^{(0)}|.$$
(19)   
We assume that each of the state vectors $`|\psi _N^{(0)}`$ is symmetric with respect to exchanges of particles. Thus the effect of the summation over the permutations $`P`$ is simply to multiply the operator $`\widehat{\rho }_N^{(0)}`$ by $`N!`$. As a result the probability of producing exactly $`N`$ particles is also multiplied by the factor $`N!`$. The Poisson distribution goes over into a geometrical distribution and after evaluating the normalization factor we get

$$P_N=(1\nu )\nu ^N.$$
(20)   
This formula makes sense only if $`\nu <1`$, because otherwise the sum of the probabilities $`P_N`$ diverges. For the average number of particles one finds

$$\overline{N}=\frac{\nu }{1\nu }$$
(21)   
with a singularity at $`\nu =1`$. From the model presented in the following section it will be seen that this singularity corresponds to Einstein’s condensation.

In order to avoid the repeated summation of series it is convenient to introduce the generating functions. The generating function for the multiplicity distribution is

$$\mathrm{\Phi }(z)=\underset{N=0}{\overset{\mathrm{}}{}}P_Nz^N=\frac{1\nu }{1z\nu }.$$
(22)   
The logarithmic derivative of this function with respect to $`z`$ at $`z=1`$ yields the average multiplicity. The second derivative of the logarithm with respect to $`z`$ at $`z=1`$ is the dispersion and in general the $`p`$-th cumulant of the multiplicity distribution is given by

$$K_p=(p1)!\left(\frac{^pLog\mathrm{\Phi }}{z^p}\right)_{z=1}=(p1)!\left(\frac{\nu }{1\nu }\right)^p.$$
(23)   
Inclusive and exclusive momentum distributions, as well as all the correlation functions in momentum space, can be calculated by functional differentiation from the generating functional

$$\mathrm{\Phi }[u]=\underset{N=0}{\overset{\mathrm{}}{}}𝒩P_N^{(0)}𝑑q\mathrm{\Omega }_N(q)\underset{i=1}{\overset{N}{}}u(q_i)=\frac{1\nu }{1\nu 𝑑q_i\mathrm{\Omega }_{01}(q_i)u(q_i)}.$$
(24)   
For instance, the single particle distribution is

$$\left(\frac{\delta \mathrm{\Phi }[u]}{\delta u}\right)_{u=1}=\frac{\nu }{1\nu }\mathrm{\Omega }_0(q).$$
(25)   
Thus symmetrization (Bose-Einstein statistics) introduces in this simple model only a change of normalization.

## 6 Independent production

Let us consider now the full independent production model. In order to find the modification of the multiplicity distribution due to Bose-Einsten statistics it is necessary to calculate the correction factors

$$\sigma _N=\underset{P}{}𝑑q\underset{i=1}{\overset{N}{}}\rho _1^{(0)}(q_i,q_{Pi}).$$
(26)   
Since each permutation can be decomposed into cycles, this integrals can be expressed in terms of the cycle integrals

$$C_{k>1}=d^{3k}q\rho _1^{(0)}(q_1,q_2)\rho _1^{(0)}(q_2,q_3)\mathrm{}\rho _1^{(0)}(q_k,q_1).$$
(27)   
It is convenient to add the definition

$$C_1=1.$$
(28)   
Similarly the integrals necessary to calculate the generating functional for the momentum distributions can be expressed in terms of the cycle integrals

$$C_k[u]=d^{3k}qu(q_1)\rho _1^{(0)}(q_1,q_2)u(q_2)\rho _1^{(0)}(q_2,q_3)\mathrm{}u(q_k)\rho _1^{(0)}(q_k,q_1)$$
(29)   
After some combinatorics, very similar to that used when deriving the linked clusters expansion familiar from quantum field theory and many body theory, one finds the generating functional

$$\mathrm{\Phi }[u]=\mathrm{exp}\left[\underset{k=1}{\overset{\mathrm{}}{}}\nu ^k\frac{C_k[u]C_k[1]}{k}\right].$$
(30)   
Substituting in this functional $`z`$ for $`u`$ one obtains the generating function for the multiplicity distribution. Without exhibiting the actual calculations we will now present some general results, obtained for this model.

* The single particle momentum distribution and all the momentum correlation functions can be expressed in terms of one function depending on two single particle momenta

$$L(q_1,q_1^{})=\underset{k=1}{\overset{\mathrm{}}{}}\nu _kd^3q_2\mathrm{}d^3q_k\rho _1^{(0)}(q_1,q_2)\rho _1^{(0)}(q_2,q_3)\mathrm{}\rho _1^{(0)}(q_k,q_1^{}).$$
(31)   
For instance the momentum distribution is

$$\mathrm{\Omega }(q)=L(q,q).$$
(32)   
The two particle cumulant is

$$K_2(q_1,q_2)=L(q_1,q_2)L(q_2,q_1).$$
(33)   
In general the $`p`$-th correlation function is

$`K_p(q_1,q_2,\mathrm{},q_p)=L(q_1,q_2)L(q_2,q_3)\mathrm{}L(q_p,q_1)+`$   
$`(\text{permutations of the indices }2,\mathrm{},p).`$ (34)
* For typical density matrices the average square of the difference of momenta between two particles $`q^2`$ decreases due to symmetrization.
* For typical density matrices the average difference between the production points of pairs of particles decreases due to symmetrization.
* For typical density matrices the size of the interaction region as evaluated from the width of the two-particle momentum correlation function decreases due to symmetrization.

We will discuss these predictions in a further section, where we will rederive them in a more intuitive way. The references to ”typical density matrices” mean that the statement is true for most density matrices, but not for all. We have not been able to find a condition defining the relevant class of density matrices.

Probably the most interesting implication is the possibility of Einstein’s condensation, but this will be discussed in the following section.

## 7 Einstein’s condensation

Using matrix notation one can rewrite the definition of the function $`L`$ given in the previous section in the form

$$L(q,q^{})=\underset{k=1}{\overset{\mathrm{}}{}}\nu ^kq|\left(\widehat{\rho }_1^{(0)}\right)^k|q^{}.$$
(35)   
Expanding the single particle density operator in terms of its eigenvectors and eigenvalues, we find

$$\widehat{\rho }_1^{(0)}=\underset{n}{}|n\lambda _nn|$$
(36)   
and for its $`k`$-th power

$$\left(\widehat{\rho }_1^{(0)}\right)^k=\underset{n}{}|n\lambda _n^kn|.$$
(37)   
Thus in the momentum representation

$$L(q,q^{})=\underset{n}{}q|nn|q^{}\underset{k=1}{\overset{\mathrm{}}{}}\nu ^k\lambda _n^k.$$
(38)   
This expression makes sense only if for all $`n`$ there is $`\lambda _n\nu <1`$. Denoting the largest eigenvalue of the density operator by $`\lambda _0`$, we expect problems when $`\nu \lambda _01`$.

Performing the summations of the geometric series, we can rewrite the expression for $`L(q,q^{})`$ in the form

$$L(q,q^{})=\underset{n}{}\frac{\psi _n(q)\psi _n^{}(q^{})\nu \lambda _n}{1\nu \lambda _n}.$$
(39)   
For $`\nu \lambda _01`$ it is convenient to use the equivalent formula

$$L(q,q^{})=\frac{\psi _0(q)\psi _0^{}(q^{})}{1\nu \lambda _0}+\stackrel{~}{L}(q,q^{}),$$
(40)   
where $`\stackrel{~}{L}`$ remains bounded in the limit. Putting $`q=q^{}`$ and integrating over $`q`$ we get the corresponding formula for the average multiplicity

$$\overline{N}=\frac{1}{1\nu \lambda _n}+\text{bounded term}.$$
(41)   
From these formulae it is clear that when $`\nu \lambda _0`$ tends to one, Einstein’s condensation occurs. Increasing $`\nu `$ corresponds to the increasing of the number of particles in the system. In all the states with indices $`n0`$ there is place only for a limited number of particles, while all the surplus, which can be arbitrarily large, gets located in the state $`|0`$. In this sense, when the number of particles becomes very large, we recover the model with the pure state discussed previously. A very interesting question is, whether experimentally it is possible to create condition, where the Einstein condensate would dominate.

Let us conclude this section with two remarks. For the Gaussian single particle density matrix

$$\rho _1^{(0)}(q,q^{})=\frac{1}{\sqrt{2\pi \mathrm{\Delta }^2}}\mathrm{exp}\left[\frac{q_+^2}{2\mathrm{\Delta }^2}\frac{R^2}{2}q_{}^2\right],$$
(42)   
where $`q_+=(q+q^{})/2`$ and $`q_{}=qq^{}`$, the eigenvalues and eigenfunction are known . Thus all the calculations can be easily performed. In fact they have been performed by various methods , , , .

The theory can be reformulated in the second quantization formalism. Then the function $`L(q,q^{})`$ appears as the Green function $`\widehat{a}_q^{}\widehat{a}_q^{}`$ and the possibility of expressing all the correlation functions in terms of $`L(q,q^{})`$ is the Wick theorem with $`L`$ as the only non zero contraction.

## 8 Statistical physics interpretation

Many results from the previous sections can be simply reinterpreted and rederived using standard statistical physics. Consider the single particle unsymmetrized density operator

$$\widehat{\rho }_1^0=\underset{n}{}|n\lambda _nn|$$
(43)   
with the condition $`_n\lambda _n=1`$. This can be reinterpreted as the density operator corresponding to the canonical ensemble, if we put

$$\lambda _n=\frac{1}{Z}e^{\beta \epsilon _n},$$
(44)   
where, as usual, $`\epsilon _n`$ is the energy of state $`|n`$, $`\beta `$ is the inverse temperature in energy units and $`Z=_n\mathrm{exp}[\beta \epsilon _n]`$ is the canonical partition function. The corresponding (single particle) Hamiltonian is

$$\widehat{H}=\underset{n}{}|n\epsilon _nn|.$$
(45)   
This Hamiltonian, when written in the coordinate representation, may look quite unusual, but some cases are simple. For instance, the Gaussian density matrix corresponds to the Hamiltonian of a harmonic oscillator.

For indistinguishable particles a single particle is not a convenient subsystem and, as suggested by Pauli long ago, it is better to choose as subsystem the open system consisting of all the particles in state $`|n`$. The state of this subsystem is defined by the number of particles $`(N)`$ in it. The probability of state $`N`$ of the subsystem is

$$P_n(N)=\frac{1}{𝒵_n}\overline{\nu }^Ne^{\beta N\epsilon _n}.$$
(46)   
The grand partition function

$$𝒵_N=\frac{1}{1\overline{\nu }e^{\beta \epsilon }}$$
(47)   
is chosen so that $`_{N=0}^{\mathrm{}}P_n(N)=1`$, the parameter $`\overline{\nu }`$ is known in statistical physics as the fugacity and is connected to the chemical potential $`\mu `$ by the formula

$$\overline{\nu }=e^{\beta \mu }.$$
(48)   
In order to reproduce the formulae from the previous sections, one puts $`Z\overline{\nu }=\nu `$. The grand partition function can be used to find the moments of the multiplicity distribution very much like the multiplicity generating function. For instance, for the average occupation of state $`n`$ we find

$$N_n=\frac{1}{\beta }\frac{Log𝒵_n}{\mu }=\frac{1}{e^{\beta (\epsilon _n\mu )}1}=\frac{\nu \lambda _n}{1\nu \lambda _n}.$$
(49)   
The probability of no particles in state $`n`$ is

$$P_n(0)=\frac{1}{𝒵_n}=1\overline{\nu }e^{\beta \epsilon _n}=\frac{1}{N_n+1}.$$
(50)   
The probability of no particle in the whole system is

$$P(0)=\underset{n}{}\frac{1}{N_n+1}.$$
(51)   
Let us consider two limiting cases. When all the occupation numbers are very small, the product equals approximately $`\mathrm{exp}[N]`$ and for large multiplicities it is very small. Very large fluctuations of the multiplicity are very unlikely to occur. When most particles are in the state $`n=0`$, the product is approximately $`(N+1)^1`$, which is much bigger than in the previous case. Thus, when there is much Einstein condensate, large multiplcity fluctuations become much more probable. Cosmic ray physicists have been reporting  observations of centauro and anticentauro events. This are high multiplcity events, where respectively either the neutral pions or the charged pions are missing. One could speculate that this phenomena are related to Einstein’s condensation.

Statistical physics gives also a simple interpretation for the function $`L(q,q^{})`$. One finds

$$L(q,q^{})=\underset{n}{}\psi _n(q)\psi _n^{}(q^{})\frac{1}{e^{\beta (\epsilon _n\mu )}1}.$$
(52)   
This is the canonical density matrix with the Maxwell-Boltzmann weights replaced by the Bose-Einstein weights. The fact that the Bose-Einstein weights fall with increasing energy $`\epsilon _n`$ faster than the Maxwell-Botzmann weights explains qualitatively most of the observations reported previously. For most Hamiltonians the wave function spreads in ordinary space and in momentum space, when energy is increased. Since the Bose-Einstein weights enhance the low energies, they reduce the average momenta and radii. Also the reduction of the effective radius of the interaction region, as determined from the width of the correlation function in momentum space, can be easily understood. If in the previous formula all the terms had equal weights, we would obtain $`L(q,q^{})=\delta ^3(qq^{})`$. The stronger the cut on the sum, the broader the peak in $`qq^{}`$ becomes. Since the Bose-Einstein weights are more peaked at low energies than the Maxwell-Boltzmann ones, they correspond to a broader peak in the correlation function. Since the width of this peak is inversely proportional to the radius of the production region, symmetrization reduces the radius of this region. All these qualitative arguments are usually true. It is, however, easy to show examples of hamiltonians, where e.g. with increasing energy the wave function shrinks either in ordinary space, or in momentum space. Additional assumptions necessary to convert these qualitative arguments into rigorous theorems are, therefore, necessary, but not yet known.