J/psi-Pair Production at Large Momenta: Indications for Double-Parton Scatterings and Large alpha_s^5 Contributions

We demonstrate that the recent studies of J/psi-pair production by CMS at the LHC and by D0 at the Tevatron reveal the presence of different production mechanisms in different kinematical regions. We find out that next-to-leading-order single parton scattering contributions at alpha_s^5 dominate the yield at large transverse momenta of the pair. Our analysis further emphasises the importance of double parton scatterings --which are expected to dominate the yield at large J/psi-rapidity differences-- at large invariant masses of the pair in the CMS acceptance, and thereby solve a large discrepancy between the theory and the CMS data. In addition, we provide the first exact --gauge-invariant and infrared-safe-- evaluation of a class of leading-P_T (P_T^-4) next-to-next-to-leading-order contributions at alpha_s^6, which can be relevant in the region of large values of P_T,min=min(P_T1,P_T2). Finally, we derive simple relations for the feed-down fractions from the production of an excited charmonium state with a J/psi in the case of the dominance of the double parton scatterings, which significantly deviate from those for single parton scatterings. Such relations can be used to discriminate these extreme scenari, either DPS or SPS dominance.


Introduction
Heavy-quarkonium production has attracted considerable interest in the high-energy physics community since the J/ψ discovery, exactly forty years ago. It indeed probes the strong interaction at the interplay of its perturbative and non-perturbative regimes [1]. It can also help to understand a new dynamics of hadron collision where multiple (hard) parton scatterings (MPS) take place. MPS are normally very rare since already a single (hard) parton scattering (SPS) is rare as compared to soft scatterings. Owing to the high parton flux at high energies, MPS should be likelier at the LHC, starting with two short-distance interactions from a single hadron-hadron collision, usually referred to as double parton scattering (DPS). These have been searched in 4-jets [2,3,4], γ + 3-jets [5,6], W + 2-jets [7,8], J/ψ + W [9], J/ψ + Z [10], 4-charm [11], J/ψ+charm [11] and J/ψ + J/ψ [12] final states.
Along these lines, J/ψ-pair hadroproduction is of great interest. First, it provides an original tool to study quarkonium production in conventional SPSs. Most of the earlier theoretical studies are based on SPSs [13,14,15,16,17,18,19,20,21,22]; some using the colour-singlet model (CSM) [23], others Nonrelativistic QCD (NRQCD) [24]. Moreover, it is widely claimed that DPSs [25,26,27,28] could indeed be a significant source of J/ψ pairs at the LHC in proton-proton collisions and in proton-nucleus/nucleusnucleus collisions [29,30]. Generally, it remains a poorly understood process. Its measurement with both J/ψ decaying into a muon pair is a clean signal, accessible to most experiments, which is complementary to the DPS studies based on open charm mesons and hadronic jets. With re-spect to the latter, it allows one to investigate the physics of DPS at lower scales and lower x where different mechanisms may be at work (see e.g. [31]).
The first observation of J/ψ-pair events dates back to that of the CERN-NA3 Collaboration [32,33]. Recently, the LHCb [34], CMS [35] and D0 [12] collaborations reported their measurements at the LHC and the Tevatron. In contrast to Kom et al. [25], we recently pointed out [21] that no definite conclusion on the presence of DPSs in LHCb data [34] should be drawn given the very large theoretical uncertainties on the SPS predictions. However, the recent D0 [12] study could provide the very first separation of the DPSs from SPSs and a measurement of σ DPS ψψ and σ SPS ψψ by using the yield dependence on the (pseudo)rapidity difference between the J/ψ pair, as it was first proposed in [25]. Although such a separation relies on a good modelling of the DPS and SPS rapidity-difference spectra, this can reasonably be considered as the first observation of a DPS signal in quarkonium-pair production, even if SPSs also contribute in a significant fraction of the D0 acceptance. Two fundamental remaining questions are whether such DPS contributions are also of importance elsewhere than at large rapidity difference, ∆y, and whether they agree with theory. In addition, the recent CMS analysis, up to large J/ψpair transverse momenta (P ψψ T ) brought to light a new striking puzzle. As pointed out in [22], the P ψψ T and invariant mass, M ψψ , spectra measured by CMS [35] severely overshoot the SPS contributions -even at next-to-leading order (NLO), i.e. α 5 s . In this Letter, we first show that the SPS yield extracted by D0 can only be reproduced thanks to the additional α 5 or feed-down contributions from J/ψ + ψ . Then, along the lines of [25], we model the DPS spectra based on 3 parametrisations of existing single J/ψ data and extractaccounting for the predicted SPS yield up to α 5 s -from a fit to the CMS results [35] the effective cross section σ eff which characterises the effective spatial area of the partonparton interactions. Our fit result is then found to be well compatible with the DPS D0 results [12] which means that we de facto provide a solution to the aforementioned puzzle, with a coherent description of CMS and D0 results.
In addition, we provide an original test of the DPS vs SPS dominance based on the yields involving excited states. Such a test can be used to validate our explanation of the CMS puzzle. Finally, we evaluate the first piece of the next-to-next-to-leading-order contributions from gg → J/ψJ/ψcc (denoted ψccψ) which is gauge invariant and infrared finite. Although it was expected to be enhanced at large P ψ T min = min(P ψ T 1 , P ψ T 2 ), we find it is dominant only when the yields are out of reach for current experiments and we conclude that an evaluation up to α 5 s accuracy is probably sufficient at present time.

