Study of the sign change of the Sivers function from STAR Collaboration W/Z production data

Recent data on the transverse single spin asymmetry $A_N$ measured by the STAR Collaboration for $p^\uparrow \, p \to W^\pm/Z^0 \, X$ reactions at RHIC allow the first investigation of the Sivers function in Drell-Yan processes and of its expected sign change with respect to SIDIS processes. A new extraction of the Sivers functions from the latest SIDIS data is performed and a critical assessment of the significance of the STAR data is attempted.


Introduction
The Transverse Momentum Dependent Partonic Distribution Functions (TMD-PDFs) encode information on the 3-dimensional structure of nucleons in momentum space; they depend on the parton intrinsic motion inside the nucleon and, in general, on the nucleon and parton spins. At leading twist there are eight independent TMD-PDFs which have been studied in Semi Inclusive Deep Inelastic Scattering (SIDIS) processes. Among them, the Sivers distribution, which describes the momentum distribution of unpolarised quarks and gluons inside a transversely polarised proton, has a clear experimental signature [1,2] and is of particular interest for several reasons; one expects it to be related to fundamental intrinsic features of the nucleon and to basic QCD properties.
In fact, the Sivers distribution ∆ N f q/p ↑ relates the motion of unpolarised quarks and gluons to the nucleon spin S; then, in order to build a scalar, parity invariant quantity, S must couple to the only other available pseudo-vector, that is the parton orbital angular momentum, L q or L g . Another peculiar feature of the Sivers distribution is that its origin at partonic level can be traced in QCD interactions between the quarks (or gluons) active in inelastic high energy interactions and the nucleon remnants [3,4]; thus, it is expected to be process dependent and have opposite signs in SIDIS and Drell-Yan (D-Y) processes [5,6]: (1.1) This important prediction remains to be tested. The Sivers distribution can be accessed through the study of azimuthal asymmetries in polarised SIDIS and Drell-Yan (D-Y) processes. These have been clearly observed in the last years, in SIDIS, by the HERMES [1], COMPASS [2] and Jefferson Lab [7] Collaborations, allowing extractions of the SIDIS Sivers function [8][9][10][11]. However, no information could be obtained on the D-Y Sivers function, as no polarised D-Y process had ever been measured.
Recently, first data from polarised D-Y processes at RHIC, p ↑ p → W ± /Z 0 X, have become available [12]. The data show an azimuthal asymmetry, A W N , which can be interpreted as due to the Sivers effect and which hints [12,13] at a sign change between the Sivers function observed in these D-Y processes and the Sivers function extracted from SIDIS processes. However, considering the importance of the sign change issue, before drawing any definite conclusion, both the SIDIS and D-Y data and their comparison, have to be critically analysed and discussed.
In this paper we perform a new extraction of the valence and sea-quark Sivers functions from the newest experimental SIDIS data. We then perform an analysis of the RHIC W ± /Z 0 D-Y data [12], based on these new functions, trying to assess the significance of A W N on the sign change of the Sivers functions.
The paper is organised as follows: in Section 2 we recall the formalism used to analyse and interpret the experimental data. In Section 3 we present a new extraction of the Sivers functions from experimental data. In Section 4 we compute the asymmetries observable in D-Y processes and based on the SIDIS extracted Sivers functions, both with and without the sign change, comparing them with the recent RHIC results, while in Section 5 we analyse the impact of the D-Y data as a possible indication of the sign change of the Sivers function. Conclusions and final comments are given in Section 6.

