Curing the unphysical behaviour of NLO quarkonium production at the LHC and its relevance to constrain the gluon PDF at low scales

We address the unphysical energy dependence of quarkonium-hadroproduction cross sections at Next-to-Leading Order (NLO) in alpha_S which we attribute to an over-subtraction in the factorisation of the collinear singularities inside the PDFs in the MSbar scheme. Such over- or under-subtractions have a limited phenomenological relevance in most of the scattering processes in particle physics. On the contrary, it is particularly harmful for P_T-integrated charmonium hadroproduction which renders a wide class of NLO results essentially unusable. Indeed, in such processes, alphaS is not so small, the PDFs are not evolved much and can be rather flat for the corresponding momentum fractions and, finally, some process-dependent NLO pieces are either too small or too large. We propose a scale-fixing criterion which avoids such an over-subtraction. We demonstrate its efficiency for eta(c,b) but also for a fictitious light elementary scalar boson. Having provided stable NLO predictions for eta(c,b) P_T-integrated cross sections, sigma^NLO(eta(Q)), and discussed the options to study eta(b) hadroproduction, we argue that their measurement at the LHC can help better determine the gluon PDF at low scales and tell whether the local minimum in conventional NLO gluon PDFs around x=0.001 at scales below 2 GeV is physical or not.


Introduction
The production of quarkonia (Q) in inclusive proton-proton and electron-proton collisions when the protons break apart is one of most often studied process at high-energy colliders. Yet one still does not agree on how these heavy quark-antiquark bound states are produced. The interested reader will find it useful to consult the following reviews [1][2][3][4] addressing HERA and Tevatron results and more recent ones [5,6] as what regards the recent advances in the field with the RHIC and LHC. Besides probing QCD at the interplay between its perturbative and nonpertubative regimes, quarkonium production -once theoretically understood-should in principle allow us to probe the proton gluon content in terms of PDFs (see e.g. [7][8][9][10]) or TMDs (see e.g. [11][12][13][14][15][16][17][18][19][20][21]).
What we have learnt in the recent years with the advent of NLO computations of P T -differential cross sections, dσ NLO /dP T , of J/ψ and Υ [22][23][24][25][26][27][28] is that the inclusion of NLO corrections in any data-theory comparison is absolutely mandatory to extract qualitatively reliable statements. This is particularly true in two of the three most used approaches, the Colour-Singlet Model (CSM) [29][30][31] and Non-Relativistic QCD (NRQCD) [32], where P T -enhanced NLO contributions notably affect observables like dσ/dP T and the yield polarisation as a function of P T . In fact, for S -wave quarkonia, the CSM is the leading NRQCD contribution in the heavy-quark velocity, v. For the Colour-Evaporation Model [33,34], the impact of NLO corrections [35,36] to dσ/dP T is limited. In the latter model, all the spin and colour contributions of the heavy-quark pair are summed over and the possible additional gluon radiations at NLO do not open new production channels at variance with the CSM and NRQCD.
When one integrates over P T , these NLO channels, which are P T -enhanced in the CSM and NRQCD and which are precisely responsible for the large impact of the NLO corrections at mid and large P T , are just suppressed by one power of α s without any P T -enhancement factor. In this context, in 2015, we studied [37] the energy dependence of the J/ψ and Υ P T -integrated cross section at NLO, σ NLO , in NRQCD to verify the coherence with the NRQCD predictions for dσ NLO /dP T . Beside the confirmation of a possible breakdown of NRQCD universality 1 -as first claimed [38] by F. Maltoni et al. based on a partial NLO study-, we found out that, for all the NRQCD contributions 2 , the energy dependence became unphysical once the α s corrections were added. The same observation was made for the η c which is at the centre of this study. More precisely, the charmonium cross sections would become negative at increasing energies for a wide class of factorisation and renormalisation scales. In the J/ψ case, the NLO corrections significantly reduce the predicted yields close to RHIC energies [39] and σ NLO J/ψ already becomes negative at a couple of hundred GeV at central rapidities, y, for µ F ≥ M ψ [37]. Such observations were already made in the 1990's regarding the η c independently by Schuler [40] and then by Mangano & Petrelli [41] but were then essentially forgotten, see e.g. [42].
For bottomonia and for some -small-µ F scale choices for charmonia (see [37] for details), σ NLO would not become negative but the K NLO factor, defined as σ NLO /σ LO , would steadily deviate from unity for increasing √ s. As discussed above, large K NLO factors have already been observed in quarkonium production at finite P T but they can then be explained by kinematical factors scaling like P T /m Q . These are absent when P T is integrated over. As we noted such an intriguing behaviour can already be observed for √ s on 200 ∼ 300 GeV [37], so not necessarily at very high energies where large logarithms of the colliding gluon momentum fraction x should be accounted for. Indeed, such energies typically corresponds to x = 0.01 and even above.
In this article, we propose a solution to this issue which we attribute to an over-subtraction in the factorisation of the collinear singularities inside the PDF in the MS scheme. As such, it may appear in any NLO computations once a couple of unfavourable factors combine. In general, such over-subtractions indeed have a limited phenomenological relevance. It is clearly not the case for charmonium production which therefore offers a neat study case. As we will discuss, we propose a simple solution which consists in a factorisation-scale choice based on the high-energy limit of the partonic cross section and we demonstrate how well it works for η c and η b production and for the production of a fictitious elementary light boson, whose production mechanism is at odds with the production of a non-relativistic pair of heavy quarks then forming a quarkonium.
Having proposed a way to get sound NLO perturbative results, we discuss the interplay between the behaviour of the gluon PDFs at low scales and, in particular, the η c production cross sections. This motivates us to encourage a vigorous experimental effort to measure it and, before data are available to fit them, we are tempted to suggest experts in PDF fits to analyse how degraded global fits would be if xg(x; µ F ) at NLO is required to be monotonous for x < 0.01 at µ F ∼ m c .
The structure of the article is as follow. In section 2, we outline the structure of the NLO η Q production cross sections and explain how to reproduce the existing results. On the way, we provide analytical expressions in terms of the partonic cross section and of the partonic luminosities needed to compute the rapidity differential cross section, dσ NLO /dy, which are not available in the literature. In section 3, we make a brief historical survey of the past phenomenology of NLO η Q hadroproduction and of the attempts to identify the origin of these negative NLO cross sections and we explain that they come from the subtraction procedure in the factorisation of the collinear singularities in the MS scheme. Section 4 is devoted to our factorisationscale choice. Section 5 gathers our resulting cross sections for η c and η b . We first demonstrate that our proposal works by discussing the behaviour of the K NLO factors for η Q and elementary scalar bosons. Then, we discuss the interplay between the gluon luminosity and our obtained cross sections and finally we present what we believe to be the best possible NLO predictions. Section 6 gathers our conclusions and an outlook at other quarkonium-production processes.
2 η Q production up to NLO in the collinear and NRQCD factorisations

