TMD Evolution Study of the $\cos 2 \phi$ Azimuthal Asymmetry in Unpolarized $J/\psi$ Production at EIC

Semi-inclusive $J/\psi$ production in electron-proton collisions is a promising process to study gluon transverse momentum distributions (TMDs) at the future Electron-Ion Collider. In this article, we improve on previous studies of the $\cos 2 \phi$ azimuthal asymmetry that arises from the linear polarization of gluons inside unpolarized protons by including TMD evolution. We find that in the TMD regime the asymmetry grows monotonically with increasing transverse momentum of the outgoing $J/\psi$, in contrast to tree level calculations with Gaussian TMDs. Our predictions for the asymmetry at EIC can become very large at larger $x$, $Q$, and transverse momenta, even larger than the positivity bound. This problem stems from the very small $b$ region and implies a range of validity of TMD factorization that is more restricted than usually expected. We also include an estimate of the nonperturbative uncertainty from the large $b$ region and we conclude that it is smaller than the largest source of uncertainty, which stems from the choice of Color Octet Long-Distance Matrix Elements.


I. INTRODUCTION
Transverse momentum dependent parton distributions functions (TMDs) describe the transverse momentum distribution of quarks or gluons inside protons, or hadrons in general. For small longitudinal momentum fraction x the gluons dominate, but hardly anything is known yet about the gluon TMDs experimentally. Many theoretical proposals have been put forward that could potentially be used for extractions of gluon TMDs by investigating transverse momentum spectra and azimuthal asymmetries for bound and open heavy quark production, both in lepton-proton and in proton-proton collisions. This is because heavy quarks are very sensitive to the gluon content of hadrons, as they are predominantly produced from gluons and not intrinsically present in hadrons at small momentum fractions. Furthermore, some quarkonium states, like the J/ψ, are relatively straightforward to detect and numerous events can be collected. Therefore, quarkonium production can been considered as a main tool to extract gluon TMDs and many studies have appeared about this topic already [1][2][3][4][5][6][7][8][9][10][11][12][13].
One aspect of gluon TMDs that is especially of interest is that noncollinear gluons inside an unpolarized proton can be linearly polarized [14]. This manifest itself in scattering processes through azimuthal asymmetries as pointed out in [1,15,16]. These are the QCD analogues of the QED asymmetries that were recently observed in ultra-peripheral heavy ion collisions at RHIC, where linearly polarized photons from the colliding gold nuclei produce electron-positron pairs with cos 2φ and cos 4φ azimuthal asymmetries [17]. In order to measure these effects for gluons, the future U.S.-based Electron-Ion Collider (EIC) is well suited, especially using heavy quarkonium production processes. In semi-inclusive electroproduction of a heavy quarkonium state, such as a J/ψ vector meson, the linear gluon polarization manifests itself through a cos 2φ azimuthal asymmetry [7,10,11]. This observable is the main subject of this paper. Other promising processes are single or double quarkonium production at the LHC, where the linear gluon polarization manifests itself through cos 2φ and cos 4φ azimuthal asymmetries [12,18]. The advantage of the EIC is that only one TMD is involved, whereas in proton-proton collisions always a convolution of two TMDs is probed.
Semi-inclusive electro-production of a heavy quarkonium state at small transverse momentum is expected to be describable within TMD factorization, just like semi-inclusive deep inelastic scattering (SIDIS) for light hadron production and the Drell-Yan (DY) process [19][20][21]. This is thanks to the presence of a large energy scale, the photon virtuality Q, which allows to factorize the cross section description in a perturbative short-distance part that can be expanded in orders of the strong coupling constant α s and a nonperturbative long-distance part that is expressed in terms of TMDs. TMDs have to be modelled, calculated using lattice QCD, or extracted from experimental data. The TMD factorization description allows to incorporate the scale evolution, commonly referred to as TMD evolution. This in turn allows to improve on predictions and subsequent extractions.
In this article, we include TMD evolution effects at leading order in α s in a similar way as discussed in [10,22]. In this way, we obtain more realistic estimates for the transverse momentum spectrum and the associated azimuthal asymmetry for semi-inclusive J/ψ production in unpolarized electron-proton collisions for EIC kinematics. The results depend on the perturbative Sudakov factor, which has not yet been derived for this specific process. The one of [23] obtained from a matching calculation seems applicable only to light hadron production, as it has been demonstrated for pp → η c X [24] and for open heavy quark pair production in ep collisions [25] that there are no double logarithms associated to heavy quark production. Here we will implement the result of [24], which leads to considerably larger asymmetry values than when following [23]. There will also be considerable dependence on how the region of very small b (∼ 1/Q) is treated, indicating the need for matching to the collinear expression earlier than usually expected, i.e. at transverse momentum values below Q/2. All this will be addressed in detail.
In this study there is significant uncertainty from the nonperturbative contributions. Although no model for the gluon TMDs will be assumed, there is uncertainty from the nonperturbative Sudakov factor in the TMD evolution expressions, as well as from the nonperturbative formation of the bound quarkonium state from a produced heavy quark pair. The latter is usually described within nonrelativistic QCD (NRQCD) [26]. The NRQCD framework involves another factorization: a separation of the perturbative short-distance contributions (expanded in α s ) from the nonpertubative Long-Distance Matrix Elements (LDMEs). The relative importance of the various LDMEs is estimated by means of the heavy quarkantiquark relative velocity v in the bound-state rest frame in the limit v 1, by the velocity scaling rules. Typically, charmonium states have v 2 ∼ 0.3, whereas bottomonium states have v 2 ∼ 0.1, which are considered small enough to expand in.
In TMD factorization the formation of the bound quarkonium state from a produced heavy quark pair has to be incorporated by introducing Shape Functions (SFs) [27,28] that can be viewed as describing the transverse momentum smearing that arises in the final state hadronization process. It also describes the formation of a colorless hadronic state by emission of soft-gluon radiation. The SFs are related to the LDMEs of NRQCD, although the relation is not yet known in detail (the relation for large transverse momenta can be obtained from matching in the way described in [23]). In this paper we will only consider the leading order relation, which means that we only consider Color Octet (CO) SFs and LDMEs. In the process ep → e [QQ]X at leading order (LO) in α s the heavy quark-antiquark will be in a CO state and CO LDMEs will dominate. Beyond LO also Color Singlet (CS) LDMEs will contribute, but as we will see, the uncertainty in the prediction from the uncertainty of the CO LDMEs is too large to necessitate inclusion of higher order corrections. As mentioned, another source of uncertainty is from the nonperturbative part of the TMD evolution expressions, for which we will estimate a reasonable range. Given these uncertainties, we can draw only qualitative conclusions on the dependence of the asymmetry on transverse momentum and how that differs from earlier results in the literature [7,10,11]. Nevertheless, the results are promising regarding future measurements of the asymmetry at the EIC, as it could become quite large.
The paper is organized as follows. In Section II, we briefly review the TMD description of the process and introduce the SFs. In Section III, we explain the TMD evolution formalism that is employed for the numerical evaluations. Especially, we discuss the explicit form of the nonperturbative in more detail. In Section IV we present our results on the azimuthal cos 2φ asymmetry, and finally conclude and discuss our findings in Section V.

