Process dependence of the gluon Sivers function in $p^\uparrow p \to J/\psi + X$ within a TMD approach in NRQCD

We consider the transverse single-spin asymmetry (SSA) for $J/\psi$ production in $p^\uparrow p \to J/\psi + X$ within a TMD approach in non-relativistic QCD. Extending a previous study [1], we employ here the color-gauge invariant generalized parton model (CGI-GPM), in which spin and intrinsic transverse momentum effects are taken into account, together with leading-order initial- and final-state interactions (ISIs and FSIs). We find that, even when the heavy-quark pair is produced in a color-octet state, ISIs and FSIs lead to a nonvanishing SSA, allowing, in principle, to test the process dependence of the gluon Sivers function (GSF). We show that of the two independent contributions, due to the so-called $f$- and $d$-type GSFs, appearing in the CGI-GPM, the $d$-type one turns out to be dynamically suppressed. Therefore, as already found adopting the Color-Singlet Model approach for the $J/\psi$ formation [2], only the $f$-type GSF could play a role in phenomenology. A comparison with the corresponding results obtained in the generalized parton model, without the inclusion of ISIs and FSIs, is also carried out.


I. INTRODUCTION
The study of the three-dimensional structure of hadrons is of fundamental importance for our comprehension of their properties. It has certainly reached a substantial level of accuracy, thanks to many efforts, carried out in the last two decades both theoretically and experimentally. Its understanding in terms of transverse momentum dependent parton distributions (TMDs) represents the main achievement in this context [3,4]. These functions, nonperturbative in nature, have been extracted from several fits to experimental data, coming from semi-inclusive deep inelastic scattering (SIDIS) and Drell-Yan (DY) processes. As a matter of fact, most of the information collected so far is mainly restricted to the quark sector, while gluon TMDs are still very poorly known.
Among the eight leading-twist nucleon TMDs, the Sivers function [5,6] plays a seminal role. It describes the asymmetric azimuthal distribution of unpolarized partons (quark and gluons) in a fast-moving transversely polarized nucleon and is related to the orbital motion of partons. It could be responsible for azimuthal and single-spin asymmetries (SSAs) in processes where one of the initial nucleons is transversely polarized w.r.t. its direction of motion. Another important feature is its expected process dependence. This can be understood in terms of initial-and finalstate interactions (ISIs and FSIs), encoded in Wilson lines (gauge links), essential to preserve gauge invariance. One of the most striking consequences is the expectation of a sign change between the Sivers function probed in SIDIS w.r.t. the one probed in DY processes [7]. This is usually referred to as modified universality of the Sivers function.
For SIDIS and DY processes TMD factorization has been proven [8][9][10], and their analyses are well consolidated, in contrast to more inclusive processes like pp → h + X, where a TMD scheme, referred to as the generalized parton model (GPM) [11], is adopted as a phenomenological Ansatz. On the other hand, its success in describing many polarization observables and its role in looking for potential factorization breaking effects make it an important tool.
For these processes, the color gauge invariant extension of the generalized parton model (GPM), known as the CGI-GPM, has been developed in Refs. [12,13] and further extended in Ref. [2]. In this approach, ISI and FSI are taken into account assuming a single eikonal gluon exchange between the struck parton and the remnants of the transversely polarized proton. This approximation is basically the leading order contribution, in a perturbative expansion of the Wilson line, in the definition of the Sivers function.
As in the quark case, the process dependence of the gluon Sivers function (GSF) can be absorbed into the corresponding partonic hard functions entering the factorized expression of the cross sections. However, two universal,

