NLO QCD corrections to exclusive electroproduction of quarkonium

TL;DR

This paper calculates next-to-leading order QCD corrections to exclusive electroproduction of quarkonium, improving theoretical predictions and providing insights into gluon distributions in nucleons, with implications for future experiments.

Contribution

It presents the first NLO QCD calculations for exclusive quarkonium electroproduction within the nonrelativistic QCD framework, enhancing prediction accuracy and constraining gluon distributions.

Findings

HERA $J/\psi$ data well explained at high $Q^2$

Predicted $\Upsilon$ production rates

Reduced uncertainties in theoretical predictions

Abstract

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.

Equations

$$Q^2=$$

$$s=(q+p{)}^2=W^2,t={\Delta }^2$$

$${\bar{p}}^{\mu }=n_{+}^{\mu }+\frac{m_N^2-t/4}{2}{n}_{-}^{\mu },{\bar{q}}^{\mu }=-\xi n_{+}^{\mu }+\frac{Q^2-M^2+t/2}{4\xi }{n}_{-}^{\mu },$$

$${\Delta }^{\mu }=-2\eta n_{+}^{\mu }+\eta (m_N^2-t/4)n_{-}^{\mu }+{\Delta }_{\perp }^{\mu }.$$

$$\xi =-\frac{\bar{q}\cdot n_{-}}{\bar{p}\cdot n_{-}}\approx \frac{Q^2-M^2}{2s+Q^2-M^2},\eta =-\frac{1}{2}\frac{\Delta \cdot n_{-}}{\bar{p}\cdot n_{-}}\approx \frac{Q^2+M^2}{2s+Q^2-M^2}.$$

$${\mathcal{M}}^{\lambda {\lambda }^{\prime }}=$$

$$\times \sum _{p=g,q}{\int }_{-1}^{1}dx(T_p^{\lambda {\lambda }^{\prime }}(x,\eta ,\xi )F^p(x,\eta ,t)+{\tilde{T}}_p^{\lambda {\lambda }^{\prime }}(x,\eta ,\xi ){\tilde{F}}^p(x,\eta ,t)),$$

$$F^p=\frac{1}{2(\bar{p}\cdot n_{-})}\left[H^p(x,\eta ,t)\bar{u}(p^{\prime })n{/}_{-}u(p)+E^p(x,\eta ,t)\bar{u}(p^{\prime })\frac{i{\sigma }^{n_{-}\Delta }}{2m_N}u(p)\right],$$

$${\tilde{F}}^p=\frac{1}{2(\bar{p}\cdot n_{-})}\left[{\tilde{H}}^p(x,\eta ,t)\bar{u}(p^{\prime })n{/}_{-}{\gamma }_{5}u(p)+{\tilde{E}}^p(x,\eta ,t)\bar{u}(p^{\prime })\frac{{\gamma }_5(\Delta \cdot n_{-})}{2m_N}u(p)\right],$$

$$H^g(x,\eta ,t)=H^g(-x,\eta ,t),{\tilde{H}}^g(x,\eta ,t)=-{\tilde{H}}^g(-x,\eta ,t),$$

$$H^{q(+)}(x,\eta ,t)=H^q(x,\eta ,t)-H^q(-x,\eta ,t),$$

$$\displaystyle\tilde{H}^{q(+)}(x,\eta,t)=\tilde{H}^{q}(x,\eta,t)+\tilde{H}^{q}(-x,\eta,t)\$$

$$H^{q(-)}(x,\eta ,t)=H^q(x,\eta ,t)+H^q(-x,\eta ,t),$$

$${\tilde{H}}^{q(-)}(x,\eta ,t)={\tilde{H}}^q(x,\eta ,t)-{\tilde{H}}^q(-x,\eta ,t)$$

$$x\to x+i\epsilon ,\text{for}x>\eta ,$$

$$x\to x-i\epsilon ,\text{for}x<-\eta .$$

$$y_1=\frac{\eta -x}{\eta -\xi },y_2=\frac{\eta +x}{\eta -\xi }.$$

$$T_g^{\lambda {\lambda }^{\prime }}(x,\eta ,\xi )=\frac{1}{x^2-{\eta }^2}\left(A_g^{(0)\lambda {\lambda }^{\prime }}(y_1,y_2)+\frac{{\alpha }_s({\mu }_R)}{\pi }{A}_g^{(1)\lambda {\lambda }^{\prime }}(y_1,y_2)\right),$$

