Unraveling the Gluon Sivers Function in Hadronic Collisions at RHIC

We study the transverse single-spin asymmetries for $p^\uparrow p\to \pi\, X$ and $p^\uparrow p\to \gamma\, X$ within the so-called color gauge invariant generalized parton model (CGI-GPM) which, in addition to spin and transverse momentum effects, includes initial and final state interactions with the polarized proton remnants. We compute all relevant contributions, focusing in particular on the process dependence of the gluon Sivers function, which, for these processes, can always be expressed as a linear combination of two independent, universal terms. This study extends and completes a previous one, where only quark initiated partonic processes were considered. We then perform a combined phenomenological analysis of RHIC data on transverse single-spin asymmetries in $p^\uparrow p\to \pi\, X$ and $p^\uparrow p\to D\, X$, putting the first preliminary constraints on these two gluon Sivers functions. We show how their size can be estimated by means of these data, and use our results to provide predictions for the process $p^\uparrow p\to J/\psi\,X$, comparing them with data, and $p^\uparrow p\to \gamma\, X$, for which experimental information will soon become available. Corresponding estimates within the simpler GPM approach, without initial and final state interactions and with a single universal gluon Sivers function, are also given, showing that a clear discrimination between these two models is, for the moment, not possible.


I. INTRODUCTION
Among the various transverse momentum dependent parton distribution and fragmentation functions (TMDs for short), the Sivers function [1,2] is of great interest, both experimentally and theoretically. It is related to the asymmetry in the azimuthal distribution of unpolarized quarks and gluons inside a high-energy proton that is transversely polarized with respect to its momentum. As such, it can in turn give rise to azimuthal asymmetries of the produced particles in high-energy scattering processes initiated by transversely polarized protons. Moreover, the Sivers function is known to be very sensitive to the color exchanges among initial and final states, and to the color flow in the scattering processes. These peculiar properties have a clear signature [3,4], providing a strong test of the TMD formalism.
A first evidence of a nonzero Sivers distribution for quarks has come from data on single spin asymmetries for semi-inclusive deep inelastic processes (SIDIS), measured by the HERMES Collaboration at DESY [5], and confirmed later by the COMPASS Collaboration at CERN [6]. Nowadays, thanks to a continuous and dedicated experimental investigation and to new phenomenological extractions, it can be considered established.
The knowledge of the quark Sivers function, quite important by itself, provides an indirect constraint on the much less known gluon Sivers function by means of the Burkardt sum rule [7], which states that the transverse momenta of all unpolarized partons inside a transversely polarized proton add up to zero. Available parameterizations for the quark Sivers function [8,9] almost fulfill, within uncertainties, the Burkardt sum rule, pointing towards a small gluon contribution. This is consistent with theoretical arguments valid in the large-N c limit of QCD [10,11], according to which the gluon Sivers function should be suppressed by a factor 1/N c with respect to the valence quark Sivers distributions at values of the light-cone momentum fraction x of the order of 1/N c .
Turning now to the discussion of direct probes of the gluon Sivers effect, we note that a first extraction of the gluon Sivers function from very precise data on single spin asymmetries in p ↑ p → π 0 X at central rapidities [12] has been attempted in the framework of the generalized parton model (GPM) [13]. In this approach, the TMD formalism is applied even to single-scale processes and transverse momentum dependent distribution and fragmentation functions are conditionally taken to be universal. Although lacking a formal proof, the GPM is phenomenologically very successful in describing many processes for which data are available, see Refs. [14][15][16][17][18][19][20].
In the meantime, a Color Gauge Invariant formulation of the GPM, named CGI-GPM [21][22][23] has been proposed, in which the effects of initial (ISI) and final (FSI) state interactions on the quark Sivers function are taken into account, within a one-gluon exchange approximation. As a result, the Sivers function for quarks becomes nonuniversal, and its process dependence can be absorbed into the partonic cross sections. Hence, in the calculation of physical observables, for example in proton-proton collisions, one can still use the quark Sivers functions obtained from SIDIS data, but they need to be convoluted with the modified partonic cross sections calculated in Ref. [21]. In particular, the CGI-GPM can reproduce the expected opposite relative sign of the quark Sivers functions in SIDIS and in the Drell-Yan processes [3,4].
In Ref. [24], the CGI-GPM has been for the first time extended to the gluon Sivers function in the study of inclusive J/ψ and D meson production in proton-proton collisions at RHIC. These processes, as compared to pion production, have the advantage of probing gluon TMDs directly, since quark induced subprocesses can be safely neglected in the kinematical regions considered. Similarly to the quark case, the process dependence of the gluon Sivers function can still be absorbed into the hard partonic cross sections. However, one needs to introduce two different classes of modified partonic cross sections, corresponding to the two different ways in which a color-singlet state can be formed out of three gluons, i.e. either through an antisymmetric or a symmetric color combination. Each one of them has to be convoluted with a different gluon Sivers distribution. These two universal and independent distributions are named, respectively, the f -type and d-type gluon Sivers functions [25], or A 1 and A 2 in the notation of Ref. [26]: the former is even under charge conjugation, while the latter is odd. It turns out that only the f -type distribution contributes to J/ψ production, at least in the analyzed kinematical region where the color-singlet mechanism is dominant, while for D-meson production the d-type is the most relevant one [24]. Corresponding studies, within the GPM framework only, have been presented in Ref. [27] and later on in Refs. [28,29].
In the present paper we extend the formalism of the CGI-GPM to the processes p ↑ p → π X and p ↑ p → γ X. We calculate all modified partonic cross sections induced by gluons, needed for a re-analysis of the RHIC pion data of Ref. [12]. These results are therefore complementary to the quark-induced ones published in Ref. [21]. Moreover, we perform a detailed phenomenological analysis and show how it is possible to disentangle and give an estimate of the size of the two gluon Sivers functions. To this end we study, in the same framework (see Ref. [24]), also the latest available data on inclusive D-meson production [30]. We then compare our new predictions for single spin asymmetries in p ↑ p → J/ψ X with the most recent RHIC data [31] and give the corresponding theoretical estimates for the kinematics reachable at LHC with a fixed polarized target. Finally, we give predictions for the process p ↑ p → γ X currently under investigation at RHIC, for which data are expected in the near future.
The paper is organized as follows: in Section II we present the leading order partonic cross sections, within the framework of the CGI-GPM, for the gluon induced subprocesses that contribute to the Sivers asymmetry in p ↑ p → h X (Sect. II A) and in p ↑ p → γ X (Sect. II B). In Section III we perform a phenomenological analysis of available data on single spin asymmetries in p ↑ p → π X and p ↑ p → D X putting some reliable constraints on the gluon Sivers function, then in Section III B we present our predictions for the same observable in p ↑ p → J/ψ X (for which a comparison with data is possible) and p ↑ p → γ X. Conclusions and final remarks are collected in Section IV. The color factors needed for the calculation of the hard functions H Inc ab→cd within the CGI-GPM are listed in the two Appendices.

