Statistical analysis of the azimuthal asymmetry in the $J/\psi$ leptoproduction in unpolarized $ep$ collisions

In this paper, we study the azimuthal asymmetry in the $J/\psi$ leptoproduction in unpolarized $ep$ collisions. There are two independent azimuthal asymmetry modulations, namely $\mathrm{cos}(\psi)$ and $\mathrm{cos}(2\psi)$, where $\psi$ is the azimuthal angle of the lepton scattering plane with respect to the hadron-interacting plane. We calculate the two modulations as functions of four kinematic variables, and find that they provide a very good laboratory to distinguish several models describing the heavy quarkonium production, including the color-singlet (CS) model, the nonrelativistic QCD (NRQCD) associated with the $^1S_0^{[8]}$ dominance picture, and the NRQCD in which the values of all the three color-octet (CO) long-distance matrix elements are of the same order. In order to make definite conclusions, we restrict our calculation in a specific kinematic region, where the CS and CO mechanisms can be distinguished by scrutinizing the values of the $\mathrm{cos}(\psi)$ modulation, while the $^1S_0^{[8]}$ dominance picture can be tested by measuring the values of the $\mathrm{cos}(2\psi)$ modulation. Calculating their values and carrying out a meticulous statistical analysis, we find that at an integrated luminosity $\mathcal{L}=1000\mathrm{pb}^{-1}$, the statistical uncertainties of the two quantities are small enough to tell the three models apart. When this experiment is implemented at the future $ep$ colliders such as the EIC, crucial information for the $J/\psi$ production mechanism might be discovered.


Introduction
It has not been realized that the transverse polarization and transverse momentum effects in nucleons could be significant until the first measurement of inclusive pion hadroproduction was carried out at Argonne synchrotron [1][2][3], the results of which were later confirmed by Fermilab [4] at a slightly higher colliding energy. In order to understand the unexpected large transverse spin asymmetries observed in those experiments, various theoretical mechanisms and new structure functions were proposed, among which are there the well-known Sivers effect [5,6] along with a new function, the Sivers function, describing the azimuthal asymmetry of unpolarized quarks in polarized hadrons, and Collins effect [7,8], considering the different fragmenting probabilities of transversely polarized partons to hadrons with transverse momentum along different directions. Since this century, the transverse spin effects have been studied experimentally in semi-inclusive deep-inelastic scattering (SIDIS). In 2004, HERMES [9] and COMPASS [10,11] published their measurements of the single-spin asymmetry in the collisions of leptons off polarized proton and deuteron targets. One of the main advantages of SIDIS is that the transverse spin and transverse momentum effects are not mixed, as in hadroproduction; they result in different azimuthal asymmetry modulations. Another interesting feature of SIDIS is that, even with unpolarized targets, there also exist two types of azimuthal asymmetry modulations, namely the cos(ψ) and cos(2ψ) modulations, where ψ is the azimuthal angle for the observed final-state hadron with respect to the virtual-photon-target beam. A few years ago, HERMES [12] and COMPASS [13] collaborations presented the data for the azimuthal distributions of hadrons produced in deep inelastic scattering off unpolarized targets and found nonzero cos(ψ) and cos(2ψ) azimuthal asymmetry modulations. These modulations arise as long as the momenta of the final-state hadrons have transverse components, which may originate from either the intrinsic transverse motion of partons inside the targets, or the hard emission of the final-state hadrons. The former one was first studied by Cahn [14,15] and thus is called the Cahn effect. Cahn effect dominates in low p t region, while in high p t region, the real emission plays a more important role.
The rest of this paper is organized as follows. The analytical framework of our calculation is discussed in Section 2. In Section 3, we present our numerical results, according to which, the statistical analysis is given in the same section. Section 4 is a concluding remark.
2 The azimuthal asymmetry in the J/ψ leptoproduction

