$J/\psi$ production in hadron scattering: three-pomeron contribution

In this paper we discuss the inclusive $J/\psi$ production in proton-proton collisions from fusion of three pomerons. We demonstrate that this mechanism gets dominant contribution from the region which can be theoretically described by CGC/Saturation approach. Numerically, it gives a substantial contribution to the $J/\psi$ production, and is able to describe the experimentally observable shapes of the rapidity, momenta and multiplicity distributions. The latter fact provides a natural explanation of the experimentally observed enhancement of multiplicity distribution in $J/\psi$ production.

The description of the charmonium hadroproduction remains one of the long-standing puzzles since its discovery.The large mass m c of the charm quark inspired applications of perturbative methods and consideration in the formal limit of infinitely heavy quark mass [1].However, in reality the coupling α s (m c ) ∼ 1/3 is not very small, so potentially some mechanisms suppressed in the large-m c limit, numerically might give a sizeable contribution.
The Color Singlet Model (CSM) [2][3][4] assumes that the dominant mechanism of charmonia production is the gluongluon fusion supplemented by emission of additional gluon.Early evaluations in the collinear factorization framework did not agree with the experimental data at large transverse momenta p T by several orders of magnitude.The failure of the expansion over α s due to milder suppression of higher order terms at large p T [5,6] and co-production of additional quark pairs [7,8] motivated introduction of the phenomenological Color Octet contributions [9,10].The modern NRQCD formulation [11][12][13][14][15] constructs a systematic expansion over the Long Distance Matrix Elements (LDMEs) of different charmonia states which can be extracted from fits of experimental data.However, at present extracted matrix elements depend significantly on the technical details of the fit [14,18,19], which contradicts expected universality of the extracted LDMEs.At the same time, it is known that at large p T , a sizeable contribution might come from other mechanisms, like for example gluon fragmentation into J/ψ or co-production together with other hadrons [20][21][22][23].The latter findings are partially supported by experimental data on multiplicity of co-produced charged particles [24][25][26][27][28] which suggest that J/ψ production might get sizeable contribution from this mechanism.
In this paper we analyze J/ψ in the CGC/Saturation approach, which incorporates the leading-log contributions from hadron production in the small-x kinematics, and takes into account saturation effects in the region of very small-x [39,59].We focus on the mechanism of J/ψ production shown in the diagram (a) of the Fig. (1), when the production of J/ψ occurs from fusion of three gluons, accompanied by the production of two parton showers 1 .In order to emphasize the role of co-production of other particles associated with J/ψ, below we will follow the terminology used in some BFKL papers and refer to this mechanism as "two-parton shower" contribution.This mechanism differs from the so-called "single-shower" mechanism shown in the diagram (b) of the Fig. (1).The latter corresponds to gluon-gluon fusion with emission of additional (soft) gluon in collinear and k T -factorization approaches, and is a counterpart of CSM mechanism in BFKL framework.We expect that the two-parton shower mechanism should be dominant for the events with large multiplicities, exceeding the average multiplicity n of the gluon production in the inclusive process.This mechanism is similar to mechanisms studied earlier in the literature for proton-proton collisions [29,30] and for proton-nucleus and nucleus-nucleus scattering [31][32][33].For collisions involving heavy nuclei, the two parton shower mechanism inside heavy nuclei with atomic number A is dominant due to factor ∼ A 1/3 , and, because of this, has been comprehensively discussed during the past decade [34][35][36][37].On the other hand, for proton-proton scattering we found only one paper which considers this process in the k T -factorization approach [38]. ) is not suppressed by any hard scale and thus does not bring a smallness proportional to ᾱS .In the DGLAP approach, it is expected that such diagrams are of higher twist and are (at least formally) suppressed by additional powers of hard scale.However, at high energy this suppression is compensated by enhanced contribution of the second parton shower (see [39] for review).For typical r ∼ 1 GeV −1 we get ᾱS 4/r 2 ≈ 0.2 so ᾱS G BFKL (s, . . . ) ∝ ᾱS s ∆BFKL ≥ 1 at high energies, where ∆ BFKL = 4 ln 2 ᾱS is the intercept of the BFKL Pomeron [40] and G BFKL (s, . . . ) is the Green function of the BFKL Pomeron.The numerical estimates of Ref. [38] suggests that the considered mechanism yields about one third of the experimental cross section in the k T factorization approach, though the final estimate suffers significantly from the uncertainty in digluon PDF modeling, namely the choice of value of the parameter σ eff ("effective double parton cross-section").On the other hand, the diagram (b) in the Fig ( 1) at high energies gets additional suppression due to growth of the saturation scale Q s , which decreases the average dipole size and suppresses the emission of the extra gluons leading to α S (Q s ) suppression.
The main motivation of this paper is to re-visit the estimates of the contribution of the mechanism of Fig.
(1-a) in the CGC/Saturation framework and check if it can reproduce the observed multiplicity distributions.We demonstrate that: i) this mechanism can be calculated in CGC/Saturation approach (see Ref. [59] for a review) since the main contribution comes from the vicinity of the saturation scale, where we know theoretically the scattering amplitude, and (ii) it on its own gives a significant contribution to experimentally observable cross sections.In contrast to Ref. [38], we use a CGC/Saturation framework, and in order to avoid uncertainties related to digluon distributions, we relate the cross-sections of the J/ψ production to the diffractive production cross-section known from DIS.
The paper is organized as follows.In the next section II we evaluate the contribution of the suggested mechanism in the CGC/Saturation framework.We re-write this contribution in the coordinate representation and relate it to the gluon double densities.In Section III we discuss the interrelation of the suggested process with the diffractive production of J/ψ in DIS.In the Section IV we make phenomenological estimates and compare results with experimental data.The subsection IV A is dedicated to the estimates of the total cross section of the process.In this subsection we point out the differences with Ref. [38].In subsections IV B, IV C and IV D we consider the momenta, rapidity and multiplicity distributions of the differential cross-sections.Finally, in the Section V we draw conclusions.
The CGC/saturation approach, being more microscopic, describes the high energy interactions in terms of colorless dipoles, their density, distribution over impact parameters, evolution in energy, etc.The distinctive feature of this approach is the appearance of saturation effects which affect the dynamics for parton momenta comparable to some saturation scale Q s , a new dimensional parameter.The studies of J/ψ production in this approach may be found in [36] and by construction include all possible multishower contributions, although with additional model-dependent assumptions.
The BFKL Pomeron calculus works with BFKL Pomerons and their interactions, and phenomenologically is similar to the old Reggeon theory [55].This approach is suitable for describing diffractive physics and correlations in multiparticle production, so we can use the Mueller diagram technique [56].The relation between different processes at high energy are very often more transparent in this approach, since in addition to the Mueller diagram technique we can use the AGK cutting rules [57], which are useful in spite of the restricted region of their application [58].In this paper for the sake of definiteness we use the BFKL Pomeron calculus as a framework for evaluations.
In  Its contribution to the total cross section for J/ψ production2 is equal to where G cut I P are the Green functions of the cut pomerons (it is related to elastic amplitudes N BFKL cut by Fourier transform), the gluon momenta k T , p T , q T , are defined in the Fig.
(the evaluation of the factor I (k T , q T ) is discussed in more detail in Appendix A).The additional factor 4 in Eq. ( 8) which comes from two sources: 2 from the AGK cutting rules [57] and 2 from the fact that gluon that produces cc-pair can come from another proton.In Eq. (3) Ψ g (p T , r) stands for the wave function of gluon with virtuality p T , transverse quark-antiquark separation r and the light-cone fraction of the quark z, while Ψ J/ψ is the wave function of J/ψ meson.The amplitude of the BFKL Pomeron G cut I P y, k T ± 1 2 q T , Q T can be simplified if we take into account that the Q T dependence of the BFKL Pomeron is determined by the size of the largest of the interacting dipoles 3 , and accounting for Eq. ( 1), may be written as The dependence on Q T is described by S h (Q T ) in Eq. ( 5), which has the non-perturbative origin and has to be taken from the experiment.
Using Eq. ( 5) we can re-write Eq. ( 2) in the form For further evaluations it is very convenient to introduce momenta p 1,2,T = ±k T + 1 2 q T , which allow to rewrite Eq. ( 6) as where we took the integral over p T ∈ (0, 2m c ) using √ s, and Y − y ≡ ln (1/x g ).Making a Fourier transform, we may rewrite Eq. ( 7) in the coordinate space as The fact that the Q T dependence is determined by the size of the largest dipole stem from the general features of the BFKL Pomeron.Indeed, the eigenfunction of the BFKL Pomeron in the coordinate space is equal to [47] N r, r where b is the conjugate variable to Q T .From Eq. ( 4) one can see that the typical value of b is of the order of the largest of r and r .In our process r is of the order of R h , where R h denotes the radius of the nucleon.The value of 1/r is of the order of the mass of the heavy quark mc, or the saturation scale Qs and, therefore, turns out to be much larger than 1/R h , and can be neglected.The dependence on Q T ≈ 1/R h is described by S h (Q T ) in Eq. ( 5), which has the non-perturbative origin and. in practice, has to be taken from the experiment.
where the amplitudes N (y; r i ) are related to the solutions of the BK equation as N (y; r i ) = d 2 b N y; r i , b , and the variable b is a Fourier conjugate of momentum p 1,T + p 2,T − q T .In what follows, we will also need an expression for the q T -integrated cross-section, which takes a simpler form Finally, we would like to mention that the expressions presented in this section implicitly assume that each parton cascade is emitted independently, namely that parton correlations are negligible.In more general case with nonzero parton correlations, the expression (6) should be replaced with where we replaced the product of pomeron propagators with the double transverse momentum densities ρ (2) defined as where (x 1 , p 1,T ) and (x 2 , p 2,T ) are the light-cone and transverse momenta of the partons in the cascade, 4 |P is the Fock state of colliding hadrons, a + and a denote the creation and annihilation operators for gluons that have longitudinal momentum x i and transverse momentum p i,T , c i are the color indexes.However, at present there is no experimental evidence that such correlations are large, for this reason in what follows we will use a simpler expressions (6).Similarly, in coordinate space the Eq. ( 8) can be extended as where ρ (2) (x, r; x, r ; Q T ) is the double parton density in the coordinate representation, r s = 1 2 (r + r ) and While for numerical estimates we could use parameterizations of the amplitude N available from the literature, we would like to minimize dependence on parameterization and make a few model-independent estimates.One of the important parameters for understanding the small-x dynamics is the saturation scale Q s and its product on characteristic size of the dipoles r in the process.The saturation has a mild dependence on energy [61] and in the kinematics of interest for our studies ( √ s ∈ (1.9, 7) TeV) its values are Q 2 s = 0.7 − 0.9 GeV 2 .The typical size of the dipole essential in our process might be estimated as 5 so the product τ ≡ r 2 Q 2 s ≈ 0.5 ... 0.7.Contrary to the large-m c expectations, this number is not very small and corresponds to the dynamics in the vicinity of the saturation scale [37,63].The scattering amplitude is well known in this region (see (15) below), so this implies minimal model-dependence in our estimates.

