Accessing Linearly Polarized Gluon Distribution in $J/\psi$ Production at the Electron-Ion Collider

We calculate the $cos~2 \phi$ asymmetry in $J/\psi$ production in electron-proton collision for the kinematics of the planned electron-ion collider (EIC). This directly probes the Weisz\"acker-Williams (WW) type linearly polarized gluon distribution. Assuming generalized factorization, we calculate the asymmetry at next-to-leading-order (NLO) when the energy fraction of the $J/\psi$ satisfies $z<1$ and the dominating subprocess is $\gamma^* +g \rightarrow c + {\bar c}+g$. We use non-relativistic QCD based color singlet model for $J/\psi$ production. We investigate the small $x$ region which will be accessible at the EIC. We present the upper bound of the asymmetry, as well as estimate it using a (i) Gaussian type parametrization for the TMDs and (ii) McLerran-Venugopalan model at small $x$. We find small but sizable asymmetry in all the three cases.


I. INTRODUCTION
J/ψ production is a direct probe of gluon transverse momentum dependent parton distributions (TMDs), as the leading process is the virtual photon-gluon fusion. Very little is known so far about the gluons TMDs [1], apart from a positivity bound. Recently, unpolarized TMD gluon pdfs have been extracted from LHCb data [2]. TMD pdfs are process dependent due to the initial and/or final state interactions, in other words due to the presence of gauge links in their operator definitions. Each gluon TMD contains two gauge links in contrast to the quark TMDs that contain one. Because of this the process dependence of gluon TMDs is more involved than quark TMDs [3]. The simplest possible configurations are both future pointing [++] or one future and one past pointing [+-] gauge links. In the literature related to small-x physics, the former is called Weizsäcker-Williams (WW) gluon distribution [4,5]. For unpolarized gluons, WW gluon distribution can be interpreted as the number density of gluons inside hadrons in light-cone gauge. The other distribution, [+-], is called dipole distribution [6]. This appears in many physical processes and is the Fourier transform of the color dipole amplitude [7,8].
In small x physics, these two types of unintegrated gluon distributions have been discussed in the literature quite extensively [9][10][11][12][13][14][15][16][17]. Apart from the unpolarized gluon TMD, the linearly polarized gluon TMD recently has attracted quite a lot of interest. This basically measures an interference between an amplitude when the active gluon is polarized along x (or y) direction and a complex conjugate amplitude with the gluon polarized in y (or x) direction in an unpolarized hadron [6]. This was introduced for the first time in [6] and calculated in a model in [18].
It has been shown that the linearly polarized gluon distribution affect the unpolarized cross section of scattering processes, as well an azimuthal asymmetry of the type cos2φ [19] . The linearly polarized gluon distribution is a time-reversal even (T even) object, and can be WW type or dipole type, depending on the gauge links. J/ψ production in ep collision probes the WW type linearly polarized gluon TMD through a virtual photon-gluon fusion process. The leading order (LO) process γ * + g → cc → J/ψ contributes to the asymmetry at z = 1 [20]. The linearly polarized gluon distribution has not been extracted from data yet. However there are quite a large amount of theoretical studies about how to probe it in different experiments. In [21] the authors have proposed to probe it in dijet imbalance in unpolarized hadronic collision and also in heavy quark pair production in ep collision and in pp collision [15,19,22,23]. It can also be probed in quarkonium pair production in pp collision [2], and in associated production of a dileplon and J/ψ [24]. Very recently, in [25] the authors have investigated the possibility to probe it in dijet imbalance in eA collision. h ⊥g 1 affects the transverse momentum distribution of final state hadron like Higgs boson [26][27][28][29] and heavy quarkonium [30][31][32] in unpolarized pp collision. Although h ⊥g 1 can be probed in pp and pA collision, initial and final state interactions may affect the factorization in such processes. Such complications are less in ep collision processes for example at the electron-ion collider (EIC). In a previous work [20] we have investigated the possibility of probing h ⊥g 1 in cos 2φ asymmetry in J/ψ production through the leading order (LO) process γ * + g → J/ψ at the future EIC. This 2 → 1 process contributes at z = 1, where z is the energy fraction of the photon carried by the J/ψ in the proton rest frame. Here, we extend our analysis to the kinematical region z < 1. We consider the unpolarized eP collision. The production mechanism of J/ψ is not yet well-understood theoreically. The most widely used approach is based on non-relativistic QCD (NRQCD) [33]. Here one assumes a factorization of the amplitude into a hard part where the cc pair is produced perturbatively in the process γ * + g → c +c + g. The heavy quark pair then hadronizes to form the J/ψ bound state. The hadronization is described in terms of the long distance matrix elements (LDMEs) which are obtained by fitting the data. For some LDMEs lattice calculations are available. They have definite scaling properties with respect to the velocity parameter v, which is assumed to be small. The cross section for the production of J/ψ is expressed as a double expansion in terms of the strong coupling constant α s as well as v [34]. For J/ψ, v ≈ 0.3.
In NRQCD, the heavy quark pair can be produced both in color singlet (CS) state [35][36][37][38] or in color octet (CO) state [39][40][41]. The former is called CS model and the latter, CO model. In the CS model, the heavy quark pair is produced in the hard process as a color singlet with the same quantum number as J/ψ. In [42] the J/ψ production rate for unpolarized pp collision at RHIC assuming a generalized TMD factorization was calculated in CS model, and it was found that the theoretical estimate reasonably explains the data for low values of p T , where p T is the transverse momentum of J/ψ. However, high p T spectra for J/ψ production needs the inclusion of CO states. As we showed in [43] both CS and CO contributions are needed to match the HERA data. However, in this work, as a first study, we calculate the cos 2φ asymmetry in J/ψ production in ep collision in CS model. All previous studies of this asymmetry in eP collision have considered the LO process. In this work, for the first time we investigate the asymmetry in the kinematical region z < 1, where the NLO process γ * + g → J/ψ + g is dominant. We restrict ourselves to the small x region where the WW gluon TMDs are probed.
In order to estimate the cos 2φ asymmetry we use three different models for the TMDs. First, we use a Gaussian parametrization [30][31][32] for both the linearly polarized gluon distribution and the unpolarized TMD. The linearly polarized gluons satisfy an upper bound and the asymmetry reaches its maximum value when this upper bound is saturated. We also calculate the upper bound of the asymmetry. Finally, in small x region, the WW type gluon distributions are calculated using a saturation model [44][45][46]. TMDs in McLerran-Venogopalan (MV) model, although expected to work better for a large nucleus, has been found to be phenomenologically successful for the nucleon [47]. We have used a regulated MV model in small x region for the WW type gluon TMDs. We have compared the asymmetry in all three cases in the kinematics of the planned electron-ion collider (EIC).

