Global analysis of the Sivers functions at NLO+NNLL in QCD

We perform global fit to the quark Sivers function within the transverse momentum dependent (TMD) factorization formalism in QCD. We simultaneously fit Sivers asymmetry data from Semi-Inclusive Deep Inelastic Scattering (SIDIS) at COMPASS, HERMES, and JLab, from Drell-Yan lepton pair production at COMPASS, and from $W/Z$ boson at RHIC. This extraction is performed at next-to-leading order (NLO) and next-to-next-to leading logarithmic (NNLL) accuracy. We find excellent agreement between our extracted asymmetry and the experimental data for SIDIS and Drell-Yan lepton pair production, while tension arises when trying to describe the spin asymmetries of $W/Z$ bosons at RHIC. We carefully assess the situation, and we study in details the impact of the RHIC data and their implications through different ways of performing the fit. In addition, we find that the quality of the description of $W/Z$ vector boson asymmetry data could be strongly sensitive to the DGLAP evolution of Qiu-Sterman function, besides the usual TMD evolution. We present discussion on this and the implications for measurements of the transverse-spin asymmetries at the future Electron Ion Collider.


Introduction
One of the most important discoveries in hadronic physics over the past decades has been the measurements of large spin asymmetries in hadronic interactions [1,2]. These experimental measurements eventually lead to the conclusions that not only are QCD dynamics important for describing experimental data; but that these experimental measurements can be used to probe the internal structure of hadrons. For the past forty years, a major focus of the hadronic physics community has been precision extractions of the distribution functions which describe this internal structure [3][4][5][6]. In particular, the Sivers function [7,8], which provides the transverse momentum distribution of unpolarized quarks in a transversely polarized proton via a correlation between the transverse momentum of the quark and the transverse spin of the proton, has received considerable attention in recent years. By studying the Sivers function, major advancements have been made in the understanding of the spin-transverse momentum correlation and factorization theorems. For instance, theoretical investigation of the Sivers function led to the discovery that this function observes modified universality between semiinclusive deep inelastic scattering (SIDIS) and Drell-Yan process [9][10][11][12][13]. Roughly speaking, this effect occurs because the phase which is produced from the re-scattering of the unpolarized quark and the color remnant field of the initial-state hadron is opposite between these two processes. A fundamental goal of the future Electron Ion Collider (EIC) [4] will be high precision determination of these so-called transverse momentum dependent distribution functions (TMDs) over a wide range of energy scales, i.e. the so-called quantum three-dimensional (3D) imaging of the hadrons.
While the extraction of TMDs is an essential ingredient in describing transverse momentum dependent observables, high precision determination of these distributions functions has remained a challenge. The Sivers function and all other TMDs are non-perturbative objects. These TMDs must then be either computed on a lattice [5,14], or fitted from spin asymmetry data with the use of TMD factorization theorems [15][16][17][18]. The TMD factorization theorems are valid in the region where q ⊥ /Q 1 where q ⊥ is the transverse momentum resolution scale and Q is the relevant hard scale of the collision. In this region, the cross section can be factorized in terms of transverse momentum dependent parton distribution functions (TMDPDFs) and/or transverse momentum dependent fragmentation functions (TMDFFs), and perturbatively calculable short distance hard coefficients. In this paper, we rely on the TMD factorization theorems for SIDIS and Drell-Yan processes.
Despite the challenges involved with fitting TMDs, tremendous progress has been made in the field over the past few years. In particular, the focus of the field has been to increase the perturbative accuracy of the extractions of the TMDs. In [19,20] global extractions of the unpolarized TMDPDFs and TMDFFs were performed from SIDIS and Drell-Yan data at leading order (LO) and next-to-leading logarithmic (NLL) accuracy. In [21] the unpolarized TMD-PDFs were extracted at next-to-next-to leading order (NNLO) and next-to-next-to leading logarithm (NNLL) accuracy. Recently in [22] the TMDPDFs were extracted at NNLO+N 3 LL accuracy from Drell-Yan data; while in [23] the TMDPDFs and TMDFFs were extracted simultaneously from SIDIS and Drell-Yan data at NNLO+N 3 LL in which the authors further include target mass corrections as well as q ⊥ /Q power corrections. Progress has also been made in understanding the predictive power of the TMD factorization formalism in different kinematic regions [24,25], and in matching with the collinear factorization [26][27][28].
In this paper, we perform the first fit at NLO+NNLL to the Sivers function, one of the most known spin-dependent TMDs. Previously, the highest precision extraction of the Sivers asymmetry has been at LO+NLL in [29,30]. While the focus of phenomenology for unpolarized TMDs is the effects of the TMD evolution, the DGLAP evolution of twist-three function, the collinear counterpart that enters the TMD evolution for spin-dependent TMDs, introduces additional complications for fits to transverse spin-asymmetry data. For example, in the study of TMD Sivers functions with TMD evolution, the collinear twistthree Qiu-Sterman functions arise. The evolution of Qiu-Sterman function has been studied extensively in the literature [31][32][33][34][35][36][37][38][39], however a method of performing the full evolution of this function has not been well established. Nevertheless in the extractions of the Sivers functions in the literature, two approximate schemes for performing this evolution have been used in the literature. For example, in [30], the DGLAP evolution of the Qiu-Sterman function is treated to be the same as the unpolarized PDF. On the other hand, in [40], the authors use a large-x approximation for the splitting kernel [34,35] in the evolution equation of the Qiu-Stermn function. In this paper, we carefully compare the impact of these two schemes on the extraction of the Sivers function.
We perform the first global extraction of the Sivers function from all different processes, including SIDIS at HERMES, COMPASS, and JLab, Drell-Yan lepton pair at COMPASS, and W/Z production at RHIC. To perform the fit, we note that a large number of experimental data are available. At HERMES, the Sivers function has been probed by measuring both pion and kaon production in SIDIS on a proton target [41]. At COMPASS, the Sivers asymmetries have been measured in [42] for unidentified charged hadron production from the proton target, with a re-analysis of this data in [43]. The measurements with a deuteron target are presented in [44]. The Sivers function has also been probed for a neutron target at JLab for pion production in [45]. To test the modified universality prediction, Drell-Yan Sivers asymmetries have been measured at COMPASS [46] for virtual photon (or lepton pair) production at relatively small energy scales of Q ∼ a few GeV, as well as RHIC [47] for W and Z production at much large energy scales, Q ∼ M W/Z . The rest of the paper is organized as follows. In Sec. 2, we summarize the relevant TMD factorization formalism for SIDIS and Drell-Yan processes. In Sec. 3, we first discuss our nonperturbative parameterizations for the unpolarized TMDPDFs and TMDFFs, and benchmark them with the SIDIS hadron multiplicity and Drell-Yan cross section data. We then present our non-perturbative parametrization for the Sivers function, and discuss how we perform the DGLAP evolution of the Qiu-Sterman function. In Sec. 4, we present our fit results, where we explore several different ways for performing the fit. In Sec. 4.1 we present the results of a simultaneous fit to the low energy data from SIDIS and the COMPASS Drell-Yan data. In Sec. 4.2 we study the impact of the high energy data from RHIC. In Sec. 4.3 we study the impact of the DGLAP evolution scheme for the Qiu-Sterm function on the fit. In Sec. 4.4 we present the global fit where we include Sivers asymmetry data from all processes. In Sec. 5 we give predictions for Sivers asymmetry at the EIC. We conclude our paper in Sec. 6.

