Exclusive $J/\psi$ and $\Upsilon$ photoproduction and the low $x$ gluon

We study exclusive vector meson photoproduction, $\gamma p \to V + p$ with $V=J/\psi$ or $\Upsilon$, at NLO in collinear factorisation, in order to examine what may be learnt about the gluon distribution at very low $x$. We examine the factorisation scale dependence of the predictions. We argue that, using knowledge of the NLO corrections, terms enhanced by a large $\ln(1/\xi)$ can be reabsorbed in the LO part by a choice of the factorisation scale. (In these exclusive processes $\xi$ takes the role of Bjorken-$x$.) Then, the scale dependence coming from the remaining NLO contributions has no $\ln(1/\xi)$ enhancements. As a result, we find that predictions for the amplitude of $\Upsilon$ production are stable to within about $\pm 15\%$. This will allow data for the exclusive process $p p \to p\Upsilon p$ at the LHC, particularly from LHCb, to be included in global parton analyses to constrain the gluon PDF down to $x\sim 10^{-5}$. Moreover, the study of exclusive $J/\psi$ photoproduction indicates that the gluon density found in the recent global PDF analyses is too small at low $x$ and low scales.


Introduction
The present global PDF analyses (e.g. NNPDF3.0 [1], MMHT2014 [2], CT14 [3]) find that there is a large uncertainty in the low x behaviour of the gluon distribution. There is a lack of very low x data, particularly at low scales. Moreover, the gluon is determined at low x mainly by the DGLAP evolution of deep inelastic scattering, that is, not from a direct measurement 1 , but rather from the derivative dF 2 /d ln Q 2 . As a result, at low scales, the uncertainty on the gluon PDF is large for x < ∼ 10 −3 . On the other hand, the HERA data on diffractive vector meson photoproduction [4,5,6,7,8,9,10,11,12,13,14,15], γp → V + p and the LHC data on the exclusive processes pp → p V p [16,17,18] where V = J/ψ or Υ, and p Pb → p V Pb [19], sample directly the gluon distribution down to x ∼ 10 −5 . For a review of exclusive vector-meson production at the LHC see, for example, Ref. [20].
However, J/ψ and Υ data are not used in the global PDF analyses. One reason is that the corresponding cross sections are described by non-diagonal (skewed) analogues of the PDFs, namely generalised parton distributions (GPDs), due to the different masses of the incoming photon and the outgoing V meson [21]. The GPDs can be related to the PDFs via the Shuvaev transform in the low ξ region [22]. 2 In this work, following [23], we make the physically motivated assumption that the input distribution has no singularities in the right-half Mellin-N plane which implies that this relation holds at NLO to accuracy O(ξ). Another reason is the dependence of the theoretical predictions on the choice of the factorisation scale, µ F . This problem has two independent parts. One part has a technical nature and the other is more physical. Let us discuss them in turn. First, the 'technical' problem, which is related to the convergence of the perturbative expansion at low ξ and low scales. A good illustration is that the NLO amplitude for the exclusive high-energy γp → J/ψ + p process [24] was shown to yield a cross section which varies by up to an order of magnitude for a reasonable variation of µ F . This problem was emphasised recently by Wagner et al. [25]. The strong scale dependence arises because in the DGLAP evolution of low ξ GPDs the probability of emitting a new gluon is strongly enhanced by the large value of ln(1/ξ). Indeed, the mean number of gluons in the interval ∆ ln µ F is [26] n α s N C π ln(1/ξ) ∆ ln µ 2 leading to a value of n up to about 8, for the case ln(1/ξ) ∼ 8 with the usual µ F scale variation interval from µ F /2 to 2µ F . In contrast, the NLO coefficient function allows for the emission of only one gluon. Therefore we cannot expect compensation between the contributions coming from the GPD and the coefficient function as we vary the scale µ F . (At large ξ the compensation is much more complete and provides reasonable stability of the predictions to variations of the scale µ F .) In Section 2, we use the NLO contribution to fix a choice of the factorisation scale for the LO part of the amplitude, which allows the ln(1/ξ) corrections to be resummed.
The second or 'physical' reason why there is a strong scale dependence of J/ψ and Υ photoproduction is due to 'defects' in the presently available PDF sets which are obtained in the DGLAP global analyses of deep inelastic and related hard scattering data. That is, our studies indicate that there is a strong scale dependence in vector-meson photoproduction caused by the unexpectedly small LO contribution in comparison with the NLO correction. This occurs due to the very low input gluon density in the low x domain obtained in the recent global PDF analyses. We conclude that at the input scale, the low x gluon density is underestimated in comparison with that for quarks. The crucial observation is that the number of input gluons (which is parametrised freely in the global analyses) is found to be much less than the number of such gluons emitted by quarksor, to be more precise, gluons associated with the quark (each quark carries a gluon field created by its colour charge). With the present global partons the NLO sea quark contribution to V meson photoproduction grows with decreasing x faster than the input gluons. As a result, at sufficiently high energy, the quark component of the NLO correction to V photoproduction approximately cancels the main LO contribution which is due only to the gluon density. Therefore a small variation of the scale in the NLO component leads to a strong variation of the predicted cross section.
Recently the LHCb Collaboration have presented data for open charm and beauty production in the forward direction [27,28]. These data sample approximately the same kinematic domain as exclusive vector meson production; actually the x values are slightly larger, and the scale for beauty production is also larger. It was shown in [29,30] that the data can be reproduced by NLO QCD using the present global PDFs. Note, however, that again there is a large sensitivity to the choice of factorization scale. Indeed, to reproduce the data one needs to take a scale approximately a factor of two larger than the natural choice µ F = Q 2 + m 2 Q . Since the optimal or appropriate scale for open heavy-quark production is not known at present, these data do not yet reliably determine the low x gluon at low scales. In the present paper we focus on exclusive vector-meson production. Moreover, our goal is not to present a new global PDF analysis or an explicit determination of the low x gluon, but rather to study the possibility to reduce the scale dependence of the prediction, to study the qualitative structure of the NLO amplitude and to give hints of the consequences for the behaviour of the low x gluon at low scales.
In Section 2 we consider the NLO contribution originating from the light quark (singlet) GPD. We show that this part of the NLO amplitude allows us to choose a factorisation scale which sums the ln(1/ξ) enhanced contributions inside the GPD. The remaining part of the NLO contribution has no large ln(1/ξ) factors, and, as shown in Section 3.2, the compensation between the scale dependence of the GPDs and of the coefficient functions for Υ production makes the NLO result sufficiently stable for the data to be included in global parton analyses (just as for the predictions of processes at larger x). The study of the J/ψ photoproduction in Section 3.1 suggests that the gluon density obtained in the existing PDF analyses is too small at low x and low scale. In Section 4 we present our conclusions.