The calculation of the cross sections
In electron-proton collisions, which, in general, can be reduced into the virtual-photon-proton fusion process as Here, X represents the proton remnant. k, P , q, k ′ , p, and p X are the momenta of the initial-state electron, proton, and virtual photon, the final-state electron, the J/ψ, and the proton remnant. Usually, we use the following invariants to describe the kinematics of this process,  Figure 1: Illustration of the process.
the first two of which can be replaced by the following two dimensionless invariants, as long as the invariant colliding energy squared, S = (P + k) 2 , is given. Here, x is the well-known Bjorken-x, and z is named the elasticity coefficient.
To study the azimuthal effects, we adopt such frames in which the virtual photon and the proton beams are along the z and anti-z directions, respectively, as illustrated in Figure 1. The azimuthal angle of the final-state J/ψ is then defined with respect to the virtual-photon-proton beams. In such frames, the physical quantities are always labelled by a superscript ⋆, in order to distinguish from those in laboratory frames. Since the laboratory frame is not concerned in this paper, we just omit the superscript ⋆. In addition to the invariants, x, y, and z (or equivalently, Q 2 , W 2 , and z), we need two more quantities, the transverse momentum of the final-state J/ψ, p t , and its azimuthal angle, ψ, to sufficiently describe the kinematics of the J/ψ leptoproduction process. We define two dimensionless quantities, where M is the J/ψ mass, so that the SIDIS process can be described by six dimensionless variables, x, y, z, ξ, η, and ψ. Having all these kinematic variables, we can write down the hadron production cross section in unpolarized SIDIS as (see e.g. Reference [72]) where α is the electromagnetic fine-structure constant, W µν h is the hadronic tensor integrated over the phase-space of the final-state hadrons other than the J/ψ, and the normalized leptonic tensor, l µν , is defined as If only one final-state hadron is observed, the leptonic tensor, l µν , can be decomposed into the linear combination of four simpler tensors, factorizing all the azimuthally dependent terms into the corresponding coefficients as and Here, ρ is defined as where the proton mass has been neglected. Substituting Equation (2.9) into Equation (2.8) and taking into account the current-conserving equation, we can rewrite the normalized leptonic tensor in a form that is more convenient for computation as where (2.14) If we define these C i 's can be expressed in a more compact form as In collinear factorization, the protons interact with the virtual photons via the on-shell partons, and thus, the hadronic tensor, W µν h , can be further factorized as the summation of the parton-level hadronic tensors convoluted with the parton distribution functions (PDFs). For cc( 3 S 1 ) production, there is only one parton-level process at QCD LO, i.e.
while for the production of the three color-octet states, namely cc( 1 S 1 ), and cc( 3 P [8] J ), there are two parton-level processes at QCD LO, they are (2.18) Here, n represents 1 S J . Denoting the squared amplitudes of the above processes as W µν i+γ ⋆ →cc(n)+i , where µ and ν are the Lorentz indices for the virtual photons in the amplitude and its complex conjugate, respectively, W µν h can be written as where n runs over the four intermediate cc states, N i is the synthesized initial-state averaging factor, X is the fraction of the proton momentum taken by the parton, f i/p is the corresponding PDF, O J/ψ (n) is the LDME for an intermediate cc state, cc(n), producing a J/ψ, and dΦ X is defined as (2.20) Here, p i and p ′ i are the momenta of the initial-and final-state parton, respectively, and p i = Xp has been implemented. Having dX and dΦ X integrated over, W µν h can be further simplified as (see e.g. Reference [73]) Note that in Equation (2.22) the value of X has been fixed by the 0-th dimension of the δ-function at .
To write the cross section in a more compact form, we define The cross section can correspondingly be expressed as Since W µν i+γ ⋆ →cc(n)+i can be easily evaluated via Feynman diagrams, we have been equipped with all the elements needed in Equation (2.22) for calculating the J/ψ leptoproduction cross sections.

