Revisiting NLO QCD corrections to total inclusive J/psi and Upsilon photoproduction cross sections in lepton-proton collisions

We revisit inclusive J/psi and Upsilon photoproduction at lepton-hadron colliders, namely in the limit when the exchanged photon is quasi real. Our computation includes the next-to-leading-order (NLO) alpha_s corrections to the leading-order contributions in v. Similarly to the case of NLO charmonium-hadroproduction processes, the resulting cross sections obtained in the MS-bar factorisation scheme are sometimes found to be negative. We show that the scale-fixing criteria which we derived in a previous study of eta(c) production successfully solves this problem from the EicC all the way up to the FCC-eh energies. We then elaborate on how to study a scale uncertainty akin to that derived by scale variations when one fixes a scale. In turn, we investigate where both J/psi and Upsilon photoproduction could be used to improve our knowledge of gluon content of the proton at scales as low as a couple of GeV.


Introduction
Inclusive 1 production of quarkonia in hadron-hadron and lepton-hadron collisions is a potential rich source of information on the hadron structure. As such, it has been thoroughly studied both experimentally and theoretically (see [1][2][3][4][5][6] for reviews). Yet, the mechanisms underlying their inclusive production are still not an object of consensus within the community. This in turn does not encourage one to employ cross-section measurements to extract information on the gluon structure of the proton.
Owing to the presence of an electromagnetic coupling, photoproduction cross sections are smaller than hadroproduction ones which calls for large luminosities to obtain large enough quarkonium data sets. As such, the P T reach of J/ψ HERA data is limited to barely 10 GeV, there is quasi no data on ψ and none on the Υ.
In the present analysis, we focus on the P T -integrated yield which was surprisingly seldom studied at HERA. Indeed most of the inclusive data set have been selected with the minimal P T of 1 GeV. This cut however introduces a strong sensitivity on the P T spectrum of the cross section in a region where it is not necessarily well controlled. By itself, such a yield is not directly connected to the total number of J/ψ produced for which we believe the theory predictions to be more robust. One reason for such a cut is probably the difficulty to obtain numerically stable NLO results when they appeared [26,27] 3 . In fact, as we will discuss, these were probably due to the appearance of large nega-tive NLO contributions which we will address here along the same lines as our recent study on η c production [29].
Moreover, such P T -integrated cross sections will be easily measurable at high energies with a very good accuracy at the planned US Electron-Ion Collider (EIC) [30], but also at future facilities such as the LHeC [31] or FCCeh [32], thus in a region where the gluon PDFs are not well constrained. Measurements at lower energies at AMBER-COMPASS++ [33] and the EicC [34] would then rather probe the valence region, which could happen to be equally interesting.
In our study, like in the previous one [7], we will focus on the aforementioned CSM [10], corresponding to the leading-v contribution in NRQCD whose NLO QCD corrections are in principle known since the mid nineties [26,27]. As we revisited in [7], the impact of QCD corrections to J/ψ inclusive photoproduction steadily grows when P T increases. This can be traced back to the more favourable P T scaling of specific real-emission contributions. The same has been observed in several quarkoniumhadroproduction processes [35][36][37][38][39][40][41][42]. On the contrary, one expects, at low P T , a more subtle interplay between the contribution of these real emissions near the collinear region and the loop corrections. This has for a long time been understudied. We thus aim at discussing it here in detail.
The structure of our Letter is as follows. Section 2 outlines our methodology to compute vector-quarkonium inclusive photoproduction cross sections at NLO accuracy including a discussion of the reason for negative NLO cross sections, our proposed solution and a discussion on how to account for scale uncertainties in a computation when one scale is fixed. Section 3 gathers our prediction for future measurements at lepton-hadron colliders along with a discussion of the corresponding theoretical uncertainties. Section 5 gathers our conclusions.
2. J/ψ and Υ photoproduction up to one-loop accuracy