CSM, NRQCD and collinear factorisation
The present study essentially bears on collinear factorisation [43] whereby the hadronic cross-section to produce a quarkonium Q is factorised into a convolution of PDFs and a partonic cross section,σ(Q). Through NRQCD factorisation [32], the latter is further factorised into short-distance perturbative parts, computable with Feynman graphs, and long distance non-perturbative parts. As a result, one starts for the production of a quarkonium Q in a collision of two hadrons A and B from: where e.g. f a/A is the PDF of the parton a inside the hadron A, dσ ab (QQ [n] + {k}) are proportional to the partonic differential cross-section to produce a QQ pair in the (spin and colour) quantum number n, possibly with other particles {k} from the scattering of the partons ab and O n Q ) is an NRQCD Long-Distance Matrix Element (LDME) for the non-perturbative hadronisation of the pair in the state n into the quarkonium Q. NRQCD factorisation stems from an expansion in the relative velocity v between the QQ pair in the quarkonium rest frame. In this work, we focus on the terms leading in v and sub-leading in α s . As such, we only need to consider the colour singlet 1 S [1] 0 state for pseudo-scalar quarkonia, which is thus equivalent to the CSM. In such a case, the sum over n in dσ ab (Q) reduces to a single term. The purpose of the next sections is to explain how dσ ab (Q) can computed up to NLO accuracy in order to explain the appearance of negative cross sections in past computations.

η Q hadroproduction at LO
At LO (α 2 s ), η Q hadroproduction proceeds through gluon fusion, g(k 1 ) + g(k 2 ) → η Q (P), which can be computed via Feynman diagrams like Fig. 1a. In the CSM [29][30][31], the matrix element to create a 1 S 0 pseudoscalar quarkonium η Q with a momentum P, possibly accompanied by other partons, noted {k}, is obtained from the product of the amplitude to create the corresponding heavy-quark pair, M(ab → QQ + {k}), a spin projector, N(Ps.|s 1 , s 2 ) and R(0), the η Q radial wave function at the origin in the configuration space. The CSM being the leading v contribution to NRQCD, R(0) can naturally be related to a NRQCD LDME as follows: By virtue of heavy-quark spin symmetry, R(0) is identical for the η c and J/ψ for instance up to v 2 corrections. It can then be obtained from the well measured leptonic width of the J/ψ computed in the CSM/NRQCD or from potential models. In what follows, we will use |R η c (0)| 2 = 1 GeV 3 and |R η b (0)| 2 = 7.5 GeV 3 [39].
Overall, one has where P = p Q + pQ, p = (p Q − pQ)/2, s 1 and s 2 are the heavy-quark spins, and δ ii / √ N c is the projector onto a CS state. For v → 0, the spin projector on a pseudoscalar state, After one sums over the quark spin, one obtains traces which can be evaluated in a standard way.
However we note here the explicit appearance of γ 5 in N(Ps.|s i , s j ) which can cause issues within the framework of dimensional regularisation. Here we employ the standard 't Hooft-Veltman scheme to deal with γ 5 in D = 4 − 2dimensions [44]. We obtain for the LO matrix element squared [45], and for the partonic cross section (where we have set M Q = 2m Q as expected within NRQCD and N c = 3), where we have defined z = M 2 Q /ŝ. The hadronic section then reads with τ 0 = M 2 Q /s = 4m 2 Q /s and τ = τ 0 /z, in terms of the following differential gluon luminosities: These fully encapsulate the energy and rapidity dependences of the η Q yields at LO.

η Q hadroproduction at NLO
Let us now outline how to compute the η Q cross section up to NLO accuracy [40,45,46] which we will then use throughout our study. NLO contributions involve both virtual(-emission) and real(-emission) corrections via gg fusion that can be represented by diagrams in Fig. 1b, Fig. 1c and Fig. 1d. In addition to the gg fusion, qg and qq channels contribute at order α 3 s as shown in Fig. 1e and Fig. 1f. Both real and virtual contributions individually exhibit singularities. In order to deal with these singularities, we employ dimensional regularisation where we define D = 4 − 2 .
As usual, the virtual contributions exhibit both Ultra-Violet (UV) and Infra-Red (IR) divergences. The former are removed via the renormalisation procedure. To do so, we apply the on-shell (OS) renormalisation scheme for the gluon/quark wave functions and the heavy-quark mass counter-term while, for the coupling, we perform the renormalisation δZ MS g within the MS-scheme and we take [47] δZ where β 0 = 11 3 C A − 4 3 T F n f with n f being the number of active light flavours. We have above made a distinction between UV and IR to label the poles coming from UV and IR divergences respectively. In the following we will only label the poles to show their UV/IR character but not the appearing in exponents. The with and without labels ultimately originate from the regulator in D = 4 − 2 . We have also absorbed a global factor of e − γ E (4π) inside the MSrenormalised α s coupling.
As what regards the virtual corrections, we are thus only left with soft IR divergences. In contrast to the virtual contributions where the singularities are already manifest in the η c -gluon form factor 3 , the divergences in the real-emission part only reveal themselves after taking the phase-space integration 4 .
For dσ/dy, the phase-space integration is slightly less straightforward to be performed analytically than forσ where one can just integrate over the full phase-space without separating out the rapidity y and the transverse momentum P T . 3 That is the contribution proportional to δ(1 − z). 4 z → 1 for soft andt,û → 0 for collinear divergences.
For bothσ and dσ/dy, after combining the virtual with the real corrections, the soft singularities vanish 5 and we are left, as usual, with the initial-state collinear divergences which originate from diagrams such as in Fig. 1d and Fig. 1e. The occurrence of these divergences is a consequence of the fact that the initial states are fixed by the kinematics and therefore not integrated over.
Under the collinear factorisation, the divergences arising from such collinear emissions from specific initial partons are subtracted in the factorised PDFs via the corresponding Altarelli-Parisi (AP) Counter-Terms (CTs) which introduce the factorisation scale µ F in the partonic cross section [48]. In the MS scheme, the AP CT forσ ag reads, where Γ [1 + ] is the Gamma function and P ag (z) are the splitting functions between parton a and a gluon. We have given their expressions in Appendix A.
Using these standard procedures, we have reproduced the expressions of the partonic cross sections up to α 3 s [40,45,46]. It does not generate any specific complications to fold these with PDFs. We have collected in Appendix B the final expressions for the integrated cross section σ NLO η Q in terms of the partonic luminosities. On the contrary, the analytical expressions needed to obtain dσ NLO η Q /dy are absent in the literature. We have gathered them in terms of the partonic luminosities for the three channels gg, qg and qq in Appendix C. The codes which we have derived from these expressions and which we have used to generate the results presented later have been successfully cross-checked versus the semi-automatic code FDC [49].
3 On the origin of unphysical η Q cross section at high energies

