$Cos(2\phi_h)$ asymmetry in $J/\psi$ production in unpolarized $ep$ collision

We present a calculation of the $cos (2 \phi_h)$ asymmetry in $J/\psi$ production in electron-proton collision at the future electron-ion collider (EIC), a useful channel to probe the gluon TMDs. We calculate the asymmetry at next-to-leading order (NLO) in $\alpha_s$ in the framework of generalized factorization. The dominating sub-process is $\gamma^* +g \rightarrow J/\psi+g$. The production of $J/\psi$ is calculated in the non-relativistic QCD(NRQCD) framework with the inclusion of both color singlet and color octet contributions. Numerical estimates of the $cos(2\phi_h)$ asymmetry are given in the kinematical region to be accessed by the future EIC. The asymmetry depends on the parameterization of the gluon TMDs, as well as on the long distance matrix elements (LDMEs). We use both Gaussian-type parameterization and McLerran-Venugopalan model for the TMDs in the kinematical region of small-$x$, where the gluons play a dominant role. We obtain sizable asymmetry, which may be useful to probe the ratio of the linearly polarized and the unpolarized gluon distribution in the proton.


I. INTRODUCTION
Quantum Chromodynamics (QCD) is an exceptionally rich and complex theory of strong interactions between quarks and gluons, the fundamental constituents of matter. However, these quarks and gluons do not exist in nature as free particles but are confined inside hadrons, and their fundamental properties can be explored only with the help of scattering processes. Hadron physics community has expanded its inquiry beyond the ordinary one dimensional collinear parton distribution functions (PDFs) in the motion of the parton and its spatial distribution in a direction perpendicular to the momentum of the parent hadron.
To account for transverse motion, the collinear PDFs were extended to transverse-momentumdependant PDFs, also referred to as TMDs [1][2][3]. TMDs have attracted an enormous amount of interest and are being investigated at the current facilities including JLAB 12 GeV upgrade, RHIC and are planned to be investigated at the future electron-ion collider (EIC). They are considered as an extension of the standard, one-dimensional, parton distribution functions (PDFs) to the three-dimensional momentum space. Due to the gauge invariance, the operator definition of TMDs requires the inclusion of gauge links or Wilson lines. Unlike colinear PDFs, which are universal, TMDs are process dependent due to their initial and final state interactions [4,5], or the gauge links. In other words, the TMDs extracted from semi-inclusive deep inelastic scattering (SIDIS) are not the same as extracted from Drell-Yan processes due to the difference in gauge links. TMDs are sensitive to the soft gluon exchanges and the color flow in the specific sub-processes in which they are probed. SIDIS and Drell-Yan processes provide the majority of experimental data for the extraction of TMDs, where the observables of most interest are the single-spin asymmetries (SSA) and azimuthal asymmetries, which have been measured and are currently under direct experimental scrutiny [6][7][8].
A lot of work has recently been done to extract quark TMDs inside a proton from low energy data from HERMES, COMPASS or JLab experiments. Contrarily a little is experimentally known about the gluon TMDs [9], as they typically require higher-energy scattering processes and are harder to isolate in comparison to quark TMDs. Each gluon TMD contains multiple gauge links while the quark TMDs contain one, due to which the process dependence of gluon TMDs is more involved than quark TMDs [10]. Gluon TMDs parameterize the transverse mo-tion of gluons inside a proton. In the parton model, the gluon correlator of the unpolarized spin 1/2 hadron is parameterized in terms of leading twist transverse momentum dependent (TMD) distribution functions [9] i.e. f g 1 (x, k 2 ⊥ ) and h ⊥g 1 (x, k 2 ⊥ ). The function f g 1 (x, k 2 ⊥ ) represents the probability of finding an unpolarized gluon, within an unpolarized hadron, with a longitudinal momentum fraction x and transverse momentum k ⊥ , while as h ⊥g 1 (x, k 2 ⊥ ) represents the distribution of linearly polarized gluons within the unpolarized hadron and is also known as Boer-Mulders function. These are the only two TMDs of the unpolarized proton that provide essential knowledge on the transverse dynamics of the gluon content of the proton.
They are also important for the proper explanation of gluon-fusion processes at all energies. At small-x, it turns out that the linearly polarized distribution may reach its maximally allowed size, bounded by the unpolarized gluon density [9]. Depending on the gauge links, there are two types of gluon TMDs, Weizsäcker-Williams (WW) type where both gauge links are either future or past pointing [11,12]; and dipole type, where one gauge link is future pointing and one past [13]. ep collision probes the WW type gluon distribution. Various methods have been proposed to measure both f g 1 (x, k 2 ⊥ ) and h ⊥g 1 (x, k 2 ⊥ ), as well as other gluon TMDs, for example [10,.
In the last few years, the distribution of linearly polarised gluons within an unpolarized proton has attracted a lot of interest. Yet they have not been extracted experimentally. A lot of theoretical investigations has been put forward to probe h ⊥g 1 (x, k 2 ⊥ ), and a model-independent theoretical upper bound can be found in Ref. [9,35]. They are time-reversal even (T even) object, affecting the unpolarized cross-section of scattering processes, as well as an azimuthal asymmetry of the type cos(2φ h ) [36]. Processes like heavy quark pair or dijet production in SIDIS [36], diphoton pair [16] and Υ(1S)+jet [37] production in pp collision have been suggested for extracting h ⊥g 1 (x, k 2 ⊥ ). It has been seen in these processes that h ⊥g 1 (x, k 2 ⊥ ) can be probed by measuring azimuthal asymmetries. It can also be probed in heavy-quark pair production in lepton-proton, proton-proton collisions [36,38,39] and associated production of dilepton and J/ψ [25] as well. Quarkonium production is a useful tool to probe the gluon TMDs (see [40]for a recent review). J/ψ production in SIDIS is a good channel, as the effect of the gauge links are simpler compared to pp collision, and one can assume the generalized factorization.
Quite a lot of theoretical work has been done recently in this direction. TMDs can be accessed through this process when the transverse momentum of the produced J/ψ (P h⊥ ) is not very large, | P h⊥ |< M , where M is the mass of J/ψ. Such kinematical region is expected to be accessible at the future electron-ion collider (EIC). In Ref. [22] the authors probed h ⊥g 1 in cos(2φ h ) asymmetry in J/ψ production through the leading order (LO) process γ * + g → J/ψ at the future EIC at z = 1, where z is the fraction of photon energy transferred to J/ψ.
The calculation of cos(2φ h ) asymmetry which directly probes the WW type linearly polarized gluon distribution at next-to-leading-order(NLO) using color singlet model (CSM) is already calculated in [30], and in J/ψ +jet production in ep collision in NRQCD based color octet model in [41]. In this work, we extend our study to investigate the cos(2φ h ) asymmetry in J/ψ production in ep collision within the NRQCD based Color Octet Model(COM) in the kinematical region z < 1.
The J/ψ meson is a charm-anti-charm quarkonium bound state with odd charge parity. It has gained a lot of interest, it not only provides knowledge on the strong interactions responsible for hadronization but also acts as a probe of the perturbative and non-perturbative aspects of QCD.
The three main models which have been developed to describe the production mechanism of J/ψ meson are, color singlet model (CSM) [42][43][44], color evaporation model (CEM) [45,46] and color octet model(COM) [47], are all successful but at different energies. Each of these models has its field of applicability, and they have also achieved considerable success in describing the experimental data, but which one is the best is still an open problem. All these models have one thing in common, the cross-section is factorized into a perturbative hard part where the quarks and gluons form the heavy quark and antiquark pair, with momentum of order M (heavy quark mass) and a non-perturbative(non-relativistic) soft part where the heavy quark pair forming a bound state at a scale of order Λ QCD . Former is also known as short-distance part and can be calculated in order α s (M ). All the information of hadronization is encoded in the long-distance matrix elements (LDME), which are usually extracted by fitting experimental data. They describe the transition probability to form the quarkonium state from the heavy quark pair. A particularly elegant approach for separating a perturbative hard part from the non-perturbative soft part is to recast the analysis in terms of NRQCD [48], an effective field theory designed precisely for this purpose.
Color Octet Model is based on NRQCD effective field theory. In COM, the heavy quark relative momentum (M Q ν) is much less than the mass (M Q ) of heavy quark, where ν is the relative velocity between the heavy quark and anti-quark in the centre-of-mass frame of the quarkonium pair. The quark potential model calculations indicate that ν 2 ≈ 0.3 for charmonium. The cross-section of quarkonium is expressed as an infinite series in the limit ν 1. Each term in the series is the product of QQ pair cross-section in a definite state n = 2S+1 L [a] J and LDMEs. Here J, L and S are total angular momentum, orbital angular momentum and spin quantum numbers, respectively, and a is the color multiplicity carrying a value of 1 for color singlet and 8 for color octet state. Both the color singlet and color octet states are included in this model for the production rate of quarkonium. A good description of J/ψ by the color octet model is given at RHIC energies [49]. Azimuthal asymmetries in J/ψ production in these processes may also be used as a tool to get information on the LDMEs in NRQCD [50]. In our study of cos(2φ h ) azimuthal asymmetry in J/ψ, we will use NRQCD formalism based on color octet model. We will be taking NLO process γ * + g → J/ψ + g into the consideration, and investigate mainly the small-x region, where the gluon TMDs play a very important role.
The rest of this paper is organized as follows. The analytical framework of our calculation is discussed in Section II. In Section III, we present our numerical results and finally, in Section IV, we conclude.

A. Calculation of the cross section and asymmetry
We consider semi-inclusive production of J/ψ in unpolarized ep collision : where the quantities within the brackets are the four-momentum of corresponding particles and X represents the proton remnant. We use the following invariants to describe the kinematics of this process, The other two dimensionless invariants are, The quantity Q 2 is the virtuality of the photon, W 2 is virtual photon-target system invariant mass squared, z is the fraction of energy transferred from photon to J/ψ in the frame where the initial proton is at rest, y is the inelasticity variable and has a physical interpretation as the fraction of the energy of the electron transferred to the proton. The variable is the known as Bjoken-x.
To study the azimuthal asymmetry, we use a frame in which the incoming proton and the virtual photon exchanged in the process move in +z and −z directions. The kinematics here is defined in terms of two light-like vectors with the help of a Sudakov decomposition, here chosen to be the momentum P (= n − ) of the incoming proton, and a second vector n(= n + ), obeying the relations n · P = 1 and n 2 + = n 2 − = 0. Thus one can have the following expressions for the momenta of the incoming proton and virtual photon(q = l − l ): Where M p is the mass of proton. Moreover, the momenta of incoming and outgoing lepton can be written as: Here s = (l + P ) 2 = 2P · l is the squared centre-of-mass energy of the electron-proton scattering and Q 2 , s, y and x B are related by the relation Q 2 = s x B y. At leading order, J/ψ is produced by the virtual-photon gluon fusion for which z = 1 [22]. Since at small-x the proton is rich in gluons, the dominant partonic sub-process at NLO for the J/ψ production is γ * (q) + g(k) → J/ψ(P h ) + g(p g ). In terms of light-like vectors defined above the four momentum of initial state gluon, final state gluon and J/ψ are expressed as: Here x = k ·n is the light cone momentum fraction of the gluon, M is the mass of J/ψ and P 2 h = −P 2 h⊥ . The Mandelstam variables for the partonic sub-process γ * (q) + g(k) → J/ψ(P h ) + g(p g ) become:ŝ here φ and φ h represent the azimuthal angles of initial gluon and the transverse momentum of J/ψ, respectively. Due to the presence of the gauge links or initial/final state interactions, factorization in many of such processes are still not proven. However, ep collision processes are less complicated than pp collisions in terms of color flow. In [51] the process ep → e + J/ψ + X is compared with the SIDIS, for which TMD factorization has been proven. It is argued that the additional scale, namely, the mass of J/ψ does not affect the gauge link structure, and TMD factorization is expected to be valid for this process as well. Following the approach of Refs. [22,52], we assume the generalized factorization in the kinematics considered. The differential scattering cross-section can be written as a convolution of leptonic tensor, a soft parton correlator for the incoming hadron and a hard part: The term M µν represents the amplitude of J/ψ production in the γ * + g → J/ψ + g partonic sub-process. The leptonic tensor is given by where e is the electronic charge. At leading twist, the gluon correlator of the unpolarized proton contains two TMD gluon distribution functions Here g νν ⊥ = g νν − P ν n ν /P · n − P ν n ν /P · n. The quantities f g 1 (x, k 2 ⊥ ) and h ⊥g 1 (x, k 2 ⊥ ) represent the unpolarized and the linearly polarized gluon distribution functions respectively.

B. J/ψ production in NRQCD based color octet model
The amplitude for the production of J/ψ can be written as follows [52,53] Here Ψ LLz (k ) is the non-relativistic bound-state wave function with orbital angular momentum L, L z and the relative momentum k of the heavy quark in the quarkonium rest frame; k is assumed to be smaller than P h . LL z ; SS z |JJ z are the usual Clebsch-Gordon coefficients, O(q, k, P h , k ) represents the amplitude for the production of the heavy quark pair QQ and is calculated from the Feynman diagrams. The Feynman diagrams relevant for this process are given in Fig. 1. The polarization vectors of the initial gluons and the heavy quark-antiquark legs are absorbed into the definitions of the gluon correlators and in the bound state wave function, respectively. The quantity plays the role of spin projection operator [52,53] with where ε Sz (P h ) is the spin vector of the QQ system. We consider here γ * + g → J/ψ + g, which is the dominating partonic sub-process. 1. Feynman diagrams for γ * + g → J/ψ + g at NLO order partonic sub-process.
By considering the contribution from all the above Feynman diagrams we can write the amplitude as where C m represents the color factor of each diagram here the summation runs over the colors of the outgoing quark and antiquark. The SU (3) Clebsch-Gordan coefficients for color singlet and color octet states, respectively, are given as and they project out the color state of the QQ pair either in the color singlet or in the color octet state. Here, N c represents the number of colors. In the fundamental representation the generators of the SU (3) group is denoted by t a following which T r(t a t b ) = δ ab /2 and T r(t a t b t c ) = 1/4(d abc + if abc ). In case of color octet state, the color factors for the production of initial QQ is The amplitudes O m (q, k, P h , k ) for the above Feynman diagrams are written as: Here M = 2m c , with m c representing the charm quark mass and the tensor T µρσ (k, p g ) = where r = (|r|, θ, ϕ) in spherical coordinates, R L (|r|) and Y LLz (θ, ϕ) is the radial wave function and spherical harmonic respectively. In particular this means, As the first term of Taylor expansion at k = 0 does not depend on k anymore, so it gives the contribution from S waves (L = 0, J = 0, 1), where O(0) = (q, k, P h , 0) and P SSz (0) = P SSz (P h , 0). For P waves (L = 1, J = 0, 1, 2), R 1 (0) = 0, so in the Taylor expansion the terms linear in k α are considered in Eq. 14. Thus, where ε α Lz (P h ) is a polarisation vector for L = 1 bound state and R 1 (0) represents the derivative of P -wave (radial) wave function evaluated at origin. Thus For P -wave we could utilize the following relations for Clebsch-Gordan coefficients and various polarisation vectors [54,55]: Here the term ε α Jz (P h ) represents the polarization vector of bound state with J = 1, obeying the following relations: while as ε αβ Jz (P h ) is the polarization tensor for J = 2 bound state and obeys the following set of relations [54,55]: The relation between the radial wave function R 0 (0) and its derivative at origin R 1 (0) with the LDME can be found in the Eqs.(36 − 38) within Ref. [27].
The numerical values of LDMEs for these color octet states are extracted from Ref. [56]. Now using the above formalism together with the symmetry relations and sum over the color factors from Ref. [27] in the Eqs. (39,40,43,44,47), we could write the amplitude for color octet states 1,2) ).

