Prompt ${J/\psi}$-pair production at the LHC: impact of loop-induced contributions and of the colour-octet mechanism

Prompt double-$J/\psi$ production at high-energy hadron colliders can be considered as a golden channel to probe double parton scatterings (DPS) --in particular to study gluon-gluon correlations inside the proton-- and, at the same time, to measure the distribution of linearly-polarised gluons inside the proton. Such studies however require a good control of both single and DPS in the respective regions where they are carried out. In this context, we have critically examined two mechanisms of single parton scatterings (SPS) that may be kinematically enhanced where DPS are thought to be dominant, even though they are either at higher orders in the strong-coupling or velocity expansion. First, we have considered a gauge-invariant and infrared-safe subset of the loop-induced contribution via Colour-Singlet (CS) transitions. We have found it to become the leading CS SPS contributions at large rapidity separation, yet too small to account for the data without invoking the presence of DPS yields. Second, we have surveyed the possible Colour-Octet (CO) contributions using both old and up-to-date non-perturbative long distance matrix elements (LDMEs). We have found that the pure CO yields crucially depend on the LDMEs. Among all the LDMEs we used, only two result into a visible modification of the NRQCD (CS+CO) yield, but only in two kinematical distributions measured by ATLAS, those of the rapidity separation and of the pair invariant mass. These modifications however do not impact the control region used for their DPS study.


Introduction
The role of multiple parton interactions in proton-proton collisions is believed to become increasingly important when one explores the energy frontier in particle physics. As such, the relevance in LHC observables of two simultaneous hard scatterings, usually referred to as Double Parton Scatterings (DPS), has attracted much attention in the last decade with significant theory advances related to perturbative QCD [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15]. Since DPS are highertwist effects in total cross sections compared to the conventional single parton scatterings (SPS), quantitative studies of DPS remain challenging though not impossible both on the theoretical and experimental sides. These are particularly interesting since they provide us with means to study parton correlations inside the proton (see e.g. [10,16,17]).
Among the possible hard probes of DPS at high-energy hadron colliders, the associated production of quarkonia (see [18] for an exhaustive review) provides unique opportunities to measure DPS in gluon-induced reactions thus to study gluon-gluon correlations in the proton. Numerous measurements of quarkonium associated processes have been performed at the Tevatron and the LHC. They can mainly be categorised as di-quarkonium production (J/ψ + J/ψ [19][20][21][22][23], J/ψ + Υ [24], Υ + Υ [25]), associated production with a vector boson (J/ψ + W ± [26], J/ψ + Z [27]) or with another heavy quark (J/ψ+open charm [28], Υ+open charm [29]). All these measurements cover different kinematical regions with different momentum transfers in the hard scattering. Their theoretical analysis is highly non-trivial, which has triggered many theoretical studies in the recent years . Very recently, the first calculation of triple-J/ψ production showed that it can help us probe both DPS and triple parton scatterings (TPS) [51].
In this context, we focus in this paper on the di-J/ψ case with the aim to improve the existing perturbative QCD calculations for the SPS. To do so, we consider higher-order corrections in both the strong coupling constant, α S , and the heavy-quark velocity, v. First, we study the impact of a gauge-invariant and infrared-safe subset of loop-induced (LI) contributions. Our analysis follows the lines of a similar study for J/ψ + Υ production [42]. Such contributions appear at next-to-next-to-leading order (NNLO) in α S but could be enhanced at large rapidity differences and high invariant masses of the J/ψ pair because of the presence of topologies with double t-channel gluon exchanges between both charmanticharm quark lines. Second, we perform a comprehensive survey of the impact of the colour-octet (CO) contributions in three kinematical domains covered by the existing LHC measurements [21][22][23] considering the various existing fits of the non-perturbative CO long-distance matrix elements (LDMEs).
In order to disentangle DPS from SPS in observables where two particles are observed, one usually relies on the analysis of specific kinematical dependences which are believed to be drastically different in both samples. Common choices of variables are the azimuthal and the rapidity separations between both observed particles, ∆φ and ∆y. The DPS contribution, coming from two a priori independent parton scatterings, is expected to be flatter than the SPS one in both distributions.
For double-J/ψ studies, the analysis of the ∆y(J/ψ, J/ψ) distributions should be preferred compared to that of ∆φ(J/ψ, J/ψ) since the ∆y(J/ψ, J/ψ) distribution of the SPS yield is much less affected by possible non-perturbative intrinsic k T of the colliding gluons [33] than the ∆φ(J/ψ, J/ψ) one, which can become as flat as the DPS ones in some cases. In general, one expects the DPS fraction to be the largest at large |∆y(J/ψ, J/ψ)|. A precise determination of the DPS yield therefore requires a good knowledge of the SPS in this region. Both the LI and CO topologies with t-channel-gluon exchanges could result into a flat dσ d∆y like in the J/ψ + Υ case [42]. Assuming α S ∼ v 2 , the colour-singlet (CS) LI contribution should be of the same magnitude as the leading order (LO) CO contribution (yet both smaller that the bulk of the CS yield in the absence of the possible kinematical enhancement which we are after here). According to the NRQCD velocity scaling rules [52], the former one is indeed O(α 6 S v 3 ) while the latter one is O(α 4 S v 7 ). This justifies why we consider both of them in this study.
This article is organised as follows. In section 2, we first quickly review the existing LHC measurements used in our comparisons 1 . Then, we discuss our theory framework in section 3. Section 4 gathers our discussion of the impact of the inclusion LI CS contribution and section 5 comprises a comprehensive analysis of complete LO CO contribution. The appendix A collects additional plots relevant for further theory-data comparisons.