NLO corrections and the choice of µ F
Our aim in this section is to show how knowledge of the NLO contribution for γp → V +p allows us to sum the ln(1/ξ) enhanced contributions by a particular choice of µ F in the LO part of the amplitude (and correspondingly in the NLO coefficient function), which, in turn, reduces the factorisation scale dependence of the predicted cross section. We outline the procedure in Section 2.2, but first we comment on the NLO formalism.

NLO formalism for γp → V + p
The NLO formalism for γp → V p has been presented in [24]. However, before we can apply it to describe the high energy V photoproduction data, we must note one important correction.
The LO partonic amplitude, A g (x/ξ), for the subprocess γ + (gg) → V , can be calculated by considering the fusion of a photon with a pair of on-shell gluons with zero transverse momentum and physical, transverse, polarisations. If we use dimensional regularisation with D = 4 + 2 (as in [24]) to regularize the ultraviolet (UV) and infrared (IR) divergences, then the result is given, before dividing by the number of physical transverse polarisations of the incoming gluon, by In D dimensions, there are D −2 = 2+2 possible transverse directions. Since the entire calculation must be performed in D dimensions, to properly average over the polarisation of the incoming gluon, we must divide (2) by this factor, and not simply by 2. Therefore for the LO partonic amplitude we obtain, contrary to [24] (prior to the first erratum), The difference of a factor of (1 + ) does not alter the finite LO result (as is taken to zero) but does lead to a different NLO result due to the differing counter-terms. Indeed, the extra factor (1 + ) presented in [24] generates extra terms ∼ / in the counter-terms, ∆, giving them a form ∆ = . . .
where 1/ˆ = 1/ + γ E − ln(4π) and γ E is the Euler-Mascheroni constant. However, with the LO partonic amplitude (3) we instead obtain counter-terms containing It is precisely this alteration to the counter-terms which leads us to a different NLO finite result. 3 This result agrees with that of [24] with both errata.
The details of our recalculation of the NLO result will be published separately [33]. Here we simply state the final result for the high-energy limit of the NLO part of the amplitude. In the high-energy approximation, W 2 M 2 V (where W is the c.m. energy of the incoming γp system and M V the mass of the vector meson), the imaginary part of the amplitude dominates and the leading contribution to the NLO correction comes from the region ξ |x| 1. The matrix element is given by Here m is the mass of the c (for J/ψ) or b (for Υ) quark and (for this unpolarised, forward process) with H g , H S the gluon and quark singlet GPDs, respectively. The GPDs can be calculated from the diagonal PDFs, q(x, µ),q(x, µ) and g(x, µ), via the Shuvaev transforms: where with an accuracy O(ξ), see Ref. [23]. Expressions (8) and (9) were used not only for the asymptotic NLO amplitude, (4), but also for the full expression for the complete NLO amplitude used in the numerics presented below. The LO coefficient function, where e i are polarisation vectors and O 1 V is the NRQCD (non-relativistic QCD) matrix element of the cc → J/ψ (or bb → Υ) transition. In order to have a small relativistic correction to O 1 V , we have to calculate the Feynman diagrams assuming that the charm/bottom-quark line has mass One can see directly from the high energy, leading ln(1/ξ), limit, given in (4), that for the choice µ F = m the asymptotic limit of the quark NLO contribution vanishes. It is this observation that will allow us (in Section 2.3) to claim that a suitable value for the factorisation scale in leading logarithm terms is µ F = m.
Before we show how the NLO contribution allows us to fix the scale µ F in the LO term, it is informative to first recall its structure and introduce the notation in terms of the LO prediction.