$${\tilde{T}}_g^{\lambda {\lambda }^{\prime }}(x,\eta ,\xi )=\frac{1}{x^2-{\eta }^2}\frac{{\alpha }_s({\mu }_R)}{\pi }{\tilde{A}}_g^{(1)\lambda {\lambda }^{\prime }}(y_1,y_2),$$

$$T_q^{\lambda {\lambda }^{\prime }}(x,\eta ,\xi )=\frac{1}{\eta -\xi }\frac{{\alpha }_s({\mu }_R)}{\pi }{A}_q^{(1)\lambda {\lambda }^{\prime }}(y_1,y_2),$$

$${\tilde{T}}_q^{\lambda {\lambda }^{\prime }}(x,\eta ,\xi )=\frac{1}{\eta -\xi }\frac{{\alpha }_s({\mu }_R)}{\pi }{\tilde{A}}_q^{(1)\lambda {\lambda }^{\prime }}(y_1,y_2).$$

$$A_g^{(0)++}=A_g^{(0)--}=\frac{-1}{\sqrt{y_1+y_2-1}}{A}_g^{(0)00}=-\frac{1}{3}\frac{1}{y_1+y_2}.$$

$$A_{g,\text{loop-pole}}^{(1),++}$$

$$\displaystyle-\left(y_{1}^{2}+\frac{y_{1}y_{2}}{2}+\frac{y_{2}^{2}}{2}\right)\text{ln}(y_{1}+y_{2})+\frac{3C_{F}}{8C_{A}}(y_{1}^{2}-3y_{1}y_{2}+2y_{2})\bigg{)}+\{y_{1}\leftrightarrow y_{2}\}\ ,$$

$$A_{g,\text{loop-pole}}^{(1),00}$$

$$\displaystyle-\frac{3C_{F}}{4C_{A}}(2y_{1}y_{2}-y_{2})\bigg{)}+\{y_{1}\leftrightarrow y_{2}\}\ ,$$

$$A_{q,\text{loop-pole}}^{(1),++}$$

$$\displaystyle=-\frac{1}{\epsilon}\frac{C_{A}C_{F}(y_{1}-y_{2})}{9y_{2}(y_{1}+y_{2})^{2}}\bigg{(}\text{ln}y_{1}-\frac{y_{1}+y_{2}}{2y_{1}}\text{ln}(y_{1}+y_{2})\bigg{)}-\{y_{1}\leftrightarrow y_{2}\}\ ,$$

$${\tilde{A}}_{g,\text{loop-pole}}^{(1),++}$$

$${\tilde{A}}_{g,\text{loop-pole}}^{(1),00}$$

$$A_{g,\text{ct-pole}}^{(1),++}$$

$$A_{g,\text{ct-pole}}^{(1),00}$$

$${\tilde{A}}_{g,\text{ct-pole}}^{(1),++}$$

$$A_{q,\text{ct-pole}}^{(1),++}$$

$$F^p(x,\eta ,{\mu }_F)=F^p(x,\eta )-\frac{1}{ϵ}\left[\frac{{\alpha }_s}{2\pi }\frac{\Gamma (1-ϵ)}{\Gamma (1-2ϵ)}{\left(\frac{4\pi {\mu }_R^2}{{\mu }_F^2}\right)}^{ϵ}\right]\sum _{p^{\prime }}{\int }_{-1}^{1}dy{V}_{p{p}^{\prime }}(x,y,\eta )F^{p^{\prime }}(y,\eta ),$$

$$A^{(1)\lambda\lambda\prime}=\bigg{(}c_{0}^{\lambda\lambda\prime}+\sum_{i=1}^{11}c_{i}^{\lambda\lambda\prime}f_{i}\bigg{)}\pm\{y_{1}\leftrightarrow y_{2}\}\ ,$$

$$f_1=\sqrt{\frac{y_1-1}{y_1}}\text{ln}(\sqrt{y_1-1}+\sqrt{y_1}),f_2=f_1^2,$$

$$f_3=\sqrt{\frac{y_1+y_2}{y_1+y_2-1}}\text{ln}(\sqrt{y_1+y_2-1}+\sqrt{y_1+y_2}),f_4=f_3^2,$$

$$f_5=\text{ln}(\frac{y_1}{y_1+y_2}),f_6=\text{ln}2+\text{ln}(y_1+y_2),f_7={\text{ln}}^2(y_1+y_2),$$

