Very low $x$ gluon density determined by LHCb exclusive $J/\psi$ data

The low $x$ behaviour of the gluon density $xg(x,\mu^2)$ at scale $\mu^2=2.4$ GeV$^2$ is determined using exclusive $J/\psi$ production data from HERA and LHCb within the framework of collinear factorisation at next-to-leading order (NLO). It is shown that in the interval $3\times 10^{-6}<x<10^{-3}$ the gluon distribution function grows as $xg(x,\mu^2)\propto x^{-\lambda}$ with $\lambda=0.135\pm 0.006 $. The impact this experimental data will have for the global parton distribution function (PDF) analyses in this low $x$ domain is quantified. No indication in favour of parton density saturation is observed.


Introduction
At moderate values of x the parton distribution functions (PDFs) of the proton are determined from global data with good accuracy by PDF analyses, see, for example, [1,2,3]. However in the low x domain, with x < 10 −3 , there are practically no data to constrain the input parton densities. In this domain the gluon PDF dominates. The global PDF predictions are based simply on extrapolation using some more-or-less arbitrary ansatz for the input distributions. For this reason the uncertainties of the PDFs in the very low x domain are huge.
There are two types of data which probe this domain. The first is open charm production and the second is exclusive J/ψ production -both processes have been measured by the LHCb collaboration in the forward region. These processes probe mainly the gluon PDF at a rather low scale O(m c ) close to the input Q 0 values of the global PDF analyses. Here m c is the mass of the charm quark. Open charm production is experimentally more complicated to measure as it is extracted from D-meson production data, but the theoretical formalism is direct. On the other hand exclusive J/ψ production is experimentally much cleaner, but the theoretical formalism needs care.
The charm and D-meson data [4] were used to restrict the uncertainty of the NNPDF gluon PDF in the low x region in [5,6,7,8,9]. However there are some inconsistencies in the energy and rapidity behaviour of the experimental results, which were discussed in [10,11]. The exclusive J/ψ data are more consistent and have better accuracy than the inclusive D-meson cross section.
Nevertheless, until now the J/ψ data have not been used in global analyses due to theoretical complications. First, the J/ψ cross section is driven by Generalised Parton Distributions (GPDs) and not directly by the conventional collinear PDFs. Second, the first calculations of the corresponding NLO coefficient functions revealed a huge scale uncertainty of the predictions [12,13,14]. Both of these problems have been overcome.
The amplitude, A, for exclusive J/ψ photoproduction can be expressed, within collinear factorisation, as where the F i and C i are, respectively, the GPDs and coefficient functions (known at NLO, see [12]) and their dependence on the renormalisation and factorisation scales, µ 2 R , µ 2 F , has been suppressed, as well as that due to the invariant momentum transfer, t. The kinematics of our set-up are shown in Fig. 1.
It was shown that in the relevant region of very small skewedness parameter ξ, depicted in Fig. 1, the GPD functions can be related to the normal PDFs via the Shuvaev transform [15]. This relation is based on the fact that, due to the polynomial condition, the Gegenbauer moments of the GPD are equal to the known Mellin moments of the non-skewed PDF up to O(ξ) accuracy at NLO [15,16]. The GPD grids are generated from PDFs supplied on grids via the LHAPDF interface [17]. Here M ψ is the mass of the J/ψ and W 2 is the photonproton energy squared for the γp → J/ψ p subprocess.
The second problem, that concerns the strong dependence on the factorization scale observed in the low x region, was essentially removed by subtracting the low k t < Q 0 contribution from the NLO coefficient functions. This subtraction is needed to avoid the double counting between the NLO coefficient function and the contribution hidden in the input PDF [14]. In addition the Figure 1: LO (left panel) and NLO (right panel) contributions to γp → V + p, where V = J/ψ (or Υ). Here the momentum P ≡ (p + p )/2, with k the loop momentum. Note that the momentum fractions of the left and right input partons are x = X + ξ and x = X − ξ respectively; for the gluons coupled directly to the on-shell heavy quark pair, we have x x and so x 2ξ.
double log terms, [α s ln µ 2 F ln(1/x)] m , can be resummed in the leading order term by choosing the optimum scale µ F = m c = M ψ /2 for our process [18]. The NLO amplitude A(µ f ), with factorisation scale µ f , can be written schematically in the form With the choice µ F = M ψ /2, the remaining NLO coefficient function, C NLO rem (µ F ), does not contain terms enhanced by ln(1/x) ln(1/ξ).
The approach was described in more detail in [19] where it was shown that the HERA data on diffractive J/ψ photoproduction [20] with energies corresponding to x > 10 −3 are well described using the present global gluons. 1 This demonstrated the efficiency of the method, which will be used in the present note to extract the behaviour of the gluon in the low x region (x < 10 −3 ) from the exclusive J/ψ LHCb data [23] (as well as HERA photoproduction data that lie in this region). 1 We should also mention the possibility of relativistic corrections to the NRQCD matrix element that we use in our approach. Recall that, strictly speaking, if we were to include relativistic corrections, see for example [21], then we must simultaneously account for the higher, cc + g, Fock component of the J/ψ wave function. As was shown in [22], these two corrections largely cancel each other, leading to a final correction of the order of a few percent provided that the NRQCD matrix element is normalized to the leptonic decay width, J/ψ → l + l − , and the charm quark mass is chosen to be m c = M ψ /2, as is kept in the present paper. Note also that the correction to the NRQCD matrix element changes the normalization of the J/ψ cross section but does not affect the x (or W ) behaviour of the low-x gluon. The fact that at x > ∼ 0.001 the data are well described by the existing global gluons is an argument in favour of the correct normalization, that is, in favour of small relativistic corrections to our approach.
As was shown in [19], after the k t < Q 0 subtraction the quark contribution to this process is negligibly small in this x region. Thus we determine just the gluon PDF and use the quark PDF from the existing global fits.
Of course, at the moment, global PDF analyses are performed to NNLO accuracy. However, as a first step, we start fitting the J/ψ data at NLO. In the future this approach can be extended to NNLO. 2 The outline of the paper is as follows. In Section 2 we describe the ansatz that we will use to parametrize the NLO gluon PDF in the collinear factorization scheme in the low x domain, x < 0.001. In Section 3, after a brief discussion of the exclusive J/ψ data, we describe how we determine the low x gluon directly from the data. In Section 4, we compare the results we find for the low x gluon with those obtained by reweighting the NNPDF gluon using the D-meson LHCb data. Finally, in Section 5, we provide a reweighting of the NNPDF3.0 gluon via the exclusive J/ψ data and compare and contrast this with the gluon obtained from the above alternative approaches. Our conclusions are briefly summarized in Section 6.