The NLO partonic cross section and its HE behaviour
In this section, we focus on the partonic high-energy (HE) limit (ŝ → ∞ or equally z → 0) and show how this limit can help us understand the origin of the unphysical cross-section results which we referred to in the introduction.
The first NLO computation for pseudo-scalar quarkonium production was done [45] by Kühn & Mirkes in 1992 for toponium. At the time, it was not known that a toponium state could not bind. Their NLO results were confirmed by G. Schuler [40] two years later who performed the first phenomenological application for charmonia. He was the first to report negative cross sections for η c production at √ s just above 1 TeV for the central scale choice. He explained this 5 But for a soft singularity proportional to β 0 that arises through renormalisation. This factor will be absorbed inside the PDFs, see later.
unphysical behaviour by the fact that the partonic gg cross section was approaching a negative constant for µ F = M Q at highŝ. When folded with PDFs, such negative contributions coming from real emissions would become larger than the Born contributions for too flat low-x gluon PDFs. However, as what regarded the reason why the gg-partonic cross section was approaching a negative constant at highŝ, he did not provide any explanation, only a suggestion of a possible side effect of the restriction in the heavy-quark kinematics for them to be at threshold to form a non-relativistic bound state like a quarkonium.
In 1996, while presenting preliminary NLO crosssection results within NRQCD, Mangano & Petrelli discussed in a proceedings contribution [41] similar issues; they then attributed these negative cross sections to a possible over-subtraction of the collinear divergences inside the PDFs, thus rendering the partonic cross section negative in the HE limit. Quoting them, "there is nothing wrong in principle with these [partonic] cross sections turning negative in the small-[z] region, as what is subtracted is partly returned to the gluon density via the evolution equations". They however also noted that, for processes like charmonium production occurring at scales near where the PDF evolution is initiated, it is insufficient in practice -hence the negative hadronic cross sections then observed by Schuler.
An important observation we would like to make here for our following reasoning is that the magnitude of these negative partonic cross sections for z → 0 is process dependent. As such, the universal PDF evolution, for a given scale, cannot thus possibly fix the issue in a global manner.
To further assess this, let us indeed focus on the small-z limit ofσ ab which we obtained in the previous section, as done by Schuler, Mangano and Petrelli. Taking this limit in Eq. (B.6) & (B.8) 6 , one gets where C g = 2C A , C q = C F and M Q is the mass of the produced quarkonium (or 2m Q ). Equivalent expressions were obtained for P-wave quarkonium [40,41]. One can also consider such a limit [50] for the production of the Brout-Englert-Higgs (BEH) scalar boson H 0 using NLO expressions [51][52][53] or a fictitious elementary scalar boson, dubbedH 0 , whose coupling to gluons also occurs through a loop of heavy quarks and get a similar limit. In all these cases, we stress, since it will be essential for our forthcoming discussion, that this limit [50] in fact yields A g = A q .
We further note the presence of the factorisation scale µ F inside log M 2 µ 2 F in these limits. This term is in fact universal and process-independent as it originates from the AP 6 The following discussion applies to bothσ and dσ dy .
CT (see Eq. (9)) to subtract the initial-state collinear divergences. On the other side, A a is clearly process-dependent as Table 1 illustrates it. As a consequence, if A a < 0, the HE limit of the partonic cross section thus gets negative for the natural scale choice µ F = M Q and above. Whether this can make the hadronic cross section turn negative is then a matter of a complex interplay between the hadronic energy, the PDFs, µ R and the size of the (process-dependent) virtual corrections.
Before discussing this interplay, let us however go back to the notion of over-subtraction to trace back the origin of these negative limits. Away fromŝ = M 2 , only the real emissions contribute. In fact, at largeŝ, the solet-channel gluonexchange topologies depicted by When integrating overt, one will encounter the aforementioned collinear divergences which are to be absorbed in the PDF via the AP CT. Anticipating this subtraction, we can exhibit the corresponding divergence and recast the un-renormalised cross sectionσ as where above we have split the collinear part fromĀ a (z) 8 which is free of divergences for any 0 ≤ z < 1. We have multiplied the first term by a factor D a = 1 + δ ag to account for the fact that one has collinear singularities for each gluon in the gg channel. Therefore one would need to take 2σ (AP-CT) gg , i.e. for each parton, to eliminate the poles. From the equation above, it follows thatσ NLO,z 1 ag is positive-definite due to the fact that the first term evaluates to positive infinite 9 as Clearly, other schemes to absorb these collinear divergences inside the PDFs would yield differentĀ a (z) 10 . In the DIS scheme for instance,Ā a (z) [45] exhibits a log z dependence, which does not create any issue once integrated over z and this different z dependence should in principle be compensated by a different evolution of the PDFs. Yet,Ā a (z 0) should remain finite.
What we wish to argue here is that, since this subtraction is the only possible source of negative numbers at z 1, if A a (z) happens to be negative in a given scheme where PDFs are supposedly positive (see [54] for MS), this signals that the AP CT have likely over-subtracted some collinear contributions from the real-emission contributions, and this can yield the observed negative hadronic cross sections. This is indeed what happens for quarkonia since the NLO threshold contributions (ŝ = M 2 ) are found to be positive-definite for η Q and several other states at least for µ F = µ R 11 . Note that σ NLO η Q also goes negative at large √ s for µ F = µ R .
Let us re-iterate at this stage that contributions of type full square |M| 2 like the real emissions are always positivedefinite by construction at any kinematical point z. The only way to render them negative is the over-subtraction via the AP-CT inside the PDFs. We agree that evolved PDFs can 8 We remark at this stage that if the form factor of the Born cross section is resolved, i.e. considering the top-quark loop with a finite mass in the case of H 0 production via gluon fusion, A a =Ā a (z = 0) is a constant. On the contrary, i.e. the coupling is tree-level type-like as in the Higgs EFT with m t → ∞, then we have an additional log z dependence and a different off set for the gg and qg channel. It is not very surprising as for z → 0, m t andŝ are both large and HEFT cannot be applied. 9 For IR poles, one has that IR < 0, while for UV poles UV > 0. 10 In principle, one could thus look for a scheme where the partonic cross sections simply do not become negative. This is left for future investigations as it would entail refitting the PDFs with different evolution equations. 11 This is also the case for H 0 andH 0 . reduce the weight of these regions in z where the partonic cross sections are negative, and eventually avoid negative hadronic cross sections. Yet, it is hard to believe that they would do so for all possible processes where this can occur as the coefficients A a are process-dependent while the DGLAP evolution is process-independent.
3.2 From negative partonic cross sections to negative (or positive) hadronic cross sections Having now identified the origin of the negative cross sections, we can discuss their relevance to the past phenomenology which we recalled in the previous subsection.
First, we note that the η b phenomenology, for which σ NLO η b remains positive in the LHC range, is less pathological. We have indeed found out [37] that σ NLO η b only slightly deviates from σ LO η b in the LHC range. It thus seems that it is less sensitive to the limit of Eq. (10). Both charmonia and bottomonia have the same partonic cross section but for three changes: the mass shift and a trivial rescaling of the LDME and n f which plays a minor role here. This mass shift however has three immediate effects : (i) a given z = M 2 Q /ŝ value for bottomonia corresponds to 3 times larger √ŝ . Considering the rescaling on the integration bounds, [M 2 Q /s, 1], when convoluted with PDFs, this effectively corresponds to a 3 times larger √ s, which is thus easily outside the range of past studies, (ii) however, even at fixed z, the results would differ since α s (µ R M Q ) is smaller and this reduces the impact of α 3 s contributions compared to the (positive) Born ones at α 2 s , (iii) the evolved gluon PDFs up to a larger µ F become steeper which reduces the relative importance of the small-z domain compared to the threshold contribution at z = 1 which remains positive. Taken together, these 3 points explain very well why the charmonium case, at low scales, is the most pathological one and that the issue of a possible over-subtraction of the collinear divergences is usually considered to be rather academical with a limited impact on other hadronic cross sections.
It is however legitimate to wonder if further aspects specific to the modelling of quarkonium production renders its phenomenology particular. Our answer tends to be negative. Indeed, as Table 1 shows, CO and CS states are equally affected, in agreement with the past phenomenology [37]. This confirms that neglecting such higher-order v corrections is not the source the issue. Since CS and CO are both computed in the non-relativistic limit, one may wonder whether that this limit is also a source of issues, as suggested by Schuler. Yet, we anticipate that the phenomenology of ã H 0 with m Q = MH0 /2 should also be affected since A a is also negative. As our numerical results will show, it is indeed the case and this will thus confirm that this is not a quarkonium issue per se.