General procedure
We start with the LO contribution to γp → V + p. It is sketched in Fig. 1(a), and the amplitude is given by the convolution where the sum over a = q, g is understood and C (0) q = 0 for this process. The coefficient function C (0) g is calculated using the non-relativistic vector meson wave function. In general, the relativistic corrections are not small for the J/ψ case. These corrections should be considered together with the three parton (cc+g) component of the wave function; that is, accounting for the rôle of gluons which provide the interaction between charm quarks. For this process, it was shown by Hoodbhoy [34] that the consistent treatment of relativistic corrections may, to good accuracy, be effectively accounted for by choosing, in the non-relativistic formula, the charm quark mass m c = M J/ψ /2. 4 After this, the remaining part of the correction is quite small (a few percent only). So, we may use the non-relativistic J/ψ wave function with m c = M J/ψ /2, to obtain a result with good accuracy.
We return to discuss V = J/ψ or Υ. At LO the coefficient function C (0) g does not depend on µ F , whereas the low x gluon distribution F g depends strongly on the scale. In other words, based on (12), the exclusive V data measure g(x, µ F ), but we do not know the value of the scale at which it has been determined. How does this scale freedom arise? To obtain the LO result the coefficient function -the upper box in Fig. 1(a) -is calculated with on-mass-shell gluons with transverse momenta l t = 0. In this collinear approach the factorisation scale µ F acts effectively as the ultraviolet (UV) cutoff of the logarithmic integral dl 2 t /l 2 t ∝ ln µ 2 F over the gluon transverse momentum in the gluon loop of Fig. 1(a). Formally, in the LO collinear factorisation approach the value of µ F is not known. In principle, the full result does not depend on µ F since the higher-order, NLO, NNLO, ..., corrections compensate the effect of variations of µ F . However, in reality the perturbative series is truncated and the compensation may not be sufficient to provide the scale stability of the theoretical prediction.
On the other hand, we can go beyond the collinear logarithmic approximation by computing the l t integral accounting for the l t dependence 5 of the hard γ + gg → V matrix element, M, shown in the upper box of Fig. 1(a). In comparison with the coefficient function (calculated with l t = 0) now the l t and l 2 dependence of M is included explicitly, and provides the UV convergence of the integral over l t . Formally, in collinear factorisation the difference between the pure logarithmic (ln µ 2 F ) evaluation and the precise calculation of the gluon l t integral is treated as part of the NLO correction. This part of the NLO correction is of kinematic origin and is usually quite large. Fortunately, it can be moved into the LO component of the amplitude, noticeably reducing the remaining NLO correction.
Instead of performing an independent calculation which accounts for the l t dependence, we can remain within the collinear approach and determine the value of the l t integral given that the NLO coefficient function C (1) q of Fig. 1(b) is known. Indeed, Fig. 1(b) is the only diagram for the quark NLO coefficient function. In this approach the incoming quarks are assumed to be on-mass-shell and with zero transverse momenta but the loop integral over l, which contains the l t dependence, is calculated exactly. Since this is the same integral as that which occurs in Fig. 1(a) we can use the result for C to obtain a precise value, J, of the corresponding integral in the LO amplitude of Fig. 1(a). After this we choose a scale µ F = µ 0 which mimics the precise l t integration. That is, with a scale choice satisfying ln µ 2 0 = J we have moved a large contribution from NLO to LO, and can continue to work in the conventional collinear approach, but now with a smaller NLO correction. Figure 1: (a) The LO contribution to γp → V + p, showing the convolution (12). (b) The NLO quark contribution. For these graphs all permutations of the parton lines and coupling of the gluon lines to the heavy-quark pair are to be understood. Here P ≡ (p + p )/2 and l is the loop momentum.

