Gluon TMDs and NRQCD matrix elements in $J/\psi$ production at an EIC

In this paper we analyze azimuthal asymmetries in the processes of unpolarized and polarized $J/\psi \,(\Upsilon)$ production at an Electron-Ion Collider. Apart from giving access to various unknown gluon transverse momentum distributions, we suggest to use them as a new method to extract specific color-octet NRQCD long-distance matrix elements, i.e.\ $\langle0\vert{\cal O}_{8}^{J/\psi}(^{1}S_{0})\vert0\rangle$ and $\langle0\vert{\cal O}_{8}^{J/\psi}(^{3}P_{0})\vert0\rangle$, whose values are still quite uncertain and for which lattice calculations are unavailable. The new method is based on combining measurements of analogous asymmetries in open heavy-quark pair production which can be performed at the same energy. To enhance the gluon contribution one can consider smaller values of $x$ and, in order to assess the impact of small-$x$ evolution, we perform a numerical study using the MV model as a starting input and evolve it with the JIMWLK equations.

( 3 P0)|0 , whose values are still quite uncertain and for which lattice calculations are unavailable. The new method is based on combining measurements of analogous asymmetries in open heavy-quark pair production which can be performed at the same energy. To enhance the gluon contribution one can consider smaller values of x and, in order to assess the impact of small-x evolution, we perform a numerical study using the MV model as a starting input and evolve it with the JIMWLK equations.

