Multiplicity dependence of quarkonia production in the CGC approach

In the paper we develop, in the CGC approach, a new mechanism for quarkonia production in events with high multiplicity of the soft produced particles. In this mechanism both J /Ψ and soft hadrons are produced at short distances r ∝ 1 / Q s , where Q s is the saturation scale, which is considered to be proportional to the multiplicity of soft hadrons. Using this approach, we analyze the vigorously growing multiplicity dependence of quarkonia production, which was recently discovered by the STAR and ALICE col-laborations. We demonstrate that this result presents a strong evidence in favor of multigluon fusion mechanisms of the quarkonia states. We analyze in detail the contribution of the 3-gluon fusion mechanism and demonstrate that it describes correctly all the available experimental data on the multiplicity dependence of J /ψ meson production. We also make predictions for other quarkonia states, such as ψ( 2 S ) and Υ ( 1 S ) , and ﬁnd that the multiplicity dependence of these states should be similar to that of the J /ψ meson. Finally, we discuss an experimental setup in which the very strong multiplicity dependence could be observed. This observation would be a strong evidence in favor of the CGC approach.


Introduction
The production mechanisms of hadrons which contain heavy quarks, especially charmonia states, remains one of the longstanding puzzles of high energy particle physics. The conventional paradigm considers that the dominant contribution comes from gluon-gluon fusion, and that the heavy quarks formed in the process emit soft gluons in order to hadronize into quarkonia states [1][2][3][4][5][6]. This picture gives reasonable estimates for the total and differential cross-sections, within a e-mail: Marat.Siddikov@usm.cl (corresponding author) uncertainty of Long Distance Matrix Elements (LDMEs) of quarkonia states [7][8][9][10][11][12][13][14]. Nevertheless, the above described picture was challenged by recent measurements of the multiplicity dependence of certain heavy mesons. As was found recently by the STAR [15,16] and ALICE [17,18] collaborations, the J/ψ yields grow rapidly as function of the multiplicity of co-produced charged particles. Since these charged hadrons are predominantly soft, the description of their production presents in general a complicated theoretical problem in the realm of non-perturbative QCD. For this reason up to now there have been few attempts to describe theoretically the multiplicity dependence of quarkonia production [19,20]. This problem has also been addressed by using Monte Carlo event generators [21][22][23][24], which combine perturbative QCD with a phenomenological treatment of the hadronization processes. However, recent data [18] clearly demonstrated that none of the above-mentioned approaches could predict correctly the observed yields of J/ψ mesons at large multiplicities. As we will show below, the large multiplicity data present also a challenge for all the approaches based on the gluon-gluon fusion mechanism.
Recently, in [25], we suggested that the experimentally observed multiplicity dependence of quarkonia production [18,26] could indicate a sizable contribution of the multipomeron mechanisms. 1 This claim was inspired by the observation that each pomeron on the average gives the same contribution to the observed multiplicity, and for this reason the events with large multiplicities potentially might get enhanced contributions from configurations with additional pomerons. While it is expected that such higher-order contri- butions should be suppressed in the heavy quark mass limit, for the 3-pomeron term such suppression does not work due to quantum numbers of the S-wave quarkonia [25,27,28]. Indeed, as we can see from Fig. 1, formally the two-pomeron and three-pomeron mechanisms contribute at the same order in the strong coupling α s . For a long time it was believed that the two-pomeron mechanism dominates, since the softness of the emitted gluons in the two-pomeron mechanism implies that the suppression in the heavy quark mass limit is milder than expected, from O(α s )-counting. However, the contribution of the three-pomeron mechanism is also enhanced at high energies, due to the additional t-channel pomeron and the associated increase of gluon density. As was illustrated in [25,27,28], the latter mechanism provides a reasonable description of the rapidity and the p T -dependence of produced charmonia, 2 as well as describes qualitatively the main features of the multiplicity dependence. The goal of this paper is to analyze in detail the multiplicity dependence in the framework of the CGC/saturation approach, compare its predictions with available experimental data and make predictions for the forthcoming measurements.
As we have mentioned, the relation between the production of quarkonia, which comes from short distances, and of soft hadrons, which originates from long distances, is a very complicated problem for theoretical approaches due to the present embryonic state of non-perturbative QCD. One possible solution is the building of Monte Carlo programs, which are based on perturbative QCD and on the phenomenological experience in the treatment of the hadronization processes. In this paper we suggest an alternative approximation, based on the CGC approach, which being an effective QCD theory at high energies, states that a new phase of QCD: a dense sys-tem of partons (gluons and quarks) is produced in collisions with a new characteristic scale: the saturation momentum Q s (W ), which increases as a function of energy W [37][38][39][40][41]. At LHC energies, the scale Q s (W ) is much larger than the scale of soft physics μ soft , and therefore the dynamics of the events at those energies can be treated theoretically in the framework of the CGC approach. Nevertheless, the transition from this system of partons to the measured state of hadrons is still an unsolved problem, and for this reason at the moment we need to use a pure phenomenological input for the long distance non-perturbative physics. Hence, our model for multi-particle production at high energies consists of a two stage scenario. The first stage is the production of a gluon with transverse momentum p T ∼ Q s (W ), which is much larger than the scale of soft physics μ soft , and therefore can be treated theoretically in the framework of the CGC approach. The second stage is the hadronization of the produced hard gluons. In what follows we will rely on the widely used hypothesis of Local Parton-Hadron Duality (LPHD) [42][43][44], which assumes that the distributions of observed hadrons differ from the distributions of partonic mini-jets only by a numerical constant. In addition, we will assume that the correlations between hadrons in rapidity are long range in nature and stem from the production of parton showers, described by BFKL Pomerons. The experimental data show that the short range correlations are small [45][46][47], in accord with the expectation from CGC models for high energy scattering [48].
The paper is structured as follows. In the next Sect. 2 we describe a framework for quarkonia production in the CGC/Saturation approach. In Sect. 3 we outline our main theoretical results and make numerical estimates of the multiplicity dependence, compare with available experimental data, and make predictions for the ψ(2S) and Υ (1S) quarkonia. Finally, in Sect. 4 we draw conclusions and discuss some open challenges in our approach.