3 S 1 Amplitude
The final expression for the 3 S 1 scattering amplitude: Due to the symmetry relations [27] Feynman diagrams 4 and 8 cancel and do not contribute to

1 S 0 Amplitude
For the 1 S 0 scattering amplitude the expression reads as: where O 1 (0), O 2 (0) and O 3 (0) are given in Eq. (31) and The amplitude for 3 P J can be written as,

C. cos(2φ h ) Azimuthal Asymmetry
To probe the TMDs in this process, one needs two well-separated scales. While Q 2 gives the hard scale, the other scale is given by P h⊥ , and we consider a kinematical region where the transverse momentum of J/ψ is small compared to the mass of J/ψ, M , i.e. P h⊥ < M .
For large transverse momentum of J/ψ collinear factorization is expected to hold. The generic structure of the TMD cross-section defined by Eq. 11 is easily obtained from the contraction of four tensors which are We get a contribution for the amplitude and their corresponding complex conjugate amplitude from all the six states 1 S  ε λa µ (p g )ε * λa µ (p g ) = −g µµ + p gµ n gµ + p gµ n gµ p g · n g − p gµ p gµ (p g · n g ) 2 (36) with n µ g = P µ h M . 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, we can write: (37) and the delta function can be written as here, the delta function sets z 2 = (1 − z). After integration over x, z 2 and p g⊥ , the final form of the differential cross section can be written as where, | M | 2 = dφ| M | 2 , and k ⊥ is the magnitude of k ⊥ . We only keep the transverse momentum of initial gluon up to O(k 2 ⊥ /M 2 p ). We expand in P h⊥ /M and keep the terms up to O(P 2 h⊥ /M 2 ). As we are interested in the small-x domain, we have expanded the amplitudes in x B and did not consider higher order terms in x B . Thus we could write the final expression for differential scattering cross-section as where The analytic expressions of the coefficients A 0 , A 1 , A 2 , B 0 , B 1 and B 2 are very large, which we have not included in this article. They are available upon request. The cos(2φ h ) asymmetry measured in experiments is defined as where φ h is the azimuthal angle of J/ψ production plane with the lepton plane. The cos(2φ h ) asymmetry as a function of P h⊥ , x B , z and y can be written as: Thus when the color octet contributions are taken into account the unpolarized gluon distribution and the linearly polarized distributions are not disentangled in this process in the kinematics considered. However, the above asymmetry depends on the ratio of the two, and can be used to extract this ratio.