Formalism
In this section, we provide the TMD factorization formalism for the Sivers asymmetry. We begin in Sec. 2.1 with the SIDIS formalism, while in Sec. 2.2 and 2.3 we present the formalism for Drell-Yan lepton pair and W/Z boson production, respectively.

Sivers Formalism in SIDIS
The differential cross section for SIDIS, e( ) + p (P, S ⊥ ) → e ( ) + h (P h ) + X, where S ⊥ is the transverse spin vector of the polarized nucleon, can be written as the following form [16,48] where the phase space dPS = dx B dQ 2 dz h d 2 P h⊥ , the electron-proton center-of-mass (CM) energy S = (P + ) 2 and the exchanged virtual photon momentum q = − with Q 2 = −q 2 , and the usual SIDIS kinematic variables are defined as As shown in Fig. 1, the plane which contains the initial and final lepton momentum vectors is the lepton plane, while the momentum vectors of the photon and final state hadron form the hadron plane. The azimuthal angle of the hadron plane with respect to the lepton plane is denoted φ h , while the azimuthal angle of the transversely polarized proton spin with respect to the lepton plane is denoted φ s . We follow the Trento conventions [49] for the definition of the azimuthal angles. In this expression, σ DIS 0 is the leading order (LO) electromagnetic scattering cross section given by where α EM is the electromagnetic fine structure constant.
in Eq. (2.1) are the unpolarized and transversely polarized structure functions, respectively. The experimentally measured quantity, the Sivers asymmetry, , for this process is given in terms of the structure functions as follows The momentum space expression for these structure functions are given by where the hard factor, H DIS (Q; µ), is given in [50,51] as follows In these expressions, we have used the short-hand notation for the convolution integrals. In these expressions e q is the fractional electric charge for the quarks. k ⊥ represents the transverse momentum of the quark relative to the nucleon, while p ⊥ is the transverse momentum of the final state hadron relative to the fragmenting quark.ĥ = P h⊥ /P h⊥ is the unit vector which points in the direction of the final-state hadron transverse momentum and M is the mass of the struck nucleon.
is the unpolarized TMDFF. In these expressions µ and ζ are the renormalization and rapidity (Collins-Soper) scales [15], which are used to regulate ultraviolet and rapidity divergences, respectively. Moreover, the rapidity scales obey the relation ζ A ζ B = Q 4 in the TMD region.
The expressions for the structure functions are simplified by going to the b-space, the Fourier conjugate space to the transverse momentum space. In the b-space, these expressions become Here the b-space TMDs are defined as At small b where 1/b Λ QCD , one can perform an operator product expansion (OPE) of these functions in terms of their collinear counterparts: (2.14) are the collinear PDF, FF and the Qiu-Sterman function, respectively. The operator ⊗ denotes the convolution over the parton momentum fractions and are given by for f i/p and likewise for D h/i . In these expressions, the sum over the index i = q, g is implicit.
The convolution in the case of the Sivers function is more complicated, since it involves two kinematic variablesx 1 andx 2 : The C functions in the above equations are the Wilson coefficient functions, and their expressions at NLO are given in Appendix. A. Several comments are in order for the case of the Sivers function. First, although the coefficient function for general scales µ and ζ are quite complicated, it becomes much simpler when one chooses the canonical scales µ = √ ζ = µ b = c 0 /b, with c 0 = 2e −γ E and γ E the Euler constant. such scales are referred to as the natural scale of the TMDs. Second, there are different conventions/normalization for the Qiu-Sterman function. In our case, we first follow the Trento convention [49] for the quark Sivers function and then the convention for the Qiu-Sterman function is such that the coefficientC function at leading order in Eq. (2.16) is a simple delta function. Our convention is related to the so-called first transverse moment of the Sivers function [11,28] Third, in principle the convolution in Eq. (2.16) receives contribution not only from the Qiu-Sterman function which is a quark-gluon-quark twist-3 correlator, but also the so-called twist-3 three-gluon correlator. Since the three-gluon correlator is not well-known at the moment in phenomenology, we neglect all contributions from gluon to quark splitting in the Sivers function [39,52]. Finally, Eq. (2.13) is only defined for the Sivers function in SIDIS. Thus if one changes to the Sivers function in Drell-Yan, one should include an additional minus sign in the last line of this expression. The large logarithms present in Wilson coefficient functions are resummed in the renormalization group evolution of TMDs from the natural scale µ 2 Such a TMD evolution is encoded in the exponential factor, exp [−S], with the so-called Sudakov form factor S. The perturbative part of the Sudakov form factor is given by where Γ cusp and γ V are the cusp and non-cusp anomalous dimensions, respectively, and D is the rapidity anomalous dimension (Collins-Soper kernel) [15,53]. In this paper, we perform the resummation of these logarithms up to NNLL. All information on the anomalous dimensions up to NNLL are given in Appendices B and C. When b becomes large and thus µ b Λ QCD , the TMD evolution runs into the nonperturbative region. We follow the usual b * -prescription [54] that introduces a cut-off value b max and allows for a smooth transition from perturbative to non-perturbative region, With the introduction of b * in the Sudakov form factor, the total Sudakov form factor can be written as the sum of perturbatively calculable part and nonperturbative contribution. The final expressions for the structure functions are given by where we have replaced µ b by µ b * = c 0 /b * , and Q 0 is the reference scale of the TMDs. The functions S f NP , S D NP , and S s NP are the corresponding non-perturbative Sudakov form factors for the unpolarized TMDPDF, TMDFF, and the Sivers function, respectively, and they will be given in the next section. Note that in these expressions we have introduced the vector q ⊥ = −P h⊥ /z h , while q ⊥ = |q ⊥ | denotes its magnitude.