Transfer of part of the NLO to the LO contribution
Since the above observation is crucial, let us demonstrate the procedure in more detail. At NLO level the LO+NLO amplitude at some factorisation scale µ f may be expressed in the form 6 where the F s are the GPDs and where the coefficient function C (0) does not depend on the factorisation scale. Note that we are free to evaluate the LO contribution at a different scale µ F , since the resulting effect can be compensated by changes in the NLO coefficient function, which then also becomes dependent on µ F . Then eq. (13) becomes Note that although the first and second terms on the right hand side depend on µ F , their sum does not (to O(α 2 s )) and is equal to the full LO+NLO amplitude calculated at the factorisation scale µ f . In (13) the NLO coefficient function C (1) is calculated from Feynman diagrams which are independent of the factorisation scale. How does the µ F dependence of C (1) rem in (14) actually arise? It occurs because we must subtract from C (1) the α s term which was already included in the LO contribution. 7 Since the LO contribution was calculated up to some scale µ F the value of C (1) after 6 For ease of understanding we omit the parton labels a = g, q on the quantities in (13) and the following equations. The matrix form of the equations is implied. 7 Simultaneously this subtraction also provides the infrared convergence of C (1) .
the subtraction depends on the value µ F chosen for the LO component. The change of scale of the LO contribution from µ f to µ F also means we have had to change the factorisation scale which enters the coefficient function C (1) from µ f to µ F . The effect of this scale change is driven by the LO DGLAP evolution, which is given by where V denotes the skewed splitting functions. That is, by choosing to evaluate A (0) at scale µ F we have moved the part of the NLO (i.e. α s ) corrections given by the last term of (15) from the NLO to the LO part of the amplitude. In this way C (1) becomes the remaining µ F -dependent coefficient function C rem (µ F ) of (14). In spite of the unusual form of (14), with two different scales µ f and µ F , it is an exact equality at NLO and could, in principle, be generalised to higher orders. 8