Ansatz for the low x gluon
It was demonstrated in [19] that the diffractive J/ψ cross section is driven by the Generalised Parton Distributions, GPD(X + ξ, X − ξ), of the gluon with X ξ, see Fig. 1. That is, to describe the LHCb data, we effectively need the gluon in the region of low x X + ξ only. So it is sufficient to parametrize the gluon in the region x < 10 −3 . On the other hand the Shuvaev transform, that relates the GPD to the conventional collinear gluon PDF, includes an integral over the whole x < 1 interval. Moreover, the transform was derived assuming that the gluon had a smooth analytical behaviour with the property that g(x) → 0 as x → 1. In order to satisfy these requirements we choose the following ansatz for the conventional gluon PDF, and where xg global is the value of the gluon PDF obtained in a global PDF analysis. The simplest low x form for the gluon would be where the normalization factor N 0 is chosen so that for n = 1 the gluon PDF has the matching at x = x 0 , The factor n in (6) is close to 1. It allows the possibility of matching to a global gluon whose normalization differs from N 0 but still lies within the global gluon error band at x = x 0 . The factor (1 − x) in (6) provides the vanishing xg → 0 as x → 1. Note that due to the smooth form of C in (5) the complete distribution (4) does not violate analyticity even for n = 1.
Alternatively, in order to compare our present collinear determination of xg new with an earlier determination of the low x gluon obtained in the k t factorization approach [24], we also use the ansatz where the parameter a now plays the role of λ. Here, with three light quarks (N f = 3) and N c = 3 we have β 0 = 9. We take Λ QCD = 200 MeV and q 2 0 = 1 GeV 2 , as in [24], with µ 2 0 = 2.4 GeV 2 fixed. The exponent in (8) resums, to all orders in m, the double logarithmic terms (α s ln(1/x) ln µ 2 ) m and hence we find that, to good accuracy, we reproduce the NLO DGLAP low x evolution in the interval of Q 2 from 2 to about 30 GeV 2 . Therefore this parametrization can be used to describe Υ photoproduction data as well.
3 Determination of the low x gluon from J/ψ data Here, we show the results of our fits to J/ψ photoproduction data for x < 10 −3 , using an ansatz for the gluon PDF as described in eqns. (4)- (7). The matching is made at x 0 = 10 −3 using the gluon PDF from three NLO parton global analyses, NNPDF3.0 [1], MMHT14 [2] and CT14 [3]. Due to the small contribution of the quark sector at NLO to the J/ψ cross section [19], we do not attempt to fit the quark PDFs but only the gluon PDF around its input scale. The quark PDFs obtained in the global NLO analyses are therefore used for all x.

