Complete NLO QCD study of single- and double-quarkonium hadroproduction in the colour-evaporation model at the Tevatron and the LHC

We study the Single-Parton-Scattering (SPS) production of double quarkonia (J/psi+J/psi, J/psi+Upsilon, and Upsilon+Upsilon) in pp and pp(bar) collisions at the LHC and the Tevatron as measured by the CMS, ATLAS, LHCb, and D0 experiments in the Colour-Evaporation Model (CEM), based on the quark-hadron-duality, including Next-to-Leading Order (NLO) QCD corrections up to alpha_s^5. To do so, we also perform the first true NLO --up to alpha_s^4-- study of the p_T-differential cross section for single-quarkonium production. This allows us to fix the non-perturbative CEM parameters at NLO accuracy in the region where quarkonium-pair data are measured. Our results show that the CEM at NLO in general significantly undershoots these experimental data and, in view of the other existing SPS studies, confirm the need for Double Parton Scattering (DPS) to account for the data. Our NLO study of single-quarkonium production at mid and large p_T also confirms the difficulty of the approach to account for the measured p_T spectra; this is reminiscent of the impossibility to fit single-quarkonium data with the sole 3S18 NRQCD contribution from gluon fragmentation. We stress that the discrepancy occurs in a kinematical region where the new features of the improved CEM are not relevant.