II. SINGLE-SPIN ASYMMETRIES IN THE CGI-GPM APPROACH
The SSA for the p ↑ p → J/ψ + X process is defined as follows: where dσ ↑(↓) is the differential cross section, E h d 3 σ ↑(↓) /d 3 P h , with one of the initial protons polarized along the transverse direction ↑ (↓) with respect to the production plane. We consider the proton-proton collision along the z axis in the center of mass frame, with the polarized proton moving along +ẑ, wherein the J/ψ is produced in the x − z plane, and the ↑ transverse polarization is along +ŷ. The numerator of the asymmetry receives a sizeable contribution only from the Sivers function [2] which is defined as [34] ∆f a/p ↑ (x a , k ⊥a ) ≡f a/p ↑ (x a , k ⊥a ) −f a/p ↓ (x a , k ⊥a ) FIG. 1: Diagrams for the dominant gluon fusion contribution to the process p ↑ p → J/ψ + X in the GPM (a) and in the CGI-GPM approaches with inclusion, at leading order, of additional effects from initial-state (b) and final-state ((c) and (d)) interactions. FSIs are effective only when the J/ψ is produced in a color-octet state. Notice that there are analogous diagrams for other 2 → 2 subprocesses as well as for the 2 → 1 channels, like g + g → J/ψ. The scattering amplitudes for the underlying partonic reaction, g + g → J/ψ + g, are represented by the central blobs, while the upper and lower ones describe the soft proton → gluon transitions.
This TMD describes the azimuthal distribution of an unpolarized parton a with light-cone momentum fraction x a and intrinsic transverse momentum k ⊥a = k ⊥a (cos φ a , sin φ a , 0) in a high-energy, transversely polarized nucleon with mass M p , moving along the z direction.
In order to proceed with the calculation of the asymmetry within the CGI-GPM framework we take into account the proper insertion of the leading order contribution, in the strong coupling constant power expansion, of the gauge links, for all diagrams relevant for J/ψ production in NRQCD (see, as an example, Fig. 1 for the gluon-gluon 2 → 2 channel). We note that we did not consider the FSIs of the unobserved particle (gluon) because they are known to vanish after summing the different cut diagrams, see for example the discussion in Ref. [12].
One has to include the 2 → 1 partonic subprocesses, namely g + g → J/ψ and q +q → J/ψ, as well as the 2 → 2 ones, that is g +g → J/ψ +g, g +q(q) → J/ψ +q(q) and q +q → J/ψ +g, and exploit the contributions from the 3 S states respectively. Here we refer to the standard notation for a heavy-quark pair state 2S+1 L (c) J , where S is the spin of the pair, L and J the orbital and total angular momentum and c the color configuration, with c = 1, 8. Notice that since in the 2 → 1 channels, at leading order, the J/ψ gets its transverse momentum only from the intrinsic ones of the two initial partons, these contributions can be relevant only in the low-P T region.
As one can see, in the framework of NRQCD we have to take into account contributions to the SSA coming also from the quark Sivers function, not present in the CSM, where only the gluon-gluon fusion channel is at work. Moreover, as already pointed out previously, in the CGI-GPM formalism we have to consider two possible independent sources for the GSF.
Following Refs. [1,2], to which we refer the reader for more details, the numerator of the asymmetry for 2 → 1 channels, a + b → J/ψ, i.e. g + g → J/ψ and q +q → J/ψ, within the CGI-GPM approach is given by where q = u, d, s,ū,d,s, and M inc are the hard scattering amplitudes modified by ISIs and FSIs, as detailed in the sequel. Moreover, at order O(k ⊥ / √ s), with M T = P 2 T + M 2 , being M and P T the mass of the J/ψ and its transverse momentum, respectively, and y its rapidity.
Similarly for 2 → 2 channels, a + b → J/ψ + c (i.e. g + g → J/ψ + g, g + q(q) → J/ψ + q(q) and q +q → J/ψ + g), the numerator in Eq. (1) is given by The denominator in Eq. (1) is just twice the unpolarized cross section, discussed in detail within a TMD scheme in Ref. [1]. For completeness we give the expressions separately for the 2 → 1 and the 2 → 2 channels: All the hard scattering amplitudes squared |M Inc | 2 , where we have omitted the final state quantum numbers, are calculated perturbatively by incorporating, at leading order, the ISIs and FSIs within the CGI-GPM approach, and are listed in Appendix A, for all color octet states. The expressions for the color-singlet contributions can be directly found in Ref. [2], where SSAs within the Color-Singlet Model were discussed.
Here we only point out that, for the 2 → 1 processes, due to cancellations between the ISIs and FSIs, |M Inc(f,d) | 2 = 0 for the g + g → J/ψ subprocess independently of the 2S+1 L (c) J state, leaving active only the qq channel (second line in Eq. (3)). Moreover, for the g + g → J/ψ + g subprocess, the hard parts corresponding to f ⊥g(d) 1T , |M Inc(d) | 2 , turn out to be zero for all states (see also Ref. [35]). This means that, as already found in the study of SSAs for J/ψ production within the CSM [2], only the f -type GSF enters the dominant gg channel. As we will see, this has important consequences in the phenomenological study.