III. INTERRELATION WITH THE DIFFRACTIVE J/ψ PRODUCTION IN DIS
The cross-section of the J/ψ production in the small-x kinematics is closely related to the diffractive J/Ψ production in DIS, and this relationship is useful to fix the unknown non-perturbative factor ∼ d 2 Q T S 2 h (Q T ) which appears in (7,8).In the BFKL picture, there are elastic and inelastic contributions to diffractive J/ψ production, shown in the left and right panels of the Figure 3. Taking into account approximate equality of the elastic and inelastic contributions [64,65], we may focus on evaluation of the diagram (a) of the Fig. (3).For fixing the prefactor , it is sufficient to consider only the q T -integrated cross-section dσ diff /dy, which is given by (see Fig. (3)): In the vicinity of the saturation scale, which gives the dominant contribution, the CGC/saturation approach predicts that the amplitude N has a form where γ ≈ 0.63 [59], x ≈ M J/ψ e −y / √ s, and Q 2 s (x) is the saturation scale.The Eq. ( 14) in this case takes the form , which allows to factorize all the energy-and rapidity dependence, and in this way facilitate scaling from HERA to Tevatron and LHC energies.Similarly, for the gluon induced elastic J/ψ production in the vicinity of the saturation scale (see Ref. [37]) we may simplify (9) to 5 See Appendix A and [62] for paramertizations of the wave functions which also has a factorizable dependence on energy and rapidity, Taking into account that the wave functions of the proton and gluon are proportional to each other in the leading order of pQCD [35], we may expect the proportionality of the diffractive and inclusive cross-sections where and for estimates of parameter R we used the parameterizations of the wave functions given in Appendix A. The data on diffractive production are available from H1 and ZEUS experiments [64][65][66], and they show that in HERA kinematics the elastic and inelastic contributions (diagrams a and b in the Fig. 3, respectively) are approximately equal.As a consequence, for W ≡ √ s ≈ 30 GeV (x g ≈ 0.01) the diffractive cross-section This allows to to use Eq. ( 21) for estimate of the hadroproduction cross-section at √ s ≈ 30 GeV, In the next section we will extrapolate it with the help of the small-x evolution up to LHC energies.Table I.Comparison of the theoretical estimates with experimental data for the cross-section dσ/dy at central rapidities.
Theoretical estimates correspond to values of parameter λ ∈ (0.2, 0.3) (lower and upper values respectively).In the last column we've quoted the data on prompt J/ψ production.The cross-section dσ/dy at Tevatron was extracted dividing the total cross-section σtot(|y| < 0.6) by the width of the bin ∆y ∼ 1.2.

