Double prompt $J/\psi$ hadroproduction in the parton Reggeization approach with BFKL resummation

We study double prompt $J/\psi$ hadroproduction within the nonrelativistic-QCD factorization formalism adopting the parton Reggeization approach to treat initial-state radiation in a gauge invariant and infrared-safe way. We present first predictions for the cross section distributions in the transverse momenta of the subleading $J/\psi$ meson and the $J/\psi$ pair. Already at leading order in $\alpha_s$, these predictions as well as those for the total cross section and its distributions in the invariant mass $m_{\psi\psi}$ and the rapidity separation $|Y|$ of the $J/\psi$ pair nicely agree with recent ATLAS and CMS measurements, except for the large-$m_{\psi\psi}$ and large-$|Y|$ regions, where the predictions substantially undershoot the data. In the latter regions, BFKL resummation is shown to enhance the cross sections by up to a factor of two and so to improve the description of the data.

We study double prompt J/ψ hadroproduction within the nonrelativistic-QCD factorization formalism adopting the parton Reggeization approach to treat initial-state radiation in a gauge invariant and infrared-safe way. We present first predictions for the cross section distributions in the transverse momenta of the subleading J/ψ meson and the J/ψ pair. Already at leading order in αs, these predictions as well as those for the total cross section and its distributions in the invariant mass m ψψ and the rapidity separation |Y | of the J/ψ pair nicely agree with recent ATLAS and CMS measurements, except for the large-m ψψ and large-|Y | regions, where the predictions substantially undershoot the data. In the latter regions, BFKL resummation is shown to enhance the cross sections by up to a factor of two and so to improve the description of the data. Despite concerted experimental and theoretical endeavors ever since the discovery of J/ψ meson more than four decades ago, its production mechanism has remained mysterious; for a recent review, see Ref. [1]. The factorization approach [2] to nonrelativistic QCD (NRQCD) [3] endowed with velocity scaling rules [4] for the longdistance matrix elements (LDMEs), which is by far the most acceptable candidate theory for heavy-quarkonium production and decay and has been elaborated at nextto-leading order (NLO), has been challenged by the longstanding J/ψ polarization puzzle [5] and by the inadequate description of η c hadroproduction data with J/ψ LDMEs converted via heavy-quark spin symmetry [6]. The high flux of incoming partons at LHC allows us to study J/ψ production more thoroughly, also in association with other charmonia, bottomonia, W or Z bosons, so as to pin down the J/ψ production mechanism. Among these production processes, double J/ψ hadroproduction is of special interest because J/ψ formation takes place there twice, making this particularly sensitive to the nonperturbative aspects of NRQCD, allowing for an independent test of the predicted universality of the LDMEs, and providing additional strong constraints on their values [7]. Moreover, this is believed to be an exquisite laboratory to study double parton scattering (DPS) and to extract its key parameter, the effective cross section σ eff [8].
In recent years, double prompt J/ψ hadroproduction has been measured extensively by the LHCb [9], CMS [10], and ATLAS [11] Collaborations at the CERN LHC, and by the D0 Collaboration [12] at the FNAL Tevatron. On the theoretical side, the complete NRQCD results at leading order (LO), including the color octet (CO) contributions due to the cc Fock states 1 S [8] 0 , 3 S and the feed-down contributions due to the χ cJ and ψ ′ mesons have been obtained recently [13].
For some channels, the relativistic corrections [14] and NLO QCD corrections [15] are also available. The available NRQCD predictions can explain the LHCb and D0 data, and to some extent also the CMS data, reasonably well. However, these single-parton scattering (SPS) predictions only amount to a few percent of the CMS data in the regions of large invariant mass m ψψ or rapidity separation |Y | = |y 1 − y 2 | of the J/ψ pair, although the CO contributions, in particular those involving t-channel gluon exchanges (see Feynman diagrams in Fig. 1 of Ref. [13]), enhance the QCD-corrected [15] color-singlet (CS) contribution of gg → 2cc( 3 S [1] 1 ) there by more than one order of magnitude. The value of σ eff extracted by fitting the DPS contribution on top of this [12] is considerably smaller than typical values from other processes, and the resulting SPS plus DPS results still undershoot the CMS data in the upper m ψψ and |Y | bins [16].
As noticed in Ref. [13], the CMS kinematic conditions [10] render 2 → 3 subprocesses predominant, which enter at NLO in the collinear parton model (CPM). However, a complete NLO NRQCD computation is presently out of reach from the technical point of view. On the conceptual side, the conventional NRQCD factorization formalism needs to be extended to cope with the double P -wave case [17]. Moreover, the perturbative expansion is spoiled for small values of the J/ψ pair transverse momentum p ψψ T . The characteristic scale µ ∼ [(4m c ) 2 + (p ψψ T ) 2 ] 1/2 of the hard-scattering processes of double J/ψ hadroproduction satisfies Λ QCD ≪ µ ≪ √ S, where Λ QCD is the asymptotic scale parameter and √ S is the centerof-mass energy. We are thus accessing a new dynamical regime, namely the high-energy Regge limit, where the NLO QCD corrections can be largely accounted for through the unintegrated parton distribution functions (unPDFs) in the parton Reggeization approach (PRA) [18] based on Lipatov's effective field theory (EFT) formulated with the non-Abelian gauge invariant action [19]. The PRA has already been successfully applied to the interpretation of measurements of single heavyquarkonium hadroproduction [20,21]. In the large-|Y | region, the two J/ψ mesons are well separated obeying multi-Regge kinematics (MRK). For subprocesses containing t-channel gluon exchange type diagrams, there will be large logarithms of form (α s ln |s/t|) n in the higher-order QCD corrections, where s and t are the Mandelstam variables of the partonic 2 → 2 Born process. Such large logarithms can be resummed by the Balitsky-Fadin-Kuraev-Lipatov (BFKL) formalism [22]. Recently, BFKL resummation has been studied for single J/ψ [23] and J/ψ plus jet [24] inclusive hadroproduction. In this Letter, we will take a crucial step towards a full-fledged NLO NRQCD study of double prompt J/ψ hadroproduction, by adopting the PRA and performing BFKL resummation.
Owing to the PRA and NRQCD factorization, the cross section of inclusive double prompt J/ψ hadroproduction can be expressed as: where dσ PRA mn is the short-distance coefficient (SDC) of the partonic subprocess if H = J/ψ. Since partonic subprocesses initiated by Reggeized quarks and antiquarks are greatly suppressed by their unPDFs, we may disregard them here. Furthermore, we may neglect the χ c0 feed-down contribution because BR(χ c0 → J/ψ + X) = 1.40% [25] is so small. Representative Feynman diagrams for the partonic subprocess R + R − → cc(m)cc(n) at LO in α s are depicted in Figs. 1(a) and (b). By the EFT Feynman rules [19], they come in the same topologies as those for gg fusion in the CPM. Because of J P C conservation, for given m and n, not all diagrams contribute. To discuss the BFKL resummation effect conveniently, we divide the partonic subprocesses into three categories, according to the order in α s where t-channel gluon exchange emerges, namely at LO, NLO, and next-to-next-to-leading order (NNLO): (i) LT, with m, n = 1 S We first compute the LO contributions to all the three categories. Due to lack of space, we relegate the details of our calculation to a separate paper. In contrast to other k T factorization approaches, the PRA yields gauge invariant SDCs with off-shell initial-state partons, which provides a strong check for our analytic calculations. In the collinear limits t 1,2 → 0, we recover the CPM formulas [13], which constitutes yet another nontrivial check. In the numerical analysis, we adopt the Kimber-Martin-Ryskin scheme [26] to generate the unPDF from the MSTW 2008 set of LO CPM PDFs [27], which comes with α s (M Z ) = 0.13939. We choose the renormalization and factorization scales to be T + p H2 T )/2, and ξ is varied between 1/2 and 2 about its default value 1 to estimate the unphysicalscale uncertainty. In the case of feed-down from H = χ c1 , χ c2 , ψ ′ , we put p [28]. For the J/ψ, χ cJ , and ψ ′ LDMEs, we use the values specified in Table I. The CS results have been derived from the Buchmüller-Tye potential in Ref. [29]. The CO results have been fitted to LHC data of inclusive single charmonium hadroproduction [30] in the very theoretical framework described above; they supersede pre-LHC results [20]. For H = J/ψ, ψ ′ , there is a strong correlation be- 0 ) , so that only a linear combination of them can be determined. We may thus put O H ( 3 P [8] 0 ) = 0. We have checked that the theoretical uncertainties in our predictions for double prompt J/ψ hadroproduction due to this freedom are negligible. The masses M J/ψ = 3.097 GeV, M χc1 = 3.511 GeV, M χc2 = 3.556 GeV, M ψ ′ = 3.686 GeV, and branching fractions Br(χ c1 → J/ψγ) = 34.3%, Br(χ c2 → J/ψγ) = 19.0%, Br(ψ ′ → J/ψ + X) = 61.4% are adopted from Ref. [25].
Since the LHCb [9] and D0 [12] measurements of prompt double J/ψ hadroproduction agree reasonably well with the LO NRQCD predictions in the CPM [13], we may focus here on the CMS [10] and ATLAS [11] ones. The CMS data were taken at √ S = 7 TeV requiring for each J/ψ meson to be in the rapidity range |y| < 2.2 and to satisfy a y-dependent minimum-p T cut, as described in Eq. (3.3) of Ref. [10]. The CMS total cross   Fig. 2(a)-(c), respectively. There is generally very good agreement, except for the upper two m ψψ bins and the upmost |Y | bin, where the predictions undershoot the data. The advancement of the PRA beyond the CPM is most striking for Fig. 2(a) because p ψψ T = 0 at LO in the latter case, but it is also significant for the m ψψ and |Y | distributions, as may be observed by comparing Figs. 2(b) and (c) with their CPM counterparts in Figs. 3 and 4 of Ref. [13]. In both cases, the predictions are substantially increased in the first three bins, so as to nicely match the data, while the K factors are of order unity in the upmost bins. It is interesting to observe that, in the upmost m ψψ and |Y | bins, the LT contributions almost exhaust the cross section, making up 88% and 96% of it, respectively, while they are comparable to the NLT and NNLT ones in the residual bins.
A LO NRQCD study of double J/ψ hadroproduction has also been performed in a non-gauge-invariant k Tfactorization formalism, including the 3 S  CO channels as well as the feed-down from the χ cJ and ψ ′ mesons [31]. The m ψψ and |Y | distributions thus predicted agree, within their uncertainty bands, with the CPM results of Ref. [13], and so fail to describe the CMS data [10] both in the lower and upper m ψψ and |Y | bins. Also the CMS p ψψ T distribution is poorly described, especially in its peak region.
ATLAS took their data at √ S = 8 TeV imposing the acceptance cuts p T > 8.5 GeV and |y| < 2.1 on each J/ψ meson [11]. They separately studied the central (I) and forward (II) y regions of the subleading J/ψ meson (J/ψ 2 ), with p 2T < p 1T , namely |y 2 | < 1.05 and 1.05 < |y 2 | < 2.1. Their respective total cross sections σ ATLAS,I = (82.2 ± 8.3 ± 6.3) pb and σ ATLAS,II = (78.3 ± 9.2 ± 6.6) pb are both compatible with our LO PRA predictions σ PRA ATLAS,I = 133.6 +89.6 −52.2 pb and σ PRA ATLAS,II = 105.2 +73.8 −41.6 pb. Their respective p 2T , p ψψ T , and m ψψ distributions are compared with our LO PRA predictions in Fig. 3. We find fairly good agreement for the p 2T and p ψψ T distributions, especially in region II, with regard to both normalization and line shape. In particular, the predictions faithfully reproduce the peaks of the measured p ψψ T distributions. As for the m ψψ distributions, there is decent agreement for m ψψ 40 GeV, while the predictions significantly undershoot the data in the upper m ψψ bins, as in the CMS case above.
As for LT dominance, we find this contribution to exceed 90% of the total prediction in the upmost m ψψ bin of region-I, where the maximum allowed value of |Y | is 3.15, while this happens in the upper four bins in region II, where the maximum allowed value of |Y | is 4.2.
Although the LO CPM relationship m ψψ = 2[(2m c ) 2 + (p J/ψ T ) 2 ] 1/2 cosh(|Y |/2) [13] is evaded by the PRA, detailed analysis of the m ψψ and |Y | distributions reveals that they are strongly correlated in the large-m ψψ and -|Y | regions. At LO, only the LT subprocesses involve t-channel gluon exchanges, and thus make up the bulk of the cross section at large values of m ψψ and |Y |, as demonstrated above. This will generate large logarithmic corrections of the type (α s ln |s/t|) n in higher orders, which can be efficiently included via BFKL resummation.
The appearance of α s (µ 2 ) in the exponent of G in Eq. (4) may produce a potentially large theoretical uncertainty. Several approaches have been proposed to remedy this. As frequently done [32], we adopt here the one [33] based on a non-Abelian physical renormalization scheme choice in connection with Brodsky-Lepage-Mackenzie (BLM) optimal scale setting [34], which yields µ 2 BLM = µ 2 0 e CBLM , where µ 0 is a reference scale and sch depends on the renormalization scheme and, in the momentum subtraction (MOM) scheme advocated in Ref. [33], takes the value −1.281. We choose µ 0 = ξ[|k 1T ||k 2T |] 1/2 and vary ξ from 1/2 to 2 about its default value 1 to estimate the residual unphysical-scale uncertainty in G.
In our numerical analysis, we evaluate Eq. (1) with Eq. (2), which we call dσ BFKL , for each value of l separately. We find that only the l = 0 mode contributes to the |Y | distribution, while all l modes contribute to the m ψψ one. However, the series in |l| rapidly converges. For example, in the upmost m ψψ bin of CMS, the |l| = 2 to |l| = 0 ratio is below 2% for all LT subprocesses. We may, therefore, safely neglect all modes with |l| > 2.
We merge the full LO PRA calculation dσ PRA , appropriate in the small-|Y | region, and the LL-resummed LT contribution dσ BFKL , appropriate in the large-|Y | region, as where the asymptotic term dσ BFKL,0 , which is obtained from dσ BFKL by substituting Eq. (4) with its initial condition, is to avoid double counting. Equation (6) smoothly interpolates from dσ PRA at perturbatively small |Y | values to dσ BFKL at asymptotically large |Y | values.
The BFKL-improved PRA predictions thus evaluated are also included in Figs. 2 and 3. Their uncertainties are obtained by combining the PRA and BFKL ones in quadrature. They exceed the PRA uncertainties only moderately, which indicates that the scale uncertainties in Eq. (4) are well under control in the MOM-BLM approach [33]. In the p 2T and p ψψ T distributions and in the lowest few bins of the m ψψ and |Y | distributions, the BFKL resummation effects are so insignificant that we refrain from displaying the BFKL-improved results. On the other hand, these effects are truly significant in the upper m ψψ and |Y | bins, where they may even double the pure PRA results, so as to improve the description of the CMS [10] and ATLAS [11] data. In the latter case, even agreement is reached in some medium m ψψ bins. Unfortunately, large gaps remain in the upmost |Y | bin of CMS and the upmost few m ψψ bins of CMS and ATLAS.
While the remaining discrepancies may be attributed to the DPS mechanism, it would be premature to do so before the complete NLO NRQCD corrections are available, for the following reasons. Firstly, conventional NRQCD factorization is known to break down at NLO for double P -wave channels due to a new type of infrared divergence [17]. The quantitative influence of this is presently unclear. Secondly, the NLO NRQCD corrections to the NLT subprocesses can be quite sizable. While t-channel gluon exchange type diagrams, which contribute at leading power in the large-|Y | region, arise for LT subprocesses already at LO [see Fig. 1   1,2 channels, the type of diagrams in Fig. 1(c) form a gauge invariant subset, but not for the 3 S 1 )cc( 1 S [8] 0 )g subprocess that our MRK approximation reproduces the exact result for the t-channel gluon exchange type diagrams in the upmost CMS |Y | bin within a factor of 1.2. In this way, we find that such partial NLO (NLO * ) results for the individual (m, n) channels among the NLT subprocesses can be up to a hundred times larger than the LO PRA results for these channels in the upper |Y | and m ψψ bins. The effect of adding the total NLO * NLT contribution on top of the central LO NRQCD prediction in the PRA with BFKL resummation is shown for the upmost m ψψ and |Y | bins in Figs. 2 and 3. In Fig. 2(c), this amounts to 45% and 16% for direct and prompt production, respectively. The total NLO * NLT contributions will, in turn, be enhanced by BFKL resummation, which we leave for future work.
To summarize, we have pushed the NRQCD factorization approach to double prompt J/ψ hadroproduction beyond LO in two important ways, so as to consolidate the theoretical basis for the interpretation of Tevatron and LHC data and, ultimately, for meaningful extractions of the DPS key parameter σ eff . On the one hand, we have incorporated multiple gluon radiation off the initial state via the PRA, which, unlike other k T factorization approaches frequently used in the literature, ensures for the SDCs to be manifestly gauge invariant, infrared safe, and devoid of artificial kinematic cuts. On the other hand, we have resummed, via BFKL evolution in |Y |, the LLs of the form (α s ln |s/t|) n arising from t-channel gluon exchanges in the LT subprocesses [see Fig. 1(b)], which would otherwise inevitably invalidate the fixed-order treatment at large |Y | and m ψψ values. Using updated CO LDMEs from a LO PRA fit to LHC data of prompt single J/ψ hadroproduction [30] (see Table I), we have predicted total cross sections as well as p 2T , p ψψ T , m ψψ , and |Y | distributions of prompt double J/ψ hadroproduction and confronted them with CMS [10] and ATLAS [11] data. Passing from the CPM to the PRA has enabled p 2T and p ψψ T distributions, which turned out to nicely agree with the CMS and ATLAS data, and has enhanced the m ψψ and |Y | distributions in the lowest few bins so as to nicely match the CMS and ATLAS data there. On the other hand, the desirable enhancements in the upper bins have failed to appear. However, performing BKFL resummation for the LT processes on top of the PRA has greatly improved the situation there, by generating K factors of up to 2 in the upmost m ψψ and |Y | bins. Using the gauge invariant MRK formalism, we have also estimated the NLO NLT contributions there. The combined BFKL and NLO NLT effects significantly reduce the gaps between SPS predictions and data in the upper m ψψ and |Y | bins, leaving less room for DPS and, in turn, bringing the anomalously small value of σ eff [12] up, in line with other determinations, as expected if σ eff is universal.