Azimuthal asymmetries in semi-inclusive $J/\psi\,+\,\mathrm{jet}$ production at an EIC

We consider transverse momentum dependent gluon distributions inside both unpolarized and transversely polarized protons and show how they can be probed by looking at azimuthal modulations in $e \, p \to e \, J/\psi \, \mathrm{jet} \, X$. We find that the contribution due to quark induced subprocesses is always suppressed in the considered kinematic regions, accessible in principle at a future Electron-Ion Collider. Our model-independent estimates of the maximal values of these asymmetries allowed by positivity bounds suggest the feasibility of their measurement. In addition, by adopting the McLerran-Venugopalan model for the unpolarized and linearly polarized gluon densities, we study the behavior of the $\cos2\phi$ asymmetries in the small-$x$ limit.


I. INTRODUCTION
Transverse momentum dependent distribution functions (TMDs) of gluons inside unpolarized and polarized protons have been defined in terms of QCD operators for the first time in Refs. [1,2]. Since then, they have received growing attention, mainly because they encode essential information on the transverse motion of gluons and their spin-orbit correlations. As such, they parameterize highly nontrivial features of the partonic structure of the proton. For instance, the distribution of linearly polarized gluons can be nonzero even if the parent proton is unpolarized and, if sizable, can modify the transverse spectra of (pseudo)scalar particles like, for example, the Higgs boson [3][4][5][6][7]. Another example is provided by the gluon Sivers function [8,9]. In general, the Sivers function describes the azimuthal distributions of unpolarized partons inside a proton that is transversely polarized with respect to its momentum; it is expected to generate observable single spin asymmetries in processes initiated by transversely polarized protons. Moreover, it can provide an indication on how much quarks and gluons contribute to the total spin of the proton through their orbital angular momentum.
Gluon TMDs are also of great interest because of their intrinsic process dependence, due to their gauge link structure. This is determined by soft gluon exchanges occurring in the particular processes in which the TMDs are probed. Unambiguous tests of such properties have been proposed for both the quark [10,11] and the gluon [12] Sivers distributions. Their verification will represent an important confirmation of our present knowledge of the TMD formalism and nonperturbative QCD in general.
Experimentally, not much is known about gluon TMDs. However, many proposals have been suggested to probe them, especially through the measurement of transverse momentum distributions and azimuthal asymmetries for heavy-quarkonium production, both in lepton-proton [13][14][15][16] and in proton-proton collisions [17][18][19][20][21][22]. In particular, in a very recent publication [16], some of us considered the inclusive production of a J/ψ or a Υ meson in deep-inelastic electron-proton scattering, namely e p → e J/ψ (Υ) X. The corresponding cross section can be calculated in the TMD formalism only when two well-separated scales are present: a soft one, sensitive to the intrinsic parton transverse momenta, and a hard one, which allows for a perturbative treatment. Hence the analyzed kinematic region was the one in which the transverse momentum q T of the J/ψ meson (the soft scale) is much smaller than its mass M (the hard scale), namely q T M . Such an analysis could be performed at the future Electron-Ion Collider (EIC) planned in the United States [23,24].
In the present paper, we study instead a somewhat complementary process: the associated production of a J/ψ or a Υ meson and a hadronic jet in deep-inelastic electron-proton scattering, e p → e J/ψ (Υ) jet X, where the initial proton can be either unpolarized or transversely polarized in the virtual photon-proton center of mass frame. In this case the soft scale is given by the total transverse momentum of the J/ψ + jet pair, required to be much smaller than its invariant mass, while the individual transverse momenta of the J/ψ and the jet need not to be small. The advantage is that, by varying the invariant mass of the pair, one can now access a range of scales, in contrast to the