III. NUMERICAL RESULTS
We proceed now with the phenomenological analysis of SSAs for J/ψ production within the CGI-GPM approach. To this aim, following Refs. [17,18], we adopt for the unpolarized TMDs a Gaussian factorized form where f a/p (x a ) is the collinear parton distribution. The Sivers function is parameterized as with where |N a | ≤ 1 and Eq. (9) can be rewritten as where ρ a = M 2 k 2 ⊥a +M 2 and 0 < ρ a < 1. With these choices, the Sivers function satisfies the model independent positivity bound for all values of x a and k ⊥a : For the collinear unpolarized parton distributions we will adopt the CTEQL1 set [36], at the factorization scale equal to M T , adopting the DGLAP evolution equations.
These parameterizations allow us to integrate analytically the expressions entering the numerator and the denominator of A N for the 2 → 1 channels, as follows: where, while in Eq. (14) (a, b) = (q,q), (g, g), in Eq. (15) only the qq channel is active. For the 2 → 2 channels we will have to proceed by a numerical integration. At this point, we have to fix all parameters entering our expressions. In order to carry out a direct and easier comparison with the corresponding analysis performed in the GPM framework and study the impact of ISIs and FSIs, we adopt the same choices made in Ref. [1]. More precisely, for the quark unpolarized Gaussian width we use k 2 ⊥q = 0.25 GeV 2 , as extracted in Ref. [37], while for the gluon TMD we use k 2 ⊥g = 1 GeV 2 , that allows for a reasonably good description of the unpolarized cross section data in the low-P T region relevant for our study [1,2].
Moving to the LDMEs, we consider the BK11 [38] and the SYY13 [39] sets, whose values are given in Appendix A. As extensively discussed in Ref. [1] (where all details and physics motivations can be found), these sets are suitable enough for a study, within a TMD framework, of the low-P T region, where also SSA data are available: see, for instance, Figs. 1 and 2, left panels, in Ref. [1], for a comparison with unpolarized cross section data.
We will start focusing on the relevance of ISIs and FSIs, and their interplay with the production mechanism, by a detailed comparison with the corresponding estimates in the GPM approach. To properly study the role of the partonic dynamics, we compute the contributions to A N by maximizing, separately, the Sivers effect both for quarks and gluons. This can be obtained by using ρ q,g = 2/3, N q (x) = +1 and N (f,d) g (x) = +1 in Eq. (12). Notice that the chosen positive normalization allows for a better understanding of the relative signs coming from the hard dynamics. Recalling that in NRQCD the heavy-quark pair can be produced in a CS or a CO state, and the latter with different angular momentum quantum numbers, we will also discuss separately each contribution. In the following, we will adopt the kinematics of the PHENIX experiment, at √ s = 200 GeV, for which SSA data are available. As we will show below, and as it happens in the GPM, also in the CGI-GPM the quark-initiated contributions are negligible. Similarly, the contribution from the d-type GSF is extremely small, due to the absence of the gg channel in the numerator of the SSA (see comment at the end of the previous section). This means that within the CGI-GPM and NRQCD frameworks, one can still concentrate on the f -type gluon Sivers function.
In Fig. 2 we show our maximized estimates for A N at √ s = 200 GeV, coming from the dominant gluon contribution (red solid lines) in the CGI-GPM, f -type, (left panel) and the GPM (right panel) approaches, at x F = 0.1 as a function of P T , together with a full wave decomposition (see the legend). Note the one order of magnitude difference in the vertical scale between the two panels. We also observe that some GPM contributions are larger than one since the denominator of the SSA includes all terms (entering with relative signs), while in the numerator we consider, term by term, only a specific wave state. The overall result (red solid lines) is, as it has to be, smaller than one.
We can further observe that, while the CGI-GPM estimates show clear oscillations, with a change of sign around P T 1 GeV for CO states, the corresponding ones within the GPM have all a definite sign. This oscillating behaviour is specifically due to the 2 → 2 channels and plays a role also in the GPM.
On the other hand, in the GPM the 2 → 1 channels, at least in the small-P T region where they are relevant, compensate for this effect, leading to an overall definite sign. In contrast, as already pointed out in the previous section, in the CGI-GPM the gluon contributions from the 2 → 1 channels are identically zero. Notice that this oscillation in sign for the 2 → 2 channels has nothing to do with the role of ISIs and FSIs (as it is clear from the fact that also the GPM estimates present this feature) and comes directly from the parton dynamics, as weighted by the Sivers azimuthal phase. This can be easily understood recalling that within the GPM the hard parts in the numerator are identical to those in the denominator, which does not manifest any oscillation in sign. We also observe ( Fig. 2) that this oscillating behaviour does not affect CS states, whose amplitudes squared present a much simpler structure in terms of their Mandelstam invariant dependence.
Another interesting feature in Fig. 2 is that among the CO contributions, two of them, the 1 S 0 and 3 P [J] (where the symbol [J] stands for a sum over J = 0, 1, 2) wave terms, are very large, almost comparable in size but opposite in sign (this is due to the sign of the corresponding LDMEs), while one, the 3 S 1 wave, is extremely small. This happens in both approaches. Moreover, the CS contribution, which as already said has a definite sign, shifts the zero in the CGI-GPM estimates to larger P T values, as compared with the CO terms.
Coming back to the size of each contribution, the smaller values in the CGI-GPM approach w.r.t. the corresponding ones in the GPM are directly due to the cancellations between different hard partonic parts, entering with proper color factors and having in some cases opposite signs.
In Fig. 3 we show the corresponding maximized estimates at √ s = 200 GeV, within the CGI-GPM, at x F = −0.1 as a function of P T (left panel) and at fixed P T = 1.65 GeV as a function of x F (right panel), as in the PHENIX analysis (see below). Concerning the backward region, the main aspect is the suppression of all contributions as compared to those at x F = 0.1 (Fig. 2, left panel), that leads to much smaller results. In fact, besides the effects already discussed, the dependence on the Sivers azimuthal phase, cos φ a (through the Mandelstam invariants, see Appendix A), is less relevant in the hard parts, and therefore the integration over it (see Eq. (5)) is more effective in reducing the effect. Moreover, at such fixed P T value (right panel), as it is evident from Fig. 2, the cancellation among the various contributions in the CGI-GPM is much more effective for all x F values, leading to very small maximized SSAs. The use of the other LDME set (SYY13) gives very similar results.
Having analysed in Figs. 2 and 3 the maximized contribution of the f -type GSF, to have a more complete view, in Fig. 4 we show for the BK11 (left panel) and the SYY13 (right panel) LDME sets a collection of results for maximized A N at √ s = 200 GeV and x F = 0.1 as a function of P T , adopting different models and approaches. PHENIX data [40] are also shown. More precisely, we present the maximized estimates obtained within NRQCD, employing both the CGI-GPM (this work) and the GPM [1], and those within the CSM in both approaches [2]. For the CGI-GPM all contributions are shown and, as already pointed out above, those from the quark (with the exception of the very small P T region for the SYY13 set) and the d-type gluon Sivers functions are negligible. We also recall that the SYY13 LDME set includes only color-octet states. This implies, still within the CGI-GPM, larger values of A N at large P T , due to the missing of the relative cancellation between the CO and the CS contributions, at work for the BK11 set. Notice that the estimates obtained within the CSM do not depend on the LDME set and, in this respect, they could appear identical also in the right panel. From these plots we see that, even if to a much lesser extent as compared to the GPM, also within CGI-GPM one can potentially put some constraints on the size of the f -type GSF already with the few data points available.
As expected from previous considerations, the corresponding results at x F = −0.1 show a quite different situation, see Fig. 5. In this configuration the maximized SSAs within the CGI-GPM approach are strongly suppressed. For the SYY13 LDME set, the situation looks only slightly different since the absence of the CS contribution prevents a   further relative cancellation. On the other hand, even maximizing the GSF, the estimates also in this case are already very close to the data. In this respect, within the CGI-GPM approach data at negative x F appear not very useful to constrain the GSF. This is in contrast to what happens in the GPM framework where, both in CSM and in NRQCD, only a strongly suppressed GSF w.r.t. its positivity bound could give estimates compatible with PHENIX data.
Complementary information can be also obtained by looking at the same quantities at fixed P T (chosen here to be 1.65 GeV as in PHENIX data) as a function of x F . This is shown in Fig. 6, where, adopting the same choices as in Figs. 4 and 5, we see that for such P T values one cannot put any constraint on the f -type GSF. A better configuration, in this respect, would be exploring larger P T values (around 2-3 GeV) and, with some care, very low-P T values (below 1 GeV), in the positive x F region (see Fig. 4). In such cases the maximized A N would be sizeable enough and any data could help in constraining the GSF within a CGI-GPM approach.
Before concluding this comparison, we have to mention that the use of the available extraction of the f -type (and a fortiori the d-type) GSF, from midrapidity pion and D-meson SSA data [18], would give results 0.05 (0.15) times smaller than the corresponding maximized estimates, that is almost compatible with zero, and in reasonable agreement with J/ψ data. For its relevance, we also present estimates for the corresponding A N in J/ψ production for the kinematics reachable at LHC in the fixed target mode with a transversely polarized target (see the AFTER [41,42] and LHCSpin [43,44] proposals at CERN). In such a configuration one could probe even larger light-cone momentum fractions in the polarized proton, accessing the gluon TMDs in a very interesting and complementary region.
In Fig. 7 we present our maximized estimates for A N for pp ↑ → J/ψ + X at √ s = 115 GeV, at fixed P T = 3 GeV, as a function of x F (left panel) and at fixed rapidity y = −2, as a function of P T (right panel), adopting the BK11 set. Notice that in such a configuration the backward rapidity region refers to the forward region for the polarized proton target. As one can see, at P T = 3 GeV (left panel) the maximized contribution from the f -type GSF at backward rapidity is around 5% and, in principle, could be accessed/constrained experimentally. The same is true at very small (< 1 GeV) or large (≥ 3 GeV) P T values at y = −2 (right panel). We notice that the corresponding estimates, from the f -type GSF, at P T around 2 GeV would be almost negligible at all rapidities. As already discussed in our previous studies, adopting the GPM or the CGI-GPM together with the CSM, the maximized A N would be much larger and, potentially, easier to constrain.

