Comparison of the Color Glass Condensate to multiplicity dependence of the average transverse momentum in p+p, p+Pb and Pb+Pb collisions at the LHC

The first momentof the charged-particle transverse momentum spectrum and its correlation with the charged-particle multiplicity N_{ch} provide vital information about the underlying particle production mechanism. The ALICE collaboration recently reported thatversus N_{ch} in Pb+Pb collisions is smaller than in p+p and p+Pb collisions. Other interesting features of data is rather flatness ofat high N_{ch} in Pb+Pb and p+Pb collisions in seemingly striking contrast to the case of p+p collisions. With a detailed calculation, we show all these peculiar features in a wide range of energies and system sizes can be well described by the idea of gluon saturation within the Color Glass Condensate framework using the k_T-factorization. This establishes an important fact that the bulk of the produced particles in heavy-ion collisions at the LHC carries signature of the initial stage of collisions. We also show that the recent scaling property seen by the CMS collaboration between the number of tracks in p+p and p+Pb collisions may provide a strong evidence in favor of geometric-scaling phenomenon and gluon saturation, indicating that the underlying dynamics of high multiplicity events in p+p and p+Pb collisions should be similar.


I. INTRODUCTION
The recent LHC measurements of particle correlations in azimuthal and pseudorapidity in proton-lead (p+Pb) collisions [1][2][3], the so-called Ridge, and its similarity to the same effect in lead-lead (Pb+Pb) collisions, have raised the question whether the same underlying dynamics is responsible in both cases. There are two very distinct pictures which equally provide a good description of this phenomenon, the Color Glass Condensate (CGC) approach based on the initial-state (and gluon saturation) effects in the nucleon or nuclear wavefunctions [4] (see also Refs. [5][6][7]) and the hydrodynamical approach based on the final-state rescattering effects [8][9][10]. In order to distinguish between different scenarios of particle production mechanisms in p+Pb collisions, it is important to investigate on the equal-footing, the multiplicity dependence of particle production, correlations, event shapes in p+p, p+Pb, and Pb+Pb collisions. Such studies have been already undertaken by the LHC experiments.
Recently the ALICE [11] and the CMS [12] collaborations reported the measurements of the first moment p T of the charged-particle transverse momentum spectrum and its correlation with the charged-particle multiplicity N ch in p+p, p+Pb, and Pb+Pb collisions. While the rise of p T with N ch in p+p collisions can be understood within a model with final-state color reconnection between strings produced in multiple parton interactions [13,14], the same model cannot describe p+Pb and Pb+Pb data. On the other hand, while the EPOS model [15] which assumed collective flow, describes p+Pb data at high-multiplicity, it over-predicts the Pb+Pb data and gives opposite trend versus N ch . Other Monte Carlo event generators such as DPMJET [16], HIJING [17] and AMPT [18] also fail to describe these data, for model comparison see Ref. [11]. Therefore, an unified description of p T versus N ch from p+p, p+Pb to Pb+Pb collisions seem to be a challenge. The main aim of this letter is to address the importance of the initial-state effect, a missing important ingredient in all above-mentioned models. Here, we provide an unified description of all these data, within the CGC approach [19][20][21] based on the gluon saturation [19] and Glasma physics [22,23]. We also show that the scaling property between the number of tracks in p+p and p+Pb collisions, recently observed by the CMS collaboration [12], indicates the geometric-scaling phenomenon [24,25] as expected in high-energy QCD.