IV. PHENOMENOLOGICAL ESTIMATES
A. Total J/Ψ production cross section As we demonstrated in Section II, the typical size r of Q Q dipole is small when the saturation scale Q s (x) m c (see Eq. ( 13)), and for this reason from (18) we expect that the cross-section at central rapidities should scale with energy as For numerical estimates we assume that the saturation scale Q 2 s scales as [63] with λ ≈ 0.2 − 0.3 [70][71][72] and x ≈ M J/ψ e −y / √ s.In case of Tevatron and LHC kinematics, this leads to the cross-section estimates given in the Table I.In these estimates we use that the evolution does not change the ratio between the elastic and inelastic contributions which correspond to two terms in the solution to the evolution equation for ρ (2) (see Appendix B and Fig. ( 9)).
As we can see, the suggested mechanism gives a significant contribution to the total cross-section, though in view of inherent uncertainties of the CGC/Saturation approach we cannot make more accurate estimate about the fraction of charmonia produced via this mechanism.
As we have mentioned earlier, our estimates exceed the result of Ref. [38], where similar contribution was found approximately twice smaller.Below we would like to analyze the reasons which might be responsible for this discrepancy.It was assumed in [38] that the digluon distribution is proportional to a direct product of independent gluon distributions, with a phenomenological multiplicative pre-factor usually described by the so-called effective crosssection σ eff .In the CGC/Saturation picture, each gluon effectively is replaced by a pomeron, so the factorized model would correspond to the diagram (a) in the Fig. (3).As we discussed earlier, the inelastic contributions (diagram (b) in the same Figure) yields a numerically comparable contribution which would correspond to corrections breaking the factorized form.This additional inelastic contribution might be one of the reasons of the strong channel dependence of the phenomenological effective cross-section σ eff 6 .In our approach, we fix the the only unknown constant d 2 Q T S 2 (Q T ) from experimental data at HERA and evolve the cross-section as prescribed by small-x evolution, thus taking into account both diagrams of the Fig. (3).The main source of uncertainty in this procedure is the energy dependence of the saturation scale Q 2 s (x), namely the choice of coefficient λ in (27).Theoretical estimates of this parameter significantly depend on the choice of the scale in the running coupling constant ᾱs , and lead to estimates λ ≈ 0.3 − 0.4 [73], whereas phenomenological estimates of this parameter restrict it to the range λ ∈ (0.2, 0.3) [59].We can see that it results in up to fifty per cent uncertainty in the theoretical estimates shown in the Table I.

