Double prompt $J/\psi$ hadroproduction in the parton Reggeization approach with high-energy resummation

TL;DR

This paper investigates double prompt $J/$ hadroproduction using the parton Reggeization approach, providing predictions that largely agree with experimental data and highlighting the importance of high-energy resummation in certain kinematic regions.

Contribution

It introduces the application of the parton Reggeization approach with high-energy resummation to double $J/$ production, improving theoretical predictions in high-energy collider physics.

Findings

Predictions agree with ATLAS and CMS data at leading order.

High-energy resummation enhances cross sections in large invariant mass and rapidity separation regions.

Resummation improves the description of data in previously underestimated regions.

Abstract

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.

Tables

Table 1. Table 1: Adopted values of LO NRQCD LDMEs in units of GeV 3 .

$⟨{𝓞}^{𝑱/𝝍}\left({}_{}{}^{\mathrm{𝟑}}𝑺_{\mathrm{𝟏}}^{\left[\mathrm{𝟏}\right]}\right)⟩$	$⟨{𝓞}^{𝑱/𝝍}\left({}_{}{}^{\mathrm{𝟏}}𝑺_{\mathrm{𝟎}}^{\left[\mathrm{𝟖}\right]}\right)⟩$	$⟨{𝓞}^{𝑱/𝝍}\left({}_{}{}^{\mathrm{𝟑}}𝑺_{\mathrm{𝟏}}^{\left[\mathrm{𝟖}\right]}\right)⟩$	$⟨{𝓞}^{{𝝍}^{\prime }}\left({}_{}{}^{\mathrm{𝟑}}𝑺_{\mathrm{𝟏}}^{\left[\mathrm{𝟏}\right]}\right)⟩$	$⟨{𝓞}^{{𝝍}^{\prime }}\left({}_{}{}^{\mathrm{𝟏}}𝑺_{\mathrm{𝟎}}^{\left[\mathrm{𝟖}\right]}\right)⟩$	$⟨{𝓞}^{{𝝍}^{\prime }}\left({}_{}{}^{\mathrm{𝟑}}𝑺_{\mathrm{𝟏}}^{\left[\mathrm{𝟖}\right]}\right)⟩$	$\frac{⟨{𝓞}^{{𝝌}_{𝒄\mathrm{𝟎}}}\left({}_{}{}^{\mathrm{𝟑}}𝑷_{\mathrm{𝟎}}^{\left[\mathrm{𝟏}\right]}\right)⟩}{{𝒎}_{𝒄}^{\mathrm{𝟐}}}$	$⟨{𝓞}^{{𝝌}_{𝒄\mathrm{𝟎}}}\left({}_{}{}^{\mathrm{𝟑}}𝑺_{\mathrm{𝟏}}^{\left[\mathrm{𝟖}\right]}\right)⟩$	$\frac{⟨{𝓞}^{𝑱/𝝍\left({𝝍}^{\prime }\right)}\left({}_{}{}^{\mathrm{𝟑}}𝑷_{𝑱}^{\left[\mathrm{𝟖}\right]}\right)⟩}{{𝒎}_{𝒄}^{\mathrm{𝟐}}}$
1.16	$3.61\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟐}}$	$1.25\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟑}}$	$0.76$	$2.19\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟐}}$	$3.41\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟒}}$	$4.77\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟐}}$	$5.29\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟒}}$	0

Equations

$$d{\sigma }^{\mathrm{PRA}}(AB\to 2J/\psi +X)=\sum _{m,n,H_1,H_2}\int \frac{d{x}_1}{x_1}\int \frac{d^2{\mathbf{k}}_{1T}}{\pi }$$

$$\frac{d{\hat{\sigma }}_{mn}^{\mathrm{BFKL}}}{d{\mathbf{p}}_{1T}d{y}_{1}d{\mathbf{p}}_{2T}d{y}_2}=\int \frac{d^2{\mathbf{q}}_T{d}^2{\mathbf{q}}_T^{\prime }}{(4\pi {)}^2{S}^2{x}_1{x}_2}$$

$$d{\sigma }^{\mathrm{PRA}+\mathrm{BFKL}}=d{\sigma }^{\mathrm{PRA}}+d{\sigma }^{\mathrm{BFKL}}-d{\sigma }^{\mathrm{BFKL},0},$$

Full text

††preprint: DESY 19–111                                                                                                                         ISSN 0418–9833              ††preprint: June 2019                                                                                                                                                              