The calculation of the azimuthal asymmetry modulations
Since W µν i+γ ⋆ →cc(n)+i and the tensor basis in Equation (2.13) are all restricted in the hadronic plane, i.e. they are independent of the azimuthal angle ψ, the azimuthal asymmetry modulations appear only in the coefficients in the leptonic tensor expansions, namely in C i 's. In any case, the cross section can be written in the following form, (2.26) Therein, A cosψ and A cos2ψ are two independent azimuthal asymmetry modulations. They are functions of x, y, z, and ξ. Explicitly, they can be accessed via the following equations, (2.27)  J states, respectively.
Since the four intermediate cc states contribute independently, one can study the azimuthal asymmetry originated from them respectively. In this sense, the cross section can be written in the following form as dσ = n σ n 0 1 + A n cosψ cos(ψ) + A n cos2ψ cos(2ψ) dxdydzdξdψ. (2.28) In the next section as we present their numerical results, we just omit the superscript n. Sometimes, it is useful to study the integrated cross sections, the azimuthal asymmetry modulations should be redefined accordingly. When we observe the azimuthal asymmetry modulations in a specific kinematic region, e.g D, where D is a 4-dimensional area in the x-y-z-ξ space, the corresponding azimuthal asymmetry modulations in the region D, A D cosψ and A D cos2ψ , are defined according to the following equation, They can also be expressed in the following explicit form as (2.30) -7 -  J states, respectively.