Sivers Formalism in Drell-Yan
For Drell-Yan scattering, p(P A , S ⊥ ) + p(P B ) → [γ * (q) →] + − + X, the differential cross section with the relevant terms is given in [55][56][57][58] by the expression where dPS = dQ 2 dy d 2 q ⊥ , y is the rapidity of the lepton pair while q ⊥ and Q are the transverse momentum and invariant mass of the virtual photon, respectively. Here, W U U and W sin(φq−φs) U T are the unpolarized and transversely polarized structure functions. Note that we have deviated from the notation in [56] by writing the Drell-Yan structure functions as W in order to differentiate them from the SIDIS structure function. The leading order electro-magnetic scattering cross section is given by where S = (P A + P B ) 2 is the center of mass energy squared and N C = 3 is the number of color. As shown in Fig. 2, the plane which is perpendicular to the spin vector S ⊥ and which also contains the initial hadrons forms the hadron plane. The plane which contains the hadron momenta and which contains the vector boson (i.e. γ * here) momentum generates the vector boson plane. We use the convention that the polarized hadron moves in the z direction while S ⊥ moves in the y direction. We note that the convention for the x and z axes must be reversed in order to compare with the COMPASS Drell-Yan data. For the Drell-Yan production, φ q , the azimuthal angle of the vector boson, and φ s , the azimuthal angle of S ⊥ generate the sin(φ q − φ s ) modulation for this process.
Analogous to the asymmetry in SIDIS, the Drell-Yan Sivers asymmetry can be written in terms of the structure function as 1 A sin(φq−φs) In the TMD formalism, these structure functions are given by the following expressions For Drell-Yan process, the above convolution in the structure functions is given by where x a and x b are the momentum fractions of the hadrons carried by the quarks and are given by The usual Feynman-x is related to x a,b as follows x F = x a − x b , which will be used in the next section. On the other hand, k a⊥ and k b⊥ are the transverse momenta of the parton relative to their corresponding nucleon. The hard function is given in [18] by The expressions for the structure functions can once again be simplified by going to the b-space. At this point, it might be important to emphasize again that the Sivers function f ⊥ 1T above for the Drell-Yan process differs by a sign from that in SIDIS in Eq. (2.6): This will lead to slightly different definition for the Sivers function in the b-space: Note that another single spin asymmetry denoted as AN for Drell-Yan process has also been frequently used in the literature, which is related to the Sivers asymmetry defined here by a minus sign: AN = −A sin(φq−φs) U T . For details, see [56].
Note the additional minus sign in the second line of the equation, in comparison with the corresponding SIDIS expression in Eq. (2.13). The final expressions for the b-space structure functions are given by Note that in the second expression, we have already taken into account the sign change in the Sivers functions between DY and SIDIS processes in Eq. (2.32).

Sivers formalism for W/Z Production
The case for W/Z boson production in the proton-proton collisions is similar to the case for virtual photon production. In this case, the hard scale Q is set equal to the mass of the produced vector boson, Q = M W, Z . The expression for the differential cross section is given by where the phase space dPS = dy d 2 q ⊥ and V = W, Z. The leading-order scattering cross sections are given by where G F is the Fermi weak coupling constant. On the other hand, the structure functions are given by where we have Here |V qq | 2 is the CKM matrix, while V q and A q are the vector and axial couplings of the Z boson to a quark of flavor q. Just like Eq. (2.26) in the last section, the asymmetry can be written as a ratio of these structure functions in the exactly same form.

Non-Perturbative Parameterization
Now that we have included all of the perturbative elements of the Sivers asymmetry, we begin discussing the non-perturbative contributions to the Sivers function. As we have seen in the previous section, the Sivers asymmetry depends not only on the Sivers functions but the unpolarized TMDs as well. Therefore, in order to isolate the fit to affect only the Sivers function from these experimental data, it is first necessary to fix the non-perturbative evolution of the unpolarized TMDs. In Sec. 3.1, we choose a parameterization for the unpolarized TMDPDF and TMDFF from a previous extraction and use this formalism to describe unpolarized SIDIS and Drell-Yan data. In Sec. 3.2 we provide the details of our numerical scheme for the Sivers function.  Figure 3. The experimental data for Drell-Yan lepton pair production measured by the E288 collaboration [59] plotted as a function of q ⊥ /Q are compared with the normalized theoretical curve. Different colors represent different invariant mass of the lepton pair from 4 < Q < 5, 5 < Q < 6, 6 < Q < 7, 7 < Q < 8, 8 < Q < 9, 11 < Q < 12, 12 < Q < 13, 13 < Q < 14 GeV, respectively. Three panels correspond to different energies for incident proton beams: 200 GeV (left), 300 GeV (middle), and 400 GeV (right).  The non-perturbative evolution functions for the unpolarized TMDs have been extracted widely in the literature. Because we perform a simultaneous fit between SIDIS and Drell-Yan data in this paper, the appropriate parameterizations for the unpolarized TMDs are those that have also been obtained in simultaneous fits. Furthermore since we perform our fit at NLO+NNLL, the optimal parameterization is one that has been obtained at the same perturbative order.

