IPPP/24/41 LTH 1381

January 29, 2025

## A method to include exclusive heavy vector-meson production data at small x in global parton analyses

C.A. Flett , A.D. Martin , M.G. Ryskin a b c and T. Teubner d

a Universit´ Paris-Saclay, CNRS, IJCLab, 91405 Orsay, France e 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 d Department of Mathematical Sciences, University of Liverpool, Liverpool, L69 3BX, U.K.

## 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 within xFitter to determine the gluon parton distribution function (PDF) at nextto-leading order (NLO) at moderate-to-low x using the measurements of exclusive J/ψ production in ep and pp collisions from HERA and LHC. We further study the constraints from the corresponding Υ 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, 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. For x &gt; ∼ 10 -3 , there is a general agreement among the results from various global fit analyses [1, 2, 3]. However, 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.

Figure 1: (Left) LO contribution to γp → V + p , where the vector meson V = J/ψ, Υ . (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 NLO g and GPD F g , dominates the NLO quark contribution. In these diagrams, the momentum P ≡ ( p + p ′ ) / 2 and l is the loop momentum. Note that the momentum fractions of the left and right partons are x = X + ξ and x ′ = X -ξ respectively; for the gluons connected to the heavy quark-antiquark pair, we have x ′ ≪ x and so x ≃ 2 ξ .

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01128v2-figure-bd432af4-66f2-47be-817a-722b09eb2f7f.jpg)

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/ψ and Υ [4, 5, 6] have allowed for the determination of the gluon Parton Distribution Function (PDF) down to x ∼ 3 × 10 -6 at factorization scales µ 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 ≈ 2 ξ at very small x , where 2 ξ 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 ( ∼ O x ( )) accuracy using the Shuvaev transform [8]. The Shuvaev transform makes use of the fact that as ξ → 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. [10]) even for ξ = 0 the Gegenbauer

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 [9]. As ξ → 0 the Gegenbauer moments become equal to the Mellin moments.

moments can be obtained from the Mellin moments to O ξ ( ) accuracy. Thus it is possible to obtain the full GPD function at small ξ from its known moments.

In principle, this allows us to include the lowx exclusive J/ψ data in the global PDF analysis. In practice the problem is that the Shuvaev transform amounts to a slowly convergent double integral and the corresponding computation is time consuming. This is troublesome for the global fit as after every iteration, a new three-dimensional grid defined over the variables X,ξ and Q 2 has to be computed to obtain the updated theory prediction.

In the present paper we describe a method which helps to 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/ψ and Υ exclusive cross sections affect the result we compare the standard xFitter [11] 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 lowx gluon PDF data

To overcome this difficulty, we propose the following procedure. We translate the experimental J/ψ 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 \rightarrow p + J / \psi + p$$

contains two components,

$$\frac { d \sigma ^ { ( p p \to p + J / \psi + p ) } } { d Y } \ = \ S ^ { 2 } ( W _ { + } ) \left ( k _ { + } \frac { d n } { d k _ { + } } \right ) \sigma _ { + } ( \gamma p \to J / \psi p ) \ + \ W _ { - } \ \text{term},$$

corresponding to energies W ± , 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 σ γp ( → J/ψp ), for each rapidity Y of the J/ψ , are weighted by the corresponding gap survival factors S 2 ( W ± ) [12] and photon fluxes dn/dk ± [13]. The cross section with the lower γp -energy ( W -) was measured at HERA [14, 15, 16, 17] 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/ψ data provide reliable values for σ + ( γ + p → J/ψ + p ), corresponding to the component with the larger γp -energy, W + , see Fig. 3. 2

These σ + data were well described in [18] by a power-like lowx gluon

$$x g ( x ) = A \cdot x ^ { - \lambda } \quad \text{ with } \quad \lambda = 0. 1 3 5 \pm 0. 0 0 6,$$

Figure 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.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01128v2-figure-d543ec03-8186-4229-a5ac-b92807227415.jpg)

also shown in Fig. 3.