II. THEORETICAL FRAMEWORK
The single-spin asymmetries (SSAs) for the processes p ↑ p → h X and p ↑ p → γ X are defined as follows where dσ ↑(↓) denotes the single-polarized cross section, in which one of the protons in the initial state is polarized along the transverse direction ↑ (↓) with respect to the production plane. As extensively studied in Ref. [16], within a TMD approach, the numerator of the asymmetry is mainly driven by only two contributions: the Sivers [1,2] and the Collins [32] effects. Furthermore, in suitable kinematical regions, as we are going to discuss below, only the Sivers effect can be sizeable. Hence, the numerator of the asymmetry is sensitive to the quantity [33] ∆f withf a/p ↑ (x a , k ⊥a ) being the number density of partons a with light-cone momentum fraction x a and transverse momentum k ⊥a = k ⊥a (cos φ a , sin φ a ) inside the transversely polarized proton with mass M p , which is taken to move along theẑ-axis. The Sivers distribution of parton a is represented either by ∆ N f a/p ↑ (x a , k ⊥a ) or f ⊥a 1T (x a , k ⊥a ) and fulfills the following positivity bound We note that, since a can be either a quark (antiquark) or a gluon, the Sivers contribution to the asymmetry can be expressed as a sum of two terms, namely where quark (gluon) refers to the parton inside the polarized proton in the numerator of A N . The quark and gluon contributions to A N cannot be directly disentangled either in p ↑ p → π X or in p ↑ p → γ X. For this reason, in our numerical studies, focused on the extraction of the gluon Sivers function, we will use all the available information on the quark Sivers functions coming from the analysis of azimuthal asymmetries in SIDIS processes.
In the next two subsections, we provide the explicit expressions of the numerators of the asymmetries for p ↑ p → π X and p ↑ p → γ X, respectively, in the CGI-GPM approach. The corresponding formulae for p ↑ p → J/ψ X and p ↑ p → D X are given in Ref. [24], where it was found that, for such processes, the gluon contribution to the asymmetry is dominant.
Within the framework of the CGI-GPM, the numerator of the asymmetry is given by where J(k ⊥π ) is a kinematical factor [15] andŝ,t,û are the usual Mandelstam variables for the partonic subprocess ab → cd. Furthermore, f b/p (x b , k ⊥b ) is the TMD distribution for an unpolarized parton b inside the unpolarized proton, while D π/c (z, k ⊥π ) is the the fragmentation function of an unpolarized parton c into a pion. Finally, H Inc ab→cd are the perturbatively calculable hard scattering functions. In particular, the ones for which a is a quark or an antiquark, are well-known and can be found in Ref. [21], while the remaining ones have been evaluated here for the first time along the lines of Ref. [24]. As already pointed out, in the CGI-GPM approach there are two independent gluon Sivers contributions: the f -and d-type. The leading order (LO) explicit expressions for the hard functions corresponding to the gluon Sivers distribution f H Inc (f ) where N c is the number of colors. For the other gluon Sivers function f ⊥g (d) 1T , one has H Inc (d) More details on their calculation are given in Appendix A. For comparison, we show the corresponding, well-known unpolarized hard functions, defined in such a way that which appear in the denominators of the asymmetries.
The numerator of the SSA for the process p ↑ p → γ X reads As for p ↑ p → π X, the partonic hard functions in which the parton a inside the polarized proton is a quark or an antiquark are given in Ref. [21]. For the gluon induced subprocesses, we find for the f -and d-type gluon Sivers functions, respectively. The unpolarized hard function is given by and is normalized such that the corresponding partonic cross section has the following form: We refer to Appendix B for further details of the calculation.

