Extraction of the Sivers Function from SIDIS, Drell-Yan, and W =Z Data at Next-to-Next-to-Next-to Leading Order

We perform the global analysis of polarized semi-inclusive deep inelastic scattering (SIDIS), pioninduced polarized Drell-Yan (DY) process, and W =Z boson production data and extract the Sivers function for u, d, s and for sea quarks. We use the framework of transverse momentum dependent factorization at next-to-next-to-next-to leading order (NLO) accuracy. The Qiu-Sterman function is determined in a model independent way from the extracted Sivers function. We also evaluate the significance of the predicted sign change of the Sivers function in DY process with respect to SIDIS.

Introduction.-The three-dimensional (3D) hadron structure is an important topic of theoretical, phenomenological, and experimental studies in nuclear physics. In the momentum space, the 3D nucleon structure is described in terms of transverse momentum dependent distributions and fragmentation functions, collectively called TMDs, which depend both on the collinear momentum fraction and the transverse momentum of parton. The TMD factorization theorem [1][2][3][4][5][6][7] provides a consistent operator definition and evolution of TMDs. Among TMDs, the Sivers function f ⊥ 1T ðx; k T Þ [8,9] is the most intriguing since it describes the distribution of unpolarized quarks inside a transversely polarized nucleon and generates single-spin asymmetries (SSAs).
The Sivers function arises from interaction of the initial or final state quark with the remnant of the nucleon and thus, many of its features reveal the gauge link structure that reflects the kinematics of the underlining process [10]. Above all, the difference between initial and final state gauge contours leads to the opposite signs for Sivers functions in semi-inclusive deep inelastic scattering (SIDIS) and Drell-Yan (DY) kinematics [11][12][13], In the limit of the large transverse momentum the Sivers function is related [14] to the key ingredient of collinear factorization of SSAs, the Qiu-Sterman (QS) function [15][16][17][18], which describes the correlation of quarks with the null-momentum gluon field. Therefore, the measurement of the Sivers function and the exploration of its properties is a crucial test of our understanding of the strong force, and one of the goals of polarized SIDIS and DY experimental programs of future and existing experimental facilities such as the Electron Ion Collider [19,20], Jefferson Lab 12 GeV Upgrade [21], RHIC [22] at BNL, and COMPASS [23,24] at CERN. In this work, we perform the global analysis of transverse spin asymmetries at next-to-next-to-next-to leading order (N 3 LO) perturbative precision within the TMD factorization approach and extract Sivers function. Several important features make our results stand out from the previous results [25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40]. First of all, we use unprecedented N 3 LO perturbative precision, together with the ζ prescription [41]. Second, we use unpolarized proton and pion TMDs extracted from the global fit of SIDIS and DY data [42,43] at the same perturbative order and scheme, which allows us for the first time to describe with very good quality, SIDIS, DY, and W AE =Z experimental data. Lastly, we use the novel model-independent approach to obtain QS function from the Sivers function. Also, we estimate the significance of the sign flip relation, Eq. (1).
The variable P hT is the transverse momentum of the produced hadron h 2 in the laboratory frame. The azimuthal angle ϕ h and ϕ S are measured relative to the lepton plane [48]. The single-spin Sivers asymmetry in SIDIS is defined as the ratio of structure functions and can be written in the TMD factorization as where M is the mass of the nucleon h 1 , and where f and D are the TMD parton distribution function (PDF) and fragmentation function (FF), J n is the Bessel function of the first kind and the summation runs over all active quarks and antiquarks q with electric charge e q . The arguments μ and ζ are the ultraviolet and the rapidity renormalization scale, correspondingly. The Q dependence of the ratio in Eq. (3) is due to the scales ζ 1;2 , which obey ζ 1 ζ 2 ¼ Q 4 [4,[49][50][51][52][53]. To respect it, we fix ζ 1 ¼ ζ 2 ¼ Q 2 , and also μ 2 ¼ Q 2 . The dependence on ðμ; ζÞ of a TMD distribution is dictated by the pair of TMD evolution equations [4,41,54], which, in turn, relate measurements made at different energies. In this work we use the ζ prescription [41] which consists in selecting the reference scale ðμ; ζÞ ¼ ðμ; ζ μ ðbÞÞ on the equipotential line of the field anomalous dimension that passes through the saddle point. In this case, the reference TMD distribution is independent on μ (by definition) and perturbatively finite in the whole range of μ and b. The solution of the evolution equations can be written [41,55] in the following simple form: and similar for other TMDs. The function f ⊥ 1T;q←h ðx; bÞ ¼ f ⊥ 1T;q←h (x; b; μ; ζ μ ðbÞ) on the right-hand side of Eq. (5) is the optimal Sivers function [55]. The function ζ μ ðbÞ is a calculable function of the universal nonperturbative Collins-Soper kernel Dðb; μÞ [56]. The N 3 LO expression used in this work is given in Ref. [42].
Drell-Yan process.-The relevant part of the differential cross section for DY reaction The variables φ and q T are the angle and the transverse momentum of the electroweak boson measured in the center-of-mass frame and y is its rapidity. The experimentally measured transverse spin asymmetry is where M is the mass of the polarized hadron h 1 , and where f 1 and f 2 are TMD PDFs for hadrons h 1 and h 2 . Often, the experiment provides measurements related to A TU [Eq. (7)]. In particular, the process h 1 ðP 1 Þ þ h 2 ðP 2 ; SÞ → l þ l − þ X (i.e., with the polarized hadron h 2 ) measured by COMPASS [58] is described by where the exchange of Sivers and unpolarized TMD PDFs takes place in the numerator of Eq. (7) and M refers to h 2 . Another important case is the asymmetry A N [59] measured by the STAR Collaboration and defined such that A N ¼ −A TU [60]. The STAR measurements are made for W AE =Z-boson production, and thus B DY n [Eq. (8)] should be updated replacing P q e 2 q by an appropriates structure, which can be found, e.g., in Ref. [42].
Nonperturbative input.-In addition to the Sivers function, SSAs (3), (7) contain nonperturbative unpolarized TMDs and the Collins-Soper kernel. We use these functions from Ref. [42] (SV19). SV19 was made by the global analysis of a large set of DY and SIDIS data, including precise measurements made by the LHC, performed with N 3 LO TMD evolution and NNLO matching to the collinear distributions. The unpolarized TMD PDFs for the pion were extracted in the same framework in Ref. [43]. In these extractions the Collins-Soper kernel is parameterized as , D resum is the resummed N 3 LO expression for the perturbative part [61], and c 0 is a free parameter. The linear behavior at large b of Eq. (9) is in agreement with the predicted nonperturbative behavior [62,63] and coefficient c 0 can be related to the gluon condensate [63].
It is customary in the TMD phenomenology to match TMDs to collinear distributions at small b [4,50,64,65]. In the present work, we do not use the matching of the Sivers to QS function [34,65,66], since it is not beneficial in the Sivers case. The reason is that QS function is not an PHYSICAL REVIEW LETTERS 126, 112002 (2021) autonomous function, but mixes with other twist-three distributions [67]. Therefore, a consistent implementation of the matching requires introduction of several unknown functions-subjects of fitting. Instead, we use the reversed procedure. We consider the optimal Sivers function as a generic nonperturbative function that is extracted directly from the data. QS function is then obtained from the smallb limit of the extracted Sivers function. For the Sivers function, we use the following ansatz: We will distinguish separate functions for u, d, s quarks, and a single sea Sivers function forū,d, ands quarks. The Sivers function does not have the probabilistic interpretation and can have nodes [68,69], which is realized by the parameter ϵ. The presence of a node is predicted by various models [68,[70][71][72]. We set β s ¼ β sea and ϵ s ¼ ϵ sea ¼ 0, since they are not restricted by the current experimental data. In total, we have 12 free parameters in our fit. Notice that the absence of the small-b matching is advantageous for our analysis as it allows us both to circumvent the difficulties of evolution of QS functions and to reach N 3 LO precision [73]. Such a strategy is allowed in the ζ prescription, and would also work in other fixed scale prescriptions [62], but is not consistent in the resummation-like schemes, e.g., used in Refs [35,38,39]. In the latter scheme, one would need to use the matching function for the Sivers function at N 3 LO, which it is not yet available [65].
Fit of the data.-The TMD factorization theorems are derived in the limit of large Q and a small relative transverse momentum δ, defined as δ ¼ jP hT j=ðzQÞ in SIDIS, δ ¼ jq T j=Q in DY process. We apply the following selection criteria [42,43] onto the experimental data: hQi > 2 GeV and δ < 0.3: The Sivers asymmetry has been measured in SIDIS and DY reaction [58,59,[74][75][76][77][78]. In total, after data selection cuts [Eq. (12)], we use 76 experimental points. We have 63 points from SIDIS measurements collected in π AE and K AE production off the polarized proton target at HERMES [74], off the deuterium target from COMPASS [76], and the 3 He target from JLab [78,79], and h AE data on the proton target from COMPASS [80]. We use 13 points from DY measurements of W AE =Z production from STAR [59] and pioninduced DY from COMPASS [58]. Let us emphasize that the recent 3D binned data [74] from HERMES allowed us to select a sufficient number of data (46 points) from SIDIS measurements. COMPASS and JLab measurements in SIDIS are done by projecting the same data onto x, z, and P hT . In order not to use the same data multiple times and for better adjustment to the TMD factorization limit, we use only P hT projections.
The evaluation of the theory prediction for a given set of model parameters is made by ARTEMIDE [81]. The estimation of uncertainties utilizes the replica method [82], which consists of the fits of data replicas generated in accordance with experimental uncertainties. From the obtained distribution of 500 replicas, we determine the values and the errors on parameters and observables, including, for the first time, propagation of the errors due to the unpolarized TMDs. We use the mean value of the resulting distributions due to SV19 uncertainty as the central fit value (CF value), which is our best estimate of the true values for the free parameters. The uncertainty is given by a 68% confidence interval (68% CI) computed by the bootstrap method. The resulting replicas are available as a part of ARTEMIDE [83].
We performed several fits with different setups. In particular, we distinguish the fits with and without the inclusion of DY data. We found that the Sivers function extracted in the SIDIS-only fit nicely describes the DY data without extra tuning. Indeed, the N 3 LO SIDIS-only fit has χ 2 =N pt ¼ 0.87 and without any adjustment describes also the DY data with χ 2 =N pt ¼ 1.23.
The combined SIDIS þ DY fit reaches a very good overall χ 2 =N pt ¼ 0.88 for all 76 DY and SIDIS data points, with χ 2 =N pt ¼ 0.88 for SIDIS and χ 2 =N pt ¼ 0.90 for DY cases.  [74], COMPASS pion-induced DY process [58], and STAR W AE =Z data [59]. Open symbols: data not used in the fit. Orange line is the CF and the blue box is 68% CI.
Parameters of the Sivers function resulting from SIDIS-only and SIDIS þ DY fits are compatible with each other [84]. The quality of data description in the SIDIS þ DY N 3 LO fit can be seen in Fig. 1.
We have performed a fit without the sign change of the Sivers function from Eq. (1) in order to estimate the significance of the sign change from the data. The resulting fit does exhibit tensions between DY and SIDIS data sets, however, the fit has χ 2 =N pt ¼ 1.0 and cannot exclude the same sign of Sivers functions in DY and SIDIS kinematics. The sign of the sea-quark Sivers function plays the central role here. Indeed, the sign of the DY cross section is mostly determined by the sea contribution due to the favored q þq → W=Z=γ subprocess, whereas the sea contribution in SIDIS is suppressed. Therefore, with the current data precision, the flip of the sign for the N sea parameter alone is sufficient to describe the data and almost compensates the effect of the overall sign-flip [Eq. (1)] at the level of the cross section. The future data from RHIC and COMPASS together with EIC and JLab will allow us to establish the confirmation of the sign change [Eq. (1)].
Extracted Sivers functions.-The extracted Sivers functions in b space for u and d quarks are shown in Fig. 2. One can see that our results confirm the signs of the u quark (negative) and d quark (positive) at middle-x range known from the previous analyses [25][26][27][28][29][30][31][32][33][34][35][36][38][39][40], and also show a node for the u quark at large x. We have not explicitly used the positivity relation [85] of Sivers functions because it is only a LO statement and can be violated in higher order calculations. However, we verified numerically that our results do not exhibit any substantial violation of positivity bounds.
The magnitude of s and sea quarks contribution in our fit is substantially different from other extractions where the biased ansatz f ⊥ 1T ðxÞ ∝ f 1 ðxÞ is used [27,[29][30][31][32][33][34][35][36]38,39] and the nonvalence contribution is artificially suppressed. In our case, the sea-and s-quark Sivers functions are comparable in size with u and d quarks, at x > 0.1 (and vanish at x < 0.1). Our unbiased extraction of the Sivers function reproduces large SSAs measured in the DY W AE =Z processes, see Fig. 1.
Determination of the Qiu-Sterman function.-The Sivers function at small b can be expressed via the operator product expansion (OPE) through the twist-three distributions [65,66,86]. At the OPE scale μ ¼ μ b ≡ 2 expð−γ E Þ=b the NLO matching expression [65] depends only on the QS function and can be inverted. We obtain the following relation: whereȳ ¼ 1 − y, α s is the strong coupling constant, and T q and G ðþÞ are QS quark and gluon functions. This expression is valid only for small (nonzero) values of b. We use  The solid black line shows the CF value and blue band shows 68% CI. The light green band shows the band obtained by adding the gluon contribution G ðþÞ . We compare our results to JAM20 [40] (gray dashed lines) and ETK20 [39] (violet dashed lines). b ≃ 0.11 GeV −1 such that μ b ¼ 10 GeV. The resulting QS functions are shown in Fig. 3. To estimate the uncertainty due to the unknown gluon contribution we allow for G ðþÞ ¼ AEðjT u j þ jT d jÞ. The resulting 68% CI uncertainty band and comparison to Refs. [39,40] are also shown in Fig. 3.
Conclusions.-In this Letter, we have presented the first extraction of the Sivers function that consistently utilizes previously extracted unpolarized proton and pion TMDs, and uses SIDIS, pion-induced Drell-Yan, and W AE =Z-bozon production experimental data. The extraction is performed at unprecedented N 3 LO perturbative precision within the ζ prescription that allows us to unambiguously relate the Sivers function and QS function. This relation has been used to obtain the QS function and to evaluate the influence of the unknown gluon QS function. We also examined the significance of the predicted sign change of Sivers functions in SIDIS and DY processes. Our results compare well in magnitude with the existing extractions [25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40] and confirm the signs of Sivers functions for u and d quarks while we also obtain non-negligible Sivers functions for antiquarks. We will study the impact of the future Electron-Ion Collider data on the knowledge of the Sivers function in a future publication.
Our results set a new benchmark and the standard of precision for studies of TMD polarized functions and are going to be important for theoretical, phenomenological, and experimental studies of the 3D nucleon structure and for the planning of experimental programs.