Energy dependence of direct-quarkonium production in $\varvec{pp}$ collisions from fixed-target to LHC energies: complete one-loop analysis

Abstract

We compute the energy dependence of the $P_T$-integrated cross section of directly produced quarkonia in pp collisions at next-to-leading order (NLO), namely up to $\alpha _S^3$, within nonrelativistic QCD (NRQCD), treating the $P_T$-integrated and the $P_T$-differential cross sections as two different observables. The colour-octet NRQCD parameters needed to predict the $P_T$-integrated yield can thus be extracted from the fits of the $P_T$-differential cross sections at mid and large $P_T$. For the first time, the total cross section is evaluated in NRQCD at full NLO accuracy using the recent NLO fits of the $P_T$-differential yields. Both the normalisation and the energy dependence of the $J/\psi$, $\psi '$ and $\Upsilon (1S)$ we obtained disagree with the data except when using the fit results of Butenschoen and Kniehl. If one disregards the colour-octet contribution, the existing data in the TeV range are well described by the $\alpha _S^3$ contribution in the colour-singlet model – which, at $\alpha _S^4$, however, shows an unphysical energy dependence. A similar observation is made for $\eta _{c,b}$. All this underlines the necessity for a resummation of initial-state radiations in both channels, which is, however, beyond the scope of this article. In any case, past claims that colour-octet transitions are dominantly responsible for low-$P_T$ quarkonium production are not supported by our results.

1 Introduction

Understanding the production mechanism of low-$P_T$ quarkonia in nucleon–nucleon collisions is of fundamental importance to properly use them as probes of deconfinement or collectivity in heavy ion collisions. Indeed, most of the analyses of quarkonium production in nucleus–nucleus collisions are carried out on the bulk of the cross section, namely at low $P_T$.

Recently, comparisons between ALICE data [1] without $P_T$ cut and CMS data [2] with $P_T$ cut in PbPb collisions at $\sqrt{s_{NN}}=2.76$ TeV showed an unexpected suppression pattern of the charmonia, at variance with the simple picture of quarkonium melting in deconfined quark matter [3]. However, to properly interpret this observation, it is essential to rule out the possibility that a part of the effect observed could be due to a difference in the production mechanism in individual nucleon–nucleon collisions at low and at larger $P_T$. The propagation of a colour-octet pair in a deconfined medium certainly differs from that of a colour-singlet pair; this can result in a different nuclear suppression (see e.g. [6]). On the contrary, as regards the bottomonia, the observation of the expected sequential-suppression pattern has been claimed by CMS [4, 5].

Further, the effect of normal nuclear matter may also significantly depend on how the pair is produced: the recently revived fractional energy loss [7, 8] would for instance act on long-lived colour-octet states and probably differently if the heavy-quark state is already produced colourless at short distances, as postulated in the CSM [9]. Saturation effects in pA collisions also do depend on the colour state of the perturbatively produced heavy-quark pair [10–12].

Despite the possibility that NRQCD factorisation would not hold at low $P_T$, several NRQCD analyses have thus been carried earlier to evaluate the impact of the colour-octet channels to the $P_T$-integrated $J/\psi$ yields [13–15]. A first study of the impact of initial-state radiation (ISR) on the very low $P_T$ $J/\psi$’s and $\Upsilon$’s was recently carried out successfully in NRQCD [16] – yet at the cost of introducing additional non-perturbative parameters.

Whereas, based on an analysis of the sole early RHIC data, Cooper et al. argued [14] that the universality of NRQCD was safe and that colour-singlet contributions to the $P_T$-integrated $J/\psi$ yields were negligible, the global analysis of Maltoni et al. at NLO showed [15] that the colour-octet Long-Distance Matrix Elements (LDMEs) required to describe the total prompt $J/\psi$ yield from fixed-target energies to RHIC were one tenth of that expected from the – leading-order – fit of the $P_T$-differential cross sections at Tevatron energies.

Such fits of the $P_T$-differential $J/\psi$ cross sections have recently been extended to NLO – i.e.  one-loop – accuracy on the prompt $J/\psi$ yields – some of them focusing on the larger $P_T$ data and explicitly including the feed-down contributions [17, 18], some enlarging the analysis beyond hadroproduction and including rather low-$P_T$ data [19] – and on the $\Upsilon (nS)$ yields [20, 21]. Thanks to these studies, we can significantly extend the existing NRQCD studies of the $P_T$-integrated cross section by combining in a coherent manner, the hard parts – or Wilson coefficients – up to $\alpha _S^3$, first derived in [22], and which we have systematically checked with FDC [23], with the NRQCD matrix elements fitted at NLO on the $P_T$ dependence of the yields. One can indeed consider the $P_T$ -integrated and the $P_T$-differential cross sections as two different observables – their Born contributions involve different diagrams – and such a procedure is not at all trivial from the physics point of view.

As we shall discuss in detail later, our results show that the data do not allow for a global description of both the $P_T$-integrated and the $P_T$-differential quarkonium yields. As a point of comparison, we also had a look at Colour-Evaporation-Model-like (CEM) predictions derived from NRQCD following the work of [26] and we found that it cannot reproduce $P_T$-integrated yields using the LDMEs obtained following the relations of [26] after identifying the minimal singlet transition to that of the CSM. A contrario, results obtained from the traditional CEM implementation at one loop do not show a similar issue.

The inability of colour-octet dominance within NRQCD to provide a global description of both low and large $P_T$ data is in line with the recent findings [28–31] that the sole LO colour-singlet contributions are sufficient to account for the magnitude of the total cross section and its dependence in rapidity, $\mathrm{d}\sigma /\mathrm{d}y$, from RHIC, Tevatron all the way to LHC energies. Any additional contribution in this energy range creates a surplus¹ as compared to data.

However, as we also study in a dedicated section, the total NLO CSM cross section shows a weird energy dependence at LHC energies. The problem is striking for the $J/\psi$, less for the $\Upsilon$. In any case, one should be very careful in interpreting these results. In particular, such NLO results cannot be considered as a improvement of the LO ones. We also observe the same issue for $\eta _c$ and $\eta _b$ production for which there is no final-state-gluon radiation at Born order. We are therefore tempted to attribute this behaviour to large loop contributions which become negative at low $P_T$, rather than to specific effects related to the $^3S_1$ production per se. A quick inspection of the rather concise one-loop results [32] for $\eta _c$ and $\eta _b$ production in the TMD factorisation formalism unfortunately does not reveal any obvious negative contributions and does not help in the understanding of this rather general issue of quarkonium production in collinear factorisation.