Three-pomeron contribution to charmonia production
The contribution of the three-pomeron mechanism to the cross-section of the S-wave quarkonia production is given by [25] dσ y, √ s dy where y is the rapidity of the produced quarkonia, measured in the center-of-mass frame of the colliding protons; Ψ M (r, z) is the light-cone wave function of the quarkonium M (M = J/ψ, ψ(2S), Υ (1S)) with transverse separation between quarks r ; and z is the light-cone fraction carried by the quark (see Appendix B for its parametrization). In the heavy quark limit the wave function Ψ g of the heavy quark-antiquark pair might be approximated by the perturbative expressions from Ref. [49] (see Appendix B for details). We also use the notation is the dipole scattering amplitude. The factors in the first line of (1) are explained in detail in [25] and are related to the gluon uPDFs and DPDFs in the proton. The nonperturbative factor ∼ d 2 Q S 2 h (Q T ) is important for the study of the absolute cross-sections considered in [25], although eventually it will cancel in the selfnormalized observables considered in this paper. The second factor includes the gluon PDF x 1 G (x 1 , μ F ) , which we take at the scale μ F ≈ 2 m c . In LHC kinematics and at central rapidities (our principal interest) this scale significantly exceeds the saturation scale Q s (x), which justifies the use of the two-gluon approximation. However, in the kinematics of small-x g (large energies) there are sizable non-perturbative (nonlinear) corrections to evolution in the CGC/Saturation approach. In this kinematics the corresponding scale μ F should be taken at the saturation momentum Q s , while the gluon PDF x 1 G (x 1 , μ F ) is in this approach closely related to the dipole scattering amplitude N (y, r ) = d 2 b N (y, r, b) as [50,51] where y = ln(1/x). Equation (3) might be inverted and gives the gluon uPDF in terms of the dipole amplitude, This result allows us to rewrite (1) in a more symmetric form, entirely in terms of the dipole amplitude N , which will be used for analysis below. While the relations (4,5) are valid only when both xG and φ satisfy the linear evolution equations, Eq. (6) is the correct generalization of Eq. (1) for the BK non-linear evolution (see Refs. [33][34][35][36]). Therefore, Eq. (6) correctly accounts for the diagrams of Fig. 2b, which is the goal of this paper.