Formalism
We consider a generalised Drell-Yan process, p ↑ p → W ± X, in which one observes a W boson, with four-momentum q, created by the annihilation of a quark and an antiquark. We define our kinematical configuration with the polarised p ↑ proton, with four-momentum p 1 , moving along the positive z-axis and the unpolarised one, with four-momentum p 2 , moving opposite to it. We adopt the usual variables: (2.1) The annihilating quarks have an intrinsic transverse motion, k ⊥1 and k ⊥2 . We fix the azimuthal angles by choosing the "up" (↑) polarisation direction as the positive y-axis (φ S = π/2). The spin "down"(↓) polarisation direction will have φ S = 3π/2. The other transverse momenta azimuthal angles are defined as: In the kinematical region using the TMD factorisation formalism at leading order, the unpolarised cross section for the p p → W X process can be written as [13][14][15][16] dσ pp→W X dy W d 2 q T =σ 0 where f q i /p (x i , k ⊥i ) are the unpolarised TMDs, V q 1 ,q 2 are the weak interaction CKM matrix elements and the q 1 ,q 2 runs over all appropriate light quark and antiquark flavours (q 1 q 2 = ud,du, us,su for W + , etc.).σ 0 is the lowest-order partonic cross section (with G F the Fermi weak coupling constant),σ 5) and the parton longitudinal momentum fractions are given, at O (k ⊥ /M W ), by Notice that, with the definition of x F adopted in Eq. (2.1), one has In such a formalism, the distribution for unpolarised quarks with transverse momentum k ⊥ inside a proton with 3-momentum p and spin S, generates a transverse Single Spin Asymmetry (SSA) where dσ stands for dσ pp→W X /(dy W d 2 q T ) and ∆ N f q/p ↑ (x, k ⊥ ) is the Sivers function. The above expression much simplifies adopting, as usual, a Gaussian factorised form both for the unpolarised distribution and the Sivers functions, as in Ref. [8]: (2.12) (2.14) where f q (x) are the unpolarised PDFs, M 1 is a parameter which allows the k ⊥ Gaussian dependence of the Sivers function to be different from that of the unpolarised TMDs and N q (x) is a function which parameterises the factorised x dependence of the Sivers function.
The following moment of the Sivers function is of importance: (2.17) With the choices of Eqs. (2.12)-(2.15) the k ⊥ integrations can be performed analytically in Eq. (2.11), obtaining: and where, in the last line, we have used, according to our kinematics, S·(p 1 ×q T ) = cos φ W . A N (y W , q T ) is the quantity measured at RHIC [12] 1 . Let us notice that the RHIC measurements of W ± production at √ s = 500 GeV [12] cover the rapidity region |y W | < 1. In particular, data are available for y W ± 0.4 and y W 0. This corresponds to: where x 1 refers to the polarised proton and x 2 to the unpolarised one. Then, although the x region is predominantly the valence one, the data at y W −0.4 are expected to be more sensitive to the sea-quark Sivers functions.

