NLO QCD corrections to exclusive electroproduction of quarkonium

The process of exclusive electroproduction of vector quarkonium (EEQ), $e p\to epV$, is per se an interesting topic in studies of quarkonium production mechanism, QCD description of diffractive interaction and nucleon structure. We investigate this process in the framework of nonrelativistic QCD and QCD collinear factorization at the next-to-leading order QCD accuracy. The perturbative convergence behavior is discussed in a large range of photon virtuality $Q^2$. The $J/\psi$ large-$Q^2$ electroproduction data at HERA can be well explained, and the $\Upsilon$ differential production rate is predicted. The uncertainties in theoretical predictions with radiative corrections are greatly reduced. Notice the EEQ process is extremely sensitive to the gluon distribution in nucleon, the generalized parton distribution, our results will constraint the gluon density with high precision while confronting to the future experimental data. For the sake of comparing convenience, the analytic expressions are provided.


I. INTRODUCTION
The processes of exclusive electro-and photo-production of quarkonium γ ( * ) p → V p with V = J/ψ or Υ, are particular interesting and important but yet not well explored, especially the former. The photon here can be highly virtual for electroproduction or real for photoproduction. These processes provide unique opportunities in studying the quarkonium production mechanism, and perturbative QCD calculation reliability.
They are experimentally adjustable in physical parameters, say for example the virtuality of the photon. Furthermore, they are gluon rich, hence extremely sensitive to the gluon distribution in the nucleon, particular to the off-diagonal effects. Stable theoretical calculations are necessary therefore to reduce the uncertainty of gluon density distribution in the still vague domain of small Bjorken variable x B .
In experiment, J/ψ production processes have been extensively studied at HERA, whereas for Υ production the data are limited to the photoproduction case. For a review, see for example [1]. For future, some projects on deep inelastic experiment are in progress or proposed, like ENC at FAIR [2], eRHIC at BNL [3], LHeC at CERN [4] and EIC in China [5], where the EEQ process will be further explored attentively. Theoretically, two main frameworks are employed in the evaluation, that is the QCD collinear factorization [6] and BFKL k T factorization [7,8]. Although the BFKL approach has a solid perturbative QCD foundation, can sum up large logarithms of energy ln(1/x) and implies the k T information of gluon in the nucleon in the description of hard diffractive processes [9][10][11][12], it is impaired by the absence of full next-to-leading order (NLO) calculation, which however tends to be tough and may yield enormous corrections. Higher order calculation in collinear factorization hopefully can explain the existing Υ exclusive photoproduction data, but is fraught with difficulties for the J/ψ case [13][14][15].
Whereas the exclusive quarkonium electroproduction processes, which in some sense are even more important in physics due to the adjustable virtuality of the intermediate photon, have not been explored properly. In this paper, we calculate the NLO QCD corrections to exclusive quarkonium electroproduction processes, and investigate their implications to the parton distribution in the nucleon.
According to QCD factorization, the EEQ processes may be allocated into three domains, where the purtabative calculable sector and nonperturbative part belonging to different energy scales are properly separated. Namely, 1) The hard partonic process γ * g(q) → QQg(q). Due to the hard scales provided by the heavy quark mass or by the photon virtuality Q 2 , this part can be described by perturbative QCD (pQCD).
2) The transition from the QQ pair to the physical quarkonium state. Herein, the transition probability can be encoded into the non-relativistic QCD (NRQCD) [16] matrix element O V .
3) The parton distribution within the nucleon. This effect may be separated from the hard process via k T factorization or collinear factorization mechanism.
It is worth mentioning that the parton distribution in the nucleon here refers to the so-called generalized parton distribution (GPD), as in the case of deeply virtual Compton scattering (DVCS) [17][18][19]. The GPD extends the usual forward PDF to the non-forward situation, and encode more richer information about the nucleon, like the nucleon spin [18]. The study of GPD is nowadays a very dynamical field, for reviews see for instance references [20][21][22].