A. Leading order TMD description of the process
We study the semi-inclusive process and we consider all particles to be unpolarized. The kinematic variables for this process are given by with q = l − l and q 2 = −Q 2 . The partonic subprocess contributing to this reaction at leading order J . The spectroscopic notation denotes that the produced QQ pair will form a bound quarkonium state with spin S, orbital angular momentum L and total angular momentum J. This transition from the QQ pair to a bound quarkonium state is commonly described in terms of NRQCD LDMEs. At this LO only the CO configuration, indicated with the superscript (8), will contribute. At order α 2 s also CS states contribute. Although the CS contribution is suppressed relatively to the CO by a perturbative coefficient of the order α s /π, according to the NRQCD scaling rules the relevant CO LDMEs 0|O |0 by a factor of order v 3 and v 4 , respectively [29]. Taken together, the CO configuration is enhanced with respect to the CS by a factor v 3 π/α s which for the smallest Q value that we consider in this paper (Q = 3 GeV) is around 2 and grows with increasing Q. Furthermore, it is known that at large z, the CS term becomes negligible [30], and z is fixed to 1 in our analysis. For these reasons we solely consider CO contributions.
The reference frame for this process is chosen such that both the virtual photon from the electron and the incoming proton move along theẑ-axis. The azimuthal angle φ T of the quarkonium transverse momentum is defined with respect to the lepton scattering plane, i.e. φ l = φ l = 0. Moreover, the scattered electron has a scattering angle θ defined in this plane. In Figure 1 a schematical setup of the reaction is shown. Because z = 1, where M denotes the J/ψ mass. In TMD factorization the differential cross section is schematically written as: where p 2 T = −p 2 T denotes the transverse component of the gluon momentum, L the leptonic tensor, Γ g the gluon correlator and A the partonic scattering amplitude. Subsequently, the correlator is parametrized at leading order in terms of two TMDs [14], with f g 1 (x, p T ), the unpolarized gluon distribution and h ⊥g 1 (x, p T ), the linearly polarized gluon distribution. To ensure we can apply TMD factorization, only the kinematic region is considered in which the transverse momentum P T of the J/ψ is small compared to the virtuality of the photon Q. Usually P T < Q/2 is considered as the range of validity for the TMD region, but we will see that this is too optimistic in the present case.
Without taking into account smearing effects in the transition from the QQ pair to a bound quarkonium state, the cross section is written in terms of TMDs and NRQCD LDMEs as (see e.g. [10], where q T ≡ P T ) with the normalization factor where e Q is the fractional electric charge of the quark and M h the mass of the proton. The explicit forms of the prefactors are where y is the inelasticity and the superscripts U + L, L and T refer to the specific polarization of the photon [1,31]. Employing heavy-quark spin symmetry relations [26], one finds in terms of CO LDMEs [10]: Taking into account evolution in the above result for the cross section requires consideration of the more general TMD factorization expression including SFs.

