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

TL;DR

This paper investigates the dominant three-pomeron fusion mechanism for inclusive $J/\psi$ production in proton-proton collisions, showing it aligns with CGC/Saturation theory and explains experimental multiplicity enhancements.

Contribution

It introduces the three-pomeron fusion as a key mechanism for $J/\psi$ production, connecting it with CGC/Saturation and experimental multiplicity observations.

Findings

Three-pomeron fusion significantly contributes to $J/\psi$ production.

The mechanism explains the shapes of rapidity, momentum, and multiplicity distributions.

The approach accounts for the observed multiplicity enhancement in experiments.

Abstract

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.

Figures

Tables

Table 1. Table 1: Comparison of the theoretical estimates with experimental data for the cross-section d ​ σ / d ​ y 𝑑 𝜎 𝑑 𝑦 d\sigma/dy at central rapidities. Theoretical estimates correspond to values of parameter λ ∈ ( 0.2 ,  0.3 ) 𝜆 0.2 0.3 \lambda\in(0.2,\,0.3) (lower and upper values respectively). In the last column we’ve quoted the data on prompt J / ψ 𝐽 𝜓 J/\psi production. The cross-section d ​ σ / d ​ y 𝑑 𝜎 𝑑 𝑦 d\sigma/dy at Tevatron was extracted dividing the total cross-section σ tot ​ ( | y | < 0.6 ) subscript 𝜎 tot 𝑦 0.6 \sigma_{{\rm tot}}(|y|<0.6) by the width of the bin Δ ​ y ∼ 1.2 similar-to Δ 𝑦 1.2 \Delta y\sim 1.2 .

	Theoretical estimates	Experiment
$\sqrt{s}\approx 1.96$ TeV	2.1-2.6 $\mu$b	2.38  $\mu$b PSICDF
$\sqrt{s}\approx 7$ TeV	3.8-5.6 $\mu$b	5.8 $\mu$b Abelev:2012gx

Equations

$$N_{\mathrm{cut}}^{\mathrm{BFKL}}\left(Y,r,b\right)\equiv 2N^{\mathrm{BFKL}}\left(Y,r,b\right)$$

$$\frac{d^2\sigma \left(Y,q_T\right)}{dy{d}^2{q}_T}=\frac{4C_F^3{\bar{\alpha }}_S^3}{(2\pi {)}^6}\int d^2{k}_T{d}^2{p}_T{d}^2{Q}_T{G}_{IP}^{\mathrm{cut}}\left(Y-y,{\mathbf{p}}_T,0\right)\times$$

$$I\left(\boldsymbol{k}_{T},\boldsymbol{q}_{T}\right)\,\,=\,\,\int_{0}^{1}dz\int d^{2}r\,\Big{(}e^{i\frac{1}{2}\boldsymbol{q}_{T}\cdot\boldsymbol{r}}\,\,-\,e^{i\boldsymbol{k}_{T}\cdot\boldsymbol{r}}\Big{)}\,\,\Psi_{g}\left(p_{T},r;z\right)\,\Psi_{J/\psi}\left(r;\,z\right)\,\,=\,\,F\left(\frac{1}{2}\boldsymbol{q}_{T}\right)\,-\,F\left(\boldsymbol{k}_{T}\right)$$

$$N\left(r,r^{\prime };b\right)={\left(\frac{r^2{r}^{\prime 2}}{{\left(\vec{b}-\frac{1}{2}\left(\vec{r}-{\vec{r}}^{\prime }\right)\right)}^2{\left(\vec{b}-\frac{1}{2}\left(\vec{r}-{\vec{r}}^{\prime }\right)\right)}^2}\right)}^{\gamma }$$

$$G_{IP}^{\mathrm{cut}}\left(y,{\mathbf{k}}_T\pm \frac{1}{2}{\mathbf{q}}_T,{\mathbf{Q}}_T\right)\approx 2G_{IP}^{\mathrm{BFKL}}\left(y,{\mathbf{k}}_T\pm \frac{1}{2}{\mathbf{q}}_T,{\mathbf{Q}}_T=0\right)S_h\left(Q_T\right)$$

$$\frac{d\sigma \left(Y,Q^2\right)}{dy{d}^2{q}_T}$$

$$\frac{d\sigma \left(Y,Q^2\right)}{dy{d}^2{q}_T}=\frac{4C_F^3{\bar{\alpha }}_S^2}{(2\pi {)}^4}\int \frac{d^2{Q}_T}{(2\pi {)}^2}{S}_h^2\left(Q_T\right)x_{g}G\left(x_g,M_{J/\psi }\right)$$

$$\times \int d^2{p}_{1,T}{d}^2{p}_{1,T}{\delta }^{(2)}\left({\mathbf{p}}_{1,T}+{\mathbf{p}}_{2,T}-{\mathbf{q}}_T\right)I^2\left({\mathbf{p}}_{1,T},{\mathbf{p}}_{2,T}\right)G_{IP}^{\mathrm{BFKL}}\left(y;{\mathbf{p}}_{1,T},0\right)G_{IP}^{\mathrm{BFKL}}\left(y,{\mathbf{p}}_{2,T},0\right)$$

$$\frac{d\sigma \left(Y,Q^2\right)}{dy{d}^2{q}_T}=4\int \frac{d^2{Q}_T}{(2\pi {)}^2}{S}_h^2\left(Q_T\right)x_{g}G\left(x_g,M_{J/\psi }\right)$$