Extraction of Sivers functions from SIDIS data
The quark flavours involved in W production include anti-quarks. Thus, in order to estimate the asymmetry A W N , it is important to have a reliable extraction of both quark and antiquark Sivers functions.
For instance, in order to produce a W + , u,d ands quarks from the polarised proton combine withd,s, u quarks from the unpolarised proton, such that the asymmetry is proportional to (3.1) Both quantities in the round brackets in the above equation contain a sea and a valence quark distribution. However, because of the numerical values 2 of |V u,d | and |V u,s |, the last two terms in Eq. (3.1) are much suppressed with respect to the first two. Thus, we expect that A W + N mainly depends on the u quark andd sea quark Sivers functions. Likewise, for W − production, the asymmetry is proportional to 2) and we expect that W − data are mainly sensitive to d quark andū sea quark Sivers function.
A previous extraction of the Sivers functions that included anti-quark distributions was reported in Ref. [8]. However, new data have become available since then and we perform here a new complete extraction of the Sivers functions. We refer to Ref. [8] for more details about the procedure.
One may notice that in our simple parameterisation of the Sivers functions as given in Eqs. (2.12)-(2.15) the knowledge of the width k 2 ⊥ of the unpolarised TMDs is important. Such a study was performed in Refs. [18,19]. We adopt here the parameters from Ref. [18], fixed by fitting the HERMES multiplicities [20]: where p 2 ⊥ is the width of unpolarised Transverse Momentum Dependent Fragmentation Functions (TMD-FFs): Notice that the study of Ref. [18] found no flavour dependence of the widths of the TMDs. The collinear distribution and fragmentation functions, f q/p (x) and D h/q (z), needed for our parameterisations are taken from the available fits of the world data: in this analysis we use the CTEQ6L set for the PDFs [21] and the DSS set for the fragmentation functions [22]. The LHAPDF [23] library is used for collinear PDFs. We fit the latest data from the HERMES Collaboration on the SIDIS Sivers asymmetries for π ± and K ± production off a proton target [1], the COMPASS Collaboration data on LiD [24] and NH 3 targets [25], and JLab data on 3 He target [26]. These available SIDIS data cover a relatively narrow region of x, typically in the socalled valence region. It suffices to use the most simple parameterisation for the anti-quark Sivers functions [see Eqs. (2.13), (2.14)]: It means that we assume the anti-quark Sivers functions to be proportional to the corresponding unpolarised PDFs; we have checked that a fit allowing for more complicated structures of Eq. (2.14) for the anti-quarks, results in undefined values of the parameters α and β. The Sivers asymmetry measured in SIDIS can be expressed using our parameterisations of TMD functions from Eqs. (2.12-2.15, 3.4) as Thus, we introduce a total of 9 free parameters for valence and sea-quark Sivers functions: N uv , N dv , Nū, Nd, α u , β u , α d , β d , and M 2 1 (GeV 2 ). In order to estimate the errors on the parameters and on the calculation of the asymmetries we follow the Monte Carlo sampling method explained in Ref. [8]. That is, we generate samples of parameters α i , where each α i is an array of random values of {N uv , N dv , Nū, Nd, α u , α d , β u , β d , M 2 1 }, in the vicinity of the minimum found by MINUIT, α 0 , that defines the minimal total χ 2 value, χ 2 min . We  generate 2 · 10 4 sets of parameters α i that satisfy with the high tolerance ∆χ 2 = 17.21 that corresponds to the 95% C.L. of coverage probability for 9 free parameters. The fit is performed with MINUIT minimisation package and the resulting parameters can be found in Table 1; the corresponding extracted Sivers functions are shown in Fig. 1. We indicate both the errors for the standard definition of ∆χ 2 = 1 and the high tolerance error with ∆χ 2 = 17.21 (the errors given in parentheses).
The main new features of the fit are the parameters N dv = −0.52 ± 0.20 and N uv = 0.18 ± 0.04. The previous extraction [8], that used different gaussian width values, k 2 ⊥ = 0.25 GeV 2 and p 2 ⊥ = 0.20 GeV 2 , yielded N d = −0.9, which almost saturated the positivity bound |N q | = 1, and N u = 0.35. Theū andd Sivers functions turn out to be both small, compared to the quark distributions, and negative. Future Electron-Ion Collider data will be crucial for the investigation of the anti-quark Sivers distributions. The parameters that control the large-x behaviour of the functions, β uv and β dv , have big errors, see Table 1. The future Jefferson Lab 12 GeV data will allow a better precision extraction in the high-x region.
The partial contributions to χ 2 from different experiments are shown in Table 2. One can see that the proton data on π + from the HERMES Collaboration and the positive hadron data from the COMPASS Collaboration show some larger χ 2 values that might be attributed to possible effects of TMD evolution [10,11,27].
Several plots showing the quality of our best fits of the data are presented in Fig. 2