B. Smearing effects and shape functions
As mentioned, the TMD factorized expressions have to take into account final state smearing effects that are encoded in the SF ∆ [n] [27,28]. This nonperturbative hadronic quantity describes the transition from the QQ pair to a bound quarkonium state, not only the formation of the bound state, but also the soft-gluon radiation required to produce a final state hadron in the CS state, which generally will change the momentum of the quarkonium. Including the SFs in the cross section becomes where we have introduced the transverse momentum convolutions: with the transverse momentum dependent weight function The SF can be thought of as a generalization of the LDMEs in collinear factorization. It is expected that the shape functions are proportional to the LDMEs, also beyond LO: for some universal SFs ∆(k 2 T ) and ∆ h (k 2 T ) that could in principle be unequal. In absence of smearing ∆(k 2 T ) = ∆ h (k 2 T ) = δ 2 (k T ), the convolutions reduce to products of an LDME with a gluon TMD and Eqs. (5)-(10) are recovered. This is the simplification we will adopt, but only after including effects from TMD evolution.
The process under investigation has an azimuthal asymmetry Adopting Eq. (15), this becomes where A and B are given in Eqs. (6) and (7), and .
Before continuing, we make a few remarks on the prefactor B/A. It turns out that B/A depends very strongly on the specific set of LDMEs adopted. Also, since | cos 2φ T | ≤ 1 and |R| ≤ 1, the LDMEs must be such that |B/A| ≤ 1 and A ≥ 0, which are important constraints to impose on LDME extractions at EIC. These constraints are not satisfied by the BK set for Q 2 < ∼ 2.5M 2 , but that seems a problem of applying LDMEs that were obtained from an NLO analysis in a LO analysis. Finally, we note that B/A vanishes in the limit y → 1 when the virtual photon is longitudinally polarized and maximizes when y → 0.