II. FRAMEWORK FOR CALCULATION
The process we have considered here is e(l) + p(P ) → e(l ) + J/ψ(P h ) + X Both the scattering electron and target proton are unpolarized. Four momentum of particles is represented within the round brackets. The dominating subprocess at NLO for quarkonium production in ep collision is photon-gluon fusion process, at leading order this process contributes at z = 1 [20]. In this work, we consider the NLO process γ * (q)+g(k) → J/ψ(P h )+g(p g ) and the kinematcal region z < 1, which will be accessible at EIC. The final state gluon is not detected. Here the variable z is defined as z = P · P h /P · q which is the energy fraction of J/ψ in the proton rest frame. We use a generalized factorization scheme taking into account the partonic transverse momenta. We consider the frame in which the virtual photon and proton are moving in +z and −z direction respectively. The incoming and outgoing electron form a lepton plane, which provides a reference for measuring azimuthal angles of other particles. The four momenta of proton and virtual photon q = l − l are given by [20]: where Q 2 = −q 2 and Bjorken variable, x B = Q 2 2P ·q . M p is the mass of proton. All four momenta are written in terms of light like vectors n − = P and n + = n = (q + x B P )/P · q, such that n + · n − = 1 and n 2 − = n 2 + = 0. The leptonic momenta can be written as here, s = (l+P ) 2 = 2P ·l, is the center of mass energy of electron-proton scattering. y = P ·q/P ·l, such that the relation Q 2 = sx B y hold. The virtual photon and target proton system invariant mass squared is defined as W 2 = (P + q) 2 . In terms of the light-like vectors defined above, the four momenta of initial state gluon is given as where, x = k · n is the light-cone momentum fraction. The four momentum of the final state J/ψ and the final state gluon are give by M is the mass of J/ψ. For the partonic level process: γ * (q)+g(k) → J/ψ(P h )+g(p g ), we can define the Mandelstam variables as followsŝ The φ and φ h are the azimuthal angles of the initial gluon and J/ψ transverse momentum vector respectively. Now, within the generalized TMD factorization framework, the differential cross section for the unpolarized process is given by [20] ; where L µµ is leptonic tensor which is given by with e is the electric charge of electron.
Φ νν is gluon correlator which can be parametrized in terms of gluon TMDs. For unpolarized proton, at leading twist, gluon correlator can be given as [1]: where f g 1 (x, k 2 ⊥ ) is the unpolarized gluon distribution and h ⊥g 1 (x, k 2 ⊥ ) is the linearly polarized gluon distribution.