III. RESULTS AND DISCUSSION
In this section, we present numerical estimates of the cos(2φ h ) asymmetry in the kinematical region to be accessed at EIC. As gluon TMDs are particularly important in the small-x domain, we are considering the small-x kinematics. At z = 1, one gets contribution from the leading order as well as diffractive contributions. Here the momentum fraction of the final gluon is (1 − z). This means as z → 1, the final gluon becomes soft. To keep the final gluon hard, we have imposed a cutoff z < 0.9. Also, to avoid the gluon fragmentation contribution to J/ψ, we have used a lower cutoff z > 0.1. We have checked that changing the lower cutoff does not affect the asymmetry much. We took mass of the proton to be M p = 1 GeV.
The contraction in the calculation above for the different states i.e., 1 S J(=0,1,2) is calculated using the FeynCalc [57,58]. In all the plots of the asymmetry, the long distance matrix elements (LDMEs) are taken from Ref. [56] except for the right panel of Fig. 3 and for the plots in Fig. 4, where we have used different sets of LDMEs. The asymmetry depends on the parameterization/model used for the gluon TMDs. In our estimate, we have used two sets of parameterization to calculate the cos(2φ h ) asymmetry: 1) Gaussian-type parameterization [20,23,52] for both the unpolarized as well as the linearly polarized TMDs and 2) McLerran-Venugopalan model [59][60][61] at small-x region.