Despite concerted experimental and theoretical endeavors ever since the discovery of $J/\psi$ meson more than four decades ago, its production mechanism has remained mysterious; for a recent review, see Ref. Brambilla:2010cs . The factorization approach Bodwin:1994jh to nonrelativistic QCD (NRQCD) Caswell:1985ui endowed with velocity scaling rules Lepage:1992tx for the long-distance matrix elements (LDMEs), which is by far the most acceptable candidate theory for heavy-quarkonium production and decay and has been elaborated at next-to-leading order (NLO), has been challenged by the long-standing $J/\psi$ polarization puzzle Butenschoen:2012px and by the inadequate description of $\eta_{c}$ hadroproduction data with $J/\psi$ LDMEs converted via heavy-quark spin symmetry Butenschoen:2014dra . The high flux of incoming partons at LHC allows us to study $J/\psi$ production more thoroughly, also in association with other charmonia, bottomonia, $W$ or $Z$ bosons, so as to pin down the $J/\psi$ production mechanism. Among these production processes, double $J/\psi$ hadroproduction is of special interest because $J/\psi$ formation takes place there twice, making this particularly sensitive to the nonperturbative aspects of NRQCD Barger:1995vx . 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 $\sigma_{\mathrm{eff}}$ Kom:2011bd .

In recent years, double prompt $J/\psi$ hadroproduction has been measured extensively by the LHCb Aaij:2011yc , CMS Khachatryan:2014iia , and ATLAS Aaboud:2016fzt Collaborations at the CERN LHC, and by the D0 Collaboration Abazov:2014qba at the FNAL Tevatron. On the theoretical side, the complete NRQCD results at leading order (LO) have been obtained recently He:2015qya . For some channels, the relativistic corrections Li:2013csa and NLO QCD corrections Sun:2014gca 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_{\psi\psi}$ or rapidity separation $|Y|=|y_{1}-y_{2}|$ of the $J/\psi$ pair, although the color-octet contributions, in particular those involving $t$-channel gluon exchanges, enhance the QCD-corrected Sun:2014gca color-singlet contribution of $gg\to 2c\bar{c}({}^{3}S_{1}^{[1]})$ there by more than one order of magnitude. The value of $\sigma_{\mathrm{eff}}$ extracted by fitting the DPS contribution on top of this Abazov:2014qba 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_{\psi\psi}$ and $|Y|$ bins Lansberg:2014swa .

As noticed in Ref. He:2015qya , the CMS kinematic conditions Khachatryan:2014iia render $2\to 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 He:2018hwb . Moreover, the perturbative expansion is spoiled for small values of the $J/\psi$ pair transverse momentum $p_{T}^{\psi\psi}$. The characteristic scale $\mu\sim[(4m_{c})^{2}+(p_{T}^{\psi\psi})^{2}]^{1/2}$ of the hard-scattering processes of double $J/\psi$ hadroproduction satisfies $\Lambda_{\mathrm{QCD}}\ll\mu\ll\sqrt{S}$, where $\Lambda_{\mathrm{QCD}}$ is the asymptotic scale parameter and $\sqrt{S}$ is the center-of-mass energy. We are thus accessing the high-energy Regge regime, where the NLO QCD corrections can be largely accounted for through the unintegrated parton distribution functions (unPDFs) in the parton Reggeization approach (PRA) Kniehl:2014qva based on Lipatov’s effective field theory formulated with the non-Abelian gauge invariant action Lipatov:1995 . The PRA has already been successfully applied to the interpretation of measurements of single heavy-quarkonium hadroproduction Kniehl:2006sk ; Kniehl:2006vm . In the large-$|Y|$ region, the two $J/\psi$ 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 $(\alpha_{s}\ln|s/t|)^{n}$ in the higher-order QCD corrections, where $s$ and $t$ are the Mandelstam variables of the partonic $2\to 2$ Born process. Such large logarithms can be resummed by the Balitsky-Fadin-Kuraev-Lipatov (BFKL) formalism Kuraev:1976ge . Recently, BFKL resummation has been studied for single $J/\psi$ Kotko:2019kma and $J/\psi$ plus jet Boussarie:2017oae inclusive hadroproduction. In this Letter, we will take a crucial step towards a full-fledged NLO NRQCD study of double prompt $J/\psi$ hadroproduction, by adopting the PRA and performing BFKL resummation.

Owing to the PRA and NRQCD factorization, the cross section of inclusive double prompt $J/\psi$ hadroproduction can be expressed as:

[TABLE]