J/ψ PRODUCTION IN NRQCD BASED COLOR SINGLET (CS) FRAMEWORK
The dominating subprocess at NLO is γ * + g → J/ψ + g. All the tree level Feynman diagram corresponding to this process are given in Fig. 1 .
The general expression of the amplitude for the bound state production of J/ψ in NRQCD framework can be written as [20,30] : As we have imposed a cutoff on z, z < 0.9, we do not need to consider the virtual diagrams as they contribute at z = 1. In the above equation, 2k is the relative momentum of heavy quarks and O(q, k, P h , k ) is calculated from the Feynman diagrams. The polarization vector of heavy quark, anti-quark legs are absorbed into the bound state wave function. By considering contribution from all the Feynman diagrams, O(q, k, P h , k ) is given by Where, O m , (m = 1, 2, ...6) are corresponding to each Feynman diagrams and C m represents the color factor of corresponding diagram. The expressions for O m are given below Here, the mass of bound state M is assumed to be twice the mass of charm quark (m c ) i.e. M = 2m c , .The charge conjugation invariance allow us to write the expressions for (O 4 , O 5 and O 6 ), from the other Feynman diagrams, by reversing the fermion line and replacing k by −k . Assuming the QQ is formed in color singlet state, the color factor of each diagram is given by The SU (3) Clebsch-Gordan coefficients for CS are given by where N c is the number of colors. The genertors of the SU(3) group satisfies the relations: . From these relations, we get the color factor for the production of QQ pair in CS state as follows; The spin projection operator, given in the equation of amplitude of the bound state, includes the spinors of heavy quark and anti-quark and is given by [30]: where Π SSz = γ 5 for singlet (S = 0) state and Π SSz = / ε sz (P h ) for triplet (S = 1) state. ε sz (P h ) is the spin polarization vector of QQ pair. Since, k << P h , one can perform Taylor expansion of the amplitude around k = 0. In that expansion, the first term gives the S-waves(L=0,J=0,1).
For the P-waves(l=1,J=0,1,2), we need to consider the linear terms in k in the expansion as the radial wavefunction R 1 (0) = 0 for P −wave. Since, J/ψ is a 3 S 1 state, in the color singlet model we calculate contribution of the CS state 3 S 1 . where, We have the following symmetry relations for 3 S 1 state The final expression for the amplitude for 3 S 1 state is given by

