$\cos2\phi_t$ azimuthal asymmetry in back-to-back $J/\psi$-jet production in $e~p\rightarrow e~J/\psi~Jet~ X$ at the EIC

In this article, we investigate the $\cos2\phi_t$ azimuthal asymmetry in $e ~p\rightarrow e ~J/\psi ~Jet~ X$, where the $J/\psi$-jet pair is almost back-to-back in the transverse plane, within the framework of the generalized parton model(GPM). We use non-relativistic QCD(NRQCD) to calculate the $J/\psi$ production amplitude and incorporate both color singlet(CS) and color octet(CO) contributions to the asymmetry. We estimate the asymmetry using different parameterizations of the gluon TMDs in the kinematics that can be accessed at the future electron-ion collider (EIC) and also investigate the impact of transverse momentum dependent (TMD) evolution on the asymmetry. We present the contributions coming from different states to the asymmetry in NRQCD.


I. INTRODUCTION
Transverse momentum dependent parton distributions (TMDs) [1][2][3][4][5][6] give a tomographic picture of the nucleon in terms of quarks and gluons in momentum space. TMDs play an important role in processes where two scales are involved; for example, in semi-inclusive deep inelastic scattering (SIDIS) where apart from the photon virtuality, one measures the transverse momentum of the outgoing particle, or Drell-Yan process where the transverse momentum of the outgoing lepton-pair provide the second scale. For such processes, one can apply the generalized factorization involving TMDs. However, the TMD factorization has not been proven for all processes.
In the kinematical limit when the collinear factorization becomes valid, the result based on TMD factorization should be matched with that obtained using collinear factorization in the same process, with the inclusion of soft factors in the TMD framework [7][8][9][10][11][12]. To ensure gauge invariance, the TMDs include gauge links or Wilson lines, which comes due to initial or final state interactions [13][14][15][16]. These introduce process dependence in them. Recent results from RHIC are in favour of the theoretical prediction of a sign change of the Sivers function observed in SIDIS and DY processes, respectively [17]. More data is needed to have a firm understanding of the process dependence of the TMDs. Gluon TMDs [18] till now are far less investigated than the quark TMDs. The positivity bound gives a constraint on them [18]. Recently, there is an extraction of unpolarized gluon TMDs from LHCb data [19]. Gauge invariance of the gluon TMDs requires the inclusion of two gauge links in the definition, which makes the process dependence more involved than the quark TMDs [20]. The most common are one past and one future pointing gauge links, [+−] or [−+] (f-type) and both past or both future pointing, [−−] or [++] (d-type). The operator structures of these two types of gluon TMDs are different [20]. In the literature on small-x physics, these two gluon TMDs are called Weiszacker-Williams (WW) [21,22] type and dipole [23] type, respectively. These contribute in different processes.
In an unpolarized proton, there is a nonzero probability of finding linearly polarized gluons.
The linearly polarized gluon distributions were introduced in [18], and first investigated in a model in [24]. Recently these have attracted quite a lot of attention, although till now they have not been extracted using data. Linearly polarized gluon distributions can be probed in ep and pp collisions [25][26][27][28][29][30][31][32][33][34][35][36][37][38]. These give an azimuthal asymmetry of the form cos 2φ t [26]; also, they affect the transverse momentum distribution of the outgoing particle. Depending on the gauge links, the linearly polarized gluon distribution can be WW or dipole type. These are time reversal even (T-even) objects. In pp scattering processes, the initial and final state interactions often affect the TMD factorization; the asymmetries and cross sections also involve both f-type and d-type gluon TMDs, and disentangling the two is difficult from the observables. The gauge link structure is simpler in ep scattering processes [39], and the upcoming electron-ion collider (EIC) at Brookhaven National Lab will play an important role in probing the gluon TMDs, including the linearly polarized gluon TMD over a wide kinematical region. cos 2φ t asymmetry in J/ψ production in unpolarized ep collision has been shown to be a useful observable to probe the linearly polarized gluon TMD. Contribution to the asymmetry comes already at the leading order (LO) through the virtual photon-gluon fusion process [40]; this contributes at z = 1, where z is the fraction of the energy of the photon carried by the J/ψ in the rest frame of the proton. In the kinematical region z < 1 [41] one has to incorporate higher order Feynman diagrams. In this process, the J/ψ produced needs to be detected in the forward region, or with its transverse momentum p T not so large; otherwise TMD factorization is not expected to hold. In this work, we investigate the cos 2φ t asymmetry in a slightly different process, namely when a J/ψ and a jet are observed almost back-to-back in ep collision. Only the WW type gluon TMDs contribute in this process [42]. Here the J/ψ produced can have large transverse momentum, as the soft scale required for the TMD factorization is provided by the total transverse momentum of the J/ψ-jet pair, which is smaller than their invariant mass as they are almost back-to-back [43]. In fact, by varying the invariant mass of the pair, one can also probe the TMDs over a wide range of scales and investigate the effect of TMD evolution on the asymmetry. In [42], the upper bound of the cos 2φ t asymmetry was investigated in this process, as well as the asymmetry in the small-x region. Here we present a calculation of the asymmetry using some recent parameterization of the gluon TMDs and also investigate the effect of TMD evolution.
A widely used approach to calculate the amplitude of J/ψ production is based on an effective field theory called non-relativistic QCD (NRQCD) [44][45][46]. Here, one assumes that the amplitude for J/ψ production process can be factorized into a hard part where the cc pair is produced perturbatively, and a soft part where the heavy quark pair hadronizes to form a J/ψ. The hadronization process is encoded in the long distance matrix elements (LDMEs) [47]; which are usually extracted using the data. The cross-section is expressed as a double expansion in terms of the strong coupling α s as well as the velocity parameter associated with the heavy quark v [48,49]; in the limit v 1. For charmonium v ≈ 0.3. The heavy quark pair in the hard process is produced in different states denoted by 2s+1 L (c) J where s denotes the spin of the pair (singlet or triplet), L is the orbital angular momentum, J is the total angular momentum and (c) denotes the color configuration, which can be singlet or octet. The heavy quark pair produced in the hard process emits soft gluons to evolve into J/ψ. For the S-wave contribution, the dominant term in the limit v ≈ 0 gives the result of the color singlet model (CSM) [50,51], where the heavy quark pair in the hard process is assumed to be produced with the same quantum numbers as the J/ψ, and in the color singlet state. In our work, we include both color singlet (CS) and color octet (CO) contributions.
The paper is arranged as follows: In section II, we present the TMD formalism adopted. In Section III, we investigate the effect of the TMD evolution on the asymmetries. In sections IV and V, we present two recent parameterizations of the gluon TMDs, based on the spectator model and a Gaussian parameterization, respectively. Numerical results are presented in section VI. The conclusion is discussed in section VII.