cos(2φ h ) asymmetry in the Gaussian parameterization
In Gaussian type of parameterization, both the TMDs, f g 1 and h ⊥g 1 are factorized into a product of collinear PDFs times exponential factor which is a function of the transverse momentum.
r(0 < r < 1) is a parameter. We took the Gaussian width k 2 ⊥ = 0.25 GeV 2 and r = 1/3 in all plots except Fig. 5. The term f g 1 (x, µ) is the collinear PDF, for which MSTW2008 [62] is used and is probed at the scale µ = M 2 + P 2 h⊥ , where M (= 3.096 GeV) is the mass of J/ψ. The linearly polarized gluon distribution in the parameterization above satisfies the positivity bound [9], but does not saturate it. integration ranges are 0 < P h⊥ < 2, 0.1 < z < 0.9. For both plots CMSWZ set of LDMEs [56] is used.
In Fig. 2-5, we have used the Gaussian parameterization of the TMDs. In Fig. 2 and 5 we calculated the total asymmetry by including the contribution from both the color singlet and color octet states in the framework of NRQCD.
In left panel of Fig. 2 we showed the asymmetry as a function of P h⊥ at center of mass energy √ s = 100 GeV and at a fixed value of Q 2 = 15 GeV 2 . The integration ranges are 0.0015 < x B < 0.1, 0.1 < z < 0.9 and y is set by the values of Q 2 , s and x B . In the right panel of Fig. 2, we showed the asymmetry as a function of y at √ s = 120 GeV and at fixed x B . The asymmetry is larger for smaller values of x B and approaches to zero as y tends to 1. We obtained negative asymmetry which is consistent with the LO calculations shown in [22].  [56]. Right: comparing the asymmetry for different LDMEs sets: CMSWZ [56], SV [63], ZSSL [64], BK [65].
In the left panel of Fig. 3, we present the contribution to the asymmetry from the color singlet(CS) and the color octet(CO) states. From the plot, we see that the color octet states are giving a significant contribution to the asymmetry, whereas the contribution from the CS is almost zero and slightly positive in higher P h⊥ region. In the right panel of Fig. 3, we showed the asymmetry for different sets of LDMEs. We see that the magnitude and the sign of the asymmetry depends on the set of LDMEs used. In fact, this is because of different states contributing to the asymmetry depending on the LDMEs. We have a larger asymmetry for LDMEs set BK [65]. Also, in Ref. [22] and Ref. [41], it has been shown that at LO and at NLO, respectively, in the color octet model, the asymmetry comes from the several states and the results depend on the choice of LDMEs. In the left panel of Fig. 4, we show the contribution coming from all the individual states to the asymmetry for the LDMEs set CMSWZ [56]. For this set of LDMEs, there is a dominance of one single state, 1 S 0 to the asymmetry. Whereas for the LDMEs set SV [63], which is shown in the right panel of the same figure, we have major contributions from two states 1 S In Fig. 5, we plotted the asymmetry as a function of P h⊥ for different values of the parameter r(left) and Gaussian width k 2 ⊥ (right). We do not see any significant change in the asymmetry over different values of Gaussian width and the parameter r.
Also, we have checked that there is no significant change in the asymmetry if the value of the centre of mass energy is changed for a fixed value of Q 2 ; but the asymmetry decreases with the increase of the Q 2 for a fixed value of s.

