How to include exclusive $J/\psi$ production data in global PDF analyses

We compare the cross section for exclusive $J/\psi$ photoproduction calculated at NLO in the collinear factorization approach with HERA and LHCb data. Using the optimum scale formalism together with the subtraction of the low $k_t<Q_0$ contribution from the NLO coefficient function to avoid double counting we show that the existing global parton distribution functions (PDFs) are consistent with the data within their uncertainties. However, at low $x$ the uncertainties of the present global PDFs are large. On the other hand, the accuracy of the LHCb data are rather good. Therefore, these data provide the possibility to directly measure the gluon PDF over the very large interval of $x$, $10^{-6}<x<10^{-2}$, at a fixed low scale.


Introduction
The parton distributions of the proton are known with good precision from global parton analyses as long as the x values are not too low. Indeed for x > ∼ 10 −3 the results of the different groups [1,2,3] agree with each other quite well. However, the uncertainty in the parton distributions strongly increases as we go to lower values of x, especially at low scales. This simply reflects the fact that no experimental data are used to directly probe this region. 1 On the other hand there exist data on diffractive J/ψ production measured by LHCb [6] which sample much smaller values of x > ∼ 10 −6 at low scales. Moreover, this process is driven dominantly by the square of the gluon PDF, a parton for which little direct experimental information exists. However, these data are not included in the conventional global analyses for two reasons. First, the corresponding cross section is not described by the usual PDFs but by the more complicated generalised parton distributions (GPDs), see [7] for a review. Next, the NLO corrections are large and the results strongly depend on the choice of scale.
In the present paper we recall how these problems can be solved within the conventional collinear approach by using the Shuvaev transform [8], which at small x allows the calculation of the GPDs from the conventional integrated PDFs. Secondly, the strong scale dependence can be reduced by choosing a factorization scale which effectively resums the double logarithmic α s ln(µ 2 ) ln(1/x) terms (which are enhanced by the large values of ln(1/x) at small x) and transfers them into the incoming PDFs. Finally, and most importantly, to avoid double counting, we have to subtract the low transverse momentum, k t , contributions below the input scale Q 0 from the NLO coefficient functions, as these contributions are already included in the input PDFs. The subtraction is of the form of a power correction which, as expected, is large.
Previously, the LHCb data for forward ultraperipheral J/ψ production were successfully described in [9] using the k t factorization framework. However, the k t factorization approach does not include the complete set of NLO corrections. Thus this approach does not allow these J/ψ data to be included in the NLO global analyses based on the collinear factorization theorems. Our formalism is based on the conventional collinear framework and includes all NLO corrections. In Section 5 we show that three existing sets of PDFs (NNPDF3.0 [1], MMHT2014 [2], CT14 [3]) taken at the optimal scale mentioned above, and convoluted with the NLO coefficient functions from which the low k t < Q 0 contribution has been subtracted, give a satisfactory description of the diffractive J/ψ HERA data, but vastly different predictions in the region of the LHCb J/ψ data.
The plan of the paper is as follows. In Section 2 we give our notation. In Section 3 we explain how our approach can be used to probe the PDFs. In Section 4 we demonstrate the stability of the analysis with respect to variations of the remaining scale dependence. In 1 Besides its intrinsic value, there are at least two further reasons to be interested in the behaviour of the gluon PDF at very small x and low scales µ ∼ 1.5 GeV. First, recall that the distribution of gluons as x → 0 governs the high-energy asymptotics of the scattering amplitude. In particular, the gluon distribution at some relatively low scale can be used as the boundary condition for the BFKL equation. This boundary condition for BFKL is needed to account for the effects of confinement. As was shown in [4,5], such a boundary condition replaces the BFKL cut (in the complex momentum j-plane) by a series of Regge poles. At very low x the boundary condition should indicate the presence of saturation effects that are needed to stop the power growth of the original BFKL amplitude. Another motivation for obtaining a reliable gluon PDF at small x is that it may be used to evaluate the production cross section of a possible new light particle at the LHC (if such a new particle exists) or to put a limit on the corresponding coupling.
Section 5 we show that the PDFs given by the existing global analyses agree with the J/ψ exclusive photoproduction data measured at HERA [10] and that they can be constrained at even smaller x ∼ 10 −6 using LHCb ultraperipheral J/ψ data. We discuss our results in Section 6 and present our conclusions in Section 7.