In the recent years, there has been an accumulation of experimental hints [2,[23][24][25][26][27][28][29][30] that quarkonium pairs can be produced in a significant amount by two simultaneous parton-parton scatterings -the DPS. This is particularly true at large rapidity separations, ∆y ψψ , where the a priori leading Single Parton Scatterings (SPS) are suppressed since they generate highly correlated quarkonium pairs, thus at low ∆y ψψ . The region of large ∆y ψψ is therefore a good candidate for a control region for DPS extraction.
At small ∆y ψψ , all the experimental data sets are in fact in good agreement with the SPS predictions from the Colour-Singlet Model (CSM), i.e. the LO in the heavyquark relative velocity, v, expansion of Non Relativistic QCD (NRQCD) [52]. These predictions are known up to NLO accuracy [41,42,46]. In addition, the NRQCD Colour-Octet (CO) contributions are found to be negligible 1 in this region [2,43,51] which is in line with the expected suppression by O(v 8 ) with respect to the CS contributions.
For increasing ∆y ψψ , the lack of complete NLO NRQCD studies is prejudicial and opens the door to some debates [2,43,51] about the possibility for unexpectedly large SPS contributions from CO contributions. Indeed, owing to the large uncertainties in the LDME determinations [51], NRQCD at LO shows a very low predictive power, e.g. in regions where the DPS is thought to be the dominant source of quarkonium pairs. Hopefully, possible future NLO stud-ies could close this debate. This is important since recent direct and indirect DPS extractions based on quarkonia in pairs [2,[26][27][28] or in association with a vector boson [53,54] seem to point [54] at an unexpected flavour or momentum dependence of the parton correlations in the proton -as encoded in the well known quantity σ eff. [11]-, when compared to other (direct and indirect) extractions [55][56][57][58][59][60][61][62][63][64][65][66].
Using the colour-evaporation model (CEM) [67,68] a model based on quark-hadron duality but which shares some features of NRQCD [69], in particular the direct production of vector quarkonia from gluon fragmentation-we wish to advance our understanding of the SPS contributions to quarkonium pairs. In addition, its implementation is very similar to that of open heavy-flavour production and can be done in MADGRAPH5 AMC@NLO [70] with some tunings. Finally, it is straightforward to treat the feed-down contributions (e.g. from χ c ) to prompt J/ψ in the CEM. Altogether, this allows us to perform the first complete NLO study of quarkonium-pair production using one of the widely used quarkonium-production models.
In the case of J/ψ + Z [53] and J/ψ + W [54] production, we have shown that the CEM provides an upper limit on the SPS contributions. This is also likely the case for J/ψ + Υ production and for the J/ψ + J/ψ case in the kinematical region where gluon fragmentation to both quarkonia is expected to be dominant. More generally, it offers an indirect way to scrutinise whether some specific configurations of the heavy-quark pair receive at NLO kinematicallyenhanced contributions, which result in large K factors (see e.g. [41]). Indeed, if we were to observe a large K factor to the di-quarkonium CEM yields, where all the pair (spin and colour) configurations are summed with the same weights, this would necessarily signal a potential large K factor to some NRQCD contributions.
Such a complete NLO study for di-quarkonium necessitates a coherent determination of the non-perturbative CEM parameters -one per particle. Therefore, an interesting side product of our study is the corresponding NLO study of the p T -differential cross section of single J/ψ, ψ(2S ) and Υ(nS ) hadroproduction. To what concerns the ψ(2S ) and Υ(nS ), this is the very first study of this kind. So far the NLO CEM studies [71,72] were held for the p T -integrated yield at α 3 s . This paper is organised as follows. In section 2, we explain the methodology of our NLO CEM calculation. In section 3, we discuss our original results for the p T distribution of single J/ψ, ψ(2S ) and Υ(nS ) in the CEM at NLO, which we use to fit the corresponding non-perturbative CEM parameters. In section 4, we then present our results for the production of di-J/ψ for the CMS, ATLAS, and LHCb acceptances, of J/ψ + Υ in the D0 acceptance and for the di-Υ in the CMS acceptance. Section 5 is devoted to our conclusions and outlook.

The CEM in a few formulae
In the CEM, a given quarkonium-production cross section is obtained from that to produce the corresponding heavy quark-antiquark pair QQ with the sole constraint that its invariant mass lies between twice the quark mass, 2m Q , and twice that of the lightest open-heavy-flavour hadron, 2m H . The same logic applies in the case of a pair of quarkonia. The cross section for single quarkonium production is then given by and that for the production of a pair, Q 1 + Q 2 , of quarkonia where the corresponding (doubly) differential cross section for QQ (Q 1Q1 +Q 2Q2 ) production as a function of the pair invariant mass(es), m QQ (m Q 1Q1 m Q 2Q2 ) and P Q i is a non-perturbative parameter encapsulating the probability for the hadronisation of the QQ pair into the quarkonium Q i . It is supposed to be universal and independent of the production of the pair.
In principle, having the heavy-quark cross section differential in the invariant mass, dσ/dm QQ is sufficient to obtain the short-distance part of the CEM for single or associated production and correspondingly for pair production. The automated tool MADGRAPH5 AMC@NLO with specific cuts can provide such cross sections up to NLO accuracy, also differential in other variables, like the rapidity or the transverse momentum of the QQ pair which translates 2 into that of the quarkonium Q. Such cross sections should just then be multiplied by the non-perturbative parameter P Q i which is usally tuned to match the single-quarkonium production data.
3. The p T -differential cross section for singlequarkonium hadroproduction at NLO The existing CEM studies of quarkonium production at RHIC, the Tevatron and the LHC rely on a hard-scattering matrix element at one loop for inclusive heavy-quark production, namely α 3 s (see [20] for a recent review). This is based on the well-known multi-differential MNR computation [74] using the aforementioned invariant-mass cut. At this order, a heavy-quark pair with a non-zero p T (irrespectively of the invariant mass of the pair) comes from realemission graphs, where a final light parton recoils against the QQ pair. The virtual-emission contributions do not contribute away from p T,QQ = 0. Such existing computations for p T,QQ 0 are effectively Born-order or tree-level computation from the partonic processes gg[qq] → (QQ)g or gq → (QQ)q, and thus not effectively at NLO accuracy. As a case in point, the renormalisation-scale dependence of the resulting cross section is simply that of the third power of α s (µ R ).
Thanks to MADGRAPH5 AMC@NLO, we are able to provide complete NLO CEM hadroproduction results for dσ/d p T,Q by computing pp → (QQ) CEM + 1 parton up to α 4 s where the subscript indicates that the pair invariant mass is integrated as in Eq. (1). A first J/ψ study was presented along with our J/ψ + Z CEM computation [53]. Here we go further and consider in addition the ψ(2S ) and Υ(nS ) cases. We also discuss in more details the resulting CEM parameter depending on whether it is fit at mid or large p T or on the p T -integrated yields.  As what regards the parameters of our computation, they remain very standard. We have used the PDF set NLO NNPDF 3.0 set [79] with α s (M Z ) = 0.118 provided by LHAPDF [80] from which we have derived the PDF uncertainty. The latter remains negligible compared to the factorisation-and renormalisation-scale uncertainties, which are evaluated by varying them independently in the interval 1 2 T . Like in [53], we use m c = 1.27 GeV for charmonium production in the CEM as suggested in [72]. It is important to note that the quark mass enters the cross section both via dσ/dm QQ and via the integration range. Results with m c = 1.5 GeV are sometimes slightly different. However, when the CEM is tuned to data, the mass dependence is mostly absorbed in the change of P Q i and the physics conclusion always remains nearly identical. For the bottomonia, we have used m b = 4.7 GeV. For the upper bounds of integrations, 2m H , we have used 3.728 GeV for cc and 10.56 GeV for bb.   Figure 1: The p T (and y) differential cross sections of single (a) J/ψ and ψ(2S) and (b) Υ(nS ) (n = 1, 2, 3) production in the CEM. The plotted data from ATLAS (8 TeV) [75] and from CMS (7 TeV) [76] were used to fit We have performed a number of fits of P Q using the experimental data of single inclusive prompt J/ψ and ψ(2S) and Υ(nS ) data from ALICE, ATLAS and CMS over different p T ranges. We could also have used the very precise LHCb data [81][82][83] but we preferred to restrict our fit to central rapidity data. Including them would not have changed our conclusions since the CEM does not describe well the p T spectrum in any case. In the J/ψ and Υ(nS ) cases, the obtained values of P Q i correspond to prompt production. For ψ(2S ), they hold for direct production. Table 1 gathers the used kinematical ranges and the corresponding fit results at LO and NLO.
Since the K factors for pp → cc (+jet) + X and pp → bb (+jet) + X near threshold are larger than unity, the P Q i at NLO are correspondingly smaller than at LO. Moreover, since the CEM p T spectra are usually too hard (see [20]), the P Q i also tend to decrease in order to match data at high p T and the fits overall worsen. This well-known (LO) trend is indeed confirmed at NLO. In the present LHC kinematics, this is particularly obvious in Fig. 1a and Fig. 1b.
The ALICE J/ψ data set and one CMS Υ data set extend to p T = 0 which allowed us to fit the p T -integrated cross section with a NLO α 3 s computation of pp → (QQ) CEM + X. These results for the Υ case for the entire p T range are comparable to those of Vogt, i.e. 2 to 3% depending on m b , the scale choice and the PDFs (see "Υ1" and "Υ4" of Table 8 of Ref. [71]). We see that the CEM parameters obtained by fitting the p T spectra are systematically smaller than those obtained by fitting p T -integrated yields. For the J/ψ case, we have quoted a range. Indeed, as can be seen in [72], σ NLO (cc) shows a very large uncertainty, which even tends to increase at large √ s ending up to be as large as one order of magnitude. The obtained lower values are systematically much smaller than the open-charm data. If we were to fit the ALICE J/ψ data with σ(cc CEM ) computed with the scale values corresponding to these lower values, the discrepancy would be absorbed in P ψ which would become anomalously large, even above unity in some cases. This would be unphysical. Since the open-charm data systematically lie between the central and upper NLO values, we thus only quote in Table 1 the corresponding range for P ψ . It is in fact in line with the values quoted in [72] for m c = 1.27 GeV but for different (fixed) scales and PDFs. The data sets used for these older fits are also obviously different.

LO and NLO CEM results for di-quarkonium hadroproduction
In this section, we present our LO and NLO CEM results for all the existing LHC and Tevatron results [24,25,[28][29][30]84], but the D0 analysis [24] for which no normalised distributions were released and the early LHCb analysis at √ s = 7 TeV [23] which we consider to be superseded by their 13 TeV analysis. The corresponding kinematical conditions are summarised in Table 2.
Like for the NLO single-quarkonium study presented above, we employ the NLO NNPDF 3.0 set [79]. The dependence of the result on the renormalisation µ R and factorisation µ F scales is quantified by varying them independently in the interval 1 2 µ 0 ≤ µ R , µ F ≤ 2µ 0 where µ 0 depends on the considered system. For charmonium and bottomonium pairs, it is fixed to be (4m Q ) 2 + p 2 T where p T is randomly selected from one of the pair members. For charmonium+bottomonium, it is the average of the transverse masses, 0.5 × (m T 1 + m T 2 ). We also do not consider uncertainties from the heavy-quark mass as they are mostly absorbed in the CEM parameters, P Q i . This is surely the case for the invariant-mass-integration region. The remaining uncertainty from the value of the hard matrix element may differ, but in view of the data-theory disagreements which we discuss next, we consider this approximation to be reasonable.

J/ψ pairs
Let us first discuss our results of the CEM calculation of J/ψ-pair production in the CMS setup. The differential cross section in the rapidity separation, |∆y ψψ |, is shown in Fig. 2a, in the invariant mass, M ψψ , in Fig. 2b, and in the transverse momentum of the J/ψ-pair, p T ψψ , in Fig. 2c. We see that the CEM results are at least an order of magnitude below the experimental data of CMS at both LO and NLO, even considering their (large) uncertainties. We note that the scale uncertainty in the NLO calculations is half of that in the LO ones which indicates the absence of kinematically-enhanced topologies. The regions of large |∆y ψψ | and large M ψψ are those where the DPS contributions are extracted and our computations confirm that SPS contributions, from the CEM for sure but likely as from other models, are negligible there. Unsurprisingly, the NLO yield populates the p T ψψ distributions but its contribution is clearly too small and does not even show the bump generated by the kinematical cut in the CMS acceptance. Such a bump is well seen in NLO CSM computations, which describes the data well all over the entire spectrum [2,51]. The very same observations can be made for the p T ψψ distributions measured by ATLAS [28].
We note that ATLAS only released the M ψψ and |∆y ψψ | distributions for their fiducial yields. For the latter, only a small fraction of the events passes the muon cuts, even after the J/ψ cuts (which are also stringent, p T ψ > 8.5 GeV). In addition, the CEM yield is computed over a tiny fraction of the possible cc invariant masses. Overall, the relevant multi-dimensional hyperspace where the integration is performed can be extremely small and complex. The result depends on an extremely small part of the physical phase for pp → cc + cc + X, especially at LHC energies. To help the integrator find the CEM 'domain', one has to enlarge the invariant-mass regions at the MC generation level and then to restrict it at the histogramming level. This is unfortunately extremely ineffective. Whereas the MAD-GRAPH5 AMC@NLO integrator manages to perform well the integration at LO in a reasonable amount of time, it be-comes highly CPU consuming at NLO, for instance O(10 8 ) CPU·seconds to get distributions where all the bins are simply populated. Unfortunately, we did not manage to obtain reliable NLO results for these distributions. As such, we only plot the LO results (see Fig. 2e for |∆y ψψ |, Fig. 2f for M ψψ ). We expect the NLO results to be similar in view of the LO/NLO ratio for the corresponding CMS distribution in a similar (inclusive) phase space.
Contrary to CMS, ATLAS also released their data as a function of |∆φ|, the azimuthal angle between both J/ψ. It is useful in quantifying the relative size of the DPS vs SPS contributions, especially when transverse-momentumsmearing effects can be neglected. Indeed, in such a case, the SPS contributions usually exhibit a peak at |∆φ| = π (both particles recoil on each other) and sometimes at |∆φ| = 0 (the pair recoils against a third particle) whereas the DPS contributions should exhibit a flat distribution if both partonic scatterings are indeed uncorrelated. This remains of course an approximation although, until now, never falsified.
Along these lines, the concave data |∆φ| distribution shown in Fig. 2g is indicative of a significant SPS contributions. According to ATLAS [28], it amounts to about 90% of the entire yield. Clearly, the CEM is unable to account for this SPS contribution. The same distribution measured by LHCb at 13 TeV shown on Fig. 2h is however much more intricate to interpret. Indeed, the LHCb measurement was performed without p T ψ cut which allows the momentumsmearing effects to be significant. As a consequence, the DPS vs SPS separation is much more involved. On a logarithmic plot, the NLO CEM yield already looks nearly flat, not very different than the shape of the data distribution. Yet, the normalisation is again off by more than an order of magnitude. Let us stress that, for the LHCb acceptance, we have used a P Q value fit on the p T -integrated yields, which is the largest of all those discussed above. One observes a similar gap on the other distributions 3 (see Fig. 2i-2l) which confirms that the CEM is unable to account for any measured di-J/ψ data sets. This is even the case in regions where at the same time the DPS contributions are expected to be mild -if not negligible-and the CSM has been found to match the data.

Υ + J/ψ pairs
We now move to the J/ψ + Υ(nS ) case as measured by the D0 Collaboration at √ s = 1.96 TeV [26]. The only 3 A T , also called the transverse momentum asymmetry, is defined as released kinematical distribution was that of |∆φ| which we have compared to our CEM computations in Fig. 3a. We note that this measurement was performed at low p T and we have used the corresponding CEM parameter values. If we had used parameters fit to the p T -differential data, the CEM predictions would have been even smaller. At |∆φ| = π, the NLO CEM is at best 5 times below the data and ends up to be 100 times lower at |∆φ| = 0. This is in line with the current interpretation of these D0 data, namely that they are dominated by DPS contributions [27].

Υ pairs
Finally, we move to Υ(1S)-pair production as measured by the CMS experiment. In a first study at 8 TeV [84] for |y Υ | < 2.0, they only measured the integrated cross section, which they found to be Very recently, CMS performed a new study at 13 TeV [30] with significantly more events which allowed them to perform differential measurements as a function of ∆y and M ΥΥ . The comparisons are shown in Fig. 3b & 3c.
As of now, there does not exist any direct or indirect DPS/SPS separation. As such, we are not able to claim that the CEM is in contradiction with the data. Yet, any reasonable estimate of the DPS yield would indicate that the SPS fraction should be significant [20]. This would mean that the CEM indeed cannot account for the SPS yield to di-Υ in the CMS acceptance.

Conclusions
We have presented the first CEM study at LO and NLO for the SPS yields in hadroproduction of quarkonium pairs. Our computation -fully accounting for contributions up to α 5 s -was performed thanks to a tuned version of MAD-GRAPH5 AMC@NLO taking into account the specificities of the CEM. Except for those kinematical distributions where the LO distributions are trivially suppressed, the K-factors we have found are systematically on the order of the unity, in particular at large ∆y QQ and M QQ . This lends support to the irrelevance of possible kinematically-enhanced contributions from QCD corrections in these regions (see also [38,47,50,51]) where the dominance of DPS contributions is sometimes questioned. Owing to the similarities between the CEM and the COM of NRQCD, we foresee a similar situation when the first NLO COM studies are performed.
On the quantitative level, we have compared our computations to a large selection of data for J/ψ + J/ψ, J/ψ + Υ, and Υ + Υ hadroproduction in pp and pp collisions at the LHC and the Tevatron as measured by the ATLAS, CMS, D0, and LHCb collaborations. In all the cases, the computed yields are one to two orders of magnitude below the experimental data. This is also the case in the kinematical regions where it is established that the SPS contributions are dominant, or equivalently that the DPS contributions cannot reasonably describe the data. As such, it provides evidence that the CEM does not encapsulate the leading production mechanism for this SPS yield. This is another case, with J/ψ + cc [20], where the CEM fails to describe the data while the CSM can.
In order to present coherent results at one-loop accuracy, we have also studied p T -differential cross sections for single quarkonium hadroproduction up to α 4 s . We have used these in order to fit the non-perturbative CEM parameters. As far as the description of p T -differential yields are concerned, our results therefore naturally supersede existing CEM results available in the literature (see e.g. [71]) which were only performed up to α 3 s . We have also updated our J/ψ NLO results made along a previous J/ψ + Z NLO CEM study [53]. Let us add that this is the first time that p Tdifferential cross sections for Υ(nS ) are computed at this order and compared to the data.
Overall, the CEM features for single quarkonium hadroproduction observed at α 3 s , i.e. at LO, are confirmed and the CEM remains unable to provide a satisfactory description of the single-quarkonium-hadroproduction data with too hard a spectrum at large p T . For this reason, we provide different values of the CEM parameters as needed makeshift if one wants to still perform other phenomenological studies similar to the present one for di-quarkonium production. Indeed, the CEM still represents a handy approach to estimate quarkonium cross sections when computations under other approaches like the CSM and NRQCD are too complex, especially beyond LO accuracy.