Bose enhancement, the Liouville effective action and the high multiplicity tail in p-A collisions

In the framework of dense-dilute CGC approach we study fluctuations in the multiplicity of produced particles in p-A collisions. We show that the leading effect that drives the fluctuations is the Bose enhancement of gluons in the proton wave function. We explicitly calculate the moment generating function that resums the effects of Bose enhancement. We show that it can be understood in terms of the Liouville effective action for the composite field which is identified with the fluctuating density, or saturation momentum of the proton. The resulting probability distribution turns out to be very close to the gamma-distribution. We also calculate the first correction to this distribution which is due to pairwise Hanbury Brown-Twiss correlations of produced gluons.


I. INTRODUCTION
Study of correlations in p-A and p-p collisions have been a very active area in the last years due to the observation of the ridge correlations at LHC. Since the ridge signal is much more pronounced in high multiplicity events, it is very important to understand the origin of the multiplicity fluctuation and especially the high multiplicity tail of the distribution, see Ref. [1][2][3][4] for the most recent experimental studies 1 .
In the present paper we address this question in the framework of the Color Glass Condensate (CGC) approach. We calculate the multiplicity momentum generating function, using McLerran-Venugopalan (MV) model for the projectile and assuming that the target is very dense. The latter assumption allows us to employ factorizable form for the averages of Wilson lines, as explained in Ref. [5]. We show that the main contribution to multiplicity fluctuations arises form the Bose enhancement (BE) of gluons in the projectile wave function (see Ref. [6] providing the interpretation of the correlations in the projectile in terms of BE; a related effect of Pauli blocking for quarks is discussed in Ref. [7]). This effect produces fluctuations which are not suppressed by the factor of the area of the projectile. We are able to resum the BE contributions exactly in the multiplicity generating function. The resulting distribution turns out to be γ-distribution, which for large moments practically coincides with the negative binomial distribution.
The other important effect in correlated gluon production is the Hanbury Brown-Twiss (HBT) effect. As discussed at length in recent literature [8], it is the leading cause for the angular correlations of produced gluons in the CGC approach. Its contribution to the total multiplicity on the other hand is suppressed relative to that of BE, as the correlated peak that it produces is very narrow. Nevertheless we identify the contributions to the multiplicity generating function due to the HBT within our calculational framework and calculate corrections induced by it.
We note that a calculation along similar lines has been undertaken some years ago in Ref. [9], see also Refs. [10][11][12] for numerical calculations and comparison to experimental data and Ref. [13] for higher order correction to an MV model. There are however significant differences between our approach and that of Ref. [9]. In particular Ref. [9] treated both the projectile and the target as dilute. It turns out that the large density effects of the target suppress half of the contributions considered as leading in Ref. [9]. Additionally the HBT contributions have not been included in the analysis of Ref. [9]. The resulting multiplicity distribution we obtain is somewhat different from that in Ref. [9].
The structure of this paper is the following. In Section II we lay out the general framework of the calculation and perform the averaging over the projectile wave function using the MV model.
In Section III we make a detour and consider specifically the second moment of the multiplicity distribution, a.k.a. inclusive two gluon cross section. We show explicitly that BE leads to the largest contribution to this moment and that this contribution is not suppressed by a power of area relative to the square of the single inclusive cross section. The argument is very similar to that in Ref. [9] except for the differences alluded to above due to the dense nature of the target. Similar observation was made also in Ref. [14].
In Section IV we calculate in closed form the momentum generating function which resums all BE contributions. We also relate it to the constrained effective potential approach proposed in Ref. [15] and show that our result for the momentum generating function is identical to the Liouville theory for the "composite field" discussed in Ref. [15].
In Section V we consider the HBT corrections to the distribution and provide a closed expression that resums the leading correction.
Finally Section VI contains a short discussion of our results.