Dipole scattering amplitude in CGC approach
In the CGC/saturation approach the dipole amplitude N (y, r, b) should satisfy the non-linear Balitsky-Kovchegov (BK) equation [52][53][54] for dipoles of small size r . In the saturation region this solution exhibits a geometric scaling, being a function of one variable τ = r 2 Q 2 s , where Q s is the saturation scale [55][56][57][58].
The main observation of Ref. [25] is that the typical r and r in Eq. (6) are in the vicinity of the saturation scale r 2 Q 2 s ∼ 1 (r 2 Q 2 s ∼ 1), where we know the general solution of the non-linear BK equation [59]. It has the form: At long distances the scattering amplitude N approaches 1, having the following form [58]: where κ is a constant (see Refs. [50,58]) Hence the difference of the amplitude in Eq. 6 at large r 2 Q 2 s 1 vanishes as where z ± = ln r ± r 2 Q 2 s /4 . Therefore, one can see that at large τ δN 2 vanishes.
However, to find the behavior of δ N 2 at moderate τ we need to solve the non-linear Balitsky-Kovchegov [52][53][54] equation in the entire kinematic region. Generally speaking, the BFKL kernel (χ(γ ) in γ -representation) takes into account all twists contributions. Nevertheless, for our estimates we restrict ourselves to the leading twist term only, which has the form [58] (see also Ref. [48]): As indicated in Eq. (10) we have two types of logs: denotes the saturation scale. To sum these logs it is necessary to modify the BFKL kernel in different ways in the two kinematic regions, as shown in Eq. (10).
In the saturation region r 2 Q 2 s > 1 it has been shown [55][56][57][58] that the scattering amplitude manifests geometric scaling behavior and the BK equation with the simplified kernel of Eq. (10) can be re-written in the form [58]: where N (z) = z dz N z and z = ln τ 2 with τ = r Q s .
In Ref. [58] it has been shown that the solution of this equation takes the form for the function φ which has been defined as The value of the constant C has to be found from matching with the region τ < 1. For study of the small φ ∼ φ 0 we can see that we should choose C = 0. Indeed, in this case the solution at small φ has a form which coincides with the general solution of Eq. (7) for the region τ < 1 at small φ 0 = N 0 1. In Fig. 2 we plot the solution of Eq. (12). One can see that the function which enters the master equation (see Eq. (6)), 3 decreases with growing τ and the largest contribution stems from the region τ ∼ τ . Therefore, the estimates from Eq. (11) supports our main idea that the typical distances in Eq. (6) are in the vicinity of the saturation scale where both τ ∼ 1 and τ ∼ 1. The solid curve in Fig. 2 gives ΔN 2 at r = r which corresponds to the main contribution to Eq. (6). One can see that at τ = 1÷2 this curve shows τγ behaviour. The red curve presents the behaviour of the overlapping integral which is a function of r = τ/Q s but not of τ . In the figure we plot this function fixing Q 2 s = 9 GeV 2 which corresponds to large n ≈ 10 ÷ 12. One can see that the main contribution stems from the region of τ , where we can use N ∝ τγ .

Saturation model
There are different phenomenological parametrizations for N available in the literature. One of such parametrizations, which we will use for our numerical estimates, is the CGC parametrization [60] (for sake of definiteness, we use the first set of parameters from Table I in [60]). As we have shown, the dominant contribution to the integrals over r, r in (6) comes from the vicinity of the saturation scale, |r | ∼ |r | ∼ Q −1 s , where the dipole amplitude is given by [59] . For this reason in what follows we will also use the parametrization for semi-analytical estimates of multiplicity dependences. Note that in Eq. (16) we restricted τ < 2 in accord with Fig. 2. As we can see from the structure of the cross-section (6), the contribution of the large r > 4 Q −1 s (x) (saturated region) is negligible, for this reason the parametrization (16) gives results close to the realistic dipole amplitude [60] (see Appendix C for more details).

