IPPP/20/23

LTH 1237

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

C.A. Fletta, A.D. Martinb, M.G. Ryskinb,c and T. Teubnera

a Department of Mathematical Sciences, University of Liverpool, Liverpool, L69
3BX, U.K.

b Institute for Particle Physics Phenomenology, Durham University, Durham, DH1
3LE, U.K.

c Petersburg Nuclear Physics Institute, NRC Kurchatov Institute, Gatchina, St.
Petersburg, 188300, Russia

The low $x$ behaviour of the gluon density $xg(x,\mu^{2})$ at scale
$\mu^{2}=2.4$ GeV2 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.

##  1 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/\psi$ 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 ${\cal 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/\psi$
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/\psi$ data are
more consistent and have better accuracy than the inclusive $D$-meson cross
section.

Nevertheless, until now the $J/\psi$ data have not been used in global
analyses due to theoretical complications. First, the $J/\psi$ 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/\psi$ photoproduction may be written,
using collinear factorization, in the form [12]

$A=\frac{4\pi\sqrt{4\pi\alpha}e_{q}(\epsilon_{V}^{*}\cdot\epsilon_{\gamma})}{N_{c}}\left(\frac{\langle
O_{1}\rangle_{V}}{m_{c}^{3}}\right)^{1/2}\int_{-1}^{1}\frac{\mathrm{d}X}{X}\,\left[C_{g}\left(X,\xi\right)F_{g}(X,\xi)+C_{q}(X,\xi)F_{q}(X,\xi)\right],$
(1)

where we have suppressed the dependence on the renormalization and
factorization scales, $\mu_{R},\mu_{F}$, and on the invariant transferred
momentum squared, $t$. Here, the non-relativistic QCD (NRQCD) matrix element
$\langle O_{1}\rangle_{V}$ describes the formation of the $J/\psi$ meson with
$m_{c}$ the charm quark mass. The quantum numbers of the photon and the
$J/\psi$ meson select the charge conjugation even exchange in the $t$-channel,
so that only the charge conjugation even quark singlet and gluon GPDs, denoted
$F_{q}$ and $F_{g}$ respectively, contribute. The quark and gluon coefficient
functions $C_{q}$ and $C_{g}$ are known at NLO [12].

The kinematics of the process are displayed in Fig. 1. The partons carry
momentum fractions $(X+\xi)$ and $(X-\xi)$ of the plus-component of the mean
of the incoming and outgoing proton momenta, $P=(p+p^{\prime})/2$, where the
skewedness parameter $\xi$ is defined by

$\xi~{}=~{}\frac{p^{+}-p^{\prime+}}{p^{+}+p^{\prime+}}~{}=~{}\frac{M^{2}_{\psi}}{2W^{2}-M^{2}_{\psi}},$
(2)

with $M_{\psi}$ the mass of the $J/\psi$ and $W^{2}$ the photon-proton energy
squared for the $\gamma p\to J/\psi~{}p$ subprocess. It was shown that in the
relevant region, $\xi~{}<~{}10^{-3}$, 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 GPDs are
equal to the known Mellin moments of the non-skewed PDFs up to ${\cal O}(\xi)$
accuracy at NLO [15, 16]. The GPD grids are generated from PDFs supplied on
grids via the LHAPDF interface [17].

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
double log terms, [$\alpha_{s}\ln\mu^{2}_{F}\ln(1/x)]^{m}$, can be resummed in
the leading order term by choosing the optimum scale
$\mu_{F}=m_{c}=M_{\psi}/2$ for our process [18]. The NLO amplitude
$A(\mu_{f})$, with factorisation scale $\mu_{f}$, can be written schematically
in the form

$A(\mu_{f})~{}=~{}C^{\rm LO}\otimes{\rm GPD}(\mu_{F})~{}+~{}C^{\rm NLO}_{\rm
rem}(\mu_{F})\otimes{\rm GPD}(\mu_{f}).$ (3)

With the choice $\mu_{F}=M_{\psi}/2$, the remaining NLO coefficient function,
$C^{\rm NLO}_{\rm rem}(\mu_{F})$, does not contain terms enhanced by
ln$(1/x)\simeq\ \ln(1/\xi)$.