Large ln(1/ξ) terms resummed in the LO contribution
The idea is to use the above procedure to reduce the scale dependence of the exclusive V photoproduction amplitude. Unfortunately the value of µ F is just one number, while C (1) (x/ξ) is a function of the ratio x/ξ. So we have no chance to nullify C (1) completely. Nevertheless we can nullify the most important NLO correction which is enhanced by the large value of ln(1/ξ). This contribution is generated by the integral (dx/x) in the C (1) rem (x/ξ, µ F )F (x, ξ, µ f ) convolution of eq. (14) over the kinematic region of 1 x ξ. That is, we choose a scale µ F = µ 0 which nullifies C (1) (x/ξ) in the limit of x ξ. It can be seen from (4) that the scale µ F = µ 0 = m ensures that the ln(1/ξ)-enhanced NLO corrections completely vanish for both the quark and the gluon components.
From the NLO viewpoint the particular choice of µ F = m in the LO part is irrelevant. The corresponding order α s effect is exactly compensated by the remaining C (1) rem term. However, the variation of µ F in the LO part affects not only the O(α s ) terms but the higher-order α s contributions as well. Therefore, in this way, we resum the important large ln(1/ξ)-enhanced part of the higherorder α s corrections inside the parton distribution convoluted with the LO coefficient function and improve the convergence of the perturbative series. 9 Actually our approach is rather close in spirit to the k t -factorisation method. Indeed, there, the value of the factorisation scale is driven by the structure of the k t = l t (or the virtuality, Q 2 ) integral 8 For example, if the NNLO contribution were known, then we will have three scales: µ f , µ F ≡ µ LO , and µ NLO , where the NNLO correction to (14) takes the form α 2 s C rem (µ LO , µ NLO ) ⊗ F (µ f ). The scale µ NLO is fixed to nullify this term in the limit x ξ, and hence further reduce the sensitivity to variations of the scale µ f . 9 In particular, the most important factorisation scale dependence, enhanced by large ln(1/ξ), is caused by the double log terms, [α s ln(µ F ) ln(1/x)] n , generated in the axial gauge by ladder-type diagrams. For GPDs this ladder contribution was studied in [36], where it was shown, in the large ln(1/x) limit, that the skewed splitting functions in the diagrams of Fig. 1. 10 In the k t -factorisation approach this k t integral is written explicitly, while the parton distribution unintegrated over k t is generated by the last step of the DGLAP evolution, similar to the prescription proposed in Refs. [37,38]. Now, using the known NLO result, we account for the exact k t integration in the last cell adjacent to the LO hard matrix element. This hard matrix element M, shown in the upper box of Fig. 1(a), provides the convergence of the integral at large k t . In this way it puts an effective upper limit of the k t integral, which plays the role of an appropriate factorisation scale.
The details of the prescription, for the case of the high-energy Drell-Yan process, were discussed in [39]. Indeed, Drell-Yan production of low-mass lepton pairs at high rapidity is another process for which the NLO prediction depends sensitively on the choice of the factorisation scale, unless the ln(1/x) enhancements are first resummed in the incoming parton distributions. It was found in Ref. [39] that, after the scale µ F = µ 0 is fixed for the LO contribution, the variation of the scale in the remaining NLO part does not noticeably change the predicted Drell-Yan cross section. In [39] it was shown how to calculate µ 0 for the Drell-Yan process and, moreover, that the NLO prediction with µ F = µ 0 is very close to the NNLO result.
Returning to the earlier discussion in this subsection, it is clear from (4) and (14)  We can also explain this pictorially. As noted above, the µ F dependence of the NLO correction is caused by the subtraction of the contribution generated by the evolution equation. In particular, the correction which (i) depends on µ F and (ii) is enhanced by ln(1/ξ), is generated by the laddertype diagrams 11 shown in Fig. 2. The choice of the factorisation scale µ F i determines which part of the diagram is attributed to the evolution of the PDF and which to the hard matrix element. By choosing the value of µ F = µ 0 = M V /2 we include all the ladder cells in the LO part, so that are proportional to 1/x. So these double log terms come from DGLAP integrals of the form where we have strong k 2 ≡ k 2 t and x ordering. Thanks to the strong k t and x ordering in these double log LO integrals, the correct upper limit, µ F , in the first integral automatically provides the exact resummation of all the terms in the double log series. 10 We stress again that, in the high energy (x ξ) contribution, the form of the integral over the gluon loop momentum l t is exactly the same in both the quark ( Fig. 1(b)) and gluon ( Fig. 1(b) but with the quark lines replaced by gluons) NLO contributions. Therefore the scale µ F = µ 0 simultaneously nullifies the high energy quark and gluon contributions. Note that this is only true after the corrections discussed in Section 2.1 are included. 11 Ladder diagrams occur in the axial gauge which is conventionally used to calculate the GPDs. the remaining NLO matrix element will not receive any contribution from ladder-type diagrams (at small ξ). Note that the ladder diagrams would contribute to the NLO matrix element (from µ F to µ 0 ) if we were to choose a scale µ F < µ 0 . For the choice µ F > µ 0 the NLO matrix element would contain the ladder diagrams with the negative sign in order to compensate the extra contributions (from µ 0 to µ F ) included in the LO part.
Since we consider an exclusive process, the mass of the final state is fixed, but the intermediate states, corresponding to the discontinuity which gives the imaginary part of the amplitude (depicted using a dotted line in Fig. 2 It should be emphasised that the asymptotics of the NLO amplitude is used only to determine the effective scale µ F . In all our further numerics we use the full expression for the complete NLO amplitude as given in [24,32,33].

Can exclusive vector meson production data be included in a global PDF fit?
In this section we present our results for the exclusive γp → J/ψ + p and γp → Υ + p processes as functions of the factorisation scale µ f , and the renormalisation scale µ R . The cross section relies sensitively on the gluon PDF at small x, in a kinematic regime where it is very poorly known, as well as depending of the choice of scales. For these quasi-elastic processes we present just the dominant imaginary part of the amplitude. The real part of the coefficient functions has been calculated exactly, both via a dispersion relation [24] and by directly computing the real part of the loop integrals [32,33], it is non-zero in the time-like region |x| < ξ. Therefore, to compute the real part of the amplitude directly, we need also the GPDs in this region. Unfortunately, the Shuvaev transform is not valid in the time-like region [23]. Nevertheless, the real part of the amplitude can be included via dispersive methods on the level of the amplitude. However, giving a prediction of the full cross section is not our objective here. Rather our aim is to study the stability of the perturbative predictions and to investigate whether or not we can determine optimum scales so that experimental data for these processes, and the related pp → p V p processes, can be included in global PDF (collinear) analyses to constrain the gluon PDF at low x, in a domain for which, at present, there are no data. For this goal it is sufficient to work with the more simple imaginary part of the amplitude.
We recall the main result of Section 2. We start from the key equation, (14), that is To obtain this result we had introduced by hand a new scale µ F in the LO term and showed that the choice µ F = m = M V /2 allowed us to resum and transfer all the enhanced large ln(1/ξ) terms from the NLO contribution to the LO term, leaving a much smaller remaining NLO contribution There is still a µ f scale dependence, but now this should be much weaker.

The process γp → J/ψ + p
To obtain predictions for exclusive J/ψ photoproduction we are working at scales close to the input GPDs used in the calculation. The results for the scale variation of the imaginary part of the amplitude are shown in the two plots of Fig. 3. In each plot we show separately the LO and NLO contributions to the amplitude. 12 The left plot shows how these contributions change if we vary all the scales µ f = µ F = µ R ≡ µ simultaneously, taking µ 2 = m 2 /2, m 2 , 2m 2 . In the right plot we fix µ F at the optimum scale µ F = m and vary µ f = µ R ≡ µ. The result is dramatic. The transfer of the ln(1/ξ) terms from the NLO to the LO contribution has significantly reduced the µ f scale dependence.
But there is another, more severe problem. The LO contribution is dominated by the NLO contribution of opposite sign. Thus the imaginary part of the quasi-elastic γp → J/ψ + p amplitude changes sign when the NLO contribution is added. As it stands, this result is in contradiction with modelling the interaction as an elastic forward scattering where the imaginary part of the amplitude is positive. What is happening? The explanation is interesting. In general, the global DGLAP PDF analyses start from input forms which are completely arbitrary and, moreover, neglecting any constraints, know nothing about the structure of the evolution at low Q 2 . It is then found that the gluon PDF tends to be small, or even negative, in the low Q 2 , low x (10 −4 < ∼ x < ∼ 10 −2 ) domain, see Fig. 4. Clearly due to the lack of data constraints in this domain the gluon is not reliably known.
Let us study the over-simplified case with an input gluon g(x) = 0 so that at the input scale we have only quark PDFs. For γp → J/ψ + p the quark contribution only appears at NLO. The imaginary part of this NLO contribution is negative with respect to the normal (g = 0) LO term. Indeed, when calculating the NLO coefficient function we must subtract the contribution which is already generated by LO evolution. However, the subtraction is performed purely formally in order to avoid the infrared divergence and does not account for the actual value of the gluon density used NNPDF3.0nlo Figure 4: The gluon distribution at Q 2 = 1.21 GeV 2 as determined by two recent global analyses [1,2]. Figure produced using [42,43].
to calculate the LO term. Therefore, at low scales, where the LO term is suppressed by much too small 'unphysical' gluon PDFs obtained from the global fits, the NLO term generated by the quark is not just a correction, but is the dominant contribution. As a consequence, the imaginary part of the quasi-elastic amplitude becomes negative (in comparison with the LO contribution) and the observed behaviour reveals a lack of stability of the perturbative series at this order. Of course, at high scales where evolution, and not the input, determine the PDFs, we obtain a sensible NLO prediction.
This simple consideration demonstrates that if γp → J/ψ + p and pp → p + J/ψ + p data could be included in a global fit, then they would put a strong constraint on the low x input gluon distribution. It is not simply that g(x) must be positive, but that actually the input gluons cannot be smaller than the density of gluons emitted by the quarks before the beginning of the evolution, when the parton virtuality Q < Q 0 . If J/ψ data would have been included in the collinear global parton analyses, this constraint on the low x input gluon would have been automatically satisfied.
Indeed, it is seen from the second plot of Fig. 3 that after the ln(1/ξ) enhanced corrections are resummed by fixing µ F = m, the stability of both the LO and the NLO components of the amplitude under the scale variations are much better. From this viewpoint the J/ψ data may be included in a global PDF analysis. The only problem is that the gluon density obtained in the existing PDF analyses is too small at input scales and low x (10 −5 < ∼ x < ∼ 10 −2 ), so that we get the wrong sign of the imaginary part of the amplitude. This means that if J/ψ data were included in the global PDF analyses, we must get a larger gluon density at low scales. 13 Larger gluons in this kinematic domain will increase the LO component of the J/ψ photoproduction amplitude and will provide the correct sign for the whole LO+NLO amplitude.

The processes γp → Υ + p and pp → pΥp
The top two plots of Fig. 5 show the results for the scale sensitivity of γp → Υ + p, analogous to those of Fig. 3 for γp → J/ψ + p. Again we see that fixing the scale µ F = M Υ /2 reduces the µ f scale uncertainty. The lower two plots show that part of this stability arises because the change caused by variation of µ R is to some extent compensated by an 'opposite' change due to the variation of µ f . If we were to choose a non-optimal scale, say for example, µ 2 F = 2m 2 , then the scale variation turns out to be about twice as large as for the optimal choice µ 2 F = m 2 . Due to the larger Υ mass we are now working at scales with more perturbative stability. The LO contribution is partly cancelled by the NLO term, but is not dominated by it. If we take the upper right or lower left plot of Fig. 5 then the µ f scale uncertainty of the amplitude is about ±15% and ±25% respectively. As a result the data for exclusive Υ production can be used in global PDF analyses to probe the gluon distribution down to x ∼ 10 −5 , in the case of LHCb kinematics. The calculation of exclusive pp → pΥp is described in [35]. It is based on the sum of the two diagrams shown in Fig. 6. For an Υ produced at large rapidity y, the dominant contribution is from the diagram with the larger γp centre-of-mass energy W + , which depends on the gluon density at x M Υ e −y / √ s. The small contribution from the other diagram, with much lower energy W − , may be estimated from the existing HERA data.
If we use the central values of the presently available PDF sets obtained from global analyses, then we find cross section predictions which are about a factor of four below the HERA exclusive Υ data [6,14,15]. This indicates that future PDF global analyses with Υ data included will, just as for J/ψ data, require a larger gluon distribution at low values of x. While the exclusive Υ process samples the PDFs at scales Q 2 ∼ M 2 Υ /4, the larger gluon required at these Q 2 will affect the gluon densities at all values of Q 2 via the evolution.

Note on the alternative k t -factorisation approach
The collinear DGLAP global analyses tend to result in gluon distributions with valence-like x distributions at low scales, see Fig. 4. On the other hand, the approach used in the JMRT paper [35], to study exclusive vector meson production, is free from this problem. There, the 'NLO' prediction was not calculated as a correction from NLO Feynman diagrams in collinear factorisation, but  Figure 5: Predictions of ImA/W 2 for γp → Υ + p as a function of the γp centre-of-mass energy W , produced using CTEQ6.6 partons [40]. The scales are  Figure 6: The two diagrams describing exclusive Υ production at the LHC. Diagram (a), the W + component, is the dominant contribution for an Υ produced at large rapidity y. Thus data for pp → pΥp allow a probe of very low x values, x M Υ e −y / √ s; recall that in the (dominant) imaginary part of the Born amplitude we have x = ξ. approximated by taking the full integral over the gluon k t in Fig. 1(a), including an ansatz for the now k t dependent gluon. The convergence of this explicit k t integral provides effectively the 'optimal' value of the factorisation scale and in this way sums all ln(1/ξ) enhanced terms, accounting for a large part of the NLO corrections. It does, of course, not include corrections from diagrams which do not have a ladder structure, however, the contribution from the light quarks, Fig. 1(b), is already included in the unintegrated gluon distribution used in JMRT [35]. Clearly, then the quasi-elastic amplitude is positive definite. Let us briefly describe how our 'NLO' predictions for γp → Υ + p data were made. First the incoming gluon distribution was fitted, within this approach and using an ansatz for the gluon based on the important double-logarithmic dependence on ln(1/x) and ln Q 2 , to reproduce the J/ψ data from HERA and the LHC. Then we proceeded to Υ production using our fitted gluon. We verified that, in the relevant kinematic domain, this gluon reproduces the NLO DGLAP evolution to good accuracy. Therefore it was not surprising that the 'NLO' predictions for the Υ cross section, as a function of W , agreed well with the LHCb data [18] when they became available.

Conclusions
Here we have been concerned about the inclusion of exclusive vector meson production data in global PDF analyses in order to probe the gluon density at small values of x; that is, in the x < ∼ 10 −4 domain. Note that in this domain the input gluon PDF is freely parametrized, and is found to have a tendency to be valence-like with large uncertainties. It was hoped that the situation would be changed as data for exclusive vector meson production, V = J/ψ and Υ, which are very sensitive to the gluon at small x, became available at HERA and the LHC; particularly measurements of the exclusive process pp → p V p at the LHC with V detected at large rapidity. Why has this not happened (so far)?
For J/ψ, the main problem is, at first sight, the very poor convergence of the perturbative expansion in the collinear approach. This can be seen, e.g., from Fig. 3, where the NLO contribution is comparable to the LO term and opposite in sign. However, we argued that this large NLO contribution, in comparison with that of the LO, reflects not the poor convergence, but rather that the global PDF analyses find a gluon density which is too small at low x and low scales. For this reason we find a LO contribution to exclusive J/ψ photoproduction amplitude which is too small. We argued that the input gluon in the collinear approach, used in the global PDF analyses, should not be parametrized freely, but should be subject to some constraints. We noted that working in the physical, k t -factorisation type scheme would avoid these problems.
For exclusive Υ production the situation is much better. We found that the optimum factorisation scale is much higher, the perturbative expansion at NLO in collinear factorisation converges well, and the remaining mild scale dependence of the predictions (∼ ±15% on amplitude level) means that data for pp → pΥp can now be included in the global PDF fits to determine the gluon in the low x regime for the first time.
It is appropriate to list the theoretical uncertainties of the present calculation. Although the leading double log terms have been resummed correctly to all orders, there still exist remaining NNLO and higher contributions, which are unknown at present. When the full NNLO amplitude is known we showed how the scale uncertainty can be further reduced. Next, we consider the accuracy of the expressions which relate GPDs to conventional PDFs, (8) and (9). These relations are based on conformal invariance and the equality of the Gegenbauer moments of GPDs to the Mellin moments of PDFs. Due to the polynominal property [21] the accuracy of this equality of the moments is O(ξ 2 ), providing an O(ξ) accuracy for (8) and (9), see [23]. However, recall that the Shuvaev transform assumes the absence of the singularities in the right-half Mellin-N plane for the input distribution. This assumption is physically reasonable since in the Regge approach there are no singularities in the right-half (j > 1) plane in the space-like (|x| > ξ) domain where we actually work. Nevertheless, whenever possible, this assumption should be checked. One check is that the predictions for the GPD/PDF ratio obtained in this way are in a good agreement (especially for low x gluons, which are the most important for exclusive J/ψ or Υ production) with the NLO results of [44] coming from a fit to deeply virtual Compton scattering (DVCS) HERA data. Finally, the relativistic correction to the vector meson wave function, which is discussed in Section 2.2, is not expected to be significant [34].
We noted that an alternative probe of the low x gluon at low scales has been considered in [29,30]. There they study the data for charm (and beauty) production obtained in the forward direction by the LHCb Collaboration [27,28]. Again it is found that the predictions strongly depend on the choice of factorization scale. Their prediction is lower than the data if the natural scale µ 2 F = m 2 Q + p 2 t is chosen, and it is found that a larger µ F is needed to reproduce the data using the existing global PDFs. Thus our expectation of a larger gluon at low x is not in contradiction with the LHCb charm (and beauty) forward data.