The same problems appear with some CO channels as well and may therefore cast doubts on the reliability of our results in the $\sqrt{s}$ region where some contributions shows a strange behaviour – in particular at LHC energies. At this stage, we are not able to conclude from our observations whether these problems are indicative of the breakdown of NRQCD factorisation at low $P_T$ and low x or at low x only. However, for sure, none of the above observations can reasonably support the idea that CO transitions are dominant at low $P_T$. Such a conclusion would at least be premature.

The structure of the paper is as follows. In Sect. 2, we detail the procedure to evaluate the $P_T$-integrated yield at one-loop accuracy in NRQCD and we explain the idea underlying this first complete one-loop analysis. In Sect. 3, we explain our selection of LDMEs determined at NLO. In Sect. 4, we briefly comment on the existing world data sets for $J/\psi$, $\psi (2S)$ and $\Upsilon (1S)$.² We also explain how we estimate the direct yields. In Sect. 5, we show and discuss our results for the first full one-loop NRQCD analysis of quarkonium hadroproduction. To go further in the interpretation of some of our results, we discuss in Sect. 6 the prediction of NRQCD using CEM-like LDMEs. This is also compared with the conventional approach based on quark–hadron duality. Section 7 focuses on CSM results both for the $^3S_1$ states considered here and for the $^1S_0$ states for which analytical results exist. Our conclusions are presented in Sect. 8.

2 A full one-loop cross-section computation

2.1 Generalities

Following the NRQCD factorisation, the cross section for quarkonium hadroproduction can be expressed from the parton densities in the colliding hadrons, f(x), a hard-part – the partonic cross section – for the production of a heavy-quark pair with zero relative velocity, v, in a definite angular-momentum, spin and colour state, and a LDME connected to the hadronisation probability of the intermediate state into the quarkonium. Namely, one has for the production of a quarkonium $\mathcal{Q}$, along with some unidentified set of particles X,

$$\begin{aligned} \sigma =\!\!\sum _{i,j,n}\!\!\int \!\! \mathrm{d}x_{1}\mathrm{d}x_{2}f_{i}f_{j}\hat{\sigma }[i+j\rightarrow (Q\overline{Q})_{n}\!+\!X]\langle \mathcal{O}_{\mathcal{Q}}(n)\rangle \end{aligned}$$

(1)

where the indices i, j run over all partonic species and n denotes the colour, spin and angular-momentum states of the intermediate $Q\overline{Q}$ pair.

For the $^3S_1$ quarkonium states, the first CO states which appear in the v expansion are the $^{1}\!S^{[8]}_{0},^{3}S^{[8]}_{1}$ and $^{3}\!P^{[8]}_{J}$ states, in addition to the leading v contribution $^{3}S^{[1]}_{1}$ from a CS transition. One, however, has to note that, for hadroproduction, whereas the CO contributions already appear at $\alpha _S^2$ (Fig. 1a, b), the CS one only appear at $\alpha _S^3$ (Fig. 1c). These $\alpha _S^2$ CO graphs nevertheless do not contribute to the production of quarkonia with a nonzero $P_T$, since they would be produced alone without any other hard particle to recoil on.

Fig. 1

Representative diagrams contributing a, b at Born order to $i+j\rightarrow \mathcal{Q}$, c–e both at Born order to $i+j\rightarrow \mathcal{Q}+$jet and at one loop to $i+j\rightarrow \mathcal{Q}$, f at one loop to $i+j\rightarrow \mathcal{Q}$, g–k at one loop to $i+j\rightarrow \mathcal{Q}+$jet

The Born contributions from CS and CO transitions are indeed different in nature: the former is the production of a quarkonium in association with a recoiling gluon, which could form a jet, while the latter is the production of a quarkonium essentially alone at low $P_T$.

Let us now have a look at the $\alpha _S^3$ CO contributions (Fig. 1d–f) which are then NLO – or one-loop – corrections to quarkonium production and which are potentially plagued by the typical divergences of radiative corrections. Yet, the real-emission $\alpha _S^3$ corrections to CO contributions (Fig. 1d, e) can also be seen as Born-order contributions to the production of a quarkonium + a jet – or, to put it otherwise, of a quarkonium with $P_T \gg \Lambda _\mathrm{QCD}$. As such, they do not show any soft divergences for $P_T \ne 0$. These are supposed to be the leading contribution to the $P_T$-differential cross section in most of the data set taken at a hadron collider (Tevatron, RHIC and LHC). These are now known up to one-loop accuracy, namely up to $\alpha _S^4$ (see e.g. [17–21, 34, 35]) (Fig. 1g, h).

It is important to note that one cannot avoid dealing with the divergences appearing at $\alpha _S^3$ if one study the $P_T$-integrated cross section.

2.2 Different contributions up to $\alpha ^3_S$

At $\alpha _S^2$, the CO partonic processes are

$$\begin{aligned} q+\bar{q}\rightarrow & {} Q \bar{Q}[^{3}S^{[8]}_{1}] \quad \hbox { (Fig.~1(a))} \nonumber \\ g+g\rightarrow & {} Q\bar{Q}[^{1}\!S^{[8]}_{0}, ^{3}\!P^{[8]}_{J=0,2}] \quad \hbox { (Fig.~1(b))} \end{aligned}$$

(2)

where q denotes u, d, s.

At $\alpha _S^3$, the QCD corrections to the aforementioned channels include real (Fig. 1d, f) and virtual (Fig. 1f) corrections. One encounters UV, IR and Coulomb singularities in the calculation of the virtual corrections. The UV-divergences from the self-energy and triangle diagrams are removed by the renormalisation procedure. Since we follow the same lines as [18, 35] where all the procedure is described, we do not repeat its description. As regards the real-emission corrections, they arise from three kinds of processes (not all drawn):

$$\begin{aligned} \begin{aligned}&g+g\rightarrow Q \bar{Q}[^{1}\!S^{[8]}_{0},^{3}\!\!S^{[8]}_{1},^{3}\!\!P^{[8]}_{J=0,2}]+g, \\&g+q(\overline{q})\rightarrow Q \bar{Q}[^{1}\!S^{[8]}_{0},^{3}\!\!S^{[8]}_{1},^{3}\!\!P^{[8]}_{J=0,2}]+q(\overline{q}), \\&q+\overline{q}\rightarrow Q \bar{Q}[^{1}\!S^{[8]}_{0},^{3}\!\!S^{[8]}_{1},^{3}\!\!P^{[8]}_{J=0,1,2}]+g. \end{aligned} \end{aligned}$$

(3)

As usual, the phase-space integrations generate IR singularities, which are either soft or collinear and can be conveniently isolated by slicing the phase space into different regions. Here we adopt the two-cutoff phase-space slicing method to deal with the problem [36].

