High energy QCD: multiplicity dependence of quarkonia production

In this paper we propose an approach which demonstrates the dependence of quarkoni production on the multiplicity of the accompanying hadrons. Our approach is based on the three gluons fusion mechanism, without assuming the multiplicity dependence of the saturation scale. We show, that we describe the experimental data, which has a dependence that is much steeper than the multiplicity of the hadrons.

Fortunately, in Ref. [21] it was shown, that these two approaches are equivalent for the description of the scattering amplitude in the rapidity range : where ∆ BFKL denotes the intercept of the BFKL Pomeron. In this paper it is also shown, that for Eq. (1) we can use the Mueller, Patel, Salam and Iancu approximation(MPSI) [25] for hadron-hadron scattering at high energies.
The recent experiments by ALICE [26][27][28][29][30] and STAR [31,32], show that the cross sections for J/Ψ production strongly depends on the multiplicity of accompanying hadrons. These data have stimulated theoretical discussions on the origin of such dependence (see Refs. [33][34][35][36][37]). In this paper, we develop an approach to this problem based on two ingredients. First, we assume ,that the production of quarkonia stems from the triple gluon fusion [35,39,40] (see Fig. 1). For the interaction with nuclei this mechanism is dominant [41][42][43][44][45]; and it has been demonstrated in Ref. [35], that this mechanism gives a substantial contribution in hadron-hadron collisions.  Fig. 1-a: The three gluon fusion mechanism of J/Ψ production. Fig. 1-b: the Mueller diagram [38] for J/Ψ production, which illustrates the inter-relation between three gluon fusion and the triple BFKL Pomeron interaction. The wavy lines describe the BFKL Pomerons, while the helical curves represent gluons.
Second, we showed in Refs. [46,47], that in spite of the fact that in different kinematic regions, the QCD cascade leads to a different energy and dipole size dependence of the mean multiplicity, the multiplicity distribution has a general form: where N denotes the average number of partons. The paper is organized as follows: In the next section we describe our approach to hadron-hadron collisions. In section III we discuss quarkonia production in simplified Reggeon Field Theory, defined in zero transverse dimensions. In section IV we generalized the result in this toy-model approach to high energy QCD, and compare our estimates with the experimental data . We summarize our results in the Conclusions.

II. HADRON-HADRON INTERACTION IN MPSI APPROACH
This section does not contain new results, and we include it in the paper for completeness of presentation, as well as a kind of an introduction to the notation, and the main ideas.