III. TMD EVOLUTION IMPLEMENTATION
Beyond tree level, the TMDs and SFs, become scale dependent which governs the TMD evolution [32]. Implementing TMD evolution is more easily done in impact parameter space, where convolutions become simple products. In general, we can write where b T is the conjugate of momentum q T .Ŵ consist of three factors: Here H denotes the hard part,Â andB the Fourier transformed TMD and SF, respectively, µ the renormalization scale and ζ A/B the rapidity variables. The latter arises due to the required regularization of lightcone divergences from lightlike Wilson lines [32], although in this case there are no such divergences associated to the SF [24,27]. The natural choice to minimize large logarithms in ζ A/B will be ζ A = Q 2 , ζ B = 1, similar as in [33] for open heavy quark pair production in electron-proton collisions.
The Fourier transforms of f g 1 , h ⊥g 1 and ∆ [n] (h) are defined as follows: In terms of these functions the convolutions can be written as: where we suppressed the dependence on ζ and µ.
The TMDs obey the Collins-Soper and Renormalization Group equations with respect to the scale parameters ζ and µ [32]. These equations can be used to evolve the TMDs from high to low scales [32,34]: with Sudakov factor S A : While the renormalization scale µ in the hard scattering part H should be set to µ ∼ Q to avoid large logarithms of µ/Q, the TMDs should be evaluated at much lower scale in order to avoid large logarithms of µ/M h or µ/Λ QCD . Instead of selecting a fixed, low (but still perturbative) scale for the TMDs, it is common to take The Sudakov factor then expresses the resummation of logarithms in µ b /Q.
In leading order in α s the perturbative expression for the above Sudakov factor can be written as [32,34] Adopting the scale µ ∼ Q in the SF is also expected to lead to large logarithms in Q/µ b , which should also be resummed. This can again be done using a Renormalization Group equation, leading to a contribution to the overall Sudakov factor at the single log level. It can be incorporated by changing S A into: The constant B CO from the soft-gluon radiation is not derived for the current process yet. However, since ∆ (h) (k 2 T ) in Eq. (15) is expected to be universal, we take for the present case of CO J/ψ production B CO = 1 as obtained in pp → η c X for CO η c production [24]. This is also consistent with recent results presented in [35] obtained within the SCET formalism.
Including the one-loop running of α s , one can preform the µ integral explicitly: Note that S A is spin independent, and thus the same for all convolutions.
We note that the perturbative Sudakov factor for our process is quite different from the one of [23], where it is assumed that the SF also leads to a factor S A given in Eq. (28), like for TMD fragmentation functions. This leads to a larger Sudakov factor and as a consequence to more Sudakov suppression. We then find smaller values for the asymmetry, but the shape of the asymmetry is quite similar.
The above perturbative expression for the Sudakov factor is valid in the region b 0 /Q ≤ b T ≤ b T,max . The lower limit is the point beyond which µ b becomes larger than Q, such that the Sudakov integral flips sign. The upper limit marks the point where perturbation theory starts to fail, which is not exactly known. It is common to take b T,max = 0.5 GeV −1 or b T,max = 1.5 GeV −1 in phenomenological analyses.
The two limits to the b T integration are implemented in different ways. One way to ensure b 0 /Q ≤ b T is to consider the replacement [22] in the Sudakov factor, which effectively boils down to a different resummation: in logarithms of µ b /Q rather than µ b /Q. For consistency this is then also the scale one should use in the TMDs and SFs. In this way the perturbative expression for S A is valid for all b T ≤ b T,max . Often a slightly different replacement is considered [36] However, we will use the µ b replacement for our predictions. We will comment on this choice below.
Different ways to ensure that b T ≤ b T,max in the perturbative expression have been employed, but the most common one is the b * T -method [37]: This is introduced in the following way:Ŵ where forŴ (b * T ) one can always use the perturbative expression for the Sudakov factor and the nonperturbative Sudakov factor S N P makes up for the difference to the realŴ (b T ). This factor can only be extracted from data and is in principle different for different convolutions. In addition, S N P depends on the prescriptions used to separate the perturbative and nonperturbative components inside the convolution. There are different parameterizations used in the literature, but typically it is chosen to be a Gaussian. In general it is Q dependent and of the form [37]: where Q N P is a parameter with the dimension of mass, that typically is chosen to be (near) the smallest scale at which perturbation theory is expected to be valid; of course, any change in Q N P can be compensated for by a change in the other terms. The specific form of S N P will be discussed in the next section. General constraints are that exp (−S N P ) should be unity at b T = 0, and it should smoothly vanish at large b T , in order to exclude contributions from (far) outside the proton. It also guarantees convergence of the convolutions.
Implementing both (31) and (33) should be done in the right order, i.e.
or similarly forμ b . In this way one ensures that S A (0) = 0. This is also the scale that one should adopt in the TMDs and SFs.
If we take all the above into account, the convolutions read where we suppressed the dependence on the scale √ ζ = µ = µ b * in the TMDs and SFs. At this perturbative scale, the b T dependence of the TMDs and SFs can be calculated. For the TMDs these "perturbative tails" are [38]: One sees that both expressions are determined by the collinear distributions f i/P , but start at different orders in α s , the reason being that there is no collinear version of h ⊥g 1 .
The leading order perturbative transverse momentum tail of∆ [n] has been studied in [23], where it is observed that the tail of ∆