where $d\hat{\sigma}_{mn}^{\mathrm{PRA}}$ is the short-distance coefficient (SDC) of the partonic subprocess $R^{+}(k_{1})R^{-}(k_{2})\to c\bar{c}(m)c\bar{c}(n)+X$, $\Phi_{R^{\pm}/A,B}(x_{1,2},t_{1,2},\mu^{2})$ are the unPDFs of the Reggeized gluons $R^{\pm}$ with four-momenta $k_{1,2}^{\mu}=x_{1,2}K_{1,2}^{\mu}+k_{1,2T}^{\mu}$ and virtualities $t_{1,2}=-k_{1,2}^{2}=\mathbf{k}_{1,2T}^{2}$, $K_{1,2}^{\mu}$ are the four-momenta of the colliding hadrons $A,B$ with light-cone components $K_{1}^{-}=K_{2}^{+}=0$ ($K_{i}^{\pm}=K_{i}^{0}\pm K_{i}^{3}$), and $\langle\overline{\mathcal{O}}^{H}(m)\rangle$ is the product of LDME $\langle\mathcal{O}^{H}(m)\rangle$ of $H=J/\psi,\chi_{cJ},\psi^{\prime}$ and the branching fraction $\mathrm{BR}(H\to J/\psi+X)$, with the understanding that $\mathrm{BR}(H\to J/\psi+X)=1$ if $H=J/\psi$. 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 $\chi_{c0}$ feed-down contribution because $\mathrm{BR}(\chi_{c0}\to J/\psi+X)=1.40\%$ Tanabashi:2018kda is so small.

Representative Feynman diagrams for the partonic subprocess $R^{+}R^{-}\to c\bar{c}(m)c\bar{c}(n)$ at LO in $\alpha_{s}$ are depicted in Figs. 1(a) and (b). By the Feynman rules of Ref. Lipatov:1995 , they come in the same topologies as those for $gg$ fusion in the CPM. To discuss the BFKL resummation effect conveniently, we divide the partonic subprocesses into three categories, according to the order in $\alpha_{s}$ where $t$-channel gluon exchanges emerge for the first time, namely (i) LO $t$-channel (LT), with $m,n={}^{1}\!S_{0}^{[8]},{}^{3}\!S_{1}^{[8]},{}^{3}\!P_{J}^{[1,8]}$; (ii) NLO $t$-channel (NLT), with $m={}^{3}\!S_{1}^{[1]}$ and $n={}^{1}\!S_{0}^{[8]},{}^{3}\!S_{1}^{[8]},{}^{3}\!P_{J}^{[1,8]}$; and (iii) NNLO $t$-channel (NNLT), with $m,n={}^{3}\!S_{1}^{[1]}$; see Figs. 1(b)–(d).

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 Baranov:2015cle , 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}\to 0$, we recover the CPM formulas He:2015qya , which constitutes yet another nontrivial check. In the numerical analysis, we adopt the Kimber-Martin-Ryskin scheme Kimber:2001sc to generate the unPDFs from the LO CPM PDFs of Ref. Martin:2009iq , which come with $\alpha_{s}(M_{Z})=0.13939$. We choose the renormalization and factorization scales to be $\mu_{r}=\mu_{f}=\xi[(4m_{c})^{2}+\bar{p}_{T}^{2}]^{1/2}$, where $m_{c}=1.5~{}\mathrm{GeV}$, $\bar{p}_{T}=(p_{T}^{H_{1}}+p_{T}^{H_{2}})/2$, and $\xi$ is varied between $1/2$ and 2 about its default value 1 to estimate the scale uncertainty. In the case of feed-down from $H=\chi_{c_{1}},\chi_{c_{2}},\psi^{\prime}$, we put $p_{T}^{J/\psi}=p_{T}^{H}M_{J/\psi}/M_{H}$, which is a good approximation because $p_{T}^{J/\psi}\gg M_{H}-M_{J/\psi}$ Ma:2010vd . For the $J/\psi$, $\chi_{cJ}$, and $\psi^{\prime}$ LDMEs, we use the values specified in Table 1. The color-singlet results have been derived from the Buchmüller-Tye potential in Ref. Eichten:1995ch . The color-octet results have been fitted to LHC data of inclusive single charmonium hadroproduction ATLAS:2014ala in the very theoretical framework described above; they supersede pre-LHC results Kniehl:2006sk . For $H=J/\psi,\psi^{\prime}$, there is a strong correlation between $\langle\mathcal{O}^{H}({}^{1}\!S_{0}^{[8]})\rangle$ and $\langle\mathcal{O}^{H}({}^{3}\!P_{0}^{[8]})\rangle$, so that only a linear combination of them can be determined. We may thus put $\langle\mathcal{O}^{H}({}^{3}\!P_{0}^{[8]})\rangle=0$. We have checked that the theoretical uncertainties in our predictions for double prompt $J/\psi$ hadroproduction due to this freedom are negligible. All other input parameters are adopted from Ref. Tanabashi:2018kda .