A. QCD parton cascade
We start with the equation for the QCD parton cascade which can be written in the following form. [1,3,17,18]: where P n (Y ; {r i , b i }) denotes the probability to have n-dipoles of size r i , at impact parameter b i , and at rapidity is a typical cascade equation in which the first term describes the depletion of the probability of n, due to one dipole decaying into two dipoles of arbitrary sizes, while the second term describes, the growth due to the splitting of (n − 1) dipoles into n dipoles.
The initial condition for the DIS scattering is which corresponds to the fact that we are discussing a dipole of definite size which develops the parton cascade. Since P n (Y ; {r i }) is the probability to find dipoles {r i }, we have the following sum rule i.e. the sum of all probabilities is equal to 1. This QCD cascade leads to the Balitsky-Kovchegov (BK) equation [1,4,5] for the amplitude, and gives the theoretical description of DIS. To see this we introduce the generating functional [3] where u (r i b i ) ≡ = u i is an arbitrary function. The initial conditions of Eq. (4) and the sum rules of Eq. (5) require the following form for the functional Z: Multiplying both terms of Eq. (3) by n i=1 u (r i b i ) and integrating over r i and b i , we obtain the following linear functional equation [18]; Searching for a solution of the form Z ([u(r i , b i , Y )]) for the initial conditions of Eq. (7a), Eq. (8a) can be re-written as the non-linear equation [3]: Therefore, the QCD parton cascade of Eq. (3) takes into account non-linear evolution. Generally speaking the scattering amplitude can be written in the form [5,18]: where N (Y 0 , r i , b i ) is the amplitude of the interaction of dipole r i with the target at low energy Y = Y 0 , and the n-dipole densities in the projectile ρ p n (r 1 , b 1 , . . . , r n , b n ) are defined as follows: For ρ n we obtain [18] : For ρ 1 we have the linear BFKL equation [8]: However, to obtain the BK equation for the scattering amplitude we need to use Eq. (10), in which we introduce the amplitude of interaction of the dipole with the target at low energies. Using Eq. (8a),Eq. (10) and Eq. (11) , we can obtain the non-linear BK equation from Eq. (9) in the following form [5] B. The interaction of two dipoles at high energies We first consider the simplest case of scattering, the high energy interactions of two dipoles with sizes r and R and with r ∼ R. In Ref. [21] it is shown that in the limited range of rapidities, which is given by Eq. (1), we can safely apply the Mueller, Patel, Salam and Iancu approach for this scattering [25](see Fig. 2-a). The scattering amplitude in this approach can be written in the following form [18]: where ρ t n and ρ p n denote the parton densities in the target and projectile, respectively. These densities can be calculated from P n using Eq. (11). N BA is the scattering amplitude of two dipoles in the Born approximation of perturbative QCD (see Fig. 2). Eq. (15) simply states that we can consider the QCD parton cascade of Eq. (3) generated by the dipole of the size r for the c.m.f. rapidities from 0 to 1 2 Y , and the same cascade for the dipole of the size R, for the rapidities from 0 to − 1 2 Y . One can see that Eq. (15) is the t-channel unitarity re-written in a form, convenient for applying the evolution of the parton cascade in the form of Eq. (12).
Generally speaking, for the dense system of partons at Y = 0 n-dipoles from upper cascade could interact with m dipoles from the lower cascade, with the amplitude N m n {r i }, {r j } [18]. In Eq. (15) we assume that the system of dipoles that has been created at Y = 0 is not very dense, at least for the range of rapidities given by Eq. (1). In this case and after integration over {r 1 } and {r j }, the scattering amplitude can be reduced to a system of enhanced BFKL Pomeron diagrams, which are shown in Fig. 2 The average number of dipoles at Y = 0 is determined by the inclusive cross section, which is given by the diagram of Fig. 2-c and which can be written at y → 0 as follows [48]: The average number of dipoles that enters the multiplicity distribution of Eq. (27) is equaln = N = Indeed, the enhanced diagrams of Fig. 2-b lead to the inelastic cross section which is constant at high energy.