[n]
h would require a study at higher order in α s . However, the Fourier transformed functions, i.e. the perturbative small-b expressions of the SFs, both start at order α 0 s , in contrast to the TMDs discussed above. This order α 0 s contribution stems from the small transverse momentum part, irrespective of whether it is the simple ∆(k 2 T ) = ∆ h (k 2 T ) = δ 2 (k T ) expression or more realistic smeared-out versions of ∆(k 2 T ) and ∆ h (k 2 T ). Given the uncertainty in the α s corrections to the tails of the SFs, here we will simply adopt the LO description in terms of the LDMEs.
In this section we have discussed all perturbative ingredients to perform TMD evolution numerically at leading order in α s . In order to make numerical predictions we need to discuss the nonperturbative Sudakov factor in greater detail. As this factor is unknown for gluons, we need to estimate the uncertainty this introduces in the predictions.

A. Nonperturbative Sudakov factor and error estimation
A parameterization for S N P was obtained from fits to low energy SIDIS data as well as higher energy DY and Z boson production data [34] S N P (b T ; Q) = g 1 ln Q 2Q N P + g 2 1 + 2g 3 ln 10xx 0 with g 1 = 0.184 GeV 2 , g 2 = 0.201 GeV 2 , g 3 = −0.129 GeV 2 , x 0 = 0.009, Q N P = 1.6 GeV and b T,max = 1.5 GeV −1 . We take this expression as our starting point for two fixed small x values and Casimir scaled (i.e. multipied by a factor C A /C F ) in order to apply to gluons: The B-term can be thought of as related to the intrinsic transverse momentum of the TMDs, which generally is x dependent. By Fourier transforming a simple Gaussian dependence for f g Of course, the B-term is not necessarily the same for the h ⊥g 1 (x, p 2 T ) convolution, but for simplicity we do not distinguish between these cases. Matching Eqs. (41) and (42) for the x values that we will consider yields the specific B values shown in Table I.  To perform an error estimate and to assess the importance of S N P for the size of the convolutions and asymmetry, we will vary A within the extreme limits following a previous study [12], instead of taking the one quoted in Eq. (42). The idea is that one expects the exp(−S N P ) term to be non-negligible (here defined as being larger than 10 −3 ) anywhere between b T,max and the charge radius of the proton. If the exp(−S N P ) term becomes negligible around b T,max already, then there will be hardly any nonperturbative contribution outside the perturbative regime, which is not realistic. On the other hand, one does not expect significant contributions beyond the charge radius of the proton, offering a generous upper bound.
To implement this range, we define a value b T,lim such that at large Q, where the A-term is dominant and the B-term can be neglected, exp[−S N P ] becomes negligible, i.e. ∼ 10 −3 . Considering b T,lim as the diameter, since it is conjugate to q T , r is defined as the characteristic radius r = b T,lim /2, that can be thought of as the range over which the interactions occur from the center of the proton. For three values of b T,lim or r we determined the A values at Q = 12 GeV, which are shown in Table I We will include the B-term in S N P as it is needed for smaller Q values and large b T,lim values. This is illustrated in Figure 2: the product exp[−S A ] exp[−S N P ] at Q = 3 GeV receives large perturbative contributions from S A at large b T , all the way up to or even beyond the proton radius. This results in an upward bump at small q T in the convolutions, in particular in C[wh ⊥g 1 ∆

[n]
h ]. Moreover, at small Q the limit of B → 0 can result in a violation of the positivity bound for the convolutions: ]. For larger Q values, the curves with and without B-term lie increasingly closer to each other, as expected.