Elements of kinematics
We will consider the process γ(P γ )+ p(P p ) → Q(P Q )+X where the photon γ(P γ ) is emitted by an electron e(P e ). Let us then define s ep = (P e + P p ) 2 ≈ 4E e E p + m 2 p (E e(p) is the electron (proton) beam energy, m p is the proton mass) and s γp = W 2 γp = (P γ + P p ) 2 . We can then introduce x γ as P γ = x γ P e such as s γp = x γ s ep . As announced, in the present study, P 2 γ 0. Diffractive contributions are suppressed for increasing P T and away from the exclusive limit, i.e. when the quarkonium carries nearly all the photon momentum. A cut on P QT is usually sufficient to get rid of them, which we do not wish to apply here. One can however cut on a variable called elasticity, defined as z = P Q ·P p P γ ·P p . z indeed corresponds to the fraction of the photon energy taken by the quarkonium in the proton rest frame, with the proton momentum defining the z axis. It can be rewritten as z = 2 E p m T W 2 γp e y in terms of the quarkonium rapidity y (with y and E p being defined in the same frame) and the quarkonium transverse mass, m T = M 2 Q + P 2 QT with M Q being the quarkonium mass. Such diffractive contributions are known to lie at z → 1 [13]. At low z, resolved-photon contributions can appear as important where only a small fraction of the photon energy is involved in the quarkonium production. At HERA, they had a limited impact [6,13]. At lower energies, like at the EIC, their impact should be further reduced. On the contrary, at the LHeC or FCC-eh, their impact might be sizable even at moderate z. Nonetheless, their modelling requires a good control of the contributions from gg and gq channels. However, our understanding of the very same channels in inclusive quarkonium hadroproduction, especially at low P T [1], is clearly limited. When comes the time for the building of these future lepton-hadron colliders, it will be needed to re-evaluate the impact of the resolvedphoton contributions at low z and high energies. For the time being, we simply note that imposing a lower bound on z would not alter our conclusions at all, precisely because it does not correspond to the low-x region in the proton.

The Colour-Singlet Model
As mentioned earlier, there is no agreement on which mechanism is dominant in quarkonium production. The most popular approaches are: the Colour-Singlet Model (CSM) [9][10][11], the Colour-Evaporation Model (CEM) [43,44] and the Non-Relativistic QCD (NRQCD) [12], whose leading-v contribution is the CSM for S -wave quarkonia. These mechanisms mainly differ in the way they describe hadronisation. The factorisation approach of the CSM is based on considering only the leading Fock states of NRQCD. In the CSM, there is no gluon emission during the hadronisation process, and, consequently, the quantum state QQ does not evolve during the binding. Thus, spin and colour remain unchanged.
The matrix element M to produce a vector state Q + {k}, where {k} is a set of final state particles, from the scattering of the partons ab in CSM is: where N(λ|s 1 , s 2 ) (resp. δ ii / √ N c ) is the projector onto a vector (resp. CS) state and M(ab → QQ + {k}) is the amplitude to create the corresponding heavy-quark pair. When one then sums over the heavy-quark spins, one obtains usual traces which can be evaluated without any specific troubles. In fact, such a computation can be automated at tree level as done by HELAC-Onia [45,46]. In the present case of inclusive photoproduction of a vector Q, there is a single partonic process at Born order, αα 2 s , namely γg → Qg (see Fig. 1a). One could also consider γQ → QQ at the same order, but we have shown it to be small at low P T [7].
The value of R(0), the Q radial wave function at the origin in the configuration space, can be in principle extracted from the leptonic decay width computed likewise in the CSM. The latter is known up to NLO [47] since the mid 70's, up to NNLO since the late 90's [48,49] and up to N 3 LO since 2014 [50]. However, as discussed in Appendix A, the short-distance amplitude receives very large QCD radiative corrections which translate into significant renormalisation and NRQCD-factorisation scale uncertainties. These essentially preclude drawing any quantitative constraints on |R(0)| 2 from the leptonic decays width.
In principle, we should thus associate to it a specific theoretical uncertainty which is however supposed to only affect the normalisation of the cross sections. In what follows, we will employ a similar value as Krämer [27], 1.25 GeV 3 for the J/ψ and 7.5 GeV 3 for the Υ(1S ). As for the masses, we will use m c = 1.5 GeV and m b = 4.75 GeV. Let us recall that within NRQCD  The hadronic cross section is readily obtained by folding the partonic cross section, dσ (0) γg , with the corresponding PDFs and, if relevant, summing over the parton species. Generically, one has: where µ F is the factorisation scale and f i (x, µ F ) is the PDF that gives the probability that a parton i carries a momentum fraction x of the parent proton. At LO, dσ γi identifies to dσ (0) γg .