Numerical Scheme for Unpolarized TMDs
Simultaneous extractions from SIDIS and Drell-Yan data have been performed in [19,20,23,62] 2 . In [62] the extraction was performed at NLO+NLL. Similarly in [19,20] the extraction was performed at LO+NLL. In [23], the authors performed the fit of the unpolarized data at NNLO+N 3 LL level, where they further included both m/Q and q ⊥ /Q power corrections. This could introduce additional complications when performing the fit to the Sivers asymmetry, since those power corrections are likely to be different for spindependent cross sections.
In view of the current status, we choose the non-perturbative parametrization in [62] for the unpolarized TMDs in our study at NLO+NNLL accuracy. We will first verify that such a parametrization describes the unpolarized experimental data well. From [62], the Figure 5. The COMPASS multiplicity in [61] for charged hadron production from a deuteron target is compared with the normalized theory curve. The triangular points represent the h + data points while the circular data points represent the h − data points. For better presentation, the h + data is offset by a factor of 0.4.
non-perturbative factors in Eqs. (2.14) and (2.15) have the following form The factors which contain g f 1 and g D 1 contain information on the Gaussian width of the TMDs in momentum space at the initial scale Q 0 , while the factor which involves g 2 controls how the TMDs evolve from Q 0 to Q. The latter is universal to all TMDs [15] and will enter into our discussion in the Sivers non-perturbative parameterization. The values of the parameters that were obtained in this reference are given by Note that in the expression of Eq. (3.1), the non-perturbative parameterization is independent of x. Thus we have dropped explicit dependence on the variable x. At this point, it is important to note that for the COMPASS Drell-Yan data in [61], the asymmetry was measured for π + p scattering. In [63] the pion TMDPDF was extracted from the experimental data in [64] and it was found that g f 1 = 0.082 for pions.
To perform numerical calculations, we choose to use HERA NLO as 118 parametrization in [65] for the collinear parton distribution functions. For the collinear pion fragmentation function, D π/q (z h , µ b * ), we use the DSS14 parameterization [66]. While for the collinear kaon fragmentation function D K/q (z h , µ b * ), we use the DSS17 parameterization in [67]. For unidentified charged hadrons, we follow the work in [23] to use the approximation To demonstrate that this parameterization describes the unpolarized TMDs, we now compare this numerical scheme with the unpolarized TMD data. We start this comparison by examining a sample of Drell-Yan data in order to check the validity of the scheme for the TMDPDF. We note that the Drell-Yan Sivers asymmetry data which enters into our fit from COMPASS and RHIC do not contain so-called fiducial cuts. In order to avoid complications associated with these cuts on Drell-Yan data, we choose to benchmark our expression for the unpolarized cross section against the E288 data [59], which also does not contain fiducial cuts, see Tab. 2 of [22]. For E288, the target nucleus is Copper. In order to describe the Copper TMDPDF, we use nuclear modification prescription in [68]. In Fig. 3, we plot the theoretical curve against the experimental data [59], as a function of q ⊥ /Q. For each bin, we have normalized the theory such that the theory and data are equal at the first point. Different colors represent different invariant mass of the lepton pair from 4 < Q < 5, 5 < Q < 6, 6 < Q < 7, 7 < Q < 8, 8 < Q < 9, 11 < Q < 12, 12 < Q < 13, 13 < Q < 14 GeV, respectively. Three panels correspond to different energies for incident proton beams: 200 GeV (left), 300 GeV (middle), and 400 GeV (right). We find that the parameterization of [62] is well-suited at describing the shape of the Drell-Yan data.
To check the validity of our scheme for the unpolarized TMDFFs, we now examine the HERMES multiplicity defined as where the superscript h denotes the species of the final state observed hadron, and the subscript "H" represents the HERMES data. We also study the COMPASS multiplicity data, which has a slightly different convention and is given by where the subscript "C" denotes the COMPASS data and M h H is defined in Eq. (3.4). On the other hand, the denominator in Eq. (3.4) is the inclusive DIS cross section and is given by where F 2 is the usual DIS structure function while F L is the longitudinal structure function. For their precise definitions see [69]. We compute the denominator at the NLO by using the APFEL library [70].
In the left panel of Fig. 4 we plot the HERMES pion multiplicity data [60] as a function of q ⊥ /Q along with the numerical results for the theory. In the right panel of this figure we plot kaon multiplicity data and theory. As shown in the figure, different colors represent different average z h values from z h = 0.15, 0.23, 0.28, 0.34, 0.34, 0.42, 0.53, respectively. In these plots, we have normalized the theory so that data is equal to the theory at the second point of each data set 3 . In Fig. 5, we plot the COMPASS multiplicity data [61] for charged hadron production from a deuteron target along with the numerical results of our scheme. The triangular points represent the h + data points while the circular data points represent the h − data points. Here again, different colors represent different z h = 0.2, 0.3, 0.4, 0.6, respectively. From these plots, we find that the presented parameterization work very well at describing the shape of the multiplicity data for both HERMES and COMPASS data, indicating that the scheme for the TMDFFs are valid.

Numerical Scheme for Sivers Function
Now that the non-perturbative evolution for the unpolarized TMDs have been fixed, we present the numerical scheme for the Sivers function in our fit. Analogous to the unpolarized TMDPDF, we take the polarized non-perturbative parameterization As we have emphasized in the previous section, the parameter g 2 is spin-independent and thus we take the same value as in the unpolarized TMDs in Eq. (3.3). On the other hand, we introduce the parameter g T 1 , which describes the Gaussian width of the momentum space distribution for the Sivers function and will be a fit parameter. We once again note that since this parameterization is independent of x, we will drop its explicit dependence in future notation.
For the Qiu-Sterman function T F q/p , we find that the parameterization in [29] is still the most economical choice, which sets T F q/p (x, x, µ 0 ) to be proportional to the unpolarized PDF f q/p (x, µ 0 ) at some initial scale µ 0 : Note that N q (x) characterizes the non-perturbative collinear physics of the Qiu-Sterman function and is to be fit from the experimental data. In this expression, the parameters α u and N u are used to fit the up quarks. α d and N d are the fit parameters for the down quarks and Nū, Nd, N s , Ns, α sea are for sea quarks and β q = β is the same for all flavors. This parameterization enforces that the form of the sea quarks is the same while the normalization of each sea quark can vary. Overall we use 11 parameters in total to perform the fit, including g T 1 .
In order to obtain a numerical result for the Sivers function in Eq. (2.16), DGLAP evolution of the Qiu-Sterman function must be performed from µ 0 to the natural scale, µ b * . As we have emphasized, the DGLAP evolution of the Sivers function has been studied extensively in the literature, see for instance [31][32][33][34][35][36][37][38][39]. However, to perform the full evolution of the Qiu-Sterman function is highly nontrivial due to its dependence on two momentum fractions x 1 , x 2 in general [31,71]. Thus in the TMD global analysis, the evolution of the Qiu-Sterman function has been implemented under certain approximations. There are two schemes that are used to perform this evolution in the literature. For both schemes that we discuss in this paper, the relevant DGLAP evolution equation for the Qiu-Sterman function is given by the expression In the first scheme that we consider, from [35], the authors show that at large x, the transverse spin dynamics leads to a modification to the quark to quark splitting kernel, P T q←q , with where P q←q (x) is the standard quark to quark splitting kernel for unpolarized PDFs, This scheme has been used for instance in [40]. In the second scheme, for phenomenological purposes, the evolution of the Qiu-Sterman function has often been treated to be the same as the unpolarized collinear PDF, with P T q←q (x) = P q←q (x). See e.g. Ref. [30]. Apparently, for both cases, we can write the relevant spitting kernel as where η is a parameter that controls the numerical scheme used to perform the DGLAP evolution. When η = N C , the evolution matches the result of [35]. On the other hand, for the second scheme that we consider, we set η = 0 so that the evolution models the standard DGLAP evolution of the unpolarized PDF.
To solve this evolution equation, it is useful to take the Mellin transform of this expression; for details on Mellin-space evolution, see Sec. 3 in [72]. After performing the Mellin transform of this expression, the evolution equation becomes (3.14) In this expression, T F q/p (N, µ) is the Mellin transforms of the Qiu-Sterman function, i.e.
Similarly γ(N ) is the Mellin transform of P T q←q (x) which can be written as Here γ u (N ) is the Mellin transform of the unpolarized splitting function P q←q (x) and is given by with S 1 (N ) the harmonic sum function.
In the region where µ b * < m b , the mass of the b quark, the solution of the evolution equation is given by .