II. MAIN FORMULATION
In the CGC approach [19][20][21], the gluon jet production in A+B collisions can be described by k T -factorization given by [26], where C F = (N 2 c − 1)/2N c with N c being the number of colors. We define x 1,2 = (q T / √ s)e ±y where q T , y are the transverse-momentum and rapidity of the produced gluon jet, and √ s is the collision energy per nucleon. The unintegrated gluon density φ G A (x i ; k T ) gives the probability to find a gluon that carries x i fraction of energy with k T transverse momentum in the projectile A (or target B). The unintegrated gluon density is related to the imaginary part of the forward quark anti-quark scattering amplitude on a proton (or nucleus) target N p(A) (x i ; r T ; b) via, with a notation where r T denotes the dipole transverse size and b T is the impact parameter of the scattering. The most important ingredient of the single inclusive hadron production cross-section in Eq. (1) which captures the saturation dynamics is the fundamental or adjoint dipole amplitude. The impact-parameter dependence of the amplitude is crucial here. We use the impact-parameter Color Glass Condensate (b-CGC) dipole saturation model [27]. In the b-CGC dipole model, the color dipole-proton amplitude is given by, with effective anomalous dimension defined as where Y = ln(1/x) and κ = χ ′′ (γ s )/χ ′ (γ s ), with χ being the LO BFKL characteristic function. The parameters A and B in Eq. (4) are determined uniquely from the matching of the dipole amplitude and its logarithmic derivatives at rQ s = 2. The dipole amplitude depends on the saturation scale Q s , namely the momentum scale at which non-linear gluon recombination effects start to become as important as the gluon radiation [19]. In the b-CGC dipole model, the saturation scale of proton explicitly depends on the impact-parameter and is given by The parameters of the b-CGC dipole amplitude (N 0 , B CGC , γ s , x 0 , λ) are determined via a fit to HERA data [27]. This model approximately incorporates all known properties of small-x regime of QCD [28] including the impactparameter dependence of the scattering amplitude [29] and it describes all existing small-x data at HERA both in the inclusive and exclusive diffractive processes including the recently released high precision combined HERA data [27]. The dipole-nuclear amplitude is constructed in the standard fashion from the dipole-proton amplitude and proton saturation scale by employing the nuclear thickness function [30,31]. For a review of k T -factorization phenomenology within different saturation models see Refs. [32,33]. In order to take account of the difference between rapidity y and the measured pseudo-rapidity η, we employ the Jacobian transformation h(p T , η, m 2 jet ) between y and η where the parameter m jet is mini-jet mass which mimics the pre-hadronization effect 1 [30]. For the value of strong-coupling α s in Eqs. (1,2) we employ the running coupling prescription used in Ref. [30,31,34]. The average transverse momentum of the mini-jet is calculated from Eq. (1), The average transverse momentum of the mini-jet can be directly related to the saturation scale, see below.
The first question one has to address here is whether via the single-inclusive gluon production, one can generate high multiplicity N ch >> N ch in p+p collisions. One can show that with centrality cut based on the single-inclusive k T -factorization supplemented with the proton-dipole amplitude, one can maximally generate upto N ch ≈ 1 ÷ 2 N ch for very central p+p collisions at √ s = 7 TeV and η = 0. Although the main underlying mechanism of highmultiplicity event production in p+p and p+A collisions is still unknown, multi-gluon production and fluctuations in geometry and color charge distributions should play significant role. In the CGC framework, the yield of multiple gluons production can be obtained from Glasma color flux tubes decay [23] which leads to the negative binomial distribution for the multiplicity, in accordance with experimental observation at RHIC and the LHC [35]. Using the multiple gluons production yield from Glasma flux tubes decay [23], one can immediately show that in leading log approximation, deep in the gluon saturation regime, the average transverse momentum of a produced gluon in multiple-gluon production is equal to computing the same quantity in the single-inclusive production, where we used the notation . . . to indicate the averaging in multiple-gluon yield while the notation . . . is used for a averaging over a single-inclusive yield. In the Glasma picture, the transverse area S ⊥ is filled with Q 2 s S ⊥ independent flux tubes of size 1/Q 2 s and each of these tubes emits gluons with the typical momentum of order the saturation scale Q s . It may then seem natural to expect that Eq. (8) to be held deep inside the saturation region. Therefore, at high-multiplicity events at the LHC, assuming that we are in the saturation region, the average transverse momentum of a gluon-jet can be directly obtained via the single-inclusive k T -factorization defined in Eq. (7), even though, the multiple-gluon production is the main source behind the production of such rare events. This may also be considered here as our working hypothesis.
In the CGC picture, the gluon saturation scale is proportional to the parton density, and since the parton density is proportional to the particle multiplicity, then the saturation scale appeared in Eq. (8) should depend on the multiplicity, its own first cause. Let us define two saturation scales: one for the projectile Q A s and another for the target Q B s , we also introduce two auxiliary variables Q s,min = min.
Using the k T -factorization Eq. (1), one can show that deep in the saturation region where q T < Q s,min we approximately have, while in the kinematic region where only one of target or projectile is in the saturation region Q s,max > q T > Q s,min , within the same approximation we have [36], where ζ = Q s,min /Q s,max . Note that the above relations are by construction built in the KLN model [37]. Here, in accordance with the underlying assumption in Eq. (8), we consistently assume that we are deep inside of saturation regime and take the relation given in Eq. (9) to relate the saturation scale to charged-particle multiplicity N ch by fixing the prefactor in Eq. (9) to the average minimum-bias charged-particle multiplicity N ch . Therefore, we replace where the integral in b is now in minimum-bias 2 . We recall that large multiplicity events with N ch >> N ch can be described by the multiple gluon production yield obtained from Glasma flux decay yield, see the discussion after Eq. (7). On the other hand, the average transverse momentum of the produced single inclusive gluon in a highmultiplicity event can be still obtained via the k T -factorization Eq. (1) due to the equality given in Eq. (8). Note that the saturation scale given in Eq. (11) by construction reproduces the charged-particle multiplicity N ch via the k T -factorization Eq. (9) upto some logarithmic corrections (consistent with the same approximation made in Eq. (8)). In this sense, deep inside of the saturation regime, one can approximately mimic the effects of Glasma flux decay (and fluctuations) by redefining the saturation scale and employ Eq. (1) to compute the average transverse momentum of the produced gluon in such rare high-multiplicity events.
It is instructive to note that from Eqs. (9,10), for symmetric collisions like p+p and A+A collisions at mid-rapidity with Q s,min = Q s,max , we approximately have while in asymmetric cases like p+A collisions where one of the sources is significantly weaker than the other one with | ln ζ| >> 1 − ζ, the multiplicity dependence of q T becomes weaker than the former cases and the simple relation in Eq. (12) is not more reliable [36]. This is consistent with the fact that based on only dimensionality consideration, the relation in Eq. (12) is approximately expected to be held in the saturation region. However, it is important to note that at a fixed multiplicity the interaction area in p+p, p+A and A+A collisions are very different and in practice there are large logarithmic corrections to Eq. (12). Our main aim here is to calculate q T via Eq. (7) without resorting to any approximation. Finally, we should also specify our hadronization scheme. It is generally an open problem how to incorporate the fragmentation processes into the CGC/saturation formalism. In order to clearly disentangle the initial from final-state effect, we employ here a simple scheme for the final-state hadronization, based on the so-called Local Parton-Hadron Duality (LHPD) principle [38]. Namely we assume that the hadronization is a soft (or semi-soft) process and cannot change the direction of the emitted radiation. Furthermore, we assume that the hadronization occurs in vacuum and it is the same for p+p, p+Pb, and Pb+Pb collisions. Note that the main contribution of the k T -factorization for the multiplicity distribution comes from small transverse momentum p T < 1 GeV. This is in accordance with the experimental observation that the average transverse momentum of the produced charged-particle is rather small p T < 1 GeV for a wide range of kinematics and system size, see Figs. 1,2. Note that the fragmentation functions based on fixed order approximation to DGLAP evolution is dubious at low transverse momentum 3 and we do not employ them here since we are only interested in low p T region. In the framework of the LHPD, the p T spectrum of the produced hadron upto a normalization, is obtained from Eq. (1) by replacing q T = p T / z where p T is the transverse momentum of the produced hadron and the parameter z is the average fraction of energy of the mini-jet carried by the hadron. The average transverse momentum of the produced hadrons can be then obtained by where p intrinsic T is the average intrinsic transverse momentum of the hadron in the mini-jet and q mini-jet T is computed via Eq. (7).