B. Transverse momentum distribution
Qualitatively the main features of q T -distribution may be understood from (8).For small momenta 1 2 q T Q s (x), we can safely use the behaviour of the Pomeron Green function in the vicinity of the saturation momentum [63] In the kinematic region 1 2 q T k T the scattering amplitude becomes sensitive to dynamics at shorter distances, which is described in perturbative QCD.In this region we may use a parameterization (15), which is a BFKL prediction for small dipoles with r Q −1 s (x) [59].The two approaches are valid for different values of q T and thus complement each other.The choice of the threshold scale q 0 between them is somewhat arbitrary, yet we expect that q 0 should be of order 0.001 0.010 0.100

ATLAS
. The shape of the qT distribution of produced J/ψ at central rapidities.The threshold scale q0 is chosen as q0 ∼5 GeV.The experimental data are taken from Refs.[74,75].
In Fig. (4) we compare model predictions for the shape of the q T -dependence with experimental data from [74,75].In order to avoid the above-mentioned global uncertainty in normalization, we consider the normalized ratio d 2 σ/dy dq T / (dσ/dy).The cusp on the curve near the threshold scale q T = q 0 can be smoothed out by more relaxed conditions.

C. Rapidity distribution
For numerical estimates in the previous sections we used the leading log approximation (LLA) for the BFKL Pomeron Green function, assuming additionally that mass m c is large, so that we could use a small-r approximation.The shape of the rapidity distribution in this approximation has a very simple form, which follows directly from Eq. ( 18): since all dependence on rapidity is concentrated in the energy dependence of the saturation scale (27) which contributes to the cross-section (18) multiplicatively, and the gluon density in prefactor.The latter in the small-x kinematics is expected to have a simple power-like behaviour However, such simple parameterization is valid only at central rapidities (near y = 0), and as could be seen from the Figure 5, outside this region (|y| 1) mismatches the experimental cross-section even on the qualitative level.This happens because at very forward or very backward kinematics we should add in (29) additional factor x = e −y * ±y per each gluon/pomeron in order to have correct endpoint behaviour of the gluon densities in the x → 1 limit [76].
As we can see from the Figure 5, the model gives a reasonable description of the rapidity dependence, which might be interpreted as confirmation of the applicability of the small-r approximation.The ratio R J/ψ (y) is defined as R J/ψ = dσ J/ψ /dy / dσ J/ψ /dy y=0 .The dashed curve corresponds to Eq. ( 29) with BFKL-style parameterization 31 for gluon densities, the solid line takes into account ∼ (1 − x) 5 -endpoint factors as described in the text.The data are taken from Ref. [75].