Numerical results and phenomenological analysis
In our numerical calculation, we dopt the following parameter choices. The charm quark mass, m c , is fixed at one half of the J/ψ mass, which is approximated to M = 3.0 GeV. The colliding energy is chosen as √ S = 318 GeV to accord with the HERA experiment. To evaluate the parton distributions in protons, we employ the GRV PDF given in Reference [74]. Therein, the factorization scale is set to be µ f = Q 2 + p 2 t = S(xy + ξ). With the above parameter choices, we present the values of A cosψ in Figure 2, 3, 4, and those of A cos2ψ in Figure 5, 6, 7, for individual cc states. The ξ dependence of the azimuthal asymmetry modulations are presented for any combination of the following parameter choices, x =0.005, 0.05, 0.5, y =0.1, 0.4, 0.7, and z=0.3, 0.6, 0.9. The value of ξ ranges from 0.00001 to 0.001, corresponding to p t ≈ 1 GeV to p t ≈ 10 GeV.
As an inspiring result, the values of the azimuthal asymmetry modulations for the four cc states are remarkably distinguished in some of these kinematic regions. To tell the CS and CO channels apart, we can measure A cosψ at x > 0.05, and z ∼ 0.9, where the values of A cosψ for all the three CO states lie around 0, while those for the CS state is as large as 1. If, as argued in some papers, the J/ψ production is dominated by the CS channel, the value of A cosψ in this region should coincide with the solid curve presented in the plot at the lower right corner in Figure 4. However, this strategy is not feasible because in large-x region, the cross sections are so small that no enough events can be produced in this region to perform reasonable analysis. Fortunately, this problem can be amended by including the  small-x region, where although the CS and 3 S [8] 1 channels cannot be distinguished, the CS and CO mechanisms are still distinguishable as the 3 S [8] 1 channel is greatly suppressed in J/ψ leptoproduction, For this reason, we perform our analysis in the region x > 0.001 and 0.75 < z < 0.9, in which the value of Q 2 is generally larger than 4 GeV 2 and the validity of the perturbative expansion is guaranteed.
Another interesting feature of our numerical results is that the 1 S [8] 0 channel can also be well distinguished in some kinematic regions. For x ∼ 0.005, y ∼ 0.1, and z ∼ 0.9, the value of A cos2ψ for cc( 1 S [8] 0 ) is almost -1, while that for all the other three states is almost 0. Since most of the extractions of the CO LDMEs lead to the same picture, i.e. the 1 S [8] 0 LDME is at least one order of magnitude greater than the other two, which however is not consistent with the NRQCD scaling, this feature provides a perfect laboratory to test the 1 S [8] 0 dominance picture. If the J/ψ leptoproduction is also dominated by the 1 S [8] 0 channel, we should observe the value of A cos2ψ at almost -1 in this kinematic region. Since at larger value of x, the cross sections are very small, we can carry out the analysis in the same kinematic region as in the above paragraph, namely x > 0.001 and 0.75 < z < 0.9.
In order to draw up a practical strategy for the experimental measurement, we need to calculate the number of events assuming a luminosity, and analyse the systematic uncertainties. For this purpose, we need to specify the value of the renormalization scale µ r , the strong coupling constant α s , and electromagnetic fine structure constant α. The value of µ r in our calculation is given by µ r = S(xy + ξ). The running of α s follows the following where at n f = 5, Λ QCD is given by Λ QCD = 0.226 GeV. It is easy to verify that at the Z 0 boson mass M Z ≈ 91 GeV, the value of α s is approximately 0.130. The electromagnetic fine structure constant α = 1/137 is adopted. In the following, we carry out our study in the kinematic region, say x > 0.001 and 0.75 < z < 0.9. Due to the limit of the capability of the detectors, we further constrain the value of y and ξ in the region 0.04 < y < 0.6 and ξ > 0.00001. In the following, we denote this kinematic region as D. This configuration correspond to Q 2 > 4 GeV 2 , 60 GeV < W < 240 GeV, and p t > 1 GeV. The lower bound of X can be easily obtained as X min ≈ 0.0013, which is a moderate value to neglect the gluon saturation effects.
Integrating over x, y, z, and ξ in the region D, we obtain the cross sections contributed from the four channels as  where and J states, respectively.
As expected, the CS and CO channels can be well separated by measuring A D cosψ , while the 1 S [8] 0 channel can be distinguished by studying A D cos2ψ . Although the short-distance coefficient for the 3 S 1 channel is almost twice of that for the CS one, the cross section for this channel, however, is negligible due to a much smaller LDME which is suppressed by two orders of magnitude relative to the CS one. On the experiment side, the J/ψ events are collected in bins of ψ. The azimuthal modulations thus can be extracted by fitting the data to Equation (2.29), where linear regression is always used as a standard technique. Since the quantities that we are interested in are ratios, most of the systematic uncertainties are cancelled when doing the data analysis. As a result, we consider only the statistical uncertainty in the following. In our analysis, the range, [0, 2π), of ψ is divided into 12 equidistant bins, namely [0, π/6), . . ., [11π/6, 2π). The integral over ψ of the three modulations, 1, cos(ψ), and cos(2ψ), in each bin are calculated, and make up three vectors each of which consists of 12 elements. Explicitly, they are It is easy to verify that the three vectors are linearly independent and not correlated. If the J/ψ leptoproduction is dominated by the CS mechanism, the integrated cross section in the region D can be expressed as where is employed [75]. Assuming the integrated luminosity at the future ep collider is L = 10 3 pb −1 , we can construct the similar vector, as in the above paragraph, for this distribution as If the CO parts are not negligible, we need to sum over their contributions. Employing the LDMEs obtained in References [55,60], which are given below, This equation leads to A cosψ = −0.246. It is easy to see that this number is well separated from that led to by the CS model in the sense of the uncertainty given in Equation (3.10).
In the rest of this paper, we will demonstrate that the azimuthal asymmetry in the J/ψ leptoproduction can also distinguish the LDMEs that are consistent with the 1 S Since with the LDMEs given in Equation (3.11), the value of A cos2ψ is -0.0155, it is obvious that one can distinguish the two sets of LDMEs by scrutinize the value of A cos2ψ .

Summary
In this paper, we calculated the azimuthal asymmetry modulations in the J/ψ leptoproduction, namely A cosψ and A cos2ψ defined in Equation (2.27), as functions of x, y, z, and ξ. By scrutinizing their behaviours, we find that the CS and CO mechanisms can be well distinguished through the measurement of the values of A cosψ . Further, the 1 S [8] 0 dominance picture can also be tested by measuring the values of A cos2ψ . Restricted in the region, 0.001 < x < 1, 0.04 < y < 0.6, 0.75 < z < 0.9, and p t > 1 GeV, we carried out the calculation of the differential cross sections with respect to the azimuthal angle for three models, and found that they lead to clearly different values of the azimuthalasymmetry modulations. Having applied rigorous statistical analysis, we found that at the integrated luminosity L = 1000pb −1 , the statistical uncertainties of A cosψ and A cos2ψ are small enough to tell the three models apart. As a conclusion, the azimuthal asymmetry in the J/ψ leptoproduction provides a good laboratory for the study of the quarkonium production mechanisms. We suggest that this experiment be implemented at the future ep colliders such as the EIC.