The exclusive J/ψ data from LHCb
The LHCb experiment, by design, does not directly measure the cross section for J/ψ photoproduction but instead that for exclusive pp → p + J/ψ + p [23]. The experiment is unable to tag forward protons accompanying the J/ψ so instead only the rapidity of the J/ψ is measured. Events are selected by ensuring a large rapidity gap on both sides of the J/ψ measurements, where the transverse momentum of the J/ψ is small, and assumed to correspond to exclusive reactions. The lack of forward proton tagging means it is also not possible to determine which of the two protons emitted the photon. The ultraperipheral amplitude for a given J/ψ rapidity is then generally the sum of two photoproduction amplitudes with different W 2 , depending Figure 2: Two leading-order (LO) diagrams describing exclusive J/ψ production at the LHC. The left diagram, the W + component, is the major contribution to the pp → p + J/ψ + p cross section for a J/ψ produced at large rapidity Y . Thus such data allow a probe of very low x values, The q T of the photon is very small and so the photon can be considered as a real on-mass-shell particle.
on which proton emitted the photon and which was the target, see Fig. 2. The interference contribution is suppressed as the photons transverse momentum, q T , is much smaller than that of the proton exchanging the gluons. 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 existing description of HERA data. The cross section for J/ψ photoproduction at the large energy, W + , may therefore be extracted from the LHCb measurements.
Additionally, at the LHC, there is a non-negligible probability of additional soft interactions between the two colliding protons that can result in secondary particles polluting the rapidity gaps used to select the exclusive events. This will suppress the number of events deemed exclusive and therefore one must account for the gap survival probability, S 2 < 1, to have no such additional interaction. The value of S 2 depends on the pp collider energy and the partonic energy W . The values of S 2 (W ) as a function of W were calculated using the eikonal model [25] which well describes the data for the differential dσ(pp)/dt cross section and lowmass diffractive dissociation. The details of the procedure to extract σ(γp → J/ψ + p) at large W + energies is described in [24]. We use the low x LHCb "data" points obtained in this way by the LHCb collaboration [23].