(3.18)
Here β 0 (µ 0 ) = 11 − 2/3 n f (µ 0 ), where n f (µ 0 ) is the number of active flavors at the scale µ 0 . In the region where µ b * > m b , the solution of the evolution equation is given by where T F q/p (N, m b ) is given by 20) and n f (µ b * ) is the number of active flavors at the scale µ b * . In order to construct the Sivers function in Eq. (2.16) at NLO, there is an additional convolution of the coefficient C function and the Qiu-Sterman function. We find that it is useful to first take its Mellin transform and thus the convolution over the momentum fraction becomes a simple product in Mellin space: where the parameter c must be taken such that all of the singularities in the function f ⊥q 1T,q/p c + ze iφ , b; µ, ζ lie to the left of the line x = c in the imaginary plane. In our code, we use c = 2 which satisfies this criteria. We also take φ = π/4 to optimize the numerical integration. fit scheme SIDIS Drell-Yan W/Z N data η in evolution  In this section we present a simultaneous fit to measurements of the Sivers asymmetry from SIDIS data sets from JLAB in [45], HERMES in [41], COMPASS in [43,44] and the COMPASS Drell-Yan data in [46]. We note that we do not include the COMPASS data set in [42] since the data set in [43] is a re-binning of this set. Furthermore the data set in [43] was projected into two sets of data z h > 0.1 and z h > 0.2. To avoid fitting correlated data sets, we choose to fit only the z h > 0.1 data set. We then compare our prediction for the RHIC asymmetry against the RHIC data.

Fit Results
While typical kinematic cuts from unpolarized SIDIS fits for instance in [23] select only data which has q ⊥ /Q < 0.25, we find that this selection process leaves very few data points for the available Sivers data. In Fig. 6 we plot a histogram of the selected data SIDIS data as a function of q ⊥ and Q. We find that the cut q ⊥ /Q < 0.25 leaves only 12 SIDIS data points, while the cut q ⊥ /Q < 0.5 leaves 97 data points. In fact, we find that the majority of the data has q ⊥ /Q > 0.5. In order to retain a large enough data set to perform a meaningful fit we perform the cut q ⊥ /Q < 0.75. Furthermore to restrict the selected data set to the TMD region, we also enforce that the SIDIS data must have P h⊥ < 1 GeV. At the same time in order to avoid the threshold resummation region, we also enforce that z h < 0.7.
In order to perform the fit, we use the MINUIT package [73,74] to minimize the χ 2 . In this section, we define the χ 2 as where E i are the central values of the experimental measurements, ∆E i are the total experimental errors, T i ({a}) is the theoretical value at the experimental kinematics, and {a} is a vector containing the fit parameters. For this section, we take η = N C to perform the DGLAP evolution of the Qiu-Sterman function, referred to as fit 1 in Tab. 1. In order to optimize the minimization process, the denominator of our asymmetry is pre-calculated at the beginning of the fit. We also perform pre-calculations for the unpolarized TMDs and use grid interpolation in the numerator of the asymmetry. For the NLO Sivers function, we find that the Mellin space prescription leads to a massive speeds compared to performing the convolution integrals. Furthermore we use the numerical method in [75] to perform all Bessel integrals.
In order to generate an uncertainty band, we follow the work in Ref. [20,76] to use the replica method. To generate one replica, we shift each of the the data points by a Gaussian noise with standard deviation corresponding to the experimental error. The fit is performed on the noisy data 200 times as well as the no noise data. This result in 201 sets of stored fit parameters. Using each of the 201 sets of stored parameters, we calculate the asymmetry for each of the included data as well as calculate the first transverse moment of the Sivers function in Eq. (2.19) for each of the quark flavors. The uncertainty band is generated at each point by retaining all contribution within the 68% region.
In Table. Table 3. The distribution of experimental after taking the kinematic cuts q ⊥ /Q < 0.75, P h⊥ < 1 GeV, and z < 0.7. The column Q avg gives the average hard scale for the measured data set. On the right column, we have included the χ 2 /N data for each set of data from the extraction in fit 1. The RHIC data was not included into the fit. Here we give the χ 2 /N data for the prediction.
function at the initial PDF scale, f ⊥ (1) 1T (x, µ 0 ) with µ 0 = √ 1.9 GeV as defined in Eq. (2.19). In this figure, we have plotted all 200 replicas for each of the extracted quark flavors. We again use the middle 68% of the data points in the plot to generate the grey uncertainty band for each of the Sivers moments. For theū, s ands-quarks, the Sivers moment have been multiplied by a factor of 5 while ford, we have multiplied by a factor of −5. We find that the Sivers d function is the largest in magnitude and is positive; while the Sivers u function is nearly as large but is negative. Furthermore we find that theū andd-quark functions are nearly equal to one another in magnitude, both are more than 5 times smaller in magnitude than the valence quarks, and are both positive. For the s-quark, we find that the magnitude is approximately 5 times smaller than the valence quarks in magnitude and is negative. Finally for thes-quark, we find that the magnitude is very small and that the sign is not well determined in this fit.
In Figs. 8, 9, and 10, we plot our theoretical curves against the SIDIS data. Fig. 8 is  Figure 8. Left: The COMPASS deuteron target measurement [44] for π + , π − , K + , K − , and K 0 from top to bottom, and as a function of x B (left), z h (middle), and P h⊥ (right). Right: HERMES proton target measurement [41] π + , π 0 , π − , K + , K − , and (π + − π − ) from top to bottom, and as a function of x B (left), z h (middle), and P h⊥ (right). The data is plotted in red along with the total experimental error. The central curve in blue as well as the uncertainty band in gray are generated using the result from fit 1 in Tab. 1. for COMPASS deuteron target (left panel) and for HERMES proton target (right panel), and for both pions and kaons. Fig. 9 is for charged hadrons from COMPASS proton target. Fig. 10 is for pion production on a neutron target from JLab. Finally in Fig. 11 we plot theoretical curves against the COMPASS Drell-Yan lepton pair data in π − + p collisions. We plot the asymmetry A sin(φq−φs) U T as a function of transverse momentum q ⊥ , invariant mass Q, Feynman x F = x π −x N , momentum fraction x N in the proton target, and momentum fraction x π in the pion target, respectively. The experimental data along with the total experimental uncertainties are plotted in red. The blue curves are the theory curves from the fit with no noise. The uncertainty band in grey is generated from the stored values of the asymmetry for each of the replicas. For each data point, the maximum and minimum value of the asymmetry within the middles 68% are used to generate these error bars. As it is indicated already in Tab. 3 and as it is evident from the figures, the agreement between our theory and SIDIS and Drell-Yan data is very good, although to a less degree with the Drell-Yan data because of the much larger experimental uncertainty. We note that very recently in [77], the HERMES collaboration provided additional experimental data for the Sivers asymmetry. Since the  Figure 9. Left: The COMPASS proton target measurement for h − for 1 GeV 2 < Q 2 < 4 GeV 2 , 4 GeV 2 < Q 2 < 6.25 GeV 2 , 6.25 GeV 2 < Q 2 < 16 GeV 2 , 16 GeV 2 < Q 2 < 81 GeV 2 from top to bottom [43]. Right: Same as the left except for h + production. The central curve as well as the uncertainty band are generated using the result from fit 1 in Tab HERMES paper has not yet been published, we cannot implement this data into our fit. However, we find that there is very strong agreement between our extracted asymmetry and this new data.
In Fig. 12, we plot the prediction for the RHIC data in p + p collisions at √ S = 500 GeV using the extracted Sivers function from this fit. In the left panel, we plot the Sivers asymmetry A N as a function of rapidity for W − (left), W + (middle), and Z 0 (right), respectively. We integrate vector boson transverse momentum over 0.5 < q ⊥ < 10 GeV. On the right  Figure 12. Prediction for the Sivers asymmetry for p + p → W/Z at √ S = 500 GeV [47] using the result of fit 1 in Tab. 1. We plot only the central curve from fit 1 here since the size of the uncertainty band is small for this prediction. Left: The y dependent data integrated in q ⊥ from 0.5 to 10 GeV. Right: The q ⊥ dependent data integrated in y from −1 to 1. panel, we plot A N as a function of q ⊥ while we integrate over the rapidity |y| < 1. We find that the asymmetry for W/Z for the central fit is at most 2%, which is more than an order of magnitude smaller than the central values recorded at RHIC. This leads to a χ 2 /N data of 2.015 for the prediction for RHIC, as shown in Tab. 3. Even if one considers the very large error bars in the RHIC data, this comparison seems to indicate some tension between our theory and the RHIC data.