Notation & collinear factorization
The exclusive J/ψ photoproduction amplitude may be written, using collinear factorization, in the form [11] where we have suppressed the dependence on the renormalization and factorization scales, µ 2 R , µ 2 F , and on the invariant transferred momentum t. Here, the non-relativistic QCD (NRQCD) matrix element O 1 V describes the formation of the J/ψ meson with m c the charm quark mass. The quark singlet and gluon GPDs are denoted F q and F g , respectively. The quark and gluon coefficient functions C q and C g are known at NLO [11] and are given at tree level by The kinematics of the process are displayed in Fig. 1. The partons carry momentum fractions (X + ξ) and (X − ξ) of the plus-component of the mean incoming/outgoing proton momenta P = (p + p )/2. The photon-proton centre of mass energy squared is given by W 2 = (q + p) 2 , where q is the photon momentum. The asymmetry between the momentum fractions carried by the partons is parametrised by the skewness parameter, Due to the vanishing of the quark coefficient function at LO the process is predominantly sensitive to the gluon GPD. At LO, the gluon coefficient function is strongly peaked for |X| ∼ ξ and so the gluon GPD is probed close to F g (ξ, ξ). In fact, for the imaginary part of the amplitude, the LO gluon coefficient function acts as a Dirac delta function and the GPD is probed at exactly |X| = ξ. 3 Connecting exclusive production to the PDFs Firstly, let us recall the advantage of using the exclusive J/ψ LHCb data in global parton analyses in the collinear factorization scheme. It offers the possibility to probe PDFs (mainly the gluon PDF) at extremely low x in a so far unexplored kinematic regime. In particular, for forward ultraperipheral production, pp → p + J/ψ + p, the LHCb experiment can reach 2 for √ s = 13 TeV and rapidity Y = 4.5. Moreover the cross section is proportional to the square of the parton density, so the uncertainty on the PDF is reduced.
However, as mentioned in the introduction, there appear to be two disadvantages. First, the description of the exclusive J/ψ process depends on Generalized Parton Distributions (GPDs), and, second, there is a strong dependence on the choice of scale, indicating a large theoretical uncertainty. Immediately below we note how the first disadvantage is overcome. Then, in the next section, we discuss the removal of the sensitivity to the scale dependence.
Though exclusive J/ψ production is described by GPDs, at very low values of x and small momentum transfer t the GPD can be related to the conventional integrated PDF, via the Shuvaev transform, with accuracy O(x) [8]. The Shuvaev transform makes use of the fact that as ξ → 0 (and at t = 0) the Gegenbauer moments of the GPD become equal to the known Mellin moments of the PDF. Due to the polynomial condition (see e.g. [12]) even for ξ = 0 the Gegenbauer moments can be obtained from the Mellin moments to O(ξ) accuracy. Thus it is possible to obtain the full GPD function at small ξ from its known moments. Based on this fact we can obtain an expression which transforms the low x PDF to the corresponding GPD [8].
The GPD function (denoted by F a (X, ξ) with a = g, q in Fig. 1) accounts for the fact that the momenta of the 'left' and 'right' partons in the diagrams of Fig. 1 are different. In particular, they carry proton momentum fractions X + ξ and X − ξ respectively. The Shuvaev transform relates the GPD F a (X, ξ) to the PDF(X +ξ). We systematically construct GPD grids from a three-dimensional parameter space in X, ξ/X and Q 2 with forward PDF grids taken from the LHAPDF interface [13]. It turns out that the values of X that are most relevant in the convolution of the GPD with the coefficient function are of the order of ξ. Thus, indeed, in this way we probe the gluon PDF at values of x close to 2ξ.
Strictly speaking, by using such a transform we assume that the amplitude has no additional singularities in the right half (j > 1) of the complex angular momentum j plane. This assumption is well motivated physically, and moreover it was shown [14] that the results agree with those obtained [15] in an independent global GPD analyses of the available data.

Overcoming the strong scale dependence
The strong sensitivity to the choice of scale in the predictions for diffractive J/ψ photoproduction was first observed in [11,16] and recently confirmed in [17]. There are two sources for this sensitivity to the scale choice. Firstly, there is the double logarithmic contribution which contains a large ln(1/x) factor. For the region of interest, x ∼ 10 −5 , this means an order of magnitude enhancement. Secondly, there is double counting in the coefficient functions for Q 2 < Q 2 0 . We discuss how these problems are overcome in turn.