II. OUTLINE OF THE CALCULATION
We consider the process depicted in Fig. 1, where the proton is polarized with polarization vector S, while we do not observe the polarization of the other particles. The kinematics is defined with the help of a Sudakov decomposition in terms of two lightlike vectors, here chosen to be the momentum P of the incoming proton, and a second vector n which fulfills the relations n · P = 1 and n 2 = 0. In a frame where the three-momenta of the incoming proton and the virtual photon exchanged in the reaction lie on theẑ-axis, one has the following expressions for the momenta of the incoming parton (p), virtual photon (q), outgoing jet (K j ) and quarkonium (K ψ ): In the above, x = p · n, M is the mass of the quarkonium, while the virtuality of the photon is defined as Q 2 = −q 2 , the Bjorken-x variable is given by x B = Q 2 /2P · q, and z = P · K ψ /P · q is the energy fraction transferred from the photon to the quarkonium in the proton rest frame. We clarify that, in our conventions, the subscript T refers to a soft transverse momentum, whereas ⊥ denotes the large transverse component of a measurable hadronic momentum. In this process, and for our specific choice of the reference frame, the two directions are the same, as can be seen from the two-dimensional delta function in Eq. (12) below. Moreover, the momentum of the incoming lepton reads: with y = P · q/P · the inelasticity variable, and from which the momentum of the scattered lepton can be obtained as = − q. Similarly, for the proton spin vector with S 2 L + S 2 T = 1 and M p denotes the proton mass. Assuming TMD factorization, the cross section can be written as [37]: where we have introduced the antisymmetric transverse projector µν T = µναβ P α n β , with 12 T = +1. In the symmetric part of the correlator, Γ {µν} T /2, three T -odd TMDs appear: the gluon Sivers function f ⊥ g 1T (x, p 2 T ), as well as the TMDs h ⊥g 1T and h g 1T which are chiral-even distributions of linearly polarized gluons inside a transversely polarized proton. In analogy with the transversity TMD for quarks, one introduces the combination which vanishes upon integration over transverse momentum [12], in contrast to its quark counterpart. Finally, the distribution g g 1T (x, p 2 T ) is T -even and represents the circular polarized gluon content of a transversely polarized proton. Note that the factor 1/2, necessary to average over the incoming gluon polarization, is already taken into account in Eqs. (7) and (8).
Integrating out the azimuthal angle of the final lepton [39], one has Moreover, we can write: and where we made use of the relation Q 2 = x B ys, with s = ( + P ) 2 the square of the center-of-mass energy. Finally, integrating overz, p T and x, the cross section in Eq. (5) can be rewritten as FIG. 2: Azimuthal angles for the process e p → e J/ψ jet X in a reference frame where φ = φ = 0. The vectors qT = K ψ⊥ + K j⊥ and K ⊥ = (K ψ⊥ − K j⊥ )/2 define the azimuthal angles φT and φ ⊥ , respectively.
In the above expression, The region of validity of our calculation is then the one where there are two strongly ordered scales: q T ≡ |q T | K ⊥ ≡ |K ⊥ |, corresponding to the setup where the J/ψ and the jet are produced almost back to back in the plane transverse to the collision axis. We can then set K ⊥ K ψ⊥ −K j⊥ in our calculation, and retain only those terms which are, at most, of the order of q T /M p , discarding the ones suppressed by powers of q T /M or q T /K ⊥ . Note that, integrating Eq. (13) over q T , one recovers the collinear result [25] dσ where we have introduced the usual Mandelstam variablesŝ andt for the subprocess γ * g → J/ψ g, and where x min = (Q 2 + M 2 )/(ys),t min = −(ŝ + Q 2 )(ŝ − M 2 )/ŝ.

III. ANGULAR STRUCTURE OF THE CROSS SECTION: ANALYTIC RESULTS
In the reference frame defined above, and measuring the azimuthal angles w.r.t. the lepton plane (φ = φ = 0) as depicted in Fig. 2, we denote by φ S , φ T and φ ⊥ the azimuthal angles of the three-vectors S T , q T and K ⊥ , respectively. In the region q T K ⊥ , we obtain that the differential cross section can be written in the form [37] dσ At leading order (LO) in perturbative QCD, we find The normalization factor N is given by where e Q is the electric charge of the heavy quark in units of the proton charge, and with the transverse mass M ⊥ defined as M ⊥ = M 2 + K 2 ⊥ . The amplitudes of the modulations A eg l , with l = 0, 1, 2, describe the interaction of an electron with an unpolarized gluon, through the exchange of a photon which can be in different polarization states [37], and can be expressed as where the subscripts i = U + L, L, I, T refer to the specific polarizations of the photon: unpolarized plus longitudinal, longitudinal, interference transverse-longitudinal and transverse, respectively [37,40]. Namely, denoting by A γ * g λγ ,λ γ , with λ γ , λ γ = 0, ±1, the helicity amplitudes squared for the process γ * g → QQ 2S+1 L (c) J g, the following relations are fulfiled where we have omitted numerical prefactors. Analogously, for the amplitudes B eg m with m = 0, 1, 2, 3, 4, one can write The amplitudes A γ * g and B γ * g have been calculated at LO in the framework of NRQCD, taking into account both the CS and CO contributions corresponding to the Feynman diagrams in Fig. 3.
According to the CS mechanism, the heavy-quark pair is produced in the hard process J ] g directly with the quantum numbers of the observed quarkonium, which in the case of the J/ψ meson are S = 1, L = 0 and J = 1. This means that the pair is in a 3 S (1) 1 state. On the other hand, the allowed CO states are 1 S J ] q should also be taken into account, in principle. However, since it turns out to be negligible in the kinematic region under study (see Sec. V), we will not consider it in the following.
In the calculation of A γ * g i and B γ * g i , all the CS and CO perturbative contributions need to be added, each one weighted by its corresponding LDME, denoted by 0|O c ( 2S+1 L J )|0 , encoding the soft hadronization of a QQ pair with angular momentum configuration 2S+1 L J and color assignment c = 1, 8, into a colorless state with the quantum numbers of the J/ψ. One can write with an analogous expression holding for the amplitudes B γ * g i .

