<figure><img src="image_1.png" /></figure>

# Including exclusive heavy vector-meson production data in a global parton analysis

A.D Martin, M.G Ryskin, T Teubner, Christopher Flett

► To cite this version:

A.D Martin, M.G Ryskin, T Teubner, Christopher Flett. Including exclusive heavy vector-meson production data in a global parton analysis. European Physical Journal C: Particles and Fields, 2025, 85 (4), pp.434. 10.1140/epjc/s10052-025-14142-9. hal-04667847

HAL Id: hal-04667847

https://hal.science/hal-04667847v1

Submitted on 16 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A method to include exclusive heavy vector-meson production data at small $x$ in global parton analyses

C. A. Flett$^{1,a}$, A. D. Martin$^{2}$, M. G. Ryskin$^{3}$, T. Teubner$^{4}$

Received: 29 November 2024 / Accepted: 29 March 2025

© The Author(s) 2025

**Abstract** We propose a method which allows the inclusion of exclusive heavy vector-meson production data at low $x$ in future global parton analyses. As an example we perform a study withinxFitter to determine the gluon parton distribution function (PDF) at next-to-leading order (NLO) at moderate-to-low $x$ using the measurements of exclusive $J/\psi$ production in ep and pp collisions from HERA and LHC. We further study the constraints from the corresponding $\Upsilon$ production process. We finish by discussing the possible effects at next-to-next-to-leading order (NNLO) through incorporation of a $K$ factor for the exclusive heavy vector-meson coefficient function at NLO.

## 1 Introduction

The precision of the parton distributions of the proton is well established through global analyses [1–3], provided the resolution scale $Q^2$ is not too low and the momentum fraction $x$ remains within a moderate range, neither too small nor too large. As $x$ decreases, particularly at lower scales, the uncertainty in these distributions increases significantly. This increase in uncertainty is due to the lack of experimental data directly probing this region.

The experiments at the Large Hadron Collider (LHC) are capable of particle detection and reconstruction over a wide rapidity range. Notably, various measurements of the differential cross sections for exclusive heavy vector mesons such as $J/\psi$ and $\Upsilon$ [4–6] have allowed for the determination of the gluon Parton Distribution Function (PDF) down to $x \sim 3 \times 10^{-6}$ at factorisation scales $\mu_F = m_q$, where $q = c, b$.

Unfortunately it is not easy to include these data in global parton analyses. The cross section for exclusive meson production is driven not by the usual (diagonal) PDFs but by the more complicated (skewed) generalised parton distributions (GPDs), and is proportional to the skewed gluon-density squared, as indicated in Fig. 1; see [7] for a review. Note from the caption of the figure that $x \approx 2\xi$ at very small $x$, where $2\xi$ is the proton momentum fraction transferred through the GPD to the vector meson.

However, in the small $x$ region of our interest, the value of the GPD can be calculated from the conventional PDF with good ($\sim O(x)$) accuracy using the Shuvaev transform [8,9]. The Shuvaev transform makes use of the fact that as $\xi \to 0$ (and transverse-momentum transfer $p_T = 0$) the Gegenbauer moments$^1$ of the GPD become equal to the known Mellin moments of the PDF. Due to the polynomial condition (see e.g. [11]) even for $\xi \neq 0$ the Gegenbauer moments can be obtained from the Mellin moments to $O(\xi)$ accuracy. Thus it is possible to obtain the full GPD function at small $\xi$ from its known moments.

In principle, this allows us to include the low-$x$ exclusive $J/\psi$ data in the global PDF analysis at NLO. In practice the problem is that already at NLO the full Shuvaev transform amounts to a slowly convergent double integral, and the corresponding computation is time consuming as the integrand depends on the derivative of a PDF obtained through interpolation on a grid. This is troublesome for the global fit as, after every iteration, the derivative of the PDF varies and so a new grid defined over the variables $X, \xi$ and

<footnote>$^{1}$ Université Paris-Saclay, CNRS, IJCLab, 91405 Orsay, France</footnote>

<footnote>$^{2}$ Institute for Particle Physics Phenomenology, Durham University, Durham DH1 3LE, UK</footnote>

<footnote>$^{3}$ Petersburg Nuclear Physics Institute, NRC Kurchatov Institute, Gatchina, St. Petersburg 188300, Russia</footnote>

<footnote>$^{4}$ Department of Mathematical Sciences, University of Liverpool, Liverpool L69 3BX, UK</footnote>

<footnote>$^a$e-mail: christopher.flett@ijclab.in2p3.fr (corresponding author)</footnote>

<footnote>$^1$Gegenbauer moments are the analogue of Mellin moments which diagonalize the $Q^2$ evolution of PDFs. The corresponding operator diagonalizes the $Q^2$ evolution of the GPDs [10]. As $\xi \to 0$ the Gegenbauer moments become equal to the Mellin moments.</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1** (Left) LO contribution to $\gamma p \to V + p$, where the vector meson $V = J/\psi$, $\Upsilon$. (Right) NLO quark contribution. For these graphs all permutations of the parton lines and couplings of the gluon lines to the heavy-quark pair are to be understood. The NLO gluon contribution, with coefficient function $C_g^{\text{NLO}}$ and GPD $F_g$, dominates the NLO quark</figcaption></figure>

<figure><img src="image_4.png" /></figure>