Treatment of double log contributions
It was shown in Ref. [16] that it is possible to find a scale (namely µ F = M ψ /2) which effectively resums all the double logarithmic corrections enhanced by large values of ln(1/ξ) into the gluon and quark PDFs, where ξ is the skewedness parameter of the Generalised Parton Distributions (GPDs). In terms of the usual (unskewed) PDFs related to GPDs via the Shuvaev transform, x 2ξ. That is, it is possible to take the (α S ln(1/ξ)ln(µ 2 F )) term from the NLO gluon (and quark) coefficient functions and to move it to the LO GPDs. This allows a resummation of all the double logarithmic, i.e (α S ln(1/ξ)ln(µ 2 F )) n , terms in the LO contribution by choosing the factorization scale to be µ F = M ψ /2. The details are given in Ref. [16], see also Ref. [18].  The result is that the γp → J/ψ p amplitudes, taken at factorization scale µ f , are schematically of the form With the choice µ F = µ 0 = M ψ /2, the remaining NLO coefficient function, C NLO rem (µ F ), does not contain terms enhanced by ln(1/x) ln(1/ξ).
Thus to summarize, eq. (4) allows us to consider different factorization scales µ f . However the scale in the first term on the right-hand-side is fixed to be µ F = m c independent of the value of µ f . Since the contribution from the second term is small we predominantly probe the gluon distribution at scale µ F = µ 0 .
Moreover, we find that after the scale µ F in (4) is fixed to µ F = µ 0 , the result (shown in the left panel in Fig. 2) becomes more stable with respect to variations of the factorization scale µ f in comparison to the huge variations seen in [11]. However note that the NLO correction is still comparable to the LO term and opposite in sign. As we discuss in Section 4.2, this is due to double counting between the NLO coefficient function and the contribution coming from DGLAP evolution. Once we avoid this double counting, we will see that the perturbative treatment is brought under control and also that we have a further reduction of the scale sensitivity.

Alternative resummation
As an aside, we note that in [20], instead of fixing µ F = µ 0 , it was proposed to resum the BFKL corrections, like α s ln(1/x), already in the coefficient function. It was stated that this would allow good scale stability to be obtained. The problem, however, is that then the coefficient function absorbs all the BFKL effects, and therefore the convolution of the GPD with the coefficient function is such that the coefficient function, C a (X, ξ), occupies almost the whole available ln(1/X) interval; that is the dominant contribution comes from X ∼ O(1) and not X ∼ ξ. Thus, we would lose the main advantage of probing the unexplored very small x regime.

Treatment of double counting power corrections
Next we consider a power correction which may further reduce the NLO contribution and, moreover, may reduce the sensitivity to the choice of scale. The correction is O(Q 2 0 /M 2 ψ ) where Q 0 denotes the input scale in the parton evolution which turns out to be important for the relatively light charm quark, m c M ψ /2. Let us explain the origin of this 'Q 0 correction' following Ref. [21]. We begin with the collinear factorization approach at LO. Here, we never consider parton distributions at low virtualities, that is for Q 2 < Q 2 0 . We start the PDF evolution from some phenomenological PDF input at Q 2 = Q 2 0 . In other words, the contribution from |l 2 | < Q 2 0 of Fig. 1(b) (which can be considered as the LO diagram, Fig. 1(a), supplemented by one step of the DGLAP evolution from quark to gluon, P gq ) is already included in the input gluon GPD at Q 0 . That is, to avoid double counting, we must exclude from the NLO diagram the contribution coming from virtualities less than Q 2 0 . At large scales, Q 2 Q 2 0 this doublecounting correction will give small power suppressed terms of O(Q 2 0 /Q 2 ), since there is no infrared divergence in the corresponding integrals. On the other hand, with Q 0 ∼ 1 GeV and µ F = m c (∼ M ψ /2), a correction of O(Q 2 0 /m 2 c ) may be crucial. Beyond NLO single logarithmic terms, ln(1/x), may again be present in the amplitude. However, we anticipate that including the Q 0 subtraction their impact will be much smaller.
In the present paper we use the NLO correction C NLO rem for J/ψ photoproduction excluding the contribution coming from the low virtuality domain 3 (< Q 2 0 ). We find that for J/ψ this procedure substantially reduces the resulting NLO contribution and, moreover, reduces the scale dependence of the predictions. It indicates the stability of the perturbative series. Indeed, as shown in the left panel of Fig. 2, before the Q 0 subtraction the NLO corrections may exceed the value of the LO contribution and, depending on the scale, even the sign of the amplitude can change. However, after the subtraction and choosing the optimal scale µ F = M ψ /2 in the leading order part of the amplitude (first term of (4)), we observe a rather good scale stability as shown in the right panel of Fig. 2.
In Fig. 3 we show the results for ImA a with a = g, q for the choice µ F = M ψ /2 = m c for two values of the factorization scale: µ 2 f = m 2 c and µ 2 f = 2m 2 c . We take µ R = µ f . Here A a=g,q are the gluon and quark contributions to the γp → J/ψ + p amplitude in the collinear factorization scheme at NLO. The plot shows the stability of the amplitude with respect to variations of µ f , and also that the Q 0 subtraction practically fully absorbs the quark contribution. With this set-up, we can therefore say that low x exclusive J/ψ photoproduction probes predominantly only the gluon distribution.