IV. NUMERICAL PREDICTIONS FOR THE AZIMUTHAL ASYMMETRY
For the numerical calculations we can use different sets of extractions of CO LDMEs for the J/ψ case. These values are obtained from fits to Tevatron, RHIC and LHC data and summarized in Table II. Most of these results are obtained from next-to-leading order analyses, except for the SV set, which is based on a leading order calculation. The mass of the J/ψ (M ) and the charm quark (m c ) are taken to be 3.1 GeV and 1.4 GeV, respectively.
The perturbative tails of the TMDs in terms of the collinear quark and gluon PDFs, are computed by using the MSTW2008LO PDF set [43]. Furthermore, we employ one-loop running of α s . To determine the QCD scale we matched α s to the PDF set: α s (M Z ) MSTW2008LO ⇒ Λ QCD = 0.255 GeV . For computations with Q = 3 GeV < m b ≈ 4 − 5 GeV, we use n f = 4 instead of 5.  [39] 8.9 ± 0.98 0.56 ± 0.21 SV [40] 1.8 ± 0.87 1.8 ± 0.87 BK [41] 4.5 ± 0.72 −0.54 ± 0.16 BCKL [42] 9.9 ± 2.2 0.49 ± 0.44 As an illustration of the typical features of the TMD evolution of the convolutions q T C[f ⊥g h ], and their ratio, in Figure 3 we show the results for one particular LDME, namely 0|O  We observe that the q T -spectrum broadens and the estimated uncertainty band from the unknown nonperturbative contributions becomes smaller with increasing Q, as one would expect. All curves are shown for q T < Q/4 in order to ensure that the positivity bound is respected. For most curves this restriction is sufficient, except for x = 10 −1 the restriction q T < Q/4 is not enough for the higher Q values. Therefore, we cut off at an even lower q T when making azimuthal asymmetry predictions for x = 10 −1 .
In our computations the ratio of convolutions R violates the positivity bound within what is usually expected to be the range of validity of TMD factorization q T < Q/2, as shown in Figure 4. This problem originates from the very small b region. This is clear from Figure 4, where we compare the results for the two ways to ensure that b T ≥ b 0 /Q in the perturbative Sudakov factor that we discussed earlier: µ b → µ b or µ b →μ b . The region where the asymmetry starts to become sensitive to how the very small b region is treated roughly corresponds to the region where the positivity bound is violated. As far as we know this is the first instance of an azimuthal asymmetry to display sensitivity to the very small b region well within what is commonly considered the TMD region. Since the replacement µ b → µ b leads to slightly smaller asymmetries, we adopt that choice in what follows, but for q T < Q/4 that choice does not really matter. Returning to the discussion of Figure 3, generally the maxima of the convolutions increase towards smaller Q, except for the small Q cases. For q T C[f g 1 ∆] this can be understood from the significant relative decrease in magnitude of the perturbative f g 1 tail for smaller Q values in combination with the Sudakov factors, where the B-term is then of large influence. The decrease of q T C[wh ⊥g 1 ∆ h ] is, on the other hand, predominantly due to the Sudakov factors as the perturbative h ⊥g 1 tail at these x values stays approximately constant with varying Q.
The differences in behavior of the two convolutions depends on various factors: the differences between the tails of the TMDs, the Sudakov factors and the type of Bessel function. The product of the tail with the Sudakov factors will go to zero at large b T and the presence of h ⊥g 1 in a convolution contributes to reducing the magnitude of the integrand and to b T -broadening. The h ⊥g 1 tail is naturally suppressed by order α s in comparison to the f g 1 tail. However, α s (µ b * ) is growing with b T , up to its upper value ≈ α s (b 0 /b T,max ), and the h ⊥g 1 tail can become larger in comparison with the f g 1 tail at large b T . This effect becomes more pronounced for smaller x.
The consequence of the b T -broadening is that more damped oscillations of the J 0 Bessel function in C[f g 1 ∆] occur before the integrand becomes zero. Each additional oscillation in the integrand brings the convolution closer to zero and more oscillations fit in a given b T -range when q T increases. Therefore, C[f g 1 ∆] with smaller Q decreases faster. The situation is different for C[wh ⊥g 1 ∆ h ] that contains the J 2 Bessel function. This function starts its damped oscillation at zero and goes up. The consequence is that the b T -integrals benefit from unsuppressed intermediate b T values. This results in a peak maximum before large b T oscillations will bring the convolution down toward zero in a similar way as for C[f g 1 ∆]. This point is at smaller q T for smaller Q taking into account the b T -broadening of the Sudakov factors. Another crucial difference is that the envelope of J 2 tends slower toward zero than the J 0 one with increasing b T . The consequence is that C[wh ⊥g 1 ∆ h ] falls slower than C[f g 1 ∆]. Hence, R and the azimuthal asymmetry, always grow with q T , but more slowly for larger Q. In addition, as the large b T values in C[wh ⊥g 1 ∆ h ] are less suppressed than in C[f g 1 ∆], the azimuthal asymmetry and C[wh ⊥g 1 ∆ h ] are more sensitive to the changes in S N P .
Varying x in the computations changes the convolutions. First, we notice that the overall magnitude of the convolutions is smaller for larger x, because the collinear gluon PDF is less prominent at larger x. On the other hand, the magnitude of R becomes higher for large Q, but lower for small Q when increasing x in the range of chosen values. Secondly, the shape of the perturbative TMD tails is different, in particular the f g 1 tail is broader in b T for larger x. On top of that, the B-term in S N P is smaller, as can be seen in Table I. Together, the broader b T integrands make these convolutions go faster to zero for larger x. This explains the behavior of the magnitude of R and that the azimuthal asymmetry become less straight at x = 10 −1 , especially visible for larger Q.
After these general qualitative observations, we present predictions for the azimuthal asymmetry at the EIC, shown in Figures 5 to 8 ( 3 P 0 )|0 from the BK set is negative which can lead to negative cross sections for certain values of Q when used in our LO expressions, and the values of BCKL are very similar to the CMSWZ values, but with larger errors. In general the asymmetry is (much) larger when using the CMSWZ set compared to the SV set, especially at small Q. TMD evolution predicts the azimuthal asymmetry to grow with q T as discussed previously, and we see that the choice of LDMEs is of large influence in the predictions. This results in an uncertainty in the predictions that is larger than that due to S N P . In the plots we restrict to the central values of the LDME fits, but of course, taking into account the uncertainties in these values would broaden the bands further, such that the bands for CMSWZ and SV start to overlap more, giving rise to a large range of possible asymmetry values at EIC. Asymmetries anywhere between 1% and 20% may thus be expected at the EIC, which seem feasible to measure. Improved constraints on the LDMEs, and more generally on the TMD shape functions, can very likely be obtained in this way.