FIG. 3: Feynman diagrams representing the partonic subprocesses underlying the reaction
The crossed diagrams of (a)-(e), obtained by reversing the heavy quark lines, are included in the calculations but not shown explicitly. Only diagrams (a)-(c) contribute to the CS production mechanism.
Making use of the heavy-quark spin symmetry relation [27] which is valid up to higher order corrections in the small expansion parameter v, and the fact that the color factors read Eq. (24) can be simplified further to yield where we defined A The explicit expressions for the hard scattering functions A γ * g i and B γ * g i we introduced above, depend on the production mechanism of the quarkonium we consider. For the CS hard parts, we obtain: and Since the CO expressions for the hard parts are very long and not very enlighting, we do not show them here. We note that our results for the angular independent part of the cross section fully agree with the ones published in Ref. [25], obtained within the framework of collinear NRQCD, while, to the best of our knowledge, the angular modulations are derived here for the first time.

IV. AZIMUTHAL MODULATIONS
The cross sections in Eqs. (17) and (18) exhibit various azimuthal modulations, and have the same structure as similar 2 → 2 processes such as γ * g → QQ [37]. These modulations can be exploited to extract ratios of specific gluon TMDs over the unpolarized one, f g 1 (x, q 2 T ). To this end, we define the following azimuthal moments where the denominator is given by As an example, to extract the ratio h ⊥ g is the distribution of linearly polarized gluons in an unpolarized proton, one could measure: or .
The gluon TMDs for a transversely polarized proton can be extracted in a similar fashion. In particular, the polarized cross section, Eq. (18), can be simplified by integrating over the angle φ ⊥ , yielding: where we have introduced the combination h g 1 of Eq. (9). Clearly, the resulting integrated cross section has only three independent azimuthal modulations left, each related to a different T-odd gluon TMD. Note that this situation is analogous to the case of quark azimuthal asymmetries in SIDIS (e p ↑ → e h X), in which the role of φ T is played by φ h [41].
Using Eq. (31), one then obtains (setting |S T | = 1) where we note that the asymmetries in Eqs. (36) and (37) vanish when y → 1, see the last line of Eq. (23), corresponding to a longitudinally polarized virtual photon. Furthermore, as already pointed out in Ref. [16] for the process e p → e J/ψ X, the following ratios of asymmetries are directly sensitive to the relative magnitude of the various gluon TMDs, without any dependence either on the different LDMEs or on the other kinematic variables in the process.

A. Use of positivity bounds
In order to assess the maximum allowed size of the above asymmetries, one can make use of the fact that the polarized gluon TMDs have to satisfy the following, model independent, positivity bounds [1]: Applying these to the asymmetries in Eqs. (32) and (33), we find that: where the maxima on the r.h.s. are always independent of the center-of-mass energy. As can be seen from Eqs. (36) and (37), the upper bounds of the asymmetries A sin(φ S +φ T ) and A sin(φ S −3φ T ) are equal to half the one for A cos 2φ T , while the bound for the Sivers asymmetry A sin(φ S −φ T ) , Eq. (35), is simply equal to one.