As we previously alluded to, the $\alpha _S^3$ CS contribution is particular since, in the limit $v=0$, it would be strictly speaking of Born order for both the production of a quarkonium and of a quarkonium + a jet. It arises from the well-known process

$$\begin{aligned}&g+g\rightarrow Q\bar{Q}[^{3}\!S^{[1]}_{1}]+g \quad \hbox { (Fig.~1(c))}. \end{aligned}$$

(4)

Our calculations are equivalent to previous work by Maltoni et al. [15, 22] and we have checked that we reproduce their results for all the relevant channels. As announced in the introduction, one of the novelties in our study resides in the use of the LDMEs fitted at the same order, i.e. one loop, to the $P_T$-differential cross sections. As such, this is the first global NLO analysis of hadroproduction.

Since we also look at data at rather low energies, we also included a CS channel via $\gamma ^\star$ exchange. Indeed, as noted in a different context in [37], the QED CS contributions via $\gamma ^\star$ are naturally as large as the CO $^3S_1^{[8]}$ transition via $g^\star$ – the $\alpha _\mathrm{em}$ suppression being compensated by the small relative size of the $^3S_1^{[8]}$ CO LDME ($\mathcal{O}(10^{-3})$) as compared to the $^3S_1^{[1]}$ CS LDME ($\mathcal{O}(1)$). The real-emission contributions arise from

$$\begin{aligned} \begin{aligned}&q+\overline{q}\rightarrow Q \bar{Q}[^{3}\!S^{[1]}_{1}]+g,\\&g+q(\overline{q})\rightarrow Q \bar{Q}[^{3}\!S^{[1]}_{1}] +q(\overline{q}), \end{aligned} \end{aligned}$$

(5)

whereas the loop contributions are only from

$$\begin{aligned}&q+\overline{q}\rightarrow Q \bar{Q}[^{3}\!S^{[1]}_{1}]. \end{aligned}$$

(6)

Figure 1a (Fig. 1f) with the s-channel gluon replaced by a $\gamma ^\star$ would depict the Born (a one-loop) contribution. We have, however, found that they do not matter in the regions which we considered.

3 Constraints on the LDMEs from the $P_T$-differential cross section

The CS LDMEs can either be extracted from the leptonic decay width at NLO or can be estimated by using a potential model result, which gives for the Buchmüller–Tye potential [38] $|R_{J/\psi }(0)|^{2}=0.81$ GeV$^{3}$, $|R_{\psi (2S)}(0)|^{2}=0.53$ GeV$^{3}$ and $|R_{\Upsilon (1S)}(0)|^{2}=6.5$ GeV$^{3}$.

As regards the CO LDMEs, they can only be extracted from data. As we discussed above, our aim is to analyse the $P_T$-integrated yield using the constraints from the $P_T$ dependence of the yields.

3.1 $J/\psi$

In the $J/\psi$ case, we will use the results of five fits of this dependence [17, 18, 39–41]. The first two were limited to pp data but explicitly took into account the effect of the feed-down.³ The latter fit was based on a wider set of data including ep and $\gamma \gamma$ systems but the feed-down effects were only implicitly included through constant fractions for these systems. The fourth one includes the recent $\eta _c$ measurement at $P_T \ge 6$ GeV by LHCb [42] by relying on the heavy-quark spin symmetry of NRQCD, which relates colour-octet matrix elements of spin-singlet and -triplet quarkonia with the same principal quantum number n. The fifth one incorporates the leading-power fragmentation corrections together with the usual NLO corrections, which results in a different short-distance coefficient and allows for different LDMEs.

Another recent fit [43] took the $\eta _c$ measurement into account. The LDME values which they found fall into the range considered for [40], therefore we do not use it separately.

In Ref. [17], Ma et al. have based their analyses on the fit of two linear combinations⁴ of LDMEs:

$$\begin{aligned} M^{J/\psi }_{0,\,r_{0}}= & {} \langle \mathcal{O}_{J/\psi }(^{1}\!S^{[8]}_{0})\rangle + \frac{r_{0}}{m^2_c}\langle \mathcal{O}_{J/\psi }(^3\!P^{[8]}_0)\rangle , \nonumber \\ M^{J/\psi }_{1,\,r_{1}}= & {} \langle \mathcal{O}_{J/\psi }(^{3}\!S^{[8]}_{1})\rangle + \frac{r_{1}}{m^2_c}\langle \mathcal{O}_{J/\psi }(^3\!P^{[8]}_0)\rangle . \end{aligned}$$

(7)

They proceeded to two fits with different $P_T$ cuts. We use that for $P_T>7$ GeV and limit ourselves to the central values they obtained: $M^{J/\psi }_{0,\,r_{0}}=7.4\times 10^{-2}$ GeV$^3$ and $M^{J/\psi }_{1,\,r_{1}}=0.05\times 10^{-2}$ GeV$^3$, since a single set of values of $M^{J/\psi }_{0,r_{0}}$ and $M^{J/\psi }_{1,r_{1}}$ translates anyhow into a wide range of values of the LDMEs. Indeed, limiting ourselves to positive values of $\langle \mathcal{O}_{J/\psi }(^{1}\!S^{[8]}_{0})\rangle$ and $\langle \mathcal{O}_{J/\psi }(^{3}\!S^{[8]}_{1})\rangle$, one can solve Eq. 7 and get the loose constraint $\langle \mathcal{O}_{J/\psi }(^3\!P^{[8]}_0)\rangle \in [-0.2, 4.3]\times 10^{-2}$ GeV$^5$. As a central value, we choose the middle of the allowed interval. The same group has, however, recently improved their analysis by taking into account the feed-down [44]. As aforementioned, they in turn performed a new fit [40] including $\eta _c$ data. The six sets of LDMEs to be used to probe the allowable parameter space of the fit are given in Table 1.

Table 1 Values of $\langle \mathcal{O}_{J/\psi }(^{3}\!P^{[8]}_{0})\rangle$, $\langle \mathcal{O}_{J/\psi }(^{1}\!S^{[8]}_{0})\rangle$ and $\langle \mathcal{O}_{J/\psi }(^{3}\!S^{[8]}_{1})\rangle$ from five NLO (i.e. at one loop) fits of the $P_T$ dependence of the yields, which we will use to compute [the energy dependence of] the $P_T$-integrated yields

As mentioned above, in [39], Butenschoen et al. proceeded to a global fit of prompt $J/\psi$ data from pp, $\gamma \gamma$, $\gamma p$ systems.⁵ Since $\gamma \gamma$, $\gamma p$ mostly lies at low $P_T$, they also considered data at rather low $P_T$ from RHIC. They did not include NLO predictions for $\chi _c$ in the fit. Rather they assumed a constant direct fraction, for instance 36 % for hadroproduction.