The NLO corrections and their divergences
At order αα 3 s , two categories of new contributions arise. Those from the real emissions represented by Fig. 1b,1c, 1d, 1e and from virtual emissions (or loops) represented by Fig. 1f, 1g, 1h. Specific topologies of the former category benefit from P T /M Q enhancement factors which make them leading at large P T . This for instance justifies to employ the NLO approximation [7]. When P T is integrated over, a priori all these contributions should be accounted for. Such a computation was first carried out by Krämer [26] in the mid 1990's. We will briefly outline now what it amounts to.
As usual, one critical feature of such NLO computations in collinear factorisation [51] is the appearance of various types of singularities. Of particular relevance for our discussion are the collinear divergences from initialstate emissions which arise when one cannot distinguish two massless particles, with an angle between them close to zero. These remain because the initial states are fixed by the kinematics of collinear factorisation and, consequently, are not fully integrated over. These singularities are absorbed inside the MS-renormalised PDFs via the processindependent Altarelli-Parisi Counter Terms (AP-CT). These AP-CT introduce a µ F dependence in the partonic cross section, which would, in an all-order computation, cancel that introduced by the PDF scale evolution governed by the DGLAP equation.

The cross section in terms of scaling functions
For our analysis, we have found useful to employ the NLO cross-section decomposition in terms of scaling functions derived by Krämer [27]. Using FDC [38,52], we have reproduced his results (scaling functions as well as hadronic cross sections) and, with the appropriate parameter choices and kinematical cuts, those of [53].
Indeed, the advantage of considering P T -and zintegrated cross sections is that the hadronic photoproduction cross sections can be recast in terms of a simple convolution of the PDF and scaling functions of a single scaling variable 4 . This allows one to outline the structure of the result to better understand some specific behaviour (like the scale dependencies discussed in the previous section, hence the importance of negative contributions to the cross section) of the NLO yield. This formulation is also useful because it allows one to economically vary parameters like the c.m. energy, the heavy-quark mass, the renormalisation and factorisation scales.
Along the lines of Krämer [27], we express the partonic cross section as 5 : γg and c (1) γg from the αα 3 s (NLO) γg contributions and c (1) γq and c (1) γq from the αα 3 s (NLO) γq contributions. c (1) γg encapsulates contributions 6 from both real-and virtual emissions. If it had contained only virtual contributions, it would scale like c (0) γg and eventually vanish at largeŝ. This implies that the asymptotic value of c (1) γg entirely comes from the real emissions. c (1) γg only includes real emissions and comes along with an explicit µ F dependence from the AP-CT. The last term, whose form is generic, comes from the renormalisation procedure. The hadronic cross section is then obtained according to Eq. (2).  Already at this stage, one can note then that, at largeŝ, the NLO cross section will be proportional to ln(M 2 Q /µ 2 F ) and a process-dependent coefficient, , which only comes on the real-emission contributions, like for η Q hadroproduction [29]. Fig. 3 shows the √ s γp -dependence of σ γp for J/ψ photoproduction integrated over z < 0.9 and P T , for different choices of µ R and µ F among M J/ψ × (0.5, 1, 2), using the CT18NLO PDF set [54] and with a 20% feed-down contribution from ψ decay (like in [7]). We expect the b feed down on the P T -integrated yields to be on the order of 5% and we do not include it as it can be experimentally removed. For Υ(1S ) photoproduction, we have estimated 7 the feed-down contributions from Υ(2S ) to be 12.5% and from Υ(3S ) to be 2.2%.  The long dashed grey curve is the LO cross section for µ R = M Q and µ F = 0.86 M Q . We have checked that it remains positive for any µ R and µ F scale choice which is expected provided that the PDFs are positive. It happens to reasonably account for the available experimental values if one notes that the theoretical uncertainties from the scales, the mass and R(0) (not shown) are significant. All other curves represent the NLO cross section for different scale choices. In two cases, the NLO cross section becomes negative as s γp increases. Anticipating the results of the next section, the same behaviour is obtained with other well-known PDF sets (MSHT20 [58], NNPDF31 [59]). As for Υ(1S ), the cross section remains positive in the considered energy range for any realistic scale choice, like for η b hadroproduction up to 100 TeV [29]. Let us now discuss the origin of such an unphysical behaviour for J/ψ photoproduction and propose a solution to it.