Impact of the RHIC data
In this section, we study the impact of the RHIC data to the fit. One possible issue which may be arising in the description of the RHIC data is that while there are a large number of experimental data at small Q, there are much less data at RHIC energies. In order to access the impact of the RHIC data, it is therefore convenient to follow the work in [78] to introduce a weighting factor to the calculation of the χ 2 . Thus in this section, the expression for the χ 2 is given by We also define the N data for this weighted fit as For the first term of Eq. (4.2), the sum is performed over all data in the previous section, i.e., all the SIDIS data plus COMPASS Drell-Yan data. In the second term, the sum is performed only over the RHIC data. In this second expression, ω is the weighting factor. In order to emphasize the contributions of the RHIC data, we choose ω = N/N R = 226/17 so that the RHIC data and the rest of the experimental data sets are equally weighed in the calculation of the χ 2 . Furthermore, in order to perform the DGLAP evolution of the Qiu-Sterman function, we take η = N C . Using this definition of the χ 2 , we perform a fit to the selected data. In Tab. 4, we provide the distribution of the χ 2 for this fit. With the addition of the weighting factor, we find that the χ 2 /N data = 1.888 for the RHIC is quite large while for the low energy data the χ 2 /N data = 0.996. This result indicates that the issue with describing the RHIC data is not that the high energy data has a small number of data points. Rather, it indicates that when using our theoretical assumptions, these sets of data disagree on the properties of the Sivers function.
In order to access which one of our theoretical assumptions is responsible for the large χ 2 of the RHIC data, we have performed several tests. Firstly, we have checked whether the quality of the description of the RHIC data was due to the cut on q ⊥ /Q. In order to check if quality of the fit is due to the value of this cut, we have performed an additional fit with the cut q ⊥ /Q < 0.5. We find that this change leads to a χ 2 /N data is 1.885 for the RHIC data. While it would be preferable to perform an fit with q ⊥ /Q < 0.25, we note that there is not enough data in this region to constrain the parameters of the fit. Because there is no strong improvement in the description of the RHIC data after applying the q ⊥ /Q < 0.5, we conclude that this cut is not responsible for the disagreement between the data sets.
Another possible assumption that could be causing the large χ 2 of the RHIC data is the assumption that the sea quarks have the same α and β parameter. To check this, we have performed a 13 parameter fit with the chosen parameter with the parameters α u , N u , β val , α d , N d , Nū, Nd, N s , Ns, α + , α − , and β sea . Here αd = αs = α + and α s = αū = α − . The introduction of the α + and α − parameterization decouples the positive and negative sea quarks from one another while the introduction of the parameters β val and β sea decouples the valance and sea quarks. However, we find that the addition of these parameters lead to a χ 2 /N data is 1.885. This implies that this assumption on the function form is not the issue.
In order to address the disagreement between the RHIC data and the rest of the data sets, in Fig. 13 we plot the profiles of the χ 2 /N data using the 13   we set all but one of the parameters equal to the values which are determined by the fit and we vary the remaining parameter about its best value. The best value determined by the fit is given by a vertical gray line. In this plot, we see that the curves for the RHIC χ 2 do not change much as the α, β, and g T 1 parameters are varied. This indicates that the RHIC data is insensitive to these parameters. On the other hand, we see that when N q parameters are varied that there are large modifications to the RHIC χ 2 . Thus, the RHIC data is sensitive to these parameters. We see from the N q plots that the RHIC data and the rest of the data sets agree on the sign of the quark-Sivers functions for N d ,Nū,Nd, and Ns while the data sets disagree strongly about the magnitude of the parameters. For N s , we see that the RHIC data appears to be insensitive to the sign of this parameter so that the disagreement is not striking. However, we find that the SIDIS and COMPASS Drell-Yan data sets indicate that the sign of the u-quark is positive while the RHIC data is indicating that the sign of the u-quark is negative. This disagreement is occurring because the fit program is attempting to describe the large positive A N asymmetry for the W + RHIC data. Thus in order to describe this data, either the Nd or Ns parameters must be large or the sign of the N u is incorrect. Since the value of the parameter N u is extremely well constrained by the SIDIS and COMPASS Drell-Yan data while the value of Nd and Ns parameters are weakly constrained, we conclude that this sign disagreement will be resolved once the magnitude of Nd and Ns parameters are addressed. Overall in Fig. 13, we see the trend that the RHIC data N q requires much larger values for the N q parameter than the SIDIS and COMPASS Drell-Yan data. Since the SIDIS and COMPASS Drell-Yan data were gathered at much lower energy scales that the RHIC data, this tension between the sets indicates that the size of the Sivers asymmetry grows as a function of the hard energy scale. This result indicates that the issue in describing the RHIC data appears because of a possible evolution effect. Since the perturbative TMD evolution of the Sivers asymmetry is known, this issue is either occurring due to the chosen non-perturbative parameterization of the Sivers function or from the choice of the DGLAP evolution of the Qiu-Sterman function. RHIC is expected to release the new measurement for W/Z Sivers asymmetry [79] in the near future in which they have much more statistics and thus smaller experimental uncertainty. The new data will be very valuable in constraining the non-perturbative component of the TMD evolution for the Sivers function. In the next section, we will study the effects of the DGLAP evolution of the Qiu-Sterman function and how they will affect the size of the asymmetry.