SPS Contribution to J/ψ + J/ψ production
In this section, we outline the computation of the SPS contribution in the CSM [23] or equally NRQCD at LO in v 2 . The amplitude to produce of a pair of S -wave quarkonia denoted Q 1 and Q 2 , of given momenta P 1,2 and of polarisation λ 1,2 accompanied by other partons -inclusive production-, noted k, is then given by the product of (i) the amplitude to create the corresponding double heavy-quark pair, in the specific kinematical configuration where the relative momenta of these heavy quarks (p 1,2 ) in each pairs is zero, (ii) two spin projectors N(λ 1,2 |s 1,3 , s 2,4 ) and (iii) R 1,2 (0), the radial wave functions at the origin in the configuration space for both quarkonia. At the partonic level, the SPS amplitude thus reads: where one defines from the heavy-quark momenta, q 1,2,3,4 , P 1,2 = q 1,3 + q 2,4 , p 1,2 = (q 1,3 − q 2,4 )/2, and where s 1,3 ,s 2,4 are the heavy-quark spin components and δ c i c j / √ 3 is the colour projector. N(λ|s i , s j ) is the spin projector, which has a simple expression in the non-relativistic limit: where m Q is the heavy-quark mass and Γ S is ε λ µ γ µ when S = 1 (e.g. J/ψ, ψ ). Such a partonic amplitude is then squared, summed over the colour and spin of external partons and convoluted with the partonic densites (PDFs) in the allowed kinematical phase space.
In the case of gg → J/ψJ/ψ + cc, there exist more than 2000 graphs (see Fig. 1 (c-f)) which are non-trivially zero even after the topologies with a single gluon connected to an individual heavy-quark lines are removed. In order to generate the amplitude for this process, and the other ones considered in this study, we use HELAC-ONIA, described in [36], which generates the amplitude based on Eq. (1) using recursion methods [37]. It can also deal with P-wave production which involves the derivative of the amplitude in the relative momentum of the heavy-quark forming the quarkonia.
HELAC-ONIA also performs the helicity-amplitude calculations and the convolution with the PDF as well as the final-state variable integration. At LO, J/ψ-pair production at colliders is from gg → J/ψ + J/ψ at α 4 s (see e.g. Fig. 1a). At α 5 s , one needs to consider real-emission contribution (see e.g. Fig. 1b) as well as loop corrections. In [22], a full NLO computation showed that for 1 P ψ T > 2 GeV, it is sufficient to rely on a NLO [21] evaluation where only the LO and NLO real-emission contributions are accounted for; the latter being regulated by an infrared cut-off. This is easily explained by the P 2 T suppression of the loop contributions. For the ψccψ contribution which we compute, there is no need to apply any infrared cut-off. Since the LO kinematics is that of a 2 → 2 process (Fig. 1a), it generates a trivial P ψψ T dependence (δ(P ψψ T )) in the collinear factorisation with conventional PDFs. We have therefore accounted for the k T 's of the initial partons with a Gaussian distribution with k T = 2 GeV as in [21] in order to obtain fairer comparisons of the P ψψ T spectra. Technical details about the implementation can be found in [38].