Description of the J/ψ data
To set the scene, we first use eq. (1) at LO and NLO to generate and compare cross section predictions using the existing LO and NLO partons from [1,2,3], respectively, for the x-range where we have used exclusive J/ψ data from H1, ZEUS and LHCb. In this way, we are able to quantify the scale dependence of the theoretical prediction as well as the size of the NLO result relative to the LO one. In Fig. 3, we show such a comparison using CT14 partons [3]. Our choice of scales is explained in [18]. The NLO scale variation is smaller than that at LO and a better description of the HERA data is obtained with the NLO result. The plot emphasises that in the region where the current PDFs are well constrained, it is still crucial to use the NLO description. It is reassuring and non-trivial that our NLO prediction, with the 'optimum' scale choice, agrees well with the HERA data.
We now determine the low-x gluon by performing a two-parameter (λ and n, as defined in eq. (6)) fit of all the σ(γp → J/ψ + p) LHCb and HERA data with x < 0.001 using, as input, NLO parton PDFs from [1,2,3]. The results are shown in Table 1 Table 1: The values of λ and n obtained from fits to the J/ψ data using three sets of global partons.
The respective values of the total χ 2 min (and χ 2 min /d.o.f) for 45 data points are also shown.
The respective values of the χ 2 min statistic were calculated accounting for the bin-to-bin correlated errors within each individual experimental data set as well as uncorrelated errors. The covariance matrix was constructed, and iterated, according to the 't 0 prescription' as outlined in [26]. We use all HERA data points [20] with W > 100 GeV and all LHCb [23] data points. [20] we allow for a fully correlated 6.5% normalisation error. For the H1 2006 data set [20] we include a fully correlated 5% normalisation error while for the H1 2013 data set [20] we use the full covariance matrix as provided by H1. For the LHCb 2014 data [23] we allow for a fully correlated ∼ 7% normalisation error. Finally, for the LHCb 2018 data [23], we use the covariance matrices supplied by the collaboration as well as a fully correlated normalisation error of ∼ 4%.

For the ZEUS 2002 and 2004 data sets
The description of the exclusive J/ψ cross section is shown in Fig. 4, while the gluons extracted from the J/ψ data at µ 2 = 2.4 GeV 2 and x < 0.001 are shown in Fig. 5. The error bands are obtained by sampling over the two parameters within their individual 1σ standard deviations, accounting for their correlation. The hatched green band in Fig. 5 in addition accounts for the uncertainty due to the choice of the global (NNPDF3.0, MMHT2014 or CT14) partons. The gluon at very small x shows no hint of the onset of saturation.
Starting from three different sets of global partons, we obtain practically the same low x gluons with the same quality (χ 2 min ) of the description. The typical errors are ±2.5% for the normalization and ±4% for λ. We see from gluon density provides an excellent description of the J/ψ data in the fitted x < 10 −3 region, irrespective of which global parton set is used. In fact, the three descriptions only visibly differ for x < 10 −5 . Note that the observed hierarchy of central cross section predictions at x ∼ 3 × 10 −6 differs from that expected given the power behaviours in Table 1. We have checked that this is due to the small x and small scale quark behaviour of the global sets. In the left hand side of Fig. 6 we compare the uncertainties of the gluon densities given at x = 0.001 and µ 2 = 2.4 GeV 2 by the global analyses, while in the right hand side we show the values that are obtained after fitting the J/ψ data. The J/ψ data are seen to greatly improve the knowledge of the gluon in the low x interval 3 × 10 −6 < x < 10 −3 . In particular, we find at x 0 = 0.001 that x 0 g(x 0 , µ 2 = 2.4 GeV 2 ) = 2.28 ± 0.06.

The alternative double-log parametrization
While the simple two parameter ansatz in (6) leads to a very good description of the J/ψ data, it is still informative to repeat the procedure using the double-log ansatz in (8). Recall that a similar form was used in [24]. The result obtained using the NNPDF3.0 NLO parton set is a = −0.046 ± 0.006, n = 0.979 ± 0.025, χ 2 min /d.o.f = 1.05. The description and the behaviour of the low x gluon are very similar to that obtained using (6). We find that the fit using the double log parametrization gives the central value x 0 g(x 0 , µ 2 = 2.4 GeV 2 ) = 2.31 in agreement with (10).
Note that the double-log parametrization gives a result close to that obtained in the k tfactorization approach [24]. However now, accounting for the complete set of NLO corrections, we find that the gluon growth with energy (1/x) is less steep than that obtained in [24]. Instead of a = −0.10 we now have a ∼ −0.05. The LHCb16 data used in [24] have been replaced by the data in [23] that is used here, but this is not accountable for the difference in a.