II. FORMALISM
We consider a semi-inclusive electroproduction of a J/ψ and a jet, where the four-momentum of the particles are given in their corresponding round brackets. Here, both the incoming electron beam and the target proton are unpolarized with their respective momenta l and P . The kinematics of the process can be described in the following variables The virtuality of the scattering photon is given by Q 2 , s is the square of the electron-proton center of mass energy whereas, W is the invariant mass of photon-proton system. x B is the Bjorken-x variable, y is the inelasticity variable that gives the fraction of the energy of the electron taken by the scattering virtual photon and the variable z defines the fraction of energy of the photon carried by the outgoing J/ψ particle in proton rest frame. We consider virtual photon-proton center of mass frame where they move along +z and −z directions, respectively.
To define the kinematics, we use the light-cone coordinate system. We use two light-like vectors, one of which is chosen to be in the direction of the proton momentum, P = n − and the other n = n + , such that P · n = 1 and n 2 − = n 2 + = 0. In terms of these vectors, the momenta of the particles involved in the process can be written as follows; Momenta of the initial proton and virtual photon can be expressed as where M p is proton mass and Q 2 = x B ys. The expressions for the incoming and outgoing lepton momenta can be written in terms of light cone coordinates using the inelasticity variable y as, At the partonic level process, J/ψ can either be produced through the gluonic channel: g +γ → J/ψ + g or through the quark (anti-quark) channel: q(q) + γ → J/ψ + q(q). However, in the small-x region, the gluonic channel dominates over the quark (anti-quark) channel [42]. Hence, we have considered the contribution from the gluonic channel only in the following estimate of the asymmetry. The above processes contribute at the next-to-leading (NLO) order in α s and in the kinematic region z < 1. The outgoing energetic gluons produced in this process gives the jet. Now, we can express the momentum of the initial gluon, p g , the final J/ψ and jet with momentum P ψ and P j , respectively, in terms of the light-like vectors as where, x is the collinear momentum fraction of the initial gluon, P ψ⊥ and P j⊥ are the transverse momenta of J/ψ and jet respectively. The incoming and the outgoing scattered lepton form the leptonic plane. All the azimuthal angles of the final state particles are defined with respect to the leptonic plane with φ l = φ l = 0.
For the process under consideration, TMD factorization has not formally been proven yet, although it is expected to be valid. In our study we have assumed TMD factorization. The total differential scattering cross-section for the unpolarized process: ep → J/ψ Jet X can be written as [26] dσ = 1 2s The function M µν represents the scattering amplitude of J/ψ production in the photon-gluon fusion process: γ * (q) + g(p g ) → QQ(P ψ ) + g(P j ) partonic subprocess. The leptonic tensor L µµ describes the electron-photon scattering and can be written as, where e is represents the electronic charge. The gluon correlator, Φ νν g (x, p 2 T ) describes the gluon content of the proton. At the leading twist, for the case of an unpolarized proton, it can be parameterized in terms of two TMD gluon distribution functions as [18] Φ νν g (x, p 2 Here, g νν ⊥ = g νν − P ν n ν /P · n − P ν n ν /P · n. The quantities f g 1 (x, p 2 T ) and h ⊥g 1 (x, p 2 T ) represent the unpolarized and linearly polarized gluon TMD, respectively.

A. J/ψ production in NRQCD framework
The Feynman diagrams for the dominant subprocess of photon-gluon fusion, which results in the production of a J/ψ and a jet, are shown in Fig. 1.
The amplitude for the production of J/ψ within the NRQCD framework can be written as follows [36,51] where k is the relative momentum of the heavy quark or the anti-quark in the rest frame of the non-relativistic quarkonium bound state, which is assumed to be very small as compared with the rest mass of the quarkonium. Here Ψ LLz (k) is the nonrelativistic bound-state wave function with orbital angular momentum L, L z . The Clebsch-Gordan coefficients LL z ; SS z |JJ z projects out their angular momentum. The mass of the quarkonium, M ψ , is taken to be twice the heavy quark mass. The O(q, p g , P ψ , k) represents the amplitude for the production of the heavy quark anti-quark pair, QQ, without the inclusion of the polarization of the quark and anti-quark. This can be calculated by considering the contributions from all the above Feynman diagrams and can be written as where, i denotes the contribution from the individual Feynman diagrams given in Fig. 1 and C i corresponds to the color factor for each diagram. The O i (q, p g , P ψ , k) for the above Feynman diagrams are written as The expressions for the remaining, O 5 , O 6 , O 7 and O 8 can be obtained by reversing the fermionic current and replacing k to −k.
In the NRQCD framework, the outgoing QQ pair can be formed in the color singlet (CS) state or in the color octet (CO) states. The color factors, C i , corresponding to the CO case are given as [52], The SU(3) Clebsch-Gordon coefficients for CS and CO states are given by, respectively, where N c is number of colors, t c is the generators of SU(3) in the fundamental representation.
Their properties are given by T r[t a t b ] = δ ab /2 and T r[t a t b t c ] = 1 4 (d abc + if abc ). By using these relations along with the relation in Eq. (17), one obtains the color factors for individual Feynman diagrams for the QQ formed in color octet states as For the case of formation of QQ pair in CS state, the color factors are given by The spin projection operator for the bound state of the J/ψ includes the spinors of the heavy quark and anti-quark, cc and is given as where, Π SSz = γ 5 for spin singlet (S = 0) and Π SSz = / ε Sz (P ψ ) for spin triplet (S = 1). The ε Sz is the spin polarization vector of the outgoing cc pair. Now, since, in the rest frame of the bound state, k P ψ , thus one can Taylor expand the amplitude given in Eq. (13) around k = 0 limit. In that expansion, the terms corresponding to k 0 give S-wave scattering amplitude (L = 0, J = 0, 1) and terms linear in k correspond to P-wave scattering (L = 1, J = 0, 1, 2) and the corresponding amplitudes are given as Where R L represents the radial wave function and the shorthand notations used in the above expressions are Following are the relations for the Clebsch Gordon coefficient and the polarization vector of J/ψ which we can use to calculate P wave amplitudes [53] LzSz The ε α Jz (P ψ ) is the polarization vector corresponding to J = 1 angular momentum state, that follows the current conservation and obeys the following relations [53], The ε αβ Jz (P ψ ) is the polarization tensor corresponding to J = 2 which is symmetric in the Lorentz indices and follows the relations [53], The radial wave function and its derivative evaluated at origin R 0 (0), R 1 (0) given in Eq. (21) and Eq. (22) are related to LDMEs by the following equations [42], The two different sets of LDMEs that we have used in our numerical calculation are given in Table. I Ref. [55] 8.9 ± 0.98 0.30 ±0.12 1 ) and CO states 1 state can be written as, .
The symmetry relations, given in Ref. [56], lead to the cancellation of contributions from Feynman diagrams 4 and 8 for 3 S The total amplitude for 1 S (8) 0 state can be written as, (36) and amplitude can be written as [56], E. Total differential cross-section and the asymmetry The structure of differential cross-section as defined in Eq. (10) has a contraction of tensors which can be schematically written as where i, j = 1, 2, 3, 4 (corresponds to the Feynman diagrams given in Fig. 1 and we have We have already defined all tensors in the above convolution and contributions to the amplitudes come from all the CS and CO states ( 3 S J(0,1,2) ). We have summed over all the polarization states of the outgoing gluon using the relation given where, n ψµ = P ψµ /M ψ . We use the frame where the incoming virtual photon and proton moves along z-axis. The azimuthal angles of the lepton scattering plane is defined as φ l = φ l = 0. We integrate out the azimuthal angle of the final lepton [57], we can write Moreover, for the other phase factors, one can write and conservation of the four-momenta can be written as where we have used the relation Q 2 = x B ys. Now, we define the sum and difference of the transverse momentum of J/ψ and jet as From Eq. (44), we havez = (1 − z) and q t = p T . We use TMD factorization in the kinematic region |q t | |K t |. This leads to the situation where the outgoing J/ψ and jet are almost back to back in the transverse plane, i.e. the xy plane, w.r.t the virtual photon-proton colliding axis, i.e. the z-axis, thus allowing us to set K t P ψ⊥ −P j⊥ . φ t and φ ⊥ are the azimuthal angles, respectively, for q t and K t and are defined with respect to the leptonic plane as illustrated in  Finally, integrating overz, p T and x, the differential cross-section, in Eq. (10), as a function of z, y, x B , q t , K t can be rewritten as In the following calculation, we have only kept the zeroth and first order terms in we obtain the total differential cross-section [42], The analytic expressions of the coefficients A i 's and B i 's are very lengthy, which we have not included in this article. They are available upon request.
The TMDs come with various azimuthal modulations. These modulations can be used to extract information about the ratio of TMDs. This can be done by defining the asymmetries as [42] A Here, we are interested in one particular asymmetry, i.e., cos 2φ t asymmetry to extract the linearly polarized gluon TMD. This can be written as, Now, by plugging the differential scattering cross-section from Eq. (47) in the above equation, we get the cos 2φ t asymmetry as function of z, y, For estimating the cos 2φ t asymmetry numerically, we have used two recent parameterizations of the gluon TMDs, and also estimated the effect of TMD evolution. The following section gives the details of the TMD evolution formalism used.

III. TMD EVOLUTION
The evolution of the TMDs [58,59] with the scale affects the asymmetries measured in the energies of different experiments, and it is important to estimate the effect of this evolution to obtain the angular asymmetries of produced hadrons measured by the HERMES, COMPASS, JLab as well as the future EIC experiments at different energies. The TMD evolution is usually studied in the impact parameter space [58]. The impact parameter dependent TMDs can be written as Fourier transforms of the TMDs, In the TMD evolution approach, TMDs not only evolve with the intrinsic transverse momentum of the parton but also evolve with the probing scale. The expression for TMD evolution at a given final scale Q f can be obtained by solving Collins-Soper evolution equation and renormalization group equation. Using this approach, the expression for the gluon TMD in the impact parameter space can be written as [7,60,61] f partons like quark/anti-quark or gluon. The exponents S A and S np are the perturbative and non-perturbative Sudakov factors, respectively. We note that the Sudakov factor S A is spin independent, and thus the same for all (un)polarized TMDs [35,62]. As stated before, a formal proof of TMD factorization for this process remains to be done, and here we study the effect of TMD evolution on the asymmetry assuming such factorization. The subsections below contain a description of the TMD evolution formalism used and the relevant formulas.

A. Coefficient functions and perturbative Sudakov factor
The coefficient function can be written as a series of strong coupling constant α s [63], The Sudakov factor as the leading order of α s can be written as [63], The running of the coupling α s is ignored in Eq. (54) because it starts at α 2 s . The C A = N c , n f denotes the number of active flavours. In b t 1/Λ QCD the Sudakov factor can be Taylor expanded. Further, by substituting the expressions from Eq. (53) and Eq. (54) the unpolarized gluon TMD at LO without non-perturbative Sudakov factor is given as [63], Scale evolution of the collinear PDFs are given by the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equation. Using this evolution equation, one can evolve f g Here, P gg and P gi are leading order splitting functions which are given as, Where, C F = (N 2 c − 1)/2N c and T R = 1/2. In Eq. (57), the first term involves the plus prescription and thus avoids an infrared divergence because of (1 −x) in denominator. The plus prescription is given as [63], The ⊗ symbol denotes convolution of the two quantities; After convolution and substitution of the Eq. (56) in the Eq. (55), we have the final equation for the unpolarized gluon TMD as, Now, we can write the above equation in the q t -space by making one-to-one correspondence between the functions in impact parameter and momentum space using a general formula [35,64,65] here, n is the rank of function in q t -space. Since the unpolarized vector-meson production generally has a rank-zero structure, we can write the unpolarized gluon TMD together with the non-perturbative Sudakov factor in terms of q t -space as Now let us write the expression for the linearly polarized gluon distribution function h ⊥g 1 (x, b 2 t ). The perturbative tail of h ⊥g 1 can be computed in the same way as the perturbative tail of f g 1 , with the key difference that its expansion in powers of the QCD coupling constant begins at O(α s ). Using Eq. (3.13) of Ref. [63] and then performing the Fourier transformation using Eq. (62), we write the linearly polarized gluon distribution at LO in terms of the unpolarized Here we have used the general formula for Bessel integral [66] ∞ where J ν (z) is the Bessel function of the first kind of order ν. The q t dependence for h ⊥g 1 in the gluon correlator has a rank-two tensor structure in the non-contracted transverse momentum, thus we could write the linearly polarized gluon TMD h ⊥g 1 in the q t -space is given as,

B. Non-perturbative Sudakov factor
In the above Eqs. 63 and 66, the perturbative part is strictly valid in the perturbative domain, which means low b t (or b t << 1/Λ QCD ). However, in order to perform the corresponding Fourier transform, we need to integrate the expression from small to large b t . As a result, the perturbative expression for the Sudakov factor given above should not be used alone rather one needs to introduce a non-perturbative Sudakov factor, this should suppress the large b t domain.
The functional form of the non-perturbative Sudakov factor is constrained by two conditions, one of which is it has to be equal to 1 for b t = 0 and for large b t it is supposed to decrease monotonically and ultimately should vanish. The functional form attributed to S np is quadratic in b t with e −Snp reaching 0 within a certain value of b t called b tlim . However, when b t gets too small, the lower scale Q i = 2e −γ E /b t will be larger than final scale Q f and hence, we expect the evolution should stop. This could be resolved by taking a b t prescription as given below [67],

This prescription constrains the
and Q f (for b t → 0). Motivated by [67] we choose the Gaussian behavior of e −Snp as , The parameter A controls the width of the non-perturbative Sudakov factor for a particular Q f .
In this calculation, we have taken the final scale as Q f = M 2 ψ + K 2 t . We have taken A = 2.3 GeV 2 , which is calculated for K t = 1 GeV. The K t dependence present in Q f does not affect the value of A, as K t increases the e −Snp remains below the convergence criteria at the given b tlim .
The value of b tmax for the following numerical study is 1.5 GeV −1 which is in consistent with [67], Collins-Soper-Sterman formalism given for Z Boson production [68] and Collins-Soper-Sterman formalism implemented by [69]. Below we present two recent parameterizations of the gluon TMDs that we have used.

IV. SPECTATOR MODEL
In this section, we discuss a recent parameterization of the gluon TMDs based on a spectator model [70]. According to this model, the remnant after the gluon emission from the nucleon is treated as a single spectator particle, which is on-shell, with mass M X . The mass can take a range of values given by a spectral function. The nucleon-gluon-spectator coupling is encoded in an effective vertex that contains two form factors. The expression for a given TMD reads as Here ρ X (M X ) is the spectral function and can be written as  Table II below. These model parameters have been fixed by fitting the NNPDF data at scale Q = 1.64 GeV [70]. We assumed the same set of parameters to probe the TMDs at Q = M 2 ψ + K 2 t . In this model, we don't have direct inference of scale dependency on the TMDs unlike the case of Gaussian parameterization. However, the longitudinal momentum fraction x depends on Q but going from some low Q = 1.64 GeV to some relatively large Q = M ψ + K 2 t (≈ 6.75 GeV for M ψ = 3.1 GeV and K t = 6 GeV) hardly changes the x range.
Parameter Replica 11 Parameter Replica 11 The leading-twist T -even unpolarized and linearly polarized gluon TMDs can be written as [18,24]f where p is the gluon momentum and The form factors given above are dipolar in nature; the main advantage of using dipolar form factors consists in the possibility of canceling gluon-propagator singularities, quenching the effects of large transverse momenta where a pure TMD description is not any more adequate, and removing logarithmic divergences emerging in p t -integrated densities.

V. GAUSSIAN PARAMETERIZATION OF THE TMDS
The most widely used parameterization of the TMDs are Gaussian in nature, Here, both the TMDs, f g 1 and h ⊥g 1 are assumed to be factorized into a product of a x dependent part given in terms of the collinear pdfs and an exponential factor which is a function of only the transverse momentum (q t ). The width of the Gaussian is usually expressed in terms of the average value of the transverse momentum, which is taken as a model parameter [36][37][38] where, r(0 < r < 1) is a parameter and in our case we take r = 1/3. The term f g 1 (x, µ) is the collinear PDF which follows the DGLAP evolution equation. The Gaussian width here is q 2 t = 0.25 GeV 2 . The linearly-polarized gluon distribution in the parameterization above satisfies the positivity bound [18], but does not saturate it.

VI. RESULTS AND DISCUSSION
In the present work, we have numerically calculated the cos 2φ t azimuthal asymmetry in the unpolarized electroproduction of J/ψ process: e p → e J/ψ Jet X, within TMD factorization approach. As depicted in Fig. 2, we have J/ψ and jet almost back-to-back in the transverse plane as we consider the kinematics, |q t | |K t | which is the required condition to assume TMD factorization. Contributions from the virtual photon-quark(anti-quark) initiated sub-processes in the unpolarized cross-section is very small in the kinematics considered [42] as compared with gluon initiated subprocess: γ * + g → J/ψ + g. Therefore, in the numerical estimate of the asymmetry, we have included only the gluon-photon fusion subprocess and neglected the contribution from the quark(anti-quark) initiated sub-processes. We have used MSTW2008 [71] set of collinear PDFs and adopted two sets of LDMEs for the study of azimuthal asymmetry as listed in Table I with the charm mass m c = 1.3 GeV. We used NRQCD framework for J/ψ production rate and included contributions from both color singlet and color octet states in the asymmetry. We also calculated the asymmetry taking into consideration contribution only from J(=0,1,2) is calculated using FeynCalc [72,73]. We have investigated the effect of TMD evolution on the asymmetry. For the gluon TMDs we use two parameterizations, Gaussian and based on the spectator model, as discussed above. The J/ψ mass is taken to be M ψ = 3.1 GeV. The collinear PDFs are evaluated at the scale Q = M 2 ψ + K 2 t . The numerical results are presented in the kinematical region that can be accessed at the future EIC.
We have imposed a cut on the variable z, namely, 0.1 < z < 0.9, to estimate the asymmetry.
As z → 1, the final state gluon becomes soft, which leads to infrared divergences. We impose the upper cut to avoid this gluon to become soft. Contribution of J/ψ production from fragmentation of final hard gluon comes from lower z region, and we impose the lower cut to minimize this contribution. The asymmetry gets maximized around z = 0.7 for the kinematics we have considered. Hence, we took z = 0.7 for all plots where z is fixed. In addition, we also show the cos 2φ t asymmetry as function of z. In our estimate, we have neglected the contribution of J/ψ production via feed-down from excited ψ(2S) and the decays of χ c states.   We have used CMSWZ set of LDMEs [55].  We have used CMSWZ set of LDMEs [55].
In the upper left panel (a) of Figs. 3-5, we have plotted the cos 2φ t asymmetry as functions of K t at √ s = 140 GeV. In these plots, we integrated q t and y in the range (0, 1) and, (0.1, 1) respectively. We see that the cos 2φ t asymmetry is maximum (negative) for lower K t and monotonically decreases as we go in the higher K t region. The maximum asymmetry we obtained is ≈ 29% in the spectator model at K t = 1 GeV followed by the Gaussian parameterization, ≈ 17%. Incorporation of the TMD evolution results in a smaller asymmetry, ≈ 1% at K t = 1 GeV.
The y dependence of cos 2φ t azimuthal asymmetry is shown in the upper right panel (b) of Figs. 3-5 at √ s = 140 GeV; we have shown the results both in NRQCD and in the CS. We plotted the asymmetry for two fixed values of K t , namely 2 GeV and 4 GeV. We have plotted the asymmetry in the range of y ∈ [0.1, 1], however in the lower y region, the magnitude of the asymmetry is similar in both the spectator and the Gaussian models, whereas in TMD evolution approach, the asymmetry is maximum around y = 0.44 at K t = 2 GeV in NRQCD.
The asymmetry is small in the TMD evolution approach; however, we obtain a significant asymmetry, ≈ 29%, in the spectator model followed with ≈ 18% in the Gaussian model at We have taken fixed values of y, namely, 0.1 and 0.8. The peak of the asymmetry is ≈ 12% at z ≈ 0.7 in NRQCD and ≈ 2% at z ≈ 0.4 in the CS at y = 0.1. In all these plots, we have used LDMEs from [55].  In Fig. 6, we have plotted both the TMDs, f g 1 and h ⊥g 1 and their ratio, as functions of q t for all three parameterizations. In all these plots we have used similar kinematics as we considered for the plots in Fig. 3-7. In this kinematics, the x value of the gluon TMDs are of the order of 10 −3 − 10 −2 . We have plotted at the probing scale which is the virtuality of photon, where K t = 3 GeV and at the fixed values of y = 0.3 and z = 0.7. This sets x ≈ 0.012. From the plots of the ratio, Fig. 6), we see that the TMDs in the spectator model indeed saturates the positivity bound, whereas, the Gaussian parameterizations and TMD evolution approach satisfy the positivity bound but do not saturate it except for q t ≈ 0.36 GeV, where Gaussian parameterization is saturating the positivity bound. Moreover, the ratio is larger in the case of Gaussian as compared with TMD evolution approach for almost whole range of q t considered. In the Spectator model, tails of the TMDs in the small-x domain, depend on the trend of spectral function at large M X [70]. We have checked that the spectator model results in Eqs. (12) and (15) of [70] for the ratio we have fixed the virtuality of photon Q = M 2 ψ + K 2 t . We can see that the asymmetry calculated with the spectator model is maximum and agrees with the upper bound. As seen from Fig. 7, the asymmetry incorporating TMD evolution is significantly smaller than that calculated using Gaussian and spectator models for the gluon TMDs. This is because the denominator of the asymmetry receives contribution from the unpolarized gluon distribution which has a leading order (LO) term Eq. (61) whereas the numerator contains the linearly polarized gluon distribution whose leading contribution comes at O(α s ).
If we exclude the LO term in the unpolarized gluon TMD, we find the asymmetry increases approximately three times. We also find that the magnitude of the asymmetry does not change much if √ s is lower.
j . One can see from (a) that for the LDME set CMSWZ [55], the dominating contribution comes from one single state, 1 S (8) 0 ; while from (b) one can see that for SV set of LDME [54], the dominating contribution comes from two states, 1 S j . So, we can conclude that the asymmetry depends on the LDME set chosen. It is worth mentioning here that our results for the upper bound of the asymmetry do not match with those presented in [42] even if we plot it using the same scale as in this reference. We have traced this mismatch to a difference in the sign of the contribution coming from the 3 P

VII. CONCLUSION
We have presented a calculation of the cos 2φ t asymmetry in almost back-to-back production of a J/ψ and a jet in ep collision, using TMD factorization and generalized parton model. This asymmetry is sensitive to the still unknown linearly polarized gluon distribution. We present a numerical estimate of the asymmetry in the kinematical region that will be accessible at the future EIC. We have used NRQCD to calculate the J/ψ production rate and two recent parameterization for the gluon TMDs, one based on a Gaussian type distribution and another based on spectator model. The asymmetry is quite sizable; in fact in spectator model the asymmetry agrees with the upper bound that is obtained by saturating the positivity condition of the gluon TMDs. TMD evolution affects the asymmetry at the energy of the EIC, making it smaller. The asymmetry also depends on the LDMEs used, and dominating contribution comes from different states. We conclude that the back-to-back production of J/ψ and a jet at the