Moreover it was shown that after the optimal factorisation scale , 3 µ F = M J/ψ / 2, was chosen, where M J/ψ is the mass of the J/ψ , and the so-called Q 0 subtraction [21] 4 was performed, the quark contribution to σ + becomes negligible [23]. That is the value of

$$\sigma _ { + } \, \in \, \left ( \, 2 \xi \, g ( 2 \xi ) \, \right ) ^ { 2 }$$

is just proportional to the gluon density squared. 5

Now for each data point, i , we calculate the predicted J/ψ cross section σ + (fit) i based on the gluons from our fit and compare it with the experimental value σ + (data) . i The ratio σ + (data) i /σ + (fit) i is 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 }.$$

In other words, in this way we can calculate the effective gluon density (PDF) corresponding to x i = 2 ξ 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 } } }.$$

2 Actually, the subtraction of the W -contribution was done already in the LHCb analyses where the results are presented as the σ + ( γ + p → V + ) cross sections. p

3 At this optimal scale the large double-logarithmic contributions of the form c n ( α s ln(1 /ξ ) ln µ 2 F ) n are resummed into the incoming PDF/GPD [20].

4 The Q 0 subtraction is necessary to exclude the double counting between the low l T part of the coefficient function and the PDFinput . This, together with the resummation of the aforementioned double-logarithmic terms [20], helps mitigate the problem of the strong factorisation scale dependence [22] of the NLO amplitude.

5 Strictly speaking, the value of σ + is proportional to the square of the corresponding GPD.

Figure 3: The description of the J/ψ photoproduction HERA [14, 15, 16, 17] and LHCb [4, 5] data based on using the central value of the global gluon PDF from the MMHT14 global analysis for x &gt; 10 -3 and a fitted power-like gluon PDF, see eqn. (3) , for x &lt; 10 -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.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01128v2-figure-4e19a076-3256-4296-ad74-ab4734032da5.jpg)

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 } }.$$

The error of the lower energy contribution, σ -( x i ), is accounted for in the experimental value of σ + (data) . i Recall that this lower energy contribution is small and its errors are negligible.

To summarize: we propose to include these 'effective' gluon PDF data points calculated according to eqns. (6), (7) in our parton analysis in order to reduce the present uncertainties in the behaviour of the lowx partons.

Since the final parton distribution may differ from eqn. (3) obtained from [18], to achieve greater accuracy several iterations may be needed such that the resulting output parameters in iteration i do not differ from those in iteration i -1 by ± 1 σ . The final determined iterated data points are given in Table 1. 6

6 Note that we do not fit to the individual σ γ ( + p → J/ψ + ) data points in our approach. p Instead, we use effective xg eff ( x ) points. Recall the values of xg fit ( x ) were obtained from a previous fit to the J/ψ cross section only, and thus may not be sufficiently accurate for a global analysis. Therefore we do not use the values

Table 1: The values of xg eff and δxg eff obtained from eqns. (6) and (7) at Q 2 = M 2 J/ψ / 4 = 2 4 GeV . 2 using the exclusive J/ψ production data at 7 TeV [4] and 13 TeV [5] obtained by the LHCb collaboration.

| 13 TeV   | x               |   xg eff |   δxg eff |
|----------|-----------------|----------|-----------|
|          | 2 . 84 × 10 - 5 |     3.83 |      0.26 |
|          | 2 . 21 × 10 - 5 |     4.08 |      0.17 |
|          | 1 . 72 × 10 - 5 |     3.96 |      0.14 |
|          | 1 . 34 × 10 - 5 |     4.05 |      0.13 |
|          | 1 . 04 × 10 - 5 |     4.14 |      0.12 |
|          | 8 . 14 × 10 - 6 |     4.32 |      0.13 |
|          | 6 . 34 × 10 - 6 |     4.35 |      0.13 |
|          | 4 . 93 × 10 - 6 |     4.34 |      0.14 |
|          | 3 . 84 × 10 - 6 |     4.6  |      0.17 |
|          | 2 . 99 × 10 - 6 |     4.91 |      0.27 |
| 7 TeV    | 5 . 28 × 10 - 5 |     3.45 |      0.24 |
|          | 4 . 11 × 10 - 5 |     3.65 |      0.2  |
|          | 3 . 20 × 10 - 5 |     3.68 |      0.19 |
|          | 2 . 49 × 10 - 5 |     3.67 |      0.19 |
|          | 1 . 94 × 10 - 5 |     3.79 |      0.19 |
|          | 1 . 51 × 10 - 5 |     3.88 |      0.19 |
|          | 1 . 18 × 10 - 5 |     3.99 |      0.2  |
|          | 9 . 18 × 10 - 6 |     4.19 |      0.2  |
|          | 7 . 15 × 10 - 6 |     4.59 |      0.23 |
|          | 5 . 57 × 10 - 6 |     4.84 |      0.28 |