I. INTRODUCTION
Transverse momentum dependent parton distributions (TMDs) are fundamental objects which encode information on the motion of partons inside hadrons and on the correlations between spin and partonic transverse momenta. As such, they can be considered as an extension of the standard, one-dimensional, parton distribution functions (PDFs) to the three-dimensional momentum space. Contrary to PDFs, TMDs are in general not universal. This is due to their sensitivity to the soft gluon exchanges and the color flow in the specific process in which they are probed. A typical example is provided by the Sivers function for quarks [1], namely the azimuthal distribution of unpolarized quarks inside a transversely polarized proton, which is expected to enter with opposite sign in the single spin asymmetries for semi-inclusive deep inelastic scattering (SIDIS) and for the Drell-Yan processes [2,3]. More recently, a similar sign change test has been proposed for the gluon Sivers function as well [4,5]. Experimental verification of these properties would strongly corroborate our present understanding of the structure of the proton and nonperturbative QCD effects.
Among gluon TMDs, the distribution of linearly polarized gluons inside an unpolarized proton [6][7][8] has attracted a lot of attention in the last few years. It corresponds to an interference between +1 and −1 gluon helicity states which, if sizable, can affect the transverse momentum distributions of final state particles like, for instance, the Higgs boson [9][10][11]. Linearly polarized gluons have been investigated theoretically in the dilute-dense regime in protonnucleus and lepton-nucleus collisions as well [12][13][14][15][16][17][18][19]. Very interestingly, it turns out that at small-x fractions of the gluons inside a nucleus, the linearly polarized distribution may reach its maximally allowed size, bounded by the unpolarized gluon density [6], although it depends on the process whether the observable effects are maximal [20].
From the experimental point of view, almost nothing is known about gluon TMDs, because they typically require higher-energy scattering processes and are harder to isolate as compared to quark TMDs. Many proposals have been put forward to access them by looking at transverse momentum distributions and azimuthal asymmetries for bound or open heavy-quark pair production, both in lepton-proton and in proton-proton collisions. The reason is that heavy quarks are very sensitive to the gluon content of hadrons, as is well known from studies of gluon PDFs. A first Gaussian shape extraction of the unpolarized TMD gluon distribution has been recently performed from LHCb data on the transverse spectra of J/ψ pairs [21].
In a series of papers [4,22,23], the process e p → e Q Q X, with Q being either a charm or a bottom quark, has been considered as a tool to extract gluon TMDs at a future Electron-Ion Collider (EIC) [24][25][26]. The observables, needed to disentangle the five different gluon TMDs contributing to the unpolarized and transversely polarized cross sections, have been properly defined, each one of them corresponding to a specific azimuthal modulation. Moreover, especially in Ref. [4], attention has been paid to the small-x behavior of all the distributions and to their process dependence, by relating them to other reactions which could be measured, for example, at the proposed fixed target experiment AFTER@LHC [27,28]. It is natural at this point to perform a similar analysis for the case in which the two heavy quarks form a bound state. We therefore consider here inclusive J/ψ and Υ production in deepinelastic lepton-proton scattering, namely e p → e J/ψ (Υ) X, where the electron is unpolarized and the proton can be either unpolarized, or polarized transversely to the electron-proton plane. In addition to unpolarized quarkonium production, we examine the cases in which the quarkonium state is polarized either longitudinally or transversely with respect to its direction of motion in the γ * p center-of-mass frame, with γ * being the virtual photon exchanged in the reaction. Analogous studies, although limited to the Sivers and linearly polarized gluon densities and to unpolarized quarkonium production, have been published recently [29,30].
In the present analysis we adopt the TMD framework in combination with nonrelativistic QCD (NRQCD) [31][32][33], which is the effective field theory that allows for a factorized treatment of the heavy-quark pair production, calculable in perturbative QCD, and the nonperturbative hadronization process leading to the binding of the pair, encoded in long-distance matrix elements (LDMEs) [34]. Since these LDMEs, which are assumed to be universal, obey specific scaling rules in the average velocity v of the heavy quark in the quarkonium rest frame [35], the corresponding cross section can be evaluated through a double expansion in the strong coupling constant α s and in the velocity v, with v 2 0.3 for charmonium and v 2 0.1 for bottomonium. The heavy quark-antiquark pair can be produced directly in a color-singlet (CS) configuration, with the same quantum numbers as the observed quarkonium, but also as a color-octet state (CO) with different quantum numbers. In the latter case, the pair becomes colorless after the emission of soft gluons. The CS LDMEs are commonly obtained from potential models [36], lattice calculations [37] or from leptonic decays [38], while the CO ones are usually determined by fits to data on J/ψ and Υ yields [39][40][41][42][43], but not from lattice calculations. As a result, at present our knowledge of the CO matrix elements is not very accurate (cf . Tables I and II below). Moreover, although NRQCD successfully explains many experimental observations, it has problems to reproduce all cross sections and polarization measurements for charmonia in a consistent way [44,45]. As a consequence, alternative approaches to NRQCD are used as well, also in TMD studies. For instance, J/ψ photoproduction as a way to access the gluon Sivers function [46][47][48] has been studied in the so-called Color Evaporation Model [49], which is based on quark-hadron duality and assumes that the probability to form a physical (colorless) quarkonium state does not depend on the color and the other quantum numbers of the hadronizing QQ pair.
In the specific kinematical configuration considered, i.e. P QT M Q , the CO production mechanism is expected to be the dominant one [50]. This is an important aspect, because some of our proposed observables, namely the single spin asymmetries, are expected to vanish in semi-inclusive deep-inelastic scattering in the CS mechanism, due to the absence of any initial or final state interactions [33]. Most of the previous studies on gluon TMDs in protonproton collisions focussed on scattering processes in which the CS production mechanism is the dominant one, such as p p → η c,b X, p p → χ 0c,b (χ 2c,b ) X [51], p p → J/ψ (Υ) γ X [52], p p → J/ψ (Υ) ¯ X [53] and p p → J/ψ J/ψ X [21]. The reason to concentrate on these CS dominated processes is to avoid the presence of final state interactions which, together with the initial state interactions present in proton-proton collisions, would lead to the breaking of TMD factorization [54]. Furthermore, as already discussed in Ref. [4], the gluon distributions extracted in e p → e J/ψ (Υ) X or in e p → e Q Q X, which correspond to the so-called Weizsäcker-Williams (WW) distributions in the small-x limit, are all related to the TMDs entering in the above mentioned proton-proton reactions and differ from them by, at most, an overall minus sign.
Investigating the process e p → e J/ψ (Υ) X can be also very helpful to improve our understanding of the mechanisms underlying quarkonium production. To this end, here we propose a new method to extract, apart from various TMDs that are at present still unknown from the experimental point of view, also the dominant CO LDMEs, namely 0|O J/ψ 8 ( 1 S 0 )|0 and 0|O J/ψ 8 ( 3 P 0 )|0 , by combining measurements of azimuthal asymmetries in e p → e J/ψ (Υ) X, with analogous ones in e p → e Q Q X. In this way, heavy-quark final states at an EIC can contribute to the determination of the CO LDMEs. We should remark that, in order to avoid having to deal with evolution in this comparison, one should consider the same value of the photon virtuality Q 2 in both processes. In the first one, e p → e J/ψ (Υ) X, we consider the transverse momentum P QT of the produced quarkonium small with respect to the hard scale, which is identified with the quarkonium mass M Q ≈ 2M Q . In order to avoid the presence of two very different hard scales, we take Q = 2M Q . Although one can take the same Q value in the second process, e p → e Q Q X, there will be another hard scale given by the transverse momentum K ⊥ of each heavy quark, which we assume to be K ⊥ = Q = 2M Q for simplicity.
The processes considered in this paper are gluon induced, and are therefore expected to be enhanced when consid-ering smaller x values. At an EIC, the smaller the x value, the smaller the Q values covered, so one has to keep a balance between the x and Q ranges. For the J/ψ case one can go to lower x values. Since we consider only a limited Q range, we will not include TMD evolution, although this can be done along the lines considered in Ref. [29]. To assess the less studied influence of evolution in x, we perform a numerical study of the implications nonlinear small-x evolution would have in the range from x ∼ 10 −2 to x ∼ 10 −4 covered by the EIC at low Q values of a few GeV. It turns out to have only a moderate suppression effect. This study is limited to the unpolarized proton case, for which nonperturbative models are available for the corresponding small-x gluon distributions [55][56][57], which we use as the initial condition for the evolution. The Color Glass Condensate effective theory [58] makes it then possible to calculate the nonlinear evolution in rapidity of these distributions, in the presence of saturation. This was done with the help of a numerical implementation of JIMWLK on the lattice in Refs. [17,18,59]. We use the results therein obtained for the unpolarized and linearly polarized WW TMDs inside an unpolarized hadron, to show predictions for our azimuthal modulations at different values of rapidity in the low-x limit. The paper is organized as follows. In Section II we provide the operator definition of gluon TMDs and discuss their process dependence. The derivation of the cross section for unpolarized quarkonium production in DIS, within the TMD framework, can be found in Section III. Further details of the calculations are relegated to Appendix A. The azimuthal moments providing direct access to gluon TMDs are defined in Section IV. Similar observables for polarized quarkonium production are discussed in Section V. Our strategy for the extraction of the CO LDMEs, based on the combination of azimuthal asymmetries for bound and open heavy-quark pair production, is described in Section VI. Upper limits of the azimuthal moments, as well as an analysis of the small-x evolution of gluon TMDs and the cos 2φ asymmetries, are presented in Section VII. Summary and conclusions are given in Section VIII.