Generalities
As was illustrated in [25], the quarkonia production crosssection (6) with CGC dipole parametrization provides a very reasonable description of the shapes of produced charmonia as a function of rapidity and transverse momentum, as well as predicts that the suggested mechanism gives the dominant contribution to J/ψ yields. The description of the multiplicity dependence presents more challenges at the conceptual level, because there are different mechanisms to produce enhanced number of charged particles N ch . The probability of multiplicity fluctuations decreases rapidly as a function of the number of produced charged particles N ch [61], and for this reason for the study of the multiplicity dependence it is more common to use a normalized ratio [18] where n = N ch / N ch is the relative enhancement of the number of charged particles in the bin, w N J/ψ / w N J/ψ and w (N ch ) / w (N ch ) are the self-normalized yields of J/ψ and charged particles (minimal bias events) in a given multiplicity class; dσ J/ψ (y, √ s, n) is the production cross-sections for J/ψ with rapidity y and N ch = Δη d N ch /dη charged particles in the pseudorapidity window (η − Δη/2, η + Δη/2).

Traditional approach
In our scenario of multi-particle production, described in the introduction, the enhanced multiplicity N ch > N ch is due to contributions from multigluon (multipomeron) configurations, as shown in Fig. 3. This mechanism leads to the multiplicity dependence (18) whereP (N IP ) is the probability to produce N IP BFKL cascades (cut pomerons). The latter parameter N IP is proportional to the number of charged particles N ch , due to the LPHD hypothesis. In the case of inclusive charged particles production in a given multiplicity class, due to the LPHD hypothesis the production cross-section is given by where c is some numerical coefficient, so the yield w (N ch ) / w (N ch ) is given byP (N IP ) . The ratio (17) evaluated with this mechanism should have a linear dependence on n = N ch / N ch ∼ N IP / N IP , which reflects the simple physics that J/Ψ can be produced by three Pomeron mechanism in each BFKL cascade. The experimental observation [26,62] that the n-dependence grows faster than linearly implies that other mechanisms must give important contributions.