cos(2φ h ) asymmetry in the McLerran-Venugopalan (MV) model
The McLerran-Venugopalan (MV) model [59][60][61], a classical model, helps us to calculate the gluon distribution in a large nucleus at small x. In this model, the Gaussian distribution of color charges is assumed which act like a static sources, producing the soft gluons by the Yang-Mills equations. Although originally proposed for a large nucleus, this model has been found to give reasonable phenomenological results in azimuthal asymmetries for a proton in small-x domain [30,36,41]. The analytical expressions for unpolarized and linearly polarized Weizsäcker-Williams gluon distribution can be written as, following [66,67] : here S ⊥ is the transverse size of the proton, Q sg (ρ) is the saturation scale for gluons, which in general is a function of x but in MV model depends on dipole size ρ logarithmically. Here 0 state. In the same Fig. 6 (right), we show the comparison of the asymmetry in the Gaussian and MV model, respectively (with two different values of Q sg0 [68,69]) within the same kinematical region as discussed above. The asymmetry in MV model depends on the saturation scale. We note that the Gaussian parameterization of the TMDs gives larger cos(2φ h ) asymmetry. In Fig. 7 we compared the contribution to the asymmetry from color singlet and color octet states. We found that larger contribution to the cos(2φ h ) asymmetry comes from color octet states, which is negative. Color singlet states give a small positive contribution to the asymmetry.

IV. CONCLUSION
In this article, we studied the cos(2φ h ) asymmetry in J/ψ production in an electron-proton collision in the kinematics of the future electron-ion collider. We calculated them in the small-x domain, where the gluon TMD, namely the unpolarized and linearly polarized gluon TMD are important in unpolarized scattering. The dominant sub-process in this kinematical region for J/ψ production is the virtual-photon-gluon fusion process γ * + g → J/ψ + g. We used the NRQCD based color octet formalism for calculating the J/ψ production rate. We obtain a small but sizable cos(2φ h ) asymmetry. The asymmetry depends on the parameterization of the gluon TMDs used. We used both the Gaussian as well as the MV model for the parameterization. The magnitude of the asymmetry was found to be larger for Gaussian parameterization. We included contributions both from CO as well as CS states. The asymmetry depends on the LDMEs used, in particular, contributions from individual states were found to depend substantially on the set of LDMEs used. Overall, our calculation shows that the cos(2φ h ) asymmetry in J/ψ production could be a very useful tool to probe the ratio of the linearly polarized gluon TMD and the unpolarized gluon TMD in the small-x region at the EIC.