2.5.
A new scale prescription to cure the unphysical behaviour of the NLO quarkonium photoproduction cross section From the above discussion, there can only be two sources of negative partonic cross sections: the loop amplitude via interference with the Born amplitude and the real emissions via the subtraction of the IR poles from the initial-emission collinear singularities. As we will argue now, the latter subtraction is the source of the negative cross section which we have just uncovered. As it was mentioned before, such divergences are removed by subtraction into the PDFs via AP-CT and the high-energy limit of the resulting partonic cross section takes the form: are the coefficients of the finite term of NLO cross section in the highenergy limit. As can be seen from Fig. 2, A γi is negative for z < 0.9, i.e. −0.29. It is also clear from Fig. 2 Unless µ F is sufficiently smaller than M Q in order to compensate A γi , lim s→∞σ NLO γi is negative, like for η Q [29] and it is another clear case of oversubtraction by the AP-CT. Indeed, in this limit, the virtual contributions are suppressed; only the real emissions contribute via their square. As such, they can only yield positive partonic cross sections before the subtraction of the initial-state collinear divergences. Since, after their subtraction, the partonic cross section is negative, it has to come from the AP-CT. This is what we refer to as oversubtraction by the AP-CT.
In principle, the negative term from the AP-CT should be compensated by the evolution of the PDFs according to the DGLAP equation. Yet, for the µ F values on the order of the natural scale of these processes, the PDFs are not evolved much and can sometimes be so flat for some PDF parametrisations that the largeŝ region still significantly contributes. This results in negative values of the hadronic cross section. Indeed, A γg and A γq are process-dependent, while the DGLAP equations are process-independent, which necessarily makes the compensation imperfect. Going to NNLO and even higher, this should naturally improve. If one has only NLO computations, this is however greatly problematic. A solution to this problem is [29] to force the partonic cross section to vanish in this limit, whose contribution should in principle be damped down by the PDFs.
According to this prescription, one needs to choose µ F such that lim s→∞σ NLO γi = 0. It happens to be possible since A γg = A γq . This amounts to consider that all the QCD corrections are in the PDFs [29]. From Eq. (4), we have: Using the scaling function of Fig. 2 when one fully integrates over P T and over z < 0.9, one getsμ F = 0.86M Q . From now on, all our NLO results will be shown with this value of the factorisation scale.  On Fig. 4, one can see the LO (in blue) and NLO (in red) µ R dependence of σ γp for J/ψ photoproduction, still using CT18NLO and integrated over P T and over z < 0.9, at two values of √ s γp = 20 GeV (short and long dashdotted lines) and √ s γp = 100 GeV (solid and dashed lines).
In both cases, the µ R sensitivity is drastically reduced at NLO. However, one notes that at the higher energy, for µ R ∼ M J/ψ , σ NLO γp is twice smaller than σ LO γp (see the arrows by the y axis). This is due to a large negative contribution from the loops (see the negative dip in the c (1) γg in Fig. 2). Since the LO and NLO cross section are however similar for µ R ∼ 2M J/ψ , the question of the natural scale of the process naturally arises. In fact, as the Born process is γg → Qg, it appears reasonable to consider √ŝ rather than M Q . A quick LO computation for the J/ψ case shows that √ ŝ ranges from 4 GeV at low hadronic energies up to even 10 GeV at high hadronic energies. In what follows, we thus consider µ R within the range [2.5 : 10] GeV for J/ψ and, for Υ(1S ), µ R ∈ [8 : 32] GeV (with µ R = 5(16) GeV being the center of this range in the two cases) at both LO and NLO and for µ F at LO. The procedure for µ F at NLO is discussed in the next section.