IV. CONCLUSIONS
In this paper we have extended, and somehow completed, a detailed analysis of SSAs for J/ψ production in pp collisions within a phenomenological TMD scheme. This study started in a previous paper, where, employing the Color-Singlet Model for quarkonium formation, we compared the Generalized Parton Model and the Color-Gauge-Invariant GPM. It has been then continued quite recently in a second work, adopting the NRQCD framework within the GPM. Here we have eventually considered its extension within the CGI-GPM. The main interest of this analysis is to see whether and to what extent one can extract information on the poorly known gluon Sivers function, focusing only on this specific process.
We have considered all relevant subprocesses in NRQCD, both for the 2 → 1 and the 2 → 2 channels, including effects of initial and final state interactions, in the one-gluon-exchange approximation. This leads to the introduction of new color factors, diagram by diagram, and the computation of modified hard scattering amplitudes. In such a way one can move the process dependence, coming from ISIs and FSIs, into the hard parts, factorizing the corresponding TMDs. One, well-known, outcome of this approach is the appearance of two independent gluon Sivers functions, referred to as the d-type and the f -type distributions.
We have then calculated the maximized contributions to A N , separately for the gluon and the quark Sivers effects, adopting the kinematics of the PHENIX experiment, for which data are available. The main findings are that the quark as well as the d-type gluon Sivers functions, even if maximized, give almost negligible contributions to the SSA, leaving at work, as in the CSM, only the f -type GSF. On the other hand, within NRQCD this contribution is also generally quite small and could be relatively sizeable only at forward rapidities and P T around 2-3 GeV, at least for the two LDME sets considered.
Therefore, while within the GPM, the GSF could be easily constrained by PHENIX SSA data for J/ψ production alone, the situation in the CGI-GPM is quite different. Indeed, if one adopts the CSM, the f -type GSF (the only one active) gives still a potentially sizeable contribution; on the contrary, in full NRQCD it could be hardly constrained, and definitely not in the backward region.
We have also presented some maximized estimates of A N , for the kinematics reachable at LHC in a fixed target mode, showing similar features as those discussed for PHENIX setup.
More data, with higher statistics, could certainly help in shedding light on the role of the gluon Sivers function, as well as on its process dependence.
Appendix A: Color factors and amplitudes squared in p ↑ p → J/ψ + X for color-octet states within the CGI-GPM approach Here we collect all color factors as well as the amplitudes squared for the relevant subprocesses in p ↑ p → J/ψ + X within the CGI-GPM approach (for color-octet states). The color factors and the corresponding amplitudes squared for the color-singlet states can be found in Ref. [2].
The modified amplitudes squared in the CGI-GPM or, more precisely, each contribution to the MM * product between any two of the Feynman diagrams for the specific subprocess, can be written as where M U are the scattering amplitudes for the unpolarized partonic processes.
The Feynman diagrams contributing to the g + q → J/ψ + q channel are shown in Fig. 9, while the corresponding color factors are collected in Table III (upper signs).
Here, as well as in the following Table IV, the label A refers to the 1 st Feynman diagram, B to the 2 nd , while C represents the grouping of the 3 rd , 4 th and 5 th diagrams in Fig. 9. For the 1 S states, the 1 st , 2 nd and 3 rd Feynman diagrams in Fig. 9 do not contribute.
Collecting together all contributions of Fig. 9 with the appropriate color factors from Table III, we find the following expressions for the amplitudes squared: FIG. 9: Feynman diagrams for the g + q → J/ψ + q process. states for the g + q → J/ψ + q (upper signs) and g +q → J/ψ +q (lower signs) channels.
The resulting amplitudes squared are the following: The color factors for the q + g → J/ψ + q channel are listed in Table IV (upper signs), and the corresponding amplitudes squared are: states for the q + g → J/ψ + q (upper signs) andq + g → J/ψ +q (lower signs) channels.