Predictions for W and Z asymmetries and comparison with data
We can now compute the asymmetry A N (y W , q T ), according to Eqs. (2.18)-(2.19), using the Sivers functions -or their opposite -as given in Eqs. (2.12)-(2.15) with the parameters, and the corresponding uncertainties, shown in Table 1.
Actually, in order to compare with data [12], we integrate both the numerator and denominator of A W N , Eq. (2.11), either over q T in the region [0.5, 10] GeV, or over y W from −1 to 1. The results, reversing the sign of the SIDIS extracted Sivers functions as in Eq. (1.1), are shown and compared with data respectively in Fig. 3 and in Fig. 4. For completeness, despite the much limited amount and quality of data, we also show our estimate of A N , integrated over q T , for Z 0 production, in Fig. 5. The results without the sign change can be easily deduced by reversing the sign of the asymmetry in Figs. 3-5. Before trying, in the next Section, a quantitative evaluation of the significance of the data regarding the issue of the sign change of the Sivers function going from SIDIS to D-Y processes, a few comments are in order.
• In general, the agreement between our estimates and the few data is rather poor, both with and without sign change. In particular, this is evident from the q T dependence  of A N , Fig. 4, and the y Z dependence of A N for Z 0 , Fig. 5. In the latter case there is only one single data point, with a big error, indicating a large positive asymmetry.
• The data on the y W dependence are given by collecting all W 's produced with q T up to 10 GeV. The simple model of D-Y TMD factorisation without evolution that we use in this analysis is expected to hold for lower values of q T ; integrating the theoretical results up to such values, in order to compare with the available data, is a somewhat ambiguous procedure. Implementation of the TMD evolution would not help to make the agreement with the data better in this case, as TMD evolution predicts a suppression of the asymmetries for higher values of Q 2 with respect to the initial lower scale [11]. This suppression might become moderate depending on the shape of the non-perturbative input of TMD evolution [28][29][30]. • Considering the q T integrated data, from a first look at Fig. 3 it appears that indeed W − data are compatible with the sign change, while W + data may be compatible with either sign of the Sivers functions.
• The shape of the TMDs and the values of the parameters here adopted allow a good description of the SIDIS data; however, they are still rather flexible, and our numerical estimates for the D-Y asymmetry might depend on the choice, for example, of the values of the Gaussian width, Eq. (3.3). A full study of combined unpolarised SIDIS, D-Y and (future) e + e − data is mandatory.