$$\times {\int }_0^{1}dz{\int }_0^{1}d{z}^{\prime }\int \frac{d^{2}r}{4\pi }\frac{d^2{r}^{\prime }}{4\pi }{d}^{2}b{e}^{-i{\mathbf{q}}_T\cdot \left(\mathbf{b}+\frac{1}{2}(\mathbf{r}-{\mathbf{r}}^{\prime })\right)}⟨{\Psi }_g\left(r,z\right){\Psi }_{J/\psi }\left(r,z\right)⟩⟨{\Psi }_g\left(r^{\prime },z^{\prime }\right){\Psi }_{J/\psi }\left(r^{\prime },z^{\prime }\right)⟩$$

$$\displaystyle\times\,\,\Bigg{(}N\left(y;\,\boldsymbol{b}\right)\,+\,N\left(y;\,\boldsymbol{b}+\boldsymbol{r}-\boldsymbol{r}^{\prime}\right)\,-\,N\left(y;\,\boldsymbol{b}+\boldsymbol{r}\right)\,-\,N\left(y;\boldsymbol{b}-\boldsymbol{r}^{\prime}\right)\Bigg{)}^{2}$$

$$=4\int \frac{d^2{Q}_T}{(2\pi {)}^2}{S}_h^2\left(Q_T\right)x_{g}G\left(x_g,M_{J/\psi }\right)$$

$$\times {\int }_0^{1}dz{\int }_0^{1}d{z}^{\prime }\int \frac{d^{2}r}{4\pi }\frac{d^2{r}^{\prime }}{4\pi }{d}^{2}b{e}^{-i{\mathbf{q}}_T\cdot \left(\mathbf{b}\right)}⟨{\Psi }_g\left(r,z\right){\Psi }_{J/\psi }\left(r,z\right)⟩⟨{\Psi }_g\left(r^{\prime },z^{\prime }\right){\Psi }_{J/\psi }\left(r^{\prime },z^{\prime }\right)⟩$$

$$\displaystyle\times\,\,\,\Bigg{(}N\left(y;b-\frac{1}{2}\left(\boldsymbol{r}-\boldsymbol{r}^{\prime}\right)\right)\,+\,N\left(y;\boldsymbol{b}+\frac{1}{2}\left(\boldsymbol{r}-\boldsymbol{r}^{\prime}\right)\right)\,$$

$$\displaystyle\,\,-\,\,N\left(y;\boldsymbol{b}+\frac{1}{2}\left(\boldsymbol{r}+\boldsymbol{r}^{\prime}\right)\right)\,\,-\,\,N\left(y;\boldsymbol{b}-\frac{1}{2}\left(\boldsymbol{r}+\boldsymbol{r}^{\prime}\right)\right)\Bigg{)}^{2}$$

$$\frac{d\sigma \left(Y,Q^2\right)}{dy}=16\int \frac{d^2{Q}_T}{(2\pi {)}^2}{S}_h^2\left(Q_T\right)x_{g}G\left(x_g,m_c\right)$$

$$\times {\int }_0^{1}dz{\int }_0^{1}d{z}^{\prime }\int \frac{d^{2}r}{4\pi }\frac{d^2{r}^{\prime }}{4\pi }⟨{\Psi }_g\left(r,z\right){\Psi }_{J/\psi }\left(r,z\right)⟩⟨{\Psi }_g\left(r^{\prime },z^{\prime }\right){\Psi }_{J/\psi }\left(r^{\prime },z^{\prime }\right)⟩$$

$$\displaystyle\times\,\,\Bigg{(}\,N\left(y;\,\frac{\vec{r}+\vec{r}^{\prime}}{2}\right)\,-\,N\left(y;\,\frac{\vec{r}-\vec{r}^{\prime}}{2}\right)\Bigg{)}^{2}$$

$$\frac{d\sigma \left(Y,Q^2\right)}{dy{d}^2{q}_T}$$

$${\rho }^{(2)}\left(x_1,{\mathbf{p}}_{1,T};x_2,{\mathbf{p}}_{2,T},{\mathbf{Q}}_T\right)=$$

$$\displaystyle\,\left\langle P\left|\Bigg{\{}a^{+}\left(x_{1},\boldsymbol{p}_{1,T}+\frac{1}{2}\boldsymbol{Q}_{T};b\right)\,a^{+}\left(x_{2},\boldsymbol{p}_{2,T}-\frac{1}{2}\boldsymbol{Q}_{T};c\right)\,a\left(x_{2},\boldsymbol{p}_{2,T}+\frac{1}{2}\boldsymbol{Q}_{T};c\right)\,a\left(x_{1},\boldsymbol{p}_{1,T}-\frac{1}{2}\boldsymbol{Q}_{T};b\right)\Bigg{\}}\right|P\right\rangle$$

$$\frac{d\sigma \left(Y,Q^2\right)}{dy{d}^2{q}_T}=4\int \frac{d^2{Q}_T}{(2\pi {)}^2}{x}_{g}G\left(x_g,2m_c\right)$$