Table 2: The values of xg eff and δxg eff obtained from eqns. (6) and (7) at Q 2 = M / 2 Υ 4 = 22 4 GeV . 2 using the exclusive Υ production data at 7, 8 TeV [6] obtained by the LHCb collaboration.

| 7, 8 TeV   | x               |   xg eff |   δxg eff |
|------------|-----------------|----------|-----------|
|            | 1 . 01 × 10 - 4 |     19.9 |       3.1 |
|            | 4 . 77 × 10 - 5 |     23.9 |       3.3 |
|            | 2 . 26 × 10 - 5 |     31.1 |       4.1 |

In Table 2, we show the final effective iterated gluon points for exclusive Υ production [6] at the corresponding 'optimal' scale 22.4 GeV . 2 The accuracy of these points is worse than

of xg fit ( x ) from this former J/ψ fit, but instead fit to the corrected g eff ( x ) = g fit ( x ) √ σ (data)( x /σ ) (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.

that obtained from J/ψ production, but they allow for a constraint on the Q 2 -evolution in the lowx region. 7

For the results presented in Tables 1 &amp; 2, we actually use the Shuvaev transform for x ≤ 10 -4 ( x ≤ 10 -3 for the data presented in Fig. 3).

## 3 The impact of the lowx gluon effective data points on a global analysis

The gluon parton distribution resulting from a global xFitter [11] NLO analysis of pure DeepInelastic-Scattering (DIS) data (red) and DIS data supplemented with the effective gluon points extracted from the exclusive J/ψ and Υ 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 &gt; 2 2 4 GeV . NLO coefficient . 2 functions and NLO PDFs were used to obtain the g fit values. The datasets used and the partial and total minimum χ 2 figure of merit per degree of freedom, χ 2 min / d.o.f, for both fits are shown in Table 3.

It is clearly seen in Fig. 4 that the inclusion of the lowx effective gluon points extracted from the exclusive heavy vector-meson production at the LHC essentially improve the accuracy of the global parton analysis, and we obtain a good overall final χ 2 min / d.o.f ∼ O (1). The blue error band is much smaller than the red one. 8 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 eff lowx 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.

7 So far the exclusive Υ data for the γp → Υ p cross section from the LHCb collaboration have been obtained using gap survival factors S 2 ( W ± ) and photon fluxes dn/dk ± that need updating or correcting, see [24]. Nevertheless, the resulting Υ data shown in Table 2 already include the corresponding corrections, which however never exceed the experimental error bars. Note that the resulting Υ points already show evidence of a rising gluon PDF power behaviour in the region 10 -5 &lt; x &lt; 10 -4 and Q 2 = m 2 b ∼ 22 GeV 2 (see Fig. 2 of [24]). Recall that as was shown in [24], at the border of the LHCb acceptance, the choice of survival factor and photon flux only amounted to 5% for J/ψ production, while for Υ production it is of the order of 25%.

8 As usual, additional theoretical uncertainties such as the correction to the Non-Relativistic QCD (NRQCD) approach used for the qq ¯ → V transition in the computation of σ γp ( → V + ), or possible higher twist effects, p 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 µ 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 [18]. Recall that, as was shown in [19], normalising the vector-meson wave function to the experimentally measured J/ψ → l + l -decay width allows us to account for relativistic corrections in the NRQCD approach with good accuracy (within a few percent).

| Dataset                 | χ 2 min / d.o.f (DIS)   | χ 2 min / d.o.f (DIS+eff. gluon pts.)   |
|-------------------------|-------------------------|-----------------------------------------|
| HERA1+2 NCep 820        | 80/73                   | 79/73                                   |
| HERA1+2 NCep 460        | 220/207                 | 220/207                                 |
| HERA1+2 CCep            | 43/39                   | 44/39                                   |
| HERA1+2 NCem            | 221/159                 | 220/159                                 |
| HERA1+2 CCem            | 54/42                   | 56/42                                   |
| HERA1+2 NCep 575        | 223/257                 | 227/257                                 |
| HERA1+2 NCep 920        | 465/391                 | 470/391                                 |
| LHC excl. J/ψ pp 7 TeV  | N/A                     | 8.95/10                                 |
| LHC excl. J/ψ pp 13 TeV | N/A                     | 3.51/10                                 |
| LHC excl. Υ pp 7,8 TeV  | N/A                     | 3.23/3                                  |
| Total χ 2 min / d.o.f   | 1412/1154 ≈ 1.22        | 1444/1177 ≈ 1.23                        |

Table 3: The partial χ 2 min / 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 &gt; 2 2 4 . GeV . The total 2 χ 2 min / d.o.f is also given.

## 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/ψ and Υ photo- and electroproduction are known only at NLO level [22, 25, 26].

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 γ + p → J/ψ + p (and/or γ + p → Υ+ ) data. In the following, we denote p 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 factorization scale and on the ratio z = ( X + ) 2 ξ / ξ of the parton momentum fractions x = ( X + ) to ξ ξ . That is the K factor should be included in the convolution for the γ + p → J/ψ + p amplitude

$$\mathcal { M } ( \gamma + p \to J / \psi + p ) \ = \ \sum _ { i = g, q } N _ { i } \int _ { - 1 } ^ { 1 } \frac { d X } { X } \ F _ { i } ( X, \xi ) \ K \ C _ { i } ^ { \text{NLO} } ( \xi / X )$$

where the constants N i provide the correct normalization, see e.g. [18], and F i and C i denote quark and gluon GPDs and coefficient functions respectively. In our case the factorization scale is fixed, µ F = M J/ψ / 2, and the four-momentum transfer t = t min = (2 ξm p ) 2 / (1 -2 ) ξ ∼ 0, so these parameters are not shown explicitly in eqn. (8). The integral in eqn. (8) converges 9 at

9 It is known that the γ + p → J/ψ + p amplitude increases as a power of the photon energy (i.e. as a power of 1 /ξ ). On the other hand the NNLO gluon coefficient function can contain only ln 2 z from the higher-order

1

-

10

1

-

10

x

x

1

1

0.1

0

10

10

10

10

10

1

5

-

4

-

3

-

2

-

1

-

x

0.1

0

10

10

10

10

10

1

5

-

4

-

3

-

2

-

1

-

x

x

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01128v2-figure-45bec1bf-0589-4e9a-ac4a-1c9f087384db.jpg)

