Particle multiplicities in the central region of high-energy collisions from $k_T$-factorization with running coupling corrections

Horowitz and Kovchegov have derived a $k_T$-factorization formula for particle production at small $x$ which includes running coupling corrections. We perform a first numerical analysis to confront the theory with data on the energy and centrality dependence of particle multiplicities at midrapidity in high-energy p+A (and A+A) collisions. Moreover, we point out a strikingly different dependence of the multiplicity per participant on $N_\text{part}$ in p+Pb vs.\ Pb+Pb collisions at LHC energies, and argue that the observed behavior follows rather naturally from the convolution of the gluon distributions of an asymmetric vs. symmetric projectile and target.

The Color Glass Condensate (CGC) approach to particle production in high-energy collisions conjectures that the energy and system size dependence of the p T -integrated multiplicity can be computed in weak coupling. The qualitative argument for this conjecture is that the running coupling at the particle production vertex would be effectively evaluated at a scale of order the semi-hard "saturation scale" Q s Λ QCD , even at low p T . McLerran and Venugopalan have shown that such a semi-hard scale indeed emerges for a large nucleus due to the high density of valence color charge per unit transverse area [1]. Furthermore, the running coupling Balitsky-Kovchegov (rcBK) equation [2] for the unintegrated gluon distribution (UGD) shows that the saturation scale grows with energy. Most gluons "in the wave function" of a hadron or nucleus have transverse momentum k T ∼ Q s , suppressing the sensitivity to the infrared, k T ∼ Λ QCD [3]; c.f. Fig. 1 below.
There have been many studies of the energy and centrality dependence of particle production in p+A and A+A collisions within the k T -factorization approach [4], using UGDs which exhibit "saturation" at k T < Q s [5][6][7][8][9]. Whenever running of the coupling has been considered, an ad-hoc choice for the scale of α s (Q) in the k T -factorization formula had to be made 1 . For example, ref [9] assumed that the coupling is evaluated at max | p T ± k T |/2 so as to avoid the infrared regime (thanks to k T ∼ Q s , as mentioned above). While in practice the sensitivity to such ad-hoc running coupling prescriptions may not be very large, it is clearly worthwhile to assess running coupling corrections from a more solid theoretical basis.
Horowitz and Kovchegov have derived a k T -factorization formula beyond LO to include running coupling corrections [11] (also see ref. [12]) to single-inclusive (small-x) gluon production in the scattering of two valence quarks. Their expression results from a resummation of the relevant one-loop corrections into the running of the coupling. They propose the following generalization to hadron-hadron or hadron-nucleus collisions: (1) (Our notation follows ref. [11]; k now denotes the transverse momentum of the produced gluon while q and k − q are the "intrinsic" transverse momenta from the gluon distributions.) This distribution of gluons in transverse momentum and rapidity has to be convoluted with a fragmentation function in order to obtain the p T -distribution of produced hadrons. Eq. (1) implicitly assumes that collinear factorization applies in fragmentation. Λ 2 coll is a collinear infrared cutoff which should match the scale of the fragmentation function typically chosen as µ 2 FF p 2 T . We have computed hadron transverse momentum distributions in p+A collisions in this way and shall report our findings elsewhere. Here, we are primarily concerned with the p T -integrated multiplicity where the most relevant regime is that around the average p T . For this regime we employ a simple model fragmentation function D(z, µ 2 FF ) ∼ δ(1 − z). For the observables considered here slight modifications of this fragmentation function mainly affect the normalization in Eq. (1) and can be absorbed into N . The normalization also absorbs "K-factors" due to higher order corrections and will be fixed by matching to data.
The unintegrated gluon distribution is given by Note that these functions do not involve a factor of 1/α s (k 2 ); instead, the factors of the inverse coupling with the appropriate scale appear explicitly in Eq. (1). N A (r, y, b) denotes the forward (adjoint) dipole scattering amplitude at impact parameter b. We assume a uniform gluon density within a proton 2 . N A (r) approaches a constant as r → ∞ and so φ(k) ∼ k 2 vanishes at low transverse momentum. For k 2 Q 2 s , on the other hand, φ(k) ∼ 1/k 2 . In the absence of the non-linear corrections to small-x evolution present in the BK equation this behavior would extend down to low k T . The numerical form of N A employed here is identical to that used previously in ref. [9]; it has been obtained in ref. [13] by solving the rcBK equation 3 . Here, we restrict to using their solution for MV model initial condition since p T -integrated observables for their other UGD sets are not much different at small x. Finally, we stress that the dipole forward scattering amplitude obtained in ref. [9] has been averaged over all BK gluon emissions without bias. A plot of φ(k) is shown in Fig. 1. This corresponds to the unintegrated gluon distribution after 3 units of rcBK rapidity evolution.
The factors of the inverse coupling in Eq. (1) are determined by the coefficient of the one-loop β-function (with N c = N f = 3), and we set Λ MS = 0.24 GeV. The Q 2 -dependence of the coupling is given by is given by the complex conjugate of this expression so that the product ln Q 2 is real, as it should be.
In the limit q k this simplifies to 4 so that Therefore, at high transverse momentum, and choosing the collinear cutoff scale Λ 2 coll = k 2 , the spectrum of produced gluons is proportional to α s (k 2 ) k −4 ln 2 k 2 .
For k q we have Here, the dominant contribution to Eq. (1) is from q ∼ Q s since φ(q) quickly decays when q is far from Q s . Eq. (6) shows that the distribution of produced gluons is well defined at low transverse momentum, k → Λ MS . The spectrum can be integrated over k > Λ MS without encountering a divergence. (It is not sensible to address the spectrum of "gluons" with k < Λ MS .) Physically, dσ/d 2 kdy ∼ α s (Λ 2 coll ∼ k 2 )/k 2 should rather level off when both collision partners are dense; k Tfactorization fails here. For proton-nucleus collisions with Q s,A Q s,p the actual contribution to dN/dy from k < Q s,p is small and we can therefore simply apply Eq. (1) down to k = Λ MS . For nucleus-nucleus collisions this is not justified since there are many particles at Λ MS < k < Q s . Here, the correct spectrum below Q s can only be obtained from a "dense-dense" computation which does not rely on k T -factorization. (To date, such calculations, for example, refs. [10,15]), have been performed only for fixed coupling, or with ad-hoc running.) On the other hand, phenomenological applications of k T -factorization have been rather successful in reproducing the dependence of the multiplicity in A+A collisions on energy and centrality [5,8,9]. This is presumably due to the fact that for large nuclei and high energies this dependence is entirely determined by the single scale Q s . In any case, given those previous applications of k T -factorization with ad-hoc scale choice to the centrality dependence of the multiplicity in A+A collisions it is certainly interesting to also see the result obtained from Eq. (1). Hence, we shall simply cut off Eq. (1) at k min of order Λ MS . We have checked the dependence of dN/dη in Pb+Pb collisions at LHC energies on N part for k min = Λ MS . . . 2Λ MS and obtained virtually identical curves.
To compute produced particle multiplicities in p+A and A+A collisions we convolute Eq. (1) with a Monte-Carlo Glauber simulation, which has been described in more detail in refs. [9]. This allows us to compute the dependence of the multiplicity on the number of participants. First results using Eq. (1) were obtained in ref. [16] for minimum bias collisions with KLN model [5,6] gluon distributions. Fig. 2 shows our results for the multiplicity per participant pair in A+A collisions at RHIC and LHC energies. We have fixed the normalization factor in Eq. (1) to match to central Pb+Pb collisions at 2.76 TeV; the same normalization has been used for all other energies, centralities, and collision systems. The data shown in Fig. 2 is from refs [17][18][19][20][21][22][23]. Our curves are very close to those published in ref. [9] using ad-hoc scale setting. The data at √ s = 200 GeV mainly probes the MV model gluon distribution at the initial x 0 = 0.01 rather than small-x rcBK evolution. The multiplicity as a function of N part shows the well known increase of dN/dη per participant towards more central collisions. It is driven by the increasing overlap in transverse coordinate space of the 2d projections of the nuclear Woods-Saxon distributions. This leads to increasingly symmetric collision partners at any given point in the transverse plane so that the convolution integral of the gluon distributions in Eq. (1) increases as both transverse momentum arguments can be near the "saturation peak".
We also show our prediction for Xe+Xe collisions at 5.44 TeV in Fig. 2. We have updated the figure to include new data for Xe+Xe collisions released by the ALICE collaboration [24].   [25][26][27]. The energy dependence of the multiplicity obtained from the r.c. k T -factorization formula with rcBK UGDs compares well to the measurements at LHC energies. The extrapolation to RHIC energy of √ s = 200 GeV, however, is significantly too low. This is not unexpected since at such energies one is sensitive mainly to the MV model initial condition imposed at x 0 = 0.01 rather than to small-x evolution. This should in fact fail, in particular for small systems, since the MV model assumes a large nucleus. Improving results for p/d+A collisions at RHIC energy (and η ∼ 0) will require an improved theoretical understanding of the unintegrated gluon distribution of a proton at x = 0.01 and greater, as well as possibly additional corrections to Eq. (1).
The dependence of the multiplicity at LHC energies on N part is rather interesting. Somewhat surprisingly perhaps, we find that beyond N part 4 the multiplicity per participant decreases slightly with N part . This is due to the fact that for increasingly asymmetric collisions the convolution in transverse momentum space of the gluon distributions does not increase in proportion to N part . Simple considerations suggest that it grows logarithmically (also see the discussion in refs. [6]). A numerical fit to the 5.02 TeV curve shown in Fig. 3(right), for 15 ≥ N part ≥ 5, gives ∼ ln 1.25 (N part )/N part . In contrast, A+A collisions become more symmetric as the impact parameter decreases and the multiplicity per participant increases with N part .
Such a feature is also seen in data, as shown in Fig. 3 (right), where we show ALICE [28] and ATLAS [29] data for (1/N part ) dN ch /dη vs. N part in p+Pb collisions at 5 TeV. In fact, in ref. [26] the ALICE collaboration already noted that the multiplicity per participant in NSD p+Pb collisions at 5 TeV (averaged over N part !) is 16% lower than in NSD p+p collisions interpolated to the same collision energy. Remarkably, this trend appears to continue beyond N mb part 8. While the distribution of multiplicity is a quantity that can be measured in a fairly direct way, the number of participants is not, and in principle depends on the method of centrality selection. We refer to the above-mentioned publications for more detailed discussions of the experimental centrality selections and their determination of N part , but here show the ALICE "ZNA N Pb-side part " and the ATLAS "GGCF ω σ = 0.2" results, which are believed to be less model dependent and may be the most suitable to compare to N part as used in our model. These data exhibit a trend similar to the calculation, but with somewhat flatter dependence on N part . This could be due to the lack of a realistic impact parameter dependence of the proton-UGD in our computations, and due to a bias on the gluon distribution introduced by the experimental centrality selection. A more accurate matching of (1/N part ) dN ch /dη to the measurements would entail accounting for such bias on the configurations of small-x gluon fields through reweighting [31].  (1) and has instead been normalized to the CMS measurement [32] of dE T /dη in minimum bias p+Pb collisions at 5 TeV. The dependence on N part and on energy is then a prediction. As before, the number of participants should be determined in the fragmentation region of the nucleus with a method that smoothly approaches the p+p limit as N part → 2. This omits the bias on E T due to fluctuations of the multiplicity at midrapidity which p+A collisions inherit from p+p [33].
Our computation predicts an increase of E T per particle for increasingly asymmetric collisions. Such "broadening" of the transverse momentum distributions of produced gluons is expected due to the increase of the saturation scale of the target nucleus [5]. We should point out that dE T /dη is, however, sensitive to final state interactions which may reduce its magnitude [34].
One may attempt to interpret the decrease of the particle multiplicity per participant with N part noted above in a simple two-component "soft + hard" model [35]. However we can see that this is not possible.
Let f ( √ s) ≥ 0 denote the fractional contribution from the hard component which is proportional to the number of binary collisions. 1 − f ( √ s) then corresponds to the soft contribution, proportional to N part : In ref. [35] Kharzeev and Nardi obtained f ( √ s = 56 GeV) = 0.22 and f ( √ s = 130 GeV) = 0.37. For √ s = 5 TeV we find that a universal fit from p+p to central Pb+Pb with Eq. (7) is impossible. Fitting to very peripheral Pb+Pb collisions only (N part ≤ 34 corresponding to ≤ 17 participants per nucleus, on average) we estimate f ≈ 0.26 ± 0.01. A similar fit of the new Xe+Xe data by ALICE [24], again for N part ≤ 34, gives f ≈ 0.34 ± 0.01. We consider this a lower bound on the value of f appropriate for p+A collisions since leading-twist perturbative processes may already experience slight "shadowing" even in rather peripheral heavy-ion collisions.
For p+A collisions we can rearrange the above equation as follows:  Fig. 3(right). The formula describes the data below the average N part 8 for minimum bias collisions fairly well. However, the trend for more central p+Pb collisions appears different from the data shown. Furthermore, since N coll is linear in N part for p+A collisions, this simple model would not predict an increase of the transverse energy per particle like in Fig. 4.
Let us summarize the main points of this paper. We have performed the first analysis of the energy and centrality dependence of particle multiplicities in the central region of high energy p+A collisions predicted by k T -factorization with running coupling corrections, and rcBK gluon distributions. We point out that the formula derived by Horowitz and Kovchegov [11] results in a well-defined gluon transverse momentum distribution 5 down to p T ∼ Λ MS . Our numerical results show that the r.c. k T -factorization formula with rcBK gluon densities provides a good description of the energy and centrality dependence of the multiplicity in p+A collisions at LHC energies, √ s > 1 TeV. For p+A collisions at RHIC energies, on the other hand, a better understanding of the unintegrated gluon distribution of a proton at x ∼ 0.01 is required. It may be worth pointing out that we have attempted to introduce as few model parameters as reasonably possible in order to exhibit where the current theory fails. In our analysis of particle multiplicities we fitted a single energy, centrality, and system independent constant: the normalization factor N in Eq. (1).
Our main observation relevant for phenomenology is to note that the convolution of the unintegrated gluon densities of a proton and of a nucleus increases more slowly than linear with the number of participants. The asymmetry of the gluon distributions in p+A collisions results from the coherence of the interaction with the dense target. As a consequence, the multiplicity per participant in increasingly asymmetric p+A collisions is found to decrease slowly. (This may flatten out somewhat if a bias on the small-x gluon distribution is taken into account.) Such behavior is markedly different from that for more and more central, and increasingly symmetric A+A collisions, as well as from expectations based on a simple two-component "soft+hard" particle production model with only energy dependent shares.