CALCULATION OF THE ASYMMETRY
As we are using a framework based on generalized TMD factorization, we consider a kinematical region in which the transverse momentum of J/ψ is small compared to the mass of J/ψ, M i.e. P h⊥ < M . The final gluon carry the momenta fraction (1 − z), as z = P ·P h P ·q , is the energy fraction transferred from the photon to J/ψ in proton rest frame. So, this means that when z → 1, the outgoing gluon is soft which leads to infrared divergence. Hence, we consider z < 0.9 to keep the final gluon hard and avoid such divergence. Moreover, gluon and heavy quark fragmentation can also contribute to quarkonium production significantly for P h⊥ > 4 GeV . We have imposed an upper limit on P h⊥ . In order to eliminate the fragmentation of the hard gluon into J/ψ we also use a lower bound on z, namely 0.3 < z.
In the differential cross section, there is a contraction of four tensors which is written as where the individual components are defined above. The summation over the transverse polarization of the final on-shell gluon is given by with n µ = P µ h M . We have three amplitudes and their corresponding conjugates that will contribute. We use the notation where i = 1, 2, 3, corresponds to the contribution from each diagram. So, the cross section will get contribution from nine terms (of the form M i M j , where i, j = 1, 2, 3) and hence, the differential cross section can be written as dσ = 1 2s where M = i M i . Out of the nine terms in | M | 2 , six are interference terms with a symmetry So, effectively we need to compute six terms.
In a frame where the virtual photon and target proton move along the z-axis, and the lepton scattering plane defines the azimuthal angles φ l = φ l = 0, then we have and the delta function can be expressed as here z 1 = z and the delta function sets z 2 = (1 − z). Hence, after integrating with respect to x, z 2 and p g⊥ , the final form of the differential cross section can be given by As we are interested in the small x region, we neglect terms containing higher powers of x B ; also as z < 1, we neglect terms containing higher powers of z. We also keep terms only upto ( . The leading terms in the numerator of the cos(2φ) asymmetry comes only from the first Feynman diagram. These terms are given in the appendix. The denominator of the cos(2φ) asymmetry, which is defined below, is simply the cross section integrated over the azimuthal angle φ h . The leading terms in the cross section comes from f g 1 term. All the terms corresponding to h ⊥g 1 are suppressed by k 2 ⊥ /M 2 p . Hence, from the approximations we mentioned above, the leading terms in the cross section in the denominator of cos(2φ) asymmetry are coming from M 1 M 1 , M 1 M 2 , M 1 M 3 . These contributions are given in the appendix.
The differential cross section then can be given as The coefficients A 0 , A 1 , B 0 , B 1 and B 2 are given in the appendix. The cos(2φ) asymmetry is defined as In order to estimate the cos(2φ) asymmetry, we need to parametrize the TMDs. We discuss two parametrization models, Gaussian parameterization of the TMDs and McLerran-Venugopalan(MV) model. We also calculate the upper bound of the asymmetry.

A. Gaussian parametrization of the TMDs
Both for the linearly polarized gluon distribution and the unpolarized gluon TMD, a Gaussian parametrization is used widely in the literature. The linearly polarized gluon distribution satisfies the model independent positivity bound [1]; The Gaussian parametrizations satisfy the positivity bound but does not saturate it. They are as follows [30][31][32]; where, r(0 < r < 1) is a parameter, in our numerical estimates we took r = 1/3. f g 1 (x, Q 2 ) is the collinear gluon pdf. The width of the Gaussian, k 2 ⊥ , depends on the energy scale of the process. Following [30], we took k 2 ⊥ = 0.25 GeV 2 .

B. Upper bound of the asymmetry
The asymmetry reaches its maximum value when the positivity bound given by Eq. (38) is saturated. Using this, we calculate the upper bound of | cos(2φ h ) | as below [19]; where S ⊥ is transverse size of the nucleus or nucleon. Q sg is the saturation scale, which in MV model, is defined as Q 2 sg = α s N c µ A ln 1 and µ A S ⊥ = α s 2πA, where A = 1 for the proton. Following the approach of [47], we have used a regulated version of the MV model in our calculation of the asymmetry. The ratio of linearly polarized and unpolarized distribution in MV model can be given by For Q 2 sg0 = (N c /C F ) × Q 2 s0 , where Q 2 s0 = 0.35 GeV 2 at x = 0.01 and Λ QCD = 0.2 GeV, the ratio is below 1 for all k ⊥ . Below we give our numerical results. We have estimated the cos(2φ h ) asymmetry in J/ψ production in the kinematics of EIC.
MSTW2008 [48] is used for collinear PDFs. As stated in the introduction, we have used cuts on z, 0.3 < z < 0.9. As we know, gluon initiated processes are enhanced at small x. In fact small x values will be accessed at EIC, and this kinematical region will be very important in determining the gluon TMDs including the linearly polarized gluon TMD. In this work, we have studied the cos 2φ asymmetry for EIC in the small x region. It is to be noted that x is related to the Bjorken variable x B through Eq. 9. Smaller x values also restrict Q 2 to be small, in this work we took Q 2 to be of the same order and bounded by M 2 , which is the mass of J/ψ. For both the parametrizations used, the asymmetry is negative, which is consistent with the LO calculation [20]. In the plots we show the magnitude of the asymmetry. We have taken the CS LDMEs from [49]. As only one state contributes in CS model, namely 3 S 1 the asymmetry does not depend on the specific set of LDME. This is different from CO model, where even at LO, contribution comes from several states [20], and the result depends on the choice of LDMEs.

IV. CONCLUSION
In this work, we have calculated the cos 2φ asymmetry in electroproduction of J/ψ at EIC, that probes the linearly polarized gluon distribution in the unpolarized proton. We calculated the asymmetry in the kinematical region z < 1, where the NLO subprocess γ * + g → J/ψ + g gives the leading contribution. The gluon TMDs probed in this process are of Weizsäcker-Williams (WW) type. As gluon distributions pay an important role in small x region, we investigate the asymmetry in the small x kinematical region, using a Gaussian parametrization of the TMDs as well as in McLerran-Venugopalan model. We also show the upper bound of the asymmetry saturating the inequality for the linearly polarized gluon distribution. At EIC, low values of x also restricts the Q 2 (virtuality of the photon) values. We have calculated the J/ψ production amplitude in NRQCD based color singlet (CS) approach. The asymmetry in the kinematical region considered is small but sizable. The magnitude of the asymmetry may depend on the production mechanism of the quarkonium. As shown in [20], CS mechanism underestimates the J/ψ production at HERA, and both CS and CO contributions are needed to describe the data. In CO formalism contribution will come from several LDMEs in the final state, which may enhance the asymmetry. We plan to see the effect of the CO mechanism on the asymmetry in a future work. Another interesting study would be the effect of small-x evolution on the asymmetry. In any case, the cos 2φ asymmetry in J/ψ production at the EIC will be an important tool to gain information on the WW type linearly polarized gluon distribution.

V. ACKNOWLEDGEMENT
We thank P. J. Mulders, E. Petreska and D. Boer for discussion. Part of the work of AM was done at NIKHEF, Amsterdam, the visit was supported by the European Research Council under the "Ideas" program QWORK (contract 320389).