$$f_8=\text{ln}(\frac{4m^2}{{\mu }_F^2})\text{ln}(\frac{y_1}{y_1+y_2})+{\text{ln}}^2{y}_1,f_9={\text{Li}}_2(1-2y_1),f_{10}={\text{Li}}_2(1-2y_1-2y_2),$$

$$f_{11}=C_0(4(1-y_1-y_2),1-2y_1,1-2y_2,1,1,0),$$

$$C_0(4(1-y_1-y_2),1-2y_1,1-2y_2,1,1,0)=$$

Full text

PACS numbers: 13.60.-r, 14.20.Dh, 14.40.Pq

I Introduction

The processes of exclusive electro- and photo-production of quarkonium $\gamma^{(*)}p\to Vp$ with $V=J/\psi$ or $\Upsilon$, 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/\psi$ production processes have been extensively studied at HERA, whereas for $\Upsilon$ production the data are limited to the photoproduction case. For a review, see for example exreview . For future, some projects on deep inelastic experiment are in progress or proposed, like ENC at FAIR ENC , eRHIC at BNL eRHIC , LHeC at CERN LHeC and EIC in China EicC , where the EEQ process will be further explored attentively. Theoretically, two main frameworks are employed in the evaluation, that is the QCD collinear factorization cfs1 and BFKL $k_{T}$ factorization cch1 ; ce1 . 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 ryskin ; ktfrthr1 ; ktfrthr2 ; ktfrthr3 , 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 $\Upsilon$ exclusive photoproduction data, but is fraught with difficulties for the $J/\psi$ case photoNLO1 ; photoNLO2 ; photoNLO3 . 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 $\gamma^{*}g(q)\to Q\bar{Q}g(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 $Q\bar{Q}$ pair to the physical quarkonium state. Herein, the transition probability can be encoded into the non-relativistic QCD (NRQCD) matrix element $\langle O^{V}\rangle$ NRQCD .

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) DVCS1 ; Jispin ; DVCS3 . The GPD extends the usual forward PDF to the non-forward situation, and encode more richer information about the nucleon, like the nucleon spin Jispin . The study of GPD is nowadays a very dynamical field, for reviews see for instance references review1 ; review2 ; review3 .

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^{\prime}$ and $p^{\prime}$ respectively. In the calculation, notations $\bar{q}=(q+q^{\prime})/2$, $\bar{p}=(p+p^{\prime})/2$ and $\Delta=p^{\prime}-p$ are also employed, and the following Lorentz invariants

[TABLE]

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 where $\bar{q}$ and $\bar{p}$ to be collinear. By introducing two light like vectors $n_{+}=(1,0,0,1)\Lambda$ and $n_{-}=(1,0,0,-1)/(2\Lambda)$, where the arbitrary value parameter $\Lambda$ has the dimension of mass, their momenta can be expanded as

[TABLE]

and hence

[TABLE]

Here,

[TABLE]

The skewedness parameter $\eta$ 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 $\eta$” in this paper.

According to the NRQCD and collinear factorization, the EEQ process amplitude can be expressed as

[TABLE]

where $m$ is the mass of heavy quark equal to $M/2$ in the leading order of relativistic expansion, $e_{q}$ is heavy quark electric charge. The superscript $\lambda$ ($\lambda^{\prime}$) denotes the helicity of incoming photon (outgoing quarkonium). $T_{p}^{\lambda\lambda^{\prime}}$ ($\tilde{T}_{p}^{\lambda\lambda^{\prime}}$) and $F^{p}$ ($\tilde{F}^{p}$) represent the hard partonic amplitude and the matrix element of light-cone parton operators respectively. Their dependence on the factorization scale $\mu_{F}$ is suppressed for brevity. Following the convention and definition of Ref.review1 , $F^{p}$ and $\tilde{F}^{p}$ may be expressed as

[TABLE]

where $H^{p}(x,\eta,t)$, $E^{p}(x,\eta,t)$, $\tilde{H}^{p}(x,\eta,t)$ and $\tilde{E}^{p}(x,\eta,t)$ are twist-2 GPDs with arguments $x$ and $\eta$. The parton momentum fractions are represented by the combination of $x$ and $\eta$ as shown in Fig.1. The effect of higher-twist GPDs are generally complicated and small in comparison with the NLO perturbative corrections MJP , which hence will not be taken account in this work.