III. NUMERICAL RESULTS AND DISCUSSIONS
In our formalism, we have only four unknown parameters, the overall normalization factor in Eq. (1), the mini-jet mass m jet , the average fragmentation parameter z and the average intrinsic transverse momentum of the hadron p intrinsic T . Note that for spectra of hadron, the free parameter p intrinsic T does not enter into the calculation while for the average transverse momentum the normalization will drop out. Therefore, for all observables considered here we practically have only 3 free parameters. All these unknown parameters are fixed via a fit to minimum-bias p+p data at low-energy [30,39]. Therefore, our results shown here at higher energies and different system sizes from p+p to p+A and A+A collisions should be considered as a free-parameter calculation. The over-all normalization, and the mini-jet mass are fixed via a fit to the experimental data on the charged-particle multiplicity at mid-rapidity. Unfortunately, we do not know how mini-jet mass changes with the system size and kinematics. Here, we focus on kinematics around mid-rapidity and assume that the mini-jet mass is independent of kinematics [39]. The parameters z and p intrinsic T are fixed via a fit to a spectra and average transverse momentum of charged-particle in minimum-bias p+p collisions at low energies (not shown here). We therefore obtain z ≈ 0.5 ÷ 0.6, and p intrinsic In Fig. 1 left panel, we show our results for the spectra of charged-hadron production in minimum bias p+p collisions at √ s = 2.36 and 7 TeV. The fact that the shape of spectra of produced hadron at low p T resembles the spectra of mini-jet spectra upto a normalization provides the first hint toward importance of the initial-state effect. Note that the curve at 7 TeV in Fig. 1 (left panel) was predicted in Ref. [30], see also Ref. [43]. In order to make prediction for the multiplicity dependence of the average charged-hadron transverse momentum, we should also know the value of the average multiplicity in the minimum bias event selection N ch at various energies for different interactions. We impose the same kinematic constrains as the ones employed in the ALICE measurements [11], namely we restrict the integrals in Eq. (7) to 0.15 < p T < 10 GeV and |η| < 0.3. In the case of p+Pb collisions, similar to the experiments at the LHC, we also take into account the rapidity shift ∆y = −0.465. As a first internal check, we reproduced the experimental values of N ch reported by the ALICE collaboration given in table 2 of Ref.
[11]. In Fig. 1 right panel, we compare our results with the ALICE recent data on the average transverse momentum of charged particles as a function of multiplicity in p+p collisions at various energies. Note that ALICE data [11] is consistent with the corresponding CMS [12] and ATLAS [44] measurements. The increase of p T with energy and multiplicity is naturally expected in the gluon saturation picture. This is simply because p T increases with the saturation scale, the only dynamical scale available in the system, and the saturation scale grows with energy and density, see Eq. (12).
In Fig. 2 left panel, we show average transverse momentum of charged-particle as a function of charged-particle multiplicity N ch in p+p, p+Pb and Pb+Pb collisions at √ s = 7, 5.02 and 2.76 TeV respectively, in the range 0.15 < p T < 10 GeV and |η| < 0.3. In Fig. 2 right panel similar to left panel, we show the corresponding average transverse momentum of the mini-jet in p+p, p+Pb and Pb+Pb collisions. Comparing both panels in Fig. 2, it is seen that remarkably the trend of the average transverse momentum of charged particles as a function of N ch at various system size and energies resembles the same quantity for the mini-jet gluon. The average transverse momentum of charged particles p T is smaller in Pb+Pb than p+Pb and p+p collisions at high-N ch . The main reason is that the effective interaction area is different in Pb+Pb compared to p+p and p+Pb collisions. Note that events with N ch < N ch , is more prepherial and less dense compared to the minimum bias collisions. We recall that the average with a factor roughly about A 1/6 (where A is atomic number). Therefore, we expect p intrinsic T for the case of lead target to be A 1/6 bigger than the same quantity for the proton, namely p intrinsic  charged-particle multiplicity measured by the ALICE collaboration [11] is N ch ≈ 259.9, 11.9 and 4.42 in Pb+Pb, p+Pb and p+p collisions respectively at the kinematics considered in Fig. 2. Therefore, a event with a given value of N ch in Fig. 2 corresponds to different density in p+p, p+Pb and Pb+Pb collisions. In other words, at moderate N ch in range of N ch < 150 shown in Fig. 2 we are in dilute region in Pb+Pb collisions given that N ch ≈ 259.9 while the same multiplicity event selection corresponds to a very rare high-density event in p+p collisions. Note that the TeV as a function of tracks charged-particle multiplicity for |η| < 2.4 obtained from the IP-Glasma model [48]. ALICE data in Pb+Pb collisions shown in Figs. 2,3 correspond to peripheral collisions with N part < 52 [45][46][47]. We recall that at 50 − 60% centrality in |η| < 0.5 we have dN ch /dη = 149 ± 6 and N part = 52.8 ± 2.0 [45].
The rise and flatness of p T in different collisions with different system sizes can be also simply understood from Eq. (12). First, we recall that in the saturation region, the interaction area is subject to rapid rise with multiplicity and then it becomes independent of multiplicity and flattens at a critical N cri which is typically larger than N ch [48]. Now, given the fact that we have at low N ch < N ch , increase in multiplicity and the interaction area, approximately cancel out each others in Eq. (12) leading to the flatness of the average transverse momentum. In the case of p+p collisions with N ch >> N ch , shown in Fig. 2, the interaction area is roughly constant and then the average transverse momentum increases with multiplicity as expected from Eq. (12), in accordance with experimental data. Similar to Fig. 2, in Fig. 3 we show the average transverse momentum of charged particles as a function of charged-particle multiplicity N ch in a wider range of N ch , upto N ch = 1400 for Pb+Pb collision at the LHC. In order to clearly see the slow rise of p T with N ch , we plotted the results in logarithmic scale. The bands labeled b-CGC is consistent with the curves shown in Fig. 2 and include the theoretical uncertainties associated to fixing the mini-jet mass m jet in Eq. (7), and the parameters z and p intrinsic T in Eq. (13). The CMS collaboration [12] has recently observed a very intruging scaling behaviour in average transverse momentum of various identified hadrons as a function of the true track multiplicity N tracks in p+p and p+Pb collisions for |η| < 2.4. Namely it was found that the p+Pb curves (in left panel Fig. 5) is approximately similar to p+p curves by taking the p+p values and multiplying their N tracks coordinate by a factor of 1.8, for all particle types, see Fig. 5. In other words, p+Pb collisions with a given N tracks is similar to a p+p collisions with 0.55 × N track . It is important to notice that this scaling is the same for all different identified hadrons indicating that perhaps its origin may be traced back to the initial stage of collisions before hadronization. The observed scaling in number of tracks is in fact expected in the CGC approach. We recall that in the CGC framework, two multiplicity events with the same saturation scale should lead to the same physics. Now, using Eq. (9), one can immeaditely relate track multiplicity in p+p and p+Pb collisions with a similar saturation scale: where N p+p tracks and N p+P b tracks denote track multiplicity in p+p and p+Pb collisions respectively. S p and S A are the effective transverse interaction areas for a given centrality or multiplicity in p+p and p+A collisions, respectively. The pre-factor K takes into account the effect of different center-of-mass energy in p+p and p+Pb collisions at the LHC. The saturation scale of proton is Q 2 s ∝ (x 0 /x) λ ∝ s λ/2 , where the parameter λ ∼ 0.22 was extracted from the LHC data in p+p collisions [34]. Therefore we have K ∼ ( √ s = 7 TeV/ √ s = 5.02 TeV) λ = 1.076. For calculating the interaction area as a function of tracks (or multiplicity), we use the recent results of the IP-Glasma initial-state model of hadrons and nuclei [48]. The IP-Glasma model [48,49], is a CGC-type approach to particle production at early stage of collisions based on the classical Yang-Mills description of initial Glasma fields which properly incorporates the impact parameter dependence via the IP-Sat dipole model for proton [50]. In the IP-Glasma [48], one can from the first principle compute the radius of interaction in term of gluon multiplicity. In order to relate the gluon multiplicity to the number of corrected tracks, following Ref. [25] we employ where for the CMS kinematics of interest with |η| < 2.4, we take ∆η = 4.8. In Fig. 4, we show the ratio of interaction areas in p+p √ s = 7 TeV and p+Pb √ s = 5.02 TeV collisions as a function of number of tracks. Note that the scaling-factor in Eq. (16) (shown in Fig. 4 up to the pre-factor K) changes slowly with N tracks with a median about 0.55, and at very large N tracks it is constant since at very high-multiplicity the interaction areas in both p+p and p+A collisions do not change further. In Fig. 5 left panel, we show first average transverse momentum of identified charged hadrons (pions, kaons, protons) in the range |y| < 1, for all particle types, as a function of the track multiplicity for |η| < 2.4 in p+p and p+Pb collisions. Using S p /S A shown in Fig. 4, we can relate the CMS number of tracks N tracks (and ALICE charged-particle multiplicity N ch ) in p+p and p+A collisions via Eq. (16). In Fig. 5, we rescaled track multiplicity in p+Pb collisions via Eq. (16), and compare the average transverse momentum of various identified charged hadrons in p+p and p+Pb collisions. In Fig. 5 (lower panel), we demonstrate that the ALICE data shown in Fig. 2, follows the same scaling property seen in the CMS data, within error bars. In Fig. 5 (lower panel), we also compare our results obtained from the k T -factorization having rescaled the results in p+Pb collisions. It is interesting to see that the rescaled full calculation results agree very well with the rescaled data. This indicates that the small discrepancy seen in Fig. 5 between p+p and rescaled p+Pb data, is mainly due to approximation employed in Eq. (16). Nevertheless, it is remarkable that simple geometric-scaling relation in Eq. (16), without introducing any new parameters or ingredients, is in accordance with data for identified charged hadrons and charged particles. Note that the key ingredient in Eq. (16) is the saturation scale Q s,min (the proton saturation scale) which relates charged-particle multiplicity in p+p and p+Pb collisions. The interaction areas and the over-all normalization factor in Eq. (16) are then computed within the saturation model.
There are, however, a number of caveats which need further study before taking the number predicted here for the scaling-factor at face value. First, note that in the above, for simplicity we have ignored the impact-parameter b dependence of the saturation scale and possible correlation between x and b. Therefore, the relation given in Eq. (16) is less reliable for prepherial collisions at low N tracks where a proper treatment of impact-parameter dependence of the collisions is indispensable. Moreover, in Eqs. (16,17), we approximately related the multiplicity to number of track with a constant factor and assumed that this relation to be the same for both p+p and p+A collisions. One should also bear in mind that the ratio of interaction areas computed here in the IP-Glasma model [48,49], is based on the IP-Sat model [50]. It will be of great interest to compute the ratio of interaction areas in the Glasma model based on the b-CGC saturation model [27].
Some words of caution are in order here. Strictly speaking our formalism is less reliable for p+p (and A+A) collisions at around mid-rapidities. This is due to the fact that the k T -factorization is valid for asymmetric dilute-dense collisions [23,26] like p+A (with a very large nuclei) or p+p (and A+A) collisions at forward rapidities.
Note that the existence of the geometric-scaling in p+Pb collisions was also discussed in Ref. [25]. The ALICE collaboration later confirmed this with rather large error bars [11]. However, the above-mentioned caveats are also presents in previous studies. Therefore further investigations on this line are needed in order to pin down the true nature of the observed scaling behaviour. It would be of great interest to see if final-state type approaches like hydrodynamic can also explain this scaling phenomenon. In principle, the scaling property in number of tracks between p+p and p+Pb collisions given in Eq. (16), should be also correct for multiplicity dependence of other observables in p+p and p+Pb collisions at high-energy (and low p T ) in the saturation region. This can be verified at the LHC.
To summerize, in this letter within the CGC framework we provided an unified description of the recent ALICE (and CMS) data on p T at various energies and system sizes from p+p to p+Pb and Pb+Pb collisions. In our approach neither final-state hadronization nor collective hydrodynamic-type effects are important to describe the main features of the current data. This clearly indicates that the underlying dynamics of particle production at small-x region is universal and the main behaviour of the bulk of the produced particles in heavy-ion collisions at the LHC can be simply described by the idea of gluon saturation. In Fig. 3, we show our predictions for p T in a wider range of N ch . In order to discriminate among difference models and pin down the importance of final-state effects like hydrodynamic, it is crucial to have experimental date at large N ch . We also showed that the observed scaling property seen in the CMS and the ALICE data between the number of tracks in p+p and p+Pb collisions may provide a strong evidence of the gluon saturation and the so-called geometric-scaling, indicating that the underlying dynamics of high multiplicity events in p+p and p+Pb should be similar.