Effects of the DGLAP evolution
In order to examine how the DGLAP evolution of the Qiu-Sterman function affects the size of the asymmetry, we begin by examining Eq. (2.41). The largest contributions to this expression should appear in the region where µ b * ∼ Q = M V [54,80]. In this region, the size of the asymmetry is roughly proportional to T F q/p (N, M V ) in Eq. (3.19). To examine how the magnitude of the Qiu-Sterman function evolves in energy, we start from the evolution equation in the moment space in Eqs. (3.14) and (3.16), and examine the ratio of this function at the two relevant scales µ 0 and M V . One can easily show that this ratio is given by , (4.4) where N (µ 0 , M V ) is given by .

(4.5)
From this expression, it becomes clear that when η > 0, the factor N in Eq. (4.5) becomes small at large scales. Thus this factor leads to a suppression of the quark-Sivers function at large scales. Thus with this suppression factor, the values of the N q parameters must be very large in order to describe the RHIC data. On the other hand, when η < 0, the factor N leads to an enhancement of the asymmetry at large scales. In order to test the sensitivity of each data set to changes in the evolution kernel due to the change in the η parameter in Eq. (3.13), we define the quantity which gives the average percent difference between the two theory calculated with η = N C and η = 0 for a given set. In this expression, {a} are the parameters obtained from the η = N C fit. In Tab. 4, we provide the value for ∆T for each data set. We find that the result of the low energy data can vary only within a few percent on the choice of the DGLAP evolution kernel. On the other hand, the high energy RHIC data sets varies by a factor of 50% when using these different kernels.
In order to explicitly demonstrate the dependence on the DGLAP evolution scheme, in Fig. 14 we plot a profile of the χ 2 /N data as a function of the parameter η, while the rest of the parameters are fixed as those from scheme fit 2a. As we can see from this plot, χ 2 /N data for RHIC data decreases as η decreases. This indicates the RHIC seems to prefer smaller η ∼ 0 or even negative values. This trend is opposite to what is seen in the SIDIS+DY data. Because of the driving from the RHIC data, the global χ 2 seems to favor the evolution scheme with η = 0 or even negative.

Global fit of the Sivers function
In Sec. 4.1, we have presented fit 1, which was performed to Sivers asymmetry for SIDIS+DY data at the low energy. The strengths of this extraction are that the theoretical uncertainties were small so that this extraction should describe very well future low energy experiments. However, as we showed in the prediction for the RHIC data, this extraction failed to describe the high energy data. In this section, we present a fit which emphasizes the contributions  Table 5. The distribution of χ 2 for each data set for the fit 2b.
of the RHIC data in order to allow future predictions for high energy measurements of the Sivers asymmetry.
To emphasize the contributions of the high energy data, we retain the weighted definition of the χ 2 in Eq. (4.2). On the other hand, as we have seen in our model, the description of the high energy data from RHIC depends strongly on the choice of the parameter η. By performing a global fit with η = N C , we found that the χ 2 /N data for RHIC was 1.888. In order to eliminate the suppression from the −N C δ(1 − x) term in the evolution kernel Eq. (3.13). In this section, we perform the fit with η = 0. This fit is referred to as fit 2b in Tab. 1.
For this fit, we recover a χ 2 /d.o.f of 1.482 with a χ 2 /N data of 1.778 for the RHIC data. The parameter values for this fit are given in Tab. 6 while the distribution of the χ 2 is given in Tab. 5. We can see from Tab. 6 that while the extraction of the Sivers function from the low energy data could not resolve the sign of the s-quark Sivers function, this fit finds that the s-quark should be positive. At the same time, the sign of all other quark functions are consistent with the previous extraction. However, we note that the central values for the N q parameters are much larger than the previous fit. This is occurring because of the large RHIC asymmetry along with the weighting used in the fit. We see also in this table that the uncertainties in the parameters are very large and tend to skewed in one direction. The magnitude of this uncertainty is due to the large experimental uncertainties in the RHIC data while the skew favors fits which increase the size of the asymmetry.
In Fig. 15, we plot the extracted transverse momentum moment of the Sivers function, f ⊥(1) 1T (x, µ 0 ) as a function of x at the scale µ 0 = √ 1.9 GeV. The blue curve is the fit to the experimental data with no Gaussian noise, while the grey uncertainty band is generated from the middle 68% of the curves. In comparison with the extracted Sivers function in Fig. 7 from SIDIS+DY data at low energy, this fit leads to much larger uncertainty band for the  Figure 16. Left: The COMPASS deuteron target measurement [44] for π + , π − , K + , K − , and K 0 from top to bottom, and as a function of x B (left), z h (middle), and P h⊥ (right). Right: HERMES proton target measurement [41] π + , π 0 , π − , K + , K − , and (π + − π − ) from top to bottom, and as a function of x B (left), z h (middle), and P h⊥ (right). The data is plotted in red along with the total experimental error. The central curve in blue as well as the uncertainty band in gray are generated using the result from fit 2b in Tab. 1.
Sivers function. The size of the Sivers functions in this fit is also significantly larger. This is of course due to the much larger asymmetries for W/Z bosons measured at RHIC. We note that we have also checked the extracted asymmetry in fit 2b against the new HERMES data in [77]. We find that there is very strong agreement between this extracted asymmetry and the new data. In Figs. 16, 17, and 18, we plot the theoretical curve of this fit against the low energy experimental data for SIDIS Sivers asymmetry. In Fig. 19, the comparison with the COM-PASS Drell-Yan data is presented. While the theoretical uncertainties are much larger than the previous extraction, the fitted asymmetry still describes the this subset of the data very well. Finally, in Fig. 20, we plot the fitted asymmetry to the RHIC data. We find that in this scheme, the size of the asymmetry for the central fit can now be up to 5%. Overall, this scheme describes the RHIC data much better than the previous extraction. The future RHIC data with much smaller experimental uncertainty will for sure help to reduce the theoretical uncertainties in the extracted Sivers functions, as well as the Sivers asymmetries computed based on these Sivers functions.