Note, the symmetric properties of GPDs may simplify the calculation. The gluon distributions satisfy

[TABLE]

and similarly are the $E^{g}$ and $\tilde{E}^{g}$. For quark distributions, different combinations are considered usually in the calculation. That is

[TABLE]

referring to the “singlet” combination, and

[TABLE]

for the “nonsinglet” case. As well, similar combinations exist for $E^{q}$ and $\tilde{E}^{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|>\eta$ (DGLAP) region to the $|x|<\eta$ (ERBL) region. In $|x|<\eta$ region, the amplitude does not contain imaginary part, and the $i\varepsilon$ prescription in propagators can be dropped. The analytic continuation is performed by restoring the $i\varepsilon$ via

[TABLE]

The manipulation of calculation may be simplified greatly by making a coordinate transformation, namely

[TABLE]

We then have

[TABLE]

Here, $\mu_{R}$ is the renormalization scale. The helicity amplitudes obey $A_{p}^{++}=A_{p}^{--}$, $\tilde{A}_{p}^{++}=-\tilde{A}_{p}^{--}$ and $A^{0+}=A^{0-}=A^{+0}=A^{-0}=\tilde{A}^{0+}=\tilde{A}^{0-}=\tilde{A}^{+0}=\tilde{A}^{-0}=0$, which in accordance with the requirement of helicity conservation. There exist also crossing symmetries under the variable exchange $y_{1}\leftrightarrow y_{2}$ (i.e. $x\leftrightarrow-x$), i.e., $A_{g}(y_{1},y_{2})=A_{g}(y_{2},y_{1})$, $\tilde{A}_{g}(y_{1},y_{2})=-\tilde{A}_{g}(y_{2},y_{1})$, $A_{q}(y_{1},y_{2})=-A_{q}(y_{2},y_{1})$, $\tilde{A}_{q}(y_{1},y_{2})=\tilde{A}_{q}(y_{2},y_{1})$.

The LO result is pretty simple,

[TABLE]

In the computation of NLO corrections, the ultraviolet (UV) and infrared (IR) singularities are regulated in dimensional regularization with $D=4-2\epsilon$ prescription adopted. The following singular terms arising from one-loop diagrams contain both UV and IR singularities:

[TABLE]

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 ($\overline{\text{MS}}$) scheme is used. The singular terms arising from counter terms are

[TABLE]

After the cancellation between (19) and (20), the remaining singularities will be absorbed into the parton distribution functions, achieved by introducing the scale dependent GPD, i.e.

[TABLE]

where $V_{pp^{\prime}}$ denotes the GPD evolution kernel.

Finally, the NLO result is finite and can be written in a general form

[TABLE]

where the plus sign corresponds to $A_{g}$ and $\tilde{A}_{q}$, while the minus to $A_{q}$ and $\tilde{A}_{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:

[TABLE]

with

[TABLE]

Note, by taking the $Q\to 0$ limit in (22), we can readily reproduce the amplitude of quarkonium photoproduction photoNLO1 .

At the high energy limit, the leading contribution to the NLO correction comes from the region $\eta\ll|x|\ll 1$, and the amplitude can be simplified to

[TABLE]

This expression suggests a suitable value of factorization scale, $\mu_{F}=\sqrt{m^{2}+\tfrac{Q^{2}}{4}}$. By taking this value, the NLO correction in (25) vanishes.

IV Numerical analysis

For full numerical evaluation, the knowledge of GPD over the full range of $x$ (or at least $x>\eta$) 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

[TABLE]

with $\mu_{0}=1$ GeV as initial scale, and MSTW08 mstw 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 evo . From the imaginary part of the amplitude, real part can be restored via the derivative dispersion relation (DDR) dispersion :

[TABLE]

In numerical evaluation, those contributions from $E^{p}$, $\tilde{E}^{p}$ and $\tilde{H}^{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:

• •

$\Lambda_{\text{QCD}}^{3}=332$ MeV , $\Lambda_{\text{QCD}}^{4}=292$ MeV , $\Lambda_{\text{QCD}}^{5}=210$ MeV ;

• •

$|R_{J/\psi}(0)|^{2}=0.903$ GeV3 , $|R_{\Upsilon}(0)|^{2}=7.76$ GeV3 ;

• •

$1.4$ GeV$\leq m_{c}\leq 1.6$ GeV , $4.8$ GeV$\leq m_{b}\leq 5.0$ GeV ;

• •

\text{max}\big{[}\frac{1}{2}\sqrt{m^{2}+\frac{Q^{2}}{4}},\mu_{0}\big{]}\leq\mu_{F}\leq 2\sqrt{m^{2}+\frac{Q^{2}}{4}} ;

• •

$\mu_{R}=\mu_{F}$ .

The renormalization scale in our evaluation is set to be equal to the factorization scale, by which contributions from terms proportional to $\beta_{0}\ln\frac{\mu_{R}}{\mu_{F}}$ are eliminated (in $c_{0}$s of Eqs. (28) and (29) in the Appendix). The theoretical uncertainties are estimated by varying the values of quark mass and factorization scale.

The imaginary parts of the amplitudes in $J/\psi$ electroproduction as function of $Q^{2}$ are shown in Figure 3(a), where the perturbative convergence exhibits. As $Q^{2}$ increases, the perturbative convergence is improved due to the increase of energy scale and the departure of small $x_{B}$ region. At large $Q^{2}$, say $Q^{2}>10$ GeV2, the convergence manifests itself well. At low $Q^{2}$, the convergence is poor, especially when $Q^{2}<2$ GeV2, where the full LO+NLO amplitude has opposite sign to the LO one. In the region $Q^{2}\approx 2$ GeV2, the absolute value of LO+NLO amplitudes are very small, and hence the cross section as well. For $\Upsilon$ production, due to the large quark mass, the convergence works well even for the photoproduction, as demonstrated in Figure 3(b).

Predictions of ${\rm Im}\mathcal{M}/W^{2}$ for $J/\psi$ electroproduction at $Q^{2}=22.4$ GeV2 as function of $W$ are shown in Figure 4 (a), based on which we can estimate the error of DDR (27). The DDR was derived from the standard integral dispersion relation (IDR) IDR1 ; IDR2 by taking the $s\to\infty$ limit and neglecting the subtraction constant term. An extended derivative dispersion relation (EDDR) given in Ref. EDDR is equivalent to IDR. In practice, when $s$ is big enough but finite the applicability and accuracy of (27) depend upon the behavior of ${\rm Im}\mathcal{M}/W^{2}$. By fitting ${\rm Im}\mathcal{M}/W^{2}$ to the form of $a_{1}W^{b_{1}}+a_{2}W^{b_{2}}$ and calculating ${\rm Re}\mathcal{M}$ via the DDR and EDDR separately, we find the discrepancy in amplitudes of these two methods is less than $1\%$. Besides, in fact the variable $s$ in (27) should be replaced by $\nu=s+\frac{Q^{2}}{2}-\frac{M^{2}}{2}$, and the corresponding corrections are then suppressed by a factor of $\tfrac{Q^{2}-M^{2}}{2s}$, which is less than $10^{-2}$ in our case. The ratios of real parts, restored via (27), to imaginary parts as function of $W$ are shown in Figure 4 (b).

Confronting to the HERA experimental condition, we calculate the differential cross section $d\sigma/dt$ at $|t|=|t|_{\text{min}}$, i.e. $\Delta_{\perp}=0$. Taking $t$ dependence measurement from H1 experiment as input H1 , we obtain the total cross section with different $W$ and $Q$. The numerical results of $J/\psi$ electroproduction are shown in Figure 5 and 6. The $\psi(2S)$ feed down contribution is not taken into account in our calculation, since in experiment H1 ; ZEUS this effect was subtracted in the data. 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 $\mu_{F}$, especially at small $\eta$. With NLO333Strictly speaking, matrix element squared $|\mathcal{M}^{(0)}+\mathcal{M}^{(1)}|^{2}$ includes some of the NNLO contributions. corrections, at $Q^{2}>10$ GeV2, the uncertainty is greatly 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, numerical results for cross section in the region of $Q^{2}<5$ GeV2 are unreliable, which casts a shadow on the perturbative QCD evaluation of $J/\psi$ photoproduction. Note, the dips on blue lines near $Q^{2}\approx 2$ GeV2 in Figure 5 are understandable because the amplitude with NLO corrections drops off at this point, as shown in Figure 3.

In comparison to the $J/\psi$ production, the theoretical evaluation of $\Upsilon$ 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 $\Upsilon$ leptoproduction yet. But hopefully it would be investigated on future lepton-nucleon colliders. For this aim, we calculate the $Q^{2}$ dependence of $d\sigma/dt$ at $|t|=|t|_{\text{min}}$ for $\Upsilon$ production in HERA experimental condition for illustration, as shown in Figure 7.

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/\psi$ electroproduction at different $Q^{2}$ and $W$ were 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$ GeV2, the NLO corrections may greatly reduce the theoretical uncertainty, and hence enable the predictions more reliable. We found a good agreement with the H1 H1 and ZEUS ZEUS data. At small $Q^{2}$, say $Q^{2}<5$ GeV2, the pQCD analysis on exclusive $J/\psi$ production tends to be dubious. To make a prediction for future experiment on $ep$ collision, we schematically calculated the exclusive $\Upsilon$ 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 GPD at DGLAP region. The input PDF and the initial scale were set to be MSTW08 mstw at $\mu_{0}=1$ GeV. As a trial, another set of input PDF, the CT14 ct14 and initial scale $\mu_{0}=1.3$ GeV were applied. We found at reasonably large $\eta$ and $\mu_{F}$ region, say $\mu_{F}>2.4$ GeV and $\eta>0.001$, the difference of the cross sections with those two sets of input PDFs and initial scales is less than $15\%$, which does not influence our conclusions on $J/\psi$ electroproduction at large-$Q^{2}$ and $\Upsilon$ electroproduction. The main reason of this is because that the discrepancy between MSTW08 and CT14 diminishes with the increase of energy scale in 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) gpdmodel , the Shuvaev transform shuvaev1 ; shuvaev2 approach and the Double Distribution (DD) model dd1 ; dd2 ; dd3 , and found:

1. 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 of the lowest Mellin moments. In this model, the real part of the amplitude can be directly calculated without employing the dispersion relation. For large-$Q^{2}$ $J/\psi$ electroproduction and $\Upsilon$ electroproduction, the FMS model may yield a similar numerical result as the Forward Model does with about $10\%$ of discrepancy.

2. The Shuvaev transform can correlate the GPD to an auxiliary function, the usual parton distribution function in our evaluation. By performing the transform, the skewed effect from evolution (at least LO evolution444Hopefully also the higher order evolution in the conformal scheme ($\overline{\rm CS}$).) would be incorporated as well. The Shuvaev transform can also be viewed approximately as the LO GPD evolution of the $\eta$-dependent initial distribution. However, with the progress of evolution, the discrepancy between the $\eta$-dependent initial distribution and the initial distribution in Forward model shrinks, while the discrepancy may be enlarged by different evolution equations. Numerically, except for the $x\to 1$ region where the parton distributions tend to be diminished, the discrepancy between GPD from Shuvaev-transform approach and what from the Forward model is less than $15\%$ at moderate scale ($\mu>2$ GeV), and the difference in amplitude is about $8\%$ and $30\%$ at LO and NLO respectively.

3. The DD model is also a widely used GPD model. It provides an natural way to implement the polynomiality, nevertheless is unstable in evolution. It may yield a huge result stemming from the enhanced quark contribution, which should be tamed somehow.

Acknowledgements

This work was supported in part by the Ministry of Science and Technology of the Peoples’ Republic of China(2015CB856703); by the Strategic Priority Research Program of the Chinese Academy of Sciences, Grant No.XDB23030100; and by the National Natural Science Foundation of China(NSFC) under the Grants 11375200 and 11635009.

Appendix

The nonzero coefficients $c_{i}$ in (22).

For $A_{g}^{++}$($A_{g}^{--}$), $c^{++}_{i}=c^{--}_{i}=c_{i}$, with

[TABLE]

For $A_{g}^{00}$, $c^{00}_{i}=-\sqrt{y_{1}+y_{2}-1}c_{i}$, with

[TABLE]

For $A_{q}^{++}$($A_{q}^{--}$), $c^{++}_{i}=c^{--}_{i}=c_{i}$, with

[TABLE]

For $A_{q}^{00}$, $c^{00}_{i}=-\sqrt{y_{1}+y_{2}-1}c_{i}$, with

[TABLE]

For $\tilde{A}_{g}$, $\tilde{c}^{00}_{i}=0,\tilde{c}^{++}_{i}=-\tilde{c}^{--}_{i}=c_{i}$, with

[TABLE]

For $\tilde{A}_{q}$, $\tilde{c}^{00}_{i}=0,\tilde{c}^{++}_{i}=-\tilde{c}^{--}_{i}=c_{i}$, with

[TABLE]