Heavy quarkonium in saturated environment of high-multiplicity pp collisions

High-multiplicity pp collisions exhibit features, traditionally associated with nuclear effects. Coherence motivates to treat high-multiplicity pp, pA and AA collisions on an equal footing. We rely on the phenomenological parametrization for mean multiplicities of light hadrons and J/ψ, assuming their linear dependence on Ncoll in pA collisions. The results of this approach underestimate the recently measured production rate of J/ψ at very high hadronic multiplicities. The linear dependence of J/ψ multiplicity on Ncoll is subject to predicted nonlinear corrections, related to mutual boosting of the saturation scales in colliding dense parton clouds. A parameter-free calculation of the non-linear corrections allows to explain data for pT -integrated yield of J/ψ at high hadronic multiplicities. Calculations are in a good accord with data binned in several pT -intervals as well. As was predicted, Υ and J/ψ are equally suppressed at forward rapidities in pA collisions. Consequently, their fractional multiplicities at forward rapidities in pp collisions are equal as well, and their magnitude agrees with data.


I. INTRODUCTION
The popular models of multiparticle production in pp and pA collisions are based on the eikonal multi-Pomeron exchange, or Glauber models. On the contrary to the wide spread believe, the Glauber model contains no clue to the multiplicity distribution. This is a model for the elastic scattering amplitude, which is related to the total cross section by the unitarity relation. One can also calculate the total inelastic cross section, however no relation between the multiple scattering terms in the elastic amplitude and inelastic processes follows from the Glauber model. The unitarity relations between each term in the Glauber elastic amplitude, expanded over multiple interactions, and the corresponding inelastic processes, was proposed by Abramovsky, Gribov and Kancheli [1], known as AGK cutting rules.
The unitarity cut of the elastic scattering amplitude can be done simultaneously through several Pomerons, while the other uncut Pomerons play role of absorptive corrections. Those corrections, also known as shadowing, lead to a substantial reduction of the total inelastic cross section. However, they cancel each other and do not affect the inclusive particle production cross section. This peculiar feature of absorption corrections to the inclusive cross section is known as AGK cancellation [1].
Relating different unitarity cuts of the elastic amplitude with multiplicity of produced particles (any species) one finds that multiplicity is proportional to the number of cut Pomerons, usually called number of collisions, N coll . Therefore, the normalised multiplicities of light hadrons and of J/ψ, defined respectively as, are expected to be equal, R J/Ψ = R h , and proportional to N coll . This expectation is in apparent contradiction with data [2][3][4][5], presented below in Fig. 1, demonstrating a significantly steeper rise of R J/Ψ with hadron multiplicity in comparison with R h . Moreover, it has been known since long time ago [6] that the simple relation R h = N coll , offered by the eikonal model, is strongly broken, and the multiplicity dependence on N coll for light hadrons produced in pA collisions, is usually parametrized as, with β h < 1 fitted to data. In the Glauber model N coll is given by The analysis of available data on pA collisions performed in [6] demonstrated consistency between different fits and no evidence for energy dependence of β h . The multiplicity measured at √ s = 5 TeV in [7], dN pA h /dy = 17.24 ± 0.66, which leads to β h ≈ 0.55. This is the value we rely upon in what follows.
The breakdown of the AGK cancellation is caused by coherence effects in gluon radiation [8]. The AGK rules [1] assume the Bethe-Heitler regime of particle production, which e.g. leads to n-times higher rapidity density of multiplicity dN h /dy for n cut Pomerons compared to arXiv:1910.09682v2 [hep-ph] 24 Oct 2019 a single cut. This regime is known to be broken by the Landau-Pomeranchuk effect of coherent radiation from multiple interactions. The reduction of gluon radiation caused by coherence is also known as gluon shadowing [9]. Only the accumulated mean transverse momentum of the radiated gluons keeps memory of multiple interactions, which lead to p T -broadening of radiated gluons, usually called saturation or color-glass condensate [10].
Besides broken N coll -dependence, Eq. (2), coherent production of quark-antiquark pairs also leads to a strong deviations from universality of R h . At first glance this might look surprising, because gluons are expected to hadronize at long distances, where no final state interaction is possible. Such a simplified interpretation of space time development is not correct.
As an example,qq production in DIS at small x exhibits shadowing due to propagation and attenuation of theqq dipole through the whole target [11]. This shadowing is a higher twist, so it is weaker for heavy than for light quarks.
Production ofqq by gluons is more involved, both the incoming gluon and producedqq interact with the target and screen each other. As a result, shadowing of quark production also is a higher twist and is described as propagation of an artificial 3-body dipole |qqg . Correspondingly, charm quarks [12], charmonia [13][14][15] are shadowed considerably less than light quarks. This explains in particular the observed non-universality The effects of coherence possess a strong scale dependence for the relation between fractional multiplicity and number of collisions. If the popular parametrization Eq. (2) is applied to J/ψ, one should expect a larger value of the parameter β J/ψ . Thus, from (2)-(4) we get, The difference between the parameters β J/ψ > β h , leads to a steeper rise of R J/Ψ in comparison with R h , as is demonstrated by data. The upper bound for the speed of rise for the ratio (5) is The processes of least shadowing corrections, like production of Drell-Yan pairs, Υ, B and D mesons, are expected to approach this bound. We do not attempt here at theoretical evaluation of β J/ψ (though is doable [15]), but prefer to rely on a fit to data. Another popular parametrization applied to measured A-dependence of J/ψ production is, where the fitted parameter α turns out to be close to unity, implying a weak nuclear suppression of J/ψ. Expanding this expression in small parameter 1 − α 1 we can relate β J/ψ in (4) and α, assuming a large number of collisions. The exponent α in (7) also does not expose a strong energy dependence. It was accurately measured at α = 0.95 in the fixed target experiment E866 at √ s = 40 GeV [16]. This value agrees with α = 0.9 − 0.98 measured at the mid-rapidity at √ s = 5.02 TeV [17,18], or interpolated between forward and backward rapidities [19][20][21].
Thus, according to (5) R J/Ψ rises linearly with increasing R h , or nearly linearly if to rely on Eq. (2), as was done in [22]. The dependence of R J/Ψ on R h shown in Fig. 1 by dashed curves was calculated at √ s = 5.02 TeV, with β h = 0.55 and the interval values of α = 0.95 − 0.98, which we consider as the corridor of the current uncertainty in measured A-dependence of J/ψ production. While these results agree with data at R h ≤ 4, available at the time of publication of [22], new measurements at √ s = 13 TeV and higher multiplicity up to R h = 7 [3][4][5] considerably exceed the predicted values of R J/Ψ , as one can see in Fig. 1. As we discussed above, neither β J/ψ , nor β h demonstate any significant energy dependence, therefore comparison of results in Fig. 1 at different energies is legitimate.

II. EFFECTS OF GLUON SATURATION
Multiple interactions naturally lead to broadening of the transverse momentum distribution, which has been calculated in [23,24] for proton-nucleus collisions in agreement with data, and predicted for high-multiplicity pp collisions in [22]. In the leading order broadening is proportional to the number of collisions with known coefficient, where N coll is related to R h by Eq. (4). The energy dependent factor C(E) in (9) controls the cross section σq q (r)of interaction of a small size r colorless dipole of energy E (in the nuclear rest frame) [23] C Notice that for broadening of gluons, which are responsible for J/ψ production, the pQCD factor 9/4 is introduced in (9). Since parton model interpretation of high-energy processes is not Lorentz invariant (only observable are), broadening of the p T -distribution, which looks like a result of multiple interactions in the nuclear rest frame, is interpreted as saturation of small-x parton density in the nucleus in its infinite momentum frame. Moreover the value of the saturation scale is directly related to the magnitude of the broadening, Q 2 s = ∆p 2 T [24]. Apparently, the additional kicks gained by the parton from multiple collisions, increase the effective scale of the process Q 2 ⇒ Q 2 + ∆p 2 T [13-15, 25, 26]. As a result of DGLAP evolution the nuclear gluon density g A (x, Q 2 + ∆p 2 T ) turns out to be suppressed at large x, but enhanced at small x [25].
Thus, the rate of J/ψ production in high-multiplicity pp collisions turns out to to be enhanced by gluon saturation and rises with R h steeper in comparison with Eq. (5).
In pA collisions only the projectile proton undergoes multiple interactions, which modify its parton distribution function (PDF), while the PDFs of bound nucleons remain unchanged. In the case of AA, or high multiplicity pp collisions, the interaction becomes symmetric, both assembles of colliding constituents are subject to multiple interactions, increasing their partonic content at small x 1,2 . On the other hand, multiple interactions in a denser gluonic medium become more intensive, leading to further increase of gluon density. Thus, a rising small-x gluon density in one of the colliding protons, induced by multiple interactions at large R h , stimulates stronger broadening of the partons in another colliding proton, which in turn leads to further increase of gluon density in the first proton. Such a mutual boosting of the gluon densities and saturation scales in the colliding protons at large R h is similar to the rise of the saturation scales in colliding nuclei is described in more detail in [25], and satisfies the bootstrap equation, where N coll is related to R h by (2). The characteristic scale Q 0 is the merging line between the nonperturbative and perturbative regimes of broadening. While the the latter is described perturbatively, the former is not calculable, but relies on the dipole phenomenology [24].
where C(E) is given by Eq. (10) and E = Q 2 0 /2m N x. Such a strategy is similar to what is used in DGLAP based analyses of data, except our evolution equation (11) is essentially nonlinear.
The saturated scale calculated with Eq. (11) for highmultiplicity pp collisions relative to the characteristic scale Q 2 = 4m 2 c vs R h is shown in Fig. 2. At high multi- plicities the boosted scale significantly exceeds the original scale, generating via DGLAP evolution more gluons at small x and a larger yield of J/ψ.
Notice that increase of the saturation scale hardly affects the hadron multiplicity R h , because hadrons with low p T , which give the main contribution to the multiplicity, are already produced with parton densities, which are at the unitarity bound. This is why β h show no evidence for energy dependence.

III. TRANSVERSE MOMENTUM DISTRIBUTION
An increase of the scale of the process at small x leads to a steep growth of the gluon density and as a result to enhanced production of J/ψ. Therefore the relative production rate R J/ψ plotted in Fig. 1 by dashed curves, should be corrected for the effects of saturation multiplying by the factor g N (x,Q 2 s + Q 2 )/g N (x, Q 2 ). The results are depicted in Fig. 1 by solid curves. Apparently, introduction of the effects of saturation improved agreement with data at high multiplicities.
It is also instructive to compare our results and data in different intervals of transverse momentum p T os J/ψ.
Here the multiplicity of J/ψ production in the nominator is taken at a given fractional multiplicity R h of light hadrons, while in the denominator the J/ψ production rate is summed over hadron multiplicity. To perform integrations in (13) one needs to know the explicit form of p T -dependence of the cross section. We rely on the popular parametrization, The mean transverse momentum squared can be extracted directly from data [20] applying a simple interpolation procedure. For √ s = 13 TeV we arrived at p 2 T = 11.72 GeV 2 . Fixing this value and the shape (14) we fitted the remaining parameter at n = 3.2.
We apply the same parametrization of the p Tdependence to the numerator of (13), but increase the mean value of p 2 T by broadening, p 2 T ⇒ p 2 T + ∆p 2 T , which was calculated with Eq. (9).
The data for R J/ψ measured in different p T intervals are depicted in Fig. 3 together with results of our calculations of Eq. (13) including the effects of saturation and mutual boosting of the saturation scale.

IV. J/ψ AND Υ AT FORWARD RAPIDITIES
Data show that J/ψ production rate in pA collisions is more suppressed at forward than at the mid-rapidity [16,19,21]. For the interval 2.5 < y < 4, measured in [19], the suppression factor A J/ψ ef f /A ≈ 0.65, which corresponds to α = 0.92.
Although the radius of Υ is smaller than of J/ψ, and Υ is less suppressed at the mid-rapidity, thebb dipole cross section rises with energy faster thancc and A Υ ef f is steeply falling towards forward rapidities (see Fig. 4 in [13]). At rapidities 2.5 < y < 4 higher-twist shadowing is predicted [13] to suppress Υ as much as J/ψ, A Υ ef f ≈ A J/ψ ef f = 0.65A. Notice that the parameters β J/ψ = β Υ = 0.57, are only slightly above the hadronic value β h = 0.55, so the originally expected relation R J/ψ = R h looks well satisfied. Of course, this does not mean restoration of the eikonal model, this is just a numerical coincidence.
Applying relation (5) we get R J/ψ = R Υ vs R h depicted by red solid curve in Fig. 4. Notice that it does not need to be corrected for the boosted saturation effects, because the correction is small at these rapidities (see Fig. 6 in [13]).

V. SUMMARY
High-multiplicity pp collisions exhibit features, which traditionally are associated with nuclear effects. Motivated by a long coherence time of interaction, we observe a close similarity between multiple interactions in pp, pA and AA collisions. In order to enhance multiple interactions in the former case one should trigger on high multiplicity of produced hadrons, while nuclei allow to reach the same multiplicity easier, by increasing the number of collisions Eq. (3). We employ the phenomenological description of the mean multiplicity in pA collisions, Eq. (2), violating the AGK cutting rules. Fitted to data, the observed nuclear effects for J/Ψ production, enable one to predict the multiplicity dependence of the J/Ψ production rate in pp collisions. However, the linear relation Eq. (5) between the fractional multiplicities R J/ψ and R h underestimate the yield of J/ψ at large hadron multiplicities [3][4][5].
Nevertheless, the linear dependence of the J/ψ production rate on N coll is subject to nonlinear corrections, predicted in [25]. They are related to mutual increase of the saturation scales in a colliding nuclei [25], or two dense parton clouds. The parameter-free calculation of the nonlinear corrections allows to explain the observed nonlinearity of R J/ψ vs R h , as well as data binned in several p T -intervals.
A sensitive test of the present model can be made with recently measured [3][4][5] yield of J/ψ at forward rapidities, which demonstrate similar fractional multiplicities R J/ψ = R h . We explain this result within our phenomenological description, while it might be a challenge for other available models. Besides, we predicted the same relation for heavy bottomia, R J/ψ = R Υ , which is also in a good accord with data [3][4][5].
Notice that the correlation of multiplicities of J/ψ vs hadrons observed at forward rapidities, are well described by Pythia [27]. However, this is not a parameter-free approach, the parameters of the generator are adjusted to data to be explained.