II. OPERATOR DEFINITION OF GLUON TMDS
The transverse momentum distribution of a gluon with four-momentum p inside a proton with four-momentum P and spin vector S can be defined as follows. We first perform a Sudakov decomposition of p and S in terms of P and a light-like vector n, conjugate to P . Namely, where M p is mass of the proton and S 2 T = −S 2 T , with 0 ≤ S 2 L , S 2 T ≤ 1, such that S 2 L + S 2 T = 1. We then introduce the following matrix element of a correlator of the gluon field strengths F µν (0) and F νσ (ξ), evaluated at fixed light-front (LF) time ξ + = ξ·n = 0, with U [0,ξ] and U [0,ξ] being two process dependent gauge links (or Wilson lines) that are needed to ensure gauge invariance. By means of the symmetric and antisymmetric transverse projectors, respectively given by the correlator in Eq. (3) can be parametrized in terms of gluon TMDs [6][7][8]. For an unpolarized proton, one has where f g 1 (x, p 2 T ) is the TMD unpolarized distribution and h ⊥ g 1 (x, p 2 T ) is the distribution of linearly polarized gluons. They are both T -even, i.e. they can be nonzero even in those processes where there are neither initial nor final state interactions. The correlator for a transversely polarized proton can be parametrized in terms of five independent gluon TMDs as follows where the symmetrization operator is defined as p {µ q ν} = p µ q ν + p ν q µ . The three gluon TMDs that appear in its symmetric part, (Γ µν is the gluon Sivers function, while the h functions are chiraleven distributions of linearly polarized gluons inside a transversely polarized proton. In analogy to the transversity function for quarks, we define the combination which however, in contrast to quark transversity, vanishes upon integration over transverse momentum [4]. Because of the definition in Eq. (3), the TMDs introduced in Eqs. (6) and (7) will depend on the gauge links, the specific structure of which is determined by the process under consideration. In this case, as in e p → e Q Q X [4], the partonic reaction γ * g → Q Q probes gluon TMDs with two future pointing Wilson lines, denoted as + links. In the small-x limit they correspond to the WW distributions. As already pointed out in Ref. [4], these TMDs can be related to the ones having two past-pointing, or − gauge links, which could be accessed in processes like p p → γ γ X in the back-to-back correlation limit [60]. More specifically, the T -even unpolarized and linearly polarized gluon TMDs are expected to be the same in the two kind of processes, while the T -odd densities, like the gluon Sivers functions, should be related by a minus sign. On the other hand, gluon TMDs with both a + and − link (future and past pointing), corresponding to the dipole distributions at small x, cannot be related to the TMDs discussed here. They could be accessed in processes like p p → γ * jet X [19], in the kinematic region where gluons in the polarized proton dominate, such that the partonic channel q g → γ * q is effectively selected [20]. However, TMD factorization for p p → γ * jet X has not been established so far.

III. OUTLINE OF THE CALCULATION
We study the process where Q is either a J/ψ or a Υ meson, the incoming proton is polarized with polarization vector S, and the other particles are unpolarized. We choose the reference frame such that both the virtual photon exchanged in the reaction and the incoming proton move along theẑ-axis, and azimuthal angles are measured w.r.t. to the lepton scattering plane, such that φ = φ = 0. Moreover, in order to apply a framework based on TMD factorization, we consider only the kinematic region in which the component of the quarkonium momentum transverse w.r.t. the lepton plane, denoted by q T ≡ P QT , is small compared to the virtuality of the photon Q and to the mass of the quarkonium M Q . The differential cross section can be written as where s = ( + P ) 2 ≈ 2 · P is the total invariant mass squared and Q 2 = −q 2 ≡ −( − ) 2 . Moreover, the gluon correlator Γ g is defined in Eq. (3) and the leptonic tensor L( , q) is given by with e the electric charge of the electron. The calculation proceeds along the same lines of Ref. [23], which we summarize for completeness in the following. We start with introducing the light-like vectors n + and n − , which obey the relations n 2 + = n 2 − = 0 and n + · n − = 1. Then we note that the four-momenta P and q can be written as where x B is the Bjorken-x variable, with x B = Q 2 /2P · q up to target mass corrections. We will thus perform a Sudakov decomposition of all the momenta in the reaction in terms of n + = P and n − = n = (q + x B P )/P · q. Therefore, the leptonic momenta can be written as where we have introduced the inelasticity variable y = P ·q/P · , such that the following relations hold: The invariant mass squared of the virtual photon-target system is defined as W 2 = (q + P ) 2 , and can be expressed in terms of the other invariants: Similarly, the gluon momentum can be expanded as where x = p · n, while for the momentum of the quarkonium state Q we have with z = P Q · P/q · P and P 2 QT = −P 2 QT . In a reference frame in which azimuthal angles are measured w.r.t. the lepton plane (φ = φ = 0), denoting by φ S , φ T the azimuthal angles of the three-vectors S T and P Q T , respectively, the phase-space elements in Eq. (10) can be written as Furthermore, using the Sudakov decomposition of the gluon momentum in Eq. (15), the δ-function in Eq. (10) can be re-expressed as Therefore, upon integration over the variables x, z and p T , the cross section takes the final form with z fixed to the value z = 1, the transverse momentum of the incoming gluon equal to that of the quarkonium, and its longitudinal momentum fraction x given by Within the framework of NRQCD, at leading order in the strong coupling constant α s , the partonic subprocess that contributes to J/ψ production is γ * g → QQ[ 2S+1 L (8) J ], as depicted in Fig. 1, where we have used a spectroscopic notation to indicate that the QQ pair forms a bound state with spin S, orbital angular momentum L and total angular momentum J. The additional superscript (8) denotes the color configuration. The relevant CO LDMEs are 0|O configuration in association with a gluon. As pointed out in Ref. [50], the CS contribution is suppressed relatively to the CO by a perturbative coefficient of the order α s /π. On the other hand, 0|O respectively. Hence, according to the NRQCD scaling rules, the CO contribution should be enhanced by about a factor v 3 π/α s ≈ 2 with respect to the CS one. This factor becomes ≈ 4 in the actual numerical analysis presented in Ref. [50] for values of Q 2 > 4 GeV 2 . A further suppression of the CS contribution can be achieved by applying a cut on the variable z, for instance by taking z ≥ 0.9, because at high z the CS term is known to become negligible [50] and will be therefore neglected it in our analysis. Within these approximations, the final unpolarized and transversely polarized cross sections read and with the normalization factor N given by where e Q is the fractional electric charge of the quark Q. Details of the derivation can be found in Appendix A.
Expressions for A U , B U and A T have also been given in Ref. [29], where some power suppressed terms were included, as well as an additional power suppressed cos φ T amplitude. The explicit expressions of the terms A U/T in Eqs. (21) and (22) read where the subscripts U + L, L, T refer to the specific polarization of the photon [23,61]. If we denote by A λγ ,λ γ , with λ γ , λ γ = 0, ±1, the helicity amplitudes squared for the process γ * g → QQ 2S+1 L J , the following relations hold (omitting numerical prefactors) Furthermore, in terms of the CO LDMEs we obtain and We conclude this section by noticing that each of the four independent azimuthal modulations in the cross section for e p → e J/ψ X, that is cos 2φ T , sin(φ S − φ T ), sin(φ S + φ T ) and sin(φ S − 3φ T ), probe a different gluon TMD. These modulations are the same as for the process e p → e Q Q X [4], after integration over the azimuthal angle φ ⊥ . As already pointed out in Ref. [4], such angular structures and the corresponding TMDs are very similar to the quark asymmetries in the SIDIS process e p → e h X, where the role of φ T is played by φ h [62].

IV. AZIMUTHAL ASYMMETRIES
In order to single out the different azimuthal modulations of the cross section dσ, given in Eq. (19) and Eqs. (21)-(22), we define the following azimuthal moments where the denominator reads with N and A U given by Eqs. (23) and (24), respectively. By taking W = cos 2φ T we obtain Moreover, assuming |S T | = 1, the other moments can be written as where the explicit expressions for (29) can be further simplified, if one employs the heavy-quark spin symmetry relations [34] 0|O J/ψ 8 Hence, at leading order in v, we obtain: and The asymmetries in Eqs. (32), (34) and (35) vanish in the limit y → 1 when the virtual photon is longitudinally polarized. Moreover, very importantly, we point out that a measurement of the ratios would directly probe the relative magnitude of the different gluon TMDs, without any dependence on the color octet LDMEs. Notice that the validity of Eqs. (40)- (42) is quite general and is not based on the heavy-quark spin symmetry relations in Eq. (36). It would be interesting to check experimentally the behavior of these ratios of asymmetries because currently there are no reliable theoretical predictions. However, for the dipole gluon TMDs we expect, from model independent considerations [8], that the observable in Eq. (42) will reach the value one in the small-x limit. It remains to be seen if this holds also for the WW gluon distributions discussed in this paper.

V. QUARKONIUM POLARIZATION
The study of J/ψ polarization is often considered as a test of NRQCD. Hence, in this section we present the cross sections for the processes e p → e Q L/T X, where the quarkonium in the final state is polarized either longitudinally (L) or transversely (T ) with respect to the direction of its three-momentum in the photon-proton center-of-mass frame. The cross sections have the same angular structure as in Eq. (19) and Eqs. (21)- (22). Namely, in terms of the kinematic variables defined in the previous section, where the superscript P = L or T denotes the polarization of the quarkonium and the superscripts U and T refer to the possible polarization states of the initial proton. Clearly, dσ = dσ L + dσ T , where dσ is the cross section for unpolarized quarkonium production given in Eq. (19). Furthermore, and with N defined in Eq. (23) and The explicit expressions for longitudinally polarized quarkonium production read in agreement with the results in Ref. [50], while is new. For completeness, the results corresponding to transverse polarization of the quarkonium read and More details on the derivation of the above cross sections, performed along the lines of Refs. [63][64][65], can be found at the end of Appendix A. We note that, also for polarized quarkonium production, it is possible to define azimuthal moments exactly as in Eq. (30), as well as their ratios in Eqs. (40)- (42). In particular, it turns out that such ratios of asymmetries depend neither on the LDMEs, nor on the polarization state of the detected quarkonium.

VI. A STRATEGY FOR THE DETERMINATION OF THE DOMINANT CO LDMEs
In this section we define novel observables that are only sensitive to the CO LDMEs 0|O J/ψ 8 ( 1 S 0 )|0 and 0|O J/ψ 8 ( 3 P 0 )|0 , and to the corresponding ones for Υ production, but not to TMDs. This is possible by combining azimuthal asymmetries in e p → e J/ψ (Υ) X with analogous quantities for open heavy-quark pair production in e p → e Q Q X [4,22,23] in the following way, where dσ Q now denotes the differential cross section for the process e p → e Q X defined in Eq. (19), and is the differential cross section for the process in which the quark-antiquark pair is almost back-to-back in the plane orthogonal to the direction of the proton and the exchanged virtual photon. Hence, in the γ * p center-of-mass frame, the difference of the transverse momenta of the outgoing quark and antiquark, conventionally specified by K ⊥ = (K Q⊥ − K Q⊥ )/2, should be large compared to their sum q T = K Q⊥ + K Q⊥ . In Eqs. (54)-(55), φ S , φ T and φ ⊥ are the azimuthal angles of the proton polarization vector S, q T and K ⊥ , respectively. Furthermore, y is the inelasticity and x B is the Bjorken variable, while z = K Q · P/q · P , with q = − as for e p → e Q X. We note that in the definitions of R cos 2φ and R the proton does not need to be polarized. The hard scale of the process e p → e Q X is identified with the quarkonium mass M Q ≈ 2M Q . Moreover, to avoid the presence of two very different hard scales in the calculation of the numerators of R cos 2φ and R, we simply take the photon virtuality Q to be Q = M Q ≈ 2M Q . Therefore, from the results for the cross section presented above, we obtain where we have introduced the shorthand notation Since we would like to have an exact cancellation of the gluons TMDs in the ratio, we need to consider the same value of the photon virtuality Q in both processes. Furthermore, the other hard scale K ⊥ ≡ |K ⊥ | in e p → e Q Q X is taken to be K ⊥ = Q to avoid any possible TMD evolution effect. From the results in Ref. [4] calculated at K ⊥ = Q = 2M Q and z = 1/2, we get   By taking the ratio of the two cos 2φ T -weighted cross sections in Eqs. (58) and (60), and the ratio of the two cross sections in Eqs. (59) and (61), it turns out that the two independent observables give access to two different combinations of the CO LDMEs O S 8 and O P 8 . Therefore the knowledge of both of them will allow to single out these two LDMEs. A similar comparison could be done by using jet pair production instead, but there quark contributions may spoil the cancellation of the gluon TMDs to some extent, depending on the kinematic region under study. This may be an option worth considering, should open heavy-quark pair production turn out not to be feasible at an EIC. In a recent Monte Carlo analysis for a future EIC it is shown that at least single spin asymmetries will be hard to study in open heavy-quark pair production (assuming 10 fb −1 ) [26], but here we propose a comparison to open heavy-quark pair production in the unpolarized proton case.
We point out that the measurement of quarkonium polarization will probe other combinations of the CO LDMEs and can be used for consistency checks. By extending the definitions in Eqs. (54) and (55) to longitudinally polarized quarkonium production, using the cross sections presented Section V, we find while, for transversely polarized quarkonium production, with R cos 2φ In particular, we notice that the measurement of R cos 2φ T , for values of the photon virtuality such that Q = M Q , will directly probe the matrix element The labels SV and CMSWZ denote the adopted LDME sets, given in Tables I and II. The bounds obtained with the BCKL set, not shown explicitly, lie always between the SV and CMSWZ results presented in the left panel.

A. Upper bounds of the asymmetries
The polarized gluon TMDs have to satisfy the following, model independent, positivity bounds [6] which can be used to calculate the upper limits of the azimuthal moments defined in the previous section. It can be easily seen that the Sivers asymmetry in Eq. (33) is bound to 1, while the asymmetries in Eqs. (34) and (35) have the same upper bound, which we denote by A W N . This is also the same upper bound of the weighted cross section cos 2φ T defined in Eq. (32).
For the numerical estimate of our asymmetries, we use different sets of extractions of the CO LDMEs for the J/ψ [39][40][41][42] (see Table I), and one set for the Υ(1S) ( Table II), all of them obtained from fits to TEVATRON, RHIC and LHC data. Note that most of these results are obtained from next-to-leading order (NLO) analyses, except for the SV set in Ref. [41], which is based on a LO calculation like our asymmetries. The negative value of 0|O J/ψ 8 ( 3 P 0 )|0 from the BK set leads to negative LO unpolarized cross sections for certain values of Q 2 . For this reason, the results obtained using the BK parametrization are not shown explicitly. The mass of the J/ψ is taken to be 3.1 GeV, while the one for Υ is 9.5 GeV. We define the charm and bottom quark masses to be equal to half of the mass of the J/ψ and the Υ, respectively.
First, in Fig. 2 we plot the maximum values of the azimuthal asymmetries A W N computed for the processes e p → e J/ψ X (left panel) and e p → e Υ X (right panel), in which the J/ψ and Υ are unpolarized, as a function of y and for different values of Q 2 . The maximum values for A W N in the case of longitudinally and transversely polarized quarkonium production are presented in Figs. 3 and 4, respectively. In general, it turns out the asymmetries depend very strongly on the specific set of LDMEs adopted. As mentioned above, they always vanish in the limit y → 1, and reach their maximum when y → 0. Alternatively, in Fig. 5 we show the Q-dependence of these maxima for the choice y = 0.1 and only for unpolarized quarkonium production. The results for polarized quarkonia do not present significative differences and, for this reason, are not shown explicitly. , for e p → e J/ψ X (left panel) and e p → e Υ X (right panel), with the J/ψ and Υ mesons longitudinally polarized along their direction of motion in the γ * p center-of-mass frame, as a function of y and for different values of the scale Q 2 . The labels SV and CMSWZ denote the adopted LDME sets, given in Tables I and II. The bounds obtained with the BCKL set, not shown explicitly, lie always between the SV and CMSWZ results presented in the left panel. In the J/ψ production plots one can see that depending on the set of LDMEs used and the quarkonium polarization state, the obtained behavior of the asymmetry as a function of Q 2 can be quite different, in some cases it increases or decreases, or even first decreases and then increases again. The Q 2 behavior of the asymmetries can thus be a further tool to determine the CO LDMEs, or at least their relative magnitude.
Finally, we point out that more stringent bounds on gluon TMDs than the ones presented in Eqs. (70) can be obtained by comparison with experiments. The COMPASS Collaboration has reported preliminary results on the Sivers asymmetry for e p → e J/ψ X [66]. The obtained value A sin(φ S −φ T ) = −0.28 ± 0.18 in the kinematical region of validity of our results, z ≥ 0.95, points towards a negative gluon Sivers function at low x and Q 2 , with a size about 1/3 of its positivity bounds. An even smaller gluon Sivers function has been found from the analyses of available data on inclusive pion and heavy-quark production in proton-proton collisions at RHIC [67,68]. However, in those analyses, TMD factorization has been assumed for single-scale processes, even if not supported by a formal proof. This phenomenological approach is known as generalized parton model (GPM) and is able to successfully describe many features of several available data. It still remains to be seen whether the effective TMDs determined within the GPM differ from the ones extracted from TMD-factorizing processes. , for e p → e J/ψ X (left panel) and e p → e Υ X (right panel), as a function of Q for y = 0.1. The labels CMSWZ, SV and BCKL denote the adopted LDME sets, given in Tables I and II. B. cos 2φ asymmetries in the MV model In the small-x limit, the nonperturbative McLerran-Venugopalan (MV) model [55][56][57] allows to calculate the gluon distributions inside an unpolarized large nucleus or energetic proton. The analytical expressions for the unpolarized and linearly polarized Weizsäcker-Williams (WW) gluon distributions in this model are given by [13,14]: where S ⊥ is the transverse size of the nucleus or proton, Λ is an infrared cutoff such as Λ QCD , and where Q sg (r) is the saturation scale for gluons, which in the MV model depends logarithmically on the dipole size r, and in general is a function of x. Factorizing the dependence on r as follows: Q 2 sg (r) = Q 2 sg0 ln (1/r 2 Λ 2 ), we can take Q 2 with Q 2 s0 = 0.35 GeV 2 at x = x 0 = 10 −2 from the fits to HERA data [69]. Adding e as a regulator for numerical convergence, in accordance with Ref. [4], we obtain for the ratio of both gluon TMDs: which is shown as a solid black line in the left panel of Fig. 6. The nonlinear evolution in rapidity of the gluon density in the presence of saturation is governed by the JIMWLK equation [58], which can be solved numerically [70][71][72]. In Refs. [18,59], an implementation of JIMWLK on a twodimensional lattice with spacing a was used to evolve the gluon TMDs f g 1 and h ⊥g 1 towards larger rapidities (or lower values of x), starting from an initial condition which corresponds to the MV model, using the prescription of Ref. [71]. We assigned a physical value to the lattice spacing a according to the relation aQ s0 = 0.015, which was obtained in Ref. [59] by studying the universal large-p T behavior of the gluon TMDs. For our choice of the saturation scale, this relation yields a = 0.025 GeV −1 . The numerical JIMWLK evolution is performed in steps δs = (α s /π 2 )δy = 10 −4 , with y = ln (x 0 /x) the rapidity. We show in Fig. 6 the initial condition at y = 0 for the ratio , as well as the result after 500 and 1000 steps in the evolution, which at a coupling α s (M J/ψ ) 0.2 corresponds to the values x 10 −3 and x 10 −4 , respectively.
In the same Fig. 6, we also compare with the analytical expression of the initial condition in Eq. (73). It can clearly be seen that the latter differs somewhat from its implementation on the lattice. Indeed, the uncertainty in the choice of the saturation scale and lattice size aside, the major source of this discrepancy is the way in which the infrared (IR) [GeV]   Table I. is regulated. On the lattice, this is taken care of by the finite lattice spacing, while in Eq. (73) we added a regulator by hand. To illustrate the freedom in the way the IR can be regulated, and the ensuing theoretical uncertainty for the ratio of the TMDs, we plot this ratio in the right panel of Fig. 6 for some different choices of the regulator.
With this ratio at hand, we show predictions for the cos 2φ T asymmetry, Eq. (32), as a function of q T in Fig. 7 and as a function of y in Fig. 8, both for the J/ψ and for the Υ mesons. Note that in Fig. 8 we took the numerical value for T ) corresponding to p T = 1.5 GeV from the numerical results in Fig. 6, and then plotted the analytical y-dependence with the help of Eq. (32). We restrict the range in p T to the region of validity of TMD factorization, which we estimate as p T = q T ∈ [0, M Q /2]. Note that in Fig. 7, the effects of the lattice discretization, which become more important for small values of q T , are clearly visible, even suggesting a rise of the ratio towards q T → 0, which is clearly an artifact of the discretization.

