QCD Corrections to the Color-singlet $J/\psi$ Production in Deeply Inelastic Scattering at HERA

TL;DR

This paper calculates QCD corrections to J/psi production in deep inelastic scattering at HERA, showing that corrections are significant but still insufficient to match experimental data, implying the need for color-octet contributions.

Contribution

First calculation of QCD corrections to color-singlet J/psi production in DIS at HERA, highlighting the importance of color-octet mechanisms.

Findings

QCD K-factors range from 0.85 to 2.38

Perturbative expansion converges well at low transverse momentum

Color-singlet cross section remains below experimental data

Abstract

We present the first study of the QCD corrections to the color-singlet (CS) $J/\psi$ production in deeply inelastic $ep$ scattering at HERA. The $K$-factor ranges from $0.85$ to $2.38$ in the kinematic regions we study. In low transverse momentum regions, the $K$-factors is even smaller, and close to 1, which indicates good convergence of the perturbative expansion. With the QCD corrections, the CS cross section is still below the data. At least at QCD next-to-leading order, the color-octet mechanism is necessary to describe the data.

Equations

$$Q^2=-q^2,W^2=(P+q{)}^2,z=\frac{P\cdot p_{\psi }}{P\cdot q},$$

$$S=2P\cdot k,s=2p\cdot q.$$

$$\mathrm{d}\hat{\sigma }(e+p\to e+c\bar{c}{[}^3{S}_1^{[1]}]+X)$$

$$=\frac{\alpha }{(4\pi {)}^3{N}_c{N}_s{S}^2}\sum _n\frac{1}{x}{f}_{a_n/p}(x,{\mu }_f)\sum _{m=1}^4{C}_m{H}_m^n$$

$$\times \frac{\mathrm{d}{Q}^2}{Q^2}\frac{\mathrm{d}{W}^2}{W^2+Q^2}\mathrm{d}{p}_t^2\frac{\mathrm{d}z}{z(1-z)}\mathrm{d}\psi \mathrm{d}{\Phi }_X^n,$$

$$\mathrm{d}{\Phi }_X^n$$

$$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),C_4=\frac{1}{p_t^{\star 2}}{A}_T,$$

$$\beta =\frac{p_t^{\star 2}+M^2+z^2{Q}^2}{2z{p}_t^{\star }Q},$$

$$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{6(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 }).$$

$$H_1=-g^{\mu \nu }{H}_{\mu \nu },H_2=p^{\mu }{p}^{\nu }{H}_{\mu \nu },$$

$$H_3=p^{\mu }{p}_{\psi }^{\nu }{H}_{\mu \nu },H_4=p_{\psi }^{\mu }{p}_{\psi }^{\nu }{H}_{\mu \nu },$$

$$\mathrm{d}\sigma =⟨\mathcal{O}⟩\mathrm{d}\hat{\sigma }(e+p\to e+c\bar{c}{[}^3{S}_1^{[1]}]+g),$$

$$H_m^S\equiv {\int }_{softregion}\mathrm{d}{\Phi }_{ij}{H}_m^{rc},$$

$$H_m^{HC}\equiv {\int }_{hardcollinearregion}\mathrm{d}{\Phi }_{ij}{H}_m^{rc}.$$

$$H_m^{(2)}=H_m^{Born}+H_m^V+H_m^S+H_m^{HC},$$

$$\mathrm{d}\hat{\sigma }=\frac{\alpha }{(4\pi {)}^3{N}_c{N}_s{S}^2}\frac{\mathrm{d}{Q}^2}{Q^2}\frac{\mathrm{d}{W}^2}{W^2+Q^2}\mathrm{d}{p}_t^2\frac{\mathrm{d}z}{z(1-z)}\mathrm{d}\psi$$

$$\times \frac{1}{x}{f}_{a_n/p}(x,{\mu }_f)\sum _{m=1}^4{C}_m[H_m^{(2)}+{\int }_{H\overline{C}}\mathrm{d}{\Phi }_{ij}{H}_m^{rc}],$$

Full text