Scale fixing and theoretical uncertainties
Just as we have discussed above, theoretical uncertainties of perturbative computations including radiative corrections arise from the appearances of unphysical scales, usually µ R and µ F . In principle, predictions of physical observables, e.g. a cross section, should not depend on them. This would only be so if all radiative corrections could be accounted for. At NLO, it is far from being the case and the scale dependences can be strong, so strong that some results are sometimes unphysical like in the process under discussion for some (reasonable) µ F values.
In general, the scale dependence is however mild (like for µ R here) and evaluating observables at natural values of the scales, i.e. those entering the kinematics of the process, yield to good predictions, which usually improve at NNLO and so forth. Since the scale dependences are meant to disappear in all-order computations, one can revert the argument and consider that the scale dependences give us some information about the impact of higher orders. This is why one varies the scale, like we have discussed at the end of the previous section, hoping to seize up some theoretical uncertainties from Missing-Higher Orders (MHO) in the jargon. This is however not done without ambiguity, to say the least, as the resulting uncertainties fully depend on the range of scale variation, conventionally chosen to be a multiplicative factor 2 for historical but unclear reasons 8 . Whatever the variation should be, one then faces an apparent impossibility with our scale-fixing prescription: how to vary a scale which is fixed?
To address this issue, it is necessary to go back to the very motivation of scale-fixing criteria. Like for our scale prescription, the other prescriptions (PMS [61], BLM-PMC [62][63][64], FAC [65,66] 9 ) stem from physical pictures: the result should be stable, the results should be maximally conformal, the convergence should be as fast as possible or, like in our case, the result should exhibit no over-subtraction of collinear singularities inside the PDFs. The question is then how much one can depart from this expectation? Presumably, the scale value from a given prescription will not be the same at NLO and NNLO for instance. One could then try to derive the NLO scale from its formal expression artificially corrected by a typical NNLO corrections scaling like α s /π, on the order of unity in our case. This would provide a range of prescribed scales that we could plug in the NLO computation to get a scale uncertainty. In our case, instead of forcing log M 2 Q µ 2 F + A γi to vanish, one could consider a range of µ F such that its absolute value is bounded below unity. As such, we would probe the robustness of the scale-fixing criteria against expected higher-order corrections. This seems a very elegant solution for which the range of µ F is 1/ √ e ≤ µ F /μ F ≤ √ e. There is however a caveat to compare it with the usual way, not because of the above solution, but because the conventional method of scale variation by 2 is purely arbitrary 10 . For meaningful comparisons, the range of variation should be accounted for. After all, what one looks after is the scale dependence. A natural way to proceed should instead be to compute dσ/d ln µ which would then be a local estimation of the scale dependence. Assuming that the variation by 2 is an approximate way to numerically evaluate this derivative, the connection between both requires to include a ln 2 factor and ln 2/ ln e for the proposal above. We refer to Appendix B for more details.
In this context, we will show a NLO µ F scale uncertainty derived from ln 2 × dσ/d ln µ, in fact very similar to that obtained from 1/ √ e ≤ µ F /μ F ≤ √ e rescaled by ln 2/ ln e.

