Towards a first estimate of the gluon Sivers function from $A_N$ data in $pp$ collisions at RHIC

Within a generalized parton model approach, with inclusion of spin and intrinsic transverse momentum effects, we show how the latest, highly precise, midrapidity data on the transverse single spin asymmetry measured in $pp\to\pi^0\, X$ by the PHENIX Collaboration at RHIC [1], can be used to get a first estimate on the still poorly known gluon Sivers distribution. To this end we also adopt the present information on the quark Sivers functions, as extracted from semi-inclusive deeply inelastic scattering data. This analysis updates a previous study by some of us where a first bound on this distribution was obtained [2].

studying the Sivers asymmetry in the production of high-p T hadron pairs in muon scattering off polarized proton and deuteron targets [64]. This process should be dominated by the photon-gluon fusion mechanism and therefore allows to get information on the GSF. First results gave an asymmetry compatible with zero for deuteron target at x G = 0.13. This fact, together with additional theoretical considerations, led Brodsky and Gardner [65] to state that the gluon contribution to parton orbital angular momentum (and the GSF) should be negligible. However, very recent preliminary measurements of the same observable for proton target give a negative gluon Sivers asymmetry, −0.26 ± 0.09 ± 0.08 at x G = 0.15 [64]. This value even if 3σ below zero is still compatible with the deuteron result.
From the phenomenological point of view, it has been suggested to study the role of the GSF in polarized protonproton collisions in several processes: SSAs in inclusive photon production in the large negative x F region (measured w.r.t. the polarized proton) [66]; back-to-back azimuthal correlations in two-jet production [67]; SSAs in inclusive D meson production at RHIC [68]; SSAs in J/ψ electroproduction with transversely polarized electron and proton beams [69]. A detailed and updated discussion on the gluon Sivers function and additional references can be found in Ref. [70].
Besides the gluon Sivers function, the role of linearly polarized gluons inside (un)polarized protons in inclusive processes in proton-proton collisions has been also actively investigated in recent years, e.g. in pion-jet production [71], heavy quark and jet-pair production at electron-ion or hadron colliders [72,73], and Higgs production at the LHC [41,74,75].
Assuming the validity of the TMD formalism for a single-scale process we show here how the analysis of highly precise midrapidity data in single polarized pp → πX processes could strongly constrain the gluon Sivers function. As shown in a series of papers [20,51,53,[76][77][78] this phenomenological approach, nowadays known as generalized parton model (GPM) is able to describe fairly well many features of several available data for such a process and it is worth to be further investigated. Even if not supported, as already said above, by a formal proof of factorization in terms of TMDs, the study of these processes can be very useful also in clarifying the role of process dependence and factorization breaking effects.
It has been already shown that, within a TMD scheme, due to strong partonic azimuthal phase cancellations, the backward hemisphere can be of little use to get information on polarized TMDs, since all effects are almost washed out [78]. On the other hand, as shown by some of us in Ref. [2], a study of midrapidity A N data at high energy can be used to constrain the gluon Sivers function. Here we will upgrade this result by using more recent information both on the phenomenological and the experimental side.
Indeed we have now at our disposal phenomenological extractions of the quark Sivers functions from SIDIS processes (not available at that time), one of them including also the sea quark contributions [45,50]. At the same time new and highly precise data from the PHENIX Collaboration at RHIC have been made available [1]. For these reasons we believe that such a reanalysis is timely.
As said, the issue of QCD evolution of TMDs, strongly related to factorization, is still an open question for such single scale processes. Concerning its potential role in the following analysis, we believe that the relatively modest range of pion transverse momentum involved (at least for the more precise data which can significantly constrain the GSF) prevents it to be effective for the asymmetries. Therefore, in the sequel we will keep including QCD evolution only in the collinear factorized component of the involved TMDs (see also below for more details).
Although in the TMD approach several terms may in principle contribute to the single spin asymmetry A N (p ↑ p → π X) ≡ (dσ ↑ − dσ ↓ )/(dσ ↑ + dσ ↓ ), in the kinematical regime of Ref. [1], as extensively discussed in Ref. [2], A N is largely dominated by the Sivers effect alone, and its numerator is given by (for details see Refs. [77,78]) where (M p denotes the proton mass) ∆ N f a/p ↑ (x a , k ⊥a ) [or f ⊥a 1T (x a , k ⊥a )] is referred to as the Sivers distribution function of parton a inside a transversely polarized proton [79]. φ a is the azimuthal angle of the intrinsic transverse momentum k ⊥a of parton a. For details and a full explanation of the notations in Eq. (1) see Ref. [78]. It suffices to notice here that J(k ⊥π ) is a kinematical factor, which at O(k ⊥π /E π ) equals 1 and dσ ab→cd /dt is the partonic differential cross section for the subprocess ab → cd.
Notice that the parton a inside the polarized proton can be a quark (or an antiquark) and a gluon, that is the Sivers contribution to the asymmetry can be expressed as a sum of two terms, that cannot be disentangled in this process. For this reason in the following analysis we will take into account all available information on the quark Sivers functions. In particular we consider two extractions. The reason for this is twofold: from one side, in the first extraction (SIDIS1 in the following) [45] only u and d flavours were considered, while in the second one (SIDIS2) [50] also the sea quark contributions were accounted for; secondly, and somehow more relevant, in SIDIS1 the set of fragmentation functions of Kretzer [80] was adopted, while in SIDIS2 the set of de Florian, Sassot and Stratmann (DSS) [81] was considered, which provides a much more important and very different leading order gluon fragmentation. This aspect, as shown in the following, could play a non negligible role in the present study.

II. ANALYSIS AND RESULTS
The presently available data on A N (p ↑ p → π 0 X) by the PHENIX Collaboration [1] are extremely precise, of the order of per mil, and with tiny errors, in particular in the region of moderate P T , where the gluon initiated processes dominate. In fact, both their central values and the error bars are at least one order of magnitude smaller than the data analysed in Ref. [2]. For this reason, while in that work a first, very conservative, upper bound on the gluon Sivers function was presented without entering into a more detailed analysis, here we want to present a more careful study aiming at a first tentative estimate of the GSF within a TMD scheme.
As stated above, in the kinematical region considered only the Sivers effect can play a relevant role, being all the other effects suppressed by strong azimuthal phase cancellations, as discussed in Ref. [2]. If one adopts the more detailed SIDIS2 parameterizations of Ref. [50], where also the sea quark Sivers functions were considered, one would find a contribution to A N compatible with zero. Taking into account that the data [1] are also almost compatible with zero, one would conclude that there is no room for the gluon Sivers effect.
In the spirit of Ref. [2] we adopt again a conservative attitude and investigate to what extent, taking into account the uncertainty on the quark Sivers distributions together with the small errors on the data, a gluon Sivers function could still play a role.
Even if the small number of data points available (ten in this case) does not allow a full statistical analysis, namely a fit, we still try to substantiate our study adopting a commonly used statistical criterium. We then define a proper χ 2 function and, using information available both on the experimental and the phenomenological side, we extract a parametrization of the gluon Sivers function by minimizing it. Since we aim at extracting the contribution to the Sivers effect from gluons, using the available phenomenological information on the corresponding contribution from quarks with its uncertainty, we define where the sum runs over the data points, σ exp is the experimental error on A exp N and σ quark the estimated theoretical uncertainty on the quark contribution A quark N , considered as a known quantity. In such a way we can constrain the gluon contribution taking into account both the theoretical (even if partially) and the experimental uncertainties.
Concerning the gluon Sivers function we adopt a somehow standard factorized functional form, analogous to the quark case [45,50], namely: where f g/p (x) is the standard unpolarized collinear gluon distribution, with |N g | ≤ 1, and With these choices, assuming that the unpolarized TMD gluon distribution is given by the Sivers function automatically fulfils its proper positivity bound for any (x, k ⊥ ) values. Consistently, for the unpolarized TMD fragmentation function (for a parton c) we use In the following, for the Gaussian width of the unpolarized TMD gluon distribution we use the same value as for the quark distribution, that is k 2 ⊥ = 0.25 GeV 2 [45]. Moreover, we define the parameter so that the k ⊥ -dependent part of the Sivers function becomes From Eq. (10) it is clear how the range of variation of the parameter ρ is 0-1. In the following analysis we will also consider ρ as a free parameter, another improvement w.r.t. Ref. [2], where a fixed value of ρ was adopted in order to maximize the gluon Sivers effect. We will then minimize the χ 2 function defined in Eq. (4) in terms of the following four parameters: N g , α, β entering Eq. (6) and ρ in Eqs. (10), (11). Since the integration in Eq. (1) is over an eight-dimensional phase space, the fit procedure over the continuous parameter phase space would be quite CPU-time consuming. Therefore, we scan the 4-dimensional parameter space over a discrete grid of values, fine enough for our purposes. More precisely, we consider the following ranges: −1 ≤ N g ≤ 1 (step value of 0.05), 0 ≤ α, β ≤ 4 (step value of 0.2), while for the ρ parameter we consider five representative values: 2/3 (as adopted in Ref. [2] maximizing the effect of the GSF), the same value as for the quark Sivers function and three more values, 0.2, 0.3, and 0.8 (lower or larger values would spoil the description of data).
As stated in the introduction, we will consider the results based on the more recent DSS-SIDIS2 parameterization [50] as well as those obtained adopting the KRE-SIDIS1 set [45], being quite representative extractions of the quark Sivers functions. In both cases, for consistency, we will adopt the GRV98-LO set [82] for the unpolarized parton distributions.
As a first step we checked that the unpolarized cross sections in the same kinematical regime, that is √ s = 200 GeV and central rapidity, can be reproduced adopting the TMD distributions and fragmentation functions discussed above. This is an important issue since these quantities enter as denominators in A N .
After that, we performed our χ 2 minimization over the discretized parameter phase-space. The best (total) χ 2 value obtained is χ 2 min = 1.93 for the SIDIS2 set and 1.86 for the SIDIS1 set. Interestingly, the corresponding best value of ρ, in both cases, is equal to the corresponding one obtained for the quark Sivers function.
Notice that since the parameters are quite correlated among them, many sets in the explored grid give χ 2 's very close to the minimum value and therefore comparable estimates (see below for a discussion on the uncertainties). An important remark is that about half of the best χ 2 value comes from the largest-P T data point, which has a very large error bar and is less sensitive to the gluon distributions (largest x). Exclusion of this point would give a total minimum χ 2 of about 1.
For completeness, even if this is not the main aim of our study, and taking them with a grain of salt (see previous comments), we give the best-fit parameter sets: Bearing in mind the caution raised above, it is nevertheless interesting to note that in both cases the gluon Sivers function turns out to be positive. This is another improvement w.r.t. the previous study [2] where no information on its sign could be extracted.
In the spirit of being conservative and including potential sources of uncertainties, we will consider also two possible uncertainty bands, generated respectively by the envelope of the A N values obtained adopting all parameter sets in the parameter-space grid leading to an increase in the χ 2 value of 2% and 10% w.r.t. the χ 2 min , that is ∆χ 2 = (2−10%) χ 2 min .  We notice here that a tolerance of 5% would give results very similar to the 10% uncertainty band. As stated above, given the limited number of experimental data, we cannot claim to have a statistically significant best fit. Therefore, it would not make sense defining and showing statistical error bands. On the other hand, it is useful to quantify the level of accuracy in the description of the data and the corresponding gluon Sivers function when the χ 2 varies within these ranges.
In Fig. 1 we present our results for A N (quark plus gluon contributions) at √ s = 200 GeV and midrapidity, compared with PHENIX data [1] and adopting the SIDIS2 [50] extraction of the quark Sivers functions. Here we show the full P T range, together with our best estimate (solid red line) and a red band corresponding to a tolerance of 10% in χ 2 , as explained above. As one can see the description of data is extremely good, even if the scale adopted in the plot and the tiny data values hide some details. Almost undistinguishable results are obtained for the SIDIS1 set.
To better visualize the data description and the differences between the two sets, in Fig. 2 we show the results for A N in the lower P T range for the SIDIS2 (left panel) and the SIDIS1 (right panel) sets. Quite importantly, this is the region that better constrains the gluon Sivers contribution. In this case, we also show the narrower tolerance green band corresponding to a 2% increase in χ 2 , together with the contribution coming from the best estimate of the gluon Sivers function (blue dotted line).
Notice that for the full-P T range, the Bjorken x explored varies, roughly, between 6 · 10 −3 and 0.6, while in the lower-P T range (up to 5 GeV) the maximum value of x is around 0.4-0.5. This has to be taken into account, together with the fact that the adopted quark Sivers functions are constrained by available SIDIS data only in the region x DSS-SIDIS2 x KRE-SIDIS1 ∆χ 2 /χ 2 min =10% ∆χ 2 /χ 2 min = 2% Ref. [2] FIG. 3: First k ⊥ -moment of the gluon Sivers function as defined in Eq. (14) for the SIDIS2 set (left panel) and SIDIS1 set (right panel) at Q 2 = 2 GeV 2 . The best estimates (red solid lines) are shown together with the tolerance bands corresponding to a 2% (narrower, green) and 10% (wider, red) variation in the χ 2 . The former bound on the gluon Sivers function (magenta dotted line), obtained in Ref. [2], is also shown.
up to x ∼ 0.3. In other words, the present analysis, which aims at constraining the gluon Sivers function adopting the information on the quark Sivers contribution and the midrapidity data in pp collisions, is sound only up to x ∼ 0.3 − 0.4. On the other hand, this is the most interesting region for a study of gluon distributions. In Fig. 3 we present the corresponding results for the first k ⊥ -moment of the gluon Sivers function, defined as More precisely we show (SIDIS2 set in the left panel and SIDIS1 set in the right panel) the best estimates, red solid line, together with the two tolerance bands of 2% (green, the narrower one) and 10% (red, the wider one) and the previous upper bound obtained in Ref. [2] (magenta dotted line). Notice that the two results (old vs. new bound) for both sets are not directly comparable due to the deep differences in the two analyses. Nevertheless from this new study one can appreciate the tiny role left to the gluon Sivers function when one tries to describe the latest A N data at midrapidity. This is confirmed even assuming a relatively large tolerance in χ 2 , like those considered here. From these results one can quantify the role played by the indeterminacy on the quark Sivers functions and on the fragmentation function sets. This is definitely an important source of uncertainty in the GSF extraction. In particular, as shown in Fig. 3, we see that the GSF is much smaller (but with larger uncertainties) for the KRE-SIDIS1 case in the low x region, while on the contrary is more constrained for the DSS-SIDIS2 case in the large-x region. For 0.05 < ∼ x < ∼ 0.3 the two extractions are almost compatible, considering the uncertainty bands, with the DSS-SIDIS2 bands narrower than the corresponding bands for the KRE-SIDIS1 set in the low-x region, while the viceversa is true in the large-x region. This is related to the fact that the SIDIS2 set has also a more constrained sea quark component and the DSS fragmentation set enhances the role played by the gluon distribution.
Another potential source of uncertainty of this analysis is related to the direct use of the quark Sivers functions as extracted from SIDIS data. As already remarked in the introduction, we do not have a proof of TMD factorization and universality of TMD PDFs for inclusive processes like the one under consideration. Nevertheless, one could speculate about the possible impact of initial and/or final state interactions.
A way to implement these effects in pp → π X processes was proposed few years ago in Ref. [83], and applied to inclusive pion jet production in pp collisions in Ref. [84], in the framework of the so-called color gauge-invariant (CGI) GPM approach. We recall here that the authors of these works focused only on the quark initiated processes and that nothing has been done so far on the gluon sector. To account also for this source of uncertainty we have reconsidered the contribution of the quark Sivers functions adopting the CGI-GPM. It is important to note that differently from what happens in the forward rapidity region, where one gets a contribution of almost the same size but opposite in sign w.r.t. the GPM, in the midrapidity region the overall effect is a strong reduction in size, but keeping the same sign. This is due to the fact that in this kinematical region many partonic channels play a comparable role, leading to relative cancellations among their contributions. The use of this result in the present analysis would imply a reduction of the GSF and a relative larger indeterminacy towards its smaller values.
In the spirit of further pursuing this issue we explore a somewhat more extreme scenario, maybe less realistic, but worth of being considered. We repeat the procedure described above, adopting a quark Sivers function reversed in  [50] for the quark Sivers functions reversed in sign. The red(green) band represents a tolerance of 10%(2%) in χ 2 (see text for details). The gluon contribution to AN , blue dotted line, is also shown. Right panel: First k ⊥ -moment of the gluon Sivers function as defined in Eq. (14) for the SIDIS2 set (with the quark parameterizations reversed in sign) at Q 2 = 2 GeV 2 . The best estimate (red solid line) is shown together with the tolerance bands corresponding to a 2% (narrower, green) and 10% (wider, red) variation in the χ 2 . The former bound on the gluon Sivers function (magenta dotted line), obtained in Ref. [2], is also shown. sign w.r.t. the one extracted from SIDIS, even if we are aware that for the process under consideration one would expect a more involved structure. We think that this attempt should give a clear indication of the most extreme variation of the GSF uncertainty.
In Fig. 4 we show the results for the SIDIS2 case, that corresponds to the most striking effect. In the left panel, one can easily see that the description of A N is now given in terms of the quark Sivers function alone and that the GSF contribution is almost negligible. This reflects also in the first moment of the GSF, right panel. Here, beside its stronger suppression (red solid curve) compared to the GPM result (left panel of Fig. 3), the uncertainty bands, even for a tolerance of 2%, extend to very low values, compatible with zero. Analogous considerations concerning the stronger suppression of the GSF in the effectively explored region apply also to the SIDIS1 case.
It is evident that such an extreme scenario would imply an even more negligible role of the GSF, that could result almost compatible with zero.
In summary, and with all the cautions discussed above, in the x region explored by pp data the gluon Sivers function can be effectively constrained. It results to be positive, strongly suppressed with respect to the previous bound [2] and much smaller than its positivity bound.

III. CONCLUSIONS
In this paper we have analyzed the impact of recent, highly precise, data for the transverse single spin asymmetry A N (p ↑ p → π 0 X) at central rapidity and moderately large transverse momentum measured by the PHENIX Collaboration at RHIC on our knowledge of the gluon Sivers function.
To this aim we have utilized the so-called transverse momentum dependent generalized parton model which takes into account intrinsic parton motion and spin effects, extending the well-known collinear leading order parton model.
Adopting the most recent phenomenological information, within the same approach, on the (sea) quark Sivers distributions, coming from SIDIS data, we have shown how the PHENIX data allow us to constrain the GSF considerably, as compared to the positivity bound as well as to a previous bound based uniquely on less precise A N data.
We have found that the new constraint is particularly significant, within theoretical uncertainties, in the region of gluon momentum fraction 0.05 < ∼ x < ∼ 0.3, that is the presently explored SIDIS region, where the quark Sivers distributions are well constrained. At lower-x values the bound is still effective but theoretical uncertainties become large, as expected. At larger-x values the bound is looser; on the other hand this is the region where the relevance of gluon contributions is small due to dominance of quark channels.
We have also considered midrapidity data measured by the STAR collaboration for pp → jet X processes [85], where one can access directly the TMD-PDFs. On the other hand these data do not improve the constraint on the GSF due to the P T region explored and their relatively large error values. We have nevertheless checked that the new bound is consistent with these data.
We can then say that the present analysis, constraining the GSF in size and sign (a new aspect w.r.t. Ref. [2]), strongly reduces the possible role of the GSF in spin and azimuthal asymmetries for processes covering x regions similar to those explored here.
Some words of caution are however required: factorization has not been proven in the context of the TMD-GPM for single inclusive processes in proton-proton collisions. In our approach TMDs keep their partonic interpretation and universality. However, initial and/or final state interactions, required to preserve color gauge invariance, might spoil the factorization for these processes, leading to process dependence and universality breaking effects. It is not easy to figure out only from theoretical considerations the phenomenological relevance of these possible effects for currently accessible processes. As a matter of fact, nowadays the TMD-GPM model is able to reproduce fairly well, within uncertainties, the majority of experimental data available on azimuthal and single spin asymmetries in SIDIS and proton-proton collision processes. Moreover, from presently available data there is no clear and unambiguous evidence of sizable universality breaking effects.
However, to investigate, even in an approximate way, the potential role of the process dependence of the quark Sivers functions, we have also considered an extreme scenario, adopting the quark Sivers function as extracted from SIDIS but reversed in sign. In such a case in the explored region one would get an even more suppressed GSF, with uncertainty bands extending to values compatible with zero.
For these reasons, we believe that, even if with some caution, the constraints on the GSF resulting from this analysis are sound within uncertainties and must be taken into account in further phenomenological analyses involving such TMD distribution.
Less involved processes from the point of view of color gauge links, where proving factorization might be easier, need to be considered to further clarify the issues concerning process dependence and test the smallness of the GSF. For example, processes like ep ↑ → QQ X, or ep ↑ → jet jet X, where only final state interactions (like in SIDIS) are involved, could be studied at future electron-ion colliders (EIC) [86]. Analogously, in p ↑ p → γγ X [87] or p ↑ p → J/ψ γ X processes only initial state interactions are involved like in DY. Moreover, factorization has been already proved, at the next-to-leading order, for p ↑ p → η c,b X [88] and A N for this process might be measurable at the proposed AFTER@LHC set-up [89,90].
In Refs. [83,84] a first attempt to compute initial and final state interactions for SSAs in hadronic collisions within the TMD-GPM approach, focusing on quark initiated subprocesses, and studying their phenomenological consequences, has been made. It would be very interesting to extend this kind of analysis to gluon initiated subprocesses, relevant to the present case, and study its effects on the bounds for the gluon Sivers function.