$$\times {\int }_0^{1}dz{\int }_0^{1}d{z}^{\prime }\int \frac{d^{2}r}{4\pi }\frac{d^2{r}^{\prime }}{4\pi }{d}^{2}b{e}^{-i{\mathbf{q}}_T\cdot \mathbf{b}},⟨{\Psi }_g\left(r,z\right){\Psi }_{J/\psi }\left(r,z\right)⟩⟨{\Psi }_g\left(r^{\prime },z\right){\Psi }_{J/\psi }\left(r^{\prime },z\right)⟩$$

$$\displaystyle\times\,\,\Bigg{(}\rho^{(2)}\left(y;\boldsymbol{b}-\boldsymbol{r}_{d};y;\boldsymbol{b}-\boldsymbol{r}_{d},Q_{T}\right)\,+\,\rho^{(2)}\left(y;\boldsymbol{b}+\boldsymbol{r}_{d},y;\boldsymbol{b}+\boldsymbol{r}_{d},Q_{T}\right)\,+\,\rho^{(2)}\left(y;\boldsymbol{b}+\boldsymbol{r}_{s},y;\boldsymbol{b}+\boldsymbol{r}_{s},Q_{T}\right)$$

$$+{\rho }^{(2)}\left(y;\mathbf{b}-{\mathbf{r}}_s,y;\mathbf{b}-{\mathbf{r}}_s,Q_T\right)+2{\rho }^{(2)}\left(y;\mathbf{b}-{\mathbf{r}}_d,y;\mathbf{b}+{\mathbf{r}}_s,Q_T\right)-2{\rho }^{(2)}\left(y;\mathbf{b}+{\mathbf{r}}_d,y;\mathbf{b}+{\mathbf{r}}_s,Q_T\right)$$

$$-{\rho }^{(2)}\left(y;\mathbf{b}-{\mathbf{r}}_d,y;\mathbf{b}-{\mathbf{r}}_s,Q_T\right)-2{\rho }^{(2)}\left(y;\mathbf{b}+{\mathbf{r}}_d,y;\mathbf{b}+{\mathbf{r}}_s,Q_T\right)-2{\rho }^{(2)}\left(y;\mathbf{b}+{\mathbf{r}}_d,y;\mathbf{b}-{\mathbf{r}}_s,Q_T\right)$$

$$\displaystyle+\,2\,\rho^{(2)}\left(y;\boldsymbol{b}+\boldsymbol{r}_{s},y;\boldsymbol{b}-\boldsymbol{r}_{s},Q_{T}\right)\Bigg{)}$$

$$⟨r^2⟩=\frac{{\int }_0^{1}dz\int \frac{d^{2}r}{4\pi }{r}^2⟨{\Psi }_g\left(r,z\right){\Psi }_{J/\psi }\left(r,z\right)⟩}{{\int }_0^{1}dz\int \frac{d^{2}r}{4\pi }⟨{\Psi }_g\left(r,z\right){\Psi }_{J/\psi }\left(r,z\right)⟩}\approx 0.76{\mathrm{GeV}}^{-2}$$

$$\frac{d{\sigma }_{\mathrm{diff}}\left(y,Q^2,\sqrt{s_{{\gamma }^{\ast }p}}\right)}{dy}=\int \frac{d^2{Q}_T}{(2\pi {)}^2}{S}_h^2\left(Q_T\right){\int }_0^{1}dz{\int }_0^{1}d{z}^{\prime }\int \frac{d^{2}r}{4\pi }\frac{d^2{r}^{\prime }}{4\pi }\times$$

$$N\left(y,r^2,0\right)\approx {\mathrm{Const}}_s{\left(r^2{Q}_s^2(x)\right)}^{\bar{\gamma }},$$

$$\frac{d{\sigma }_{\mathrm{diff}}\left(y,Q^2,\sqrt{s_{{\gamma }^{\ast }p}}\right)}{dy}=\text{mbox}{Const}_{\mathrm{s}}^2\int \frac{d^2{Q}_T}{(2\pi {)}^2}{S}_h^2\left(Q_T\right){\int }_0^{1}dz{\int }_0^{1}d{z}^{\prime }\int \frac{d^{2}r}{4\pi }\frac{d^2{r}^{\prime }}{4\pi }$$

$$\frac{d{\sigma }_{\mathrm{diff}}\left(y,Q^2,\sqrt{s_{{\gamma }^{\ast }p}}\right)}{dy}\sim {\left[Q_s^2\left(x\right)\right]}^{2\bar{\gamma }}$$

$$\frac{d\sigma \left(y,Q^2,\sqrt{s_{{\gamma }^{\ast }p}}\right)}{dy}=16\text{mbox}{Const}_{\mathrm{s}}^2\int \frac{d^2{Q}_T}{(2\pi {)}^2}{S}_h^2\left(Q_T\right)x_{g}G\left(x_g,m_{J/\Psi }\right)$$

$$\times C_F^2{\int }_0^{1}dz{\int }_0^{1}d{z}^{\prime }\int \frac{d^{2}r}{4\pi }\frac{d^2{r}^{\prime }}{4\pi }⟨{\Psi }_g\left(r,z\right){\Psi }_{J/\psi }\left(r,z\right)⟩⟨{\Psi }_g\left(r^{\prime },z^{\prime }\right){\Psi }_{J/\psi }\left(r^{\prime },z^{\prime }\right)⟩$$