State
Diagram The color factors for theq + g → J/ψ +q channel are listed in Table IV (lower signs) and the corresponding amplitudes squared are: The Feynman diagrams contributing to the q +q → J/ψ + g channel are shown in Fig. 10, while the corresponding color factors are collected in Table V (column qq and upper signs). Here the label A represents the grouping of the 1 st , 2 nd and 3 rd Feynman diagrams, B refers to the 4 th diagram, while C refers to the 5 th diagram in Fig. 10. For 1 S  Fig. 10 do not contribute.  states for the q +q → J/ψ + g (column qq, upper signs) andq + q → J/ψ + g (columnqq, lower signs) channels.

State
qq Diagramqq Diagram CU CI CF C Inc 3 S The resulting amplitudes squared are: (A32) 7.q + q → J/ψ + g channel The color factors for theq + q → J/ψ + g channel are listed in Table V (columnqq and lower signs) and the corresponding amplitudes squared are: The Feynman diagrams for the g + g → J/ψ process are shown in Fig. 11, while in Table VI we give the color factors for the sum of them. The amplitude for the 3 S    state for the q +q → J/ψ (upper signs) and theq + q → J/ψ (lower signs) channel.

CU CI CF C Inc
. q +q → J/ψ andq + q → J/ψ channels The Feynman diagram for the q(q) +q(q) → J/ψ channel is shown in Fig. 12, while in Table VII we give the corresponding color factors: qq (upper signs),qq (lower signs). Here only the 3 S 1 ]| 2 qq . (A36)

Long Distance Matrix Elements
In Table VIII we collect the values of the LDMEs adopted in the present study: