Revisiting the production of $J/\psi$ pairs at the LHC

We consider the prompt double $J/\psi$ production in $pp$ collisions at the LHC in the framework of $k_T$-factorization QCD approach. Using the fragmentation mechanism, we evaluate the color octet contributions to the production cross sections taking into account the combinatorial effects of multiple gluon radiation in the initial state driven by the Ciafaloni-Catani-Fiorani-Marchesini evolution equation. We demonstrate the importance of these contributions in a certain kinematical region covered by the CMS and ATLAS measurements. On the other hand, the experimental data taken by the LHCb Collaboration at forward rapidities and moderate transverse momenta can be described well by ${\cal O}(\alpha_s^4)$ color singlet terms and contributions from the double parton scattering mechanism. The extracted value of the effective cross section $\sigma_{\rm eff} = 17.5 \pm 4.1$~mb is compatible with many other estimations based on different final states.


Introduction
: Contribution to the J/ψ pair production from the fragmentation of gluon cascade. The dashed line encloses the hard subprocess g * g * → g * , the rest of the diagram describes the initial state radiation cascade.
Our other goal is connected with the investigation of additional production mechanism, double parton scattering (DPS), which is widely discussed in the literature at present (see, for example, [30][31][32][33][34] and references therein). Apart from the single parton scattering (SPS), where J/ψ meson pair is produced in a single gluon-gluon collision, DPS events originate from two independent parton interactions. Studying the DPS mechanism is of great importance since it can help in understanding various backgrounds in searches for new physics at the collider experiments. Despite the relative low total production rate, the DPS mechanism is expected to be important for double J/ψ production at forward rapidities [35][36][37]. Therefore, the latter can be used to determine the DPS key parameter, the effective cross section σ eff , which is related to the transverse overlap function between partons in the proton and supposed to be universal for all processes with different kinematics and energy scales. Most of the measured values of σ eff lie between 12 and 20 mb (see, for example, [38,39]). However, somewhat lower value σ eff = 8.8 − 12.5 mb was extracted from the latest LHCb data on J/ψ pair production within the NRQCD [4]. Moreover, the values σ eff = 8.2 ± 2.2 mb [40], σ eff = 6.3 ± 1.9 mb [41], σ eff = 4.8 ± 2.5 mb [42] and even σ eff = 2.2 ± 1.1 mb [43], σ eff = 2.2 − 6.6 mb [44] were obtained from recent Tevatron and LHC experiments. Below we will try to extract the effective cross section σ eff from combined analysis of the LHCb data [3,4] on double J/ψ production taken at √ s = 7 and 13 TeV.
To calculate the physical cross sections we use the k T -factorization approach [45,46]. We see certain technical advantages in the fact that, even with the LO amplitudes for hard subprocesses, one can include a large piece of higher-order pQCD corrections (NLO + NNLO + ...) taking them into account in the form of CCFM-evolved Transverse Momentum Dependent (TMD) gluon densities in a proton 3 . In this way we preserve consistency with our previous studies [18][19][20][21] and automatically incorporate the wanted effects of initial state gluon radiation. To reconstruct the CCFM evolution ladder, that is the key point of our consideration, we employ the TMD parton shower routine implemented into the Monte-Carlo event generator cascade [48]. The k T -factorization approach can be considered as a convenient alternative to explicit high-order calculations in the collinear DGLAP-based scheme. The situation in J/ψ pair production is specific since calculating even the LO hard scattering amplitudes is already complicated enough, so that extending to higher orders seems to be a rather cumbersome task. Thus, the k T -factorization remains the only way open to potentially important higher-order effects. To evaluate the DPS contributions to the double J/ψ production we will use the results of our previous analysis [19].
The outline of the paper is the following. In Section 2 we briefly describe the basic steps of our calculations. In Section 3 we present the numerical results and discussion. Our conclusions are summarised in Section 4.

The model
The neccessary starting point of our consideration is related with CS contribution to the double J/ψ production, that refers to O(α 4 s ) gluon-gluon fusion subprocess where the four-momenta of all particles are indicated in the parentheses. Some typical Feynman diagrams are depicted in Fig. 2. It is important that both initial gluons are off mass shell. That means that they have non-zero transverse four-momenta k 2 1T = −k 2 1T = 0 and k 2 2T = −k 2 2T = 0 and an admixture of longitudinal component in the polarization fourvectors (see [45,46] for more information). The corresponding off-shell (k T -dependent) production amplitude contains widely used projection operators for spin and color [49] which guarantee the proper quantum numbers of final state charmonia. Below we apply the gauge invariant expression obtained earlier [50]. The derivation steps are explained in detail there. The respective cross section can be written as where ψ 1 is the azimuthal angle of outgoing J/ψ meson, φ 1 and φ 2 are the azimuthal angles of initial gluons having the longitudinal momentum fractions x 1 and x 2 , y 1 and y 2 are the center of mass rapidities of produced particles and f g (x, k 2 T , µ 2 ) is the TMD gluon density in a proton taken at the scale µ 2 .
In addition to the CS terms above, we have considered some of CO contributions using the fragmentation approach. At high transverse momenta, p T m ψ , large logaritmic corrections proportional to α n s ln n p T /m ψ occur and, therefore, description in terms of fragmentation functions (FFs), evolving with the energy scale µ 2 , appears to be appropriate. In general, the FF D H a (z, µ 2 ) describing the transition of parton a into the charmonium state H can be expressed as follows (see, for example, [51] and references therein): where n labels the intermediate (CS or CO) state of charmed quark pair produced in the hard parton interaction and O H [n] are the corresponding LDMEs. In the leading logarithmic approximation, g * → cc[ 3 S 1 ] transition is the only one giving a sizeble contribution to S-wave charmonia production at p T m ψ [51], so that the cross section of inclusive single J/ψ production in pp collisions could be approximately calculated as where p = zp (g * ) and p (g * ) are the outgoing J/ψ meson and intermediate gluon momenta. One can easily obtain where p (g * ) = k 1 + k 2 and λ(m 2 ψ , k 2 1 , k 2 2 ) is the known kinematical function [52]. Evaluation of the off-shell production amplitude |Ā(g * g * → g * )| 2 = (3/2)πα s (µ 2 )|p (g * ) T | 2 is an extremely straightforward and, in our opinion, needs no explanation. We only note that, according to the k T -factorization prescription [45,46], the summation over the polarizations of initial off-shell gluons is carried out with µ * ν = k µ T k ν T /k 2 T . In the collinear limit k T → 0 this expression converges to the ordinary one after averaging on the azimuthal angle.
The key point of our consideration is that the gluon, produced in the hard scattering and fragmented into the J/ψ meson according to main formula (4), is accompanied by a number of gluons radiated during the non-collinear QCD evolution, which also give rise to final J/ψ mesons. Thus, taking into account all their possible combinations into the meson pairs, one can calculate corresponding gluon fragmentation contribution to the double J/ψ production up to all orders in the pQCD expansion. At high energies, the QCD evolution of gluon cascade can be described by the CCFM equation [29], which smoothly interpolates between the small-x BFKL gluon dynamics and high-x DGLAP one, and, therefore, provides us with the suitable tool for our phenomenological study. To reconstruct the CCFM evolution ladder, we generate a Les Houches Event file [53] in the numerical calculations according to (4) and (5) and then process the file with a TMD shower tool implemented into the Monte-Carlo event generator cascade [48]. This approach gives us the possibility to take into account the contributions from initial state gluon emissions in a consistent way (see also [54]).
Of course, the same scenario can be applied to fragmentation of charmed quark pairs into J/ψ mesons. So, one can first simulate the perturbative production of cc pair in the off-shell gluon-gluon fusion and then reconstruct the CCFM evolution ladder using the cascade tool. After that, one can easily produce J/ψ pairs by taking into account all possible combinations of mesons originating from the charmed quarks and/or cascade gluon fragmentation. Unlike the conventional (collinear) QCD factorization, where only fragmentation of both charmed quarks into J/ψ mesons gives contribution, the model above can lead to increase in the double J/ψ production cross section due to additional combinatorial contributions from gluons and quarks.
The charm and gluon FFs at the any scale µ 2 , D J/ψ c (z, µ 2 ) and D J/ψ g (z, µ 2 ), can be obtained by solving the LO DGLAP evolution equations: where P ab are the standard LO DGLAP splitting functions. The initial conditions for these FFs are calculated with [51] where starting scale µ 2 0 = m 2 ψ . As it was noted above, we keep only the leading contributions to corresponding FFs (see, for example, [51] and references therein). According to the non-relativistic QCD approximation, we set the charmed quark mass to m c = m ψ /2 and then solve the DGLAP equations (6) numerically. The obtained charm and gluon FFs, D J/ψ c (z, µ 2 ) and D J/ψ g (z, µ 2 ), are shown in Fig. 3 as functions of z for several values of scale µ 2 . Using these FFs, we reproduce well the results of calculations performed with the Monte Carlo event generator pegasus [55] (see Fig. 4).
Finally, we turn to the DPS contribution to the double J/ψ production. We apply a commonly used factorization formula (for details see reviews [30][31][32][33][34] and references therein): where factor 1/2 accounts for two identical particles in the final state. The effective cross section σ eff can be considered as a normalization constant which incorporates all "DPS unknowns" in to a single phenomenological parameter. Derivation of the factorization formula (9) relies on the two approximations: first, the double parton distribution function can be decomposed into longitudinal and transverse components and, second, the longitudinal component reduces to the diagonal product of two independent single parton densities. The latter is generally acceptable for such collider experiments where small-x values are probed. The typical values of the variable x in the considered process are of order x ∼ (2m 2 ψ + p 2 T ) 1/2 / √ s ∼ 10 −3 , that approximately corresponds to the kinematical region of CMS [1], ATLAS [2] and even LHCb [3,4] measurements (due to relatively small invariant mass of produced J/ψ pair, see discussion below). Therefore, one can safely omit the kinematical constraint [56,57] often applied at the edge of phase space 4 . Detailed description of evaluation of inclusive cross section σ(pp → J/ψ + X) in the k T -factorization approach supplemented with the NRQCD formalism can be found [19].
In the numerical calculations below we will use TMD gluon density in a proton obtained [59] from the numerical solution of CCFM evolution equation (namely, A0 set), where the input parameters have been fitted to the proton structure function F 2 (x, Q 2 ). At present, the A0 gluon distribution function is widely used in the phenomenological applications 5 (see, for example, [19][20][21]). The renormalization and factorization scales, µ R and µ F , were set to µ 2 R = µ 2 F =ŝ + Q 2 T , whereŝ = (k 1 + k 2 ) 2 and Q 2 T is the transverse momentum of initial off-shell gluon pair. This choice is dictated mainly by the CCFM evolution algorithm (see [59] for more information). As it is often done, the fragmentation scale µ fr is choosen to be equal to µ fr = m T , the transverse mass of fragmenting parton. We use the one-loop formula for the QCD coupling α s with n f = 4 active quark flavors at Λ (4) QCD = 250 MeV. Following [61], we set the J/ψ meson mass m ψ = 3.097 GeV. We take corresponding CS LDME from the known J/ψ → µ + µ − decay width:

Numerical results and discussion
We are now in a position to present the results of our simulations. First we discuss the role of cascade gluon fragmentation in different kinematical regimes, which correspond to the CMS, ATLAS and LHCb experiments.
In Fig. 5 we show the differential cross sections of double J/ψ production calculated as a functions of J/ψ pair invariant mass m(J/ψ, J/ψ) and difference in rapidity between the J/ψ mesons |∆y(J/ψ, J/ψ)| at √ s = 13 TeV. We have required p T (J/ψ) > 10 GeV for both produced mesons, that ensures the validity of the fragmentation approach used. Moreover, this restriction close to the CMS or ATLAS conditions. One can see that an accurate account of combinatorial contributions, originated from the cascade gluon fragmentation into the J/ψ mesons (labeled as "fragm. comb."), significantly (up to an order of magnitude) increase the cross section compared to the single gluon fragmentation, governed by the LO gluon-gluon fusion subprocess 6 (labeled as "fragm. coll."). For the latter, we reproduce the results [27]. To highlight the importance of the combinatorial gluon fragmentation, we show the results obtained using the simplified selection of J/ψ pair in each event, where one of the J/ψ mesons is originated from the gluon produced in the hard scattering subprocess and another one is produced from the leading cascade gluon (labeled as "fragm. lead."). This selection criterion almost corresponds to the collinear limit, as it is clearly demonstrated in Fig. 5. Next, we find that the cascade gluon fragmentation plays a dominant role at large invariant masses m(J/ψ, J/ψ) ≥ 25 GeV and |∆y(J/ψ, J/ψ)| ≥ 1, where it greatly overestimates the CS contributions. Taking into account these combinatorial contributions results in the drastical rise of the double J/ψ production cross sections at large m(J/ψ, J/ψ), where the strong discrepancy between the NRQCD estimations (including both the CS and CO terms) and experimental data, taken by the CMS and ATLAS Collaborations, is observed. Contrary, the combinatorial fragmentation effects should be significantly less pronounced at forward rapidities, which are covered by the LHCb measurements. To demonstrate it, we have repeated the calculations under the requirements 4.5 < p T (J/ψ) < 10 GeV and 2 < y(J/ψ) < 4.5. The upper limit of p T (J/ψ) is set to be the same as in LHCb analyses [3,4] while lower limit corresponds to the region, where the fragmentation approach is valid. Our results for distributions in m(J/ψ, J/ψ) and |∆y(J/ψ, J/ψ)| are shown in Fig. 6. One can see that the cascade gluon fragmentation gives only small contribution to the forward J/ψ pair production and, in principle, can be safely neglected. It can be easily understood since at large rapidities (or, equivalently, at large momentum fraction x of one of the interacting gluons) the gluon emissions in the initial state are insufficient.
Concerning the contributions from charm fragmentation, their role (compared to the LO predictions of conventional pQCD) is also a bit enhanced due to the multiple gluon emissions in the initial state. We find that these processes amount several percent of the J/ψ pair production cross section (see Figs. 5 and 6) and, of course, can be considered as additional non-leading terms 7 .
Thus, we have shown that taking into account the combinatorial contributions from the cascade gluon fragmentation could fill the gap between the NRQCD predictions and experimental data. However, to perform the quantitative comparison with the available CMS [1] and ATLAS [2] measurements one has to include into the analysis a number of other possible fragmentation channels playing role at low and moderate transverse momenta. Moreover, additional feeddown contributions to the double J/ψ production from the χ c and ψ decays should be taken into account. An accurate theoretical description requires a rather long-time numerical calculations. So, here we only claim the possible importance of the combinatorial fragmentation terms above and left their further cumbersome analysis for a forthcoming dedicated study. Now we turn to available LHCb data collected at √ s = 7 and 13 TeV [3,4]. These data refer to p T (J/ψ) < 10 GeV, m(J/ψ, J/ψ) < 15 GeV and forward rapidity region, 2 < y(J/ψ) < 4.5. Since the combinatorial contributions from gluon and/or charmed quark fragmentation are almost negligible there, only the CS terms and DPS production mechanism play the role. The latter give us the possibility to easily extract the key parameter of DPS mechanism, the effective cross section σ eff , from the LHCb measurements. The feeddown contributions from radiative χ c and ψ decays to the SPS cross section, which is governed by the subprocess (1), are also unimportant at small transverse momenta and invariant mass m(J/ψ, J/ψ), see discussions [37,63]. Thus, we neglect below all these terms for simplicity. To evaluate the DPS contribution to the J/ψ pair production we use the results of our previous studies and strictly follow the approach [19] for the inclusive cross section σ(pp → J/ψ + X), entering to the DPS factorization formula (9). So, the determination of σ eff can be performed in a self-consistent way.
The following kinematical variables have been investigated in the LHCb analyses [3,4]: transverse momentum p T (J/ψ, J/ψ), rapidity y(J/ψ, J/ψ) and invariant mass of the J/ψ pair, transverse momentum and rapidity of J/ψ mesons, differences in the azimuthal angle |∆φ(J/ψ, J/ψ)| and rapidity |∆y(J/ψ, J/ψ)| between the produced mesons and transverse momentum asymmetry A T , defined as The measurements have been performed for p T (J/ψ, J/ψ) > 1 GeV, p T (J/ψ, J/ψ) > 3 GeV and in the whole p T (J/ψ, J/ψ) range. We consider σ eff as an independent parameter and perform a simultaneous fit to the LHCb data. The fitting procedure was separately done for each of the measured kinematical distributions employing the fitting algorithm as implemented in the commonly used gnuplot package [64]. Not all of the existing data sets are equally informative for the σ eff extraction. Using the data where the DPS contribution is smaller than the uncertainty of the "main" contribution would only increase the total error. So, our fit is based on the following distributions (all measured at √ s = 7 TeV and 13 TeV): single J/ψ transverse momentum p T (J/ψ); single J/ψ rapidity y(J/ψ); invariant mass m(J/ψ, J/ψ); transverse momentum of J/ψ pair; rapidity of J/ψ pair; transverse momentum asymmetry A T ; rapidity separation between the two J/ψ mesons |∆y(J/ψ, J/ψ)|. For all observables except |∆y(J/ψ, J/ψ)| we used the data without cuts on p T (J/ψ, J/ψ) and with p T (J/ψ, J/ψ) > 1 GeV, while for the rapidity separation |∆y(J/ψ, J/ψ)| we used the sets without cuts on p T (J/ψ, J/ψ), with p T (J/ψ, J/ψ) > 1 GeV, and with p T (J/ψ, J/ψ) > 3 GeV.
A comparison of our predictions with the LHCb experimental data is displayed in Figs. 7 -10. The theoretical uncertainty bands include both scale uncertainties and uncertainties coming from the σ eff fitting procedure. First of them have been estimated in a usual way, by varying the µ R scale around its default value by a factor of 2. This was accompanied with using the A0+ and A0− gluon densities instead of default A0 distribution, in accordance with [59]. As one can see, we achieved a reasonably good agreement between the results of our calculations and LHCb measurements, both for √ s = 7 and 13 TeV. There is only exception in the threshold region, m(J/ψ, J/ψ) ≤ 9 GeV, where our predictions systematically overshoot the data. However, at such low m(J/ψ, J/ψ) an accurate treatment of multiple soft gluon emissions, relativistic corrections and other nonperturbative effects becomes necessary to produce the theoretical estimations. All these issues are out from our present consideration. Next, we find that neither the SPS terms, nor the DPS contributions alone are able to describe the LHCb data, but only their sum. In particular, the DPS contributions are essential to reproduce the measured |∆y(J/ψ, J/ψ)| distributions at |∆y(J/ψ, J/ψ)| ≥ 1 or 1.5, that confirms the previous expectations [35][36][37]. They are important to describe also the normalization of J/ψ rapidity distributions and shape of transverse momentum asymmetry A T at A T ≤ 0.4, see Figs. 8 -10.
The presented results, being considered altogether with the ones for inclusive single production of charmonia states [19], can give a significant impact on the understanding of charmonia production within the NRQCD framework and, in particular, on the further understanding of DPS mechanism. The most interesting outcome of our study is that the extremely low value of DPS effective cross section, σ eff ∼ 2 − 5 mb, obtained in earlier analyses of double J/ψ production at the LHC, is not confirmed.

Conclusion
We have considered the prompt production of J/ψ meson pairs in pp collisions at the LHC using the k T -factorization approach of QCD. We employ the fragmentation mechanism to evaluate the color octet contributions to the production cross sections and take into account the combinatorial effects of multiple gluon radiation in the initial state using the CCFM evolution equation. The latter could be essential in the kinematical region covered by the CMS and ATLAS measurements. On the other hand, we have demonstrated that the experimental data taken by the LHCb Collaboration at forward rapidities can be described well by the color singlet terms and contributions from the double parton scattering mechanism. We determine the DPS effective cross section σ eff = 17.5 ± 4.1 mb from the combined analysis of the LHCb data collected at √ s = 7 and 13 TeV. The extracted value is compatible with many other estimations based on essentially different final states. The extremely low σ eff ∼ 2 − 5 mb, obtained earlier from the double J/ψ production data, is not confirmed.