Dependence of the saturation scale on the multiplicity of the soft hadrons
In the 3-pomeron mechanism analyzed in this paper we expect that the multiplicity dependence is enhanced due to a larger average number of particles produced from each pomeron (see the right panel of the Fig. 3). We expect that each such cascade ("pomeron") should satisfy the nonlinear Balitsky-Kovchegov equation. The increased multiplicity in individual pomerons leads to a modification of the dipole amplitude and the saturation scale Q s As was shown in [37][38][39][40][41]50,[63][64][65][66], the saturation scale Q s is closely related to the gluon density of the target, For high multiplicity events, in view of the LPHD hypothesis, we expect that the density of gluons in the rhs of (22) should grow linearly as a function of multiplicity [50,[63][64][65][66][67][68][69][70], so we expect that the multiplicity dependence of the saturation scale is given by While at LHC energies it is expected that the typical values of saturation scale Q s (x, b) fall into the range 0.5-1 GeV, from (24) we can see that in events with enhanced multiplicity this parameter might exceed the values of heavy quark mass m Q and lead to an interplay of large-Q s and large-m Q limits. Thus the study of the high-multiplicity events gives us access to a new regime which otherwise would require significantly higher energies.
To the best of our knowledge, currently there is no microscopic first-principle evaluation of the dipole amplitude   N (y, r, b, n) . For this reason in what follows we assume that the modification of the dipole amplitude on n is properly taken into account only by modification of the saturation scale (24).
We can derive Eq. (24) directly from our scenario for multi-particle production in CGC approach (see the introduction). Indeed, the multiplicity of soft hadrons in this scenario stems from the diagrams of Fig. 2, and is equal to where N IP is the number of cut pomerons (BFKL cascades), and c is a numerical coefficient of order one. Taking into account that the distribution d N ch /dy is almost flat, we may approximate n = N ch / N ch ≈ (d N ch /dy)/ d N ch /dy , so the assumption (24) agrees with (25) up to possible logarithmic corrections, which stem from the scale dependence of α S in the denominator of (25). Note that in Eq.24 we replace Concluding, we would like to emphasize that we suggest a new mechanism for hard particle ( quarkonia) production in events with high multiplicity of soft produced particles. In this mechanism both J/Ψ and soft hadrons are produced at short distances r ∝ 1/Q s and the difference between them is originated in different dependences of the hard cross sections on the saturation scales. This mechanism is quite different from other approaches, such as the percolation approach [19] or the modification of the slope of the elastic amplitude [20]. It can be applied both to the production of soft and hard particles. In the following subsection we will apply this approach to the multiplicity dependence of quarkonia production. It is worthwhile mentioning that the suggested approach has been successfully used for description of the experimental data in Refs. [69,70].

Charmonia production outside of the rapidity window of produced hadrons
In measuring the multiplicity there are two different situations, when the bins used for collecting charmonia and co- In the former case, we may unambiguously assign all the produced particles either to one-gluon or to the two-gluon ladders in the t-channel. Since we assume that each gluon reggeizes independently and gives equal contribution to the observed average yields of charged particles, for the two-pomeron case the enhanced multiplicity should be assigned in equal parts to both pomerons. For this reason, the cross-section (6) for n = 1 modifies as dσ y, √ s, n dy where we introduced the shorthand notation added explicitly the dependence on the relative multiplicity enhancement parameter n, and for the sake of definiteness assumed that the co-produced charged particles are collected at higher rapidity η than the rapidity y of J/ψ, as seen in the Diagram (A) of the Fig. 4 (the opposite case corresponds to an inversion of sign of the rapidity y).
The evaluation of (26) requires knowledge of the probability distribution P(n) for a pomeron to have relative enhancement of multiplicity n. The theoretical evaluation of this quantity is very involved. As was found experimentally [61], the distribution of charged particles in pp collisions for n 1 is close to exponential behavior ∼ exp(−const n) (modulo logarithmic corrections). 4 However, we would like to emphasize that this result cannot be applied to single pomerons, because a convolution of any two pomeron parts does not satisfy the convolution identity On the other hand, it was predicted long ago [72] that the distribution of particles inside a single pomeron P(N ) is described by a Poisson distribution, which clearly satisfies (28). However, in order to reproduce the experimentally observed behavior in this approach, we would have to resum an infinite number of pomerons with additional model assumptions for hadron-hadron collisions. In view of this ambiguity, in what follows we will assume that the convolution of the distributions P (n i ) in (26) just cancels in the ratio (17), and in the case when more than one pomeron spans through the region covered by the detector, the multiplicity is shared equally between pomerons, as we would get with Poisson distributions. For this reason, instead of the ratio (17) we will work with a simplified ratio where we use the notation dσ instead of dσ for the crosssection, in order to emphasize that we took out of the normalization the probability distribution of charged particles (factor P(n)). According to our assumption, the integral over n 1 in (26) gets the dominant contribution from the region when the enhanced multiplicity is shared equally between pomerons, i.e. n 1 = n/2, so we can rewrite the normalized cross-section as dσ y, √ s, n dy × Ψ g r , z Ψ J/ψ r , z 4 Such an exponential behavior is an inherent feature of the Balitsky-Kovchegov parton cascade, as it has been demonstrated in [71].
In the small-n limit (1 n m 2 Q /Q 2 s (x)) we may use the dipole amplitude (16) and get for the ratio (17) 5 At the same time, for the gluon-gluon fusion mechanism we expect that in this kinematics the corresponding n-dependence would be simply ∼ nγ , since we have only one pomeron which might contribute to charged particles in a given rapidity window. This expectation does not depend on a particular choice of the NRQCD LDMEs for a given quarkonium state. As can be seen from the left panel of Fig. 6, the three-gluon mechanism (26,17) describes the experimental data [73] reasonably well, whereas the 2-gluon contribution clearly underestimates the n-dependence.

Coinciding rapidity windows for charmonia and soft hadrons productions
The situation when the J/ψ and charged particle bins partially overlap or coincide is certainly more complicated. In the overlap region we cannot attribute the enhanced multiplicity just to a single or double gluon ladder, and instead we have to average over all possible distributions of particles n 1 , n 2 , n 3 among pomerons, such that n 1 + n 2 + n 3 = n. For this reason the cross-section for this case is given by (see Appendix A for more details) 5 We would like to emphasize that we would get a similar dependence from any quarkonium mechanism which proceeds via 3-gluon fusion. E.g., we would get the same n-dependence if we considered that the fusion of three gluons produces a color octet state which later hadronizes into quarkonium.

Fig. 6
Left: multiplicity dependence at forward rapidities (2.5 < y J/ψ < 4), when charged particles are collected at central rapidities (|η| < 1). Experimental data are from [73]. Right: evaluation at central rapidities (|y J/ψ | < 0.9). Due to overlap of the bins of charmonium and charged particles, we have an additional contribution of the nonlinear term 33, which leads to a more pronounced n-dependence. Experimental data are from [26,62]. In both plots the dashed curves labeled "Small-n approx." correspond to evaluation with a dipole amplitude (16). The curves labeled "2-gluon" in both plots correspond to n-dependence of the gluon-gluon fusion mechanism, as explained in the text. In both plots we use the LC Gauss wave function [49] (see Appendix B) for more details where we integrate over the rapidity of J/ψ inside the bin, using the variable t = (y − η min )/Δη, and over relative multiplicities in each pomeron. In the case when the rapidity bins used to collect charged particles and the J/ψ fully overlap (y min = η min , y max = η max ), the variable t ∈ [0, 1]. As we discussed earlier, currently we do not have a reliable first-principle parametrization for P(z), and for this reason in what follows we will use the simplified assumption that the full cross-section might be approximated by a sum of three contributions: when the quarkonium is produced at the center of the rapidity bin or when it is produced a either of the borders, as shown in Fig. 5. The contribution of the border regions is given by (30,31). In the center of the bin, we may assume that the observed multiplicity n is shared equally between the pomerons, so the contribution of this region is given by With the dipole amplitude (16) the expected large-n dependence of this additional term (33) is given by considerably more pronounced than (31). In the experimentally measured ratio (17) for the cross-section dσ we should take a weighted sum of the cross-sections (30) and (33). We expect that with the parametrization (16) the n-dependence should have a form  at the same time the 2-gluon fusion mechanism underestimates the data. In Fig. 7 we show the dependence of our predictions on the choice of the wave function, using for comparison boosted rest frame wave functions described in Appendix B. As expected, this dependence is very mild. In the right panel of the same Fig. 7, we show predictions for other quarkonia states which might be within the reach of experimental searches in the near future. We expect that both the ψ(2S) and the Υ (1S) states should have an ndependence close to that of the J/ψ meson. Technically this happens because the typical dipole sizes r ∼ r ∼ min m (33) are significantly smaller than the sizes of the quarkonia states (see Appendix C for details), and in this approximation the wave function Ψ M of the latter might be approximated with its value at r = 0. This leads to an independence on quarkonium size or mass for the selfnormalized ratio (17,29) . The residual dependence seen in the data is due to corrections beyond the heavy quark mass limit.
In Fig. 8 we have shown results for the J/ψ multiplicity dependence, in the RHIC kinematics, which are in a very good agreement with experimental data [15,16]. We also show predictions for Υ (1S) and ψ(2S), which are close to the results for J/ψ. We hope that these predictions will be checked very soon.

Charmonia production outside of the rapidity window of produced hadrons
Finally, we would like to mention that the contribution of the last term of Eq. (32) (see also Eq. (35)) might be singled out if the width of the pseudorapidity bin Δy used for collection of quarkonia is significantly smaller than the width of the bin Δη used for collection of charged particles and it is included inside (Δy Δη), (y min , y max )⊂(η min , η max ). In this case the contribution of the border region is negligible, and we expect that the n-dependence is given only by Eq. (34). We hope that in the future the experimentalists will be able to check this prediction. Also, we would like to remind that the other approaches [19,74] do not predict such a rapid n-dependence (Eq. (34)), and for this reason its possible observation could be a decisive test in favor of the three-gluon mechanism of quarkonia production.

Conclusions
In this paper we have developed, in the framework of the CGC approach, a new mechanism for hard particle ( quarkonia) production in events with high multiplicity of the soft produced particles. In this mechanism both J/Ψ and soft pions are produced at short distances r ∝ 1/Q s and the difference between them originates in different dependences of the hard cross sections on the saturation scale. Based on this formalism we analyzed the multiplicity dependence of the twogluon and three-gluon fusion mechanisms, assuming that the observed multiplicity dependence is due to an increase of the average number of particles produced from each CGC pomeron. We found that the two-gluon mechanism significantly underestimates the multiplicity dependence, whereas the predictions of the three-gluon mechanism are in reasonable agreement with available experimental data in a wide energy range, from RHIC to LHC. This conclusion does not depend on the particular choice of the NRQCD LDMEs for a given quarkonium. We believe that this argument is a strong evidence that the contribution of the three-gluon mechanism might be substantial in J/ψ production. We call the experimentalists to measure the multiplicity dependence of other charmonia (e.g. ψ(2S) and Υ (1S)), for which we expect approximately the same multiplicity dependence as for J/ψ. Such measurements would be an important argument in favor of 3-gluon mechanism in the CGC approach. We also predict that the suggested mechanism should provide a strong ∼ n 3γ dependence, in the kinematics when J/ψ is produced at the center of the charged particle bin, i.e. the width of the bin Δy Δη. Our predictions are sensitive to the large-distance behavior of the dipole amplitude, and going to higher values of n potentially could allow to test the saturation model implemented in a dipole amplitude. However, the experimental study of this regime is challenging due to a rapid decrease of probability (yields) of such high-multiplicity events [61], which puts it out of reach of modern and forthcoming accelerators.
The mechanism of multiplicity dependence suggested in this paper is not unique: multiplicity might also be enhanced due to large number of pomerons, as shown in the left panel of Fig. 3 and in Fig. 9. We should also take into account that each of the ladders (pomerons) shown in Fig. 9 might lead to development of a Balitsky-Kovchegov cascade, as shown in the right panel of Fig. 3. We plan to address these issues in a future publication [75]. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .

Appendix A: Derivation of the Eq. (32)
In this section we will evaluate the cross-section of the pion production case, when the bins used for collection of quarkonia and charged particles overlap. In this setup we cannot attribute the observed increase of multiplicity to any of the pomerons and have to evaluate probability-weighted sum over all possible partitions of the number of charged particles N ch into three integers. Since the quarkonium might be produced at any rapidity y (not necessarily at the center of the bin), this implies that the average contribution of each pomeron to the enhanced multiplicity seen in the bin Δη will also depend on y. In what follows it is convenient to work with a variable t = (y − η min )/Δη, which for the case of the full overlap (y min = η min , y max = η max ) remains in the interval t ∈ [0, 1].
In what follows we will use the notations N (1) ch , N (2) ch and N (3) ch for the number of charged particles produced from each of the pomerons, with N (1) ch + N (2) ch + N (3) is the total number of charged particles produced inside the bin Δη. The probability that a given pomeron of length η tot units of rapidity produces N (tot) ch particles is given by P N (tot) is the average number of particles produced by the whole pomeron. Since the distribution d N ch /dη almost does not depend on rapidity, we may assume that inside each interval Δη i (i = 1, 2, 3) the distribution of produced charged particles is given by the same function, P N , where the average number Δη i /η tot , and the length of the pomeron is given by According to (A.1), the probability that the three pomerons contribute to the observed multiplicity enhancement is given by where we introduced the variables n i = N If the width of the bin Δη is sufficiently large (so that N (Δη) ch 1), we may replace the summation over N (1) ch , N ch with an integral over n 1 , n 2 and obtain the Eq. (32).

Appendix B: Wave functions
In our evaluations, unless stated otherwise, we use the socalled light-cone wave functions of 1S quarkonia states [49], where N and R are some numerical constants. In the heavy quark mass limit, as well as in the deeply saturated regime, we expect that the results should not depend at all on the choice of the wave function, since only small dipole configurations in (6) give the dominant result. In the small-r approximation we may replace the wave function Ψ M (z, r ) with its value at the point r = 0, which eventually cancels in the ratio (17). However, in order to assess the accuracy of this approximation, as well as to have the wave functions of other quarkonia (primarily ψ(2S) and Υ (1S) are of interest for experimental searches), we also used wave functions evaluated in rest frame potential models and boosted them to a moving frame in momentum space. Technically, this prescription is given by the relation [49,76] Ψ (LC) where ψ RF is the Fourier image of the eigenfunction of the corresponding Schroedinger equation. For comparison we considered three different choices of the potential: -The linearly growing Cornell potential which is discussed in detail in [77,78] -The potential with logarithmic large-r behaviour suggested in [79], We checked that the wave functions obtained with all three potentials have similar shapes, and in the case of J/ψ are close to the light-cone Gauss potential (B.6). In the heavy Fig. 10 Left: Function K r, r ∼ Ψ g · Ψ M (r ) Ψ g · Ψ M (r ) dφ ΔN r, r which appears in the integrand of (32,33) for x = 10 −3 and n = 8. Right: The same integrand for the hypothetic case of light meson production, with wave function (B.6) with parameter (meson size) R ∼ 1 fm and light quark mass m q ≈ 0.3 GeV. In both plots the semitransparent pink surface corresponds to evaluation with CGC parametrization, and the semitransparent green surface corresponds to evaluation with parametrization (16). The latter has a visible cusp due to the discontinuity of the dipole amplitude (16) quark mass limit the wave functions of the gluon Ψ g in (1) can be approximated with the perturbative expressions from [49]. Their overlap with the wave function of quarkonium are given by

Appendix C: Small-n approximation
As was illustrated in Sect. 3.1.1, the approximation (16), which was derived for r in the vicinity of the saturation scale (r Q s ∼ 1), works up to surprisingly large values of n. In this section we would like to discuss briefly the technical reasons why this happens.
In LHC kinematics, we expect that the typical sizes of r and r in the integrand of (26) are controlled by the wave functions Ψ g , for events with normal multiplicity (n ≈ 1). In the heavy quark mass limit this gives an estimate r ∼ r ∼ m −1 Q . However, the situation changes for events with elevated multiplicity (n 1), because the saturation scale Q s (x, n) increases and at some point becomes comparable and even exceeds the heavy quark mass m Q . In the deeply saturated regime Q s (x, n) m Q , so the factor ΔN defined in (27) leads to a strong suppression of dipoles with a size larger than Q −1 s (x, n). For this reason we can see that the integrand of (6,26) always gets the dominant contribution from the region r Q s , which explains why the predictions of (16) for self-normalized ratios are close to the more realistic CGC parametrization [60].
In order to illustrate these findings numerically, in the left panel of Fig. 10 we have plotted the function K r, r ∼ Ψ g · Ψ M (r ) × Ψ g · Ψ M (r ) dφ ΔN r, r , which appears in the integrand of (32,33). For the sake of definiteness we have chosen the point n ≈ 8, where the expected deviations of the parametrization (16) from CGC are the largest. We can see that for charmonia the difference between the integrands is of order 10%, in agreement with our earlier results. For the sake of comparison, in the right panel of the Fig. 10 we have also plotted a similar integrand for the hypothetic case of light meson (e.g. ρ-meson) production, assuming that its wave function is given by (B.6) with an adjusted parameter R ∼ 1 fm (size) and using the light quark mass m q ≈ 0.3 GeV. We can see that in this case the discrepancy is significantly larger, i.e, the small-n approximation (parametrization (16)) cannot be used for estimates of multiplicity dependence.