The CMS data of prompt double $J/\psi$ hadroproduction were taken at $\sqrt{S}=7~{}\mathrm{TeV}$ requiring for each $J/\psi$ 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. Khachatryan:2014iia . The CMS total cross section $\sigma_{\mathrm{CMS}}=(1.49\pm 0.07\pm 0.13)~{}\mathrm{nb}$ agrees with our LO NRQCD prediction $\sigma_{\mathrm{CMS}}^{\mathrm{PRA}}=1.68^{+1.32}_{-0.78}{}~{}\mathrm{nb}$ within errors. The CMS $p_{T}^{\psi\psi}$, $m_{\psi\psi}$, and $|Y|$ distributions are compared with our predictions in Fig. 2(a)–(c), respectively. There is generally very good agreement, except for the upper two $m_{\psi\psi}$ 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}^{\psi\psi}=0$ at LO in the latter case, but it is also significant for the $m_{\psi\psi}$ 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. He:2015qya . 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.

The unPDF uncertainty may be assessed from Fig. 2(a), which also shows the evaluations using the set produced from our default PDF set Martin:2009iq as described in Ref. Blumlein:1995eu and the set of Ref. Jung:2000hk . The latter rapidly falls off with increasing $p_{T}^{\psi\psi}$ and significantly undershoots the data in the upper $p_{T}^{\psi\psi}$ bins. This is in line with Ref. Maciula:2018bex , where the unPDFs of Ref. Jung:2000hk were found to yield a poor description of LHCb data of single prompt $J/\psi$ production Aaij:2011jh at large $p_{T}^{\psi}$. This opens a novel perspective to constrain unPDFs. To estimate the LDME uncertainty, we repeat the unresummed LO PRA evaluation in the upmost $|Y|$ bin of Fig. 2(c), which is most sensitive to the color-octet LDMEs, using in turn the NLO CPM sets of Refs. Butenschoen:2011yh ; Gong:2012ug , albeit this is slightly inconsistent. Since Ref. Butenschoen:2011yh does not provide $\chi_{cJ}$ and $\psi^{\prime}$ LDMEs, we use those of Ref. Gong:2012ug also here. We thus find an enhancement by 57% and a reduction by 7% w.r.t. our default result, respectively. For a more detailed LDME analysis, see Ref. Lansberg:2019fgm .

ATLAS took their data at $\sqrt{S}=8~{}\mathrm{TeV}$ imposing the acceptance cuts $p_{T}>8.5~{}\mathrm{GeV}$ and $|y|<2.1$ on each $J/\psi$ meson Aaboud:2016fzt . They separately studied the central (I) and forward (II) $y$ regions of the subleading $J/\psi$ meson ($J/\psi_{2}$), with $p_{2T}<p_{1T}$, namely $|y_{2}|<1.05$ and $1.05<|y_{2}|<2.1$. Their respective total cross sections $\sigma_{\mathrm{ATLAS,I}}=(82.2\pm 8.3\pm 6.3)~{}\mathrm{pb}$ and $\sigma_{\mathrm{ATLAS,II}}=(78.3\pm 9.2\pm 6.6)~{}\mathrm{pb}$ are both compatible with our LO PRA predictions $\sigma_{\mathrm{ATLAS,I}}^{\mathrm{PRA}}=133.6^{+89.6}_{-52.2}{}~{}\mathrm{pb}$ and $\sigma_{\mathrm{ATLAS,II}}^{\mathrm{PRA}}=105.2^{+73.8}_{-41.6}{}~{}\mathrm{pb}$. Their respective $p_{2T}$, $p_{T}^{\psi\psi}$, and $m_{\psi\psi}$ distributions are compared with our LO PRA predictions in Fig. 4. We find fairly good agreement for the $p_{2T}$ and $p_{T}^{\psi\psi}$ 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}^{\psi\psi}$ distributions. As for the $m_{\psi\psi}$ distributions, there is decent agreement for $m_{\psi\psi}\lesssim 40~{}\mathrm{GeV}$, while the predictions significantly undershoot the data in the upper $m_{\psi\psi}$ bins, as in the CMS case above.