V. DISCUSSION AND CONCLUSIONS
We have discussed the process of semi-inclusive electroproduction of J/ψ's. Although the initial proton and electron beams and the final state J/ψ are all unpolarized, there is an effect from the linear polarization of gluons inside unpolarized protons. It gives rise to a cos 2φ T azimuthal asymmetry for which we obtained predictions for EIC kinematics. We have included the effect of TMD evolution in order to obtain more realistic estimates. In contrast to a leading order Gaussian TMD model prediction [7], we find that the azimuthal asymmetry grows monotonically in the TMD regime, more in line with the predictions from the Generalized Parton Model approach with additional gluon radiation [44], except that the magnitude can be much larger by as much as an order of magnitude depending on the LDME set considered. Of course, the asymmetry cannot grow beyond 1, therefore, the expectation is that a maximum will be reached outside the TMD region. We found the latter region to be significantly smaller than usually expected. For this particular process and asymmetry our computations indicate that a bound of q T < Q/4 is required in order to respect the positivity bound and to not become sensitive to very small b values in the TMD region.
Larger and monotonically rising asymmetries have also been obtained in [10] for the McLerran-Venugopalan model including nonlinear evolution x, showing decreasing asymmetries with decreasing x values. In our results the x dependence of the asymmetries is less systematic in the EIC kinematic range, as it depends on the considered Q value and on the LDMEs. Especially SV at Q = 3 GeV does not follow the observed trends and is exceptionally small due to a cancellation of the S-and P -wave LDMEs in Eq. (10) at Q 2 near M 2 .
In our calculation we have restricted to leading order expressions for the small-b T expressions for the TMDs and the SFs. Assuming universality of the SFs we have adopted a perturbative Sudakov factor that was derived for pp → η c X [24], where there is no double logarithm associated to the SF of the heavy quarkonium state, only a single logarithm. The prefactor of the latter still needs to be confirmed by direct computation. Especially the absence of a double logarithmic contribution from the heavy quarkonium state reduces the amount of Sudakov suppression and brings the violation of the positivity bound down to smaller values of q T .
We have estimated the uncertainty from the nonperturbative Sudakov factor and found that the uncertainty in the CO LDMEs forms the dominant source of uncertainty. In some cases this resulted in an order of magnitude uncertainty in the predicted asymmetry. In [45] it is pointed out how one can use the polarization of the J/ψ or the comparison to open heavy quark pair production to reduce this uncertainty through measurements at the EIC, but given this uncertainty, going to the next order in α s in the TMD evolution calculation will not lead to much more precise predictions at the current stage. Nevertheless, the conclusion from our numerical study is that asymmetries are expected to be measurably large, especially at the larger center of mass energy of 140 GeV.