3.2 $\psi (2S)$

Buttenschoen et al. did not provide a fit of $\psi (2S)$ in [39] due to the lack of data besides those from pp collisions. The LDMEs which we consider for $\psi (2S)$ are therefore only from [17, 18]. For the former fit, the values are obtained in the same way as for the $J/\psi$, where $M^{\psi (2S)}_{0,\,r_{0}}=2.0\times 10^{-2}$ GeV$^3$ and $M^{\psi (2S)}_{1,\,r_{1}}=0.12\times 10^{-2}$ GeV$^3$. The resulting values as well as those from [18] are gathered in Table 2.

Table 2 Same as Table 1 for $\psi (2S)$

In [44], the authors of [17] tried to refit the existing data with a larger $P_T$ cutoff. Such a fit already badly overshoots mid-$P_T$ data. We therefore do not consider it in this work. For the same reason, we have not considered the fit of [45] since it only reproduces the $\psi (2S)$ data in an admittedly narrow – high $P_T$ – range.

3.3 $\Upsilon (1S)$

As regards the $\Upsilon (1S)$, there are two NLO analyses from [20, 21]. However, Wang et al.  used in [20] a different value of the NRQCD factorisation scale $\mu _\Lambda$, which we use in the present evaluation, that is, $\mu _\Lambda =m_b$. To perform a correct comparison would have required a new evaluation of the hard coefficients with their choice of $\mu _\Lambda$ to use their LDME values. In addition, although they did consider the effects of excited feed-down, they have not disentangled the direct contribution to that of the feed-down in their LDME extraction. The central values of [21] are gathered in Table 3.

Table 3 Same as Table 1 for $\Upsilon (1S)$

4 World data and feed-down effects

As regards the data for $J/\psi$ and $\psi (2S)$, we drew on the extensive set used in [15] with the exception that we only kept data:

• derived from more than 100 events at a given $\sqrt{s}$;

• from pp or $p\bar{p}$ collisions only in order to avoid dealing with nuclear effects;

• where $\mathrm{d}\sigma /\mathrm{d}y$ was derived at $y=0$.

To this set, we have added data published later than 2006, which includes data from the LHC. We have also added one point from the CDF Collaboration at the Tevatron.⁶ All the quoted uncertainties are combined in quadrature together with that of the direct fraction,⁷ which we assumed to be energy independent and $F^\mathrm{direct}_{J/\psi } = 60\pm 10\,\%$ [28].

As regards the $\psi (2S)$, the data sets are very scarce, especially if one focuses on $P_T$-integrated yields at $y=0$. In fact, there is only data from ISR-Clark et al. [47] averaged over $\sqrt{s}=52.4$ and 62.7 GeV and from PHENIX at $\sqrt{s}=200$ GeV. CDF measured the cross section at $\sqrt{s}=1.96$ TeV for $|y| <0.6$, but only for $P_T > 2$ GeV [53]. In order to use this precise measurement, we have extrapolated it by assuming the same ratio $\frac{\mathrm{d} \sigma (P_T < 2 \mathrm{GeV})}{\mathrm{d}y}|_{y=0}/\frac{\mathrm{d} \sigma (P_T > 2 \mathrm{GeV})}{\mathrm{d}y}|_{y=0}=0.82$ as for the $J/\psi$ [52]. As for now, there does not exist a measurement at LHC energies in the central rapidities down to small enough $P_T$ to perform a model-independent enough extrapolation.⁸

As regards the $\Upsilon (1S)$, the data set is surprisingly wider than that of $\psi (2S)$ despite a significantly smaller production cross section. It is certainly due to the larger energy of the decay leptons and to the smaller background. For a long time, it was considered that only half of the (low-$P_T$) $\Upsilon (1S)$ were directly produced ($F^\mathrm{direct}_{\Upsilon (1S)} = 50\pm 10\,\%$) based on an early CDF measurement [61]. Recent LHCb studies of $\chi _b$ production [33, 62], along with $\Upsilon (2S,3S)$ cross section measurements [59, 60, 63, 64], rather indicate that two thirds of the $\Upsilon (1S)$ are directly produced, and we will therefore opt for $F^\mathrm{direct}_{\Upsilon (1S)} = 66\pm 10\,\%$. Yet, a number of experiments could not resolve the three $\Upsilon$ states. In this case, one should apply [28] a slightly smaller direct fraction, which we take to be $F^\mathrm{direct}_{\Upsilon (1S+2S+3S)} = 60\pm 10\,\%$. As we take this fraction to be energy independent, we choose a conservative estimate of their uncertainty.

Table 4 shows the $J/\psi$ data set, Table 5 that of $\psi (2S)$ and Table 6 that of $\Upsilon$.

Table 4 $J/\psi$ data set used in our data–theory comparison. The experimental values quoted in the experimental papers may have been multiplied by feed-down factors (see text)

Table 5 $\psi (2S)$ data set used in our data–theory comparison

Table 6 $\Upsilon (1S)$ data set used in our data–theory comparison. [A star indicates that the measurement could not resolve the 1S, 2S and 3S states.]

5 Complete NLO results within NRQCD

In the numerical computation at NLO, the CTEQ6M PDF [65],⁹ and the corresponding two-loop QCD coupling constant $\alpha _{s}$ are used.¹⁰ The charm quark mass, $m_c$, is set by default to 1.5 GeV and the bottom quark one, $m_b$, to 4.5 GeV. Our default choices for the renormalisation, factorisation and NRQCD scales are $\mu _{R}=\mu _{F}=\mu _{0}$ with $\mu _{0}=2m_{Q}$ and $\mu _{\Lambda }=m_{Q}$, respectively. When other choices are made, in particular to estimate the theoretical uncertainty due to the lack of knowledge of corrections beyond NLO, they are indicated on the corresponding plots. We have taken $\delta _{s}=10^{-3}$ and $\delta _{c}=\delta _{s}/50$ for the two phase-space cutoffs – the insensitivity of the result on the chosen values for these cutoff has been checked. Our results for direct $J/\psi$, $\psi (2S)$ and $\Upsilon (1S)$ are shown on, respectively, Fig. 2a–c.

Fig. 2