$$\displaystyle\,\Bigg{\{}\left(\int d^{2}b\,Q_{s}^{2}\left(x,\,b\right)\right)^{2\bar{\gamma}}\left(\left(\frac{\vec{r}+\vec{r}^{\prime}}{2}\right)^{2\,\bar{\gamma}}\,-\,\left(\frac{\vec{r}-\vec{r}^{\prime}}{2}\right)^{2\,\bar{\gamma}}\right)^{2}\Bigg{\}}$$

$$\frac{d\sigma \left(y,Q^2\right)}{dy}\sim x_{g}G\left(x_g{M}_{J/\psi }\right){\left[Q_s^2\left(x\right)\right]}^{2\bar{\gamma }}.$$

$${\Psi }_{{\gamma }^{\ast }}\left(r,z\right)=\frac{\sqrt{2N_c}}{g}\frac{2}{3}e{N}_c{\Psi }_g\left(r,z\right)$$

$$\frac{d\sigma\left(y,\,Q^{2}\right)/dy}{d\sigma_{{\rm diff}}\left(y,\,Q^{2}\right)/dy}\Bigg{|}_{Q=m_{J/\psi},y\to Y_{0}}\,\,\approx\frac{9}{2}\,x_{g}G\left(x_{g},\,M_{J/\psi}\right)\,\mathcal{R}\,C_{F}^{2}\,\frac{\alpha_{S}}{\alpha_{{\rm e.m,.}}N_{c}^{3}}$$

$$\mathcal{R}$$

Full text

DGLAP and BFKL evolution, double parton distributions, Bose-Einstein correlations, shadowing corrections, non-linear evolution equation, CGC approach.

††preprint: TAUP - 3030/18, USM-TH-359

1. I Introduction

2. II Charmonia production in the BFKL approach

3. III Interrelation with the diffractive $J/\psi$ production in DIS

4. IV Phenomenological estimates

5. IV.1 Total $J/\Psi$ production cross section

6. IV.2 Transverse momentum distribution

7. IV.3 Rapidity distribution

8. IV.4 Multiplicity dependence

9. V Conclusions

10. VI Acknowledgements

11. A Evaluation of the factor $I\left(\boldsymbol{p}_{1,T},\boldsymbol{p}_{2,T}\right)$

12. B Evolution equation for the double gluon density

I Introduction

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 Korner:1991kf . However, in reality the coupling $\alpha_{s}\left(m_{c}\right)\sim 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)  Chang:1979nn ; Baier:1981uk ; Berger:1980ni assumes that the dominant mechanism of charmonia production is the gluon-gluon 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 $\alpha_{s}$ due to milder suppression of higher order terms at large $p_{T}$ Brodsky:2009cf ; Artoisenet:2008fc and co-production of additional quark pairs Artoisenet:2007xi ; Karpishkov:2017kph motivated introduction of the phenomenological Color Octet contributions Cho:1995vh ; Cho:1995ce . The modern NRQCD formulation Bodwin:1994jh ; Maltoni:1997pt ; Brambilla:2008zg ; Feng:2015cba ; Brambilla:2010cs 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 Feng:2015cba ; Baranov:2016clx ; Baranov:2015laa , 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/\psi$ or co-production together with other hadrons Bodwin:2012xc ; Bodwin:2014gia ; Braaten:1994xb ; Braaten:1995cj . The latter findings are partially supported by experimental data on multiplicity of co-produced charged particles PSIMULT ; Porteboeuf:2010dw ; Sjorstrand:1987 ; Bartalini:2010 ; Alice:2012Mult which suggest that $J/\psi$ production might get sizeable contribution from this mechanism.

In this paper we analyze $J/\psi$ 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$ Kovchegov:2012mbw ; KOLEB . We focus on the mechanism of $J/\psi$ production shown in the diagram $(a)$ of the Fig. (1), when the production of $J/\psi$ occurs from fusion of three gluons, accompanied by the production of two parton showers 111We would like to clarify what we mean, saying “the production of two parton showers.” The parton shower is the initial state “cut ladder,” which generates the production of gluons that are almost uniformly distributed in the entire region of rapidities. As we will see below, such parton showers correspond to the cut BFKL Pomeron. . In order to emphasize the role of co-production of other particles associated with $J/\psi$, 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 $\bar{n}$ of the gluon production in the inclusive process. This mechanism is similar to mechanisms studied earlier in the literature for proton-proton collisions KMRS ; Schmidt:2018gep and for proton-nucleus and nucleus-nucleus scattering KHTU ; Schmidt:2018rkw ; Kopeliovich:2017jpy . For collisions involving heavy nuclei, the two parton shower mechanism inside heavy nuclei with atomic number $A$ is dominant due to factor $\sim A^{1/3}$, and, because of this, has been comprehensively discussed during the past decade KLNT ; DKLMT ; KMV ; GOLEPSI . On the other hand, for proton-proton scattering we found only one paper which considers this process in the $k_{T}$-factorization approach MOSA .