II. KINEMATICS AND FACTORIZATION
The kinematics of exclusive quarkonium production is schematically shown in Figure 1. The momenta of incident photon and proton, outgoing quarkonium and proton, are denoted by q, p, q ′ and p ′ respectively. In the calculation, notationsq = (q + q ′ )/2, p = (p + p ′ )/2 and ∆ = p ′ − p are also employed, and the following Lorentz invariants are then appeared in the analytic expression, where M and m N are the mass of quarkonium and proton respectively. Noticing that working in light-cone coordinate is convenient, we choose a frame whereq andp to be collinear. By introducing two light like vectors n + = (1, 0, 0, 1)/ √ 2 and n − = (1, 0, 0, −1)/ √ 2, their momenta can be expanded as and hence Here, The skewedness parameter η here plays a similar role as the Bjorken variable x B in deep-inelastic scattering. Following common usage, we use the term "small x B " instead of "small η" in this paper.
According to the NRQCD and collinear factorization, the EEQ process amplitude can be expressed as where m is the mass of heavy quark which is equal to M/2, e q is heavy quark electric charge. The superscript λ (λ ′ ) denotes the helicity of incoming photon (outgoing quarkonium). T λλ ′ p (T λλ ′ p ) and F p (F p ) represent the hard partonic amplitude and the matrix element of light-cone parton operators respectively. Their dependence on the factorization scale µ F is suppressed for brevity. Following the convention and definition of Ref. [20], F p andF p may be expressed as where H p (x, η, t), E p (x, η, t),H p (x, η, t) andẼ p (x, η, t) are twist-2 GPDs with arguments of x and η. The parton momentum fractions are represented by the combination of x and η as shown in Fig.1. The effect of higher-twist GPDs are generally complicated and small in comparison with the NLO perturbative corrections [23], which hence will not be taken account in this work.
Note, the symmetric properties of GPDs may simplify the calculation. The gluon distributions satisfy and similarly are the E g andẼ g . For quark distributions, different combinations are considered usually in the calculation. That is referring to the "singlet" combination, and for the "nonsinglet" case. As well, similar combinations exist for E q andẼ q . Since photon and vector quarkonium have the same C parities, only the C even (singlet) component of quarks in GPD contributes to (5).

III. ANALYTICAL RESULTS
The typical leading order (LO) and NLO Feynman diagrams are shown in Fig.2.
At LO, only gluon involved process contributes, whereas the light quark induced process may appear at NLO. The effects of intrinsic heavy quark inside the nucleon are reasonably small and will be neglected.
The convolution integration in (5) stretches from |x| > η (DGLAP) region to the |x| < η (ERBL) region. In the |x| < η region, the amplitude does not contain imaginary part, and the iε prescription in propagators can be dropped. The analytic continuation is performed by restoring the iε via The manipulation of calculation may be simplified greatly by making a coordinate transformation, namely We then have Here, µ R is the renormalization scale. The helicity amplitudes obey A ++ in accordance with the requirement of helicity conservation. There exist also crossing symmetries under the variable exchange The LO result is pretty simple, In the computation of NLO corrections, the ultraviolet (UV) and infrared (IR) singularities are regulated in dimensional regularization with D = 4 − 2ǫ prescription adopted. The following singular terms arising from one-loop diagrams contain both UV and IR singularities: The UV singularities are removed by renormalization. For the renormalization of heavy quark field, heavy quark mass and gluon field, we take the on-shell (OS) scheme; for the renormalization of coupling constant, the modified minimal-subtraction (MS) scheme is used. The singular terms arising from counter terms are After the cancellation between (19) and (20), the remaining singularities will be absorbed into the parton distribution functions, achieved by introducing the scale de- where V pp ′ denotes the GPD evolution kernel.
Finally, the NLO result is finite and can be written in a general form where the plus sign corresponds to A g andÃ q , while the minus to A q andÃ g . The coefficients c 0 and c i in (22) are a bit lengthy and are given in Appendix. The f i are either logrithm or polylogarithm functions yielded from the Feynman integration: with Note, by taking the Q → 0 limit in (22), we can readily reproduce the amplitude of quarkonium photoproduction [13].
At the high energy limit, the leading contribution to the NLO correction comes from the region η ≪ |x| ≪ 1, and the amplitude can be simplified to This expression suggests a suitable value of factorization scale, µ F = m 2 + Q 2 4 . By taking this value, the NLO correction in (25) vanishes.