A scale choice as solution
We now come to our proposal to solve this unnatural behaviour of the cross section. In fact, it simply amounts to set the factorisation scale µ F such that partonic cross-section vanish at largeŝ, instead of risking it to become negative. Of course, such a scale choice is only possible provided that it is the same for all the partonic channels. Dubbing our scale choiceμ F , we just define it aŝ having in mind that that A q = A g . It is clear, from our definition that, since A a is a process-dependent quantity,μ F will be process-dependent. We have listed some values of µ F in Table 1 for the different particles we considered. It is important to note that theμ F values we have found are within or close to the usual ranges of values anyway taken in phenomenological studies.
Let us now turn to the physical picture of our reasoning. Our motivation is clear as it amounts to avoid negative cross sections which we attribute to an over-subtraction in the MS factorisation scheme. Our scale choice avoids thatσ ab be negative at small z. It makes sense to base its construction from this limit as it becomes more and more relevant at large s, precisely where σ NLO η Q can become negative. Even ifĀ a (z) becomes more negative than its limiting value, A a , our results will show that cancellingσ ab (z → 0) with µ F =μ F will be sufficient to get much more sound results, in particular to avoid σ NLO η Q < 0. Going further, we stress thatĀ g (z) also contains real emissions from the heavyquark line (see Fig. 1c) and thus differs fromĀ q (z). Working at finite z whereĀ g (z) Ā q (z) does not allow us to derive an equally simple gauge-invariant solution based on a scale choice. In the quarkonium case, the latter contributions tô σ gg are relatively suppressed by M 2 Q /ŝ = z and thus their effect disappear at small z.
Another reason to focus on the small-z limit is that when foldingσ ab (z) with the PDFs (see Eq. (1) or Eq. (B.6)), the Jacobian to transform the integration measure from dx 1 dx 2 to a measure involving dz will comprise a multiplicative factor 1/z 2 . As a result, the impact of the small-z region certainly depends much on whetherσ ab (z → 0) is zero or not, even though the z range has a lower bound set by τ 0 = M 2 /s. Indeed, for nonzeroσ ab (z → 0), the PDFs are the key element regulating the integral. By virtue of evolution, they should become steep enough as to essentially damp down the contribution of the small-z region. However, at low µ F , the PDFs can be rather flat. This can give a large weight to this small-z region where the real-emission contributions are negative for large µ F /M, hence the possibility that σ NLO η Q < 0. Now, ifσ ab (z → 0) = 0 as our µ F choice entails, the PDF shape at low scales is suddenly much less crucial to damp down a contribution which should not be the leading one in any case.
Physics wise, our scale choice essentially amounts, in the partonic HE limit, to reshuffle the entirety of the real emissions inside the PDFs 12 . From a HE viewpoint, such contributions are expected to be important, in particular at smallt since they are supposed to be enhanced by logarithms ofŝ, which should eventually be resummed. In fact, as our discussion has illustrated, the prominent effect of such contributions is a source of issues in a fixed-order computation as it jeopardises the convergence of the perturbative series with NLO contributions being more important than the Born ones. In this sense, our scale setting amounts to include all these possible HE effects in the PDFs. This makes sense as the PDFs are ultimately determined by fitting data -containing all type of higher order corrections. In fact, recent PDF analyses have been made taking into account HE effects in their evolution [56,57].