Is there evidence of saturation from exclusive J/ψ data?
High energy exclusive J/ψ production was recently described in [27] based on a BFKL approach. The authors claim that "there are strong hints for the presence of the saturation effects in exclusive photo-production of J/ψ at small x". We have to emphasize that actually the authors of [27] refer to absorptive corrections rather than saturation. Indeed, saturation means that the gluon density tends to a constant value, xg(x, µ 2 ) → const as x → 0 and at a fixed scale µ [28]. That is, the power λ in (6) behaves as λ → 0. A first hint of saturation would be to observe that the power λ (measured in some small-x interval) starts to decrease with decreasing x. The data, as shown in Fig. 3, do not indicate such behaviour.
What is actually shown in [27] is that the LO BFKL intercept, α BFKL = 1+ω 0 = 1+λ is too large to describe the high energy J/ψ data and that absorptive corrections (which are included into the non-linear BK [29] equation) are needed to tame the growth of the gluon density (6), that is to decrease the value of λ.
It is well known that the LO BFKL intercept is too large. It becomes smaller in the next-toleading (NLL) approximation. Indeed, it is seen from [27] (the short dashed green curve of their Fig. 1) that the HSS gluons [30], based on the NLO BFKL linear equation, are in agreement with the exclusive J/ψ data.
Therefore the growth of the gluon density with a smaller but non-zero λ is not evidence for 'saturation'. At the moment no hint of saturation is observed in exclusive J/ψ data at the scale µ 2 = 2.4 GeV 2 and x down to 10 −5 .

Comparison with low x gluons from D-meson data
As mentioned in the introduction, it is also possible to determine the low x gluon density from the data for various modes of inclusive open charm production of D-mesons and their excited states. In this section, we provide a comparison of the results obtained from the data for inclusive D-meson production and exclusive J/ψ production.
Inclusive D-meson production data via pp collisions at the LHC are available at centre of mass energies 5, 7 and 13 TeV [4]. The kinematics of the different modes of production of the D-mesons allow for a coverage down to x ∼ few × 10 −6 . In [5] the authors studied the impact these data for {D 0 , D + , D + s } final states would have on the small x NLO gluon within the NNPDF3.0 global analysis through a Bayesian reweighting. While the corresponding NLO calculation for D-meson production suffers from large theory uncertainties attributed to the dependence on the factorization scale and large higher order corrections, construction of ratios of the double-differential cross section in rapidity and transverse momentum bins provides a means to combat this residual scale dependence and thereby quantitatively assess the impact the data would have in the PDF fit. Of course, the overall normalisation is forfeited but the sensitivity to the x dependence of the gluon is maintained in this approach. In Fig. 7 we show the NNPDF3.0 global gluon reweighted using the ratios of inclusive D-meson cross section data at √ s = 5, 7, 13 TeV and evolved down to the J/ψ scale µ 2 = 2.4 GeV 2 (the lower grey band).
As shown and explained in [5], the data favour a decreasing gluon at the lowest value of x which the D-meson data may probe. This is to be contrasted with the same analysis performed for NNPDF3.1 supplemented with the inclusive D-meson data but now together with small x resummation [9]. In this case, the reweighting favours a much higher gluon, as shown by the upper grey band in Fig. 7. It is known that including the BFKL (small x) resummation (without a k t < Q 0 subtraction) the low scale gluons extrapolated into the low x < 0.001 region are too large and grow too fast (see e.g. [31]). That is, as shown in Fig. 8, the cross section prediction using NNPDF3.1 together with the resummation strongly overshoots the exclusive J/ψ data while the prediction using NNPDF3.0 is too low.
The comparison of these two (based on NNPDF3.0 and on NNPDF3.1) bands, together with the inconsistencies of D-meson data mentioned in [10,11], demonstrates that the quality and accuracy of D-meson data are not sufficient to get an unambiguous result and to obtain accurate low x gluons.