The cross section for direct a $J/\psi$, b $\psi (2S)$ and c $\Upsilon (1S)$ as a function of $\sqrt{s}$. The blue dot-dashed curve is the central CS curve. Its relative uncertainty is shown in the lower panels; the light green (light blue) band shows the scale (mass) uncertainty. The dashed red curve is the total CO contribution from three channels: $^{3}\!P^{[8]}_{0}$ (thin dot-dashed orange), $^{1}\!S^{[8]}_{0}$ (thin dotted magenta) and $^{3}\!S^{[8]}_{1}$ (thin dashed green). The total CO uncertainty relative to the CS central curve is shown in the lower panels; the light pink (purple) band shows the scale (mass) uncertainty. The black is the total contribution (CS + CO) at one loop. These are compared to experimental data (see text) multiplied by a direct fraction factor (when applicable) and normalised to the central CS curve in the lower panels. [Negative CO contributions are indicated by arrows.]

We first discuss the comparison between the five fits and the $J/\psi$ data (Fig. 2a). Unsurprisingly, our study shows that the global fits including rather low $P_T/m_\mathcal{Q}$ data, that is, the one of Butenschoen et al. [39], provides the only acceptable description of the $P_T$-integrated cross section. We, however, note that the latter fit does not provide a good description of the $J/\psi$ polarisation data and, as recently noted [67], it yields a negative cross section for $J/\psi +\gamma$ at large $P_T$. Finally, it does not allow [68] to describe the $\eta _c$ data. The fits of Gong et al. [18] and Ma et al. [17, 40] greatly overshoot the data in the energy range between RHIC and the Tevatron, whereas these fits a priori provide a good description of the $P_T$-differential cross section at these energies.

The fit of Bodwin et al. [41] gives the worse account of the $P_T$-integrated $J/\psi$ data in the whole energy range. Indeed, the new ingredient of [41] allows one to describe high-$P_T$ data with a large $^1S_0^{[8]}$ LDME (see Table 1) – as for [18] but without negative values for the other octet LDMEs – which results in too large a yield at low $P_T$.

In addition, we also note the strange energy dependence of at least the P-wave octet channel. We postpone its discussion to Sect. 5.1, where this is analysed in more detail for the $^1S_0^{[8]}$ transition, and to Sect. 7, where we discuss a similar observation for the CSM yield at NLO.

As regards the $\psi (2S)$ (Fig. 2b), our NLO NRQCD results do not reproduce the data at all at RHIC energies and, since both fits are dominated by the P-wave octet channel, show a nearly unphysical behaviour at LHC energies.

The comparison for the $\Upsilon (1S)$ (Fig. 2c) is more encouraging. At RHIC energies and below, the agreement is even quite good, while at Tevatron and LHC energies, the NLO NRQCD curves only overshoot the data by a factor of 2.

We finally note that from RHIC to LHC energies, the LO CSM contributions (the blue in all the plots) accounts well for the data. The agreement is a bit less good for $\psi (2S)$ if we stick only to the default/central value. This is not at all a surprise and is in line with the previous conclusions made in [28–31]. In fact, strangely enough, it seems that it is only at low energies (below $\sqrt{s}=100$ GeV) that the CO contributions would be needed to describe the data. The more recent data from the LHC and the Tevatron tend to agree more with the LO CSM.

Overall, this shows – unless the resummation of ISR modifies our predictions by a factor of 10 – that it would be difficult to achieve a global description of the total and $P_T$-differential yield and its polarisation at least for the charmonia.

As we discussed in the introduction, a first resummation study has recently been performed within NRQCD [16]. When combined with the results of [17], this resummation yields [16] a good description of the low-$P_T$ data. It should, however, be stressed that this study introduces three new parameters $g_{1,2,3}$ to parametrise the so-called $W^\mathrm{NP}$ function used the CSS resummation procedure. Moreover, we stress that such values of the CO LDMEs would result in a negative NLO $P_T$-differential cross section for $J/\psi +\gamma$ at large $P_T$ where NRQCD factorisation should normally hold.

To close the discussion of the theory–data comparison, let us note that the $^3\!S_1^{[8]}$ channel alone would provide a reasonable energy dependence. If we were to refit the low-$P_T$ data and thus obtain dominance of the $^3\!S_1^{[8]}$ channel, the yield at large $P_T$ would nevertheless dominantly be transversely polarised in disagreement with existing data (e.g. [69]). Yet, the better energy dependence of $^3\!S_1^{[8]}$ at NLO with respect to the other octet channels, which shows a flat energy dependence in the TeV region, is probably a fortunate “accident”. Indeed, most of the $^3\!S_1^{[8]}$ yield up $\alpha _s^4$ is in fact not at one loop, but for the – suppressed – $q\bar{q}$ contribution, since the $gg\rightarrow ^3\!S_1^{[8]}$ is zero. The curve for $^3\!S_1^{[8]}$ – as well as $^3\!S_1^{[1]}$ – shown on Fig. 2a is effectively a Born-order one. Clearly, these behave better than the channels where the loop contributions are allowed.

5.1 Behaviour at high $\sqrt{s}$ and scale dependence of the $^1S_0^{[8]}$ contribution

In view of the observations just made in the previous paragraphs, we have found it useful to analyse more carefully the behaviour of one specific channel. We have decided to look more carefully at the simplest one, that is, from the $^1S_0^{[8]}$ transition, in particular we will look at how the scale choices influence the behaviour of the yield at large $\sqrt{s}$. Analytical results for the hard-scattering partonic amplitude squared can be found in Appendix C of [22]. We have used them in a small numerical code to convolve it with PDFs and checked them with the results of FDC. The advantage of using FDC is that we can easily cut on $P_T$ and y.

As one could have anticipated from the band of the lower panel of Fig. 2a, one observes (Fig. 3) that for a wide range of scale choices, the energy dependence remains extremely flat. For some of the choices where $\mu _F > \mu _R$, one even sees that the cross section clearly decreases and becomes negative – when the yield becomes negative the curve stops. This is of course not satisfactory. At this stage, we are not able to conclude from our observations – normalisation and high-energy issue – whether these features point at the break down of NRQCD factorisation, NRQCD universality or should force us to continue questioning our understanding of the mid- and high-$P_T$ quarkonium-production mechanisms. To investigate this a bit further, we will look at the predictions of another approach – the colour evaporation model – in the next section and later in the colour-singlet model both for $^3S_1$ and $^1S_0$ quarkonia.

Fig. 3

The cross section for the production of a $J/\psi$ from only a colour-octet $^1S^{[8]}_0$ $c\bar{c}$ state as a function of the cms energy for various choice of the mass and scales

6 Colour-Evaporation-Model-like NRQCD evaluation