A word on our PDF choice
Given the importance of the PDF shape at low scales in the previous discussions, we have employed on purpose, thanks to LHAPDF6 [58], 3 NLO sets which show rather different features: : 1. a representative 13 set of the conventional NLO PDFs, PDF4LHC15 nlo 30 [59], 2. a dynamical PDF set, JR14NLO08VF [60], where gluons are radiatively generated from a valence-like positive input distributions at a low scale which is optimally chosen, and 3. a set taking into account HE effects in the evolution, NNPDF31sx nlonllx as 0118 [56], in order to perform our NLO cross-section evaluations.
These are plotted on Fig. 2 along with CT14nlo [61], MMHT14nlo [62], NNPDF31 nlo as 0118 [63] for comparison, for two scale choices, 1.55 GeV and 3 GeV. We note that x g(x) from PDF4LHC15 nlo 30, MMHT14nlo, CT14nlo and NNPDF31 nlo as 0118 all show a maximum around 0.02 and then a local mininum below 0.001. Such features are absent in both JR14NLO08VF and NNPDF31sx nlonllx as 0118 whereas they have a significant impact on the phenomenology as we will show later on. However, we stress that the local minimum has already disappeared once the gluon PDFs are evolved up to 3 GeV, where the 3 sets we have used display similar features but for the size of the uncertainties. At µ F = 1.55 GeV, we note that both for PDF4LHC15 nlo 30 and NNPDF31sx nlonllx as 0118, the shape can be very different within the uncertainty spanned by their PDF eigensets. At µ F = 3 GeV, this only remains the case for PDF4LHC15 nlo 30. These different behaviours will in fact be very useful to study the interplay between the scale and the PDF choices.
We further note that a recent study by Flett et al. has shown that one could extract, under specific assumptions [64], actual constraints on the gluon PDF at low scales from J/ψ exclusive production data. These are represented by the solid black lines in Fig. 2 (a) when applied to NNPDF30 nlo. Under such assumptions, the presence of a local minimum below 0.001 is unlikely as it would yield to a decrease in the J/ψ exclusive production cross section, which is absent in the data. At this stage, since it is not an actual PDF fit and since the exclusive cross sections are not directly related to the PDFs, we consider this finding as a guidance, yet a very interesting one anticipating our results. 5.2 Assessing the perturbative convergence withμ F using the K NLO factors We have found so far that significant NLO contributions to η Q production are expected to appear if the hadronic cross section becomes sensitive to z = M 2 Q /ŝ values far away from threshold. This would result in a significant √ s dependence of the NLO/LO hadronic cross-section ratio (K NLO ). As we explained, it is due to an over subtraction, in the MS scheme, of collinear contributions from the real-emission NLO contributions inside the PDFs. A relative constant offset is however expected from the virtual corrections at z = 1, like for the decay widths, in particular for reactions where α s is not very small.
To mitigate this fixed-order treatment shortcoming, we have thus proposed a specific scale choice which corresponds to the inclusion, in the z → 0 limit, of the entirety of such NLO contributions in the PDFs. The logic behind is that PDFs are fit to data which incorporate all such emissions. This is probably not a perfect solution but, beside of corresponding to perfectly acceptable µ F values, it indeed avoids erratically varying K NLO factors and negative and unphysical hadronic cross section, as the results of this section show.
Before discussing our results for K NLO , let us describe our set-up. We have evaluated them at y = 0 using Eq. (6) for LO and Eq. (C.10), (C.11) & (C.12) for NLO. The same PDF have been used for both. We have set m c = 1.5 GeV for the η c and m b = 4.75 GeV for η b . We have used the α s corresponding to our PDF choice thanks to LHAPDF6.
As for a fictitiousH 0 , we have set its mass at 3 GeV, close to that of η c . Having at our disposal, the small-z limit for different MH 0 /m Q ratio, we have chosen three values for the mass of the heavy-quark active in the loop, namely 0.5×MH 0 , MH 0 and (m t /m H 0 )×MH 0 . As can be seen from Table 1, m Q = 0.5 × MH 0 renders A a slightly negative, −0.147, whereas it is large and positive, 2.28, for the SM H 0 with m H = 125 GeV and m t = 173 GeV. The plotted K NLO factors have been computed with the publicly available code ggHiggs by Bonvini [67][68][69] based on [50,70] for the NLO result with a finite heavy-quark mass in the loop. Other than this, we have run with its default setup.
Let us first discuss the η Q results. Fig. 3 gathers our result for the K NLO factor computed at y = 0 for η c (top) and η b  (down) and for the central eigenset of our 3 NLO PDF sets, namely PDF4LHC15 nlo 30 (left), JR14NLO08VF (middle) and NNPDF31sx nlonllx as 0118 (right). We have used the conventional 7-point scale-choice values obtained by independently varying µ R and µ F by a factor of 2 about a default value which we simply chose here to be the mass of the quarkonium, M Q . We stress that LO cross sections used to compute K NLO were obtained with the same PDF and scale as those used for the NLO cross sections. In addition, we have plotted K NLO for µ F =μ F which we expect to provide the best behaviour. We have only plotted it for µ R = µ F .
We now discuss the qualitative features of the results. First, we note that, for PDF4LHC15 nlo 30, negative cross sections (K NLO < 0) appear as expected as early as 1 TeV. This happens first for µ F = M Q and µ R = 0.5M Q , then for µ F = 2M Q and µ R = M Q , and then for µ F = 2M Q and µ R = 2M Q , while K NLO essentially converges to 0 for µ F = 2M Q and µ R = M Q , which is also not acceptable. In short, all the results with µ F equal or larger than the default choice are pathological and the situation is worsened by a lower value of µ R which comes along with a larger size of α s . On the other hand, for µ F = 0.5M Q , K NLO does not get negative, neither particularly small, but shows a peak at the top LHC energies which is related to the peak in the low scale gluon distribution at low scales encoded in PDF4LHC15 nlo 30. So far, these results confirm the aforementioned past phenomenology.
On the other hand, when adopting our scale choice, µ F =μ F (= µ R ), the behaviour is smooth and seems to slowly converge towards a constant value slightly above the unity. In fact, we have checked that varying µ R at fixed µ F =μ F would simply shift the curve without changing its shape. This is connected to the fact that the limiting values is also driven by the threshold contribution at z = 1 where the virtual corrections which are sensitive on µ R sit. Now, summarising our result for a conventional PDF like PDF4LHC15 nlo 30, we can claim that the scale choice which we advocate provides a very simple solution to avoid pathological behaviour of the P T -integrated η c cross sections at NLO.
All the above observations can be made again for η b (see Fig. 3d) but for the fact that the K NLO factor does not get negative. Nonetheless, it gets so small for the large scale choices that the results remain meaningless. Presumably at √ s above those of a FCC, the cross section for µ F = 2M Q and µ R = M Q would turn negative. However, at such √ s, we admit that it is rather an academic example. Yet, we stress that K NLO varying by a factor of 10 from fixed-target en-    ergies to FCC energies is the sign of a bad convergence of the NLO computation for such scales. Besides, we do not observe any more a peak µ F = 0.5M Q , for which the energy dependence of K NLO starts to be acceptable. Choosing µ F =μ F gives the best trend with a quasi constant value, close to 1, for 1 TeV and above. Such a choice completely stabilises the K NLO energy dependence as it results that going to higher energies does not give an artificial importance to the α 3 s corrections. From the early studies of Schuler, Mangano and Petrelli one expects a strong sensitivity of the PDF shape on the impact of the NLO corrections (see also [42]). We have checked that the PDF uncertainty on K NLO derived from the 30 PDF4LHC15 nlo 30 eigensets is indeed smaller for µ F =μ F than for larger scales, despite the fact that the PDF uncertainty themselves usually decrease for growing scales. Actually to assess the PDF-shape sensitivity, it can more insightful to compare the trend with the central set of JR14NLO08VF and NNPDF31sx nlonllx which show a clear different shape in particular close to 1.5 GeV (see Fig. 2). These are respectively shown on Fig. 3b & 3e and Fig. 3c & 3f. For η c , the trend is very similar compared to what we obtained with PDF4LHC15 nlo 30 except for the absence of the peak for µ F = 0.5M Q . As we wrote above, such a peak resulted from the local maximum and minimum in the central PDF4LHC15 nlo 30 eigenset 14 .
14 Two effects can come into play here. First, the average momentum fraction of the gluons in the NLO contributions is slightly larger than for the LO one. As such, if the gluon PDF oscillates, it could happen that the PDFs product multiplying the gg NLO partonic cross section could be larger than the LO one. Second, as we previously discussed, In conclusion, even with a priori the steepest possible gluon PDFs compatible with a global NLO PDF analysis, one gets negative or strongly suppressed NLO η c cross sections for a majority of the conventional scale choices (5 out of 7), whereas that obtained with our scale choice µ F =μ F is remarkably stable. For η b , the 3 PDFs essentially yield the same K NLO factors which also shows the most stable behaviour for µ F =μ F .
Having demonstrated the efficiency of our scale choice to avoid anomalously large NLO corrections to pseudoscalar quarkonium production attributed to an over-subtraction of the collinear divergence inside the PDFs, let us now investigate whether it works for elementary scalar bosons coupling to gluons via heavy quarks. If our argumentation is correct, K NLO should, first, be rather µ F -and √ s-dependent for a scalar boson of similar mass than the η c and, second, become stable for µ F =μ F . Our results, shown on Fig. 4a, exactly confirms our expectation forH 0 with MH0 = 3 GeV and m Q = 1.5 GeV. If we were to work at even larger √ s, the larger µ F choices would eventually yield very small K NLO . They would probably not become negative but we recall that |A a | is smaller for H 0 than for η c rendering the HE limit slightly less harmful. On the other hand, too small scales yield strongly growing with flatter PDFs, σ is in principle more sensitive to the largeŝ behaviour ofσ gg . For the considered scale, µ F = 0.5M Q , this limit is positive, thus K NLO is expected to get larger, precisely right after the bump in the PDF luminosity. It is likely that the latter effect actually dominates. Indeed, for µ F =μ F , the limiting value ofσ is set to 0 while the bump in K NLO has nearly disappeared although there is still a slight bump in the PDF.
K NLO at large √ s. Finally, setting µ F =μ F , or close to it with µ F = MH0 sinceμ F = 0.93MH0 , gives remarkably stable K NLO . We consider this to be a confirmation that the instabilities in NLO computations of quarkonium production are not connected to the modelling of quarkonium production. The situation is equally good if we set consider, as an academic example, M H 0 = 125 GeV and m Q = MH0/2. Indeed, µ F =μ F yields the most stable results on Fig. 4d. On the other side, for the real H 0 case with m t = 173 GeV -for which the situation is of course not problematic-, we observe on Fig. 4f that K NLO tends to clearly increase with √ s for the smaller scale choices, like 3M H 0 yield to more stable trends. In principle, we would expect µ F =μ F to yield the most stable curve. More investigations are needed to understand why µ F = 2M H 0 shows the best trend. One observes the same for MH0 = 3 GeV and the same m Q /MH0 ratio on Fig. 4c. Fig. 4b and Fig. 4e, in between both cases, hint at an effect which would scale like m Q . A possible explanation could come from contributions of the box diagrams which would yieldĀ g (z) Ā g (z) down to very low z.
Yet, the phenomenology of H 0 being now made at N 3 LO accuracy [71,72] in the infinite top-quark-mass limit, it is rather an academical question. For the validation of our scale proposal, the success of the lightH 0 case with MH0 = 3 GeV and m Q = 1.5 GeV, with strongly decreasing K NLO for large µ F and a very stable one for µ F =μ F is much more telling since it perfectly confirms what we observed with the η c .