Also the LHCb Aaij:2011yc and D0 Abazov:2014qba measurements reasonably agree with our PRA predictions. The LHCb Aaij:2011yc data at $\sqrt{S}=7~{}\mathrm{TeV}$, with acceptance cuts $p_{T}<10~{}\mathrm{GeV}$ and $2<y<4.5$ on each $J/\psi$ meson, yield $\chi^{2}/\mathrm{d.o.f.}=12.8/5$ with respect to our LO PRA prediction with $\xi=1$ in the perturbatively safe region $m_{\psi\psi}>9~{}\mathrm{GeV}$ He:2015qya , with experimental and theoretical errors overlapping. The fiducial cross section $\sigma_{\mathrm{SPS}}=70\pm 6(\mathrm{stat})\pm 22(\mathrm{syst})~{}\mathrm{fb}$ determined by D0 Abazov:2014qba at $\sqrt{S}=7~{}\mathrm{TeV}$, with cuts $p_{T}<10~{}\mathrm{GeV}$ and $2<y<4.5$, nicely agrees with our central LO PRA prediction $\sigma_{\mathrm{SPS}}^{\mathrm{PRA}}=81.1~{}\mathrm{fb}$, where we have included the reduction factor due the acceptance cuts on the decay muons Abazov:2014qba determined in Ref. Qiao:2012wc .

Although the LO CPM relationship $m_{\psi\psi}=2[(2m_{c})^{2}+(p_{T}^{J/\psi})^{2}]^{1/2}\cosh(|Y|/2)$ He:2015qya is evaded by the PRA, detailed analysis of the $m_{\psi\psi}$ and $|Y|$ distributions reveals that they are strongly correlated in the large-$m_{\psi\psi}$ and -$|Y|$ regions. As expected, the LT contributions greatly dominate there. Specifically, they make up about 90% or more of the total PRA predictions in the upmost $m_{\psi\psi}$ bin in Fig. 2(b), the upmost $|Y|$ bin in Figs. 2(c), the upmost $m_{\psi\psi}$ bin in Fig. 4(c), and the upper four $m_{\psi\psi}$ bins in Fig. 4(f). $t$-channel gluon exchanges, which appear in the LT contributions already at LO, generate large logarithmic corrections of the type $(\alpha_{s}\ln|s/t|)^{n}$ in higher orders, which can be efficiently included via BFKL resummation.

BFKL resummation in the leading-logarithmic (LL) approximation is implemented by replacing the SDCs in Eq. (1) with

[TABLE]

where $\Psi_{+}^{(m)}(\mathbf{k}_{1T},\mathbf{p}_{1T},y_{1})$ and $\Psi_{-}^{(n)}(\mathbf{k}_{2T},\mathbf{p}_{2T},y_{2})$ are the impact factors describing the partonic subprocesses $R^{+}(\mathbf{k}_{1T})R^{-}(\mathbf{q}_{T})\to[c\bar{c}(m)](\mathbf{p}_{1T})$ and $R^{-}(\mathbf{k}_{2T})R^{+}(-\mathbf{q}_{T}^{\prime})\to[c\bar{c}(n)](\mathbf{p}_{2T})$, obtained from the appropriate $2\to 1$ PRA matrix elements in Ref. Kniehl:2006sk as explained in Ref. Kovchegov:2012mbw , and $G(\mathbf{q}_{T},\mathbf{q}^{\prime}_{T},Y)$ is the BFKL Green function given by Eq. (3.80) in Ref. Kovchegov:2012mbw , generated from the initial condition $G(\mathbf{q}_{T},\mathbf{q}^{\prime}_{T},0)=\delta^{(2)}(\mathbf{q}_{T}-\mathbf{q}^{\prime}_{T})$ via LL BFKL evolution in $Y$ Kuraev:1976ge ; see Fig. 3 for a schematic representation of Eq. (2). The resulting hadronic cross section is denoted as $d\sigma^{\mathrm{BFKL}}$.

$G(\mathbf{q}_{T},\mathbf{q}^{\prime}_{T},Y)$ depends exponentially on $\alpha_{s}(\mu^{2})$, which may produce a potentially large theoretical uncertainty. Several approaches have been proposed to remedy this. As frequently done Ducloue:2013bva , we adopt here the one Brodsky:1998kn based on a non-Abelian physical renormalization scheme choice in connection with Brodsky-Lepage-Mackenzie optimal scale setting Brodsky:1982gc . There remains a reference scale $\mu_{0}$. We choose $\mu_{0}=\xi[|\mathbf{k}_{1T}||\mathbf{k}_{2T}|]^{1/2}$ and vary $\xi$ from $1/2$ to 2 about its default value 1 to estimate the residual scale uncertainty in $G(\mathbf{q}_{T},\mathbf{q}^{\prime}_{T},Y)$.