kinematical variables
We start by introducing the kinematical variables relevant for di-quarkonium production.
On the experimental side, the second LHCb analysis [23] bears on the largest set of the kinematical variables whose distribution is used for comparisons between the experimental measurements and the theoretical calculations. Since some of these variables may not be very common, we summarise the description of their names or labels in Table 1. In particular, the transverse momentum asymmetry is defined as where J/ψ 1 and J/ψ 2 are respectively denoted as the first and second hardest J/ψ with ordered in the transverse momentum.

Available data sets
Four LHC studies of double prompt J/ψ production have so far been performed [19,[21][22][23]. LHCb performed two measurements in the same kinematical region, one at √ s = 7 TeV and another at √ s = 13 TeV; we will focus on the latter which is more precise [23]. The various kinematical cuts used in the ATLAS, CMS and LHCb analyses are summarised in Table. 2 along with the corresponding centre-of-mass energy √ s . It is useful to note that due to the different trigger and acceptance constraints on the ATLAS, CMS and LHCb data taking, the 3 samples cover complementary domains in P T and y. In particular, ATLAS [22] imposes the largest P T (J/ψ) cut (as large as 8.5 GeV), while LHCb [23] does not impose any lower P T cut on the observed J/ψ. As such, LHCb events are mostly located at low P T (J/ψ). CMS [21] imposes varying cuts from P T (J/ψ) > 4.5 GeV to P T (J/ψ) > 6.5 GeV depending on the rapidity. Moreover, LHCb can only detect forward particles whereas ATLAS/CMS have a generally larger rapidity coverage but in the central-rapidity region.
In section 5, we will discuss how these kinematical coverages can be relevant to determine the proper CO LDMEs.

Theory framework
In this section, we briefly address some specificities of our theoretical computations, which however remain very standard.

Intrinsic initial-k T smearing
An important effect for an accurate description of double-J/ψ hadroproduction is known to be the smearing of the kinematics arising from the intrinsic k T of the gluons [53]. It is in principle a non-perturbative effect which cannot properly be accounted for by the collinear factorisation. In fact, double-J/ψ production can provide new insights in the transverse dynamics of the gluons as it was shown [48] using the transverse-momentum dependent (TMD) factorisation. Clearly, a collinear computation is not meant to encapsulate such effects. As a makeshift, we simply rely on an empirical procedure to deal with them which we believe to be sufficient for our phenomenological purpose. In particular, the whole k T smearing is assumed to be factorised out by where the phase-space mapping Φ → Φ k T is determined by boosting the whole event according to the generated transverse-momentum imbalance | − → k T | = k T with a uniform distribution of the azimuthal angle in the transverse plane. Other forms are of course possible. In the present study, we assume k T to be the same for all three experimental coverages and fix its value to be 3.0 GeV. The distributions with other k T values are also not shown but can easily be obtained with the help of HELAC-Onia [54,55]. In fact, the NLO distributions with k T = 0.5 GeV and 2.0 GeV can be found in a theory-data comparison made by LHCb [23]. The k T -smearing effect is only visible for the P T (J/ψ + J/ψ), ∆φ(J/ψ, J/ψ) and A T (J/ψ, J/ψ) distributions.

Parameters entering our calculations
We now quickly describe our set-up for the present calculations before discussing the numerical results. We have fixed the charm quark mass to be 1.5 GeV and only the light u, d, s (anti)quarks and the gluons are allowed in the initial states. In order to be compatible with our previous NLO calculations, we have used the NLO parton-distribution functions (PDFs) CTEQ6M [56] for the calculations in the ATLAS and CMS acceptances and NNPDF3.0 [57] for those in the LHCb acceptance. We have explicitly checked that the PDF dependence is less than 20% and is thus a minor source of uncertainty compared to the (dominant) scale uncertainty which we discuss below. The missing higher-order terms in α S are estimated in the usual way by independently varying the factorisation and renormalisation scales as where the wave function at the origin R H QQ (0) can be determined by solving the Schrödinger equation with a given QCD potential. We will use the numerical values R J/ψ (0) 2 = 0.8 GeV 3 and R ψ(2S) (0) 2 = 0.5 GeV 3 derived in Ref. [58] using the QCD-motivated Buchmüller-and-Tye potential [59]. For the CS SPS yield, the feed-down contribution from the ψ(2S) decays is as large as the direct double J/ψ production. In practice, we take it to be equal to 2. It is thus mandatory to take it into account.

Colour-singlet contributions: partial loop-induced corrections
In principle, considering the square of a one-loop amplitude by itself should give divergent results from both the infrared and ultraviolet regions. Such one-loop amplitudes squared are part of the NNLO contributions, at O(α 6 S ) in the case of double J/ψ hadroproduction. The cancellation of the aforementioned infrared divergences would be achieved as usual by considering two-loop, one-loop single-real-emission and double-real-emission amplitudes. Such a computation is obviously beyond the scope of this study -it is not even available for single J/ψ.
However, a subset of such one-loop diagrams, restricted to the sole topologies with two separate charm-quark lines forming each a J/ψ, happens to be free of any divergence and is, in addition, gauge invariant. Correspondingly, the possible double-real emissions which could develop infrared divergences do not contribute when one of the external gluon becomes soft. This is akin to the absence of any infrared divergences at P T → 0 for gg → J/ψg. Such a subset is in fact that of the LI contribution to pp → J/ψ +Υ considered in Ref. [42] The square of the amplitude from these one-loop diagrams is what we refer here to as the (partial) LI corrections. Their computation is included in the HELAC-Onia code [54,55] and is thus available to everybody. In fact, another gauge-invariant O(α 6 S ) part, namely from pp → J/ψ +J/ψ +cc, is known [37]. It turns out to be small and can safely be ignored for our purposes. However, we wish to point out that the process pp → J/ψ + J/ψ + cc has its own interest as it can be a potential probe of the TPS at the LHC.
Let us add that we do not expect any specific kinematical enhancement of other NNLO topologies, in particular that of the double-real-gluon emission in view of the results of pp → J/ψ + Υ [42]. This is partly explained by the vanishing of these contributions when one gluon becomes soft, precisely where one can minimise the off-shellness of the other particles involved in the scattering and thus where these contributions could have been the largest. The cross-sections differential in the absolute rapidity difference between the J/ψ pair |∆y(J/ψ, J/ψ)| are shown in the left panels of Figures 1, 2 and 3 and are compared to the CMS, ATLAS and LHCb data. The NLO CS calculations are displayed by the red hatched bands in the figures. The partial LI contributions are represented by the green bands. As expected, the (partial) LI is significant at large |∆y(J/ψ, J/ψ)| region but negligible at small and intermediate |∆y|. An order of magnitude enhancement to the CS cross section is expected when |∆y(J/ψ, J/ψ)| ≥ 3.0. Nonetheless, despite the very large theoretical uncertainties from the scale variations, a discrepancy between the CS SPS and the experimental data is clearly visible at large |∆y(J/ψ, J/ψ)|, that is exactly where the DPS is expected to be important.
The invariant mass of the meson pair is also closely related to the rapidity gap |∆y(J/ψ, J/ψ)| (see the discussion in Ref. We have collected additional data-theory-comparison plots between the SPS CS yields and the LHC measurements for other observables in the appendix A.1. The data are compatible with the CS theoretical predictions but the LI contributions are found to be negligible for all the other distributions.

Comprehensive assessment of the colour-octet contributions
The whole LO CO contributions to di-ψ hadroproduction at the LHC up to O(v 7 ) in NRQCD have recently been computed by He and Kniehl [40]. Their study however bears on a single CO LDME set from an out-of-date LO single J/ψ hadroproduction fit [60] which was made with the early Tevatron data. Yet, their calculation seems to indicate that the CO contributions might be relevant at large |∆y(J/ψ, J/ψ)| and large M (J/ψ + J/ψ) due to similar t-channel gluon exchange diagrams than for the CS LI contributions. The aforementioned remaining discrepancy between this full SPS LO NRQCD calculation and the CMS data at large |∆y(J/ψ, J/ψ)| was then attributed to unknown missing higherorder QCD corrections to the CO contributions.
We however note that we do not anticipate any such so-called "giant" K factors in this region. Currently, no complete NLO CO calculation exists. Since it is important to deal with a complete set of CO channels in order to guarantee the large cancellation between Swave and P -wave contributions involved the NLO LDME fits of hadroproduction data, we consider that to rely on a LO -but complete-perturbative calculation and then to estimate the size of the missing higher-order corrections via the scale uncertainty is probably the most reasonable procedure to adopt.
An alternative approach to investigate the presence of possible "giant" K factors from new fragmentation topologies -if some are indeed relevant-without performing a full computation is that recently proposed by one of us in Ref. [61]. It has been proved useful for the single J/ψ hadroproduction case. The method is in principle general and applicable for the double J/ψ hadroproduction as well, although a new infrared divergence in double P -wave channels emerges [62]. We leave it for future studies since it may not apply to the whole phase space which we wish to consider here. Finally, we note that a similar enhancement from t-channel gluon exchange was expected for di-χ c production but its feed down was also found to be insignificant in the di-ψ yield [49].

Status and issues with the colour-octet transitions
Although the possibility for CO transitions is a robust prediction from NRQCD, their actual impact in the phenomenology has been the subject of debates for decades. The most glaring observations for the necessity of their presence are twofold. First, CO provide a natural solution for the infra-red divergence issue in P -wave production. Second, the LO v 2 NRQCD calculation involving only CS transitions still underestimates -even after including NLO QCD corrections-the yields of single J/ψ and ψ(2S) hadroproduction at large P T at the Tevatron and the LHC.
However, NRQCD computations even including CO contributions are unable to coherently describe -i.e. with the same CO LDMEs-the world data for pp, ep, γp, γγ and e + e − collisions. For a recent review, we guide the reader to Ref. [18].
The CO LDMEs are predicted to be universal non-perturbative objects by NRQCD, which should yield predictions compatible with all the data. The current status of their extractions is very confusing as their numerical values and their uncertainties are very disparate. The results of the fits of different groups disagree with each others. As long as the situation is not clarified, we believe that it is necessary to comprehensively consider these analyses instead of drawing conclusions based on a single CO LDME set as it is often done in the analysis of associated production of quarkonium (see Ref. [18] for some examples).
As such, we will use different LDME sets of which we briefly review the status and the possible limitations. As we said above, the available CO LDMEs for prompt J/ψ production are extracted from fits. According to the QCD accuracy of the short-distance coefficients (SDCs), we will categorise them in the 4 groups shown in Table 3 3. one fit based on a low-P T leading-logarithm (LL) resummed SDC [69], 4. one fit using a SDC using leading-power (LP) fragmentation matched to NLO SDC [70].
All of them have shortcomings and/or limitations. We enumerate them below: 1. First of all, we wish to emphasise that the LO fits are out-of-date and should be viewed as a pure tuning of the normalisation of the single J/ψ data. Since all of the LO fits are mainly performed with the help of intermediate and large P T hadroproduction data, where the "giant" K factors from NLO QCD corrections emerge, it is very hard to imagine that these values will give correct predictions for independent observables, like the double-J/ψ hadroproduction in our case, for which K factors would be different. We will therefore use them here for a pure illustrative purpose.
2. The LL fit in Ref. [69] concentrates on the P T (J/ψ) < m c region. The authors performed a small-P T resummation but without considering the contribution from the CS channel which is however known to saturate the data in this region [71,72]. The values of these LDMEs have never been used for the single J/ψ production at intermediate and large P T regions. They are included in our discussion like the LO fits in order to be exhaustive.
3. The NLO fit in Ref. [65] used the world data before 2011 without subtracting the feeddown contributions. The fit seems to yield a good agreement with the P T (J/ψ) < 30 GeV J/ψ yields data at different colliders but for γγ and e + e − collisions. However, it overshoots the P T > 30 GeV yields and fails to reproduce the polarisations of J/ψ, the energy-fraction distribution of the J/ψ in jets [73] and the yields of η c (by using heavy-quark spin symmetry). In addition, the SPS P T -differential cross section of J/ψ + γ [74] turns out to be negative at NLO with these values of CO LDMEs.
4. The NLO fit by Gong et al. [66] focus on the P T (J/ψ) > 7 GeV data. The feed-down contributions are subtracted. This LDME set is however not compatible with the yields (e.g. pp, γp and e + e − ) when P T (J/ψ) < 7 GeV, the polarisation of forward J/ψ [75] and the η c production. In addition, it yields to -unphysical-negative cross sections in pp → J/ψ + γ. In principle, this set is only applicable to J/ψ production with P T (J/ψ) > 7 GeV, i.e. only to the ATLAS fiducial region for our forthcoming discussion of double J/ψ production.
5. The two sets denoted sets 7 and 8 in Table 3 are two extreme cases of the PKU fit [67,68] after including the constraints from LHC η c data [67,68]. They supersede the fits including the P T (J/ψ) > 7 GeV hadroproduction data described in Refs. [76,77]. These LDME sets cannot reproduce the CDF polarisation measurement [78] -like all the other sets in fact-and are not applicable to P T (J/ψ) < 7 GeV. Both sets should only be used to di-ψ production in the ATLAS fiducial region.
6. The NLO+LP fit of Ref. [70] -as well as its update [79]-has been presented by its authors as the only fit able to reproduce the J/ψ data (both yields and polarisations) above 10 GeV after including the LP fragmentation contributions on top of the NLO calculations. However, it does not yield the correct η c cross section in the same P T region under the heavy-quark spin symmetry. As what concerns predictions for double J/ψ production, it is marginally applicable only in the ATLAS fiducial region with P T (J/ψ) > 8.5 GeV instead of 10 GeV.
Since we aim at a comprehensive analysis, we have considered all of the 9 sets listed in Table 3 to show how strongly the CO contributions depend on the available CO LDMEs. We should however recall during the discussion what we believe to be the region of applicability in P T (J/ψ) for the NLO(+LP) fits.

Numerical results
In this section, we will present our numerical results with the LO CO channels summed to the pure NLO CS channel 3 S 1 ) vary from set to set in Table 3, we will fix these values for the NLO CS channel 3 S 1 to those used in section 3.2. The uncertainty from these LDMEs is systematically subdominant compared to the scale uncertainty. All the feed-down contributions are properly taken into account as well.

LHCb data at √ s = 13 TeV
We start our discussion with the LHCb acceptance [23], where the P T (J/ψ) can be as low as zero. We have compared the CSM NLO +COM LO SPS (the green bands) with the data in Figure 4 for the ∆y(J/ψ, J/ψ) distribution and in Figure 5 for the invariant mass of the pair M (J/ψ + J/ψ) distribution. Like we have found for the CS LI contributions, the CO contributions are not relevant in the invariant mass distribution of LHCb. They start to be slightly visible in the tail of the ∆y(J/ψ, J/ψ) distribution.
This observation however very much depends on the set of CO LDMEs used. In particular, the only plausible set, i.e. set 5, in the small P T (J/ψ) region does not yield any significant contribution to the cross section. It also seems clear that none of the sets can fully account for the discrepancy between SPS and LHCb data in the last bins of dσ d∆y . Additional plots for the comparisons between CS NLO +CO LO SPS and data can be found in appendix A.2. The impact of the CO contributions on these additional distributions is in general minor.

CMS data at
The events analysed by CMS have larger P T (J/ψ) above 4.5 GeV to 6.5 GeV depending on the rapidity. In this region, the only applicable NLO fit is still set 5 taken from Ref. [65]. As opposed to the conclusion made in Ref. [40], the CO SPS contribution is either much suppressed compared to the CS SPS contributions or much smaller than the experimental data as shown in Figure 6 and Figure 7. Given that the LO fits (like that used in Ref. [40] (i.e. set 2)) are not plausible any more and that the only applicable fit is the NLO fit given by set 5, we draw the conclusion that our extraction of DPS in Ref.
[37] -made by neglecting the CO contributions-is still sound, which actually has been shown to be consistent with the ATLAS measurement thanks to a completely different method to disentangle the DPS from the SPS contributions and in a different kinematical region.

ATLAS data at
The transverse momentum cut on single J/ψ is largest in the ATLAS data sample with selected events satisfying P T (J/ψ) > 8.5 GeV. This leaves the LDME sets 5-8 as possible good

Conclusions
We have examined two SPS production mechanisms for di-J/ψ production at the LHC, which can be relevant in the control region used to determine the DPS. These are the We have also extensively compared our new SPS calculations with the existing LHC data. Our study indeed shows that the LI corrections can enhance the NLO SPS cross section at large |∆y(J/ψ, J/ψ)| and large invariant mass M (J/ψ + J/ψ). However, they are not sufficient to explain the discrepancy between SPS theoretical results and the LHC data in these regions. The inclusion of the DPS in the predictions is still crucial to account for the measurements.
On the other hand, the relevance of the CO contributions in the SPS yield strongly depends on the considered LDME set, thus with a very low predictive power -given the current status of understanding of the COM. It is in any case confined to the large |∆y(J/ψ, J/ψ)| region. We anyhow conclude that the CO contributions can only be important when compared to the ATLAS data but that the ATLAS DPS extraction via a 2D data-driven fit is very likely free of any bias due to a possibly underestimated CO contribution in their control region. Such a conclusion is backed up by studies [80,81] made in the colourevaporation model which offers a complementary framework to study the impact of CO transitions.

A Additional plots: further comparisons with data
This appendix gathers additional plots of comparisons between our SPS calculations and experimental data collected by the ATLAS, CMS and LHCb experiments.

A.1 Further comparisons with theory including partial CS LI corrections
We compare below our SPS CS NLO +LI calculation to the experimental data for other observables than the rapidity difference and the invariant mass. The transverse-momentum distributions of the pair P T (J/ψ + J/ψ) are shown in Figure 10 (CMS), in the left panel of Figure 11 (ATLAS) and in the top-right panel of Figure 12 (LHCb). The NLO +LI SPS green bands almost overlay at the red bands (NLO SPS), which implies that these LI corrections are not important for these distributions. It is interesting to note that the initial k T -smearing effect is important in the low P T (J/ψ + J/ψ) region, which illustrates that this distribution is indeed ideal to extract the transverse-momentum dependent information from the colliding partons inside the protons

A.2 Further comparisons with theory including CO contributions
Further comparisons between CS NLO +CO LO SPS results and LHCb data are shown in Figure 13 for P T (J/ψ+J/ψ), Figure 14 for ∆φ(J/ψ, J/ψ), Figure 15 for P T (J/ψ), Figure 16 for y(J/ψ), Figure 17 for y(J/ψ + J/ψ) and Figure 18 for A T (J/ψ, J/ψ) respectively. The inclusion of CO channels only slightly changes the corresponding predicted distributions of the SPS yield regardless of the set of LDMEs. Similar conclusions can be drawn for P T (J/ψ + J/ψ) and ∆φ(J/ψ, J/ψ) distributions in the CMS and ATLAS acceptances, which is clearly seen in Figures 19, 20 and 21.