contribution. In these diagrams, the momentum $P \equiv (p + p')/2$ and $l$ is the loop momentum. Note that the momentum fractions of the left and right partons are $x = X + \xi$ and $x' = X - \xi$ respectively; for the gluons connected to the heavy quark-antiquark pair, we have $x' \ll x$ and so $x \simeq 2\xi$.

$Q^2$ has to be computed to obtain the updated theory prediction. That is, after every iteration, the functional form of the PDF evolves, necessitating the computation of a new GPD grid using the Shuvaev transform. Since the functional forms of PDFs depend on various parameters that can all vary, it is impractical to precompute a grid exhaustively sampling a vast parameter space for every possible form. Therefore, an iterative computation is essential. It may be possible to employ a fast-gridding approach similar to that used for computing inclusive cross sections, e.g. as incorporated in FastNLO [12], APPLgrid [13] or PineAPPL [14], but here we recall the dependence on the PDF for exclusive observables is at the level of the amplitude rather than the total cross section, and so such an approach is less tractable.

In the present paper we describe a method which can overcome this problem by replacing the exclusive vector-meson data at small $x$ by 'effective gluon points'. To demonstrate the efficiency of the proposed method and to show how the inclusion of the LHCb $J/\psi$ and $\Upsilon$ exclusive cross sections affect the result we compare the standard xFitter [15] NLO parton analysis of DIS data with that supplemented by the vector-meson data. Note that our goal is not to present a new set of parton distributions but to emphasise the possible role of the vector-meson data in global parton analyses. The new 'effective gluon points' presented in Tables 1 and 2 can be used in future global analyses.

## 2 Effective low-x gluon PDF data

To overcome this difficulty, we propose the following procedure. We translate the experimental $J/\psi$ cross section data into a set of “effective” values for the gluon PDF.

As seen from Fig. 2, the cross section for the exclusive process

$$p + p \to p + J/\psi + p \qquad (1)$$

contains two components,

$$\frac{d\sigma^{(pp \to p+J/\psi+p)}}{dY} = S^2(W_+) \left(k_+ \frac{dn}{dk_+}\right) \sigma_+(\gamma p \to J/\psi p) + W_- \text{ term,} \qquad (2)$$

corresponding to energies $W_\pm$, where $W_+$ denotes the energy when the vector meson travels in the same direction as the photon, and $W_-$ when it travels in the opposite direction. We see that the values of the subprocess cross sections $\sigma(\gamma p \to J/\psi p)$, for each rapidity $Y$ of the $J/\psi$, are weighted by the corresponding gap survival factors $S^2(W_\pm)$ [16] and photon fluxes $dn/dk_\pm$ [17]. The cross section with the lower $\gamma p$-energy ($W_-$) was measured at HERA [18–21] and can be calculated with sufficient accuracy since it corresponds to relatively large $x$ where the uncertainties in the parton distributions are small. Hence the exclusive $J/\psi$ data provide reliable values for $\sigma_+(\gamma + p \to J/\psi + p)$, corresponding to the component with the larger $\gamma p$-energy, $W_+$, see Fig. 3.<fnref n="2" />

These $\sigma_+$ data were well described in [22] by a power-like low-$x$ gluon

$$xg(x) = A \cdot x^{-\lambda} \quad \text{with} \quad \lambda = 0.135 \pm 0.006, \qquad (3)$$

<footnote n="2">² Actually, the subtraction of the $W_-$ contribution was done already in the LHCb analyses where the results are presented as the $\sigma_+(\gamma + p \to V + p)$ cross sections.</footnote>
<figure><img src="image_5.png" /><figcaption>Fig. 2 The two diagrams describing exclusive J/ψ production at the LHC. The vertical lines represent two-gluon exchange. 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, x ~ M$_{J/ψ}$ exp(-Y)/√s, where √s is the centre-of-mass energy of the pp system; recall that for two-gluon exchange we have x ≫ x'. The q$_{T}$ of the photon is very small and so the photon can be considered as a real on-mass-shell particle</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 3 The description of the J/ψ photoproduction HERA [18–21] and LHCb [4,5] data based on using the central value of the global gluon PDF from the MMHT14 global analysis for x > 10⁻<fnref n="3" /> and a fitted power-like gluon PDF, see Eq. (3), for x < 10⁻<fnref n="3" />. The blue solid line and shaded region show, respectively, the central prediction and propagation of the ±1σ errors on the fitted parameters to the cross-section level</figcaption></figure>

also shown in Fig. 3. This provides the initialisation values for xg$_{fit}$(x), described later.

In collinear factorisation at NLO accuracy, the amplitude for the exclusive photoproduction of heavy quarkonium off unpolarised protons, γp → Vp, has the form

$$
\begin{equation}
\begin{split}
& \mathcal{M}(\gamma + p \rightarrow V + p) \\
& = \sum_{i=g,q} N_i \int_{-1}^{1} \frac{dX}{X} F_i(X, \xi) C_i^{\text{NLO}}(\xi/X),
\end{split}
\tag{4}
\end{equation}
$$

where the constants N_i provide the correct normalisation, see e.g. [22], and F_i and C_i denote quark and gluon GPDs and coefficient functions respectively. Additional dependences on the renormalisation and factorisation scales, and the four-momentum transfer t = t_min = (2ξm_p)²/(1 - ξ²) ~ 0, are not shown explicitly in Eq. (4).

We note here the sensitivity to the PDF or GPD at the amplitude level. In our approach, the small ξ GPD is obtained from the PDF using the Shuvaev transform [8,9]

$$
\begin{align}
F_q(X, \xi) &= \sqrt{1 - \xi^2} \int_{-1}^{1} dx' \left[ \frac{2}{\pi} \operatorname{Im} \right. \nonumber \\
&\qquad \left. \times \int_{0}^{1} \frac{ds}{y(s)\sqrt{1 - y(s)x'}} \right] \frac{d}{dx'} \left( \frac{q(x')}{|x'|} \right), \nonumber \\
F_g(X, \xi) &= \sqrt{1 - \xi^2} \int_{-1}^{1} dx' \nonumber \\
&\qquad \times \left[ \frac{2}{\pi} \operatorname{Im} \int_{0}^{1} \frac{ds (X + \xi(1-2s))}{y(s)\sqrt{1 - y(s)x'}} \right] \nonumber \\
&\qquad \frac{d}{dx'} \left( \frac{g(x')}{|x'|} \right), \tag{5}
\end{align}
$$

where the kernel of the transform,

$$
y(s) = \frac{4s(1-s)}{X + \xi(1 - 2s)}. \qquad (6)
$$

The subprocess cross sections

$$
\begin{align}
\sigma(\gamma + p \rightarrow V + p) &= \frac{1}{B(W)} \left. \frac{d\sigma}{dt} \right|_{t=0}, \quad \text{where} \nonumber \\
\left. \frac{d\sigma}{dt} \right|_{t=0} &= \frac{|M|^2}{16\pi W^4}, \tag{7}
\end{align}
$$

and B-slope B(W) GeV² = B₀ + 4α'ₚ ln(W/90). The inter-
cept B₀ = 4.9 (4.63) for J/ψ (Υ) production and α'ₚ =
0.06 [23].

Moreover it was shown that after the optimal factorisa-
tion scale,<fnref n="3" /> µ<sub>F</sub> = M<sub>J/ψ</sub>/2, was chosen, where M<sub>J/ψ</sub> is the
mass of the J/ψ, and the so-called Q<sub>0</sub> subtraction [26]<sup>4</sup> was
performed, the quark contribution to σ<sub>+</sub> becomes negligi-
ble [28]. That is the value of

$$
\sigma_+ \propto (2\xi g(2\xi))^2 \quad (8)
$$

is just proportional to the gluon density squared.<fnref n="5" />

The imaginary part of the LO coefficient function obtained
through the on-shell Cutkosky-cut of the LO Feynman dia-
grams (an example of which is shown in Fig. 1), enforces
X = ξ. At NLO, this is no longer true and there is a smearing

<footnote n="3">³ At this optimal scale the large double-logarithmic contributions of the form $c_n(\alpha_s \ln(1/\xi) \ln \mu_F^2)^n$ are resummed into the incoming PDF/GPD [24,25].</footnote>

<footnote n="4">⁴ The $Q_0$ subtraction is necessary to exclude the double counting between the low $l_T$ part of the coefficient function and the PDF-input. This, together with the resummation of the aforementioned double-logarithmic terms [24,25], helps mitigate the problem of the strong factorisation scale dependence [27] of the NLO amplitude.</footnote>

<footnote n="5">⁵ Strictly speaking, the value of σ⁺ is proportional to the square of the corresponding GPD.</footnote>
effect (i.e. a tail around $X = \xi$), but the integrand remains highly-peaked around this point, particularly for the gluon coefficient function [28]. In principle, while each value of the skewness parameter $\xi$ entails an integration over $X$, the dominant contribution comes from $X \approx \xi$ and so for each experimental data point at a given photon-proton centre-of-mass energy $W$, we can associate a skewness $\xi \approx M_{J/\psi}^2/(2W^2)$ and the value $x = X + \xi \approx 2\xi$. It is in this sense that Eq. (8) should be understood.

Now for each data point, $i$, we calculate the predicted $J/\psi$ cross section $\sigma_+(\text{fit})_i$ using the gluon GPDs obtained via the Shuvaev transform of our gluon PDF fit in Eq. (3) and compare it to the experimental value $\sigma_+(\text{data})_i$.<fnref n="6" /> The ratio $\sigma_+(\text{data})_i/\sigma_+(\text{fit})_i$ is approximately equal to the square of the ratio of the gluon densities

$$ \frac{\sigma_+(\text{data})_i}{\sigma_+(\text{fit})_i} = \left( \frac{g_{\text{eff}}(2\xi_i)}{g_{\text{fit}}(2\xi_i)} \right)^2 . \qquad (9) $$

In other words, in this way we can calculate the effective gluon density (PDF) corresponding to $x_i = 2\xi_i$ as

$$ g_{\text{eff}}(x_i, \mu_{\text{opt}}) = g_{\text{fit}}(x_i, \mu_{\text{opt}}) \sqrt{\frac{\sigma_+(\text{data})_i}{\sigma_+(\text{fit})_i}} . \qquad (10) $$

Correspondingly the error

$$ \delta g_{\text{eff}}(x_i) = \frac{1}{2} g_{\text{eff}}(x_i) \frac{\delta \sigma_+(\text{data})_i}{\sigma_+(\text{data})_i} . \qquad (11) $$

The error of the lower-energy contribution, $\sigma_-(x_i)$, is accounted for in the experimental value of $\sigma_+(\text{data})_i$. Recall that this lower-energy contribution is small and its errors are negligible.

We therefore propose to include these “effective” gluon PDF data points calculated according to Eqs. (10) and (11) in a parton fit analysis in order to reduce the present uncertainties in the behaviour of the low-$x$ partons. Since the final parton distribution may differ from Eq. (3) obtained from [22], to achieve greater accuracy in the final determined values of $xg_{\text{eff}}$, several additional iterations in the $xg_{\text{eff}}$ data may be needed. We therefore proceed as follows. We first calculate the cross section $\sigma_+(\text{fit})_i$ using the values of $x_i g_{\text{fit}}(x_i)$ in Eq. (3) to obtain starting values of $x_i g_{\text{eff}}(x_i)$. Such values are then used in our parton fit analysis, however they may not necessarily be sufficiently accurate to provide to the community for use in other global fits. Consequently, the outputted, new gluon PDFs $xg_{\text{fit}}$ from our parton fit analysis are then passed to our Shuvaev-integral transform routine to

obtain updated GPDs via Eq. (5) and subprocess cross sections $\sigma_+(\text{fit})$ via Eq. (7). These are then used to re-calculate $xg_{\text{eff}}$ via Eq. (10) and used in a second parton fit analysis. The procedure continues until the resulting $xg_{\text{eff}}$ values from parton analysis $n$ do not differ from those in analysis $n-1$ by $\pm 1\sigma$. Fortunately, such a procedure is well convergent. The final determined iterated data points are given in Table 1.<fnref n="7" />

In Table 2, we show the final effective iterated gluon points for exclusive $\Upsilon$ production [6] at the corresponding ‘optimal’ scale 22.4 GeV$^2$. The accuracy of these points is worse than that obtained from $J/\psi$ production, but they allow for a constraint on the $Q^2$-evolution in the low-$x$ region. So far the exclusive $\Upsilon$ data for the $\gamma p \rightarrow \Upsilon p$ cross section from the LHCb Collaboration have been obtained using gap survival factors $S^2(W_\pm)$ and photon fluxes $dn/dk_\pm$ that need updating or correcting, see [29]. Nevertheless, the resulting $\Upsilon$ data shown in Table 2 already include the corresponding corrections, which however never exceed the experimental error bars. Note that the resulting $\Upsilon$ points already show evidence of a rising gluon PDF power behaviour in the region $10^{-5} < x < 10^{-4}$ and $Q^2 = m_b^2 \approx 22 \text{ GeV}^2$ (see Fig. 2 of [29]). Recall that as was shown in [29], at the border of the LHCb acceptance, the choice of survival factor and photon flux only amounted to 5% for $J/\psi$ production, while for $\Upsilon$ production it is of the order of 25%.

### 3 The impact of the low-$x$ gluon effective data points on a global analysis

The gluon parton distribution resulting from a global xFitter [15] NLO analysis of pure Deep-Inelastic-Scattering (DIS) data (red) and DIS data supplemented with the effective gluon points extracted from the exclusive $J/\psi$ and $\Upsilon$ production at the LHC (blue) are shown in Fig. 4. Here we have used the standard xFitter ensemble of DIS data and the standard xFitter ansatz for the input parton parametrizations with the kinematic cut $Q^2 > 2.4 \text{ GeV}^2$. The datasets used and the partial and total minimum $\chi^2$ figure of merit per degree of freedom, $\chi^2_{\text{min}}/\text{d.o.f}$, for both fits are shown in Table 3.

It is clearly seen in Fig. 4 that the inclusion of the low-$x$ effective gluon points extracted from the exclusive heavy

<footnote n="7">⁷ Note that we do not fit to the individual $\sigma(\gamma + p \rightarrow J/\psi + p)$ data points in our approach. Instead, we use effective $xg_{\text{eff}}(x)$ points. Recall the values of $xg_{\text{fit}}(x)$ were obtained from a previous fit to the $J/\psi$ cross section only, and thus may not be sufficiently accurate for a global analysis. Therefore we do not use the values of $xg_{\text{fit}}(x)$ from this former $J/\psi$ fit, but instead fit to the corrected $g_{\text{eff}}(x) = g_{\text{fit}}(x)\sqrt{\sigma(\text{data})(x)/\sigma(\text{fit})(x)}$ data points. These effective values will provide a better overall description of the data and a better stability, so that the true values will not be changed too much in the fit iteration.</footnote>

<footnote n="6">⁶ For the results presented later in Tables 1 and 2, the Shuvaev transform is used in the region $x \le 10^{-4}$.</footnote>
<table><caption>**Table 1** The starting values of $xg_{\text{fit}}$ and the final values of $xg_{\text{eff}}$ and $\delta xg_{\text{eff}}$ obtained from eqns. (10) and (11) at $Q^2 = M^2_{J/\psi}/4 = 2.4 \text{ GeV}^2$ using the exclusive $J/\psi$ production data at 7 TeV [4] and 13 TeV [5] obtained by the LHCb Collaboration</caption><thead><tr><th>13 TeV</th><th>x</th><th>xg<sub>fit</sub></th><th>xg<sub>eff</sub> &plusmn; &delta;xg<sub>eff</sub></th><th>7 TeV</th><th>x</th><th>xg<sub>fit</sub></th><th>xg<sub>eff</sub> &plusmn; &delta;xg<sub>eff</sub></th></tr></thead><tbody><tr><td></td><td>2.84 &times; 10<sup>-5</sup></td><td>3.71</td><td>3.83 &plusmn; 0.26</td><td></td><td>5.28 &times; 10<sup>-5</sup></td><td>3.41</td><td>3.45 &plusmn; 0.24</td></tr><tr><td></td><td>2.21 &times; 10<sup>-5</sup></td><td>3.83</td><td>4.08 &plusmn; 0.17</td><td></td><td>4.11 &times; 10<sup>-5</sup></td><td>3.52</td><td>3.65 &plusmn; 0.20</td></tr><tr><td></td><td>1.72 &times; 10<sup>-5</sup></td><td>3.97</td><td>3.96 &plusmn; 0.14</td><td></td><td>3.20 &times; 10<sup>-5</sup></td><td>3.64</td><td>3.68 &plusmn; 0.19</td></tr><tr><td></td><td>1.34 &times; 10<sup>-5</sup></td><td>4.10</td><td>4.05 &plusmn; 0.13</td><td></td><td>2.49 &times; 10<sup>-5</sup></td><td>3.77</td><td>3.67 &plusmn; 0.19</td></tr><tr><td></td><td>1.04 &times; 10<sup>-5</sup></td><td>4.25</td><td>4.14 &plusmn; 0.12</td><td></td><td>1.94 &times; 10<sup>-5</sup></td><td>3.90</td><td>3.79 &plusmn; 0.19</td></tr><tr><td></td><td>8.14 &times; 10<sup>-6</sup></td><td>4.39</td><td>4.32 &plusmn; 0.13</td><td></td><td>1.51 &times; 10<sup>-5</sup></td><td>4.04</td><td>3.88 &plusmn; 0.19</td></tr><tr><td></td><td>6.34 &times; 10<sup>-6</sup></td><td>4.54</td><td>4.35 &plusmn; 0.13</td><td></td><td>1.18 &times; 10<sup>-5</sup></td><td>4.18</td><td>3.99 &plusmn; 0.20</td></tr><tr><td></td><td>4.93 &times; 10<sup>-6</sup></td><td>4.70</td><td>4.34 &plusmn; 0.14</td><td></td><td>9.18 &times; 10<sup>-6</sup></td><td>4.32</td><td>4.19 &plusmn; 0.20</td></tr><tr><td></td><td>3.84 &times; 10<sup>-6</sup></td><td>4.86</td><td>4.60 &plusmn; 0.17</td><td></td><td>7.15 &times; 10<sup>-6</sup></td><td>4.47</td><td>4.59 &plusmn; 0.23</td></tr><tr><td></td><td>2.99 &times; 10<sup>-6</sup></td><td>5.03</td><td>4.91 &plusmn; 0.27</td><td></td><td>5.57 &times; 10<sup>-6</sup></td><td>4.62</td><td>4.84 &plusmn; 0.28</td></tr></tbody></table>

<table><caption>**Table 2** The starting values of $xg_{\text{fit}}$ and the final values of $xg_{\text{eff}}$ and $\delta xg_{\text{eff}}$ obtained from eqns. (10) and (11) at $Q^2 = M^2_{\Upsilon}/4 = 22.4 \text{ GeV}^2$ using the exclusive $\Upsilon$ production data at 7, 8 TeV [6] obtained by the LHCb Collaboration</caption><thead><tr><th>7, 8 TeV</th><th>x</th><th>xg<sub>fit</sub></th><th>xg<sub>eff</sub> &plusmn; &delta;xg<sub>eff</sub></th></tr></thead><tbody><tr><td></td><td>1.01 &times; 10<sup>-4</sup></td><td>23.5</td><td>19.9 &plusmn; 3.1</td></tr><tr><td></td><td>4.77 &times; 10<sup>-5</sup></td><td>28.3</td><td>23.9 &plusmn; 3.3</td></tr><tr><td></td><td>2.26 &times; 10<sup>-5</sup></td><td>34.1</td><td>31.1 &plusmn; 4.1</td></tr></tbody></table>

vector-meson production at the LHC essentially improve the accuracy of the global parton analysis, and we obtain a good overall final $\chi^2_{\text{min}}/\text{d.o.f} \sim O(1)$. The blue error band is much smaller than the red one. Moreover, the inclusion of new information obtained from the heavy vector-meson ultraperipheral production affects (albeit expectedly not so strongly) even the quark distribution as demonstrated, for example, in Fig. 5. The new $g_{\text{eff}}$ low-x points affect the quark distribution even at a rather large $x$ due to the energy-momentum sum rules and the particular choice of the fixed form of the input ansatz. In Fig. 6, we show the ratio of our fitted ‘DIS + effective gluon points’ gluon and quark singlet PDFs to the central value of the DIS-only ones at a much larger scale, $Q^2 = M^2_Z$. Since during the $Q^2$ evolution the value of $x$ decreases, the parton densities at such large $Q^2$ originate mainly from the relatively large-$x$ input PDF, where the distribution is already well fixed by DIS data. Therefore, the effect of adding the new $J/\psi$ data is less pronounced at $Q^2 = M^2_Z$, but still visible, as shown.

As usual, additional theoretical uncertainties such as the correction to the Non-Relativistic QCD (NRQCD) approach used for the $q\bar{q} \to V$ transition in the computation of $\sigma(\gamma p \to V + p)$, or possible higher-twist effects, or effects of absorptive corrections, or other higher-order corrections, are not shown. Note that, apart from the absorptive effects, all other corrections (e.g. NRQCD) do not affect the energy behaviour of the amplitude, but just change its normalisation. Indeed, to ‘know’ the initial energy we need an additional

gluon, which connects the upper (photon block in Fig. 1) with the lower/target parts of the diagram. This would correspond to a higher-twist contribution and should be suppressed by a large factorisation scale $\mu_F$. Typically, these corrections act ‘locally’ in rapidity space and therefore do not affect the energy behaviour of the amplitude, only its normalisation. On the other hand, our approach, in its present form, satisfactorily describes the available HERA data [22]. Recall that, as was shown in [30], normalising the vector-meson wave function to the experimentally measured $J/\psi \to l^+l^-$ decay width allows us to account for relativistic corrections in the NRQCD approach with good accuracy (within a few percent).

## 4 Working to NNLO via a K factor

In principle, the method of “effective gluon points” described above can be extended to NNLO. Unfortunately at present the coefficient functions for $J/\psi$ and $\Upsilon$ photo- and electroproduction are known only at NLO level [27,31,32].

However, it is possible to include these points in a NNLO global analysis in an approximate way by introducing a NNLO/NLO K factor, whose value can be extracted from the description of existing HERA $\gamma + p \to J/\psi + p$ (and/or $\gamma + p \to \Upsilon + p$) data. In the following, we denote this analysis NNLO* to distinguish it from an analysis in which the full NNLO coefficient functions would be used when they become available. In general the K factor may depend on the factorisation scale and on the ratio $z = (X + \xi)/2\xi$ of the parton momentum fractions $x = (X + \xi)$ to $\xi$. That is the K factor should be included in the convolution for the $\gamma + p \to J/\psi + p$ amplitude

$$ \begin{aligned} & \mathcal{M}(\gamma + p \to J/\psi + p) \\ &= \sum_{i=g,q} N_i \int_{-1}^{1} \frac{dX}{X} F_i(X, \xi) K C_i^{\text{NLO}}(\xi/X), \end{aligned} \quad (12) $$
<table><caption>**Table 3** The partial $\chi^2_{min}/\text{d.o.f}$ for each dataset included in the baseline fit (DIS) and with the effective gluon data points added (DIS + eff. gluon pts.). Here, we use the DIS data in the kinematic range $Q^2 > 2.4 \text{ GeV}^2$. The total $\chi^2_{min}/\text{d.o.f}$ is also given</caption><thead><tr><th>Dataset</th><th>χ<sup>2</sup><sub>min</sub>/d.o.f (DIS)</th><th>χ<sup>2</sup><sub>min</sub>/d.o.f (DIS+eff. gluon pts.)</th></tr></thead><tbody><tr><td>HERA1+2 NCep 820</td><td>80/73</td><td>79/73</td></tr><tr><td>HERA1+2 NCep 460</td><td>220/207</td><td>220/207</td></tr><tr><td>HERA1+2 CCep</td><td>43/39</td><td>44/39</td></tr><tr><td>HERA1+2 NCem</td><td>221/159</td><td>220/159</td></tr><tr><td>HERA1+2 CCem</td><td>54/42</td><td>56/42</td></tr><tr><td>HERA1+2 NCep 575</td><td>223/257</td><td>227/257</td></tr><tr><td>HERA1+2 NCep 920</td><td>465/391</td><td>470/391</td></tr><tr><td>LHC excl. J/ψ pp 7 TeV</td><td>N/A</td><td>8.95/10</td></tr><tr><td>LHC excl. J/ψ pp 13 TeV</td><td>N/A</td><td>3.51/10</td></tr><tr><td>LHC excl. γ pp 7,8 TeV</td><td>N/A</td><td>3.23/3</td></tr><tr><td>Total χ<sup>2</sup><sub>min</sub>/d.o.f</td><td>1412/1154 ≈ 1.22</td><td>1444/1177 ≈ 1.23</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 4** The gluon distributions $xg$ at scales $\mu_F^2 = Q^2 = 2.4 \text{ GeV}^2$ (left) and $22 \text{ GeV}^2$ (right) obtained in xFilter fitting DIS (red) and DIS+effective gluon data points (blue)</figcaption></figure>

given previously. The integral in Eq. (12) converges<fnref n="8" /> at $X \sim \xi$. This means that actually after the convolution of the coefficient function with the gluon distribution the effective NNLO/NLO K factor may be considered as a constant assuming that the power $\lambda$ of Eq. (3) does not vary in the interval of interest.

Thus the value of this K factor can be estimated by describing the exclusive $J/\psi$ photoproduction data from HERA with the NNLO partons and NLO coefficient functions in the region $0.01 > x > 0.001$ since here the uncertainties in the gluon distribution given by the present global analysis are relatively small. In other words the square of the K factor can be obtained as the ratio of the measured $\gamma + p \rightarrow J/\psi + p$ cross section, $\sigma(\text{data})$, to $\sigma^{\text{NNLO/NNLO}}$ calculated using the NLO

coefficient function but NNLO input partons (PDF/GPD). That is,

$$K_i = \sqrt{\left(\frac{\sigma(\text{data})}{\sigma^{\text{NNLO/NNLO}}}\right)_i}, \qquad (13)$$

and then make the trivial average of the K factors for each data point. From Fig. 3 we may hope that the power $\lambda$ which describes the energy behaviour of the $\gamma + p \rightarrow J/\psi + p$ cross section is more or less constant. However, as it is seen from Eq. (12), this would require the pure power behaviour for the gluons and not for the $J/\psi$ cross section. Unfortunately in the present analysis just in the region of interest $0.01 > x > 0.001$ we observe that the gluon PDF is too flat (see Fig. 4 left). The power growth starts only at $x < 0.001$, while it is well known that the $J/\psi$ photoproduction data from HERA are all well described by a power behaviour starting at $x < 0.01$ [19,21]. Most probably this is explained by the fact that the DGLAP evolution does not account for the (higher twist) absorptive corrections which are not negligible in this low-$x$ region. To mimic the role of these corrections in DIS data, the fit chooses gluons which slightly decrease at $x < 0.01$.

<footnote n="8">⁸ It is known that the $\gamma + p \rightarrow J/\psi + p$ amplitude increases as a power of the photon energy (i.e. as a power of $1/\xi$). On the other hand the NNLO gluon coefficient function can contain only $\ln^2 z$ from the higher-order gluon loop insertions. Thus at large energy (i.e. very small $\xi$) the dominant contribution to (12) comes from $z \sim O(1)$ (i.e. $X \sim \xi$) while the power growth of the amplitude is provided by the GPD power behaviour. In this case we have to calculate the $\gamma + p \rightarrow J/\psi + p$ NNLO amplitude, $M$, by multiplying the NLO coefficient functions, $C_i^{\text{NNLO}}$ with $i = g, q$, by this K factor.</footnote>
<figure><img src="image_8.png" /><figcaption>**Fig. 6** The ratio of the fitted 'DIS + effective gluon points' gluon $xg$ (left) and quark singlet $x\Sigma$ (right) PDFs to the DIS-only ones, $xg_{\text{ref}}$ and $x\Sigma_{\text{ref}}$, respectively, at $Q^2 = M_Z^2$</figcaption></figure>

<figure><img src="image_9.png" /></figure>

In a future analysis, we have two possibilities. Either to include these absorptive effects into the DGLAP evolution following the GLR-MQ scheme [33,34] (with the possibility that this will generate power increasing gluons already in the HERA domain, $x < 0.01$), or to have at hand the full NNLO coefficient function to reach complete NNLO accuracy and go beyond the NNLO* approach discussed here.

In any case in order to obtain realistic partons at such low $x$ and scales, the absorptive effect should be accounted for.

## 5 Conclusion and outlook

We have described the method of “effective gluon points” which readily allows for the inclusion of exclusive heavy vector-meson $V$ photoproduction data, where $V = J/\psi$, $\Upsilon$ (as well as the exclusive production in ultraperipheral events at the LHC) in a conventional global parton analysis. Using xFitter, we fit DIS data together with these effective gluon points extracted from the $J/\psi$ and $\Upsilon$ exclusive production data from LHCb and demonstrate that this crucially improves

the accuracy of the obtained gluon distribution in the low-$x$ domain. The values of $xg_{\text{eff}}(x)$ extracted from the available exclusive $J/\psi$ and $\Upsilon$ data from LHCb are presented in Tables 1 and 2.

Since at present the NNLO coefficient functions for the photon to vector meson ($\gamma \to V$) transition are not available, our computations were performed at NLO level. However, there is the possibility to reach NNLO* accuracy, by extracting the NNLO/NLO $K$ factor from the existing HERA $\gamma + p \to J/\psi + p$ data as discussed in Sect. 4. Moreover, exclusive vector-meson production data have already been collected by the LHC and more should follow. This merits further theoretical study of the NNLO* analysis, which we reserve for future work.

The analysis presented in this work has demonstrated the potential of existing exclusive data to significantly improve the NLO global parton distributions at small $x$.

**Acknowledgements** This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No. 824093 in order to contribute to the EU Virtual Access "NLOAccess". This project has also received funding from the
Agence Nationale de la Recherche (ANR) via the grant ANR-20-CE31-0015 ("PrecisOnium") and via the IDEX Paris-Saclay "Investissements d'Avenir" (ANR-11-IDEX-0003-01) through the GLUODYNAMICS project funded by the "P2IO LabEx (ANR-10-LABX-0038)". This work was also partly supported by the French CNRS via the IN2P3 Master Project "QCDFactorisation@NLO". The work of TT was supported by the STFC Consolidated Grant ST/T000988/1 and currently by ST/X000699/1.

**Data Availability Statement** Data will be made available on reasonable request. [Author's comment: The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.]

**Code Availability Statement** Code/software will be made available on reasonable request. [Author's comment: The code/software generated during and/or analysed during the current study is available from the corresponding author on reasonable request.]

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Funded by SCOAP³.

## References

1. S. Bailey, T. Cridge, L.A. Harland-Lang, A.D. Martin, R.S. Thorne, Eur. Phys. J. C 81(4), 341 (2021). [arXiv:2012.04684 [hep-ph]]

2. R.D. Ball et al. [NNPDF], Eur. Phys. J. C 82(5), 428. [arXiv:2109.02653 [hep-ph]]

3. T.J. Hou, J. Gao, T.J. Hobbs, K. Xie, S. Dulat, M. Guzzi, J. Huston, P. Nadolsky, J. Pumplin, C. Schmidt et al., Phys. Rev. D 103(1), 014013 (2021). [arXiv:1912.10053 [hep-ph]]

4. R. Aaij et al. [LHCb], J. Phys. G 41, 055002 (2014). [arXiv:1401.3288 [hep-ex]]

5. R. Aaij et al. [LHCb], JHEP 10, 167 (2018). [arXiv:1806.04079 [hep-ex]]

6. R. Aaij et al. [LHCb], JHEP 09, 084 (2015). [arXiv:1505.08139 [hep-ex]]

7. M. Diehl, Phys. Rep. 388, 41–277 (2003). [arXiv:hep-ph/0307382]

8. A.G. Shuvaev, K.J. Golec-Biernat, A.D. Martin, M.G. Ryskin, Phys. Rev. D 60, 014015 (1999). [arXiv:hep-ph/9902410]

9. A. Shuvaev, Phys. Rev. D 60, 116005 (1999). [arXiv:hep-ph/9902318]

10. T. Ohrndorf, Nucl. Phys. B 198, 26–44 (1982)

11. X.D. Ji, J. Phys. G 24, 1181–1205 (1998). [arXiv:hep-ph/9807358]

12. T. Kluge, K. Rabbertz and M. Wobisch, [arXiv:hep-ph/0609285]

13. T. Carli, D. Clements, A. Cooper-Sarkar, C. Gwenlan, G.P. Salam, F. Siegert, P. Starovoitov, M. Sutton, Eur. Phys. J. C 66, 503–524 (2010). [arXiv:0911.2985 [hep-ph]]

14. S. Carrazza, E.R. Nocera, C. Schwan, M. Zaro, JHEP 12, 108 (2020). [arXiv:2008.12789 [hep-ph]]

15. S. Alekhin, O. Behnke, P. Belov, S. Borroni, M. Botje, D. Britzger, S. Camarda, A.M. Cooper-Sarkar, K. Daum, C. Diaconu et al., Eur. Phys. J. C 75(7), 304 (2015). [arXiv:1410.4412 [hep-ph]]

16. V.A. Khoze, A.D. Martin, M.G. Ryskin, Eur. Phys. J. C 74(2), 2756 (2014). [arXiv:1312.3851 [hep-ph]]

17. V.M. Budnev, I.F. Ginzburg, G.V. Meledin, V.G. Serbo, Phys. Rep. 15, 181–281 (1975)

18. S. Chekanov et al. [ZEUS], Eur. Phys. J. C 24, 345–360 (2002). [arXiv:hep-ex/0201043]

19. S. Chekanov et al. [ZEUS], Nucl. Phys. B 695, 3–37 (2004). [arXiv:hep-ex/0404008]

20. A. Aktas et al. [H1], Eur. Phys. J. C 46, 585–603 (2006). [arXiv:hep-ex/0510016]

21. C. Alexa et al. [H1], Eur. Phys. J. C 73(6), 2466. [arXiv:1304.5162 [hep-ex]]

22. C.A. Flett, A.D. Martin, M.G. Ryskin, T. Teubner, Phys. Rev. D 102, 114021 (2020). [arXiv:2006.13857 [hep-ph]]

23. V.A. Khoze, A.D. Martin, M.G. Ryskin, Eur. Phys. J. C 73, 2503 (2013). [arXiv:1306.2149 [hep-ph]]

24. S.P. Jones, A.D. Martin, M.G. Ryskin, T. Teubner, J. Phys. G 43(3), 035002 (2016). [arXiv:1507.06942 [hep-ph]]

25. E.G. de Oliveira, A.D. Martin, M.G. Ryskin, Eur. Phys. J. C 72, 2069 (2012). [arXiv:1205.6108 [hep-ph]]

26. S.P. Jones, A.D. Martin, M.G. Ryskin, T. Teubner, Eur. Phys. J. C 76(11), 633 (2016). [arXiv:1610.02272 [hep-ph]]

27. D. Y. Ivanov, A. Schafer, L. Szymanowski, G. Krasnikov, Eur. Phys. J. C 34(3), 297–316 (2004) [Erratum: Eur. Phys. J. C 75 (2015) no.2, 75]. [arXiv:hep-ph/0401131]

28. C.A. Flett, S.P. Jones, A.D. Martin, M.G. Ryskin, T. Teubner, Phys. Rev. D 101(9), 094011 (2020). [arXiv:1908.08398 [hep-ph]]

29. C.A. Flett, S.P. Jones, A.D. Martin, M.G. Ryskin, T. Teubner, Phys. Rev. D 105(3), 034008 (2022). [arXiv:2110.15575 [hep-ph]]

30. P. Hoodbhoy, Phys. Rev. D 56, 388–393 (1997). [arXiv:hep-ph/9611207]

31. S.P. Jones, Ph.D. thesis, University of Liverpool (2014) (unpublished)

32. C.A. Flett, J.A. Gracey, S.P. Jones, T. Teubner, JHEP 08, 150 (2021). [arXiv:2105.07657 [hep-ph]]

33. L.V. Gribov, E.M. Levin, M.G. Ryskin, Phys. Rep. 100, 1–150 (1983)

34. A.H. Mueller, J.W. Qiu, Nucl. Phys. B 268, 427–452 (1986)