Figure 1:  LO (left panel) and NLO (right panel) contributions to $\gamma p\to
V+p$, where $V=J/\psi$ (or $\Upsilon$). Here the momentum
$P\equiv(p+p^{\prime})/2$, with $k$ the loop momentum. Note that the momentum
fractions of the left and right input partons are $x=X+\xi$ and
$x^{\prime}=X-\xi$ respectively; for the gluons coupled directly to the on-
shell heavy quark pair, we have $x^{\prime}\ll x$ and so $x\simeq 2\xi$.

The approach was described in more detail in [19] where it was shown that the
HERA data on diffractive $J/\psi$ photoproduction [20] with energies
corresponding to $x>10^{-3}$ are well described using the present global
gluons.111We 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, $c\bar{c}+g$, Fock
component of the $J/\psi$ 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/\psi\to l^{+}l^{-}$, and the charm quark mass is
chosen to be $m_{c}=M_{\psi}/2$, as is kept in the present paper. Note also
that the correction to the NRQCD matrix element changes the normalization of
the $J/\psi$ cross section but does not affect the $x$ (or $W$) behaviour of
the low-$x$ gluon. The fact that at $x\lower
3.01385pt\hbox{$\;\stackrel{{\scriptstyle\textstyle>}}{{\sim}}\;$}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. 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/\psi$ LHCb data [23]
(as well as HERA photoproduction data that lie in this region).

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/\psi$ data at NLO. In the
future this approach can be extended to NNLO.222This would require knowledge
of the 2-loop hard scattering coefficient function.

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/\psi$ 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/\psi$ data and compare
and contrast this with the gluon obtained from the above alternative
approaches. Our conclusions are briefly summarized in Section 6.

##  2 Ansatz for the low $x$ gluon

It was demonstrated in [19] that the diffractive $J/\psi$ cross section is
driven by the Generalised Parton Distributions, GPD$(X+\xi,X-\xi)$, of the
gluon with $X\simeq\xi$, see Fig. 1. That is, to describe the LHCb data, we
effectively need the gluon in the region of low $x\simeq X+\xi$ 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)\to 0$ as $x\to 1$. In order to satisfy
these requirements we choose the following ansatz for the conventional gluon
PDF,

$xg(x,\mu_{0}^{2})~{}=~{}C~{}xg^{\rm
global}(x,\mu_{0}^{2})~{}+~{}(1-C)~{}xg^{\rm new}(x,\mu_{0}^{2})$ (4) ${\rm
with}~{}~{}~{}~{}~{}C~{}=~{}\frac{x^{2}}{x^{2}+x^{2}_{0}}~{},$ (5)

and where $xg^{\text{global}}$ is the value of the gluon PDF obtained in a
global PDF analysis. The simplest low $x$ form for the gluon would be

$xg^{\rm new}(x,\mu_{0}^{2})~{}=~{}nN_{0}~{}(1-x)~{}x^{-\lambda},$ (6)

where the normalization factor $N_{0}$ is chosen so that for $n=1$ the gluon
PDF has the matching at $x=x_{0}$,

$x_{0}g^{\rm new}(x_{0},\mu_{0}^{2})~{}=~{}x_{0}g^{\rm
global}(x_{0},\mu_{0}^{2}).$ (7)

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\to 0$ as $x\to 1$.333Note that this factor was added to
satisfy the formal conditions for the validity of the Shuvaev transform.
Practically, the results do not depend on the behaviour of the gluon at
relatively large $x$. The corresponding effects are not visible in our Figs.
4, 5, 7. Due to the smooth form of $C$ in (5) the complete distribution (4)
does not violate analyticity even for $n\neq 1$.

Alternatively, in order to compare our present collinear determination of
$xg^{\rm new}$ with an earlier determination of the low $x$ gluon obtained in
the $k_{t}$ factorization approach [24], we also use an ansatz inspired by the
double logarithm approximation,

$xg^{\rm
new}(x,\mu_{0}^{2})~{}=~{}nN_{0}~{}(1-x)~{}x^{-a}\left(\frac{\mu_{0}^{2}}{q^{2}_{0}}\right)^{-0.2}~{}{\rm
exp}\left[\sqrt{16(N_{c}/\beta_{0})\ln(1/x)\ln G)}\right]$ (8) ${\rm
with}~{}~{}~{}~{}G~{}=~{}\frac{\ln(\mu_{0}^{2}/\Lambda^{2}_{\rm
QCD})}{\ln(q_{0}^{2}/\Lambda^{2}_{\rm QCD})}~{},$ (9)