B. McLerran-Venugopalan model
As an alternative to the model-independent positivity bounds in the previous section, the gluon TMDs can be evaluated in a model, in order to show predictions for the azimuthal asymmetries in Eqs. (32) and (33   expression of x, which can be read off the delta function in Eq. (12), it follows that for realistic EIC center-of-mass energies √ s ∼ 100 GeV, and for values of y that are not too small, its value is of the order of x ∼ 10 −2 . In this regime, the nonperturbative McLerran-Venugopalan (MV) model [42][43][44] for the gluon distributions inside an unpolarized nucleus is phenomenologically very successful. In view of the total lack of information on gluon TMDs, following Ref. [12], we will apply it also to the proton case. The unpolarized and linearly polarized gluon TMDs in the MV model read [31,45]: In the above formulas, S ⊥ is the transverse size of the proton and Λ is an infrared cutoff (we take Λ = Λ QCD = 0.2 GeV). Furthermore, Q sg (r) is the gluon saturation scale, which we parameterize, as is usually done in the MV model, as Q 2 sg (r) = Q 2 sg0 ln (1/r 2 Λ 2 ). Finally, the gluon saturation scale is related to the one felt by the quarks by a color factor, Q 2 sg0 = (N c /C F )Q 2 sq0 , and we take the numerical value Q 2 sq0 = 0.35 GeV 2 at x = 10 −2 from the fit of Ref. [46].
The ratio of h ⊥g 1 over f g 1 is therefore: where e was added as a regulator to guarantee numerical convergence. In the following section, the above model is used to show the q T dependence of the asymmetries in Eqs. (32) and (33).

V. NUMERICAL RESULTS
Before presenting our results, we come back to the assumption of negligible quark contribution. To this aim, we compute the ratio of the collinear quark-over total unpolarized cross section in NRQCD, using for the parton distribution functions (PDFs) the MSTW2008LO set [49] and adopting the same LDME sets as for the study of the azimuthal asymmetries, see Table I (with the charm mass m c = 1.3 GeV) for J/ψ and Table II (with the bottom mass m b = 4.2 GeV) for Υ. Following Ref. [25], we choose the hard scale for these PDFs to be equal to ξ M 2 + Q 2 , where ξ varies between 1/2 and 2, and take a conservative estimate for the center-of-mass energy at the EIC, i.e. √ s = 65 GeV. As can be clearly seen in every plot we present in what follows, in the kinematic regions where the asymmetries can be sizable, the estimated quark contribution to the unpolarized cross section is always practically negligible, at least at our current level of accuracy.
To start, in Fig. 4 we show the upper bounds of the absolute values of the asymmetries A cos2φ T (left panel) and A cos2(φ T −φ ⊥ ) (right panel), for the J/ψ meson, as a function of the transverse momentum K ⊥ . We identified two regions in z and y where either of them is very large, corresponding to, respectively, z = 0.7 and y = 0.3, or z = 0.3 and y = 0.7. The J/ψ mass is taken to be M = 3.1 GeV, and the photon virtuality is fixed at Q 2 = 10 GeV 2 in order to keep all the large scales, i.e. K ⊥ , M and Q, in the same ballpark. The full lines correspond to the complete NRQCD calculation, while the CS results are shown as dashed lines. For these asymmetries, NRQCD consistently leads to larger upper bounds as compared to the CSM, and in certain kinematic regions even to a drastically different behavior  [47] and CMSWZ [48]. The ratio of the collinear quark contribution over the total unpolarized cross section is shown in the form of a band, representing the scale uncertainty and the results obtained with the two LDME sets, for the center-of-mass energy √ s = 65 GeV.
(see right panel). Indeed, in contrast to the CSM, the Color Octet mechanism generates a large upper bound for the A cos2(φ T −φ ⊥ ) asymmetry. Furthermore, the NRQCD results exhibit a strong dependence on the choice of the different LDME sets. Clearly, our process can lead to large asymmetries, with upper bounds in certain regions reaching almost 60%. Notice that we have selected specific kinematic ranges where the asymmetries could be potentially very large and where, consequently, any measurement could unambiguously allow to put strong constraints on the TMDs under study. . Two sets of LDMEs are used: SV [47] and CMSWZ [48]. The ratio of the collinear quark contribution over the total unpolarized cross section is shown in the form of a band, representing the scale uncertainty and the results obtained with the two LDME sets, for the center-of-mass energy √ s = 65 GeV.       The y dependence of the above asymmetries for J/ψ production is presented in Fig. 5, for a fixed value of K ⊥ = 2 GeV, and again both for the full NRQCD calculation (solid lines) and the CSM (dashed lines). Since √ s is fixed, y is consistently cut at y = 0.1 to fulfill all kinematic constraints. In the opposite limit, y → 1, the A cos 2φ T asymmetries vanish, as expected from the expression in Eq. (32). The maximum of A cos 2(φ T −φ ⊥ ) , on the other hand, is almost independent of y.
In Fig. 6, the upper bounds for the same asymmetries are presented for the Υ(1S) quarkonium state, whose mass is taken to be M = 9.5 GeV, and for which the single LDME set available, see Table II, is used (taking m b = 4.2 GeV). Moreover, to avoid potentially large logarithmic effects, we choose a virtuality of the same order as the Υ mass:   contribution to the unpolarized cross section could become important and no more negligible. As is evident from Eqs. (32) and (33), the asymmetries A cos 2φ T and A cos 2(φ T −φ ⊥ ) are sensitive to the ratio of linearly polarized and unpolarized gluon TMDs. As discussed in the previous section, in the case of a large nucleus (and at low x), analytical expressions for these TMDs are available within the MV model. Adopting this model also for a proton target and at x 10 −2 , allows us to calculate the q T dependence of the above asymmetries, for fixed values of the other variables. As can be seen in Fig. 7 for J/ψ production, using Q 2 = 10 GeV 2 everywhere, and K ⊥ = 2 GeV, z = 0.7, and y = 0.3 (left panel) or K ⊥ = 6 GeV, z = 0.3, and y = 0.7 (right panel), the corresponding asymmetries in the full NRQCD calculation can still be large, particularly A cos 2(φ T −φ ⊥ ) .
Finally, in Fig. 8 we illustrate the NRQCD decomposition in its different waves for the upper bound of A cos 2φ T for J/ψ. As one can see, for the CMSWZ set (left panel) there is a clear dominance of one single wave: 1 S

VI. CONCLUSIONS AND OUTLOOK
In this work, we have studied the quarkonium + jet electroproduction under the assumption of TMD factorization, within NRQCD. We have been able to identify broad kinematic regions where the quark contribution is negligible, and where the cross section can be analyzed only in terms of gluon TMDs. The specific azimuthal modulations entering the cross sections could allow for a direct access to important TMDs that are still completely unknown, like those involving linearly polarized gluons in unpolarized or transversely polarized nucleons.
With the help of positivity bounds for these TMDs we have demonstrated that, over a range of accessible regions of the phase space, the arising azimuthal asymmetries are potentially very large. Consequently, this process could be an excellent way to experimentally access the ratio of the linearly polarized gluon TMD over the unpolarized one at a future Electron-Ion Collider.
Besides these model independent estimates, we have also considered a nonperturbative model, valid in the small-x regime, namely the MV model, and presented some predictions which could be tested against experimental data.
Like any leading-order calculation, our work can be improved upon in different ways. Notably, should one aim to reliably extract gluon TMDs from precision data (at present not available) on the process under study, our computation of the cross section could be extended to take the following aspects into account. Firstly, the hadronization of the heavy-quark pair state into the quarkonium involves the emission of soft gluons, which can smear its transverse momentum over a range of the order of 2m c v, an effect which can be encoded in the so-called quarkonium TMD shape functions recently introduced in [50,51]. Likewise, the precision with which the transverse momentum of the gluon can be reconstructed is unavoidably limited by a factor of order Λ QCD due to hadronization. Moreover, an additional uncertainty can be expected due to wide-angle radiation that escapes the jet algorithm, and hence causes a mismatch between the momentum of the reconstructed jet and the initiating gluon. As was pointed out in Refs. [52][53][54], a part of this mismatch can be accounted for with the inclusion of TMD jet functions in the cross section.
In order to be able to include the above mentioned improvements to our computation in a systematic way, in addition to the Sudakov resummations and the higher order contributions in perturbation theory, it would be highly desirable to have a proof of TMD factorization for this process (e.g. in analogy with the one recently formulated for p p → η c X [50]). Such a factorization theorem could also shed light on the precise role of the CO long-distance matrix elements at small transverse momenta, and pave the way to a dedicated extraction of them in the TMD regime. Indeed, at our present precision, the largest source of theoretical uncertainty comes from the LDMEs which, although they are expected to be universal, vary significantly among different (collinear) fits.
We believe that this analysis represents an important step towards a better understanding of gluon TMDs. Moreover, our findings could be relevant for the study of their process dependence, in particular when compared to similar analyses for proton-proton collisions, which could be performed at various ongoing or planned experiments at RHIC and LHC [34][35][36].