QCD leading order study of the $J/\psi$ leptoproduction at HERA within the nonrelativistic QCD framework

TL;DR

This paper presents a leading-order QCD analysis of $J/ar{J}$ production in deep inelastic scattering at HERA within the nonrelativistic QCD framework, highlighting the importance of color-octet contributions for better experimental agreement.

Contribution

It provides renewed leading-order predictions for $J/\psi$ leptoproduction at HERA using nonrelativistic QCD, comparing different long-distance matrix element sets and emphasizing the role of color-octet states.

Findings

Color-singlet model predictions are below experimental data.

Color-octet contributions are comparable to singlet and improve data agreement.

The study highlights the limitations of previous collinear factorization approaches.

Abstract

As indicated in our previous paper (Phys.Rev.D96.034002), the existing literatures studying the $J/\psi$ production in deeply inelastic scattering (DIS) in collinear factorization are on the basis of a formalism that will lead to wrong results when the ranges of the transverse momentum or the rapidity of the $J/\psi$ in the laboratory frame do not cover all values possible for them. In this paper, we present the renewed results for the $J/\psi$ production in DIS at HERA within the nonrelativistic QCD framework at QCD leading order (LO). Three different sets of the long-distance matrix elements are employed for comparison. The predictions via the colour-singlet (CS) model at QCD LO are generally below the experimental data especially in the regions where perturbation theory are expected to work well, while the colour-octet contributions are of the same order of magnitude as the CS ones, however, in general make the agreement between theory and experiment better.

Tables

Table 1. Table 1: Comparison of our results with those of Reference  Kniehl and Zwirner ( 2002 ) .

Type	Cuts	Reference Kniehl and Zwirner (2002)	Our results I	Our results II
CS	$Q^2>4{\mathrm{GeV}}^2$	107pb	103.4pb	103.4pb
CS	$Q^2$, $p_t^2>4{\mathrm{GeV}}^2$	62pb	60.0pb	54.5pb
CS	$Q^2>4{\mathrm{GeV}}^2$, $p_t^{\star 2}>2{\mathrm{GeV}}^2$,	24pb	23.5pb	23.5pb
CO	$Q^2>4{\mathrm{GeV}}^2$, $p_t^{\star 2}>2{\mathrm{GeV}}^2$,	16pb	15pb	15pb

Table 2. Table 2: Three sets of the LDMEs used in our calculation.

	$⟨O^{J/\psi }{(}^3{S}_1^{\left[1\right]})⟩$	$⟨O^{J/\psi }{(}^1{S}_0^{\left[8\right]})⟩$	$⟨O^{J/\psi }{(}^3{S}_1^{\left[8\right]})⟩$	$⟨O^{J/\psi }{(}^3{P}_0^{\left[8\right]})⟩$
	${\mathrm{GeV}}^3$	$\times {10}^{-2}{\mathrm{GeV}}^3$	$\times {10}^{-2}{\mathrm{GeV}}^3$	$\times {10}^{-2}{\mathrm{GeV}}^5$
Reference Chao et al. (2012)	$1.16$	$8.9\pm 0.98$	$0.3\pm 0.12$	$1.26\pm 0.47$
Reference Butenschoen and Kniehl (2011c)	$1.32$	$3.04\pm 0.35$	$0.168\pm 0.046$	$-0.908\pm 0.161$
Reference Zhang et al. (2015)	$0.645\pm 0.405$	$0.785\pm 0.42$	$1.0\pm 0.3$	$3.8\pm 1.1$

Table 3. Table 3: χ 2 superscript 𝜒 2 \chi^{2} /d.o.f. for the LDMEs taken from References  Chao et al. ( 2012 ) ,  Butenschoen and Kniehl ( 2011c ) , and  Zhang et al. ( 2015 ) , respectively.

Condition	Reference Chao et al. (2012)	Reference Butenschoen and Kniehl (2011c)	Reference Zhang et al. (2015)
	11.54	2.07	7.03
	2.35	1.90	1.82

Equations

$$Q^2=-q^2\equiv -(k-k^{\prime }{)}^2,W^2=(P+q{)}^2,$$

$$x_B=\frac{Q^2}{2P\cdot q},y=\frac{P\cdot q}{P\cdot k},$$

$$z=\frac{P\cdot p_{\psi }}{P\cdot q},$$

$$\mathrm{d}\sigma (e+p\to J/\psi +e+X)=\int \mathrm{d}x\sum _a{f}_{a/p}(x,{\mu }_f)\mathrm{d}\sigma (e+a\to J/\psi +e+a),$$