Discussion
In this work, we too have performed a Bayesian reweighting of the NNPDF3.0 gluon but this time constrained by the exclusive J/ψ cross section. As discussed in [19] these data are in a position to be readily included in a collinear NLO global analysis due to alleviation of the large scale dependence through implementation of a Q 0 cut and resummation of a class of large  [5,9]. The latter includes low x resummation effects. The shaded blue band is the cross section prediction obtained based on our reweighting of the NNPDF3.0 NLO global gluon via the exclusive J/ψ data. The experimental data points are presented as in Fig. 4. logarithms. We have performed the reweighting using the J/ψ data in the region x < 0.01 for the NNPDF3.0 NLO set with N rep = 1000 replicas. Since the central NNPDF3.0 low x gluons are too large to describe the J/ψ data (see Fig. 4), the Shannon entropy (or effective number of contributing replicas), N eff ≈ 40 N rep . Therefore, the reweighting approach is not fully adequate. Still, the obtained gluons (hatched blue band in Fig. 7) are rather close to that obtained within the fit using ansatz (6). Since the NNPDF input distribution is mainly driven by other data at larger x ∼ 0.01 (where the effective value of λ is noticeably smaller), the reweighted NNPDF3.0 gluon has a slightly less steep growth at x < 0.001 in comparison with that coming from the power fit (6). Correspondingly, the J/ψ reweighted gluon density overshoots our (power fit) result at x = x 0 = 0.001 while undershooting it at the smallest x = 3 × 10 −6 . 3 On the other hand our J/ψ reweighting result demonstrates that the additional J/ψ data adds a lot of new information, which is to be expected as there were no data in the previous PDF analyses in this domain. The small value of the Shannon entropy means it would be desirable for the reweighting procedure to be backed up by a full new global fit. This quantifies the statements in [19] about the utility of the J/ψ data. The closeness of our reweighted gluon with the fitted gluons we have obtained provides further support for this claim. Considering all data points with W > 100 GeV, the effective χ 2 min /N dat ∼ 1.07 for the reweighted central cross section prediction.
Thus exploiting the J/ψ exclusive data we reach a much better accuracy. Now, down to x = 3 × 10 −6 , the low scale gluons (near the input Q 0 value) are known to better than 5-7% uncertainty.
An interesting observation is that in the low x < 0.001 region, the low scale fitted gluons start to grow (with 1/x) even faster (as xg(x) ∝ x −λ with λ 0.14) than the low scale global gluons do in the interval 0.001 < x < 0.01. We are able to fit a low x gluon power ansatz for the large range x < 0.001 with a single slope but find that we cannot extend this same description to 0.001 < x < 0.01. Attempting to do so results in a worsened fit and a much smaller λ. Indeed, this reflects the differing behaviour of the NLO global gluons in the intervals 0.001 < x < 0.01 and x < 0.001. The fact that the effective power λ increases with 1/x (within the 10 −2 − 10 −5 interval) is in contradiction with the assumption of saturation for which one would expect a decreasing λ → 0 as x → 0. The data with x < 0.01, therefore, cannot be described by a single power behaviour, indicative of non-trivial non-perturbative effects in the input proton wave function.
On the other hand note that the power λ 0.14 (that we obtained in the description of the J/ψ data with x < 0.001) is close to that predicted by the NLL BFKL re-summed with the optimal (BLM [32]) scale renormalization [33]. Moreover, contrary to the common expectation, even at x ∼ 10 −5 and µ 2 = 2.4 GeV 2 we see no hint for the beginning of parton density saturation.

Conclusion
High energy HERA and LHCb data on exclusive J/ψ production were described using a consistent collinear factorization approach at NLO. We fix the 'optimal' factorization scale µ F = M ψ /2, which allows for the resummation of the double-logarithmic (α s ln(1/x) ln µ F ) m corrections into the incoming PDF, and subtract the low k t < Q 0 contribution from the coefficient function to avoid double counting between the NLO coefficient function and the contribution hidden in the input PDF (or GPD) at Q = Q 0 . This provides good stability of the results with respect to variations of µ f . The generalized GPD distribution was related to the conventional (non-skewed) PDF via the Shuvaev transform. The renormalization scale is µ R = µ f .
With this, we find collinear NLO gluons at µ 2 = 2.4 GeV 2 which give an excellent description of all available accurate J/ψ data throughout the very low x interval, 3 × 10 −6 < x < 10 −3 , to about ± 5-7 % accuracy at the lowest x. The gluon PDF xg(x, µ 2 ) ∝ x −λ increases with 1/x with λ = 0.135 ± 0.006 without any hint in favour of parton density saturation.
A Bayesian reweighting approach leads to a similar behaviour of the small x gluon, emphasising the utility and constraining power of the exclusive J/ψ data. This work therefore clearly demonstrates the gains which will be achieved once these data are included in the global PDF fits.