III. PHENOMENOLOGY
We are now able to devise a possible strategy to put the first reliable constraints on the two independent gluon Sivers functions within the CGI-GPM approach. To this aim, in Section III A we will present a detailed analysis of SSA data in p ↑ p → π X and p ↑ p → D X. We will compare our findings with the available data, as well as with the corresponding results in the GPM scheme, as obtained in Ref. [13]. Finally, in Section III B we will show new predictions for SSAs in p ↑ p → J/ψ X and p ↑ p → γ X.
A. Constraints on the gluon Sivers functions from available data As discussed in the previous Section, in the CGI-GPM framework there are two universal and independent gluon Sivers functions (GSFs), the f -and d-type, and the phenomenological analysis appears more difficult with respect to the one in the GPM scheme. The reason is that, in principle, different combinations of these two contributions could lead to similar results and describe equally well the same set of data. Therefore, in order to carry out this analysis we will have to use at least two independent sets of data. In particular, we will use the extremely precise and accurate data on SSAs in pp collisions for inclusive pion production at mid-rapidity [12] and those for D-meson production [30] by the PHENIX Collaboration. They also collected SSA data for J/ψ production [31], which we will compare against our estimates. From the phenomenological point of view, it is worth noticing that for the latter process, in the CGI-GPM approach, only the f -type contribution appears. Therefore, as it will become more clear in the following, it is important to consider additional processes, where also the d-type GSF plays a role.
All these processes have a common feature: the gluon initiated subprocesses dominate over the quark ones. As was already pointed out in Refs. [13,34], the SSA for inclusive pion production in pp collisions at mid-rapidity is directly sensitive to the gluon Sivers distribution. In fact, the contribution involving the quark Sivers functions, as extracted from SIDIS azimuthal asymmetry data, is totally negligible -this is true also in the CGI-GPM approach, as we will show in the following -and all other effects, like the one driven by the Collins function, are washed out by integrations over the azimuthal phases. Concerning the SSAs in D-meson production, as discussed in Ref. [24], one has a clear and direct access to the GSF, due to the dominance of the gg → cc channel.
Within our strategy, the first issue we address is to which extent the f -and d-type contributions are effectively relevant in the process under consideration. More precisely, we start with the observation that the numerators of the SSAs, Eqs. (5) and (18), contain three fundamental quantities: the azimuthal factor of the gluon Sivers function, cos φ a (with φ a to be integrated over), the perturbatively calculable hard partonic parts, H ab→cd , and the unknown GSF, f ⊥g 1T . In order to explore the role played by the first two factors, we calculate the SSAs by maximizing the corresponding GSFs. To do this we adopt the well-known Gaussian-like and factorized parametrization for the GSF, as follows: where f g/p (x) is the standard unpolarized collinear gluon distribution, with |N g | ≤ 1, and Alternatively, if we define the parameter such that 0 < ρ < 1, then Eq. (23) becomes With these choices, assuming that the unpolarized TMD gluon distribution is given by the Sivers function automatically fulfills its proper positivity bound for any (x, k ⊥ ) values (see Eq. (3)). Analogously, for the unpolarized TMD fragmentation function (for a parton c) we use [35] D π/c (z, k ⊥π ) = D π/c (z) x F = 0 √s = 200 GeV GPM Data are from Ref. [12].
In this analysis we adopt the CTEQ6-LO parametrization [36] for the unpolarized gluon distribution, f g/p (x), with the factorization scale equal to the pion transverse momentum, p T , and the leading-order DSS set for the collinear fragmentation functions [37]. Notice that all TMDs defined above evolve with the hard scale through the scale dependence of the collinear distributions entering in their parameterizations, that is following a DGLAP evolution.
The first k ⊥ -moment of the Sivers function is also of relevance: Adopting the parameterization of Eqs. (23)- (25), In Ref. [13] a single value k 2 ⊥ = 0.25 GeV 2 [35] was adopted, the same for the unpolarized quark and gluon TMDs, while the parameters N g , α, β, ρ were fitted to the data, within the GPM scheme. Here, following Ref. [24], for the unpolarized gluon TMD we use a different value, k 2 ⊥ = 1 GeV 2 . This, indeed, gives a better account of the unpolarized cross sections for J/ψ production at not so large p T values, still allowing a good description, for instance, of the inclusive pion production. For this reason, we have reanalysed the same set of data within the GPM approach, getting results very similar to those reported in Ref. [13], although with slightly different parameters: Notice that an equally good description of pion SSA data can be obtained even with different sets of the above parameters, that are strongly correlated among each other. While this could imply very different k ⊥ dependences of the GSF, its first k ⊥ -moment remains almost unchanged in the range of x probed by data (10 −3 ≤ x ≤ 0.4).
Moving to the CGI-GPM approach, in order to maximize the effects of the GSFs, we saturate the positivity bound for their x-dependent parts (i.e. we take N g (x) = ±1) and adopt the value ρ = 2/3 [38] in Eq. (27).
For the x-dependent part of the GSF one can also use the following notation which, for N g (x) = ±1, implies ∆ N f g/p ↑ (x) = ±2f g/p (x). In Fig. 1 (left panel) we present the maximized (N g (x) = +1) gluon Sivers contributions to A N for the process p ↑ p → π 0 X at √ s = 200 GeV and mid-rapidity as a function of p T , together with PHENIX data [12], for the f -type (red solid line) and d-type (blue dot-dashed line) pieces. For completeness we also show the maximized gluon Sivers term in the GPM (green dashed line). As mentioned above, the quark Sivers contribution, also within the CGI-GPM scheme and adopting the parametrization as extracted from SIDIS data [9], is totally negligible (red dotted line). From this plot we realize that while the d-type contribution, for this process and in this kinematical region, is dynamically suppressed, the f -type one can be potentially large. The reason is that for the d-type term the hard partonic cross sections for the processes initiated by gq and gq pairs enter with a relative sign (see Eqs. (10) and (11)) and at mid-rapidity the quark and anti-quark unpolarized TMD parton distributions are equally important. On top of that, there is no gg → gg contribution (see Eq. (13)), the dominant channel at moderate values of p T . This is in contrast with the f -type term, which indeed could be potentially very large. We also notice that the corresponding effect in the GPM approach is even larger: the reason is that its partonic contributions are exactly those entering the unpolarized cross section, all positive and unsuppressed.
These considerations lead us to the second step of our strategy: the attempt to describe reasonably well the A N data for π 0 production at mid-rapidity within the CGI-GPM approach, by adopting at the same time the most conservative (that is less stringent) bounds on the f -and d-type GSFs. Notice that in the region where they are more precise (p T 5 GeV), the data are tiny, of the order of per mille, and positive. It is then clear that the most conservative scenario that could give SSAs comparable to the data implies a cancellation between the two contributions, with a strongly suppressed and positive f -type GSF and a saturated, negative d-type one (supposed totally unknown). The corresponding results, for N g . Notice that a smaller d-type GSF (in size, that is either positive or negative) would imply an even smaller f -type GSF. This issue will be addressed in the following.
Let us now consider A N for D 0 production at √ s = 200 GeV in the kinematical region relevant to carry out the corresponding analysis for its muon decays, for which data are available [30]. Actually, to be more general, we consider an even larger region both in x F = 2p L / √ s (where p L is the D meson longitudinal momentum) and p T . In Fig. 2 we show the results for A N as a function of x F and for different p T values, obtained by separately maximizing the d-(left panel) and f -type (right panel) contributions, as explained above. One can see that in the forward region, while the d-type term could be sizeable, the f -type one is relatively small. This is in contrast to what was discussed above for the case of π 0 production. The reason is that, since for D 0 production at leading order we consistently consider only the dominant fragmentation of the charm quark into the heavy meson, the cancellations between the gq and gq initiated processes, affecting the previous case, are not present anymore. Moreover, the hard partonic parts favor the d-type w.r.t. the f -type term: as one can see from Eq. (41) of Ref. [24], besides some common factors, the hard part for the f -type GSF contains a factort 2 /ŝ 2 , whilst that for the d-type GSF contains a term (t 2 − 2û 2 )/ŝ 2 . Since |t| becomes smaller and smaller as x F increases, the first piece is relatively suppressed w.r.t. the second one. On the other hand, in the backward region, where the two hard parts are similar, both contributions are relatively suppressed by the integration over the Sivers azimuthal phase, which for x F < 0 is less effective in the hard parts.
If we now use the information extracted from the analysis of π 0 SSA data, the f -type contribution in Fig. 2 should be accordingly reduced by a factor of about 0.1 (coming from the corresponding GSF), thus becoming practically negligible. This implies that for D 0 production only the d-type GSF could be considered active and therefore constrained by a comparison with the available data. Similar considerations apply also toD 0 production. In this case, as discussed in Ref. [24], within the CGI-GPM approach the f -type contribution to A N is the same as the one for D 0 production, while the d-type gets an opposite sign.
At this point, one has to convert the estimates for D-meson production to the corresponding SSAs for its muon decay products 1 , for which data are available [30]. Notice that in our LO approach the SSAs for D 0 and D + production (leading to the µ + results) are equal, as are those forD 0 and D − production (µ − results).
Since the muon SSA data are still very few and with large error bars, we refrain from performing a fit, and will consider a simple x-independent N . In the following we discuss different possible scenarios for the d-and f -type GSFs, taking into account the complementary information coming from π 0 SSAs. As we will see in a moment, even from this very conservative approach we can extract some important information.
As one can see from 1, was adopted in order to reasonably reproduce the π 0 SSA data (see Fig. 1, right panel). On the other hand, to get a fair account of the muon SSA data, one has to take indicatively |N production. Taking into account this new piece of information on the d-type GSF, we can reconsider the pion SSA data more accurately. We find that by varying N (d) g in the range −0.15 ÷ +0.15, while keeping ρ = 2/3, a very good description of both the µ ± and π 0 data can be obtained by taking N (f ) g in the corresponding range +0.05 ÷ −0.01, that is: In other words, a stronger suppression of the f -type GSF is required by the combined analysis of muon and pion SSA data. On the contrary, in the GPM approach the parametrization of the GSF extracted from the π 0 SSA data, see Eq. (32), leads to SSAs for µ ± leptons in very good agreement with available data (Fig. 3, green dashed lines). For completeness, in Fig. 4 we also show the corresponding SSA estimates as a function of p T in the positive and negative x F regions. It is worth recalling that a similar analysis of SSAs for D-meson production, within the twist-three approach, was carried out in Ref. [39]. The corresponding predictions for muon production were compared against the data in Ref. [30], showing a similar fairly good agreement.
A few comments on the above procedure are in order. The use of a fixed ρ value implies a fixed k ⊥ dependence of the GSF, therefore no such information has been extracted within the CGI-GPM approach. On the other hand, the adopted value leaves the size of the GSF practically unconstrained. Then, by tuning the parameter N g against the data we can control and estimate its size. We have also to remind that there are strong correlations between these parameters, but the amount and the precision of available data, as already stated above, prevent us from performing a true fit.
For all these reasons, in the following we will show only the first k ⊥ -moment of the GSF, which better represents its size in an almost unbiased form (at least in the x region probed by the data, 10 −3 ≤ x ≤ 0.4), without speculating on its detailed k ⊥ or x dependences. Further studies in this respect will be necessary.
In Fig. 5 we show the results for the absolute value of the first k ⊥ -moment of the GSFs as extracted from our analysis for the GPM (green dashed line) and the CGI-GPM approaches, d-type (blue dot-dashed line) and f -type (N (f ) g = 0.05, red solid line), together with the positivity bound (black dotted line). The most stringent bound is the one for the GPM approach, since in this case there are no relative cancellations between the hard partonic parts, being them all positive. In contrast, the d-type GSF within the CGI-GPM scheme is the less bounded (see comments above).
B. Predictions for SSAs in p ↑ p → J/ψ X and p ↑ p → γ X As discussed in Ref. [24], A N for J/ψ production is directly sensitive to the gluon Sivers function. Moreover, within the CGI-GPM approach and the Color Singlet model, only the f -type distribution contributes to the Sivers asymmetry. In Figs. 6 and 7 we show a comparison of our estimates, evaluated adopting M 2 T = M 2 J/ψ + p 2 T as factorization scale, with PHENIX data [31] for A N in p ↑ p → J/ψ X. In particular, in Fig. 6, left panel, we show the maximized (N With the exception of the experimental point in the most backward rapidity region, data are compatible with zero and our estimates describe them fairly well. Notice that, in principle, by using a suitable x-dependent factor, N (f ) g (x) (namely something like N g (1 − x) β , with N g −1 and a large β), also the data points at x F < 0 could be accounted  34)). Data are from Ref. [31]. Notice the different scales in the two panels.
for. On the other hand, this would prevent a description of pion SSAs at small p T , which require a strong suppression of the f -type GSF, in particular in the small-x region (see Fig. 1, left panel). If J/ψ measurements would be confirmed even in future higher statistics samples, this would definitely represent a tension with the pion SSAs, at least within a TMD approach. In this respect, more data, on a wider kinematical range and with better statistics, would be very helpful.
It is worth considering the corresponding analysis for A N in J/ψ production for the kinematics reachable at LHC in the fixed target mode with a transversely polarized target (see the AFTER [40,41] and LHCb [42,43] 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. 8 we present our estimates for A N for pp ↑ → J/ψ X at √ s = 115 GeV, at fixed p T = 2 GeV, as a function of x F (left panel) and at fixed rapidity y = −2, as a function of p T (right panel). Notice that in such a configuration the backward rapidity region refers to the forward region for the polarized proton target. In particular, we show our predictions within the GPM (thick green dashed lines) and the CGI-GPM (red bands) approaches, together with the corresponding upper/lower positivity bounds (thin lines). From these results we see that any further experimental information would be extremely useful. Another interesting observable, where the gluon Sivers function could be directly accessed, is the SSA in p ↑ p → γ X, for which we have given the complete expressions in the CGI-GPM scheme in the previous Section. We present here some estimates, both in the GPM and CGI-GPM approaches, saturating their contributions as well as adopting the results of the phenomenological analysis presented above. As for the case of SSAs in π 0 production, the most interesting regions are those at mid-and slightly backward-rapidity and not so large values of |x F |. The reason is that, at large negative values of x F , the integration over the Sivers azimuthal phase washes out the effect. This would partially spoil the analysis proposed in Ref. [44], where the authors discussed this process as a clear tool to access the GSF, also in this kinematical region.
In Fig. 9 (upper panels) we show the maximized contributions to A N at x F = 0 (left) and x F = −0.1 (right). As one can see, the d-type term at x F = 0 is dynamically suppressed, as for the π 0 production case: the reason is indeed the same, that is the partial cancellation between the hard gq → γq and gq → γq processes, see Eq. (20). Indeed, this suppression is less pronounced at x F = −0.1, where the unpolarized quark and anti-quark TMDs inside the unpolarized proton are probed at larger x values and therefore are not equally important. Moreover, in the small p T range (up to 3 GeV) the maximized estimates at x F = −0.1 are more suppressed w.r.t. those at x F = 0, due, once again, to the integration over the Sivers azimuthal phase. In the lower panels we show our estimates adopting the results discussed in the previous subsection. In all cases the values are very small and compatible with zero. Despite of this, a measure of A N for direct photon production would be extremely important to test the consistency of the whole approach.

IV. CONCLUSIONS
In this paper we have performed a study of the gluon Sivers function through a combined analysis of data on transverse single-spin asymmetries for the processes p ↑ p → π 0 X [12] and p ↑ p → D X → µ X [30], measured by the PHENIX Collaboration at RHIC. The theoretical framework adopted is the so-called transverse momentum dependent generalized parton model (GPM), in which intrinsic parton motion and spin effects are considered. In addition, we have used the color gauge invariant version of this model (CGI-GPM), which takes into account also, in the one-gluon exchange approximation, the initial and final state interactions of the active parton with the remnants of the polarized proton, leading to a process dependent Sivers function.
From a theoretical point of view, we have extended the calculation of the expressions for the single-spin asymmetries in p ↑ p → π X and p ↑ p → γ X, within the CGI-GPM approach, to the gluon sector. In this way, we completed the study of Ref. [21], in which only the corresponding quark-induced subprocesses were studied. As a byproduct, we have also shown that the one-gluon approximation employed here is sufficient to recover the exact gluonic pole strengths in any partonic process calculated at LO in perturbative QCD [25] (see the Appendices).
The analogous formulae for the single spin asymmetries in p ↑ p → D X and p ↑ p → J/ψ X were derived in Ref. [24]. It turns out that for these processes the gluon Sivers function can be re-expressed as a linear combination of two independent, universal (and so far unknown) contributions, namely the f -type and d-type Sivers distributions.
On the phenomenological side, using available knowledge of the quark and antiquark Sivers functions from SIDIS measurements, we have shown how the PHENIX data on inclusive pion and D-meson production allow us to partially disentangle and considerably constrain the size of these two gluon Sivers functions, which should be much smaller than their positivity bounds. This can be considered the first significant attempt towards a quantitative extraction of these process dependent gluon Sivers functions. On the other hand, since the number and the precision of the available data is not very high, our findings have still to be considered as preliminary. Furthermore, we have compared the extractions of the gluon Sivers function in the two approaches, with (CGI-GPM scheme) and without (GPM scheme) initial/final state interactions. The results are encouraging, even if it is not yet possible to clearly discriminate between the GPM and the CGI-GPM frameworks.
Our results have been used to predict the single-spin asymmetry for the processes p ↑ p → J/ψ X, which only depends on the f -type Sivers function. Comparison with existing PHENIX data [31], compatible with zero at forward rapidities, shows a good agreement. Predictions for the same processes have been presented in a kinematic region accessible at LHC with a fixed polarized target, and for the process p ↑ p → γ X at RHIC kinematics as well, for which data are not yet available. These will certainly help in shedding light on the still poorly known gluon Sivers function and towards our understanding of the three-dimensional structure of the nucleons.
Note added: at the very last stage of this work we have become aware of a similar study on SSAs in p ↑ p → γ X within the CGI-GPM approach [45]. While the theoretical findings are in perfect agreement, see Eqs. (19)-(20), the phenomenological analysis presents some differences, which deserve further attention. A possible explanation could be the different way of handling the role of the azimuthal phases (to be integrated over in the final observable) in the hard partonic pieces.
, gives H (f /d) I and H (f /d) Fc , respectively, and (A1) Notice that the C F d factors sum up to zero and do not play any role in the single-inclusive hadron production. Alternatively, H Inc (f /d) can be obtained directly by summing the diagrams with the color factors Finally, we have checked that, for each diagram D, the gluonic pole strengths defined by are in full agreement with the ones given in Ref. [25] for less inclusive processes like p ↑ p → π π X, for which the FSIs of parton d need to be taken into account as well.
where ab → γd is a generic partonic subprocess contributing to p ↑ p → γ X. Our results for the color factors relevant for the gluon induced subprocesses gq → γq and gq → γq are summarized in Table IV. Due to their simple color structures, all diagrams D have the same color factors. As before, C U is the unpolarized one, while C I (C F d ) is the color factor obtained when an extra gluon is attached in D to parton b (parton d). Since the photon does not interact with the remnant of the polarized nucleon, C (f /d) Fc = 0. Finally, we point out that our gluonic pole strengths, defined as are in full agreement with the ones given in Table B.4 of Ref. [46] for gq → γq, namely C (f ) Notice that the results in Ref. [46] have been derived adopting a different method, i.e. by looking at the full gauge link structure and taking the derivative of the gauge link.
For completeness, the hard functions for the quark induced subprocesses, calculated in Ref. [21], are presented below:  IV: Color factors for the process gq → γq. For the process gq → γq, the f -type color factors are the same, while the d-ones have an overall minus sign. Notation is the same as in Tab. I.