IV. NUMERICAL ANALYSIS
For the full numerical evaluation, the knowledge of GPD over the full range of x (or at least x > η) is required. Unfortunately, the available models for GPD are fraught with uncertainties, especially in the ERBL region. In order to minimize the uncertainties, we take the Forward Model (FM) in the DGLAP region, based on which the imaginary sector of amplitude can be calculated. The FM tells H g (x, η, µ 0 ) = xg(x, µ 0 ), for x > η , with µ 0 = 1 GeV as initial scale, and MSTW08 [24] as input PDF. The GPDs at the energy scale of our concern are obtained through skewed evolution equation, where the NLO evolution kernels employed come from [25]. From the imaginary part of the amplitude, real part can be restored via the dispersion relation [26]: In numerical evaluation, those contributions from E p ,Ẽ p andH p are neglected, due to the following reasons. First, these terms are kinematic or helicity suppressed.
Secondly, a rough numerical estimation tells that those terms contribute less than 1% of the total. Moreover, there still lack convinced models to evaluating these GPDs.
Other parameters taken in the calculation go as follows: • Λ 3 QCD = 332 MeV , Λ 4 QCD = 292 MeV , Λ 5 QCD = 210 MeV ; • |R J/ψ (0)| 2 = 0.903 GeV 3 , |R Υ (0)| 2 = 7.76 GeV 3 ; • 1.4 GeV≤ m c ≤ 1.6 GeV , 4.8 GeV≤ m b ≤ 5.0 GeV ; The renormalization scale in our evaluation is set to be equal to the factorization scale, the so called BLM scale [27], by which contributions from terms proportional to Q 2 ≈ 2 GeV 2 , the absolute value of LO+NLO amplitudes are very small, which will lead to small predicted cross section. For Υ production, due to the large quark mass, the convergence is good even for the photoproduction case, as demonstrated in Figure   3(b).
Confronting to the HERA experimental condition, we calculate the differential cross section dσ/dt at |t| = |t| min , i.e. ∆ ⊥ = 0. Taking t dependence measurement from H1 experiment as input [28], we obtain the total cross section with different W and Q.
The numerical results of J/ψ electroproduction are shown in Figure 4 and 5. Although the LO prediction may somehow cover the experimental data, the uncertainties are dubiously large, which mainly comes from the strong dependence of GPD on the factorization scale µ F , especially at small η. With NLO 1 corrections, at Q 2 > 10 GeV 2 , the uncertainty is largely suppressed as expected, and the prediction agrees with the experimental measurement well. At low Q 2 , there are not many experimental data and the theoretical predicability is impaired by the large uncertainty, which remains even with NLO corrections. Agreeing with above amplitude convergence analysis, numeri- cal results for cross section in the region of Q 2 < 5 GeV 2 are unreliable, which casts a shadow on the perturbative QCD evaluation of J/ψ photoproduction. Note, the dips on blue lines near Q 2 ≈ 2 GeV 2 in Figure 4 are understandable because the amplitude with NLO corrections drops off at this point, as shown in Figure 3.
In comparison to the J/ψ production, the theoretical evaluation of Υ production will definitely be more reliable because of the higher energy scale herein. However, unfortunately, to date there still has been no experimental result on the Υ leptoproduction yet. But hopefully it would be investigated on future lepton-nucleon colliders.
For this aim, we calculate the Q 2 dependence of dσ/dt at |t| = |t| min for Υ production in HERA experimental condition for illustration, as shown in Figure 6.

V. SUMMARY AND CONCLUSIONS
We calculated analytically the exclusive electroproduction of quarkonium in the NRQCD framework and collinear factorization scheme up to the NLO QCD accuracy.
In order to compare with experimental measurements, numerical evaluation about the cross section of exclusive J/ψ electroproduction at different Q 2 and W are performed.
We estimated the theoretical uncertainties by varying the magnitudes of heavy quark mass and factorization scale.
At large Q 2 , say Q 2 > 10 GeV 2 , the NLO corrections may greatly reduce the theoretical uncertainty, and hence enable the predictions more reliable. We find a good agreement with the H1 [28] and ZEUS [29] data. At small Q 2 , say Q 2 < 5 GeV 2 , the pQCD analysis on exclusive J/ψ production tends to be dubious. To make a prediction for future experiment on ep collision, we schematically calculate the exclusive Υ electroproduction in HERA condition. In this case the quark mass is much heavier and guarantees the legitimacy of pQCD use even in the photoproduction.
To understand more about the parton distribution in nucleon, the GPD here, is one of the motivations of our study. With reliable theoretical calculation, people expect to have some feedbacks on GPD by confronting to the experimental measurement. About the GPD, in the calculation we made use of the Forward Model (26) together with the NLO GPD evolution equation to evaluate GPDs at DGLAP region. The input PDFs and the initial scale were set to be MSTW08 [24] and µ 0 = 1 GeV. As a trial, another set of input PDFs, the CT14 [30] and initial scale µ 0 = 1.3 GeV were applied. We found at reasonably large η and µ F regions, say µ F > 2.4 GeV and η > 0.001, the results of cross section in two different choices are less than 15%, which does not influence our conclusions for J/ψ electroproduction at large-Q 2 and the Υ electroproduction. The main reason of this is that the discrepancy between MSTW08 and CT14 diminishes with the increase of energy scale in off small x B region. We therefore conclude that the Forward Model is simpler, parameter-free and adequate to explain the data.
Last, as an one step further investigation, we also calculated the concerned processes by using the GPD model proposed by Freund, McDermott, and Strikeman (FMS) [31] and the Double Distribution Model [32][33][34]. The former yields a similar numerical result as Forward Model does at large-Q 2 region within 10% of discrepancy, while the latter will produce a very large result stemming from the enhanced quark contribution, which needs to be tamed somehow. Note, the FMS model agrees with the forward model (26) in the DGLAP region, while in the ERBL region it is simply an ansatz based on the polynomiality for the lowest Mellin moments. In this model, the real part of the amplitude can be directly calculated without employing the dispersion relation.  (y 1 −5y 2 )(y 2 −5y 1 ) (y 1 −y 2 ) 2 (y 1 +y 2 ) 2 (  (y 1 −5y 2 )(y 2 −5y 1 ) (y 1 −y 2 ) 2 (y 1 +y 2 ) 2 (−y 2 1 + 3y 1 y 2 + y 1 y 2 − 4y 1 + 1) .