VIII. SUMMARY AND CONCLUSIONS
We have presented expressions for the azimuthal asymmetries for J/ψ and Υ production in DIS processes valid at LO in NRQCD, when the quarkonium in the final state is produced with a transverse momentum smaller than its invariant mass. Such observables can be used to extract information on the so far poorly known gluon TMDs, as well as to better understand the mechanism underlying quarkonium production by directly probing the two color-octet The labels SV and CMSWZ denote the adopted LDME sets, given in Tables I  and II. matrix elements 0|O We proposed ratios of asymmetries (Eqs. (40)- (42)) in which the LDMEs cancel out, but also ratios (Eqs. (54)- (55) and their polarized quarkonium analogues) in which the TMDs cancel out. The latter offer novel ways to extract the two mentioned CO LDMEs, which are still poorly known. This method requires one to consider comparisons of the processes e p → e Q X and e p → e Q Q X at the same hard scale, in order to establish a cancellation of the TMDs and to avoid having to include TMD evolution, although that in principle can be done with known methods as well. Since the proposed asymmetries are always given by the ratio of two cross sections, they have the further advantage of being less sensitive to the normalizations of the cross sections, to the effects of higher order corrections and to other sources of uncertainties, like the exact value of the charm and bottom mass. Moreover, as discussed in Ref. [50] to which we refer for details, we mention that for inclusive quarkonium production in DIS, nonperturbative effects like higher twist corrections and diffractive background can be more effectively suppressed as compared to quarkonium photoproduction, by looking at events with sufficiently high values of the photon virtuality Q 2 .
In our phenomenological analysis, we have used the positivity bounds for the gluon TMDs in order to estimate the maximal values of the asymmetries, and identify the kinematical regions in which they are sizable and where they could possibly be measured in future experiments at an EIC. For the gluon Sivers function this kind of studies provides little guidance because the corresponding maximal asymmetry is always 100%. A more stringent constraint comes from preliminary COMPASS results which seem to suggest a nonzero Sivers effect. An investigation by the COMPASS Collaboration of the other asymmetries as well, in a kinematical region complementary to that of the EIC, would be of course very beneficial. Finding nonzero single spin asymmetries in quarkonium production would be a sign of the CO mechanism, because in the CS mechanism they have to vanish due to the absence of any initial or final state interactions [33].
In addition to the upper bounds of the gluon TMDs, we have studied also their behavior at small values of x, which is expected to be probed at an EIC if Q 2 is not high, but of the order of Q 2 10 GeV 2 (also the higher √ s the better). The TMDs accessible in this process correspond to the so-called WW distributions at small x. This means that the T -odd distributions f ⊥ g 1T , h g 1T and h ⊥ g 1T are suppressed by a factor of x with respect to the unpolarized gluon density and, as such, cannot be described by current saturation models. On the other hand, the T -even distribution of linearly polarized gluons inside an unpolarized proton, h ⊥ g 1 , is not suppressed. We have therefore used the MV saturation model to show that it can give rise to sizable cos 2φ T modulations. Notably, we have explicitly checked that these modulations are only slightly reduced by small-x evolution effects and could be in principle measured even down to x = 10 −4 . The scattering amplitudes for the partonic processes γ * (q) + g(p) → QQ[ 1 S 1 ](P Q ) contributing to e p → e J/ψ (Υ) X can be written in the form [73] A µν [ 1 S (8) 0 ](q, p) = 1 4 1 ](q, p) = 1 4 where C S,J = 1 2J + 1 0|O M Q ≈ 2M Q is the mass of the quarkonium, and ε Sz is the polarization vector for the spin-1 quarkonium wave function. The operator O(q, p; k, P Q ) is calculated at leading order in perturbative QCD from the Feynman diagrams in Fig. 9, with k being half the relative momentum of the outgoing quark-antiquark pair, including the SU (3) color-octet projector 3i,3j|8a = √ 2 t a ij .
We obtain O µν (0) = − √ 2 δ ab ee c g s 2(M 2 Q + Q 2 ) where we have used the shorthand notation O(0) ≡ O(q, p; 0, P Q ). Similarly, the amplitudes for the subprocesses γ * (q) + g(p) → QQ[ 1 P J ](P Q ), in which a P -wave bound state is formed, read LzÔ µν A µν [ 3 P A µν [ 3 P 2 ](q, p) = i 3 C P,2 8 N c M Q ε αβ Jz (P Q ) Tr γ αÔ µν β (0) where we have definedÔ ( 1 S 0 )|0 in Eq. (48), which corresponds to a quarkonium that is longitudinally polarized, i.e. with helicity λ = 0, with respect to the one in Eq. (27), which corresponds to an unpolarized quarkonium. For QQ pairs in P -wave intermediate states, with L = S = 1, the method described above for the calculation of the unpolarized cross sections, consisting in the projection of the hard scattering amplitudes onto states of definite quantum numbers J and J z , are not useful when the final quarkonium is polarized. Instead, one can project the amplitudes onto states of definite L z and λ ≡ S z , square them and then sum over L z . The final results for the longitudinal and transversely polarized cross sections are obtained by using, respectively, the following relations for the polarization vectors ε λ (P Q ) of the quarkonium [76], ε α 0 (P Q ) ε * β 0 (P Q ) = P α Q P β Q M 2 Q − P α Q n β + P β Q n α P Q · n + M 2 Q n α n β (P Q · n) 2 , (A20) λ=±1 ε α λ (P Q ) ε * β λ (P Q ) = −g αβ + P α Q n β + P β Q n α P Q · n − M 2 Q n α n β (P Q · n) 2 , where n is any four-vector such that n 2 = 0 and P · n = 0. Obviously, by summing Eqs. (A20) and (A21), we obtain the second relation in Eq. (A18) with J z = λ.