To go further in our investigations of QCD one-loop effects on the energy dependence of quarkonium production, we have found it useful to compare our results with those of the Colour-Evaporation Model (CEM) which directly follows from the quark–hadron duality [70, 71]. The quarkonium production cross section is obtained by considering the cross section to produce a $Q \bar{Q}$ pair in an invariant mass region compatible with its hadronisation into a quarkonium, namely between $2m_Q$ and the threshold to produce open heavy-flavour hadrons, $2m_{H}$. One should multiply this by a phenomenological factor accounting for the probability that the pair eventually hadronises into a given quarkonium state. Overall, one considers

$$\begin{aligned} \sigma ^\mathrm{(N)LO,\ direct }_\mathcal{Q}= \mathcal{P}^\mathrm{direct}_\mathcal{Q}\int _{2m_Q}^{2m_H} \frac{\mathrm{d}\sigma _{Q\bar{Q}}^\mathrm{(N)LO}}{\mathrm{d} m_{Q\bar{Q}}}\mathrm{d} m_{Q\bar{Q}}. \end{aligned}$$

(8)

In a sense, the factor $\mathcal{P}^\mathrm{direct}_\mathcal{Q}$, i.e. the probability for (or the fraction of) the $Q\bar{Q}$ pairs in the relevant invariant mass region to directly hadronise into $\mathcal Q$, plays a similar role as the LDMEs in NRQCD, except that its size can be guessed. Indeed, it is expected [72] that one ninth – one colour-singlet $Q\bar{Q}$ configuration out of nine possible – of the open charm cross section in this invariant mass region eventually hadronise into a “stable” quarkonium. Taking into account this factor 9, in the case of $J/\psi$, it was argued [72] that a simple statistical counting, which would give

$$\begin{aligned} \mathcal{P}^\mathrm{direct}_{J/\psi }= \frac{1}{9} \frac{2 J_\psi +1}{\sum _i (2 J_i +1)} = \frac{1}{45}, \end{aligned}$$

(9)

where the sum over i runs over all the charmonium states below the $D\bar{D}$ threshold, could describe the existing data in the late 1990s. The solid turquoise curve – computed at NLO as opposed to the analysis of [72] – on Fig. 4a illustrates this. ¹¹

Fig. 4

The cross section for direct a $J/\psi$ and b $\Upsilon (1S)$ as a function of $\sqrt{s}$ from NLO NRQCD using the CEM-like constrained LDMEs assuming a minimal singlet transition. It is compared to the existing experimental measurements (see text)

Following the fit of Vogt in [74], $\mathcal{P}^\mathrm{direct}_{J/\psi }$ lies between 1.5 and 2.5 %. This indeed remarkably coincides with the simple statistical counting.¹² For the $\Upsilon$, the corresponding quantity is of similar size, between 2 and 5 %, although following the state-counting argument, one may expect a smaller number than for $J/\psi$. Let us nevertheless stress that a violation of Eq. 9 cannot be used to invalidate the CEM, since this relation completely ignores phase-space constraints. What CEM predicts is that $\mathcal{P}^\mathrm{direct}_\mathcal{Q}$ is process independent.

In [26], Bodwin et al. studied the connexion between the CEM and NRQCD. Following [26] up to $v^2$ corrections, only four intermediate $Q\bar{Q}$ states contribute to $^3\!S_1$ quarkonium production in a CEM-like implementation of NRQCD. One indeed has

$$\begin{aligned} \langle \mathcal{O}_{^3\!S_1}(^{3}\!S^{[1]}_{1})\rangle= & {} 3 \times \langle \mathcal{O}_{^3\!S_1}(^{1}\!S^{[1]}_{0})\rangle ,\nonumber \\ \langle \mathcal{O}_{^3\!S_1}(^{1}\!S^{[8]}_{0})\rangle= & {} \frac{4}{3} \times \langle \mathcal{O}_{^3\!S_1}(^{1}\!S^{[1]}_{0})\rangle ,\nonumber \\ \langle \mathcal{O}_{^3\!S_1}(^{3}\!S^{[8]}_{1})\rangle= & {} 4 \times \langle \mathcal{O}_{^3\!S_1}(^{1}\!S^{[1]}_{0})\rangle . \end{aligned}$$

(10)

All these nonvanishing LDMEs are then fixed if one makes the reasonable assumption that $\langle \mathcal{O}_{^3\!S_1}(^{3}\!S^{[1]}_{1})\rangle$ is indeed the usual CS LDME, i.e. $\frac{2 N_C}{4 \pi }~(2J+1)~|R(0)|^2$. As compared to the results presented in the previous section, the only additional piece to perform a full one-loop analysis is the hard part for $^{1}\!S^{[1]}_{0}$, which normally does not need to be considered for $^3\!S_1$ production at this level of accuracy in v. Computing it with FDC [23] does not present any difficulty.

Figures 4a, b show the resulting cross sections of $J/\psi$ and $\Upsilon (1S)$ production for the relevant channels and their sum, to be compared to the world data set used in the previous section. By construction, the $^{3}\!S^{[1]}_{1}$ curve is the same as in the previous section. One notes that, as for the P-wave octet, the $^{1}\!S^{[1]}_{0}$ curve strangely flattens out in the $J/\psi$ case at high energies. We will come back to this in the next section.

The total CEM-like contribution greatly overshoots the data, by a factor as large as 100. This was to be expected, since (i) following Eq. 10, all the LDMEs are roughly of the same size, (ii) the $^{3}\!S^{[1]}_{1}$ is roughly compatible with the data and (iii) the hard parts for the other transitions appear at $\alpha _S^2$ and are not suppressed as $P_T \rightarrow 0$ and are thus expected to give a larger contribution than the $^{3}\!S^{[1]}_{1}$ if one disregards the LDME.

Of course, one could question our assumption that $\langle \mathcal{O}_{^3\!S_1}(^{3}\!S^{[1]}_{1})\rangle =\frac{2 N_C}{4 \pi }~(2J+1)~|R(0)|^2$ and rather fit $\langle \mathcal{O}_{^3\!S_1}(^{1}\!S^{[1]}_{0})\rangle$. In both the $J/\psi$ and the $\Upsilon (1S)$ cases, the corresponding LDMEs would then approximately be 100 times smaller. In particular, the singlet transition would be 100 times less probable than what one expects from the leptonic decay. This would be an unlikely and dramatic violation of factorisation, which should have implications elsewhere. In particular, a pair produced at short distances with the same quantum number as the physical state, among these the colour, would have a much larger probability than expected to be broken up before eventually hadronising.

Although it is not as obvious as in the NRQCD formulae of Eq. 10, where the hypotheses of the CEM are translated into direct relations between CO and CS transition probabilities, the same should happen in Eq. 8 where all the colour configurations are summed over and then considered on the same footage. In a process where the CS configurations dominate, such as $q\bar{q} \rightarrow \gamma ^\star \rightarrow Q \bar{Q}$, CSM and CEM predictions necessarily differ. Contrary to NRQCD, which encompasses the CSM, the CEM does not encompass the CSM. If one agrees with the data, the other cannot. The matter is then how precise the predictions and the data are to rule out one approach or the other.