$$\mathrm{d}\sigma (e+p\to J/\psi +e+X)$$

$$=\int \mathrm{d}x\sum _{a,n}{f}_{a/p}(x,{\mu }_f)\mathrm{d}\hat{\sigma }(e+a\to c\bar{c}[n]+e+a)⟨O^{J/\psi }(n)⟩,$$

$$\mathrm{d}\hat{\sigma }(e+a\to c\bar{c}{[}^3{P}_J^{[8]}]+e+a)\equiv \mathrm{d}\hat{\sigma }(e+a\to c\bar{c}{[}^3{P}_0^{[8]}]+e+a)$$

$$+3\mathrm{d}\hat{\sigma }(e+a\to c\bar{c}{[}^3{P}_1^{[8]}]+e+a)+5\mathrm{d}\hat{\sigma }(e+a\to c\bar{c}{[}^3{P}_2^{[8]}]+e+a).$$

$$\mathrm{d}\hat{\sigma }(e+a\to e+c\bar{c}[n]+a)$$

$$=\frac{1}{4xP\cdot k}\frac{1}{N_c{N}_s}\frac{1}{(Q^2{)}^2}{L}_{\mu \nu }{H}^{\mu \nu }[n]\mathrm{d}\Phi ,$$

$$\mathrm{d}\Phi$$

$$S=(P+k{)}^2=2P\cdot k,\hat{s}=(p+q{)}^2,$$

$$s=\hat{s}+Q^2=2p\cdot q.$$

$$p^{\mu }=(x{E}_p^{\star },0,0,-x{E}_p^{\star }),$$

$$q^{\mu }=(q_0^{\star },0,0,E_p^{\star }),$$

$$p_{\psi }^{\mu }=(\frac{z{W}^2+m_t^{\star 2}/z}{2W},p_t^{\star },0,\frac{z{W}^2-m_t^{\star 2}/z}{2W}),$$

$$k^{\mu }=(E_k^{\star },k_t^{\star }\cos {\psi }^{\star },k_t^{\star }\sin {\psi }^{\star },k_l^{\star }),$$

$$E_p^{\star }=\frac{W^2+Q^2}{2W},q_0^{\star }=\frac{W^2-Q^2}{2W},$$

$$m_t^{\star }=\sqrt{p_t^{\star 2}+M^2},x=\frac{s}{W^2+Q^2},$$

$$E_k^{\star }=\frac{S-Q^2}{2W},k_l^{\star }=\frac{1}{2W}(Q^2+\frac{W^2-Q^2}{W^2+Q^2}S),k_t^{\star }=\frac{Q}{y}\sqrt{1-y}.$$

$$z=\frac{m_t^{\star }{e}^{y_{\psi }^{\star }}}{W},$$

$$m_t^2\equiv p_t^2+M^2=\frac{4(P\cdot p_{\psi })(k\cdot p_{\psi })}{S}=2yz(k\cdot p_{\psi }),$$

$$z=\frac{2E_p{m}_t}{W^2+Q^2}{e}^{-y_{\psi }},$$

$$k\cdot p_{\psi }=\frac{1}{2yz}[m_t^{\star 2}+(1-y)z^2{Q}^2-2z\sqrt{1-y}Q{p}_t^{\star }\cos {\psi }^{\star }].$$

$$p_t^2=p_t^{\star 2}+z^2(1-y)Q^2-2z\sqrt{1-y}Q{p}_t^{\star }\cos {\psi }^{\star },$$

$$y_{\psi }=\frac{1}{2}\ln [\frac{m_t^{\star 2}+z^2(1-y)Q^2-2z\sqrt{1-y}Q{p}_t^{\star }\cos {\psi }^{\star }}{4y^2{z}^2{E}_l^2}],$$

$$L_{\mu \nu }$$

$$l^{\mu \nu }=A_g(-g^{\mu \nu }-\frac{q^{\mu }{q}^{\nu }}{Q^2})+A_L{ϵ}_L^{\mu }{ϵ}_L^{\nu }+A_{LT}({ϵ}_L^{\mu }{ϵ}_T^{\nu }+{ϵ}_T^{\mu }{ϵ}_L^{\nu })+A_T{ϵ}_T^{\mu }{ϵ}_T^{\nu },$$

$${ϵ}_L=\frac{1}{Q}(q+\frac{2Q^2}{s}p),$$

$${ϵ}_T=\frac{1}{p_t^{\star }}(p_{\psi }-\rho p-zq),$$

$$\rho =\frac{m_t^{\star 2}/z+z{Q}^2}{s},$$

$$A_g=1+\frac{2(1-y)}{y^2}-\frac{2(1-y)}{y^2}\cos (2{\psi }^{\star }),$$

$$A_L=1+\frac{6(1-y)}{y^2}-\frac{2(1-y)}{y^2}\cos (2{\psi }^{\star }),$$

$$A_{LT}=\frac{2(2-y)}{y^2}\sqrt{1-y}\cos ({\psi }^{\star }),$$

$$A_T=\frac{4(1-y)}{y^2}\cos (2{\psi }^{\star }).$$

$$l^{\mu \nu }=C_1(-g^{\mu \nu })+C_2{p}^{\mu }{p}^{\nu }+C_3\frac{p^{\mu }{p}_{\psi }^{\nu }+p_{\psi }^{\mu }{p}^{\nu }}{2}+C_4{p}_{\psi }^{\mu }{p}_{\psi }^{\nu },$$

$$C_1=A_g,$$

$$C_2=\frac{4Q^2}{s^2}(A_L-2\beta A_{LT}+{\beta }^2{A}_T),$$

$$C_3=\frac{4Q}{p_t^{\star }s}(A_{LT}-\beta A_T),$$

Full text

I Introduction

Since the nonrelativistic QCD (NRQCD) framework was proposed to solve the $J/\psi$ and $\psi(2S)$ surplus puzzle in 1994 Bodwin et al. (1995), the $J/\psi$ production in various processes has been studied within this framework. For some of these processes, QCD next-to-leading order (NLO) results have also been achieved, including the $J/\psi$ production in $e^{+}e^{-}$ annihilation Zhang et al. (2006); Zhang and Chao (2007); Gong and Wang (2008a, b); Ma et al. (2009); Gong and Wang (2009a, b); Zhang et al. (2010); Wang et al. (2011); Feng et al. (2017), the $J/\psi$ photoproduction in $e^{+}e^{-}$ Klasen et al. (2005) and $ep$ Kramer (1996); Maltoni et al. (1998); Artoisenet et al. (2009); Chang et al. (2009); Li and Chao (2009); Butenschoen and Kniehl (2010); Butenschoen and Kniehl (2011a); Bodwin et al. (2015) collisions, and the $J/\psi$ hadroproduction Campbell et al. (2007); Gong and Wang (2008c, d); Gong et al. (2009); Ma et al. (2011a); Butenschoen and Kniehl (2011b); Ma et al. (2011b); Butenschoen and Kniehl (2011c); Butenschoen and Kniehl (2012); Chao et al. (2012); Gong et al. (2013a, b); Lansberg and Shao (2013); Li and Wang (2014); Bodwin et al. (2014); Lansberg and Shao (2015); Shao et al. (2015); Bodwin et al. (2016). The production of some other species of charmonia, such as $\eta_{c}$ Butenschoen et al. (2015); Han et al. (2015); Zhang et al. (2015); Gong et al. (2016) and $\chi_{c}$ Ma et al. (2011c); Li et al. (2011); Shao and Chao (2014); Shao et al. (2014); Zhang et al. (2016), was also investigated at QCD NLO level, which provided an alternative phenomenological test of NRQCD. Among the $J/\psi$ production processes for which experimental data are available, deeply inelastic scattering (DIS) is an interesting yet complicated one. Abundant data have been released by the H1 Adloff et al. (1999, 2002); Aaron et al. (2010) and ZEUS Chekanov et al. (2005) Collaborations, however, only the QCD leading order (LO) results have been given Baier and Ruckl (1982); Korner et al. (1982); Guillet (1988); Merabet et al. (1994); Fleming and Mehen (1998); Yuan and Chao (2001); Kniehl and Zwirner (2002). Even at QCD LO, the existing phenomenological results cannot coincide with each other (see e.g. Brambilla et al. (2004)). However, the $J/\psi$ production in deeply inelastic scattering is an excellent laboratory for the study of the $J/\psi$ production mechanism. For one thing, perturbation theory works better and the resolved photon contributions are less important in large $Q^{2}$ region, relative to the $J/\psi$ photoproduction. For another, multiple distributions are measured, which can provide reference to distinguish different models.

In Reference Zhang and Sun (2017) we pointed out that most of the existing calculations on the $J/\psi$ production in DIS are based on a formalism that is not valid when the ranges of the physical observables, such as the transverse momentum ($p_{t}$) or the rapidity $y_{\psi}$ of the $J/\psi$, in the laboratory frame do not cover all their possible values. Actually, In early 1980s, the azimuthal dependence of the $J/\psi$ production in DIS has already been studied in both the colour-evaporation Duke and Owens (1981) and the colour-singlet (CS) Korner et al. (1982) models. Unfortunately, some of the succeeding investigations were not aware of this, and abused the azimuthal symmetry under an unsuitable circumstance. Besides, after the most recent phenomenological studies, three independent measurements were published. All these facts suggest that an up-to-date theoretical study of the $J/\psi$ production in DIS in the NRQCD framework is needed.

This paper is devoted to the phenomenological investigation of the $J/\psi$ production in DIS within the NRQCD framework at QCD LO. By exploiting the equations presented in Reference Zhang and Sun (2017), we obtain the comprehensive analytical formalism for the calculation of the $J/\psi$ production in DIS, which is also valid for calculating the $p_{t}^{2}$ and $y_{\psi}$ distributions of the $J/\psi$ production cross sections. One of the crucial issue of the NRQCD research is the determination of the long-distance matrix elements (LDMEs), of which there are several independent extractions, with different strategy, having obtained different results. In order to study the aftermath of the LDME uncertainties, three typical sets of these parameters are employed to present our numerical results. The rest of this paper is organised as follows. In Section II, we briefly describe the analytic formalism of our computation, following which the numerical results are given in Section III. At the end of this paper is a concluding remark in Section IV.

II Analytical Framework

II.1 General Formalism

In electron-proton ($ep$) DIS, the kinematics of the scattered lepton can be described in terms of any two of the following variables,

[TABLE]

where $P$, $k$, $k^{\prime}$ and $q$ are, as illustrated in Figure 1, the momenta of the initial proton, the initial and final lepton, and the virtual boson, respectively. Throughout this paper, we work in the limit $m_{l}\rightarrow 0$ and $m_{p}\rightarrow 0$, where $m_{l}$ and $m_{p}$ are the mass of the initial lepton and proton, respectively. To describe the kinematics of the $J/\psi$, we need to introduce two additional variables, namely the $J/\psi$ transverse momentum ($p_{t}$) and the inelasticity coefficient,

[TABLE]

where $p_{\psi}$ is the momentum of the $J/\psi$ meson.

In collinear factorization, the $J/\psi$ producton cross section at QCD LO in $ep$ DIS can be expressed as

[TABLE]

where $a$ runs over all the species of partons, the mass of which is below the $c$-quark mass, namely $g$, $u$, $d$, $s$ and the corresponding anti-quarks. $f_{a/p}(x,\mu_{f})$ is the parton distribution function (PDF) evaluated at the factorization scale $\mu_{f}$, where the momentum of the parton is $p=xP$. Note that only for the inclusive DIS at QCD LO, when the invariant mass of the hadronic final states is zero, $x$ is identical to the Bjorken-$x$, $x_{B}$, defined in Equation 1.

Within the NRQCD framework, the partonic cross section can be further factorized, accordingly, the cross section defined in Equation 3 can be written as Bodwin et al. (1995)

[TABLE]

where $\langle O^{J/\psi}(n)\rangle$ is the LDME, which describes the hadronisation of a $c\bar{c}$ pair with quantum number $n$, and $\mathrm{d}\hat{\sigma}$ is the corresponding short-distance coefficient (SDC). For the $J/\psi$ production, $n$ can be ${}^{3}S_{1}^{[1]}$, ${}^{1}S_{0}^{[8]}$, ${}^{3}S_{1}^{[8]}$ and ${}^{3}P_{J}^{[8]}$ up to the order of $v^{4}$, where $v$ denotes the typical relative velocity of the $c\bar{c}$ pair inside the $J/\psi$ meson. Taking advantage of the relations $\langle O^{J/\psi}(^{3}P_{2}^{[8]})\rangle=5\langle O^{J/\psi}(^{3}P_{0}^{[8]})\rangle$, and $\langle O^{J/\psi}(^{3}P_{1}^{[8]})\rangle=3\langle O^{J/\psi}(^{3}P_{0}^{[8]})\rangle$, we can synthesise the three SDCs for $n=^{3}P_{0}^{[8]}$, $n=^{3}P_{1}^{[8]}$, and $n=^{3}P_{2}^{[8]}$ by defining

[TABLE]

The corresponding LDME for this synthesised SDC thus should be $\langle O^{J/\psi}(^{3}P_{0}^{[8]})\rangle$.

The partonic SDCs can in general be written as

[TABLE]

where $1/(N_{c}N_{s})$ is the color and spin average factor, $L_{\mu\nu}$ and $H^{\mu\nu}[n]$ are the leptonic and hadronic tensors, respectively, and

[TABLE]

where $p_{a}$ is the momentum of the final state parton. Note that, here, the LDMEs have been eliminated from the hadronic tensors.

In the following subsections, we will present the explicit form of the elements needed in our calculation in Equation 6.

II.2 The Parameterisation of the Physical Variables

In this subsection, we provide the expressions for the physical variables needed in our calculation.

The squared colliding energy at the lepto-hadronic level and the partonic level are defined as

[TABLE]

respectively. However, in most of the cases, $\hat{s}+Q^{2}$ emerges as a unity, thus, we define

[TABLE]

Our calculations are carried out in the virtual-boson-proton ($\gamma p$) centre-of-mass frame. All the quantities measured in this frame are labeled by a superscript $\star$ in this paper. Assigning the $J/\psi$ mass as $M$, the momenta involved can be parameterised as

[TABLE]

where

[TABLE]

The momentum of the $J/\psi$ can be obtained through the following relation:

[TABLE]

where $y_{\psi}^{\star}$ is the rapidity of the $J/\psi$ measured in the $\gamma p$ rest frame. Regarding the sign of the longitudinal component of $q$, the forward $z$ direction is defined as that of the incident virtual photon, which is consistent with the HERA convention.

For some of the available data, the ranges of the transverse momentum ($p_{t}$) or the rapidity ($y_{\psi}$) of the $J/\psi$ meson in the laboratory frame are also specified, thus, we need to find the relations between these two variables and $\psi^{\star}$111 The relation between $p_{t}$ and $\psi^{\star}$ can also be obtained in the proton rest frame by applying a rotation. However, here we introduce an alternative approach, which can avoid implementing the limit $m_{p}\rightarrow 0$. . In the laboratory frame, the rapidity of the $J/\psi$ is calculated with respect to the incident proton direction. Accordingly, we can obtain the following equations:

[TABLE]

where $E_{p}$ is the energy of the incident proton in the laboratory frame. After a short calculation, we arrive at

[TABLE]

Then we can obtain the following relations:

[TABLE]

where $E_{l}$ is the energy of the incident lepton in the laboratory frame.

II.3 The Calculation of the Leptonic and Hadronic Tensors

In HERA experimental condition, $Q$ is much smaller than the mass of the $Z_{0}$ boson, thus, we can neglect the contributions from the $Z_{0}$ propagator. The leptonic tensor can be obtained as

[TABLE]

The normalised leptonic tensor, $l_{\mu\nu}$, can be decomposed into the linear combination of four independent Lorentz invariant structures as Zhang and Sun (2017)

[TABLE]

where222We have replaced $P$ in the formalism presented in Reference Zhang and Sun (2017) by $p$.

[TABLE]

and

[TABLE]

The leptonic tensor can be modified to be more convenient for our computation as

[TABLE]

where

[TABLE]

with

[TABLE]

Then the contraction of the leptonic tensor with the hadronic one can be expressed as

[TABLE]

where

[TABLE]

and

[TABLE]

In the derivation of the above equations, we have exploited the relations, $q_{\mu}H^{\mu\nu}=0$ and $H_{\mu\nu}=H_{\nu\mu}$.

The analytical results for $H_{i}[n]$ ($i=$1, 2, 3, 4) are presented in Appendix A.

II.4 Phase Space

The only missing element for calculating the SDCs in Equation 6 is the expressions for the phase space in Equation 7.

The phase space for the scattered lepton can be obtained as

[TABLE]

while that for the hadrons is

[TABLE]

Note that Equation 27 is valid in any frame, thus, we do not restrict it in the $\gamma p$ rest frame, and erase all the superscript $\star$. To keep the Lorentz invariance of our formalism, we integrate over $\mathrm{d}x$ to eliminate the last dimension of the $\delta$ function and obtain

[TABLE]

where the value of $x$ can be found in Equation 11. Then we arrive at the final expression of the phase space,

[TABLE]

II.5 The Cross Section

The cross section for the $J/\psi$ leptoproduction can be obtained by combining Eqs. (4), (6), (23) and (29), and reduces to

[TABLE]

III Numerical Results

Before our work, three papers Fleming and Mehen (1998); Yuan and Chao (2001); Kniehl and Zwirner (2002) studying the $J/\psi$ production in DIS within the NRQCD framework have already been published. It is worth noting that they are mutually incompatible with each other. With the same leptonic tensors, our results agree with those in Reference Kniehl and Zwirner (2002). However, as we pointed out in Reference Zhang and Sun (2017), Reference Kniehl and Zwirner (2002) used a form of the leptonic tensor identical to

[TABLE]

which will lead to wrong results when $p_{t}$ or $y_{\psi}$ in the laboratory frame do not cover all their possible values. Another important issue to address is that, to calculate the polarisation of the $J/\psi$ meson, an additional momentum other than $q$ and $p$ will also emerge in the hadronic tensor, thus, with the leptonic tensor in Equation 31, one will also obtain wrong results. To illuminate the difference between Equation 20 and Equation 31, we present the comparison of our results with those in Table 1 of Reference Kniehl and Zwirner (2002). Applying the same parameter choices, namely, $m_{c}=1.5\mathrm{~{}GeV}$, $E_{l}=27.5\mathrm{~{}GeV}$, $E_{p}=820\mathrm{~{}GeV}$, $30\mathrm{~{}GeV}<W<150\mathrm{~{}GeV}$, $\alpha=1/137$, $\mu_{r}=\mu_{f}=\sqrt{Q^{2}+M^{2}}$, $\langle O^{J/\psi}(^{3}S_{1}^{[1]})\rangle=1.1\mathrm{~{}GeV}^{3}$, $\langle O^{J/\psi}(^{1}S_{0}^{[8]})\rangle=1\times 10^{-2}\mathrm{~{}GeV}^{3}$, $\langle O^{J/\psi}(^{3}S_{1}^{[8]})\rangle=1.12\times 10^{-2}\mathrm{~{}GeV}^{3}$, $\langle O^{J/\psi}(^{3}P_{0}^{[8]})\rangle/m_{c}^{2}=5\times 10^{-3}\mathrm{~{}GeV}^{3}$, the LO proton PDF given in Reference Gluck et al. (1995), and the one-loop $\alpha_{s}$ running with $n_{f}=4$ and $\Lambda^{(4)}=0.13\mathrm{~{}GeV}$, we present the corresponding results in Table 1, where our results for employing the leptonic tensors presented in Equation 31 and Equation 20 are labeled by I and II, respectively. The slight difference between our results I and those in Reference Kniehl and Zwirner (2002) might be due to the difference in the precision of the parameter choices kept in the two independent evaluations. However, when the $p_{t}^{2}$ cut is applied, the difference between our result II and that in Reference Kniehl and Zwirner (2002) is significant, which manifests the necessity of the inclusion of the azimuthal-asymmetric terms in the calculation.

Now we are in a position to confront our phenomenological results with the H1 Adloff et al. (1999, 2002); Aaron et al. (2010) and ZEUS Chekanov et al. (2005) data. In the rest part of this paper, we take $m_{c}=1.5$ GeV, $\alpha=1/137$ and $E_{l}=27.5\mathrm{~{}GeV}$. For the experimental conditions in Reference Adloff et al. (1999), $E_{p}=820\mathrm{~{}GeV}$, while for those in References Adloff et al. (2002); Chekanov et al. (2005); Aaron et al. (2010), $E_{p}=920\mathrm{~{}GeV}$. The one-loop $\alpha_{s}$ running is employed, where its value at the $Z_{0}$ boson mass is set to be $\alpha_{s}(M_{Z})=0.13$. Correspondingly, we employ CTEQ6L1 Pumplin et al. (2002) as the PDF for the proton. Our default set of the renormalisation and factorisation scales is $\mu_{r}=\mu_{f}=\mu_{0}\equiv\sqrt{Q^{2}+M^{2}}$. To present the uncertainties of our theoretical predictions, the bands in the plots cover all the results for simultaneously varying these two scales from $0.5\mu_{0}$ to $2\mu_{0}$. The NRQCD predictions suffer severely from the uncertainties of the LDMEs, thus, we will adopt three sets of the LDMEs on the market, the detailed values of which can be found in Table 2, to present our numerical results333 With the recent data of the hadroporduction of $\eta_{c}$ from the LHCb Collaboration, the authors of Reference Han et al. (2015) have imposed a strong upper limit on $\langle O^{J/\psi}(^{1}S_{0}^{[8]})\rangle$. Since the values of their new set of LDMEs almost coincide with the ones in Reference Zhang et al. (2015), in our current paper, only the latter set is adopted. For the same reason, the LDMEs presented in Reference Sun and Zhang (2015), which is an update of Reference Zhang et al. (2015) to adapt to the $J/\psi$ polarization measurements, are also ignored in the presentation of our results. To understand the $J/\psi$ polarisation puzzle, the authors of Reference Chao et al. (2012) proposed the ${}^{1}S_{0}^{[8]}$ dominance picture. There are some other works Gong et al. (2013a); Faccioli et al. (2014); Bodwin et al. (2016) studying the $J/\psi$ polarisation, which, although with different fitting strategies, also arrived at similar results, namely the ${}^{1}S_{0}^{[8]}$ LDME being almost one order of magnitude larger than the other two. Applying these LDMEs, the results for the $J/\psi$ production in DIS at HERA do not differ significantly from those obtained by employing the LDMEs in Reference Chao et al. (2012). Therefore, we just adopt the three representative sets of the LDMEs presenting our results.. The other origin of the uncertainties in our results is those arising from the LDMEs. Note that the $c$-quark mass ($m_{c}$) dependence of the SDCs can generally be balanced out by that of the LDMEs (see e.g. Ma et al. (2011b)). In this paper, we do not vary the value of $m_{c}$.

We notice that the data for the $J/\psi$ leptoproduction at HERA is plentiful. In some of the experimental conditions, the perturbative calculation based on the collinear factorisation is not expected to work well, therefore, we need first to address the ranges of the variables for good perturbative expansion. The cross sections diverge at $Q^{2}=0$, thus, when $Q^{2}$ is too small, the results will suffer from a divergence factor, $1/Q^{2}$. The newest HERA measurements applied the cut, $Q^{2}<2.5\mathrm{~{}GeV}^{2}$, for photoproduction, according to which, $Q^{2}=2.5\mathrm{~{}GeV}^{2}$ can be considered as a moderate boundary to separate regions for good and bad perturbative calculations. For a specific value of $p_{t}^{\star 2}$ and $p_{t}^{2}$, the minimum values of $x$ are approximately $[(m_{t}^{\star}+p_{t}^{\star})^{2}+Q^{2}]/(W^{2}+Q^{2})$ and $m_{t}^{2}/[z(W^{2}+Q^{2})]$, respectively. Thus, in small $p_{t}$ and $p_{t}^{\star}$ regions, the gluon saturation effect might ruin our results. The perturbation theory works better in the region, $p_{t}^{2}>6.4\mathrm{~{}GeV}^{2}$ or $p_{t}^{\star 2}>6.4\mathrm{~{}GeV}^{2}$, where the value of $x$ is generally larger than 0.001. Since the selection rules do not forbid the processes $\gamma g\rightarrow c\bar{c}[n]$ for $n=^{1}S_{0}^{[8]}$ and $n=^{3}P_{J}^{[8]}$, both soft and collinear singularities will arise for the $c\bar{c}[^{1}S_{0}^{[8]}]$ and $c\bar{c}[^{3}P_{J}^{[8]}]$ production at $z=1$. One will suffer from a divergence factor, $1/(1-z)$, in the region, $z\rightarrow 1$, as a result, only when $z$ is much smaller than 1, the NRQCD predictions are expected to be reliable. In the region, $z>0.6$, this divergence factor is larger than 2.5, which might significantly enhances the terms which ought to be small. For this reason, we consider the perturbative calculation in the region, $z>0.6$, to be not reliable. In addition, around $z=1$, the diffractively produced $J/\psi$ will dominate the prompt $J/\psi$ production, which, however, cannot be effectively eliminated from the $J/\psi$ inelastic production events. In the region $M_{X}>10\mathrm{~{}GeV}$, the value of $z$ can be as large as 0.997, thus, we completely omit the data corresponding to this condition.

As follows, we give the results for the differential cross sections with respect to $p_{t}^{2}$, $p_{t}^{\star 2}$, $Q^{2}$, $W$, $y_{\psi}$, $y_{\psi}^{\star}$, and $z$, and the double differential cross sections with respect to $p_{t}^{\star 2}$ and $z$, and $Q^{2}$ and $z$. The figures for each observed physical variable are grouped according to the experiment papers that the corresponding data are taken from. All these numerical results are presented in Figures 2-18. For each experimental condition, there are three plots, corresponding to the three sets of the LDMEs given in References Chao et al. (2012); Butenschoen and Kniehl (2011c); Zhang et al. (2015).

One can easily find that the theoretical results via the CS mechanism at LO are generally below the data. However, the discrepancy is not so large as that in the $J/\psi$ hadroproduction cases. The largest discrepancy emerges in high $p_{t}^{2}$, $p_{t}^{\star 2}$ and $Q^{2}$ regions, where the perturbative calculations are credible. In these regions, the CS results are one order of magnitude smaller than the data. According to Reference Sun and Zhang (2017), the QCD corrections to the CS channel is minor, to this end, the CO contributions are necessary for understanding the $J/\psi$ leptoproduction data.

To quantitatively investigate the agreement between the theoretical results and data, we present the $\chi^{2}$/d.o.f. for each set of the LDMEs in Table 3, where two cut conditions are considered. For the condition, $z<0.9$, all the data points calculated in this paper are included, while for $z<0.6$, only those satisfying $z<0.6$ are counted. We can see that the LDMEs given in Reference Butenschoen and Kniehl (2011c) work equally well in both large $z$ and small $z$ regions, while the other two sets fail in the region, $0.6<z<0.9$, where a divergence factor $1/(1-z)$ will make the theoretical results larger than the data. In the region, $z<0.6$, the LDMEs given in Reference Zhang et al. (2015) work best. However, the differences of the $\chi^{2}$/d.o.f. among the results obtained by using the three sets of the LDMEs are not significant. Besides, the small-$x$ effects and threshold resummation are not considered in our calculation. Therefore, it is too early to judge which of the three is better. Once these two challenging works are finished, the $J/\psi$ production in DIS can provide a good reference to distinguish the quarkonium production mechanisms.

IV Summary

In this paper, we studied the $J/\psi$ production within the NRQCD framework at QCD LO. Although this process has already been investigated in three previous papers Fleming and Mehen (1998); Yuan and Chao (2001); Kniehl and Zwirner (2002), their calculations employed a form of the leptonic tensor which will lead to wrong results when the $p_{t}^{2}$ and $y_{\psi}$ in the laboratory frame do not cover all values possible for them. Since many of the existing data are measured at specific values of $p_{t}^{2}$ or $y_{\psi}$, and even for the $Q^{2}$, $W$, $y_{\psi}^{\star}$ and $z$ distributions, some measurements applied a $p_{t}^{2}$ cut, our renewed investigation is necessary. The calculation in this paper are based on a formalism proposed in Reference Zhang and Sun (2017), which significantly reduces the complexity of the computation by reducing the number of the involved momenta.

We presented in this paper the numerical results comparing with all the existing data released from HERA collaborations. Three representative sets of the LDMEs were employed for comparison. We found that the CS contributions at QCD LO are generally below the data. In the regions where the perturbative calculations are expected to work well, the CS results are almost one order of magnitude smaller than the HERA data. The NRQCD results obtained by using the three sets of the LDMEs have their own advantages in the description of the data in specific kinematical regions. In general, the results using the LDMEs in Reference Chao et al. (2012) are the largest among the three, while those for Reference Butenschoen and Kniehl (2011c) are the smallest. Since the gluon saturation effect might be important, and the behaviour of the $z$ distribution of the ${}^{1}S_{0}^{[8]}$ and ${}^{3}P_{J}^{[8]}$ SDCs might be completely changed if the $1/(1-z)$ singularities are smeared, it is not proper to make definite conclusions until these two works are accomplished. Besides, the LDMEs employed in our calculation are all extracted with the QCD NLO SDCs. When the QCD corrections to all the CO channels are achieved, the comparison of the theoretical results to data can provide better references for the universality of the LDMEs.

In any sense, the $J/\psi$ production in DIS provides an alternative device for the study of the $J/\psi$ production mechanism. Benefiting from the various observables, they can provide new references for the theoretical studies. They may also help to fix the LDMEs for the $J/\psi$ production when new progress in the phenomenology are made, and new experiments in the future colliders, such as the EIC, are carried out.

Acknowledgements.

This work is supported by the National Natural Science Foundation of China (Grant No. 11405268, 11647113 and 11705034).

Appendix A Analytical expressions for $H_{1}$, $H_{2}$, $H_{3}$ and $H_{4}$

In this appendix, we present the analytical expressions for $H_{1}$, $H_{2}$, $H_{3}$ and $H_{4}$. The average on the initial spin and colour are not implemented. To make the expression more compact, we define

[TABLE]

the relations between the variables defined in Equation 32 and the ordinarily used Mandelstam variables,

[TABLE]

are

[TABLE]

Thus we have

[TABLE]

which we have employed to eliminate $M^{2}$ in our expressions. In the following, $e_{q}$ denotes the fractional electric charge of quark $q$.

The expressions are listed below.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]