C. Hadron -hadron collisions
In this paper we view a hadron as a dilute system of dipoles and use Eq. (17) for the average multiplicity, together with the multiplicity distribution of Eq. (1). In particular, we assume that Eq. (16) is correct, and the system of partons that is produced at c.m. rapidity y * =0 is a dilute system. However, we are aware that Eq. (17) does not describe the experimental increase of the average multiplicity, which from Eq. (17) isn ∝ exp (∆ BFKL Y ). The experimental data can be described in the framework of the CGC/saturation approach in which N BFKL were replaced by N BK [51]. Hence, we cannot view hadrons as a dilute system of dipoles, but rather have to consider them as a dense system of dipoles. For such a situation we expect thatn ∝ Q 2 s (Y )/ᾱ S (see Refs. [1,[49][50][51][52].

A. Multiplicity distribution -a recap
In the parton model [53][54][55] all partons have average transverse momentum which does not depend on energy. Therefore, we can obtain the parton model from the QCD cascade assuming that the unknown confinement of gluons leads to the QCD cascade for a dipole of fixed size. In this case the cascade equation (see Eq. (3)) takes the following simple form: where P n (Y ) denotes the probability to find n dipoles (of a fixed size in our model) at rapidity Y , and ∆ is the intercept of the BFKL Pomeron. Instead of the generating functional of Eq. (6), we can introduce the generating function: where u are numbers. At the initial rapidity Y = 0, we have only one dipole, so P 1 (Y = 0) = 1 and P n>1 = 0 (so the state is the only one dipole); at u = 1, Z (Y, u = 1) = n P (y) = 1. These two properties determine the initial and the boundary conditions for the generating function which simplify Eq. (7a) and Eq. (7b) Eq. (18) takes the following form for the generating function: The general solution to Eq. (21) is an arbitrary function (Z (z)) of the new variable: z = ∆ Y + f (u), with f(u) from the following equation: The form of arbitrary function stems from the initial condition of Eq. (20) Note, that Z (Y, u = 1) = 1, as it should be from Eq. (20).
On the other hand, we can re-write Eq. (21) in the form of the non-linear equation using Eq. (22): viz.
Comparing Eq. (24) with Eq. (19) one can see that Since from Eq. (26) it follows that the average n = N is equal to N = exp (∆ Y ) Eq. (26) can be re-written in the form of Eq. (2):

B. Quarkonia production
As we have discussed in the introduction, we assume the production of heavy quakonia stems from three gluon fusion (see Fig. 1), and it is intimately related to the triple Pomeron interaction. The Mueller diagrams for inclusive J/Ψ production are shown in Fig. 3, where the wavy lines denote the Pomeron Green's function (G I P ) which is equal to The MPSI approach for inclusive production with fixed multiplicity of produced hadrons is shown in Fig. 4 Therefore, from this figure we see that the structure of the parton cascade for the quarkonia production is quite different. In particular, for the part of the events whose weight is determined by the contribution of Fig. 3-b, the initial conditions for the parton cascade is not the ones of Eq. (20) but they have the form: This means, that for this cascade we need to find the arbitrary function Z (z) from the following equation The solution is Eq. (31) leads to a different multiplicity distribution in comparison with Eq. (27): viz.
Finally, the cross section for the quarkonia production is proportional to where κ is the weight of the contribution of Fig. 3-b to the contribution of Fig. 3-a, which is equal to In Eq. (33) P  The cross section of produced gluons (hadrons) is proportional to From Fig. 5, one can see that P (2) n (N ) and P (1) n (N ) have different dependance on z = n/N . The average number of gluons is chosen to be the mean multiplicity of hadrons in the rapidity window |η| ≤ 0.9 measured at W=13TeV. Therefore, dσ J/Ψ dy is not proportional to dσ prod. gl.i dy , but shows non-linear dependance, which we will discuss below. It should be stressed that such dependance stems from triple Pomeron mechanism of quarkonia production, as noted in Refs. [35,37]. At large N we have leading to It is worth mentioning that the average multiplicity of P (2) distribution is equal to Hence, for the multiplicity distribution P (2) the average number of accompanying gluons(hadrons) is twice larger than in distribution P (1) . Eq. (38) means that the ratio n (2) <n (2) > = n 2 N .

A. Multiplicity distribution
In QCD, to find the multiplicity distribution for hadron-hadron scattering in QCD using MPSI approach [25], we need to evaluate (see Eq. (6) However, in the MPSI approach it is more natural to introduce moments (see Eq. (12)): The integration over r i depends on the size of the initial dipole, which generates the cascade. In DIS the natural integration stems from r i > r.

Several first iterations.
We start from the first several iteration of Eq. (12), which can be re-written forρ p n (r 1 , b 1 . . . , r n , b n ) in the form For the first iterationρ 1 (Y ; r 1 , b 1 ), we obtain the BFKL equation: where ξ 1 = ln r 2 1 Λ 2 QCD and diffusion approximation where ψ(z) is Euler gamma function (see [57] formula 8.36).ρ p in,1 (γ, b) has to be found from the initial conditions. From Eq. (43) we can obtain M p 1 (Y, r, b) (see Eq. (40)) which has the following form: which satisfies the following equation: The equation for the next iteration:ρ 2 , takes the form: Equation for ρ p 2 can be re-written in the following form for ρ p 2 : For simplicity we re-write Eq. (47a) in the log approximation following Ref. [46] (see Eq. (47b)). Rewriting Eq. (47b) for M p 2 we obtain: In the last term of Eq. (48) we used Eq. (46). The solution to Eq. (48) takes the form: with χ (γ) = 1/γ. One can check that M p 2 (Y, r) which is the correct initial condition for one dipole of size r at Y = 0 , which generates the parton cascade. Actually, Eq. (48) describes M p 2 (Y, r) for the full BFKL kernel. Indeed, considering Eq. (47a) we can integrate this equation over r 1 and r 2 , to obtain the equation for M p 2 . The last term has the following form Note that r 2 2 = (r 1 − r 12 ) 2 can never approach zero, since r 2 > r. Removing this restriction, we can re-write The reggeization term in Eq. (51) describes the contribution of r 2 → 0. Plugging Eq. (51) in the last term of Eq. (47a) integrated over r 1 and r 2 , one can see that we reproduce Eq. (48) for the full BFKL kernel. Hence Eq. (49) is the solution with χ (γ) which is given by Eq. (44).

General solution
The equation for M p n (Y, r) has the following general form: The solution to this equation [46], which gives M p 1 (Y = 0, r) = 1, but all other M p n with n ≥ 2=0, are equal to which leads to the multiplicity distribution, which takes the form (see Eq. (39) and Ref. [46]): Taking the integral over γ using the method of steepest descent, using the diffusion approximation for the BFKL kernel (see Eq. (44)). The equation for γ SP has the following form: and the integral over γ is considering ξ 2 2 D z 1 2 Y 1. After normalization, we obtain that where Ψ KNO denotes the KNO function (see Ref. [58]). It is worthwhile mentioning that the multiplicity distribution of Eq.
In Fig. 6 we compare the ALICE data [59] on multiplicity distribution with Eq. (57) and with Eq. (27). On can see that the agreement is good, and the difference between the above equations can be seen at large n. In describing the experimental data we use Eq. (57), which is derived at large z, for z ≥ 3. It worthwhile mentioning that the data of CMS [60] we have discussed in our paper [46].

B. P
(2) n distribution for quarkonia production The P (1) n distribution , which we have discussed in the previous section, can be derived, using the double Laplace transform representation, both for P where ξ = ln r 2 Λ 2 QCD . Eq. (52) in the ω-representation, has the following form: with the solution: This solution gives M p 1 (Y = 0, r) = 0, while M p n (Y = 0, r) = 0 for n ≥ 2. The inverse Laplace transform leads to Eq. (53) for M p n (Y, r) and Eq. (54) for P (1) n (Y, r). To find P (2) n distribution we need to take into account that at Y = 0: One can see that the following m 1 (ω, γ) satisfies this conditions: The inverse Laplace transform with respect to ω leads to which gives the initial conditions of Eq. (63). For P (2) n we obtain with P 2 (Y = 0, r) = 1 and P (2) n (Y = 0, r) = 0 for n = 2 at Y = 0.
Repeating the same estimates as in Eq. (55), we obtain the KNO function, w with the normalization dz Ψ KN O (z) = 1.

V. COMPARISON WITH EXPERIMENTAL DATA
In both ALICE [26][27][28][29][30] and STAR [31,32] experiments the following ratio is measured: where N = n is the average number of charged hadrons in the fixed rapidity window, and n J/Ψ the average number of J/Ψ which are measured generally speaking in a different rapidity window. It turns out that F n N = n N , but it is close to this when the rapidity windows are different. When both rapidity windows are the same, F n N shows much steeper dependence than n N . In Ref. [34] the J/Ψ production is considered as being proportional to the number of collisions, since it comes from short distances, while the production of hadrons is proportional to the number of participants (see Ref. [49].) However, in the framework of the CGC approach, the J/Ψ production at high energies is proportional to the number of participants [35,[41][42][43] as it can be seen from Fig. 1.
The main ingredients for describing the experimental data are Eq. (32) and Eq. (35). As one can see from Fig On the other hand, the multiplicity of the quarkonia in the cascade of Fig. 4-b is equal to In Eq. (70) we took into account that the average number of produced partons (hadrons) is equal to 2 N . For large n/N = z one can see that Eq. (70) leads to In Fig. 7 we compare the experimental data with Eq. (72).
One can see that this simple formula provides a fairly good description of the experimental data, for central production where κ = 1. However, the experimental data for forward production of J/Ψ [26][27][28][29][30][31][32] show almost a linear dependance: n J/Ψ n J/Ψ = z. Indeed, in Eq. (72) the quadratic term is suppressed, since it is proportional to the value of κ, which is equal to (see Fig. 4 and Ref. [35] for the estimates).
pp W 13 TeV , |y|<0.9 where y * is the rapidity of the produced quarkonia in c.m.f. In leading order of perturbative QCD , in which we made all our previous estimates [1],γ = 0.63 and λ =ᾱ S χ(γ) γ ≈ 4.8ᾱ S . As is shown in Fig. 7-b, the estimate in leading order describes the data quite well. However, we need to remember that the NLO corrections both toγ and to λ are large.

VI. CONCLUSIONS
In this paper we re-visited the problem of multiplicity distributions in high energy QCD, which we have discussed in Ref. [46] and found the distribution of Eq. (57). This distribution provides a better description of the experimental data at large multiplicities n, than Eq. (27), which has been discussed previously. We also suggest a different approach to the multiplicity dependence of quarkonia production. It should be stressed that our approach is based on the three gluons fusion mechanism of Fig. 1, and it differs from the description of Refs. [36,37], since we did not assume the multiplicity dependence of the saturation scale. In our approach we assume, that the production of J/Ψ, which occurs at rapidity y, and the central production of charged hadrons, stem from the production of the same n-parton cascades, which are pictured in Fig. 8 as the production of n-gluon ladders.
In Fig. 8 one can see that J/Ψ can be produced from each of n-ladders, leading to the cross section, which is proportional to n (see Fig. 8-a). This mechanism is shown in Fig. 3-a. However, J/Ψ can be created from merging of two ladders (see Fig. 8-b) , which gives a cross section ∝ n 2 , and corresponds to Fig. 3-b. Note, that the production of the hadrons in both cases are proportional to n. Taking into account that the average number of gluons (hadrons) for the mechanism of Fig. 8-b is two time larger than for Fig. 8-a, we infer that Fig. 8 leads to the simple Eq. (72).
It should be stressed that this equation is heavily dependent on the three gluon fusion mechanism, but does not depend on the details of the cross section of quakonia production. In particular, as we have mentioned above, we do not use the dependence of the saturation scale on the multiplicity of the produced gluons. This means that the non-linear dependence of J/Ψ-production on multiplicity of charged hadrons, can stem from sources other than the dependence of the cross section on the saturation scale. Actually, this statement follows directly from the fact that 1+1 RFT generates the non-linear dependence on n.
As an aside, we note in [36,37] it was assumed that the J/Ψ-production results from the inclusive diagrams of Fig. 3-a and Fig. 3-b. This is erroneous, as at arbitrary n it is necessary to include the production of many partonic showers as illustrated in Fig. 4-a and Fig. 4-b. The difference to the percolation approach [34], lies in our hypothesis that both the production of J/Ψ and the charged pion stem from short distances of the order of r ∝ 1/Q s , and are determined by physics controlled by the CGC effective theory. The gluon jets with transverse momentum Q s , decay into charged pions (see Ref. [37] for details). The non-linear dependence of production of J/Ψ is due to the three Pomeron fusion mechanism.
In spite of the good description of the experimental data for the quarkonia production integrated over the transverse momenta (p T ), we cannot explain at the moment, why the data at fixed p T [26], shows a steeper dependence on n than the integrated data. Certainly, this problem will be the main subject of our further attempts to understand the multiplicity dependence of quarkonia production.

VII. ACKNOWLEDGEMENTS
We thank our colleagues at Tel Aviv university and UTFSM for encouraging discussions. This research was supported by ANID PIA/APOYO AFB180002 (Chile) and Fondecyt (Chile) grants 1180118.