A word on gluon luminosities at NLO and low scales
Before moving to our NLO predictions for the cross sections, we find it useful to make a short digression on the gluon luminosities. Indeed, we would not want the reader to be confused by some unnatural cross-section behaviours as a function of √ s or y which we will show and to attribute these behaviours to QCD corrections to the hard scattering. Contrary to K NLO where the PDF impact is indirect 15 because of a large cancellation of their effects in the ratio, the gg contribution to the hadronic cross section will essentially be proportional to the square of gluon PDFs at low scales. As we have seen when discussing Fig. 2, the conventional PDF sets typically exhibit a local minimum below 0.001 at scales below 2 GeV. In this region, however, the gluon PDFs are only poorly constrained by scarce data sensitive to gluons at smaller x and larger scales.
To illustrate the typical effects that such a shape can induce on both the √ s or y distribution on a low-scale process, we have plotted on Fig. 5 and Fig. 6 the corresponding differential gluon luminosity which would normally multiply a simple gg fusion process at LO. We have done so for our 3 chosen PDF sets, for 2 masses M (3 and 9.5 GeV) and 3 corresponding scales µ F (0.5M, M and 2M) as a function of √ s and y in the ranges which will correspond to the NLO cross section plots which we will show in the next sections. It clearly appears that for PDF4LHC15 nlo 30, which is representative of usual PDF sets, and to a lesser extent for NNPDF31sx nlonllx as 0118, both the √ s dependence at y = 0 and the y dependence at √ s = 14 TeV will strongly be distorted for µ F = 1.5 GeV. Not only one observes a strong scale sensitivity in the luminosity magnitude 16 , but the distributions are very different. More importantly, it is very improbable that any measured differential cross sections, even at low scales, would follow the trend of Fig. 5f (and Fig. 5d) with a yield showing a global maximum around y = 5 at the LHC or a differential cross section at y = 0 essentially constant between the Tevatron and the top LHC energy like on Fig. 5a. As we will see, the corresponding NLO cross sections will be driven by this behaviour of the gluon luminosity.