II. THE GENERATING FUNCTION
Our calculations will be performed within the dense-dilute CGC framework. In this approach the number of produced gluons for a given configuration of the projectile (proton) and a target (nucleus) is given by [?
where the square of Lipatov vertex is Here ρ p is a given configuration of the color charged density in the projectile, and U is the eikonal scattering matrixthe adjoint Wilson line -for scattering of a single gluon on the target. The target Wilson lines depend on the target color sources, ρ t ; we suppress this in our notation. The single inclusive and double inclusive production in this approach are given by and where the averaging is performed over the projectile and target color charge configurations: and The normalization factors, Z p,t , are fixed so that In general m-particle production is Instead of computing the moments of inclusive particle number fluctuations we evaluate the moment generating function (see e.g. Ref. [16]) where we introduced an arbitrary k min ≫ Λ QCD . The moments of the distribution are obviously obtained from G(t) by differentiating with respect to t at t = 0. To calculate the generating function we have to specify the distribution of the sources in the projectile and in the target. For the projectile we will use the simple Gaussian McLerran-Venugopalan (MV) model specified by which corresponds to the weight functional Note that the structure of the ρ p correlator means translational invariance of the projectile wave function in the transverse plane. This assumption is only reasonable if we concentrate on momenta of produced particles larger than the inverse radius of the projectile. Thus in the following we will always assume k min > 1 R.
The averaging over the Wilson lines of the target will be performed in the approximation articulated recently in Ref. [8]. Any product of Wilson lines is factored into pairs with the basic Wick contraction Here the adjoint dipole amplitude is defined as As explained in Ref. [8] this approximation is appropriate for the dense regime. It collects all terms in the n-particle cross section which have the leading dependence on the area of the projectile. The approximation only includes terms which contain "small size" color singlets in the projectile propagating through the target. Any non singlet state that in the transverse plane is removed by more than 1 Q s from other propagating partons must have a vanishing S-matrix on the dense (black disk) target. On the other hand if the singlet state contains more than two partons, one looses a power of the area when integrating over the coordinates of the partons. Thus the leading contribution in the black disk limit is the one where only dipole contribution to the S-matrix should be accounted for. The same approximation for the quadrupole amplitude has been used previously in Ref. [14] where its consistency with explicit modeling of the Wilson line correlators via MV model has been verified. Note that this averaging procedure for the target is formally (disregarding subtleties related to the definition of the Haar measure) equivalent to the following form of the weight functional: A. Projectile averaging We now consider the calculation of the generating function We notice immediately that, since dN d 2 kdy is quadratic in both ρ p and U , and both the corresponding weight functionals are Gaussian as well, we can integrate over one of these "fields" exactly. We choose to integrate first over ρ p . The result of this integration is where the operator M is defined by its matrix elements At this point we need to make some approximations in order to perform the remaining functional integral over U . In the next section we will consider more closely the two gluon production cross section, which will help us understand the systematics of the leading contributions to the multiplicity fluctuations, and to devise the appropriate approximation that sums these leading contributions.

III. DOUBLE INCLUSIVE PRODUCTION: DISSECTING DIFFERENT CONTRIBUTIONS
To get some insight of which effects contribute the most to the multiplicity fluctuations we make a short detour and consider the double inclusive gluon production.

A. Dipole contribution
First consider the average with respect to projectile inside each factor dN d 2 kdy ρp,ρt , that is The projectile averaging gives The subsequent target averages give two distinct contributions.The first one involves target contractions inside each single inclusive factor and is disconnected. It reproduces the square of the single gluon production probability, and is of no interest to us. The connected contribution is The result is of order N 0 c . As we will see in the following the contractions of ρ that break the factors of dN d 2 kdy ρp,ρt are of order N 2 c and thus are more important. The contribution from Eq. (21) can therefore be neglected at large N c , and we will not try to include it and analogous contributions for higher moments in the generating function. The physics of this contribution was discussed in Ref. [17]. There it was shown that it corresponds to the Bose enhancement of gluons in the target wave function. Note that in the framework of the dilute target expansion utilized in Ref. [9] this contribution is leading order in N c . The different N c counting for the same quantity in the dense and dilute limits is not uncommon in saturation approaches.

B. The quadrupole contributions
We now concentrate on the other two contractions of the projectile color charges. There are two such contractions, and they both lead to a single trace "quadrupole" contribution the production probability: Each term has two contractions with respect to U of order N 2 c . We organize them according to their physical meaning [17].

The HBT term
The following contraction leads to the HBT contribution (cyclic property of trace was used) Substituting into double inclusive production we get

Bose enhancement in the projectile
The remaining contraction reflects Bose enhancement of gluons in the projectile wave function: Again substituting into double inclusive production

C. Conclusions on double inclusive
Comparing the HBT term Eq. (24) and the BE term Eq. (26), we see immediately that Eq. (26) gives the leading contribution to the multiplicity fluctuations.
The dominant contribution to the integral over q in Eq. (26) comes from small q due to the IR divergence of the Lipatov vertex Thus approximately we have For the MV model, µ 2 p = const we get a strong IR divergence of the integral ∫ d 2 q µ 4 p q 4 . This divergence is regularized by the momentum scale inversely proportional to the projectile size Λ = 1 Rp , that is modulo constant factors where S ⊥ is the area of the projectile. Thus the area dependence of the BE contribution to the double inclusive cross section is the same as that of the single inclusive cross section squared. This feature has been noted and discussed earlier in Ref. [9] and Ref. [14], see also Ref. [18] where the first saturation correction in the projectile to the double inclusive production was derived. The behavior of the HBT contribution is different. There is no quadratic IR divergence for the integration over q or p in Eq. (24). As long as k min ≫ Λ; no extra factor of area arises, and thus the HBT term is subleading. We note that if we include the soft scales in the k integral, that is take k min ∼ Λ the situation changes and the HBT effect becomes as important as the Bose enhancement. We will not consider this situation in the present work.
Although we have concentrated on the double inclusive cross section, it is easy to see that the analysis generalizes to the higher gluon production as well [9]. For production of n gluons, the BE term which has all contractions of the Wilson loops within the same single inclusive operator yields the highest power of area S n−1 ⊥ . Thus in the leading approximation we will only keep these terms. We will discuss the first corrections to this approximation later on.

A. Leading contribution from BE
We now return to Eq. (16). We do not know how to perform the integration over U in full generality. However in view of the conclusion of the previous section, we will first only keep the leading BE contributions. It is quite clear how to do that. The leading BE contribution corresponds to contracting the two Wilson lines within the same single inclusive gluon operator. In the context of Eq. (16) this corresponds simply to contractions between the two Wilson lines inside the trace of the logarithm in the exponent: where we have defined the integrated dipole operator: In Eq. (30) we have approximated the square of the Lipatov vertex by its leading term at low q, Γ(k, q, q) ≈ 1 q 2 , and have restricted the integration over q by k min from above. These restrictions can be in principle lifted, but as long as k min ≫ Λ Eq. (30) faithfully represents all leading area terms in any moment of the probability distribution.
For the MV model µ 2 p = const, the integral can be performed analytically, as follows ALEX: Here when computing the momentum integrals we assumed that IR divergences are regulated on the momentum scale of order Λ. Although this is reasonable in our dilute-dense approach, one may expect that higher order density correction may lead to modification of this IR momentum scale to the projectile saturation momentum Q p s , if Q p s > Λ.
The generating function for the cumulants is This is our result for the multiplicity moment generating function which resums all contributions due to the Bose enhancement of gluons in the projectile wave function. We now note some of its properties.

B. Some properties of the distribution
First of all, the mean number of particles being the first moment of the distribution, is given by the first derivative of G(t) at t = 0 and is proportional to the projectile µ 2 p and the logarithm of ratio of the momentum scales as expected The higher order cumulants are given by Also note that owing to presence of the increasing powers of 1 Λ 2 , κ n+1 ≫ κ n , and thus the factorial cumulants defined by are approximately equal to the cumulants, c n ≈ κ n .
The cumulants are very close to those of the γ-distribution Since the γ-distribution is known to exhibit the KNO scaling, we naturally expect to have KNO scaling at this order as well. We can check the KNO scaling by plotting the scaling function Eq. (37) for different values of k min since varying k min one varies the mean multiplicity. This is plotted in Fig. 1, right panel. Note that what is plotted here is the result of full numerical evaluation of the probability distribution (see Appendix for detail), and not just the leading approximation. The numerics is performed by using the MV model also for the target fields with the scale µ t ≫ µ p . One observes that indeed for large enough k min the quality of the KNO scaling is very good. There is one subtlety here. Our distribution Eq. (32) is not exactly the γ-distribution. In particular the first moment Eq. (34) has an additional logarithmic factor relative to the parameter that determines the higher moments Eq. (35). Since this logarithmic factor depends on k min , it could potentially affect the KNO scaling. Nevertheless this additional logarithmic dependence is very slowly varying, and the scaling is clearly seen in the numerical results in Fig. 1. The left panel in Fig. 1 illustrates what happens when the soft scale Λ is raised and becomes comparable with the target saturation momentum. One clearly observes that as Λ grows, the probability distribution becomes narrower. Note that in this regime we cannot neglect the effects of HBT as well as correction due to exact form of the Lipatov vertex. These are the effects that drive the narrowing of the distribution for larger Λ.
The leading order probability distribution we obtained is similar in many respects to that obtained in Ref. [9]. There are however some significant differences. First, in the case of dense target the contribution of the Bose enhancement in the target wave function is suppressed by the factor 1 N 2 c , whereas in Ref. [9] the target was treated as dilute, this contribution was of order unity and contributed to the probability distribution on par with the projectile BE. Second, the cumulants in our case are very close to those given by the gamma distribution, whereas Ref. [9] finds negative binomial distribution (NBD). The n-th cumulants of the two distributions differ by the factor n − 1. This difference can be traced to a different way the function µ 2 p (p) is treated at small p in the two approaches. Our calculation corresponds to taking constant µ 2 p and cutting off the putative infrared divergence in the integrals by a finite area of the projectile. The IR regulator therefore does not arise from making the correlation between the sources nonlocal in coordinate space, but rather imposing an impact parameter profile on the source density. On the other hand the authors of Ref. [9] regulate IR divergences by taking µ 2 p (p) ∼ p 2 for momenta smaller than the projectile saturation momentum. Physically this corresponds to taking the correlation function between the two sources to be nonlocal in coordinate space.
This different treatment of the IR behavior leads to different integrals with respect to the incoming gluon from the projectile wave function q. In this paper we keep µ p = const and thus while the authors of Ref. [9] approximate µ 2 p ≈ q 2 ; in this case the integral in Eq. (39) brings no additional factors dependent on n.
Related to this point is also the different energy dependence that we expect from the distribution. In our case the parameter of the distribution is α = This parameter does not depend on energy for the dilute projectile. Thus we expect α to be approximately energy independent. On the other hand the treatment of the IR in Ref. [9] leads to replacement of Λ 2 by the saturation momentum of the projectile Q 2 s . This grows with energy, and thus the probability distribution has rather strong dependence on energy.

C. Constraint action formalism and the Liouville potential
One interesting property of our derivation is that the probability distribution of produced particles that we find is deeply related with the probability distribution of particles in the projectile wave function. Note that we have approximated the Lipatov vertex by its part that involves only the momentum of the gluon coming from the incoming wave function. Thus if in the rest of the calculation we simply set U (x) = 1, that would correspond to calculating the particle distribution in the projectile wave function. This last move would simply correspond to setting k min = 0 in Eq. (31), since ∫ d 2 k (2π) 2 D(k) = D(x = 0) = 1 which is equivalent to setting U = 1. This then gives Clearly we should not set k min = 0 in the integration limits in Eq. (30), but instead keep it fixed. Thus we conclude that taking the limit of Eq. (40) in Eq. (30) corresponds to calculating the probability distribution for particles with transverse momenta smaller than momentum k min in the projectile wave function.
Interestingly, such a calculation can be performed independently using an alternative formulation based on the framework of the constraint effective potential, see Refs. [15,19,20].
The main idea of this approach is to integrate out fluctuations of ρ p (q) that do not affect a specific operator defined through ρ p (q). In Ref. [15] the constraint effective potential for the gluon distribution defined by the covariant field, A + , was derived where A + (q) = g q 2 ρ p (q), and This potential corresponds to a Liouville potential with negative Ricci scalar 2 . The generating function for gluon multiplicity in the projectile is defined as Using the effective potential (43) we have For large S ⊥ this integral can be computed using the saddle point approximation to yield which reproduces the result obtained previously. Interestingly, the origin of the logarithm in this equation is owing to the presence of the Liouville logarithm in the effective constraint action for η(q), see Eq. (43). Note that the form of the integral appearing in Eq. (45) is quite suggestive, Here t can be viewed as an external field, while the moments and the cumulants of the gluon number can be viewed as the moments and the cumulants for fluctuations of the "composite field" ∫ d 2 q η(q) q 2 . Thus the fluctuations of the number of particles provide a direct measurement of the Liouville potential. It is interesting to note that our derivation provides a concrete realization of early ideas in the literature about relevance of Liouville action to multiplicity fluctuations in CGC. Ref. [22] postulated ad hoc such a Liouville potential for saturation momentum fluctuations. In the present paper we instead derive it form the constrained effective potential for the MV model 3 .
Although the argument of our potential is not the saturation momentum, but rather the composite filed η, they are closely related. Consider the effective potential for η close to its saddle point value at zero external field. Expanding in ln η we obtain Recall that the operator definition of the composite field η is given by Eq. (41) which is interpreted naturally as the scaled fluctuating saturation momentum of the projectile, since on a configuration by configuration basis the saturation momentum is indeed determined by the square of the color charge density. Thus we indeed can interpret Eq. (49) as the probability distribution for the fluctuating saturation momentum, which has exactly the same form as assumed in Refs. [23,24].

V. CORRECTIONS TO THE GENERATING FUNCTION
In this section we consider corrections to the cumulant generating function. There are two basic types of corrections: those attributed to BE with lower power of IR divergence when integrating with respect to the incoming gluon from the projectile wave function and those attributed to the HBT contribution, see the section for two particle gluon production.

A. Subleading BE terms
The subleading IR terms due to Bose enhancement can be easily resummed. We simply have to allow for the full Lipatov vertex and for the q-dependence of the dipole amplitude in Eq. (30). Thus formally we get whereD The explicit form of the distribution now depends on the dipole amplitude D(k), and can be calculated once this amplitude is known.

B. The HBT contributions
We now concentrate on the corrections of the second type. We will keep only terms leading in N c .
We rewrite our basic expression for the generating function as The BE contribution is represented by the first line of Eq. (54) is discussed in the detail in previous section. The last two terms are the corrections we are after. We have introduced the normal ordered product of the Wilson lines The normal ordering insures that apart from the BE term, no other terms contain contractions between two Wilson lines that belong to the same "vertex". These contractions have been completely resummed into the propagator of the color charge density We can expand the functional integral into series in the "interaction term".
where the interaction part is given by To compute the cumulant generating function it is useful to introduce the Feynman rules depicted in Fig. 2. Note that due to normal ordering in the vertex, no diagram with contraction of the two U 's belonging to the same vertex are allowed. Those have already been resummed in the propagator of ρ p .
To understand which type of diagrams give the leading corrections, let us first consider a particular example: the connected contribution involving four vertexes (correction to inclusive four particle production). The leading N c contributions have two different topologies, see Fig. 3.
The diagram a) in Fig. 3 is given by (modulo combinatorial and kinematic factors) The diagram b) in Fig. 3 is given by As before we consider the IR dominant contribution from the Lipatov vertices. For the diagram b) we get Γ(k, q 1 , q 2 )Γ(k, q 2 , q 3 )Γ(k, q 3 , q 4 )Γ(k, q 4 , q 1 ) ≈ q 1 ⋅ q 2 which after angular integration has logarithmic divergence for each integral of the form ∫ dq 2 i q 2 i . For the diagram a) we get and in this case the integral with respect to q has quadratic divergence in IR. This divergence is of course regulated by Λ 2 ∝ 1 S ⊥ which leads to the extra power of S ⊥ . It is now clear what are the leading diagrams due to the HBT corrections that contribute most to the generating function. Those are the diagrams that contain maximal number of the ρ p propagators at the same momentum q, since each such propagator is accompanied by a product of two Lipatov vertexes containing the same momentum q, thereby leading to one extra power of area for each additional ρ p propagator. These diagrams are of the type of Fig. 3 a, where every the vertexes are organized into pairs with two vertexes of the pair connected by two propagators of the U field, and one propagator of ρ p . Physically this corresponds to contributions to the n-gluon inclusive production, where gluon pairs are emitted independently, but the HBT correlations are present between two gluons in each pair.
We are now going to ignore diagrams of type B but resum all diagrams of type A, corresponding to pairwise HBT correlations.
For a diagram with 2n vertexes we have the following combinatorial coefficient 1 (2n)! from expanding the exponential; (2n − 1)!! -the number of ways to organize the 2n vertexes into pairs, which then will be contracted over U ; 2 n -the number of contractions of U -two possibilities within each pair of vertexes; 2 n -the number of contraction between two ρ p 's within the pair of vertexes. Although there are 4 possible of contractions in general, once a particular contraction of U 's is chosen, only two contractions of ρ p 's are leading order in N c ; (n−1)! 2 ways of ordering the n vertex pairs along a circle. This is not just a number of permutations n!, since the position of the first vertex should be fixed (periodicity along the circle)-that given (n − 1)!, and the factor 1 2 is due to identical nature of a permutation and its reflection, since in both cases every vertex has identical neighbors, and that's all that matters; 2 n way to contract the factors of ρ p along the circle, since for each vertex each one of two ρ p 's can either be contracted with its right or left neighbor; the color factors all cancel out except for the overall N 2 c − 1 in front for any n.
Thus at the end of the day the diagram with 2n vertexes which are contracted pairwise into the "daisy" is 2 2n−1 1 n . Resumming these terms we obtain and Finally we obtain the cumulant generating function in the form which has the same form as Eq. (32); thus the momentum integral can be performed analytically leading to Eq. (33) with the substitution D → (D + tZ). This is a rather simple expression, and one can analyze the effects of the correction given a model for the dipole amplitude D(p). These corrections may be important in the regime where k min is not significantly greater than the soft scale Λ. Although we do not consider such a situation in the present paper, at high enough energy the projectile wave function itself will acquire a saturation momentum scale Q s,p significantly larger than Λ. In this case it is quite conceivable that in our expressions the soft scale Λ will be replaced by this semi soft scale Q s,p . In this case it is perfectly sensible to consider k min < Q s,p . The HBT contributions in this regime will become significant, and the relative significance of the BE and HBT contributions to the multiplicity fluctuations has to be reanalyzed. We will not attempt to do it in the present paper.