Overall, one has to acknowledge that the conventional CEM central curves – as simplistic as the underlying idea of the model can be – give an account (Fig. 4a, b) of the world data points as satisfactory as the LO CSM. The latter seems to underestimate the data at low energies, while the former only has trouble to account for the TeV $J/\psi$ points; the slope being more problematic than the normalisation which can be adjusted. All this is qualitative, since the theoretical uncertainties on the CEM are as large as that on the open heavy-flavour production, which are admittedly large (see [27] for an up-to-date discussion of the $c\bar{c}$ case).

7 Energy dependence of the colour-singlet channel at Born and one-loop accuracy

As we just stated, the LO CSM curves provide a surprisingly good description of the $J/\psi$ and $\Upsilon$ data at high energies without adjusting – and even less twisting – any parameters. Although the central LO CSM curves agree with the data, the conventional theoretical uncertainties – from the arbitrary scales and the heavy-quark mass – are large (see the lower panels of Fig. 2). It is therefore very natural to look whether these uncertainties are reduced at one-loop accuracy. Such an observation was already made for the $\Upsilon$ case in [28], but this study was limited to a single $\sqrt{s}$, i.e. 200 GeV.

7.1 Spin-triplet quarkonia: $J/\psi$ and $\Upsilon$

Contrary the CO channels, the one-loop corrections to the CS channels only arise at $\alpha _S^4$ (see e.g. Fig. 1j, k). Nevertheless, these are known since 2007 [75] and can also be computed with FDC as done in [76]. In particular, there is no specific difficulty to integrate the $\alpha _S^3$ and $\alpha _S^4$ contributions in $P_T$, since they are finite at $P_T=0$.

However, as already noted in [29], such an NLO result tends to show negative values at low $P_T$, which can have a non-negligible impact on the total (i.e. $P_T$-integrated) cross section. To our knowledge, the energy dependence of the CSM at NLO has never been studied in detail. This is done below.

Figure 5 shows the energy dependence of the NLO CSM (7 curves). Note that, if a curve is not shown until 14 TeV, this indicates that the total yield got negative. The three red curves correspond to the default scale choices ($\mu _R=\mu _F=2 m_Q$) and are indicative of the heavy-quark mass uncertainty, on the order of a factor of 4 for the $J/\psi$ and 2.5 for the $\Upsilon$. In the former case, all three curves end up negative somewhere between 500 GeV and 2 TeV. Note also that the upper curve at low energies, i.e. for $m_c=1.4$ GeV, is the first to get negative and crosses the other ones as if the negative contribution were more important for lighter systems.¹³ In the $\Upsilon$ case, these three curves do not become negative at high energies – we have checked it up to $\sqrt{s}=100$ TeV. Nonetheless, they start to significantly differ from the LO curves (three blue curves) above 1 TeV, contrary to the good LO vs. NLO convergence found at RHIC energies in [28]. One might thus be tempted to identify this weird energy behaviour with a low- x effect.

Fig. 5

The cross section for direct a $J/\psi$ and b $\Upsilon (1S)$ as a function of cms energy in the CSM at LO and NLO for various choices of the mass and scales compared with the existing experimental measurements (see text)

Going further in the $J/\psi$ case, one can vary the factorisation and renormalisation scales about the default choice. Doing so, one obtains two classes of curves. For $\mu _R>\mu _F$ (pink and orange), the yield remains positive, but it is not less unphysical for it to be practically constant as the energy increases between 1 and 10 TeV! The only way to recover a semblance of increase is to take a large value of $\mu _R$ – and seemingly also a small value for $\mu _F$. Obviously, whatever the reason for this behaviour is, for large enough $\mu _R$, the QCD corrections which are proportional to $\alpha _s(\mu _R)$ necessarily get smaller and any difference between LO and NLO results should decrease. In the $\Upsilon$ case, as for the $J/\psi$ case, both curves with $\mu _R>\mu _F$ (pink and orange) correspond to the highest yields at high energies and the lowest at low energies. When one chooses $\mu _F>\mu _R$ (purple and green), the high-energy yields become negative both for $J/\psi$ and $\Upsilon$. In many respects, these observations are very similar to those made on the $^1S_0^{[8]}$ case. Such a pathological behaviour may thus be unrelated to the nature to the final state (see also next section).

Large NNLO corrections are expected to show up at large energies (low x) as discussed in [77]. It is not clear if they could provide a solution to this issue. Another way to solve this might be to resum initial-state radiation as done in the CEM [78] and for some CO channels [16].

In the light of such results, the most that one can reasonably say is that the NLO CSM results may be reliable for $\Upsilon$ up to 200 GeV and for $J/\psi$ up to 60 GeV, that is, up to $\sqrt{s}$ about 20 times the quarkonium mass. Above these values, the best that we have is the Born-order results.

7.2 Spin-singlet pseudo-scalar quarkonia: $\eta _c$ and $\eta _b$

Contrary to the spin-triplet case, one can obtain analytical formulae [22, 79] for the spin-singlet pseudo-scalar production cross section such as that of $\eta _c$ and $\eta _b$. This might in principle be of some help to understand the weird energy behaviour of the CS $^3S_1$ yield and of some CO channels. Indeed, the LO production occurs as for some CO channels without final-state-gluon radiation. In fact, the final state is simply colourless.

As can be seen on Fig. 6, the issue is similar in many respects but for the fact that one does not obtain negative yields for the $\eta _b$. For the $\eta _c$, the curves for $\mu _R>\mu _F$ remains positive at high energies – as for the $J/\psi$. One also sees that the crossing of the central LO and NLO curves occurs at larger $\sqrt{s}$ than for the $^3S_1$ states. However, such small quantitative differences may be due to the fact that we computed the y-integrated cross sections using the analytical expressions of [22] instead of the y-differential cross section at $y=0$.

Fig. 6

The cross section for direct a $\eta _c$ and b $\eta _b$ as a function of cms energy in the CSM at LO and NLO for various choices of the mass and scales

We have investigated this in more detail by looking at the different NLO contributions (the real emissions from gg and qg fusion as well as the virtual (loop) contributions) in order to see which channels induce the negative contributions and for which scale/mass values. However, it must be stressed that the decomposition between these different contributions depend on the regularisation method used. For instance, the decomposition is drastically different when using FDC – with sometimes a very large cancellation between the positive real-emission gg contribution and the negative sum of the Born and loop gg contributions – and the formulae of [22, 79] – where all the gg contributions are gathered. Yet, we checked that we obtain exactly the same results with both methods; the regularisation method or numerical instabilities cannot be the source of the issues observed above.