Impact of the asymmetries on the extraction of Sivers functions
In this Section we take at face value the RHIC data on A N for W ± production and, in order to quantify their significance, we calculate the deviation between the data and our estimates, separately for W + and W − : where [theory] n (α) corresponds to the calculation of the W asymmetry using the phenomenological extraction of the Sivers function performed in this paper, with model parameters α, with and without the sign change of Eq. (1.1); [exp] n are the data for W + or W − asymmetries and [∆exp] n are the corresponding experimental errors. As we explained in Section 3, in order to estimate the error on the extraction of the Sivers functions, we generate 2·10 4 sets of parameters α according to Eq. (3.7). Thus, we calculate 2·10 4 values of χ 2 using Eq. The green histogram corresponds to χ 2 with no sign change of the Sivers function, while the blue histogram corresponds to χ 2 with the sign change of the Sivers functions.
One can see from the upper left panel of Fig. 6 that W − data favour the sign change: in this case the values of χ 2 /dof are around 1.1, while without the sign change they are around 2.7. The W + data on the other hand are slightly better with no sign change, as can be seen from the upper right panel of Fig. 6. For either scenarios the χ 2 per number of data are rather large: these large values are due to the single point at y W = 0 (see Fig. 3, left panel) and the two points at large q T > 5 GeV (see Fig. 4).
If we combine both W + and W − data, then the two data sets globally favour a sign change of the Sivers functions according to Eq. (1.1). The histogram of the combined data sets is presented in the lower panel of Fig. 6. If one assumes no sign change, then χ 2 /dof = 2.35 and σ(χ 2 /dof) = 0.1, where dof = 16, while the sign change yields a lower value, χ 2 /dof = 1.75 and σ(χ 2 /dof) = 0.05. Notice that both scenarios have some disagreement with our estimates: indeed the values of χ 2 /dof are well above one. Using our results from Fig. 6 we can at most conclude that W ± data hint at an indication of the sign change according to Eq. (1.1).
Another interesting question that we would like to investigate in this paper is whether the data on the W ± asymmetries have any significant impact on the parameters of our model. Notice that we do not include W ± data in our fit of the Sivers functions. Bayes theorem allows to incorporate information from new data by applying a re-weighting of the probability densities for the model parameters. The details of application of the reweighting are explained in Ref. [31]. The probability density function for model parameters α, P(α), is going to be modified in presence of new data, and the Bayes theorem states that where P(α|D) is the so-called posterior density, that is the updated probability density function from the prior density P(α). The quantity P(D|α), called the likelihood function, represents the conditional probability for a data set D given the parameters α of the model. The quantity P(D) ensures the normalisation of the posterior density to unity. For a particular observable O one can write the expectation value with the new data as, In order to estimate the integral of Eq. (5.3) we will use a Monte Carlo approximation of the integral, such that the integral over continuous values of α will be substituted with a sum over discrete values of α k . These α k are the generated 2·10 4 sets of parameters according to Eq. (3.7). We obtain where the quantities w k are called weights and are proportional to P(D|α k ). Their normalisation is fixed by demanding E[1] = 1, that is, k w k = N . Similarly, the variance of an observable O is given by The re-weighting procedure depends on the form assumed for the likelihood function. We use χ 2 minimisation in the fits, then our weights have the following form: where all values of χ 2 (α k ) can be readily obtained from our results and Eq. (5.1). We have checked that the form of the prior density P(α) for our parameters turns out to be very well approximated by the normal distribution. Initial Reweighted Figure 7. The re-weighting procedure applied to N uv , Nd, N dv , and Nū parameters. The green hatched histogram corresponds to re-weighting, while the blue histogram corresponds to the prior distribution.
The re-weighting can be readily applied to all parameters of the model, as well as to observables: in Fig. 7 we show, as an example, the application of the re-weighting procedure to the parameters that describe the normalisation of the u-valence, d-valence, anti-d, and anti-u quark Sivers functions. One can see that the re-weighted densities (green hatched histograms) are only slightly shifted from the prior distributions (blue histograms), and the new mean value is within 1σ of the current values, given in Table 1. The same observation is true for all other parameters.

Comments and conclusions
We have analysed the recent data on the single spin asymmetry A W N measured by the STAR Collaboration at RHIC [12]; it is the first ever spin asymmetry measured in Drell-Yan processes and it might originate from the fundamental Sivers distribution of polarised quarks in an unpolarised proton. Then, it could help in testing the validity of the widely expected sign change of the Sivers function when extracted in SIDIS and D-Y processes, Eq. (1.1).
In order to perform an unbiased analysis we have re-derived, by best fitting the latest SIDIS data [1,[24][25][26], the Sivers functions, including the anti-quark ones which might play a role in the D-Y production of W s and Z 0 s. Our results are shown in Tables 1 and 2 and in Fig. 1.
Using the newly extracted Sivers SIDIS functions we have computed the D-Y SSA A N for W ± and Z 0 production, both with and without a sign change of the Sivers functions. Then, we have compared our results with the STAR data, trying to assess their significance with respect to the sign change issue. Our quantitative results, according to Eq. (5.1), can be seen in Fig. 6.
As commented throughout the paper, our simple model of D-Y TMD factorisation without evolution, Eqs. (2.9)-(2.15), is, in general, in poor agreement with the data. A more refined analysis, using the TMD evolution, would probably worsen the agreement [11]. One should add that the data, although important and pioneering, are still scarce, with large errors, and gathered in different kinematical regions.
With all the necessary caution, from our analysis of the data, one can at most conclude that, only from W − production, there is an indication in favour of the sign change of the Sivers function, which, however, is still far from being considered as proven. Soon expected data from COMPASS polarised D-Y processes, π − p ↑ → + − X, and higher statistics data from STAR Collaboration on W and Z production should add important information.