Renormalization scale
The renormalization scale is taken to be µ R = µ f . The reasons for this are as follows. Firstly, this corresponds to the BLM prescription [22]; such a choice eliminates the contribution proportional to β 0 (i.e. the term β 0 ln(µ 2 R /µ 2 f ) from the NLO terms in eq. (3.95) of [11]). Secondly, following the discussion in [23] for the analogous QED case, we note that the new quark loop insertion into the gluon propagator appears twice in the calculation. The part with scales µ < µ f is generated by the virtual component (∝ δ (1 − z)) of the LO splitting during DGLAP evolution, while the part with scales µ > µ R accounts for the running α s behaviour obtained after the regularization of the ultraviolet divergence. In order not to miss some contribution and/or to avoid double counting we take the renormalization scale equal to the factorization scale, µ R = µ f .

Description of J/ψ photoproduction data
All of the calculations presented so far are performed for the imaginary part of the production amplitude. The real part is obtained via a dispersion relation, which in the high energy limit (for the even signature amplitude) can be written in the simplified form [24] ReA ImA = tan π 2 Next we have used the nonrelativistic J/ψ wave function. As was shown by Hoodhboy [25], this provides an accuracy of a few percent, which is sufficient for our purposes. Actually, we calculate the value of ImA at t = 0 and then restore the total γp → J/ψ + p cross section assuming an exponential t behaviour with a slope B = 4.9 + 4α P ln(W/W 0 ) GeV −2 with W 0 = 90 GeV and α P = 0.06 GeV −2 . This parametrisation grows more slowly with W than the formula used by H1 [26], but is still compatible with the HERA data. We have chosen the slope parameter α P to be compatible with Model 4 of [27] which fits a wider variety of data.

HERA data
As can be seen from Fig. 4, the J/ψ photoproduction data obtained at HERA [10] are described reasonably well by all three sets of global partons [1, 2, 3] within our collinear approach. These data sample x values in the interval 4 In our approach we are free to choose the starting scale Q 0 and the µ F scale in the NLO correction. We work at LO in NRQCD and the description used for the results shown in Fig. 4 corresponds to the choices Recall that the choice µ F = m c provides the complete summation of the double log terms [16]. Besides giving a good description of the HERA data, the above choice of Q 0 and µ F give a stable theoretical prediction also when the scales µ f and µ R are varied, see Figs. 3 and 4.  Figure 4: The γp → J/ψ + p data obtained at HERA [10] and LHCb [6] compared with the predictions obtained using the PDFs taken from three different sets of global partons [1,2,3] with µ f = m c (solid lines). The dashed line for the CT14 prediction, corresponding to µ 2 f = 2m 2 c , is added to demonstrate the scale stability of our NLO predictions; but note that our optimal choice µ 2 f = m 2 c agrees better with the HERA data.