D. Multiplicity dependence
The J/ψ production accompanied by two parton showers occurs in the events with larger multiplicity than the average multiplicity n in the inclusive hadron production [24,28].The mechanism suggested in this paper provides a natural explanation of this enhancement.The cross-section includes contributions of additional parton shower as shown in  6)-b, which enhances the observed multiplicity.Technically, the dependence on multiplicity n of the produced particles affects the cross-section through the increase of the number of the particles in a parton shower and change of the value of the saturation scale, which is proportional to the density of produced gluons (see.[39,77] for more details).Assuming that for hadron-hadron collisions the area of interaction does not depend on multiplicity, the saturation scale is linearly proportional to the number of particles in the shower [78], where n is the average multiplicity at y * = 0 (n ≈6.5 at W = 13 TeV [79]).As we discussed earlier, the product of dipole size on saturation scale r 2 Q 2 (y, n) is close to one, for this reason the multiplicity events with n/n 1 probe a deep saturation regime, where approximation (18) is not valid, and we should use (9) [24].Solid line: result of evaluations with (9) using CGC parameterization for the dipole amplitude and saturation scale adjusted according to (32).The upper band marked "2-showers" stands for the estimates with approximate expression (33) and values of γ varied in the range γ ∈ (0.67, 0.76), as implemented in phenomenological parameterizations.Similarly, the lower color band marked with label "1-shower" stands for multiplicity distribution evaluated in single-shower mechanism shown in the Fig.
(1)-b.Right plot: large-multiplicity extension of the solid curve from the left plot (result of evaluations with (9) using CGC parameterization).
In the Figure 7 we plot the self-normalized multiplicity distribution evaluated with (9), using CGC parameterization for the dipole amplitude and saturation scale adjusted according to (32).We can see that agreement with experimental data from ALICE [24] is reasonable.At sufficiently small multiplicities, dN J/ψ /dy is increasing as a function of dN ch /dη, however, as can be seen from the right plot, at larger multiplicities in deep saturation regime it starts decreasing.This behaviour might be understood from the structure of the last line in (9): the dipole amplitude N saturates (approaches asymptotically a constant) [80], for this reason the difference N y; r+ r 2 , 0 − N y; r− r 2 , 0 2 gets suppressed.For the sake of completeness in the left panel of the Figure 7 we also plotted the ratio evaluated with the simplified parameterization (15).Taking into account that contributions of left and right diagrams in the Figure (6) where the coefficient κ = Q 2 s (Y − y) /Q 2 s (y) γ ≈ e γ λ(Y −2y) reflects the relative suppression of the contribution of the right diagram in the Fig. (6) with respect to the left diagram for different rapidities (at central rapidities used for comparison with data κ ≈ 1).While the LO CGC/Saturation predicts a value γ ≈ 0.63 , the phenomenological fits [70,72] favor slightly higher values of γ, for this reason we varied this parameter in the range γ ∈ (0.67, 0.76) .We also plotted the estimates of single-shower mechanism shown in the right part of the Figure (1).Within model uncertainty, we can see that the experimental data support the main hypothesis of the paper that the J/ψ production accompanied by the production of two parton showers gives a significant contribution.