Predictions for the EIC
As we have seen in the previous sections, the choice of DGLAP evolution scheme used for the evolution of the Qiu-Sterman function greatly affects the quality of the fit when considering data at large hard scales. While this issue currently presents difficulties for performing a global extraction of the Sivers function, this effect also presents an opportunity at the future EIC. The EIC will be capable of performing high precision measurements of transverse spin  Figure 19. COMPASS Drell-Yan measurement for π − -p collision [46] as a function of q ⊥ , Q, x F , x N , and x π from left to right. The central curve and uncertainty band are generated using the result from fit 2b in Tab. 1. asymmetries at a large range of scales. Experimental data which are collected over these large range of scales can be used to study DGLAP evolution effects of the Qiu-Sterman function. On the left side of Fig. 21, we plot our prediction for the Sivers asymmetry in SIDIS on a proton target as a function of x B at √ S = 105 GeV, z h = 0.25, q ⊥ /Q = 0.2 at Q 2 = 5, 50, 500 GeV 2 for π + , π − , K + , and K − production. In this figure, we have plotted our prediction for the low energy fit (fit 1 in Tab. 1) in blue, and the global fit (fit 2b in Tab. 1) in gray. While this prediction demonstrates the x-dependence of our fits, in order to demonstrate the k ⊥ -dependence of our fitted Sivers function, we also make a prediction as a function of P h⊥ on the right side of Fig. 21. In this figure, we have used the same kinematics as the left side except that we take x B = 0.2. We see from these curves that the predicted asymmetry for π − and K − production is small. This behavior is expected because of the suppression by the fractional charge e 2 d for the d-quark Sivers function, as well as the cancellation that occurs between the u and d-quarks. On the other hand, we predict an asymmetry of a few percent for π + and K + production in this kinematic region.
We see in these plots that the theoretical curves generated from fit 1 and fit 2b are very similar at Q 2 = 5 GeV 2 . This behaviour occurs because the suppression factor, N in Eq. (4.5), is close to one at small energies. However, at Q 2 = 500 GeV 2 , the theoretical curves generated from fit 1 and fit 2b can differ by a few percent. This effect presents a great opportunity at the future EIC. Since measurements at large values of Q 2 are sensitive to the DGLAP evolution effects of the Qiu-Sterman function, these data may prove useful in phenomenological studies of this evolution. At the same time, these future measurements at the EIC could provide additional statistics for high energy data which will prove useful in reducing the theoretical uncertainties for the extraction of the Sivers asymmetry at large energy scales.

Conclusions
In this paper, we have performed extractions of the Sivers function for the first time at the NLO+NNLL order. We first perform an extraction from the Sivers asymmetry data measured in SIDIS at HERMES, COMPASS and JLab, and in Drell-Yan lepton pair production at COMPASS. Using this first extraction, we generate a prediction for the Sivers asymmetry of W/Z boson at RHIC kinematics and compare with the experimental data. We find that while the SIDIS and COMPASS Drell-Yan lepton pair production data is very well described by our extraction, that our theoretical curve is much smaller than the RHIC data. We study in great detail the impact of the RHIC data and their implications. For such a purpose, we perform a fit in which we introduce a weighting factor of ∼ 13 for the RHIC data, so that the RHIC data and the rest of the experimental data sets are equally weighed in the calculation of the χ 2 . We study how RHIC data are sensitive or insensitive to the non-perturbative parameters in the Sivers function parameterization. In addition, we study in detail the dependence on the choice of the scheme used to perform the DGLAP evolution of the Qiu-Sterman function, the collinear counterpart that enters the TMD evolution formalism for the Sivers function.
We investigate the impact of two DGLAP evolution schemes which are commonly used in the extraction of the Sivers function. We find that the scheme which treats the evolution of the Qiu-Sterman function the same as the unpolarized parton distribution function, is better suited for describing the experimental data at RHIC. Using DGLAP evolution scheme, we perform for the first time a global extraction of the Sivers function and find that this scheme improves the description of the RHIC data. While our first fit describes the low energy data extremely well, our second fit describes the RHIC data much better than the first. However, due to the large experimental uncertainties at RHIC, we find that the globally extracted Sivers function has large theoretical uncertainties. We expect the forthcoming RHIC experimental data on W/Z Sivers asymmetry with large statistics and reduced experimental uncertainties would help us better constrain the Sivers function and its evolution. In addition, we make predictions for Sivers asymmetry at the future Electron Ion Collider (EIC). We find that with large range of hard scale Q to be probed at the EIC, the effects due to the DGLAP evolution of the Qiu-Sterman function can be extremely pronounced. Such measurements would present a great opportunity for testing such effects.
Upon publication, the extracted Sivers functions from both fits in this paper will be made available open source at the following link: https://github.com/UCLA-TMD/TMD-GRIDS/tree/EKT2020.

A Wilson Coefficient Functions
The scale dependent TMDPDF quark to quark and gluon to quark Wilson coefficient function is given by [81][82][83] where in these expressions, we have used the short-hand The quark to quark coefficient function for the TMDFF is given by the relation while the quark to gluon Wilson coefficient function for the TMD FF is given bŷ In these expressions, we have introduced the standard collinear splitting kernels Finally, the coefficient function for the quark-Sivers function is given bȳ

B TMD evolution ingredients
The following expansions, numbers, etc, can be found in the 2013 PDG [84]. First of all, we need the expansion of the strong coupling in terms of Λ QCD : where x = ln µ 2 /Λ 2 QCD , and the coefficients of the beta-function are given as Since we want the resummation up to NNLL, we take the expansion of α s with β 0 , β 1 and where the coefficients up to NNLL are given by On the other hand, in order to describe the perturbative TMD evolution, we want to analytically solve the integral where the coefficients of the perturbative expansions of the anomalous dimensions can be found in the below.

B.1 Integration at NLL accuracy
For this order we take γ 0 , Γ 0 , Γ 1 , β 0 and β 1 . Thus we have: The final result is then Be careful with the number of active flavors. The number of flavors for the x U that appears inside the integrand is fixed and depends on the value of µ U . However, depending on the hierarchy between µ L , µ U and m b we might have to split the integral in several pieces, and in that case, when we substitute the limits of the integral, x L and x U , they would have different numbers of active flavors (still the x U that already appeared in the integrand before the substitutions just depends on the value of µ U ).

C Evolution of the Hard Matching Coefficient
The evolution of the hard matching coefficient C V , which is related to the usual hard function as H = |C V | 2 , is given by where the cusp term is related to the evolution of the Sudakov double logarithms and the remaining term with the evolution of single logarithms. The exact solution of this equation is where we have used that d/dlnµ = β(α s ) d/dα s , where β(α s ) = dα s /dlnµ is the QCD βfunction.
Below we give the expressions for the anomalous dimensions and the QCD β-function, in the MS renormalization scheme. We use the following expansions: (C.6) The coefficients for the cusp anomalous dimension Γ cusp are Γ 0 =4C F , T F n f , The anomalous dimension γ V can be determined up to three-loop order from the partial three-loop expression for the on-shell quark form factor in QCD. We have (C.8)