VI. CONCLUSIONS
In this paper we studied the multiplicity fluctuations in p-A collisions within the framework of dense-dilute CGC formalism using the MV model for the wave function of the proton. Our approach is similar in many aspects to that of Ref. [9]. There are however some significant differences, and our results are quite different as well. As opposed to Ref. [9] we treat the target as very dense, while on the projectile side we do not assume any dynamical "correlation length" associated with the saturation momentum. We rather treat the IR physics of the proton wave function as genuinely non perturbative, governed by a soft scale of order of the inverse proton size. As a result we obtain the probability distribution which within large range of energies is energy independent.
We identified two sources of multiplicity fluctuations: those due to the Bose enhancement of gluons in the proton wave function and the HBT effect in the initial stages of scattering. Interestingly, in the dense-dilute framework the Bose enhancement in the nucleus wave function leads only to a (N 2 c − 1) −1 suppressed contribution to any cumulant of particle number, and is thus a subleading effect. We demonstrated that as long as the low momentum structure of the proton wave function is dominated by the genuine soft scale (the "proton size"), the dominant origin for the multiplicity fluctuations is the Bose enhancement.
We have calculated explicitly the moment generating function for the multiplicity distribution due to BE. The distribution we obtain is very close to the γ-distribution. Just like the γ-distribution it satisfies the KNO scaling with very good precision. Interestingly, the leading term in the generating function for multiplicity of produced particles is practically identical to the generating function for multiplicity distribution in the projectile wave function. This latter quantity can be calculated using the effective action approach as suggested in Ref. [15]. We have shown that this effective action is nothing but the Liouville action for the composite field, which can be thought of as fluctuating density (or saturation momentum).
The authors of Ref. [9] obtained the negative binomial distribution for the multiplicity. Our result, as mentioned above is somewhat different, although for large moments the NBD and the γ-distribution are quite similar. The main difference, as explained in the text is in the physics of the scale that regulates the formal IR divergences. In the dense -dilute calculation performed in the present paper this role is played by the soft scale of the proton radius, and not by the semi soft scale of the projectile saturation momentum. If one assumes that the proton wave function itself is characterized by a finite correlation length, much smaller than the proton size one would have to reanalyze to what extent the dominance of the BE persists. It may well happen that the HBT contributions become equally important and have to be included in the leading order calculation. In fact this is precisely what happens on the target side, where the presence of large saturation momentum strongly suppresses the BE effect, as we noted above. It is thus not clear to us that the approximation of BE dominance and finite saturation momentum of the projectile are mutually compatible.
We note that both NBD, and γ-distribution have rather long tails for large values of produced multiplicity. It is very natural that these tails are associated with the quantum Bose enhancement effect of identical gluons, just like the Bose-Einstein distribution of identical noninteracting bosons. Thus we believe that although the details of the distribution are model dependent (MV model in our case), the main feature of large fluctuations is universal as long as the fluctuations are dominated by BE.
Finally, we have also calculated the correction due to the generating function due to pairwise HBT correlations. In the regime studied in the present paper this correction is small. However it is bound to become important in the regime of saturated projectile, and therefore in itself would be an interesting object of study. and compute the fundamental Wilson lines Further details can be found in Refs. [26,27].