At first sight, the two parton showers mechanism should be suppressed in comparison with the production of $J/\psi$ in one parton shower shown in the diagram $(b)$ of the Fig. (1-b): in contrast to $\alpha_{s}$ coupling which appears in the shower-quark vertex, the emission of soft gluon in Fig. (1-b) is not suppressed by any hard scale and thus does not bring a smallness proportional to $\bar{\alpha}_{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 KOLEB for review). For typical $\langle r\rangle\,\sim\,1\,{\rm GeV}^{-1}$ we get $\bar{\alpha}_{S}\left(4/r^{2}\right)\approx 0.2$ so $\bar{\alpha}_{S}G^{{\rm BFKL}}\left(s,\dots\right)\,\propto\bar{\alpha}_{S}\,s^{\Delta_{\small{\rm BFKL}}}\geq 1$ at high energies, where $\Delta_{{\rm BFKL}}=4\ln 2\bar{\alpha}_{S}$ is the intercept of the BFKL PomeronBFKL and $G^{{\rm BFKL}}\left(s,\dots\right)$ is the Green function of the BFKL Pomeron. The numerical estimates of Ref.MOSA 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 $\sigma_{{\rm 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 $\alpha_{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. Kovchegov:2012mbw 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. MOSA , we use a CGC/Saturation framework, and in order to avoid uncertainties related to digluon distributions, we relate the cross-sections of the $J/\psi$ 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/\psi$ in DIS. In the Section IV we make phenomenological estimates and compare results with experimental data. The subsection IV.1 is dedicated to the estimates of the total cross section of the process. In this subsection we point out the differences with Ref. MOSA . In subsections IV.2, IV.3 and IV.4 we consider the momenta, rapidity and multiplicity distributions of the differential cross-sections. Finally, in the Section V we draw conclusions.

II Charmonia production in the BFKL approach

At present, the effective theory of QCD at high energies exists in two forms: the CGC/saturation approach GLR ; MUQI ; MV ; MUCD ; BK ; JIMWLK and the BFKL Pomeron calculus BFKL ; LIREV ; GLR ; MUPA ; BART ; BRN ; KOLE ; LELUR ; LMP . It has been proven that in general these two approaches are equivalent in a limited region of the rapidities AKLL . The interpretation of processes at high energy appears quite different in each approach, since they have different structural elements.

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/\psi$ production in this approach may be found in KMV 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 COL . This approach is suitable for describing diffractive physics and correlations in multi-particle production, so we can use the Mueller diagram technique MUDI . 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 AGK , which are useful in spite of the restricted region of their application AGKQCD . In this paper for the sake of definiteness we use the BFKL Pomeron calculus as a framework for evaluations.

In 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/\psi$ production, the Pomerons in Fig. (2) are cut Pomerons in which all gluons are produced. From the unitarity constraints for the elastic amplitude $N^{{\rm BFKL}}$ of the dipole of size $r$, rapidity $Y$ and at the impact parameter $b$, we have Kovchegov:2012mbw

[TABLE]

Its contribution to the total cross section for $J/\psi$ production222The color coefficient has been calculated in Ref.KLNT and at large $N_{c}$ it turns out to be equal $C_{F}^{3}$. is equal to

[TABLE]

where $G_{I\!\!P}^{{\rm cut}}$ are the Green functions of the cut pomerons (it is related to elastic amplitudes $N_{{\rm cut}}^{{\rm BFKL}}$ by Fourier transform), the gluon momenta $\boldsymbol{k}_{T},\,\boldsymbol{p}_{T},\,\boldsymbol{q}_{T}$, are defined in the Fig. (2), $\boldsymbol{Q}_{T}\,=\,\boldsymbol{k}_{T}\,-\,\boldsymbol{k}^{\prime}_{T}$, and

[TABLE]

(the evaluation of the factor $I\left(\boldsymbol{k}_{T},\boldsymbol{q}_{T}\right)$ 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 AGK and 2 from the fact that gluon that produces $c\bar{c}$-pair can come from another proton. In Eq. (3) $\Psi_{g}\left(p_{T},\,r\right)$ 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 $\Psi_{J/\psi}$ is the wave function of $J/\psi$ meson. The amplitude of the BFKL Pomeron $G_{{I\!\!P}}^{{\rm cut}}\left(y,\,\boldsymbol{k}_{T}\pm\frac{1}{2}\boldsymbol{q}_{T},\,\boldsymbol{Q}_{T}\right)$ can be simplified if we take into account that the $\boldsymbol{Q}_{T}$ dependence of the BFKL Pomeron is determined by the size of the largest of the interacting dipoles 333The 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 toLIREV

$N\left(r,r^{\prime};b\right)\,\,=\,\,\left(\frac{r^{2}\,\,r^{\prime 2}}{\left(\vec{b}-\frac{1}{2}\left(\vec{r}-\vec{r}^{\prime}\right)\right)^{2}\,\left(\vec{b}-\frac{1}{2}\left(\vec{r}-\vec{r}^{\prime}\right)\right)^{2}}\right)^{\gamma}$

(4)

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^{\prime}$. In our process $r^{\prime}$ 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 $m_{c}$, or the saturation scale $Q_{s}$ and, therefore, turns out to be much larger than $1/R_{h}$, and can be neglected. The dependence on $Q_{T}\approx 1/R_{h}$ is described by $S_{h}\left(Q_{T}\right)$ in Eq. (5), which has the non-perturbative origin and. in practice, has to be taken from the experiment., and accounting for Eq. (1), may be written as

[TABLE]

The dependence on $Q_{T}$ is described by $S_{h}\left(Q_{T}\right)$ 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. (II) in the form

[TABLE]

For further evaluations it is very convenient to introduce momenta $\boldsymbol{p}_{1,2,T}=\pm\boldsymbol{k}_{T}+\frac{1}{2}\boldsymbol{q}_{T}$, which allow to rewrite Eq. (II) as

[TABLE]

where we took the integral over $\boldsymbol{p}_{T}\in\left(0,\,2m_{c}\right)$ using G_{I\!\!P}^{{\rm cut}}\left(Y-y,\,\boldsymbol{p}_{T},\,0\right)\,\,=\,\,dx\,G\left(x,p_{T}^{2}\right)\Big{/}dp_{T}^{2} , $x_{g}G\left(x_{g},\,M_{J/\psi}\right)$ is the gluon structure function, $x_{g}=M_{J/\psi}e^{y}/\sqrt{s}$, and $Y-y\equiv\ln\left(1/x_{g}\right)$.

Making a Fourier transform, we may rewrite Eq. (7) in the coordinate space as

[TABLE]

where the amplitudes $N\left(y;\,\boldsymbol{r}_{i}\right)$ are related to the solutions of the BK equation as $N\left(y;\,\boldsymbol{r}_{i}\right)\,\,=\,\,\int d^{2}b^{\prime}\,N\left(y;\,\boldsymbol{r}_{i},\,\boldsymbol{b}^{\prime}\right)$, and the variable $\boldsymbol{b}$ is a Fourier conjugate of momentum $\boldsymbol{p}_{1,T}+\text{~}\boldsymbol{p}_{2,T}-\boldsymbol{q}_{T}.$ In what follows, we will also need an expression for the $q_{T}$-integrated cross-section, which takes a simpler form

[TABLE]

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 (II) should be replaced with

[TABLE]

where we replaced the product of pomeron propagators with the double transverse momentum densities $\rho^{(2)}$ defined as

[TABLE]

where $\left(x_{1},p_{1,T}\right)$ and $\left(x_{2},p_{2,T}\right)$ are the light-cone and transverse momenta of the partons in the cascade, 444The definition (11) is closely related to digluon PDFs introduced in Diehl:2011yj , though we have to mention that the probabilistic interpretation strictly speaking can be discussed only for $Q_{T}=0$. As we will see below, in final expressions for the $J/\psi$ production we will need only this case. $|P\rangle$ 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 (II). Similarly, in coordinate space the Eq. (8) can be extended as

[TABLE]

where $\rho^{(2)}\left(x,\,\boldsymbol{r};\,x,\,\boldsymbol{r}^{\prime};\,\boldsymbol{Q}_{T}\right)$ is the double parton density in the coordinate representation, $\boldsymbol{r}_{s}\,=\,\frac{1}{2}(\boldsymbol{r}+\boldsymbol{r}^{\prime})$ and $\boldsymbol{r}_{d}\,=\,\frac{1}{2}(\boldsymbol{r}-\boldsymbol{r}^{\prime})$.

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 $\langle r\rangle$ in the process. The saturation has a mild dependence on energy DKLN and in the kinematics of interest for our studies ($\sqrt{s}\in\left(1.9,\,7\right)$ TeV) its values are $Q_{s}^{2}=0.7-0.9\,{\rm GeV}^{2}$. The typical size of the dipole essential in our process might be estimated as 555See Appendix A and PSIWF for paramertizations of the wave functions

[TABLE]

so the product $\tau\,\,\equiv\,\,\left\langle r^{2}\right\rangle\,Q_{s}^{2}\approx 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 MUT ; GOLEPSI . 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/\psi$ production in DIS

The cross-section of the $J/\psi$ production in the small-$x$ kinematics is closely related to the diffractive $J/\Psi$ production in DIS, and this relationship is useful to fix the unknown non-perturbative factor $\sim\int d^{2}Q_{T}\,S_{h}^{2}(Q_{T})$ which appears in (7,8). In the BFKL picture, there are elastic and inelastic contributions to diffractive $J/\psi$ production, shown in the left and right panels of the Figure 3. Taking into account approximate equality of the elastic and inelastic contributions H1PSI ; H1ELPSI , we may focus on evaluation of the diagram $(a)$ of the Fig. (3). For fixing the prefactor $\sim\int d^{2}Q_{T}\,S_{h}^{2}(Q_{T})$, it is sufficient to consider only the $q_{T}$-integrated cross-section $d\sigma_{{\rm diff}}/dy$, which is given by (see Fig. (3)):

[TABLE]

In the vicinity of the saturation scale, which gives the dominant contribution, the CGC/saturation approach predicts that the amplitude $N$ has a form

[TABLE]

where $\bar{\gamma}\approx 0.63$ Kovchegov:2012mbw , $x\approx M_{J/\psi}e^{-y}/\sqrt{s}$, and $Q_{s}^{2}(x)$ is the saturation scale. The Eq. (III) in this case takes the form

[TABLE]

which allows to factorize all the energy- and rapidity dependence,

[TABLE]

and in this way facilitate scaling from HERA to Tevatron and LHC energies. Similarly, for the gluon induced elastic $J/\psi$ production in the vicinity of the saturation scale (see Ref.GOLEPSI ) we may simplify (9) to

[TABLE]

which also has a factorizable dependence on energy and rapidity,

[TABLE]

Taking into account that the wave functions of the proton and gluon are proportional to each other in the leading order of pQCD DKLMT ,

[TABLE]

we may expect the proportionality of the diffractive and inclusive cross-sections

[TABLE]

where

[TABLE]

and for estimates of parameter $\mathcal{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 H1PSI ; ZEUSPSI ; H1ELPSI , 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\equiv\sqrt{s}\approx 30\,{\rm GeV}$ ($x_{g}\approx 0.01$) the diffractive cross-section

[TABLE]

This allows to to use Eq. (21) for estimate of the hadroproduction cross-section at $\sqrt{s}\approx 30$ GeV,

[TABLE]

In the next section we will extrapolate it with the help of the small-$x$ evolution up to LHC energies.

IV Phenomenological estimates

IV.1 Total $J/\Psi$ production cross section

As we demonstrated in Section II, the typical size $\langle r\rangle$ of $Q\bar{Q}$ dipole is small when the saturation scale $Q_{s}(x)\lesssim 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

[TABLE]

For numerical estimates we assume that the saturation scale $Q_{s}^{2}$ scales as MUT

[TABLE]

with $\lambda\approx 0.2-0.3$ Watt:2007nr ; Kowalski:2003hm ; RESH and $x\approx M_{J/\psi}e^{-y}/\sqrt{s}$. In case of Tevatron and LHC kinematics, this leads to the cross-section estimates given in the Table  1. 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 $\rho^{(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. MOSA , 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 MOSA 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 cross-section $\sigma_{{\rm 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 $\sigma_{{\rm eff}}$ 666This phenomenological parameter manifests significant dependence on channel used to extract it and varies between 5 and 25 ${\rm mb}$ (see Ref.DPI for a review).. In our approach, we fix the the only unknown constant $\int d^{2}Q_{T}\,S^{2}\left(Q_{T}\right)$ 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_{s}^{2}(x)$, namely the choice of coefficient $\lambda$ in (27). Theoretical estimates of this parameter significantly depend on the choice of the scale in the running coupling constant $\bar{\alpha}_{s}$, and lead to estimates $\lambda\approx 0.3-0.4$ CLMP1 , whereas phenomenological estimates of this parameter restrict it to the range $\lambda\in(0.2,\,0.3)$ Kovchegov:2012mbw . We can see that it results in up to fifty per cent uncertainty in the theoretical estimates shown in the Table 1.

IV.2 Transverse momentum distribution

Qualitatively the main features of $q_{T}$-distribution may be understood from (8). For small momenta $\frac{1}{2}q_{T}\lesssim Q_{s}\left(x\right)$, we can safely use the behaviour of the Pomeron Green function in the vicinity of the saturation momentum MUT ****

[TABLE]

In the kinematic region $\frac{1}{2}q_{T}\gtrsim 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\lesssim Q_{s}^{-1}\left(x\right)$ Kovchegov:2012mbw .

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 $\sim(1-2)\,m_{J/\psi}$.

In Fig. (4) we compare model predictions for the shape of the $q_{T}$-dependence with experimental data from PSIATLAS ; PSIALICE . In order to avoid the above-mentioned global uncertainty in normalization, we consider the normalized ratio $\left(d^{2}\sigma/dy\,dq_{T}\right)/\left(d\sigma/dy\right)$. The cusp on the curve near the threshold scale $q_{T}=q_{0}$ can be smoothed out by more relaxed conditions.

IV.3 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):

[TABLE]

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

[TABLE]

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|\gtrsim 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

[TABLE]

per each gluon/pomeron in order to have correct endpoint behaviour of the gluon densities in the $x\to 1$ limit BBSH . 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.

IV.4 Multiplicity dependence

The $J/\psi$ production accompanied by two parton showers occurs in the events with larger multiplicity than the average multiplicity $\bar{n}$ in the inclusive hadron production  PSIMULT ; Alice:2012Mult . 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 Fig. (6)-a and Fig. (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.KOLEB ; KLN 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 LERE ,

[TABLE]

where $\bar{n}$ is the average multiplicity at $y^{*}=0$ ($\bar{n}\approx$6.5 at $W=13\,{\rm TeV}$ ALICEINCL ). As we discussed earlier, the product of dipole size on saturation scale $\left\langle r^{2}\right\rangle Q^{2}\left(y,\,\bar{n}\right)$ is close to one, for this reason the multiplicity events with $n/\bar{n}\gg 1$ probe a deep saturation regime, where approximation (18) is not valid, and we should use (9) with phenomenological parameterization of the dipole amplitude $N$.

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 PSIMULT is reasonable. At sufficiently small multiplicities, $dN_{J/\psi}/dy$ is increasing as a function of $dN_{{\rm ch}}/d\eta$, 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) LETU , for this reason the difference $\left[N\left(y;\,\frac{\vec{r}+\vec{r}^{\prime}}{2},\,0\right)\,-\,N\left(y;\,\frac{\vec{r}-\vec{r}^{\prime}}{2},\,0\right)\right]^{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) contribute with relative weights $\sim\left(Q_{s}^{2}\left(y,n\right)\right)^{2\bar{\gamma}}$ and  $\sim\left(Q_{s}^{2}\left(y,n\right)\right)^{\bar{\gamma}}$ respectively, we may obtain for the self-normalized multiplicity dependence of contribution in Fig. (6) a simple expression

[TABLE]

where the coefficient $\kappa=\left(Q_{s}^{2}\left(Y-y\right)/Q_{s}^{2}\left(y\right)\right)^{\bar{\gamma}}\approx\,\,e^{\bar{\gamma}\,\lambda\left(Y-2y\right)}$ 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 $\kappa\approx 1$). While the LO CGC/Saturation predicts a value $\bar{\gamma}\approx 0.63$ , the phenomenological fits Watt:2007nr ; RESH favor slightly higher values of $\bar{\gamma}$, for this reason we varied this parameter in the range $\bar{\gamma}\in\left(0.67,\,0.76\right).$ 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/\psi$ 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/\psi$ hadroproduction, accompanied by production of two parton showers. We demonstrated that this mechanism gives a large contribution to the total cross section of $J/\psi$ production, as well as to the transverse momenta and multiplicity distribution of produced $J/\psi$. The experimental data PSIMULT on multiplicity distributions of produced $J/\psi$ seem to favor this hypothesis (compared to conventional mechanisms estimated in the CGC/Saturation approach). This mechanism has no suppression of the order of $\bar{\alpha}_{S}$ at high energies compared to one parton shower contribution. We show that it alone can describe the shape of the $p_{T}$- and 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.MOSA 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/\psi$.

VI Acknowledgements

We thank our colleagues at Tel Aviv university and UTFSM for encouraging discussions. This research was supported by the BSF grant 2012124, by Proyecto Basal FB 0821(Chile), Fondecyt (Chile) grant 1180118 and by CONICYT PIA grants ACT1406 and ACT1413.

Appendix A Evaluation of the factor $I\left(\boldsymbol{p}_{1,T},\boldsymbol{p}_{2,T}\right)$

In this section we would like to comment briefly the evaluation of the function $I\left(\boldsymbol{k}_{T},\boldsymbol{q}_{T}\right)$ in (3) in the light cone perturbative QCD KOLEB .

We choose the polarization vectors of gluons $\epsilon_{\mu}^{\lambda}(p)$ in the light-cone gauge as

[TABLE]

so according to the general rules of light cone perturbation theory (LCPT)  MUDI ; KOLEB the wave function of the $\bar{Q}Q$ pair after interaction with a pair of gluons with momenta ($y_{1,2},\,\boldsymbol{p}_{1,2;\,T}$) is given by

[TABLE]

where $\lambda^{a}$ is the Gell-Mann matrix, $\boldsymbol{\epsilon}^{\alpha}$ is the polarization vector of the gluon with the helicity $\alpha$, and $\Psi_{g}\left(\boldsymbol{l}_{T},\,z\right)$ is the wave function of $\bar{Q}Q$ pair formed after gluon splitting into quark-antiquark pair with transverse momentum $\boldsymbol{l}_{T}$ and light-cone fraction $z$.

In the coordinate space the convolution of (35) with the wave function of $J/\psi$ yields

[TABLE]

which agrees with (3) if we change notations of momenta to $\boldsymbol{p}_{1,2\,T}=\frac{\boldsymbol{q}_{T}}{2}\pm\boldsymbol{k}_{T}$. In view of the heavy mass of the charm quark, in what follows we will use leading order perturbative results for the gluon wave functions $\Psi_{g}$ , as well as phenomenological parameterization of the $J/\psi$ wave function available from the literature PSIWF , so the convolution of the two objects explicitly takes the form

[TABLE]

Appendix B Evolution equation for the double gluon density

In Ref.LELU it is proven that the evolution equations for all partonic densities $\rho^{(n)}\left(\{\vec{r}_{i},\vec{b}_{i}\}\right)$ are the linear BFKL evolution equations. The non-linear corrections are essential for the scattering amplitude (see, for example, the Balitsky-Kovchegov equation for dipole scattering amplitudeBK with nuclei) but they do not give a contribution to the evolution equation for multi-gluon densities. Referring our reader to Ref.LELU for the proof and details, we present here the resulting evolution equation for $r_{1}^{2}\,r_{2}^{2}\,\int d^{2}b\,\rho_{A}^{(2)}\left(x,\vec{r};x,\vec{r}^{\prime};\vec{b}_{T}\right)\,=\,\tilde{\rho}^{(2)}\left(Y-Y_{0},\vec{r}_{1},\vec{r}_{2}\right)$:

[TABLE]

where $\left(\vec{r}_{1}+\vec{r}_{2}\right)^{2}\int d^{2}b\,\rho^{(1)}\left(Y-Y_{0},\vec{r}_{1}+\vec{r}_{2},\,\vec{b}\right)\,\,=\,\,\tilde{\rho}^{(1)}\left(Y-Y_{0},\vec{r}_{1}+\vec{r}_{2}\right)$. The two terms in Eq. (41) have clear physical meaning: the first one is the evolution of two parton showers, while the second describes the production of two gluon in one parton shower (see Fig. (9)).

The solution takes the following form (see Ref.GOLEPSI and references therein):

[TABLE]

where BFKL ; KOLEB $\xi_{1}\,=\,\ln\left(r_{1}^{2}\right)$, $\xi_{2}\,=\,\ln\left(r_{2}^{2}\right)$ and

[TABLE]

The functions$\tilde{\rho}_{{\rm in}}^{(2)}\left(\gamma_{1},\gamma_{2}\right)$ and $\tilde{\rho}_{{\rm in}}^{(1)}\left(\gamma_{1}+\gamma_{2}\right)$ are determined by the initial condition for the two and one parton shower productions (see Fig. (9). The explicit form of the function $h\left(\gamma_{1},\,\gamma_{2}\right)$ can be found in Ref.GOLEPSI . Different contributions in the solution (42) might be visualized as shown in the Figure 9.