Results
Having discussed our methodology, let us now present and analyse our results for J/ψ and Υ(1S ) photoproduction cross sections computed at NLO with theμ F prescription. certainty has essentially no physical meaning: a variation by √ 2, or 3, or whatever else not far from 2, would be equally acceptable while yielding different uncertainty estimates. We are in fact surprised how much this issue is underdiscussed when theory is confronted to experimental data. 11 We note here that the mass and R(0) uncertainites are highly kinematically correlated and essentially translate into a quasi global offset. This thus why we focus on the µ R , µ F and PDF uncertainties. NNPDF31 nlo as 0118 hessian [59]. Let us first discuss Fig. 5a. The LO cross section (2 solid and 2 dashed grey lines) relatively well describes the experimental data points in red with somewhat large uncertainties at large √ s γp . The NLO cross section is systematically smaller and one notes that the NLO µ R uncertainty (red hatched band) is reduced compared to the LO one, as expected from Fig. 4. The NLO µ F uncertainty we obtained is well-behaved without any negative values anymore and is smaller than the LO one at large √ s γp .
The PDF uncertainties at NLO from CT18NLO, MSHT20nlo as118 and NNPDF31 nlo as 0118 hessian are shown by respectively the blue, magenta and orange hatched bands. At large √ s γp , which corresponds to the low-x region in the proton, they naturally grow and eventually become larger than the µ R uncertainty. Even though it is not an observable physical quantity, we note that with our present set-up (scheme and scale choice) the relative contribution from the γq fusion channel is relatively constant and close to 5% from 20 GeV and above, about 95% then comes from γg fusion. The three PDFs we have chosen are representative of what is available from fixed-order analyses. However, they are known [67] to be artificially suppressed at low x and low scales and can even show a local mininum which then distorts in the energy dependence of η c hadroproduction cross section [29]. As such, the possibility remains, that even within their increasing uncertainties, such present PDF sets are possibly unsuited to reliably describe quarkonium production data. Using the latter to redetermine the former is then certainly something which should be attempted. The increase of the PDF uncertainty is even more visible on the relative uncertainty plots 12 Fig. 5a (lower panel) for µ R = 5 GeV. Above 300 GeV, these are clearly larger than the µ R one which slightly grows above 50 GeV due to the on-set of the negative contributions from the loop corrections (see below) and the µ F one which also slightly grows above 20 GeV, where it happens to coincidentally vanish for a wide range of values. As for the Υ(1S ) case, shown on Fig. 5b, the reduction of the µ R uncertainty at NLO is further pronounced while the PDF and µ F uncertainties remain similar.
In the J/ψ case, it is clear that it will be important to have at our disposal computations at NNLO accuracy. As Krämer noted [27] long ago, the "virtual+soft" contributions, encapsulated in c (1) , are significantly more negative than for open heavy-flavour production [69]. He suggested that this destructive interference with the Born order amplitude could be due to the momentum transfer of the ex- 12 The LO relative scale uncertainties are computed as ±(σ max γp − σ min γp )/(2σ cen γp ), where σ max/min γp is the maximum/minimum values of σ γp obtained by varying µ R , µ F in the quoted range and σ cen γp is the cross section evaluated with µ F = µ R = µ 0 , with µ 0 = 5(16) GeV for J/ψ(Υ(1S )). At NLO, the µ R uncertainty is calculated as ±(σ max γp − σ min γp )/(σ max γp + σ min γp ), while the µ F uncertainty is estimated with the local method (see Sec. 2.5 and Appendix B), normalised to σ cen γp (µ F =μ F , µ R = µ 0 ). For the PDF uncertainties, we used the normalised upper and lower PDF uncertainties [68] for µ R = µ 0 , where again µ 0 = 5(16) GeV for J/ψ (Υ(1S )).   changed virtual gluon, more likely to scatter the QQ pair outside the static limit (p 0). At NNLO, these one-loop amplitudes will be squared, the two-loop amplitudes will interfere with the Born amplitudes and the amplitudes of the one-loop corrections to the real-emission graphs will also interfere with the real-emission amplitudes. Unless the latter two are subject to the same strong destructive interference effect, one might expect relatively large positive NNLO corrections bringing the cross section close to the upper limit of the LO range and then in better agreement with existing, yet old, data. At NNLO, we also expect a further reduction of the µ R and µ F 13 uncertainties. This is particularly relevant especially around 50 − 100 GeV, which corresponds to the EIC region. This would likely allow us to better probe gluon PDFs using photoproduction data. Going further, differential measurements in the elasticity or the rapidity could provide a complementary leverage in x to fit the gluon PDF, even in the presence of sub-leading v colour-octet contributions. Indeed, these would likely exhibit a very similar dependence on x. As we will see now, the expected yields at future facilities, in particular for charmonia, are clearly large enough to perform such differential measurements. 13 It is legitimate to expect the oversubtraction by the AP counter terms to be reduced at NNLO associated with a reduction of the sensitivity on µ F . It might also be less sensitive on the PDF shape [29,70].
Let us now look at electron-proton cross sections as functions of √ s ep for J/ψ (Fig. 6a) and Υ(1S ) photoproduction ( Fig. 6b). To obtain them, Eq. (3) was convoluted with the corresponding proton PDFs and a photon flux from the electron. We have used the same photon flux as in [7]. On Fig. 6a and Fig. 6b, the same colour code and the same parameters as for Fig. 5a and Fig. 5b have been used. For σ ep , one can see the same trends for the µ R and PDF uncertainties as for σ γp . It is only at the LHeC energies and above that one could expect to constrain better the PDF uncertainty with such total cross section measurements unless we have at our disposal NNLO computations with yet smaller scale uncertainties.  (in GeV) corresponding to future facilities (using CT18NLO, µ R = 5 GeV for J/ψ and µ R = 16 GeV for Υ(1S ), µ F =μ F ) for detect = 85% via the decay channels to µ + µ − and e + e − , namely In Table 1, we provide estimations of the expected number of J/ψ and Υ(1S ) possibly detected at the different ep c.m. energies of planned experiments. As it can be seen, the expected yields are always very large for J/ψ which will clearly allow for a number of differential measurements in z, y or √ s γp . These could then be used to reduce the impact of partially correlated theoretical uncertainties, from the scales and the heavy-quark mass affecting these photoproduction cross sections, in order to bring about some additional constraints on the PDFs at low scales, in particular the gluon one. For Υ(1S ), the yields should be sufficient to extract cross sections at the EIC, LHeC and FCC-eh even below their nominal luminosities. One can also estimate the expected number of detected ψ , Υ(2S ), Υ(3S ) using the following relations N ψ 0.07 , derived from the values of 14 |R Q (0)| 2 and of the branching fractions to leptons. Using the above relations and the values in Table 1, one can see that the yield of ψ should be measurable everywhere and the yields of Υ(2S ) and Υ(3S ) are close to about half of that of Υ(1S ) and should be measurable at the EIC, LHeC and FCC-eh. The proximity between the Υ(nS ) yields follows from their similar |R Q (0)| 2 and leptonic branchings.

A note on P T -differential cross sections
As announced, the present study focuses on the fully P T -integrated case. If, instead, one is interested in P Tdifferential cross sections, our prescription for a single scale process would not work as it stands. Indeed, one has to 14 The relation for ψ was estimated using |R ψ (0)| 2 = 0.8 GeV 3 . be cautious that a class of real-emission NLO corrections (Fig. 1c) is kinematically enhanced by (P T /M Q ) 2 with respect to the Born contributions. Considering the expected P T behaviour of the scaling functions entering the definition ofμ F , one finds that c (1) γi (which contains these contributions) is enhanced with respect to c (1) γi (which only arises from the collinear contributions).
This would result inμ F (P T ) scaling as M Q × exp(P T /M Q ). This can easily be explained if one reminds that our scale prescription effectively amounts to recast, at large partonic energies, all the NLO corrections in the PDF folded with the Born cross section. This would then be done, at large P T , by trying to get those P T -enhanced contributions entirely from the PDF evolution via an unphysically large µ F . At large P T , the natural scale should instead be on the order of P T and not M Q × exp(P T /M Q ).
Common dynamical scale choices for P T -differential cross sections are (0.5, 1, 2) × m T with m T = M 2 Q + P 2 T . Two possible choices scaling like m T at large P T and compatible withμ F when integrated over P T would be T . At HERA energies [13], P 2 T 2.5 GeV 2 for J/ψ, this gives α = 0.77 and β = 0.7. The former choice is shown on Fig. 7 (red boxes) and compared to cross section obtained with a fixed µ F =μ F for different µ R (hashed red histogram). All choices give similar results which are compatible with the latest H1 data [18]. Predictions for the EIC at 140 GeV are also given. They confirm predictions given in our previous study [7] with an approximate NLO computation.