DPS Contribution to J/ψ + J/ψ production
Quarkonium-pair (Q 1 +Q 2 ) production from DPSs is usually assumed to come from 2 independent SPSs which create each a single quarkonia. One therefore assumes where m = 1 for identical final-state particles and m = 2 otherwise, σ Q is the cross section for single Q production. σ eff is expected to account for the effective size of the parton-parton interaction and should thus be universal -that is process-independent-as well as √ s-independent if the factorisation holds as in Eq. (1). Yet, there does not exist proofs of such a factorisation. Factorisation-breaking effects have been discussed in a number of recent studies (see e.g. [31,40,41]). As of today, data is needed to test it case by case. Finally σ eff cannot be determined from first principles or from perturbative methods.
Anticipating the discussion of the D0 results in the next sections, which extracted a smaller σ eff from J/ψ-pair production than usually found in previous studies involving jet observables, we note that σ eff could very well depend on the flavour of the initial partons (see e.g. [31]). In the present case, these are gluons only, whereas for high-P T jets with W, light-quark initiated processes give a significant contribution. In addition, the process considered here occurs at a rather low momentum scale. Finally, we stress that if the 1v2 contributions 2 discussed in [42] matter, this would result in larger DPS contributions, thus in a smaller extracted σ eff .
Following the common practice, DPS yields are simply computed from the corresponding measured single-J/ψ yields using a Monte Carlo code with as input a parametrisation of σ ψ , see e.g. [25]. Let us stress here that if one uses a parametrisation of prompt, i.e. excluding b-decay J/ψ, or direct quarkonium data, one would predict DPS yields for prompt or direct quarkonium pairs. Given the level of understanding of quarkonium-production mechanisms, using theoretical models to compute σ ψ entering Eq. (1) would only inflate the theoretical uncertainties. As for now, the objective of DPS studies is to quantify their impact and to verify the factorisation hypothesis in a given kinematical domain. In the present study, the objective is for instance to address the apparent discrepancy between the predicted SPS yield and the CMS data.
We thus use the setup proposed in [25] for the single J/ψ cross section, σ ψ , (Eqs.(2-4) of [25] slightly improved since we used 3 σ ψ fits in order to assess the systematic uncertainties attached to the parametrisation of σ ψ . Details regarding these fits are given in the Appendix A. As an illustrative comparison, the fits to LHC and Tevatron data are shown on Fig. 2. We stress that the data used for the 3 fits are for prompt J/ψ. The corresponding short-distance matrix LHCb, √s=7 TeV, 2 < y ψ < 2.5 CDF (x0.5), √s=1.96 TeV, |y ψ | < 0.6 Crystal-Ball fit 1 Crystal-Ball fit 2 Crystal-Ball fit 3 Figure 2: Illustrative comparison between 3 fits of σ ψ and LHCb [43] and CDF [44] data for prompt J/ψ production. element has been added to a specific branch of HELAC-ONIA [36] with as inputs the fit parameters. This branch is separate from that used to compute the SPS contributions (for technical details, the reader is referred to [38]).
We stress that the purpose of using an event generator such as HELAC-ONIA for DPS computations is to perform the spin-entangled decay of the J/ψ's under different polarisation hypotheses to apply the fiducial cuts (on the muons) of a given analysis if the muon acceptance was not corrected, as for the D0 analysis [12]. A simple combination of 2σ ψ would not allow for this. In the D0 case, a variation of λ θ within −0.45 < λ θ < 0.45, which represents a reasonable envelope of the existing experimental measurements in similar conditions (see e.g. [45]), induces a systematical 20% uncertainty. This can be compared to the 25 % systematical uncertainty on the corrected muon acceptance quoted by CMS [35] due to the unknown J/ψ polarisation.
As in [25], we use the MSTW 2008 NLO PDF set [46] and the factorisation scale µ F = m ψ T = (m 2 ψ + (P ψ T ) 2 ) 1/2 . The fit uncertainties attached to our DPS evaluation are discussed in the next sections. Finally, let us mention that we have studied one factorisation-breaking effect which is the possible correlation between two partons from a single proton as encoded in the double PDF (dPDF) of [47]. We did not find any relevant difference in the region considered in this study.

Feed-down fractions under DPS dominance
If one sticks to a simplistic -although widely used-view of the DPS production mechanisms as the one presented above, it is possible to derive general relations between the feed-down fractions of the DPS yields for double and single J/ψ production. These can be used to evaluate the feeddown impact, but also, by returning the argument, to test a possible DPS-dominance hypothesis by directly measuring pair productions involving the excited states.
Just as we define the fractions, F direct ψ , F χ c ψ and F ψ ψ , of single J/ψ produced directly, from χ c decay or from ψ decay, one can define various feed-down fractions for J/ψ + J/ψ. However, one should keep in mind that it would probably be experimentally very challenging to measure (and subtract) 3 the yield of χ c + χ c or even χ c + ψ . We therefore limit ourselves to define F χ c ψψ (resp. F ψ ψψ ) as the fraction of J/ψ + J/ψ events from the feed-down of at least a χ c (or resp. a ψ ) decay. In other words, F χ c ψψ is the fraction of events including a prompt J/ψ (direct or from χ c and ψ feed-down) plus a J/ψ identified as from a χ c . Although it is probably very difficult to measure it, we also define F direct ψψ as being the pure direct component, excluding all the possible feed-downs, which can be easier to predict theoretically.
Assuming Eq. (1) holds for all charmonia, one gets 3 In order to obtain numbers, let us recall that the world data tell us that F direct ψ , F χ c ψ and F ψ ψ are close to 60%, 30% and 10%. We then obtain F χ c ψψ 50%, F ψ ψψ 20% and F direct ψψ 35%. Although F χ c ψψ and F ψ ψψ are experimentally accessible via σ((χ c → J/ψ) + J/ψ) and σ((ψ → J/ψ) + J/ψ), they are not sufficient to determine the pure direct yield since Its extraction would require the measurement of σ(χ c + ψ ).

Feed-down fractions under SPS dominance
On the contrary, one expects a larger feed-down from ψ if SPSs dominate. In the CSM or NRQCD at LO in v 2 , the hard part for ψ + J/ψ and J/ψ + J/ψ is identical; only |R(0)| 2 differ. Taking |R ψ (0)| 2 = 0.53 GeV 3 [49], whereas |R J/ψ (0)| 2 = 0.81 GeV 3 , and B(ψ → J/ψ) = 55% [50] as well as accounting for a factor 2 from the final-state symmetry, the ratio of F ψ ψψ /F direct ψψ -defined as in section 3.1-is expected to be as large as 0.53/0.81×0.55×2+(0.53/0.81× 0.55) 2 0.85. It may even be a bit larger since we neglected σ(χ c + ψ ) in this evaluation of F ψ ψψ . The latter approximation is justified, since we checked that neither σ(χ c + J/ψ) nor σ(χ c + ψ ) are significant under SPS dominance. In the CSM they are absent at α 4 s . The colour-octet (CO) contributions for the production of these pairs are small because of the small size of the CO non-perturbative parameters (also called LDMEs) [51] and the absence of any kinematical enhancement. In the remaining our this work, we will thus consider that σ prompt SPS In turn, we also have F ψ ψψ 0.85/(1 + 0.85) 46% at any order in α s .

DPS vs. SPS
To summarise, in the SPS case, F ψ ψψ can be as large as 46% whereas F χ c ψψ is expected to be small. In the DPS case, 3 The derivation of Eq. (2) for χ c follows from the decomposition of the different sources of a prompt J/ψ + a J/ψ from a χ c . Namely, one has : direct J/ψ + χ c , χ c + χ c and ψ + χ c . Their cross section with the relevant branchings can then be decomposed in terms of single quarkonium cross sections using Eq. (1) taking care of not double counting χ c + χ c (m = 1). Their sum divided by the cross section for a pair of prompt J/ψ decomposed likewise then reads as Eq. (2)  F ψ ψψ is half as small, around 20%, and F χ c ψψ large, around 50%. This clearly means that the relative measurements of charmonium-pair production of different states can serve as a clear test to pin down DPS or SPS dominance since they correspond to rather opposite predictions. This can reliably be done provided that the single-charmonium yields are known in the same kinematical region. We stress that, for such a test, we do not need to know the value of σ eff which does not appear in Eq. (2).

The early LHCb data at low transverse momenta
Let us first look at the LHCb data. We claimed in a recent work [21] that there was no compelling reason to call for significant DPS contributions in order to describe the J/ψpair measurement by LHCb at 7 TeV in the forward rapidity region (2 < y ψ < 4.5). In particular, there is absolutely no difficulty to reproduce the measured yield with the SPS contributions alone, see the first line of Table (1). In fact, the LO [21] and NLO [22] prompt SPS values even tends to be above the LHCb one, leaving room for a possible DPS yield only when the uncertainties are accounted for. We stress that this measurement was performed without any lower P T cuts and that, in this case, the LO and NLO SPS predictions are in very good agreement, showing a good convergence of the perturbative series. We will comment later on the DPS predictions.

The D0 data up to large ∆y and the observation of DPS
contributions We now discuss the recent D0 measurement [12]. Thanks to a wide (pseudo)rapidity coverage of about 4 units (|η| ≤ 2), the D0 detector made the first extraction of the DPS contributions to J/ψ-pair production possible. As neatly discussed in [25], the yield as a function of the rapidity difference ∆y between both J/ψ's should be a good observable to distinguish DPS and SPS events. The DPS events have a broader distribution in ∆y than the SPS ones. For the latter, large values of ∆y imply large momentum transfers, thus highly off-shell particles, and are strongly suppressed. It is not the case for DPSs where the rapidity of both J/ψ is independent. Large rapidity differences are only suppressed because the individual yields are suppressed for increasing rapidities.
By fitting the ∆y distribution of their data, D0 managed [12] to separate out the DPS and the SPS yields, i.e. σ DPS ψψ and σ SPS ψψ . They found that about half of the (prompt) yield was from SPSs, the remaining half from DPSs, about 60 fb (see Table (1)).
In [21], we discussed the relevance of taking into account P T -enhanced topologies (e.g. Fig. 1b) Table 1: Comparison for σ(pp(p) → J/ψ + J/ψ + X) × B 2 (J/ψ → µµ) between the LHCb, CMS and D0 data and our predictions in the relevant kinematical regions (+ that of the forthcoming ATLAS analysis). The theory predictions are: the SPS prompt yields at LO and NLO [For LHCb, the evaluation is a complete NLO [22]], the DPS prompt yields with σ eff fitted to the CMS differential distributions (see section 4.3) and the χ 2 between the sum of DPS+SPS (resp. DPS) yield and CMS and LHCb (resp. D0 DPS) data. For the DPS yields, the first uncertainty is from σ eff (see Table (2)) and the second in parenthesis is a systematical certainty from the 3 fits (alike the variation of the central value of σ eff in Table (2)). The range of the χ 2 also comes from the 3 fits. is almost one order of magnitude larger than LO one when P ψ T = 5 GeV and that one must use a NLO (or NLO ) evaluation when dealing with data sets with a P T cut as it is the case here for all but LHCb data.
Since both J/ψ should have their P T > 4 GeV, it is interesting to look at the impact of the α 5 s corrections (NLO ). Whereas the LO SPS yield is a bit below the SPS D0 yield, the NLO yield, which is about 3 times larger, is above. Both agree with the data within the large theoretical uncertainties -mainly from m c . As we shall see later, the need for α 5 s corrections is far more obvious in the CMS acceptance with a higher P T cut.
Injecting their measured σ ψ and σ DPS ψψ in Eq. (1), D0 has found σ eff 5.0 ± 2.75 mb. This value is 3 times lower than those extracted with jet observables, which means that the DPS yield seems to be 3 times higher than what could have been naively expected -or at least twice higher accounting from their uncertainty. Yet, as next section will show, a low value of σ eff allows one to solve the CMS puzzle.

The CMS data at large momenta
LO and NLO SPS cross sections for prompt J/ψ pair production in the CMS acceptance [35] are given in Table (1). As expected because of the higher P T cut, one observes a larger NLO /LO ratio than in the D0 acceptance. Yet the NLO SPS yield is significantly below the CMS data [35] and hint at the presence of another source of J/ψ pairs. As we whall see, the discrepancy is much more evident when one looks at differential distributions.
Indeed, besides this integrated yield, CMS measured differential distributions [35] which further indicate the importance of both NLO SPSs and DPSs but in different regions. Overall CMS released in addition 17 data points of differential cross sections as a function of P ψψ T , |∆y| and M ψψ . To  quantify the impact of the DPS contributions, we have used these experimental data to fit σ eff via Eq. (1) using the 3 fits of σ ψ discussed in section 2.2 and subtracting our theoretical evaluations of the SPS NLO yield acounting for their uncertainties (green band in the plot of Fig. 4). Table (2) summarises the fit result: the χ 2 d.o.f. and the values of σ eff along with their uncertainites coming from (a) the CMS experimental uncertainties 4 and (b) the theoretical uncertainties on the SPS yield. We have also given the χ 2 d.o.f. when no DPS contribution is considered. We note that the goodness of the 3 fits is similar. The dispersion of the central values of σ eff thus allows us to assess a systematical uncertainty due to the paramatrisation of σ ψ . Fig. 4 shows the DPS distributions with the Fit 2. Comparison plots using the Fit 1 & 3 are given a supplementary materials. In addition, we note that σ eff fitted using the Fit 3 (only Tevatron data) is very close to the D0 value. Whereas the p-values 5 of our DPS fits of σ eff are about 2% (see also below), the one without DPS is below 0.03% (even less without α 5 s contributions).
As regards the J/ψ-pair P T , P ψψ T , distribution, to the SPS yield. Phenomenologically accounting for the initial parton k T is not sufficient to reproduce the data as the smeared LO curve shows. At NLO, hard real-emissions tend to generate larger momentum imbalances and configurations with P ψψ T 0. In fact, the real-emission topologies (Fig. 1b) tend to produce, at large P T , two near J/ψ -as opposed to back-to-back-with a large P ψψ T . In the case of DPSs, correlations are absent and configurations at low P ψψ T are favoured. There is also no reason for large P ψψ T configurations to be -relatively-enhanced. It is thus not surprising that the DPS band drops faster than the NLO SPS one at large P ψψ T . The "bump" around P ψψ T 12 GeV simply reflects the kinematic cuts in the CMS acceptance. Overall, one obtains a good agreement (χ 2 d.o.f.

1.1) with the P ψψ
T distribution when DPS and (NLO ) SPS contributions are considered together, but it also confirms the dominance of (NLO ) SPS contributions at large P ψψ T . In addition, CMS analysed the relative-rapidity spectrum, dσ/d|∆y|. Along the lines of the D0 data discussion, the SPS contribution dominate when |∆y| → 0, while the DPS ones are several orders of magnitude larger than the SPS ones at large |∆y|. A comparison with the CMS data is shown in Fig. 3b. Most of the data are consistent with our results, except for the last bin, which probably explains the low p-value which we obtained with the DPS fit. For this data set, χ 2 d.o.f. 2.1 with DPS for the 3 fits and about 2.6 without (only NLO ). Note the large NLO SPS uncertainty which de facto reduces the corresponding χ 2 . The CMS acceptance with a rapidity-dependent P T cut renders dσ/d|∆y| flatter but this effect is apparently not marked enough in our theory curves. More data are however needed to confirm the absence of binning effects which could have generated a dip in the distribution. As another possible cross-check, we provide predictions for the ATLAS, D0, and LHCb acceptances in Table (1) and as supplementary materials.
At P ψψ T = 0 -where the bulk of the yield lies-, the J/ψpair invariant mass, M ψψ , is closely related to |∆y| and provides similar information (see Fig. 3c). One indeed has M ψψ = 2m ψ T cosh ∆y 2 . Large ∆y -i.e. large relative longitudinal momenta-correspond to large M ψψ . [At ∆y = 3.5 and P T = 6 GeV, M ψψ 40 GeV.] Without additional cuts, the M ψψ and dσ/d|∆y| spectra of the CMS do reveal the same conclusion: the DPS contributions dominate the region of large momentum differences. At small M ψψ , SPS contributions dominate and NLO corrections are large -essentially because CMS data do not cover low P ψ T . For this data set,

Back to D0 and LHCb data
As we have just seen, the inclusion of the DPS contributions with σ eff ranging from 5 to 11 mb, depending on the fit used for σ ψ , solves the so-called CMS puzzle found (with SPS contributions only) in [22] with a χ 2 d.o.f. reduced for all the distributions. Using these fits of σ eff , we should also reproduce the D0 DPS rate [12]. The agreement (see Table (1)) is quite good (χ 2 < 1) and gives us confidence that our proposed solution to the CMS puzzle is indeed correct. The comparison with the LHCb result is less instructive owing to the uncertainties from the data and from the SPS.

A first step toward a NNLO evaluation of the SPS contributions
At high P ψ T , O(α 6 s ) (NNLO) contributions like gg → cc followed by c(→ J/ψ + J/ψ + c) (Fig. 1c) or twice by c(→ J/ψ + c) (Fig. 1d) and gg → g g → (J/ψ + J/ψ + g)g (Fig. 1g) are expected to be enhanced by factors of P ψ T w.r.t NLO [21,22]. A two-loop computation is however needed to evaluate them, which is beyond the state-of-the-art. Yet, such leading-P T contributions can be evaluated via the fragmentation approximation, as done in [22] only for topologies like Fig. 1d, which were expected to dominate at large M ψψ . However, such an approximation has been shown [52] to be unjustified for the similar process gg → J/ψcc unless P T is much larger than m ψ . As discussed in section 2, the process gg → J/ψJ/ψ + cc is infrared safe and can be computed by itself using HELAC-ONIA "out-of-the-box" [36]. As just said, it includes a class of likely dominant NNLO corrections depicted in Fig. 1c and Fig. 1d. On the contrary, the contributions from topologies Fig. 1g could be evaluated using HELAC-ONIA but only with an infrared cutoff as in [53,21] for the NNLO . In the present study, we prefer to limit ourselves to gg → J/ψJ/ψ + cc which does not require any ad-hoc prescription.
The band labelled ψccψ in Fig. 4 shows its full contribution, which is computed for the first time. This partial α 6 s contribution is as large as the NLO ones only at the highest M ψψ and ∆y, where the DPS ones are anyhow dominant. The case of another variable, the sub-leading P T , P T min , is however particular since the DPS spectrum is expected to scale as P −2 × n T min , n being the scaling power of the single J/ψ yield (P −n T ). We thus found J/ψJ/ψ + cc to be dominant (see Fig. B.5 in the Appendix B) at very high P T min , which corresponds to back-to-back production as in Fig. 1d.
Overall, the aforementioned missing fragmentation contributions (Fig. 1g) at O(α 6 s ) are expected to be of similar sizes. In general, a full NNLO computation is thus expected to be similar to one at NLO accuracy, except in kinematical regions which are not currently accessible. Corresponding predictions for the forthcoming ATLAS and LHCb analyses as well as the current D0 acceptance are given as supplementary material for comparison with forthcoming data.

Possible impact of colour-octet transitions
We have also investigated the possible impact of CO channels as discussed at LO in [16,17]. We found that, because of the double suppression of the CO LDMEs, CO+CO channels are nowhere important when P ψ T < 50 GeV, as we found out [21]. We have evaluated the contribution from 3 S [8] 1 + 3 S [8] 1 and 1 S [8] 0 + 1 S [8] 0 (see Fig. 4) using the 1-σ upper value of the 3 S [8] 1 (resp. 1 S [8] 0 ) LDMEs of the NLO prompt fit of [54], i.e. 0.00283 GeV 3 (resp. 0.0541 GeV 3 ), these are compatible with the LO direct fit of [55] and are the only ones not dramatically overshooting the low-P T single J/ψ data [56]. As we look for an upper value, we disregard the 3 P [8] J + 3 P [8] J contribution which is negative. A complete CO study is beyond the scope of this Letter and will be the object a dedicated publication.
In any case, this upper limit of the CO contributions is always smaller than the CS ones except for |∆y| > 2.5 (last bin in Fig. 3b) and M ψψ > 40 GeV (last two bins in Fig. 3c). In these regions, these SPS contributions are anyhow extremely suppressed as compared to DPS ones and the CS SPS can receive significant α 6 s contributions. The only distribution where the CO contributions might show up is for the P T min (Fig. B.5); it has a similar size and the same de-pendence as our partial NNLO evaluation. They are larger than the NLO SPS and DPS yields only where the cross sections are on the order of 10 −8 nb before accounting for the branchings. As regards the mixed CO+CS channels, there is no P ψ T enhancement to be expected and these are simply suppressed by the LDME.

Conclusion
In this Letter, we have focused on the explanation of the recent observations of (prompt) J/ψ-pair production made by D0 and CMS. The measurements by CMS [35], which severely overshoot the theory if one solely considers SPS contributions [22], indicate significant DPS contributions, which we find to agree with the magnitude measured by D0 [12]. For the first time, our study shows that both DPSs and the NLO QCD corrections to SPSs are crucial to account for the existing data.
If these experimental results are confirmed, this would provide evidence for We have also derived generic formulae predicting feeddown contributions or, equally speaking, charmonium-pairproduction rates involving excited states, in case DPSs dominate. These can be checked by measuring J/ψ + ψ or J/ψ + χ c production. Such data can also therefore check a possible DPS dominance as found by D0 and CMS at large momenta. The relatively small value of σ eff (see Fig. 4) compared to jet-related extractions we obtained to describe the CMS data -also compatible with the D0 DPS yield and their extracted value of σ eff -may be a first hint at the flavour dependence of this effective cross section. Finally, we have carried out the first exact evaluation of leading-P T NNLO (α 6 s ) contributions, i.e. J/ψ-pair production with a cc pair, which could contribute at large P T min . On the way, our study of the impact of all the real emissions at α 5 s and some at α 6 s also demonstrates the absence of a significant colour-octet contribution contrary to what was found at LO in [16,17].

Supplementary materials
As supplementary materials, we gather (i) comparisons of the DPS cross section obtained with the Fit 1, 2 & 3 and (ii) predictions for forthcoming analyses of data on tapes. We present our predictions (i) for the ATLAS acceptance in Fig. .