A digression on the η b detectability
Having now at our disposal reliable NLO predictions for η b cross sections at hadron colliders, let us address the question of the feasibility of such studies and in particular of the extraction of cross sections. Whereas prompt η c production at the LHC has now been the object of two experimental studies [73,74] by the LHCb collaboration, the prospects for η b production studies are however less clear. Compared to the η c , we do not know much on the η b properties which was only discovered in 2008 by BaBar [75] and whose mass is 9.4 GeV, while its 2S excitation, the η b (2S ), was discovered in 2012 by Belle [76] with a mass of 10.0 GeV. The η b in fact has so far been observed only in e + e − annihilations, by BaBar [75,77], CLEO [78] and Belle [76,79]. Most likely, future measurements will be carried out by Belle II [80]. Its width has been measured to be on the order of 10 MeV. We note the good agreement with the LO estimate for Γ(η b → gg) = 8α 2 The agreement is confirmed up to NNLO [81]. So far no measurement of any branching fractions are reported in the PDG [82] and it is thus not clear in which decay channel it could be measured in hadroproduction. In addition, one needs to know the branching value with an acceptable uncertainty to derive a cross section to μF=1.5GeV μF=3.0GeV μF=6.0GeV  test state-of-the-art computations which do not address the decay but only the production.
Different theoretical ideas have been pushed forward about the usable decay channels at the Tevatron and the LHC. For a long time, the decay into a pair of J/ψ potentially clearly signalled by 4-muon events was considered to be a discovery channel in the busy environment of pp collisions. Even though, from the beginning, physicists were aware that the branching fraction into this channel could be very small, we should stress here that the production cross sections for η b at colliders are not small at all, as they are comparable to those for Υ(nS ) which are routinely studied at the LHC. As such, small branching fractions could still yield observable rates. First estimations [83] reported 17 B(η b → J/ψ + J/ψ) = 7 × 10 −4±1 , using massrescaling arguments applied to η c → φφ. However, this estimate was then questioned and searches via the detection of 2 charmed mesons were suggested [84]. Actual computations based on NRQCD [85,86] later yielded a much smaller B(η b → J/ψ + J/ψ), as low as 5 × 10 −8 . It was however suggested that final-state interactions, beyond the effects included in NRQCD computations, could enhance the di-J/ψ decay width by up 2 orders of magnitude [87]. It is thus clear that until B(η b → J/ψ+J/ψ) is actually measured elsewhere, it could not provide a way to derive cross-section measurements. Yet, given the current intense activity in J/ψ + J/ψ studies with the observation [88] of a di-J/ψ resonance, a search for this decay channel at the LHC could still be an option as it may also be looked for at Belle II. Another channel subject of debates is that into 2 charmed mesons as the first branching-fraction estimate [84], on the order of a few per cent, was then also questioned [89]. Finally, let us mention the inclusive channel η b → J/ψX which however suffers from large uncertainties owing to possible CO contributions [90]. Other ideas could be found by looking at the many χ b decay channels which have been analyses so far, see [82].
An alternative to these hadronic decays is the exclusive radiative decay, η b → J/ψγ, whose branching is computed to be on the order of 2×10 −7 [85], which remains admittedly small. In principle, it suffers from smaller theory uncertainties compared to those above. An even simpler branching to predict is that of η b → γγ, whose partial width is in fact known up to two loops [91,92] in NRQCD 18 . In addition, some theoretical uncertainties cancel with those of the total width computed in NRQCD by assuming the dominance of the decay into two gluons, which is also known up to NNLO [81] (see above), yielding We thus believe that this channel should seriously be considered for a first extraction of the η b hadroproduction cross section even though channels involving J/ψ are probably 17 Note that this should then be multiplied by the square of B(J/ψ → + − ), i.e. 6 %. 18 Other approaches may sometimes give different results, but not by factors differing by more than 2 [93]. easier to deal with as what regards the combinatoric background. Tools like EtabFDC [94] should also definitely be helpful for experimental prospective studies.
As what regards the experimental setups where η b hadroproduction cross sections could be measured in pp collisions, let us cite the LHC in the collider and fixedtarget modes, in particular with the LHCb detector. The latter mode has been studied in details in [95][96][97][98] as what regards quarkonium production. Its nominal √ s for 7 TeV proton beams reaches 114.6 GeV. Another possibility is the SPD detector at the NICA facility up to √ s = 27 GeV [99].
We will provide predictions for these 3 setups.