1

Figure 4: The gluon distributions xg at scales µ 2 F = Q 2 = 2.4 GeV 2 (left) and 22 GeV 2 (right) obtained in xFitter fitting DIS (red) and DIS+effective gluon data points (blue).

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01128v2-figure-24a6f697-2630-4182-9de2-07683324db15.jpg)

X ∼ ξ . 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 λ of eqn. (3) does not vary in the interval of interest .

Thus the value of this K factor can be estimated by describing the exclusive J/ψ photoproduction data from HERA with the NNLO partons and NLO coefficient functions in the region 0 01 . &gt; x &gt; 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 γ + p → J/ψ + p cross section, σ (data), to σ NLO 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{NLO} / \text{NNLO} } } \right ) _ { i } } \,,$$

and then make the trivial average of the K factors for each data point. From Fig. 3 we may hope that the power λ which describes the energy behaviour of the γ + p → J/ψ + p cross section is more or less constant. However, as it is seen from eqn. (8), this would require the pure power behaviour for the gluons and not for the J/ψ cross section. Unfortunately in the present analysis just in the region of interest 0 01 . &gt; x &gt; 0 001 we observe that the gluon PDF . is too flat (see Fig. 4 left). The power growth starts only at x &lt; 0 001. . Most probably this gluon loop insertions. Thus at large energy (i.e. very small ξ ) the dominant contribution to (8) comes from z ∼ O (1) (i.e. X ∼ ξ ) while the power growth of the amplitude is provided by the GPD power behaviour. In this case we have to calculate the γ + p → J/ψ + p NNLO amplitude, M , by multiplying the NLO coefficient functions, C i NLO with i = g, q , by this K factor.

)