We merge the full LO PRA calculation $d\sigma^{\mathrm{PRA}}$, appropriate in the small-$|Y|$ region, and the LL-resummed LT contribution $d\sigma^{\mathrm{BFKL}}$, appropriate in the large-$|Y|$ region, as

[TABLE]

where the asymptotic term $d\sigma^{\mathrm{BFKL},0}$, which is obtained from $d\sigma^{\mathrm{BFKL}}$ by replacing $G(\mathbf{q}_{T},\mathbf{q}^{\prime}_{T},Y)$ with $G(\mathbf{q}_{T},\mathbf{q}^{\prime}_{T},0)$, is to avoid double counting. Equation (3) smoothly interpolates from $d\sigma^{\mathrm{PRA}}$ at small $|Y|$ values to $d\sigma^{\mathrm{BFKL}}$ at large $|Y|$ values.

The BFKL-improved PRA predictions thus evaluated are also included in Figs. 2 and 4. 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 $G(\mathbf{q}_{T},\mathbf{q}^{\prime}_{T},Y)$ are well under control. In the $p_{2T}$ and $p_{T}^{\psi\psi}$ distributions and in the lowest few bins of the $m_{\psi\psi}$ 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 significant in the upper $m_{\psi\psi}$ and $|Y|$ bins, where they may even double the pure PRA results, so as to reduce the shortfall with respect to the CMS Khachatryan:2014iia and ATLAS Aaboud:2016fzt data. In the latter case, even agreement is reached in some medium $m_{\psi\psi}$ bins. However, large gaps remain in the upmost $|Y|$ bin of CMS and the upmost few $m_{\psi\psi}$ bins of CMS and ATLAS, to be explained by DPS. Contrary to naïve expectations, the optimal scale turns out to be larger than $\mu_{0}$, so that using the latter instead leads to an enhancement of the BFKL-improved results, by factors of 1.0, 1.0, 1.1, 1.5, and 4.1 in the second to sixth $|Y|$ bin of CMS.

Quantitative extractions of DPS contributions from the remaining discrepancies are likely to be meaningful only after 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 He:2018hwb . The quantitative influence of this is presently unclear. Secondly, the NLO NRQCD corrections to the NLT subprocesses can be quite sizable because their $\mathcal{O}(\alpha_{s})$ suppression is expected to be compensated by the relatively large values of the color-singlet LDMEs $\langle\mathcal{O}^{J/\psi}({}^{3}\!S_{1}^{[1]})\rangle$ and $\langle\mathcal{O}^{\psi^{\prime}}({}^{3}\!S_{1}^{[1]})\rangle$. For the ${}^{3}\!S_{1}^{[1]}+{}^{1}\!S_{0}^{[8]}/{}^{3}\!P_{1,2}^{[1]}$ channels, the type of diagrams in Fig. 1(c) form a gauge invariant subset, but not for the ${}^{3}\!S_{1}^{[1]}+{}^{3}\!S_{1}^{[8]}$ channel because of $g\to c\bar{c}({}^{3}\!S_{1}^{[8]})$ formation. Their leading large-$|Y|$ contributions can be estimated via the gauge invariant MRK-asymptotic formalism, already used to evaluate $d\sigma^{\mathrm{BFKL},0}$ for the LT subprocesses in Eq. (3). We have checked for the $R^{+}R^{-}\to c\bar{c}({}^{3}\!S_{1}^{[1]})c\bar{c}({}^{1}\!S_{0}^{[8]})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_{\psi\psi}$ bins. The effect of adding the total ${\rm NLO}^{\ast}$ NLT contribution on top of the central LO NRQCD prediction in the PRA with BFKL resummation is shown for the upmost $m_{\psi\psi}$ and $|Y|$ bins in Figs. 2 and 4. In Fig. 2(c), this amounts to 45% and 16% for direct and prompt production, respectively. The total ${\rm NLO}^{\ast}$ 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/\psi$ hadroproduction beyond LO in two important ways. 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 Baranov:2015cle , 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 $(\alpha_{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_{\psi\psi}$ values. This consolidates the theoretical basis for meaningful extractions of the DPS key parameter $\sigma_{\mathrm{eff}}$.