Since the discovery of the $J/\psi$ meson Aubert et al. (1974); Augustin et al. (1974), the production mechanism of heavy quarkonia has been a hot issue in high energy physics. Before the first measurement of the $J/\psi$ and $\psi(2s)$ hadroproduction carried out by the CDF Collaboration Abe et al. (1992), the color-singlet (CS) model Einhorn and Ellis (1975); Ellis et al. (1976); Carlson and Suaya (1976); Chang (1980); Berger and Jones (1981); Baier and Ruckl (1981, 1982) was generally accepted as a natural description of the heavy quarkonia production and decay mechanism. In 1994, the nonrelativistic QCD (NRQCD) effective theory Bodwin et al. (1995) was proposed, which successfully filled the huge gap between the QCD leading order (LO) predictions via the CS model and the CDF measurement of the $J/\psi$ and $\psi(2s)$ hadroproduction Braaten and Fleming (1995); Cho and Leibovich (1996a, b). At QCD LO in the NRQCD framework, the dominant hadroproduction mechanism of the $J/\psi$ mesons is the gluon fragmentation into a color-octet (CO) ${}^{3}S_{1}$ $c\bar{c}$ pair, which produces a $J/\psi$ by emitting soft gluons in a long-distance process. Although the CS model also permits a gluon fragmenting into a $J/\psi$ with additional two hard gluons emitted, according to NRQCD, this mechanism substantially underestimates the production rate. In addition to the $J/\psi$ hadroproduction, NRQCD also worked well in many other processes, including the $J/\psi$ production in $\gamma\gamma$ fusion Klasen et al. (2002), the $J/\psi$ photoproduction at HERA Butenschoen and Kniehl (2010), the $\chi_{c}$ meson hadroproduction Ma et al. (2011); Zhang et al. (2016), and etc. However, one cannot overlook the controversies NRQCD is facing, among which the $J/\psi$ polarization puzzle is the most well known and challenging one. Many independent studies Butenschoen and Kniehl (2012); Chao et al. (2012); Gong et al. (2013); Shao and Chao (2014); Shao et al. (2014); Bodwin et al. (2014); Shao et al. (2015); Bodwin et al. (2016) at QCD next-to-leading order (NLO) level have been performed. Those who achieved good description of the $J/\psi$ polarization data are generally consistent with the ${}^{1}S_{0}^{[8]}$ dominance picture. However, the theoretical studies on the $\eta_{c}$ hadroproduction indicate that Butenschoen et al. (2014); Han et al. (2015); Zhang et al. (2014), regarding the heavy quark spin symmetry, this picture violates the recent LHCb measurement Aaij et al. (2015). This paradox was remedied in Reference Sun and Zhang (2015) which found that, even without the ${}^{1}S_{0}^{[8]}$ dominance picture, the $J/\psi$ polarization data can also be understood within the NRQCD framework.

To solve these problems, a new factorization theory was proposed in Reference Ma and Chao (2017). From another angle of view, some researchers challenge the significance of the CO contributions, and seek the way of describing the data within the CS framework. The $J/\psi$ production in $e^{+}e^{-}$ annihilation provides an example of success of this idea. With the QCD and relativistic corrections Zhang and Chao (2007); Ma et al. (2009); Gong and Wang (2009a, b); He et al. (2010), the CS contributions almost saturate the Belle data Pakhlov et al. (2009), while the inclusion of the CO ones will generally ruin the agreement between theory and experiment. Looking back at the $J/\psi$ hadroproduction, QCD NLO corrections Campbell et al. (2007) enhance the differential cross sections for the CS $J/\psi$ hadroproduction in medium and high transverse momentum ($p_{t}$) regions by one to two orders of magnitude, which reduces the discrepancy between theory and data, at the same time, change the $J/\psi$ polarization from transverse to longitudinal Gong and Wang (2008). Due to the lack of the complete next-to-next-to-leading order results, one cannot yet draw definite conclusions on the significance of the CO contributions implied by the $J/\psi$ hadroproduction data (as a review, see e.g. Reference Brambilla et al. (2011)). Such large $K$-factors are due to the fact that at QCD LO, both the leading power (LP) and next-to-leading power (NLP) terms vanish, and at QCD NLO, the NLP behaviour arises; this new behaviour enhances the cross sections significantly, especially in high $p_{t}$ regions.

The $J/\psi$ production in deeply inelastic scattering (DIS), which is also called the $J/\psi$ leptoproduction, can serve as another test of the quarkonium production mechanisms. The CS $J/\psi$ photoproduction has been studied at $\mathcal{O}(\alpha\alpha_{s}^{3})$ in Refs. Kramer et al. (1995); Kramer (1996); Artoisenet et al. (2009); Chang et al. (2009); Li and Chao (2009), which found that with the QCD corrections the CS contributions are still below the data, especially when $p_{t}$ is large. For the $J/\psi$ leptoproduction, the deflection angle of the scattered lepton is larger, accordingly the virtuality of the incident photon, which is emitted by the incident lepton and will interact in the hadronic process, need to be taken into account in the perturbative calculation. We define $Q^{2}=-q^{2}$, where $q$ is the momentum of the incident photon. Due to larger $Q^{2}$, the $J/\psi$ leptoproduction shows better features than the $J/\psi$ photoproduction. The $p_{t}$ of the $J/\psi$ yield data measured at HERA is relatively low, thus the ratio ($x$) of the proton momentum taken by the interacting parton might be very small for the $J/\psi$ photoproduction, and the gluon saturation effects Mueller and Qiu (1986) would be important. Larger $Q^{2}$ can increase the value of $x$, correspondingly suppress the gluon saturation effects. Moreover, as $Q^{2}$ increases, the contributions from the resolved photons are greatly suppressed, and the perturbative expansion in $\alpha_{s}$ will be better as well. Another interesting feature of the $J/\psi$ leptoproduction is that the NLP behaviour arises at QCD LO. Unlike the $J/\psi$ hadroproduction case, no new behaviour emerges in the $J/\psi$ production in DIS at QCD NLO, thus, the convergence of the $\alpha_{s}$ expansion should be better.

However, the emergence of $Q^{2}$ at the same time makes the computation much more complicated. We notice that although the HERA collaborations have published abundant data, the one-loop level phenomenological study of the $J/\psi$ production in DIS is still lacking. The $J/\psi$ production in DIS has been studied at QCD LO in many papers Duke and Owens (1981); 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), however, as indicated in our recent work Zhang and Sun (2017), all those Fleming and Mehen (1998); Yuan and Chao (2001); Kniehl and Zwirner (2002) under the NRQCD framework are based on a formalism that will lead to wrong results when the ranges of the $J/\psi$ $p_{t}$ or rapidity ($y_{\psi}$) are not taken to cover all their possible values. For this reason, we provided a renewed QCD LO study of the $J/\psi$ production in Reference Sun and Zhang (2017). In this paper, we study the QCD corrections to the CS $J/\psi$ production in DIS, following the calculation framework provided in Refs. Zhang and Sun (2017); Sun and Zhang (2017).

The process for the $J/\psi$ leptoproduction is illustrated in Figure 1. $k$, $k^{\prime}$, $P$, $p$ and $p_{\psi}$ are the momenta of the incident and scattered electron, the proton, the parton generated from the proton, and the produced $J/\psi$, respectively. The generally used invariants are defined as

[TABLE]

In our calculation, all the processes up to $\alpha^{2}\alpha_{s}^{3}$ are counted, including $eg\rightarrow eJ/\psi g$ at both tree and one-loop level, and $ea\rightarrow eJ/\psi ij$, in which $aij$ runs over $ggg$, $gq\bar{q}$, $qqg$, $\bar{q}\bar{q}g$, and $gc\bar{c}$. Here we use $g$ and $q(\bar{q})$ to denote a gluon and a light quark (anti-quark), respectively. Since all the singularities cancel in the hadronic process, namely $\gamma^{\star}+p\rightarrow J/\psi+X$, where $\gamma^{\star}$ denotes a virtual photon, we will directly employ the form of the leptonic tensor given in Reference Zhang and Sun (2017), which eventually leads to the same results as those by adopting the $d$-dimensional form of the leptonic tensor. The contraction between the leptonic tensor and the hadronic one can be carried out in 4-dimensions. Thus, the form of the leptonic tensor given in Reference Zhang and Sun (2017) can be directly employed without extension to its $d$-dimensional form. If we use the two-cutoff phase space slicing method Harris and Owens (2002) to separate the divergences in the real-correction processes, the phase space for the scattered electron and the $J/\psi$ can also be written in 4-dimensions. Then the short-distance coefficient (SDC) for the CS $J/\psi$ production in DIS can be expressed as Sun and Zhang (2017)

[TABLE]

where $n$ runs over all the partonic processes, $\alpha$ is the fine structure constant, $1/(N_{c}N_{s})$ is the spin and color average factor, $f_{a_{n}/p}(x,\mu_{f})$ is the parton distribution function (PDF) of a parton $a_{n}$ in a proton with $x$ and $\mu_{f}$ being the fraction of the proton momentum taken by the parton and the factorization scale, respectively, $p_{t}$ and $\psi$ are the transverse momentum of the $J/\psi$ and the azimuthal angle of the lepton plane around the $z$ axis, respectively, and for the processes $eg\rightarrow eJ/\psi g$, $\mathrm{d}\Phi_{X}^{n}=\mathrm{d}\Phi_{g}=1$, while for the real-correction processes, namely $\gamma^{\star}a\rightarrow J/\psi ij$,

[TABLE]

Note that Equation 2 is valid in any frame of reference. Since $\mathrm{d}x$ has been integrated over to eliminate the $\delta$ function which keeps the energy conservation, the value of $x$ has been fixed in Equation 2 as $x=(s+Q^{2})/(W^{2}+Q^{2})$.

The expressions for $C_{m}$ can be found in Reference Sun and Zhang (2017) and are duplicated in the following as

[TABLE]

where

[TABLE]

and

[TABLE]

Hereinafter, we denote all the physical quantities in the $\gamma^{\star}p$ center-of-mass frame by a superscript $\star$. $H_{m}$ are defined as

[TABLE]

where $H_{\mu\nu}$ is the hadronic tensor. Note that here the long-distance matrix elements (LDMEs) have been eliminated from the hadronic tensors. Then the cross section for the CS $J/\psi$ production can be expressed as

[TABLE]

where the CS LDME for the $J/\psi$ production, $\langle\mathcal{O}^{J/\psi}(^{3}S_{1}^{[1]})\rangle$, is abbreviated as $\langle\mathcal{O}\rangle$.

All the singularities are contained in $H^{\mu\nu}$’s and $\mathrm{d}\Phi_{X}$, thus, we evaluate them in $d$-dimension.

We denote $H_{m}^{Born}$ as the $H_{m}$ for the process $\gamma^{\star}g\rightarrow c\bar{c}[^{3}S_{1}^{[1]}]g$ at tree-level, $H_{m}^{V}$ as the virtual corrections to $H_{m}^{Born}$, and $H_{m}^{rc}$ as the sum of all the $H_{m}$’s for the real-correction processes and define

[TABLE]

One need to be more careful in identifying the soft and hard collinear regions, since the squared invariant mass of the incident photon is negative. The soft region is defined in terms of the energy of the final-state gluon, $E_{g}$, in the $\gamma^{\star}g$ rest frame by $E_{g}\leq\delta_{s}\sqrt{s}/2$, while the collinear region is defined by the inequality that the inner product of two massless momenta is smaller than $\delta_{c}s/2$, where $s$ is defined as in Equation 1. Here, $\delta_{s}$ and $\delta_{c}$ are two arbitrary real numbers, however, small enough to make sure that integrals of a finite function in the soft and collinear regions are negligible. Under these definitions, most of the formulas in Reference Harris and Owens (2002) remain unchanged. To avoid double counting, the gluon-soft regions need to be excluded in the calculation of the cross sections in the hard collinear regions. Accordingly, the integral domain of the ratio $z$ defined in Reference Harris and Owens (2002) should be properly determined. To distinguish this ratio and the inelasticity coefficient, we assign another symbol, $\xi$, to this ratio. When $p_{i}$ and $p_{j}$ are collinear, $\xi$ is defined as $\xi=p_{i0}/(p_{i0}+p_{j0})$. For $i=g$ and $j=q$, the integral domain of $\xi$ should be $\delta_{s}^{\prime}<\xi<1$, where $\delta_{s}^{\prime}=\delta_{s}\sqrt{s(s-Q^{2})}/(s-Q^{2}-M^{2})$, and $M$ is the invariant mass of all the hadronic final states other than $i$ and $j$ in the partonic process. In the process we study, $M$ is the $J/\psi$ mass, which is set to be twice of the $c$-quark mass, $m_{c}$, namely $M=2m_{c}$. For $i,j=g$, the integral domain of $\xi$ should be $\delta_{s}^{\prime}<\xi<1-\delta_{s}^{\prime}$. When $p$ and $p_{i}$ is collinear, the definition of $\xi$ is different, which should be $\xi=1-p_{i0}/p_{0}$. The integral domain of $\xi$ should be $x<\xi<1$ for $i=q$, and $x<\xi<1-\delta_{s}^{\prime\prime}$ for $i=g$, where $\delta_{s}^{\prime\prime}=\sqrt{(s-Q^{2})/s}\delta_{s}$. Here $x$ is the ratio of the parton momentum, $p-p_{i}$, to the incident proton momentum. The expressions of $\delta_{s}^{\prime}$ and $\delta_{s}^{\prime\prime}$ are different from those given in Reference Harris and Owens (2002) due to the nonzero $Q^{2}$. With the above configurations, the sum of $H_{m}^{V}$, $H_{m}^{S}$ and $H_{m}^{HC}$ are divergence free. Defining

[TABLE]

one can rewrite the SDC for the CS $J/\psi$ leptoproduction as

[TABLE]

where the subscript $H\overline{C}$ means the integral is carried out in the hard noncollinear region, and $\mathrm{d}\Phi_{ij}$ and $H_{m}^{rc}$ are evaluated in the limit $d\rightarrow 4$.

To evaluate the $H_{m}$’s for each process, we employ our new Mathematica package, Malt@FDC. This package can automatically reduce the loop amplitudes into linear combination of master integrals, which will be computed numerically with Looptools Hahn and Perez-Victoria (1999), and simplify the expressions of the squared amplitudes. Our final expressions contain the $\mathcal{O}(\epsilon)$ terms of the $A_{0}$ and $B_{0}$ functions, which is not given in the Looptools library. These functions are computed with a new FORTRAN package, which will be discussed elsewhere Zhang and Sun (tion). Before working on the current process, we have applied our Malt@FDC in tens of other processes. All our results are consistent with those obtained by the FDC system Wang (2004) and/or those given in published papers. As an indispensable check, we studied the $Q^{2}\rightarrow 0$ limit and compared the results with the $J/\psi$ photoproduction. Setting $Q^{2}$ to be zero, replacing the leptonic tensor associated with the virtual photon propagator by $-g_{\mu\nu}$, and implementing proper phase-space integration, we can reproduce the $J/\psi$ photoproduction results in Refs. Chang et al. (2009); Li and Chao (2009). Another check to mention is that our results are independent of $\delta_{s}$ and $\delta_{c}$ and gauge invariant.

Abundant $J/\psi$ leptoproduction data Adloff et al. (1999, 2002); Aaron et al. (2010); Chekanov et al. (2005) have been measured at HERA. However, most of them lie in the kinematic regions where the perturbative calculation might not be good. In order to make the perturbation theory work better, we constrain our concerns in the regions where $p_{t}^{2}$ and $Q^{2}$ are not too small and $z$ is not too large. Adopting the selection criteria, $Q^{2}>12\mathrm{~{}GeV}^{2}$, $z<0.9$, and $p_{t}^{2}>6.4\mathrm{~{}GeV}^{2}$, we find that only one set of data is availabe, which is presented in Reference Adloff et al. (2002). Therein, the differential cross sections with respect to $p_{t}^{2}$, $p_{t}^{\star 2}$, $Q^{2}$, $W$, $y_{\psi}^{\star}$ and $z$ are measured. To be consistent with the HERA convention, the forward direction in the $\gamma^{\star}p$ rest frame are defined as that of the incident virtual photon.

Due to the $C$-parity and color conservation, $c\bar{c}(^{3}P_{J}^{[1]})$ cannot be produced at $\mathcal{O}(\alpha^{2}\alpha_{s}^{2})$. In our calculation, we completely omit the feed down contributions from $\chi_{c}$. The contributions to the $J/\psi$ production from the $\psi(2S)$ feed down has been estimated in Reference Adloff et al. (2002). The diffractive $\psi(2S)$ are produced basically in large $z$ regions. In the region $0.75<z<0.9$, about $6\%-20\%$ of the $J/\psi$ events come from this resource. The cross sections of the inelastic $\psi(2S)$ production have the same behavior with that of the $J/\psi$ production. Throughout the kinematic region we study, the contributions from the inelastic $\psi(2S)$ feed down are estimated to be about $15\%$. The $J/\psi$ production from $b$ decay is expected in low $z$ region, and is estimated to be about $17\%$ in the region $0.3<z<0.45$ Adloff et al. (2002). All these contributions are not included in our calculation. We will address them when comparing our results to data.

To present the numerical results, we adopt the following parameter choices. The $c$-quark mass $m_{c}=1.5\mathrm{~{}GeV}$, and the electromagnetic coupling constant $\alpha=1/137$. The energy of the incident lepton and proton in the laboratory frame are $E_{l}=27.5\mathrm{~{}GeV}$ and $E_{p}=920\mathrm{~{}GeV}$, respectively. For the calculation of the $J/\psi$ production in DIS at QCD LO, the one-loop $\alpha_{s}$ running and the parton distribution function (PDF) CTEQ6L1 Pumplin et al. (2002) are employed, while for the calculation up to QCD NLO, the two-loop $\alpha_{s}$ running and the PDF CTEQ6M Pumplin et al. (2002) are employed. The value of $\alpha_{s}$ at the $Z_{0}$ boson mass is set to be $\alpha_{s}(M_{Z})=0.13$. The renormalization and factorization scales are set to be $\mu_{r}=\mu_{f}=\mu_{0}\equiv\sqrt{M^{2}+Q^{2}}$ as our default choice. The $J/\psi$ wave function at the origin is determined in terms of $|R(0)|^{2}=0.81\mathrm{~{}GeV}^{3}$ Eichten and Quigg (1995), and the CS LDME can be evaluated as $\langle\mathcal{O}\rangle=9|R(0)|^{2}/(2\pi)\approx 1.16\mathrm{~{}GeV}^{3}$. We notice that there are many parallel extractions of the CS LDME with different strategies and data. To be consistent with our calculation, we list here the value obtained through $J/\psi\rightarrow e^{+}e^{-}$ at QCD NLO in Reference Braaten and Lee (2003), i.e. $\langle\mathcal{O}\rangle=(1.005\pm 0.072)\mathrm{~{}GeV}^{3}$, that extracted from the $\eta_{c}$ hadroproduction data, $\langle\mathcal{O}\rangle=(0.645\pm 0.405)\mathrm{~{}GeV}^{3}$ Zhang et al. (2014), and that obtained from the potential model, $\langle\mathcal{O}\rangle=1.32\mathrm{~{}GeV}^{3}$ Eichten and Quigg (1995). We will find that the uncertainties in the LDMEs do not affect our phenomenological conclusions.

Our numerical results are presented in Figure 2, including the differential cross sections with respect to $p_{t}^{2}$, $p_{t}^{\star 2}$, $Q^{2}$, $W$, $y_{\psi}^{\star}$ and $z$. In the kinematic regions we study, the $K$-factor ranges from $0.85$ to $2.38$. As $p_{t}$ increases, the NLO corrections become more significant. In the largest $p_{t}^{2}$ bin, $40\mathrm{~{}GeV}^{2}<p_{t}^{2}<100\mathrm{~{}GeV}^{2}$, the $K$-factor reaches its maximum value, $2.38$. For the $Q^{2}$ bins, $12\mathrm{~{}GeV}^{2}<Q^{2}<20\mathrm{~{}GeV}^{2}$, $20\mathrm{~{}GeV}^{2}<Q^{2}<40\mathrm{~{}GeV}^{2}$ and $40\mathrm{~{}GeV}^{2}<Q^{2}<100\mathrm{~{}GeV}^{2}$, the $K$-factors are $1.25$, $1.14$ and $1.13$, respectively. As we expected, the convergence of the perturbative expansion becomes better as $Q^{2}$ increases. One can easily see from the plots that, when $p_{t}$ and $p_{t}^{\star}$ is not too large, the $K$-factor is quite close to 1, which indicates good convergence of the perturbative expansion in $\alpha_{s}$. To study the uncertainties brought in by the uncertainties of the LDMEs, we present the bands in Figure 2 covering the results for $\mu_{r}=0.5\mu_{0}$ and $\mu_{r}=2\mu_{0}$. The upper bound of the band is generally less than 1.5 times of the central curve.

Since the $\psi(2S)$ feed down contributions are not taken into account in our calculations, we need to estimate their magnitude while comparing our theoretical results to data. According to our earlier discussions, the feed down contributions are about 30% of the total $J/\psi$ production cross sections. In contrast, the ratio of the central value of the experimental data to our NLO results range from 1.4 to 4.9, among which, only 5 out of the 21 data points are smaller than 2. We can conclude that even with the feed down contributions included, the theoretical results still cannot describe the data. At least at the NLO precision, the CO mechanism is important and necessary for understanding the $J/\psi$ production in DIS.

In summary, we studied the QCD corrections to the $J/\psi$ production in DIS. This process is much more complicated than the $J/\psi$ photoproduction and hadroproduction ones due to the nonzero $Q^{2}$. In the kinematic regions that HERA experiment concerns, we found that the $K$-factors are close to 1, which indicates good convergence of the perturbative expansion. With the NLO corrections included, the CS contributions are still much smaller than the experimental data. To this end, our study iterated the importance of the color-octet mechanism.

We are indebted to Professor Geoffrey Bodwin for helpful discussions. This work is supported by National Natural Science Foundation of China (Nos. 11405268, 11647113 and 11705034).