Cross section predictions
We are now finally in position to show our results for the cross sections as a function of √ s and of y for selected values of √ s which correspond to experimental setups where we believe the challenging measurement of P T -integrated η Q yield could be performed in the future. These are the LHC at 14 TeV in the collider mode and in the fixed-target mode at 114.6 GeV, having in mind in particular the LHCb detector, and to the SPD experiment at NICA which could run up to 27 GeV. As our study mainly addresses the interplay between the size of the NLO corrections and the choice of scales, mainly that of µ F , we start by showing, dσ NLO η c /dy| y=0 as a function of √ s, for the 7-point scale choice (in black) andμ F (in green) for the central eigensets of our 3 PDF sets. These plots are the exact counterpart of the K NLO | y=0 plots of Fig. 3 which allowed us to assess the much better convergence of the NLO results at high energies when taking µ F =μ F . The essential difference here is that the PDF effects do not cancel.
As what regards the η c results shown on the first row of Fig. 7, we first note that 3 out of the 7 scale choices leads to negative cross section (see the incomplete curves) irrespective of the PDF choice. No matter what the PDF shape is, too "large" a value of µ F inevitably leads to unacceptable results. This observation is obviously in line with previous studies [37,40,41]. We now know that it is related to large negative contributions away from threshold due to a process-dependent oversubtraction of the collinear divergences which cannot be compensated by the PDF evolution, in particular at low scales where they are not much evolved. Now, as anticipated during the discussion of K NLO | y=0 plots, choosing µ F =μ F provides much more sound results. The results are particularly good up to the FCC energies with the JR14NLO08VF PDF as can be seen on Fig. 7b.
Yet, as we discussed in section 5.3, most of the conventional PDFs exhibit, at low scales, a local minimum for x around or below 0.001. This results in gluon luminosities  -which correspond the expected kinematical dependencies for a simple LO gg fusion-essentially constant in the TeV range. Without any surprise, this what we observe for theμ F curve using PDF4LHC15 nlo 30 (Fig. 7a) and to a lesser extent for NNPDF31sx nlonllx as 0118 (Fig. 7c). As expected, we observe the same with MMHT14nlo, CT14nlo and NNPDF31 nlo as 0118 on Fig. 8 both at LO and NLO.
The η c NLO energy dependence admittedly does not make sense when it remains a constant between √ s = 10 GeV and 10 TeV ! We urge the global fitters to examine whether global NLO fits cannot in fact be slightly amended in order to yield monotonous gluon PDFs at scale below 2 GeV. It is very important to realise that such a distorted shape is not at all due to the NLO corrections, but entirely due to the PDFs, as the JR14NLO08VF with monotonous gluon PDF case shows ( Fig. 7b and Fig. 8a).
Let us now turn to the η b case for which we know that the issue of over-subtraction is less problematic. Indeed, one only sees, on the second row of Fig. 7, a slight deviation at the FCC energies for the curves for µ F = 2M and µ F = 0.5M which admittedly is the most critical one according to our K NLO | y=0 analysis (see Fig. 3  K NLO η b | y=0 are significantly µ F -dependent, the dσ NLO η b /dy| y=0 are much less µ F -dependent as the difference in the - process-independent-gluon evolution induced by the different chosen µ F values efficiently compensates the explicit µ F dependence of the -process-dependent-hard scattering for η b . In such a case, the conventional scale choices would range between 5 and 20 GeV. Clearly, for η c , such a compensation does not work at NLO and the only solution to the issue remains our scale choice µ F =μ F . This is why for the following plots we stick to this scale choice which we consider to give, at the moment, the best possible NLO predictions for η c production in the TeV range. On Fig. 9, we show again dσ NLO η Q /dy| y=0 but with the PDF uncertainty associated with our PDF set choices for µ R =μ F (green band) and the µ R uncertainty for µ R ∈ [M/2 : 2M] (red band). The same observation as above can be made. Sinceμ F = 1.82 GeV for η c , PDF4LHC15 nlo 30 results show a large distortion due the PDF eigenset shapes and, clearly, nobody would expect to see it in any experimental data in the future 19 . The situation is better for NNPDF31sx nlonllx as 0118 and remains very good for JR14NLO08VF. For η b , the 3 PDF sets yield similar results. The smaller uncertainty band for JR14NLO08VF (Fig. 9e) simply comes from the very small PDF uncertainty 19 We stress that dσ NLO ηc /dy| y=0 for µ F =μ F only get negative on Fig. 9a and Fig. 9c because some gluon NLO PDF eigensets get negative. Negative cross sections would disappear at LO since the gluon PDF are squared at y = 0. of this set. It is admittedly much smaller than the conventional ones. Possible experimental data should be able to test such PDFs rather straightforwardly.
To go further in the analysis of our improved NLO results with µ F =μ F , we have plotted on Fig. 10 the relative uncertainties, dubbed as ∆σ/σ, from the PDF uncertainties and from µ R variations, obtained by normalising the upper and lower values of dσ η Q /dy| y=0 by the respective central ones, i.e. from the central eigenset for the PDF uncertainty and from µ R = M Q for the µ R uncertainty. The resulting (green and red) bands are compared to the LO µ R uncertainty (horizontal dotted blue line)-which is obviously a constant. Our first observation is that the renormalisation scale uncertainty is clearly reduced at NLO, which is a good sign of the α s convergence even at these low scales. Second we note that the scale uncertainty for η c is smaller, for √ s > 2 TeV, than that from PDF4LHC15 nlo 30, on the order of 30 % and then steadily growing; thich is representative of what NLO global fits would give. This means that forthcoming η c data at the LHC with a precision of 10 %, or lower, should already be enough to improve PDF fits, even taking into account the µ R uncertainty. For η b , both are of similar sizes and one should probably look at dσ NLO η Q /dy at fixed √ s or at various √ s if available to get more discriminating power on the PDF along the lines of [100][101][102]. The same holds for η c at lower energies, which should then be differential in y. At this point, we wish to stress that we have decided not to vary the charm and beauty quark masses 20 which are also usually considered to yield additional theoretical uncertainties. We however note that the induced variations are highly correlated in y (and √ s) and such correlations can expediently be used to alleviate their effect in a possible PDF fit.
In view of the above observations, we have decided to show dσ NLO η Q /dy at fixed √ s on Fig. 11 only for our µ F scale choice, only using JR14NLO08VF and its uncertainty (green band), along with the µ R uncertainty (red band) compared to the LO results (gray and beige bands). Our objective with these plots is to provide NLO predictions to motivate prospects for measurements and then the extraction of constraints on PDFs, rather than to test quarkonium-production models. From a kinematical viewpoint, dσ NLO η c /dy measurement should offer reliable constraints on the gluon PDF x dependence in the approximate range [5 × 10 −2 : 1] for SPD, [10 −2 : 1] for AFTER@LHCb and [10 −6 : 5 × 10 −2 ] for LHCb (at 14 TeV). At SPD, the η c production cross section is expected to be on the order of 1 pb. At AFTER@LHC, it would grow to 10 pb to reach 200 pb at the LHC. Of course, to these, one should apply rapidity-acceptance cuts, besides the appropriate branching fractions. As for the η b , the yintegrated cross section, σ NLO η b (s) is respectively expected to be 0.5 nb, 60 nb and 10 pb.

Conclusion and outlook
In this work, we have addressed the unphysical predictions of the collinear and NRQCD factorisations for the P Tintegrated quarkonium production, whereby negative cross sections are obtained for most of the conventional factorisation and renormalisation scale choices at LHC energies down to RHIC energies in some cases. In particular, we have focused on the pseudoscalar case, which is by far the simplest to tackle with analytical results available for the total cross section available since the mid nineties. On the way, we have provided analytical results for the rapidity differential cross sections, which were not available elsewhere.
We have shown that this unphysical behaviour can be explained through the high-energy limit of the partonic cross section, which is negative unless µ F is chosen to be relatively small compared to the quarkonium mass. This negative limit can be ultimately traced back to an oversubtraction of the Altarelli-Parisi counterterm in the MS scheme to absorb the collinear divergences inside the PDFs. This over-subtraction should usually be returned via the DGLAP evolution with steeper PDFs. However, the highenergy-limit values A a are process-dependent while the DGLAP evolution is clearly process-independent at fixed 20  scales. This over-subtraction cannot be returned in a global manner by the PDFs and this mismatch badly affects the charmonium phenomenology as α s is not very small and the PDF evolution is limited which results in flat PDF shapes in the mid and low-x region.
Our solution to this issue is to propose a new scale setting,μ F , which is based on the simple criterion that the partonic cross-section vanish at largeŝ (or small z). This is to be understood as that the real-emission contributions coming from the initial partons are entirely absorbed in the PDF. Although less ambitious, this is somewhat equivalent to a resummation picture, yet much simpler to implement.
We have demonstrated the efficiency of this new scale setting for a wide class of rapidity and energy shapes. We have applied this scale setting to η c,b production, a fictitious light elementary scalar bosonH 0 and also for the real BEH boson H 0 with different mass values of a fictitious heavyquark active in the loop. The success of this scale demonstrates that the issue we tackle is not in principle limited to quarkonium but rather to processes occurring at low scales, in particular when some unfavourable effects add up.
Having cured the NLO η c,b phenomenology from these negative cross sections, we have then provided predictions for SPD, AFTER@LHC, LHC up to FCC energies. Naturally, our NLO η c predictions bore on PDFs at low scales. Indeed, our scale choice for η c amounts to 2m c / √ e, thus 1.82 GeV which is well within the usual scale range for a produced system whose mass is 3 GeV. However, at such a scale, conventional PDFs such as PDF4LHC15 nlo 30 which we used exhibit a local minimum for x close or below 0.001. This translates into a distortion of the energy and rapidity dependence of the cross-section in the TeV range that nobody would expect to see in experimental data. We thus encourage the PDF community to see what would be the outcome of a global NLO PDF fit preventing a local minimum in xg(x), which is by the way absent in NNPDF31sx nlonllx as 0118, where resummation effects have been taken into account in the PDF fit, and in JR14NLO08, where gluons are evolved from a lower scale. This local minimum also seems to contradict the recent analysis of Flett et al. [64].
Similarly, we encourage the LHC experimental community to study such P T -integrated quarkonium cross section, despite the likely challenging decay channels which should be used. η c has already been studied by LHCb at finite P T (P T > 6 GeV), they can definitely push down to 0 with a targeted effort, hopefully motivated by our study. As for η b , which remains unobserved in hadroproduction, we have gathered some suggestions on how to extract its production cross section at the LHC. Indeed, these cross sections are definitely very large. Once they are available, we believe them to be ideal to better constrain the gluon PDF at low scales, and thus the gluon content of the proton in general, despite the remaining theoretical uncertainties inherent to quarkonium production. Extension to nuclei could then be considered along the lines of [103,104].
In the near future, it remains to be investigated the effect of such a scale setting for χ c and J/ψ production. The situation for the χ c and J/ψ is much worse [37,39] than that for η c as one already encounters negative yields in hadroproduction at √ s as low as a couple of hundreds of GeV. With our scale setting, we expect that both χ c and J/ψ NLO cross section will stabilise and give physical cross-section results which can then be used in NLO PDF fits.