Conclusions
Like for other charmonium production processes [29,71,72], we have observed the appearance at NLO of nega-tive total cross section which we attribute to an oversubtraction of collinear divergences into the PDF via AP-CT in the MS scheme. We applied theμ F prescription proposed in [29], which up to NLO corresponds to a resummation of such collinear divergences in High-Energy Factorisation (HEF) [73]. Expressing this integrated cross section in terms of scaling functions exhibiting its explicit µ R and µ F scale dependencies, we have found that, for z < 0.9, the optimal factorisation scale isμ F = 0.86M Q which falls well within the usual ranges of used values. Like for η c hadroproduction, such a factorisation-scale prescription indeed allows one to avoid negative NLO cross sections, but it apparently prevents one from studying the corresponding factorisation-scale uncertainties. We have thus elaborated on two approaches to study such scale uncertainties when one scale is fixed by a physical argument. We were not aware of such ideas before and these clearly deserve dedicated studies in the future.
We have seen that the NLO µ R uncertainties get reduced compared to the LO ones but slightly increase around 50 GeV, because of rather large (negative 15 ) interferences between the one-loop and Born amplitudes. As mentioned before, these "virtual+soft" contributions are significantly more negative than for open heavy-flavour production. While Krämer suggested that the difference could stem from the static limit (p 0) specific to the non-relativistic quarkonia, from which one easily departs when gluon exchanges occur, it will certainly be very instructive to have NNLO computations to see whether such one-loop amplitudes squared would bring the cross section back up close to the LO one or whether the interference between the twoloop and the Born amplitudes and between the real-virtual and the real amplitudes would be also negative and large.
Our evaluated µ F uncertainties at NLO are also reduced compared to those at LO at large √ s γp which were totally out of control before the application of theμ F prescription. We see that as a very encouraging sign. In any case, at NNLO, it is reasonable to expect a further reduction of the scale uncertainties compared to the NLO results. We have also briefly addressed the application of our scale prescription to P T -differentical cross sections and proposed choices scaling like m T at large P T while compatible withμ F = 0.86M Q when P T is integrated over.
We have also qualitatively investigated the possibility to constrain PDFs using future J/ψ and Υ(1S ) photoproduction data. Indeed, present PDF sets are possibly [29] unsuited to reliably describe high-energy quarkonium production data. Restricting to conventional sets, we have seen that PDF uncertainties get larger than the (NLO) µ R and µ F uncertainties with the growth of the γp c.m. energy, in practice from around 300 GeV, i.e. for x below 0.01. Although this is above the reach of the future EIC, we hope that with NNLO predictions at our disposal in the future, with yet smaller µ R and µ F uncertainties, one could set novel constraints on PDFs with such EIC measurements. Given our estimated counting rates for 100 fb −1 of ep collisions, we expect that a number of differential measurements will be possible to reduce the impact of highly or even partially correlated theoretical uncertainties, including the contamination of higherv corrections such as the colour-octet contributions. These, along with the forthcoming HL-LHC measurements [21], should also definitely help to improve our understanding of the quarkonium-production mechanisms.
Strictly speaking our predictions for J/ψ and ψ only regard the prompt yields. An evaluation at NLO of the beauty production cross section points at a feed-down fraction at the 5% level. Given the larger size of the other uncertainties and the possibility to remove it experimentally, we have neglected it. In general though, it will be useful to have a dedicated experimental measurements at the EIC at least to measure the beauty feed down. It may become more significant at low z where the resolved-photon contribution could set in at high √ s ep , like we have seen [7] it to become the dominant source of J/ψ at large P T . where n l f is the number of active light flavours, α is the electromagnetic coupling constant, α s is the strong interaction coupling, M Q is the Q mass, e Q is the magnitude of the heavy-quark charge (in units of the electron charge). In Table A Scale uncertainties are historically evaluated in pQCD by a multiplicative variation by a factor of 2. Assuming that the corresponding µ uncertainty, defined as σ±∆σ, is evaluated via

Acknowledgements
it is easy to note that wider variations necessarily lead to larger uncertainties and narrower ones to smaller uncertainties. However, as pointed out above, the common practice is to use such a factor 2 which then lead to the usual 7-or 9-point variation technique to estimate the envelope corresponding to the factorisation and renormalisation scale uncertainty.