where the parameter $a$ now plays the role of $\lambda$. Here, with three
light quarks $(N_{f}=3)$ and $N_{c}=3$ we have $\beta_{0}=9$. We take
$\Lambda_{\text{QCD}}=200\,\text{MeV}$ and $q_{0}^{2}=1\,\text{GeV}^{2}$, as
in [24], with $\mu_{0}^{2}=2.4\,\text{GeV}^{2}$ fixed. The exponent in (8)
resums, to all orders in $m$, the double logarithmic terms
$(\alpha_{s}\ln(1/x)\ln\mu^{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 GeV2. Therefore this parametrization can be used to describe
$\Upsilon$ photoproduction data as well.

##  3 Determination of the low $x$ gluon from $J/\psi$ data

Here, we show the results of our fits to $J/\psi$ 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/\psi$ 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$.

###  3.1 The exclusive $J/\psi$ data from LHCb

The LHCb experiment, by design, does not directly measure the cross section
for $J/\psi$ photoproduction but instead that for exclusive $pp\to p+J/\psi+p$
[23]. The experiment is unable to tag forward protons accompanying the
$J/\psi$ so instead only the rapidity of the $J/\psi$ is measured. Events are
selected by ensuring a large rapidity gap on both sides of the $J/\psi$
measurements, where the transverse momentum of the $J/\psi$ 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.

Figure 2: Two leading-order (LO) diagrams describing exclusive $J/\psi$
production at the LHC. The left diagram, the $W_{+}$ component, is the major
contribution to the $pp\to p+J/\psi+p$ cross section for a $J/\psi$ produced
at large rapidity $Y$. Thus such data allow a probe of very low $x$ values,
$x\sim M_{\psi}{\rm exp}(-Y)/\sqrt{s}\,$; recall that for two-gluon exchange
we have $x\gg x^{\prime}$. The $q_{T}$ of the photon is very small and so the
photon can be considered as a real on-mass-shell particle.

The ultraperipheral amplitude for a given $J/\psi$ rapidity is then generally
the sum of two photoproduction amplitudes with different $W^{2}$, depending on
which proton emitted the photon and which was the target, see Fig. 2. The
interference contribution is suppressed as the photon’s 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/\psi$
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 $\mathrm{d}\sigma(pp)/\mathrm{d}t$ cross section and low-
mass diffractive dissociation. The details of the procedure to extract
$\sigma(\gamma p\to J/\psi+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].

Figure 3:  LO and NLO cross section predictions obtained using the central
values of the existing global partons from [3]. Dashed (dotted) lines
correspond to the scale choices $\mu_{f}^{2}=\mu_{R}^{2}=m_{c}^{2}$
($\mu_{f}^{2}=\mu_{R}^{2}=2m_{c}^{2}$) with $\mu_{F}=Q_{0}=m_{c}$ fixed.
Figure 4: The description of the $J/\psi$ photoproduction HERA [20] and LHCb
[23] data based on using the central value of the global gluon PDF from the
three global parton analyses [1, 2, 3] for $x>0.001$. The red, blue and green
solid and dotted lines show the spread of descriptions based on the power fit,
using the $\pm\,1\sigma$ errors for the parameters. We also show by dashed
lines the cross section predictions obtained using the current central values
of the global gluons for all $x$.

###  3.2 Description of the $J/\psi$ data

The set-up of the cross section prediction follows [19]. Only the imaginary
part of the amplitude is computed using eq. (1). In this way we need only the
GPDs in the DGLAP region $|X|>\xi$. In the ERBL region, $|X|<\xi$, the
imaginary part of the coefficient function is zero. The real part is then
restored via a dispersion relation, which in the high energy limit (for the
even signature amplitude) can be written in the simplified form [26]

$\frac{{\rm Re}A}{{\rm Im}A}~{}~{}=~{}~{}{\rm
tan}\left(\frac{\pi}{2}~{}\frac{\partial(\ln{\rm Im}A/W^{2})}{\partial(\ln
W^{2})}\right).$ (10)

Next, we use NRQCD to describe the formation of the $J/\psi$ wave function. We
project the open heavy $c\bar{c}$ quark pair onto the colour singlet
configuration with the corresponding transition matrix element $\langle
O_{1}\rangle_{V}$, which is fixed by the experimentally measured leptonic
decay width of the $J/\psi$. The exclusive final state requires a colourless
high energy scattering (modelled by the two-gluon exchange) and does not allow
for an octet contribution, as this would populate the rapidity gap and destroy
the exclusivity of the final state.

Note that actually we calculate the value of Im$A$ at $t=0$ and then restore
the total $\gamma p\to J/\psi+p$ cross section assuming an exponential $t$
behaviour with a slope

$B=4.9+4\alpha^{\prime}_{P}\ln(W/W_{0})~{}~{}\mbox{GeV}^{-2}$

with $W_{0}=90$ GeV and $\alpha^{\prime}_{P}=0.06$ GeV-2. This parametrisation
grows more slowly with $W$ than the formula used by H1 [20], but is still
compatible with the HERA data. We have chosen the slope parameter
$\alpha^{\prime}_{P}$ to be compatible with Model 4 of [27] which fits a wider
variety of 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/\psi$ 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 ($\lambda$
and $n$, as defined in eq. (6)) fit of all the $\sigma(\gamma p\to J/\psi+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 and Fig. 4.

| $\lambda$ | $n$ | $\chi_{\text{min}}^{2}$ | $\chi_{\text{min}}^{2}/\text{d.o.f}$
---|---|---|---|---
NNPDF3.0 | 0.136 | 0.966 | 44.51 | 1.04
MMHT14 | 0.136 | 1.082 | 47.00 | 1.09
CT14 | 0.132 | 0.946 | 48.25 | 1.12
Table 1: The values of $\lambda$ and $n$ obtained from fits to the $J/\psi$
data using three sets of global partons. The respective values of the total
$\chi_{\text{min}}^{2}$ (and $\chi_{\text{min}}^{2}/\text{d.o.f}$) for 45 data
points are also shown.

The respective values of the $\chi_{\text{min}}^{2}$ 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 [28]. We use all HERA data points [20] with $W>100$ GeV and all
LHCb [23] data points.

For the ZEUS 2002 and 2004 data sets [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 $\sim 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 $\sim 4\%$.

The description of the exclusive $J/\psi$ cross section is shown in Fig. 4,
while the gluons extracted from the $J/\psi$ data at $\mu^{2}=2.4$ GeV2 and
$x<0.001$ are shown in Fig. 5. The error bands are obtained by sampling over
the two parameters within their individual 1$\sigma$ 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. As is seen from Fig. 5, the resulting 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 ($\chi_{\text{min}}^{2}$) of the
description. The typical errors are $\pm 2.5$% for the normalization and $\pm
4$% for $\lambda$. We see from Fig. 4 that the simple two-parameter form of
the gluon density provides an excellent description of the $J/\psi$ 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\sim 3\times 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.

Figure 4 also shows the cross section predictions obtained using the central
values of the gluon from the global parton sets extrapolated into the low $x$
region. Clearly here the global analyses have no predictive power and in each
case they have huge uncertainty bands (shown in Fig. 5 for NNPDF3.0 only)
which cover the (unfitted) $J/\psi$ data. The value of including the $J/\psi$
data is apparent.

Figure 5: The cross-hatched region shows the range of behaviour of the low $x$
NLO gluon determined by fitting to exclusive $J/\psi$ data using ansatz (6)
with $xg^{\rm global}$ taken from NNPDF3.0 [1], MMHT14 [2] or CT14 [3] parton
sets.  Figure 6:  (a) The global gluon PDF, $xg(x,\mu^{2})$, at the matching
point $x=0.001$ and $\mu^{2}=2.4$ GeV2, (b) the global gluon PDF,
$xg(x,\mu^{2})$, at the matching point $x=0.001$ and $\mu^{2}=2.4$ GeV2 after
fitting to HERA+LHCb exclusive $J/\psi$ data. Note that the errors shown on
the right hand side are those obtained by propagating the $1\sigma$
experimental data errors to our result, but do not account for theoretical
uncertainties.

In the left hand side of Fig. 6 we compare the uncertainties of the gluon
densities given at $x=0.001$ and $\mu^{2}=2.4$ GeV2 by the global analyses,
while in the right hand side we show the values that are obtained after
fitting the $J/\psi$ data. The $J/\psi$ data are seen to greatly improve the
knowledge of the gluon in the low $x$ interval $3\times 10^{-6}<x<10^{-3}$. In
particular, we find at $x_{0}=0.001$ that

$x_{0}g(x_{0},\mu^{2}=2.4\,\text{GeV}^{2})=2.28\pm 0.06.$ (11)

###  3.3 The alternative double-log parametrization

While the simple two parameter ansatz in (6) leads to a very good description
of the $J/\psi$ 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\pm 0.006,\,\,\,\,\,n=0.979\pm
0.025,\,\,\,\,\,\chi^{2}_{\text{min}}/\text{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},\mu^{2}=2.4\,\text{GeV}^{2})=2.31$ in agreement with (11).

Note that the double-log parametrization gives a result close to that obtained
in the $k_{t}$-factorization 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\sim-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$.

###  3.4 Is there evidence of saturation from exclusive $J/\psi$ data?

High energy exclusive $J/\psi$ production was recently described in [29] 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/\psi$
at small $x$”. We have to emphasize that actually the authors of [29] refer to
absorptive corrections rather than saturation. Indeed, saturation means that
the gluon density tends to a constant value, $xg(x,\mu^{2})\to const$ as $x\to
0$ and at a fixed scale $\mu$ [30]. That is, the power $\lambda$ in (6)
behaves as $\lambda\to 0$. A first hint of saturation would be to observe that
the power $\lambda$ (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 [29] is that the LO BFKL intercept,
$\alpha_{\text{BFKL}}=1+\omega_{0}=1+\lambda$ is too large to describe the
high energy $J/\psi$ data and that absorptive corrections (which are included
into the non-linear BK [31] equation) are needed to tame the growth of the
gluon density (6), that is to decrease the value of $\lambda$.

It is well known that the LO BFKL intercept is too large. It becomes smaller
in the next-to-leading (NLL) approximation. Indeed, it is seen from [29] (the
short dashed green curve of their Fig. 1) that the HSS gluons [32], based on
the NLO BFKL linear equation, are in agreement with the exclusive $J/\psi$
data.

Therefore the growth of the gluon density with a smaller but non-zero
$\lambda$ is not evidence for ‘saturation’. At the moment no hint of
saturation is observed in exclusive $J/\psi$ data at the scale $\mu^{2}=2.4$
GeV2 and $x$ down to 10-5.

###  3.5 Note on higher-twist contributions

Recall that absorptive corrections, which provide the saturation, are
described by higher-twist operators. Formally, within the collinear
factorization approach, we do not know the value of these higher-twist terms.
They have their own evolution and input conditions/functions that must be
fitted from experiment. In other words, only experiment can give us the values
of the higher-twist operator contributions. Nevertheless, let us estimate the
possible role of the higher-twist absorptive effects in the $J/\psi$
photoproduction amplitude.

The relative size of the contribution of the next twist absorptive correction
(in our $\mu^{2}$ region of interest) is driven by the parameter (see
[30])444In our approach everything below $Q_{0}$ (i.e. at scales $\mu^{2}<2.4$
GeV2) is considered as a phenomenological input distribution which is formed
mainly by non-perturbative interactions inside the proton. We never go below
$Q_{0}$; we subtract all the contributions with $k_{t}<Q_{0}$. One therefore
cannot use our higher-twist estimate (of perturbative origin) at lower scales.

$c=\alpha_{s}\frac{xg(x)}{R^{2}\mu_{0}^{2}},$ (12)

where $R$ can be as large as the proton radius $(R\sim 0.84~{}\text{fm})$. If
we consider the value of $R$ as the ‘hot spot’ radius555It may be assumed that
the low-$x$ partons group together in so-called ‘hot-spots’, with a radius
smaller than that of the proton., then we have to take a smaller gluon
density, $xg$, corresponding to only one hot spot. With $\alpha_{s}$ = 1/3 and
$\mu_{0}=M_{\psi}/2$ we obtain $c=0.008\,xg\sim 0.04$ for our gluon density
$xg\leq 5$.666A relatively large value of $xg=5$ includes/accounts for the
power growth of gluon densities at low $x$. However, actually this result is
overestimated. Indeed, the cross section of an additional high energy (gluon)
interaction is proportional to the $c$-quark separation $\langle
r^{2}\rangle.$ This means that we have to replace in (12) the factor
$1/(R^{2}\mu_{0}^{2})$ by the ratio $\langle r^{2}\rangle/R^{2}$. At the
beginning of the photoproduction process, the photon produces a point-like
$c\bar{c}$ pair. The lifetime of this pair is about
$2E_{\gamma}/M_{\psi}^{2}$, where $E_{\gamma}$ is the energy of the photon.
Accounting in addition for the Lorentz factor of the $J/\psi$, the quarks have
their ‘own’ time $\tau\sim 2/M_{\psi}=1/\mu_{0}$ to separate from each other.
However, the $J/\psi$ meson is a non-relativistic system and the heavy quark
velocity $\langle v^{2}\rangle\propto\alpha_{s}$ is small. That is we expect
the higher-twist contribution to be suppressed by an additional power of
$\alpha_{s}$ and, correspondingly, actually $c<0.015.$ Accounting for the
velocity $\langle v^{2}\rangle\propto\alpha_{s}$ can be considered as a NNLO
contribution.

##  4 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/\psi$ production.

Inclusive $D$-meson production data via $pp$ collisions at the LHC are
available at centre of mass energies $5,7$ and $13\;\text{TeV}$ [4]. The
kinematics of the different modes of production of the $D$-mesons allow for a
coverage down to $x\sim{\rm few}\times 10^{-6}$. In [5] the authors studied
the impact these data for $\left\\{D^{0},D^{+},D^{+}_{s}\right\\}$ 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
$\sqrt{s}=5,7,13~{}\text{TeV}$ and evolved down to the $J/\psi$ scale
$\mu^{2}=2.4\,\text{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.

Figure 7: Comparison of the low $x$ behaviour of the NLO gluon density
$xg(x,\mu^{2})$ at $\mu^{2}=2.4$ GeV2 obtained from exclusive $J/\psi$ data
and from inclusive $D$-meson data, see text for details.

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. [33]). That is, as shown in Fig. 8, the cross section
prediction using NNPDF3.1 together with the resummation strongly overshoots
the exclusive $J/\psi$ data while the prediction using NNPDF3.0 is too low.

Figure 8: The lower and upper bands are, respectively, the cross section
predictions obtained using NNPDF3.0 and NNPDF3.1 global partons constrained by
the $D$-meson LHCb data [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/\psi$
data. The experimental data points are presented as in Fig. 4.

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.

##  5 Discussion

In this work, we too have performed a Bayesian reweighting of the NNPDF3.0
gluon but this time constrained by the exclusive $J/\psi$ 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
logarithms. We have performed the reweighting using the $J/\psi$ data in the
region $x<0.01$ for the NNPDF3.0 NLO set with $N_{\text{rep}}=1000$ replicas.
Since the central NNPDF3.0 low $x$ gluons are too large to describe the
$J/\psi$ data (see Fig. 4), the Shannon entropy (or effective number of
contributing replicas), $N_{\text{eff}}\approx 40\ll N_{\text{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\sim 0.01$ (where the effective value of $\lambda$
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/\psi$ reweighted gluon density overshoots our
(power fit) result at $x=x_{0}=0.001$ while undershooting it at the smallest
$x=3\times 10^{-6}$.777The slightly larger normalization, at $x=10^{-3}$, of
the prediction based on the reweighting procedure is due to the greater number
of data points that are fitted in this region in the global analysis. For
smaller $x$, where the only constraining power comes from the exclusive
$J/\psi$ data in both the reweighting and power fit approaches, the
predictions are in better agreement. On the other hand our $J/\psi$
reweighting result demonstrates that the additional $J/\psi$ 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/\psi$ 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\,\text{GeV},$ the effective
$\chi^{2}_{\text{min}}/N_{\text{dat}}\sim 1.07$ for the reweighted central
cross section prediction.

Thus exploiting the $J/\psi$ exclusive data we reach a much better accuracy.
Now, down to $x=3\times 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)\propto
x^{-\lambda}$ with $\lambda\simeq 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 $\lambda$. 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 $\lambda$
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 $\lambda\to 0$ as $x\to 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 $\lambda\simeq 0.14$ (that we obtained
in the description of the $J/\psi$ data with $x<0.001$) is close to that
predicted by the NLL BFKL re-summed with the optimal (BLM [34]) scale
renormalization [35]. Moreover, contrary to the common expectation, even at
$x\sim 10^{-5}$ and $\mu^{2}=2.4$ GeV2, in our approach we see no hint in the
exclusive $J/\psi$ data for the onset of parton density saturation.

##  6 Conclusion

High energy HERA and LHCb data on exclusive $J/\psi$ production were described
using a consistent collinear factorization approach at NLO. We fix the
‘optimal’ factorization scale $\mu_{F}=M_{\psi}/2$, which allows for the
resummation of the double-logarithmic $(\alpha_{s}\ln(1/x)\ln\mu_{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 $\mu_{f}$. The generalized GPD distribution was related to the
conventional (non-skewed) PDF via the Shuvaev transform. The renormalization
scale is $\mu_{R}=\mu_{f}$.

With this, we find collinear NLO gluons at $\mu^{2}=2.4$ GeV2 which give an
excellent description of all available accurate $J/\psi$ data throughout the
very low $x$ interval, $3\times 10^{-6}<x<10^{-3}$, to about $\pm$ 5-7 %
accuracy at the lowest $x$. The gluon PDF $xg(x,\mu^{2})\propto x^{-\lambda}$
increases with $1/x$ with $\lambda=0.135\pm 0.006$ without any hint in favour
of parton density saturation at $\mu^{2}=2.4\,\text{GeV}^{2}$ and $x$ down to
$10^{-5}$. We emphasize this does not mean that the data cannot be described
by a more complicated expression which ultimately (at very small $x$) will
provide 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/\psi$ data. This work therefore clearly demonstrates the gains which will
be achieved once these data are included in the global PDF fits.

## Acknowledgements

We thank Stephen P. Jones for his major contributions (both theoretical and
phenomenological) to our previous analyses of exclusive $J/\psi$ production
which form much of the basis for the present study, and for his careful
reading of our manuscript. The authors would also like to thank Valerio
Bertone for a useful discussion and for providing the $D-$meson constrained
NNPDF parton sets in LHAPDF6 format. C.A.F and M.G.R thank the IPPP at Durham
University for hospitality. The work of C.A.F is supported by an STFC award
grant ST/N504130/1 and that of T.T is supported by STFC under the consolidated
grants ST/P000290/1 and ST/S000879/1.

## References

  * [1] R.D. Ball et al. [NNPDF Collaboration], JHEP 1504 (2015) 040 [arXiv:1410.8849].
  * [2] L.A. Harland-Lang, A.D. Martin, P. Motylinski, R.S. Thorne, Eur. Phys. J. C75 (2015) 204 [arXiv:1412.3989].
  * [3] S. Dulat et al., Phys. Rev. D93 (2016) 033006 [arXiv:1506.07443].
  * [4] LHCb Collaboration: R. Aaij et al., Nucl. Phys. B871 (2013) 1; JHEP 1603 (2016) 159, erratum: JHEP 1609 (2016) 013; JHEP 1705 (2017) 074; JHEP 1706 (2017) 147.
  * [5] R. Gauld, J. Rojo, Phys. Rev. Lett. 118 (2017) 072001 [arXiv:1610.09373].
  * [6] O. Zenaiev et al. [PROSA Collaboration], Eur. Phys. J. C75 (2015) 396 [arXiv:1503.04581].
  * [7] R. Gauld, J. Rojo, L. Rottoli, J. Talbert, JHEP 1511 (2015) 009 [arXiv:1506.0802].
  * [8] M. Cacciari, M. L. Mangano, P. Nason, Eur. Phys. J. C75 (2015) 610 [arXiv:1507.06197].
  * [9] V. Bertone, R. Gauld, J. Rojo, JHEP 01 (2019) 217 [arXiv:1808.02034].
  * [10] R. Gauld, JHEP 05 (2017) 084 [arXiv:1703.03636].
  * [11] E.G. de Oliveira, A.D. Martin, M.G. Ryskin, Phys. Rev. D97 (2018) 074021 [arXiv:1712.06834].
  * [12] D.Yu. Ivanov, A. Schafer, L. Szymanowski, G. Krasnikov, Eur. Phys. J. C34 (2004) 297, Erratum: Eur. Phys. J. C75 (2015) 75 [hep-ph/0401131].
  * [13] M. Diehl, W. Kugler, Eur. Phys. J. C52 (2007) 933 [arXiv:0708.1121].
  * [14] S.P. Jones, A.D. Martin, M.G. Ryskin, T. Teubner, Eur. Phys. J. C76 (2016) 633 [arXiv:1610.02272].
  * [15] A.G. Shuvaev, K.J. Golec-Biernat, A.D. Martin, M.G. Ryskin, Phys. Rev. D60 (1999) 014015 [hep-ph/9902410];
A.G. Shuvaev, Phys. Rev. D60 (1999) 116005 [hep-ph/9902318].

  * [16] A.D. Martin, C. Nockles, M.G. Ryskin, A.G. Shuvaev, T. Teubner, Eur. Phys. J. C63 (2009) 57 [arXiv:0812.3558].
  * [17] A. Buckley, J. Ferrando, S. Lloyd, K. Nordstr$\ddot{\text{o}}$m, B. Page, M. R$\ddot{\text{u}}$fenacht, M. Sch$\ddot{\text{o}}$nherr, G. Watt, Eur. Phys. J. C75 (2015) 132 [arXiv:1412.7420].
  * [18] S.P. Jones, A.D. Martin, M.G. Ryskin, T. Teubner, J. Phys. G43 (2016) 035002 [arXiv:1507.06942].
  * [19] C.A. Flett, S.P. Jones, A.D. Martin, M.G. Ryskin, T. Teubner, Phys. Rev. D101 (2020) 094011 [arXiv:1908.08398].
  * [20] ZEUS Collaboration (S. Chekanov et al.) Eur. Phys. J. C24 (2002) 345 [hep-ex/0201043]; Nucl. Phys. B695 (2004) 3 [hep-ex/0404008];
H1 Collaboration (A. Aktas et al.), Eur. Phys. J. C46 (2006) 585 [hep-
ex/0510016]; (C. Alexa et al.) Eur. Phys. J. C73 (2013) 2466
[arXiv:1304.5162].

  * [21] T. Lappi, H. M$\ddot{\text{a}}$ntysaari, J. Penttala, arXiv:2006.02830.
  * [22] P. Hoodbhoy, Phys. Rev. D56 (1997) 388 [hep-ph/9611207].
  * [23] LHCb Collaboration: R. Aaij et al., J. Phys. G41 (2014) 055002 [arXiv:1401.3288]; JHEP 1810 (2018) 167 [arXiv:1806.04079].
  * [24] S.P. Jones, A.D. Martin, M.G. Ryskin, T. Teubner, J. Phys. G44 (2017) 03TL01 [arXiv:1611.03711].
  * [25] V.A. Khoze, A.D. Martin, M.G. Ryskin, Eur. Phys. J. C74 (2014) 2756 [arXiv:1312.3851].
  * [26] M. G. Ryskin, R. G. Roberts, A. D. Martin and E. M. Levin, Z. Phys. C 76 (1997) 231 [hep-ph/9511228].
  * [27] V.A. Khoze, A.D. Martin, M.G. Ryskin, Eur. Phys. J. C73 (2013) 2503 [arXiv:1306.2149].
  * [28] R. D. Ball et al. [NNPDF], JHEP 05 (2010) 075 [arXiv:0912.2276].
  * [29] A. Arroyo Garcia, M. Hentschinski, K. Kutak, arXiv:1904.04394.
  * [30] L.V. Gribov, E.M. Levin, M.G. Ryskin, Phys. Rept. 100 (1983) 1.
  * [31] I. Balitsky, Nucl. Phys. B463 (1996) 99 [hep-ph/9509348];
Y. V. Kovchegov, Phys. Rev. D60 (1999) 034008 [hep-ph/9901281].

  * [32] M. Hentschinski, A. Sabio Vera, C. Salas, Phys. Rev. Lett. 110 (2013) 041601 [arXiv:1209.1353]; Phys. Rev. D87 (2013) 076005 [arXiv:1301.5283].
  * [33] H. Abdolmaleki et al. [xFitter Developers’ Team], Eur. Phys. J. C78 (2018) 621 [arXiv:1802.00064].
  * [34] S.J. Brodsky, G.P. Lepage, P.B. Mackenzie, Phys. Rev. D28 (1983) 228.
  * [35] S. J. Brodsky, V. S. Fadin, V. T. Kim, L. N. Lipatov and G. B. Pivovarov, JETP Lett. 70 (1999) 155-160 [hep-ph/9901229].