V. CONCLUSIONS
In this paper we analyzed in the CGC/Saturation framework the J/ψ hadroproduction, accompanied by production of two parton showers.We demonstrated that this mechanism gives a large contribution to the total cross section of J/ψ production, as well as to the transverse momenta and multiplicity distribution of produced J/ψ.The experimental data [24] on multiplicity distributions of produced J/ψ seem to favor this hypothesis (compared to conventional mechanisms estimated in the CGC/Saturation approach).This mechanism has no suppression of the order of ᾱS at high energies compared to one parton shower contribution.We show that it alone can describe the shape of the p Tand rapidity dependence of the cross section.As we commented in detail in section 5, our results are approximately twice larger than similar evaluation of Ref. [38] and for this reason give a significant contribution to the experimental data.However, inherent uncertainties of the CGC/Saturation approach preclude more precise estimates of its fraction.Since the main contribution to the cross section of this process stems from the vicinity of the saturation scale, the behaviour of the scattering amplitude in this region is largely under theoretical control, and thus has minor dependence on phenomenological parameterizations of the dipole amplitude.
We believe that our studies will bring attention to the contribution of multiparton densities in production of J/ψ.

CONTENTS I. Introduction 2 II. 7 IV
Charmonia production in the BFKL approach 3 III.Interrelation with the diffractive J/ψ production in DIS of the factor I p 1,T , p 2,T 13 B. Evolution equation for the double gluon density 14 References 15 I. INTRODUCTION

Figure 1 .
Figure 1.Two parton showers contribution to J/Ψ production in hadron-hadron collisions Fig. (1-a) and production of J/Ψ and even number of soft (non-perturbative) gluons from one parton shower.
the framework of the BFKL Pomeron calculus, the cross-section of the process shown in the diagram (a) of the Fig. (1) is described by the exchange of two BFKL Pomerons as shown in Fig. (2).Since we are interested in inelastic J/ψ production, the Pomerons in Fig. (2) are cut Pomerons in which all gluons are produced.From the unitarity constraints for the elastic amplitude N BFKL of the dipole of size r, rapidity Y and at the impact parameter b, we have[59]

Figure 2 .
Figure 2. The cross-section corresponding to the first diagram of the Fig. (1) in the BFKL Pomeron calculus..The vertical wavy lines of different colors and shape passing through unitarity cut are BFKL Pomerons (described by the Green functions GIP (y, kT ) in Eq. (2)), helical lines correspond to the gluons.

Figure 3 .
Figure 3. Two main contributions to the diffractive production of J/Ψ meson.The elastic contribution (diagram (a)) contributes mainly to production of hadrons with small total mass, while the inelastic contribution (diagram (b)) is the source of hadrons with large total mass.

Figure 6 .
Figure 6.Multiplicities in the J/ψ production, accompanied by two parton showers.