(x,Q

0.05

0

2

Σ

x

12

10

8

6

4

2

0

-

10

-

10

5

2

Q

DIS

DIS

5

)

xg(x,Q

2

x

0

-

5

-

4

-

3

-

2

-

1

-

5

-

4

-

3

-

2

-

1

-

5

4

-

3

-

2

-

1

-

5

0.1

0

0

0.05

0

10

10 10

10 10

10 10

10 10

10

1

-

4

-

3

-

2

-

1

-

x

x

10 10

10 10

10 10

10 10

10

1

x

1

10

1

x

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01128v2-figure-6f7dcd8b-0023-419e-9909-a2e97831f2ff.jpg)

1

Figure 5: The valence xd v quark distribution at scale 2.4 GeV 2 (left) and the singlet x Σ quark distribution at 22 GeV 2 (right) obtained in xFitter fitting DIS data (red) and DIS+effective gluon data points (blue).

is explained by the fact that the DGLAP evolution does not account for the (higher twist) absorptive corrections which are not negligible in this lowx region. To mimic the role of these corrections in DIS data the fit chooses gluons which slightly decrease at x &lt; 0 01. .

In a future analysis, we have two possibilities. Either to include these absorptive effects into the DGLAP evolution following the GLR-MQ scheme [27, 28] (with the possibility that this will generate power increasing gluons already in the HERA domain, x &lt; 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/ψ, Υ (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/ψ and Υ exclusive production data from LHCb and demonstrate that this crucially improves the accuracy of the obtained gluon distribution in the lowx domain. The values of

xg eff ( x ) extracted from the available exclusive J/ψ and Υ data from LHCb are presented in Tables 1 and 2.

Since at present the NNLO coefficient functions for the photon to vector meson ( γ → 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 γ + p → J/ψ + p data as discussed in Section 4. Moreover, exclusive vectormeson 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 .

## 6 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 ParisSaclay '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.

## References

- [1] S. Bailey, T. Cridge, L. A. Harland-Lang, A. D. Martin and R. S. Thorne, Eur. Phys. J. C 81 (2021) no.4, 341 [arXiv:2012.04684 [hep-ph]].
- [2] R. D. Ball et al. [NNPDF], Eur. Phys. J. C 82 (2022) no.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 and C. Schmidt, et al. Phys. Rev. D 103 (2021) no.1, 014013 [arXiv:1912.10053 [hep-ph]].
- [4] R. Aaij et al. [LHCb], J. Phys. G 41 (2014), 055002 [arXiv:1401.3288 [hep-ex]].
- [5] R. Aaij et al. [LHCb], JHEP 10 (2018), 167 [arXiv:1806.04079 [hep-ex]].
- [6] R. Aaij et al. [LHCb], JHEP 09 (2015), 084 [arXiv:1505.08139 [hep-ex]].
- [7] M. Diehl, Phys. Rept. 388 (2003), 41-277 [arXiv:hep-ph/0307382 [hep-ph]].

- [8] A. G. Shuvaev, K. J. Golec-Biernat, A. D. Martin and M. G. Ryskin, Phys. Rev. D 60 (1999), 014015 [arXiv:hep-ph/9902410 [hep-ph]]; A. Shuvaev, Phys. Rev. D 60 (1999), 116005 [arXiv:hep-ph/9902318 [hep-ph]].
- [9] T. Ohrndorf, Nucl. Phys. B 198 (1982), 26-44.
- [10] X. D. Ji, J. Phys. G 24 (1998), 1181-1205 [arXiv:hep-ph/9807358 [hep-ph]].
- [11] S. Alekhin, O. Behnke, P. Belov, S. Borroni, M. Botje, D. Britzger, S. Camarda, A. M. Cooper-Sarkar, K. Daum and C. Diaconu, et al. Eur. Phys. J. C 75 (2015) no.7, 304 [arXiv:1410.4412 [hep-ph]].
- [12] V. A. Khoze, A. D. Martin and M. G. Ryskin, Eur. Phys. J. C 74 (2014) no.2, 2756 [arXiv:1312.3851 [hep-ph]].
- [13] V. M. Budnev, I. F. Ginzburg, G. V. Meledin and V. G. Serbo, Phys. Rept. 15 (1975), 181-281
- [14] S. Chekanov et al. [ZEUS], Eur. Phys. J. C 24 (2002), 345-360 [arXiv:hep-ex/0201043 [hep-ex]].
- [15] S. Chekanov et al. [ZEUS], Nucl. Phys. B 695 (2004), 3-37 [arXiv:hep-ex/0404008 [hep-ex]].
- [16] A. Aktas et al. [H1], Eur. Phys. J. C 46 (2006), 585-603 [arXiv:hep-ex/0510016 [hep-ex]].
- [17] C. Alexa et al. [H1], Eur. Phys. J. C 73 (2013) no.6, 2466 [arXiv:1304.5162 [hep-ex]].
- [18] C. A. Flett, A. D. Martin, M. G. Ryskin and T. Teubner, Phys. Rev. D 102 (2020), 114021 [arXiv:2006.13857 [hep-ph]].
- [19] P. Hoodbhoy, Phys. Rev. D 56 (1997), 388-393 [arXiv:hep-ph/9611207 [hep-ph]].
- [20] S. P. Jones, A. D. Martin, M. G. Ryskin and T. Teubner, J. Phys. G 43 (2016) no.3, 035002 [arXiv:1507.06942 [hep-ph]]; E. G. de Oliveira, A. D. Martin and M. G. Ryskin, Eur. Phys. J. C 72 (2012), 2069 [arXiv:1205.6108 [hep-ph]].
- [21] S. P. Jones, A. D. Martin, M. G. Ryskin and T. Teubner, Eur. Phys. J. C 76 (2016) no.11, 633 [arXiv:1610.02272 [hep-ph]].
- [22] D. Y. Ivanov, A. Schafer, L. Szymanowski and G. Krasnikov, Eur. Phys. J. C 34 (2004) no.3, 297-316 [erratum: Eur. Phys. J. C 75 (2015) no.2, 75] [arXiv:hep-ph/0401131 [hep-ph]].
- [23] C. A. Flett, S. P. Jones, A. D. Martin, M. G. Ryskin and T. Teubner, Phys. Rev. D 101 (2020) no.9, 094011 [arXiv:1908.08398 [hep-ph]].
- [24] C. A. Flett, S. P. Jones, A. D. Martin, M. G. Ryskin and T. Teubner, Phys. Rev. D 105 (2022) no.3, 034008 [arXiv:2110.15575 [hep-ph]].

- [25] S. P. Jones, PhD thesis, University of Liverpool, 2014, Unpublished.
- [26] C. A. Flett, J. A. Gracey, S. P. Jones and T. Teubner, JHEP 08 (2021), 150 [arXiv:2105.07657 [hep-ph]].
- [27] L. V. Gribov, E. M. Levin and M. G. Ryskin, Phys. Rept. 100 (1983), 1-150.
- [28] A. H. Mueller and J. W. Qiu, Nucl. Phys. B 268 (1986), 427-452.
C. A. Flett A. D. Martin M. G. Ryskin T. Teubner
A method to include exclusive heavy vector-meson production data at small $x$ in global parton analyses
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 within xFitter 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.