LHCb data
The LHC experiments do not directly measure the cross section of photoproduction, but instead the exclusive pp → p + J/ψ + p reaction. At small transverse momentum of the J/ψ meson this process is described by the two diagrams shown in Fig. 5. The photon can be emitted either by the upper or by the lower beam protons. Since the photon's transverse momentum, q T , is much smaller than that transferred through the strong interaction amplitude (shown by the double vertical lines in Fig. 5) the interference between these two diagrams is negligible.
The contribution corresponding to the right graph, with a smaller photon-proton energy W − , comes from relatively large x, and can be subtracted using the description of HERA data. Thus the cross section for J/ψ photoproduction at the large energy, W + , may be extracted from the LHC measurements.
The last point is that in dealing with proton-proton interactions we must account for the possibility of an additional soft interaction between the two colliding protons. This interaction will generate new secondaries which will populate the rapidity gap and destroy the exclusivity of the event. The probability to have no such additional interaction is called the gap survival probability S 2 < 1. The value of S 2 depends on the pp collider energy and the partonic energy x . The q T of the photon is very small and so the photon can be considered as a real on-mass-shell particle.
W . The values of S 2 (W ) as a function of W were calculated using the eikonal model [28] which well describes the data on differential dσ(pp)/dt cross section and low mass diffractive dissociation. The details of the procedure to extract σ(γp → J/ψ + p) at large W + energies is described in reference [9]. Actually, in our figures we plot the low x LHCb data points obtained in this way and presented in [6].

Discussion of the results
The theoretical predictions, obtained by using the approach described above, are presented in Fig. 4. There we compare our predictions for the cross section for J/ψ photoproduction obtained using three different sets of global partons [1,2,3] with the HERA and LHCb data. The curves correspond to using the central values of the global PDFs. At the lower energy of the HERA data, where the global gluon PDF uncertainty is not too large, the predictions agree with the experimental values reasonably well. In the kinematic region covered by the LHCb experiment the present global PDF analyses do not sample any data, and hence they have almost no predictive power in this low x regime.
On the other hand, as is seen from Fig. 4, by exploiting the LHCb data for exclusive J/ψ production we have the possibility to greatly improve our knowledge of the gluon PDF down to x ∼ 3 × 10 −6 . The GPD(X, ξ) obtained via the Shuvaev transform is driven dominantly by the value of x = X + ξ 2ξ, while x = X − ξ x is small. Recall that in the LO contribution (given by the first term of eq. (4)) we sample the gluon PDF at x = X + ξ = 2ξ while in the NLO contribution (the second term) the momentum fraction carried by the gluon may be larger. As a check we have calculated the median value, med(X), of the corresponding X, defined in such a way that X > med(X) gives 0.5 of the NLO contribution. In the convolution of the coefficient function with the GPD (see eq. (5)) the X distribution is sharply peaked at X ξ for the gluon contribution while for the quark NLO contribution the value of med(X) 1.18 ξ. However, as it is seen from Fig. 3, the quark term is practically negligible. Thus we can say that the exclusive J/ψ production indeed probes the gluons at x = X + ξ 2ξ.

Conclusions
We show that the J/ψ meson photoproduction process and exclusive J/ψ production, pp → p + J/ψ + p, at the LHC, can be consistently described in the collinear factorization framework at NLO. The choice of the optimal scale µ F = µ 0 = M ψ /2 effectively resums the large double logarithmic terms, i.e (α s ln µ 2 F ln(1/ξ)) n . This, together with the Q 0 subtraction (needed to avoid double counting between the NLO coefficient function and the DGLAP input PDFs), leads to a largely improved scale stability of the theoretical prediction. In other words, this framework overcomes the extremely large scale uncertainties found in the existing NLO predictions [11,16,17] of diffractive J/ψ photoproduction in the collinear factorization approach. It is not surprising that at these low scales the power correction arising from the Q 0 subtraction is crucial. Another power correction coming from absorptive effects should reveal itself as the saturation of the gluon density. At the moment this is not noticeable; for small x the data appear to be compatible with the gluon PDF parametrization xg ∝ x −λ .
Huge uncertainties in the low x gluon PDF found in the existing global PDF analyses reflect the fact that no corresponding low x data were included in the fitting procedure. This is well illustrated in Fig. 6 which shows the prediction of, for example, the NNPDF3.0 [1] parton set together with its 1σ error band. However, using the proposed approach the good accuracy of the exclusive J/ψ cross section presented by LHCb will allow the determination of the NLO gluon PDF down to x ∼ 3 × 10 −6 , and the HERA data will improve the determination of the gluon for 10 −4 < ∼ x < ∼ 10 −3 .