As regards the qg contribution, only $\mu _F/m_Q$ matters to tell whether it will change sign. For $\mu _F$ close to $m_Q$ and below, it will be positive (negative) at small (large) $\sqrt{s}$. Otherwise, it remains negative for any $\sqrt{s}$. The value of $\mu _R$ only influences the normalisation.

As to the gg contributions, which are expected to be dominant at high energies, both $\mu _F/m_Q$ and $\mu _R/\mu _F$ matter. For $\mu _F \simeq m_Q$, the gg contribution monotonously increases as a function of $\sqrt{s}$ irrespective of $\mu _R/\mu _F$. For $\mu _F \gtrsim m_Q$, the gq one is rather small and, despite being negative, does not induce a turn-over in the increase of the cross section. For $\mu _F \simeq 2 m_Q$, the gg contribution gets negative at large $\sqrt{s}$ for $\mu _R \le \mu _F$. For $\mu _R > \mu _F$, it remains always positive. Yet the sum $gg +gq$ can still become negative since, in some cases, gq increases faster with $\sqrt{s}$.

All this seems in line, for the gg [gq] part, with the formulae (C.25) and (C.26) [(C.32) and (C.35)] in Appendix C of [22]. Both contributions indeed exhibit the logarithms of $\mu _F/m_Q$ multiplied by a factor function of $m_Q/\hat{s}$.

As we can see, the results are already difficult to interpret for the spin-singlet case. For the spin triplet, we do not have similar analytical results. We can only guess that the structure is similar.

What seems surprising is that when one inspects similar expressions for the $\eta _Q$ production at NLO in the TMD factorisation [32], such negative terms do not appear as obvious. We are thus entitled to wonder whether such a formalism, which automatically resums the logarithm of the transverse momenta, may provide the solution to this issue. Another possible solution may be the consideration of NNLO corrections, which may show opposite signs to that at NLO. This is obviously beyond the scope of our analysis.

8 Conclusion

We have performed the first full analysis of the energy dependence of the quarkonium-production cross section at one-loop accuracy both in NRQCD and in the CSM. Taken at face value, our results would indicate a severe breakdown of NRQCD universality – in line with the previous analysis of Maltoni et al. [15] – unless one keeps the LDMEs close to the fit of Butenschoen and Kniehl, which, however, disagrees with the $J/\psi$ polarisation measurements and the $\eta _c$ cross sections.

The situation is, however, slightly more intricate, since we have uncovered a weird – sometimes unphysical – behaviour at large energies where one approaches the small-x regime where non-linear effects in the parton densities may be relevant. This certainly casts doubts on the numerical values we obtained at LHC energies, since collinear factorisation, on which we based our analysis, could break down.

Yet, up to a few hundred GeV, the energy dependence of the various octet channels at NLO seems well behaved and there is no reason to doubt this. In this region, the NLO yield prediction by NRQCD after fitting the mid- and high-$P_T$ quarkonium data – i.e. the yield and its polarisation – would overshoot the data by a factor ranging from 4 to more than 10. The same holds true at LO (see Appendix A). To reproduce the data, the CO LDMEs should be much smaller than what they are found to be in order to reproduce the Tevatron and LHC $P_T$-differential cross section in the case of the $J/\psi$ and $\psi (2S)$.

On the other hand, the LO CSM provides a reasonable energy dependence, in agreement with the existing data, except for $\sqrt{s} \le 40$ GeV, and is therefore up to ten times below the – data-overshooting – CO contributions. At one loop, the results are, however, ill-behaved for the charmonia and the cross sections can even become negative at large $\sqrt{s}$ for some – reasonable – scale choices. The situation is a bit better for $\Upsilon (1S)$. The same occurs for the spin-singlet quarkonia. In this case, one-loop results exist in the TMD factorisation approach and do not seem to be prone to such an issue. In the case of double $J/\psi$ production, the energy dependence at one loop seems well behaved [80]. Finally, we investigated the energy dependence of the yield in the CEM where the final states are treated rather differently and we did not find any specific problems.¹⁴ We are therefore tempted to attribute this problem to initial-state effects.

In [81], Ma and Venugopalan obtained a good description of the low-$P_T$ $J/\psi$ data over a wide range of energy by, on the one hand, using the LDMEs from [17] – our second set – and, on the other, a CGC-based computation of the low-$P_T$ dependence. In reproducing the data, they found that the CS contribution is only 10 % of the total yield. This 10 % is reminiscent of the factor 10 between the CS and CO in our “collinear” study. From our viewpoint, it looks as if the specific ingredient of this CGC-based computation would correspond to an effective reduction of the two-gluon flux¹⁵ by a factor of 10. It is therefore very interesting to find new processes which would be sensitive to this physics.

The negative yields obtained in the collinear case – observed for the $^1S_0^{[8]}$, $^1S_0^{[1]}$, $^3S_1^{[1]}$, $^3P_J^{[8]}$ channels – could also be cured by adding the large contributions of the one-loop amplitude squared – thus positive. This may look like an ad hoc solution which certainly questions the convergence of the perturbative series in $\alpha _s$. However, large NNLO corrections have already been discussed 10 years ago in [77]. Another path towards a solution may be higher-twist contributions [82], where two gluons come from a single proton as recently rediscussed in [83]. In any case, whatever the explanation for this situation may be, past claims that colour-octet transitions are dominantly responsible for low-$P_T$ quarkonium production were premature in the light of the results presented here.

Appendix A: LO NRQCD predictions

Sharma and Vitev recently performed [84] a LO fit of the CO LDMEs using RHIC, Tevatron and LHC $J/\psi$ data. Setting $\langle \mathcal{O}_{J/\psi }(^{1}\!S^{[8]}_{0})\rangle =\langle \mathcal{O}_{J/\psi }(^{3}\!P^{[8]}_{0})\rangle /m^2_c$ and accounting for the possible feed-downs, they obtained:

• $\langle \mathcal{O}_{J/\psi }(^{1}\!S^{[8]}_{0})\rangle =0.018$ GeV$^3$,

• $\langle \mathcal{O}_{J/\psi }(^{3}\!S^{[8]}_{1})\rangle = 0.0013$ GeV$^3$.

Based on this fit, we derive the energy dependence for the direct $J/\psi$, which is shown on Fig. 7. Unsurprisingly, the results badly overshoot the world data.

Fig. 7

The cross section for the production of a $J/\psi$ at LO from only colour-octet states as a function of the cms energy for various choice of the mass and scales