Assessing signals of TMD physics in SIDIS azimuthal asymmetries and in the extraction of the Sivers function

New data on the Sivers azimuthal asymmetry measured in semi-inclusive deep-inelastic scattering processes have recently been released by the COMPASS Collaboration at CERN. Their increased precision and their particular binning, in terms of $Q^2$ as well as $x$, motivates a new extraction of the Sivers function, within the framework of a simple and transparent parametrization. Signals of TMD effects visible in the Sivers asymmetries are critically assessed. A thorough study of the uncertainties affecting the extracted Sivers function is presented, including the low-$x$ and large-$x$ regions.


INTRODUCTION
A successful study of the 3D-structure of nucleons depends on our ability to efficiently extract transverse momentum dependent (TMD) parton densities from experimental data. These functions encode nonperturbative information on the inner composition of hadrons in terms of their elementary constituents, and on the dynamical mechanisms which confine partons inside hadronic states.
The extraction of TMD parton distribution functions (PDFs), however, is a complex task that involves a series of steps not often free of pitfalls. On one side we have a "theory" (QCD in our case) which in principle provides a full description of the underlying physics relevant to the dynamical processes considered. Often, however, theory cannot be applied directly because it is not exactly solvable, or incomplete, or simply impractical.
Beyond the unpolarized TMD, the most interesting and studied polarized TMD-PDF is perhaps the Sivers function [1,2], which correlates the motion of unpolarized partons with the spin of the parent nucleon, and can be accessed through azimuthal asymmetries in polarized Drell-Yan (DY) and Semi-Inclusive Deep-Inelastic Scattering (SIDIS) processes. Remarkably, the Sivers function is predicted to have opposite signs in these two processes [3,4]. This sign change, in fact, has been the focus of several phenomenological analyses, although none has been totally conclusive [5][6][7][8].
Very recently the COMPASS Collaboration has presented a new re-analysis of their SIDIS measurements [9], based on a two-dimensional binning: the Sivers asymmetries are presented as functions of the kinematic variables, x, P T and z, one at a time, for four regions of the photon virtuality Q 2 . These Q 2 ranges correspond to the four regions of the di-muon mass explored in the ongoing analyses of the COMPASS Drell-Yan measurements [8].
The large number of data, together with a considerably increased precision and a finer binning in Q 2 as well as in x, poses the question of whether one can extract the Sivers function within a full QCD scheme, as that defined in Ref. [10]. In this theoretical framework, one must determine, based on the data, an input function which can be interpreted as the Sivers function at a given initial scale, and the non-perturbative function g K , which is responsible for the broadening of the TMDs as a function of the scale. Moreover, it requires a full knowledge of the unpolarized TMD PDF and FF, for which studies are still at a very early stage [11][12][13]. Further complications arise from considerations as those discussed in Refs. [14,15], where it has been suggested that at the kinematics of the current data, the errors of factorization may not be completely under control. Note that when performing a fit on experimental data, it may happen that large theoretical errors are "absorbed" by the model-dependent parts of the TMDs, which can make the interpretation of the analysis results more problematic.
For a reliable extraction of the Sivers function, it is also crucial to understand the extent to which different aspects associated with TMDs are visible in the data. Given the complications discussed above, it makes sense to use a bottom-up approach, that addresses questions that regard only the data and the information which can be inferred from them. In this paper, we will focus on these particular issues and study the extent to which effects relevant to TMD physics are likely to be observed in the existing sets of SIDIS experimental data. For this analysis, in fact, we will model the Sivers function using a parameterization similar to that used in our past work [7,16,17], but we will relax the assumption that the Sivers functions should be parameterized in terms of the corresponding unpolarized TMD PDF. The chosen parametrization will be simple but flexible enough to allow for a realistic evaluation of the uncertainties affecting the extracted functions.
Finally, as the new COMPASS experimental data are separated in different Q 2 regions, it will be interesting to compare the results obtained by using different Q 2 evolution schemes. We will analyze the scale dependence of the Sivers function predicted by three assumptions: the no-evolution case, where the Sivers function does not depend at all on the scale Q, the collinear twist-three approach, where the Sivers function varies with Q only through the kinematic variable x, and a TMD-like scheme, in which the Q 2 evolution proceeds through a modification of the width of the Sivers k ⊥ -distribution with varying Q. The new functional form of the parameterization we have introduced, independent of the unpolarized TMD PDFs, is particularly suited to be applied to the full TMD-evolution scheme.
The paper is organized as follows. In Section II we will describe the general framework used in our analysis and the parametrization which will be adopted. In Section III the two main best fits performed to extract the Sivers function will be presented and illustrated in detail, together with a thorough analysis of the corresponding uncertainties and comparisons to the experimental measurements. In Section IV some results based on the error projections of a new run of the COMPASS II experiment with polarized deuterium targets [18] will be presented. The uncertainties obtained using present data and new projected errors will be compared. In Section V we will comment on how possible signals of scale dependence can be detected in the examined SIDIS data. Final remarks and conclusions will be drawn in Section VI.

II. GENERAL STRATEGY
For the current study, we adopt a model for the Sivers function similar to that of Refs. [7,16,19,20]. We assume a factorized form for the x and k ⊥ dependences, and use a Gaussian model for the latter As mentioned above, the main difference between this parametrization and those used in previous analyses [16,19,21] is that in Eq. (1) the x-dependent part of the Sivers function, for each flavour, is no longer parametrized in terms of the corresponding unpolarized PDF. In the past, when data were scarce and affected by rather large experimental uncertainties, this parametrization provided a useful input to allow for a successful extraction of the Sivers function even though the information contained in the experimental data was quite incomplete. It also had the advantage of ensuring the automatic fulfillment of the required positivity bounds. In the current study, however, we relax this assumption in order to test in the most agnostic possible way aspects of the data related to TMD physics, like flavour separation and scale dependence. Our approach is also flexible enough to allow for a realistic determination of the uncertainties in the extraction of the Sivers function.
Furthermore, in the new model the width of the Sivers function is not written in terms of the width of the unpolarized TMD PDF; instead, we parametrize the Sivers function directly in terms of its TMD width, k 2 ⊥ S . Note also that the parameterization of Eq. (1) has been arranged in such a way that its first moment assumes a much simpler form, namely where the rightmost equation provides the relation of the first moment with the Amsterdam notation. For the unpolarized TMD PDFs and FFs we use the same functional forms as that adopted in Ref. [7], namely where f q/p (x) and D h/q (z) are the usual unpolarized PDFs and FFs, which we will take from the CTEQ6l [22] and DSS [23] leading order (LO) sets, respectively; k 2 ⊥ and p 2 ⊥ are the widths of the corresponding TMD distributions, which will be fixed according to the values extracted in Ref. [11], as we will explain in detail below. Although not explicitly indicated, our model for the unpolarized TMD PDFs and FFs depend on Q 2 , according to Dokshitzer, Gribov, Lipatov, Altarelli, Parisi (DGLAP) equations [24][25][26].
The SIDIS Sivers asymmetry, defined as can be expressed, within this framework, through the following relations We will examine different possible scenarios: our starting point is to set α q = 0 and k 2 ⊥ S = constant in Eq. (1) and to consider only the contributions from u and d flavours. This provides a reference best fit that will be used as a baseline for comparison. Then, we analyze the data with different modifications of the above reference parametrization, each of them properly devised to address different aspects regarding the sensitivity of the data to some chosen features.
Specifically, we will investigate to which extent the present experimental data support the flavour separation of the Sivers function and, in turn, how we can estimate its uncertainties in the low-x region, where the sea contributions are expected to become dominant. Moreover we will explore the sensitivity of the experimental measurements to Q 2 and x correlations, to Q 2 dependence and, possibly, to TMDevolution effects.

A. A closer look to data
Our fits will include all experimental data presently available on the Sivers asymmetries in SIDIS processes: from the HERMES Collaboration for π ± , π 0 and K + SIDIS production off a proton target [27], from the COMPASS Collaboration for π ± , K 0 and K + on LiD [28] and for h ± on NH 3 targets [9] with z > 0.2, which correspond to a very recent reanalysis of COMPASS 2010 measurements using a novel Q 2 binning, and finally from JLab data on 3 He target [29]. We will not include the K − data, as they are mainly driven by the sea contributions of the Sivers function; as explained below in this analysis sea and valence will not be separated.
For all experiments, these data are provided as functions of x, P T and z kinematic variables, with the exception of JLab data which provides only x dependent asymmetries. We will not include the z-distributions, as in our model the z dependence of the asymmetries is essentially fixed by the FFs, and it has essentially no sensitivity to our free parameters. In order to estimate the uncertainties in our extractions, we carefully explore the parameter space and consider a 2σ confidence level (C.L.), corresponding to a coverage probability of 95.4%. We then consistently accept parameter configurations that render a value of χ 2 in the range [χ 2 min , χ 2 min + ∆χ 2 ]. Note that the value of ∆χ 2 depends on the number of parameters considered and will be reported in the tables corresponding to each fit.
In order to extract the Sivers function, the first necessary ingredients are the unpolarized TMD PDFs and TMD FFs. This poses a big complication, since knowledge of these TMD functions from SIDIS data is very limited. In our study, we use the Gaussian functional forms of Eqs. (3) and (4) for the unpolarized TMDs, with the minimal parameters obtained in [11], where HERMES and COMPASS multiplicities were analyzed and fitted within the same scheme adopted here. There, it was found that COMPASS and HERMES multiplicities could be well reproduced by a simple Gaussian model, like that of Eq. (7), using only two free parameters k 2 ⊥ and p 2 ⊥ , i.e. the widths of the TMD PDF and the TMD FF, respectively. We stress that in Ref. [11] data sets from the two experiments had to be fitted separately and no simultaneous extraction was possible. In fact, it is likely that a simultaneous extraction can only be achieved by fully accounting for the highly non-trivial dynamics encoded in TMD-evolution equations, and properly dealing with the delicate interplay between perturbative and non-perturbative regimes. This is indeed a topic of ongoing research [14,15,[30][31][32]. For our current purposes all what is needed is that we consistently use the results of Ref. [11] for each individual experiment. Furthermore, we use the unpolarized widths from the HERMES extraction for the JLab Sivers asymmetries, since these two experiments were shown to be compatible in Ref. [11].
In order to illustrate this last point, we performed a simple test, where we evaluated the effects of using different Gaussian widths for the computation of unpolarized HERMES and COMPASS cross sections, i.e. our asymmetry denominators. We compared two hypotheses: i) using the same unpolarized widths for both HERMES and COMPASS asymmetries (namely the widths extracted from HERMES), ii) using different unpolarized widths corresponding to each experiment. In each case, we performed three fits on π + production from a proton target, considering only the u contributions (all other flavours being set to zero): HERMES only, COMPASS only, HERMES+COMPASS simultaneously. Results are shown in Fig. 1, where scatter plots for the parameter space are displayed, at 2σ C.L. The left panel shows that choosing the same unpolarized widths for the HERMES and COMPASS data sets results in fits that populate different regions of the parameter space (red and blue areas). In fact, the almost completely disjoint sets signal some tension. As a consequence, the combined fit (black area), although still giving a good value of the χ 2 , will have to compromise, rendering values of the normalization parameter N u which end up being "half way" between the red and the blue regions. In contrast, the regions in the parameter n. of data points = 220 One flavour fits ( . The left column indicates the flavour contribution considered in each fit (all others being set to zero). The top panel shows how the u flavour contribution dominates the effects visible in the data. The bottom panel shows the improvement on the description of the data when including one more flavour to the leading u contribution. We highlight the chosen configuration for our study: u and d contributions only. Note that adding more parameters to disentangle the sea, would put the analysis procedure at risk of over-fitting.
space explored in the right panel all overlap, visibly reducing the tension. This supports our choice to use the appropriate unpolarized Gaussian widths for each experimental set in our analysis.
This preliminary investigation illustrates how having a good knowledge of the unpolarized TMD distribution and fragmentation functions is of crucial importance for the analysis and extraction of any polarized observables. In this particular case, the fact that two different experimental data sets seem to point to different unpolarized Gaussian widths could be attributed to many different origins; possibly, a signal of TMD evolution effects. We will not get into this here, but clearly this is an issue which deserves further investigation.

B. Reference fit
The baseline of our analysis is given by Eq. (1), in which we set α u = α d = 0, so that the first moment of the Sivers function simply reduces to Furthermore, we assume the width of the Sivers function, k 2 ⊥ S , to be independent of other kinematic variables and of flavour. This introduces only one extra free parameter.
For all of the cases considered in this article, we will not attempt a flavour separation of sea and valence contributions. In fact, we have tested different hypotheses regarding the flavour content of F our results are shown in Table I, where the left column indicates which flavour component has been included in each fit (all other components being set to zero). As it can be seen in the upper panel of Table I, the u flavour Sivers function represents the leading contribution to the asymmetries. The total χ 2 improves significantly if one more flavour is added to the fit, as shown on the lower panel of Table I. Any further addition of different flavour contributions will not improve the quality of the fit, making convergence to the minimum more cumbersome and exposing us to the risk of over-fitting. For our analysis, we use the configuration that renders the smallest minimal χ 2 , i.e. we directly parametrize the total u and d flavours as follows: This results in a fit with a total of 5 free parameters, N u , N d , β u , β d and k 2 ⊥ S . We call this "reference fit". The role of the sea contributions, which are expected to be relevant at small x where the behaviour of the Sivers function is mainly driven by the α q parameters, will be addressed in Section III C.
As explained above, the k ⊥ widths of the unpolarized TMDs are fixed according to the values extracted in Ref. [11], to make sure that the unpolarized cross sections appearing in the asymmetry denominator reproduce well the measured multiplicities for both HERMES and COMPASS experiments. In this first, simple fit no Q 2 evolution is applied to the Sivers function, and the corresponding plots are labeled by "no-evolution"; the extracted function will therefore represent the Sivers function at the average Q 2 scale of the experimental data. Table II shows the values of the free parameters as determined by our best fit, together with the minimal values of the χ 2 and the total number of data points included. The errors reported in Table II are MINUIT errors, corresponding to 2σ C.L., i.e. to a coverage probability of 95.4%.
The top panels of Fig. 2 show the χ 2 tot profiles as functions of the parameters N u and N d . In these two plots the correlations between the parameters N u and N d are colour-coded: yellow corresponds to lowest, green to intermediate and purple to highest allowed values of N u (top left panel) and N d (top right panel). It is evident that these profiles are quite well approximated by a quadratic function, confirming that the Hessian method adopted to evaluate the errors on the parameters is reliable. In fact, the errors reported in Table II are well in agreement with the uncertainties on the free parameters that can easily be inferred by looking at the scatter plots. The lower panels of Fig. 2 represent the correlations between the parameters N u , β u (lower left) and N d , β d (lower right). Here it is the corresponding χ 2 tot which is colour-coded: yellow corresponds to the lowest χ 2 tot values, green to intermediate and purple to the highest χ 2 tot values. As expected from a fit consistent with a Hessian approximation, the correlations among parameters cover regions of reasonably regular, ellipsoidal shapes. Fig. 3 shows the χ 2 tot profile of the Sivers k ⊥ width, k 2 ⊥ S , and its correlation with N u (color coded). Also in this case the uncertainties indicated by the scatter plots of the parameter space are perfectly consistent with the errors estimated by adopting the Hessian approximation, reported in Table II.
Plots showing the u and d Sivers functions and their estimated uncertainty bands, as extracted in this fit will be shown below. They will be extensively discussed in Section III C.
The reference fit, with 5 free parameters, is able to reproduce all the existing SIDIS experimental measurements. Moreover, it provides a successful extraction of the Sivers function as well as a reliable estimate of the uncertainties, over the kinematic region covered by the bulk of experimental data (i.e. approximately 0.03 < x < 0.3). Below this region, where only very few data points are present, the error bands from the reference fit are at risk of being artificially small. In section III C, we will consider the case where the α parameters in Eq. (1), which regulate the low-x behaviour of the Sivers function, may be different from zero. As we will discuss, this provides error bands that better reflect the amount of information which can be inferred from data.

C. Low-x uncertainties
Starting from the reference fit described above, which represents the basis for all further studies presented in this paper, we now move on to explore in more detail the low-x kinematic region. To do this, we perform a different fit in which we allow our parametrization to become more flexible at small x by including two extra parameters, α u and α d , in the following way This best fit will be referred to as the "α-fit". Table III shows the χ 2 and the values of the parameters obtained in this case. As it is immediately evident, the value of χ 2 dof is unchanged, therefore the overall α fit -no evolution χ 2 tot = 211.5 n. of points = 220 χ 2 dof = 0.99 n. of free parameters = 7  Table II, the value of χ 2 tot remains practically unchanged. However, the uncertainty on the free parameters increases considerably. This generates much larger uncertainty bands in the low-x region, as shown in Fig. 5.
FIG. 4: Parameter space scatter plots for the α-fit, which includes the αu and α d free parameters. The regions displayed correspond to our estimate of 2σ C.L. error band. Notice that the uncertainties on the parameters which can be inferred from the scatter plots are much larger than the errors reported in Table III. quality of the fit does not improve. Moreover, the central values of the free parameters are extremely close to those obtained in the reference fit. This suggests that the experimental data currently used are not sensitive to the particular choice of value for α u and α d . Consequently, further constraining the low-x behaviour of the Sivers function seems at the moment unlikely.
One can notice, however, a sizable increase of the parameter errors. This has an effect on the uncertainty bands of the extracted Sivers function, as shown in Fig. 5, where the light-blue bands correspond to the reference fit, while the wider, gray bands refer to the "α-fit".
The modifications on the parameter space induced by adding the two α parameters is shown in detail in (top-right). Contrary to the reference fit, these profiles are very far from resembling quadratic functions, and therefore the Hessian method adopted to evaluate the errors on the parameters cannot be trusted. The MINUIT errors reported in Table II, in fact, largely underestimate the uncertainties on the free parameter determination: by looking at the plots in Fig. 4, one can easily see that, to 2σ C.L., N u can go as low as 0.1 and as large as 4.0, over a very asymmetric range. Similarly for N d , which can span over an even larger range, from 0 to −45, on an extremely asymmetric range. A clear indication, however, is given on the sign: N u is positive and N d is negative, signaling a preference of the data for a positive u and a negative d Sivers function.
In the upper panels of Fig. 4 the correlations (N u ,α u ) and (N d , α d ) are colour-coded: the very evident structure in bands of the same colour points to extremely strong correlations. This becomes even more explicit in the bottom panels of Fig. 4, where we show N u vs. α u (bottom left) and N d vs. α d (bottom right). In these scatter plots, the expected ellipsoidal shapes are replaced by very thin and stretched distributions, which indicate that an extremely large number of equally good fits can be obtained provided The results obtained from the reference fit and the α-fit are compared to the HERMES measurements of the SIDIS Sivers asymmetry for π ± production off a proton target [27] (upper panels), to the COMPASS measurements of the SIDIS Sivers asymmetry on a LiD target [28] for π ± production (middle panels), and to the JLab data for π ± production on a 3 He target [29] (bottom panel). Here we show the x dependence only. The shaded region corresponds to our estimate of 2σ C.L. error band. The light-blue bands correspond to the uncertainties of the reference fit (only N u(d) and β u(d) free parameters), while the (larger) gray bands correspond to the uncertainties of the α-fit, which includes also the α u(d) parameters.
N , α (and β) are in the appropriate ratio among each other. In other words, even very large values of the α and β parameters can result in an acceptable χ 2 , provided the corresponding N parameter is adequately large in size. Conversely, low values of α and β are also equally appropriate if N is small enough. The strong correlations introduced by α, in fact, make it cumbersome to find a good fit by a simple minimization procedure. Nonetheless, the study of the parameter space including the α parameters allows for a more realistic estimate of the uncertainty bands in the small-x region. This is shown in Fig. 5, where the gray shaded areas represent the uncertainty bands corresponding to the α-fit, while the light-blue bands represent the uncertainties corresponding to the reference best fit. Clearly, the two fits have very similar bands over the . The shaded bands correspond to our estimate of 2σ C.L. error for α-fit case. In both panels, the gray bands correspond to the fit which includes the new COMPASS-2017 data [9] for h + and h − production off NH3 target, while the meshed areas correspond to the uncertainties obtained when the COMPASS-2015 data from the older analysis [33] are included. The lower panels show the relative errors, given by the ratio between the upper/lower border of the uncertainty bands and the best-fit curve for the reference fit. region 0.03 < x < 0.3, while the α-fit uncertainties grow larger outside this range, where experimental data are less dense. Notice that the Sivers width, k 2 ⊥ S , is not significantly affected by this strong broadening of the uncertainty bands: its central value remains the same (see Tables II and III), while error bands show no significant change, as it is clearly evident in the bottom panels of Fig. 5, where the u and d Sivers functions are plotted vs. k ⊥ at x = 0.1. Fig. 6 shows the results obtained from this reference fit compared to older data, from HERMES-proton (top panels) and COMPASS-deuterium (middle panels), which have historically been present in our previous fits, together with JLab-neutron measurements (bottom panels).
The bands corresponding to the reference best fit are shown in light-blue. The enlargement of the gray bands for the α-fit provides a more sensible estimate of the uncertainties at low x. In fact, as seen in the central panels of Fig. 6, the agreement of the light blue bands with the deuteron data seem to deteriorate at small x, while the gray bands corresponding to the α-fit improve the compatibility with these experimental measurements. Note that, since separating valence and sea contributions is not possible with the current data, the effect on the uncertainties introduced by allowing α = 0 also reflects our ignorance about the sea contributions.
This supports the need to learn more about the Sivers function in the low-x region and, in turn, about its sea contributions. In fact, this is one of the main tasks of the future Electron-Ion-Collider (EIC) [34], which is planned to be built in the next few years in the USA. Besides the clear benefits of an EIC to resolve the sea of the Sivers function, we stress the importance of the deuteron target measurements as those performed by COMPASS [28]. Recall that the α-fit uncertainties only have a significant impact in the description of these data in the low-x regime, where the reference fit delivered uncomfortably small uncertainties. Thus, an improvement on the statistics for these measurements may prove very useful in constraining the low-x regime, rendering information about the sea distributions (see section IV).
In Fig. 7 we compare the results obtained from the reference fit to the newest COMPASS data on the SIDIS Sivers asymmetry on a NH 3 target [9], for h + production, binned in four values of Q 2 (the average value of Q 2 corresponding to the bin is indicated on each panel). Only the x and the P T dependences are shown, as the z dependences are not included in the fit. However, we have checked that all z-distributions are always successfully reproduced. The shaded regions correspond to our estimate of 95.4% C.L. error band.
Finally, we compare the results obtained from our fit which uses the newly re-analyzed data by the COMPASS Collaboration [9] to those obtained using the older set of COMPASS-proton data [33]. These newly released data sets belong to the same measurements (2010 run), but they differ in the way they are binned. In fact, in their more recent analysis the data are binned in x and P T as well as in four bins of Q 2 , the same bins in which the Drell-Yan analysis is being carried out. As shown in Fig. 8, some reduction of the uncertainty bands is obtained when using the COMPASS-2017 data set, indicating that an increased degree of information is reached by applying the new binning. An important feature of the new binning is the separation of different ranges of Q 2 . This, for the first time, allows to explore the possibility of scale dependence in the Sivers function. We will address this in Section V.

D. Large-x uncertainties
The study of large-x uncertainties is indeed very delicate. At present, as shown in this analysis, the Sivers function is largely unconstrained in the range from x ∼ 0.3 up to x = 1.0.
In this region the Sivers function should approach zero with the only theoretical constraint given by the positivity bound, |∆ N f q/p ↑ | ≤ 2f q/p , which should hold for any flavour q and at every value of x and k ⊥ . Notice, however, that also the integrated unpolarized PDFs are largely undetermined at large-x, undermining the significance of any phenomenological application of the positivity bound itself.
To make the large-x uncertainties more visible, in the middle panel of Fig. 5, we show the bands corresponding to the relative uncertainties, i.e. the ratios between the upper/lower border of the uncertainty bands and the best-fit curve for the reference fit, at each value of x and for any given flavour.
In a similar way, in Fig. 9, we show the ratio |∆ N f For comparison, we also display the central line of the previous extraction of the Sivers functions from Ref. [7]. As expected, in the range 0.03 < x < 0.3, the agreement between the two extractions is acceptable. As one goes outside of this region, however, the two extractions exhibit more distinct behaviours. While this does not compromise compatibility at large x, as both central lines fall within the error bands, differences at low x are more dramatic. This should serve as a warning that model dependence has an important effect outside the bounds of experimental information.
Future measurements at JLab12 [35] will be able to shed some light on the large-x kinematic region. These measurements will give a crucial contribution in the extraction of the unpolarized TMD parton distribution functions as well as the Sivers functions in the large-x range, and will give us a much clearer signature for the flavour separation of the valence contributions. The shaded bands correspond to our estimate of 2σ C.L. error. The light-blue bands show the uncertainties corresponding to our reference fit (see Table II). The red (meshed) bands correspond to the uncertainties estimated by using the same model, with the projected experimetal errors of the future COMPASS run on deuteron target [18].

IV. COMMENTS ON THE PRECISION OF DEUTERON TARGET DATA
As mentioned above, the d Sivers function is poorly determined by the existing Sivers data, in spite of the fact that it should be constrained by the identification of the final state hadrons. Possibly this can be traced back to the large u-dominance of SIDIS on proton targets.
The COMPASS collaboration has recently proposed a new run of their SIDIS measurements on a deuteron target, with increased statistics and precision [18]. It is therefore interesting to evaluate the impact of such a measurement on the extraction of the first moment of the u and d Sivers function, using the projected errors for the proposed 2021 deuteron run, as reported in Ref. [18].
Our results are shown in Fig. 10. The 2σ error bands marking the 95.45% C.L. for the first moments of the u and d Sivers functions, obtained with the reference fit (here labeled by "current") are shown in light-blue. The bands obtained when adding to the data set the projected errors on the asymmetries of the new deuteron run are shown in red, and labeled "projected". The plots for the first moments (bottom panels) show the relative uncertainty, i.e. the ratio between the upper/lower border of the uncertainty bands and the best-fit curve for the reference fit. As expected, the new deuteron run will have a small impact on the u-quark first moment of the Sivers function. On the contrary, the reduction in the error band for the first moment of the Sivers function for the d-quark is considerable, and is about a factor 2 for x < 0.1.
This new COMPASS run will therefore lead to a remarkable improvement of our knowledge on the other flavour contributions of the Sivers function, besides the already well constrained u.

V. SIGNALS OF SCALE DEPENDENCE
In all the results presented above, no Q 2 dependence of the Sivers function was considered. An important aspect of the new COMPASS binning is that it separates different ranges of Q 2 , which poses the question of whether one can distinguish different assumptions about scale dependence. To test how this scale dependence can affect our analysis, we will consider two different approaches which will serve us as comparisons: on one side, we will adopt a collinear, twist-3 evolution scheme based on Refs. [36][37][38][39]; on the other side, we will apply to the Sivers function a TMD-like Q 2 evolution similar to that described in Ref. [17].
In the collinear higher-twist evolution framework, the correlation between spin and transverse momentum is included into the higher-twist collinear parton distributions or fragmentation functions. These functions have no probabilistic interpretation: they are generated as quantum interferences between a collinear active quark state in the scattering amplitude and a collinear quark-gluon composite state in its complex conjugate amplitude. There are no intrinsic k ⊥ in this case, which are integrated over, and the Collinear twist-3 evolution χ 2 tot = 201.5 n. of points = 220 χ 2 dof = 0.94 n. of free parameters = 5 ∆χ 2 = 11.3 HERMES k 2 ⊥ = 0.57 GeV 2 p 2 ⊥ = 0.12 GeV 2 COMPASS k 2 ⊥ = 0.60 GeV 2 p 2 ⊥ = 0.20 GeV 2 Nu = 0.39 ± 0.08 βu = 3.55 ± 1.26 N d = −0.65 ± 0.27 β d = 4.77 ± 3.41 k 2 ⊥ S = 0.33 ± 0.14 GeV 2 TABLE IV: Best fit parameters and χ 2 values for the collinear twist-3 evolution case. The parameter errors correspond to 2σ C.L. Notice the reduced value of χ 2 tot w.r.t. that of the reference fit in Table II.
evolution in Q 2 occurs only through x. In other words, twist-3 PDFs and FFs evolve in Q 2 by changing shape in x.
In the TMD factorization approach, spin asymmetries are generated by spin and transverse momentum correlations between the identified hadron and the active parton. This correlations are embedded in the TMD parton distribution or fragmentation functions, which can be interpreted as probability densities.
Here the Q 2 evolution affects the x dependence as well as the shape in k ⊥ .
Although they are defined in different contexts, TMD and collinear quark-gluon correlation functions are closely related to each other. In particular, the first k ⊥ -moment of the Sivers function is related to the collinear, twist-3 quark-gluon correlation function T q,F (x, x) [39]. As the evolution equations for T q,F (x, x) are known, we can adopt them in our study to render the Q 2 dependence of the Sivers function, from the initial scale Q 2 0 = 1.2 GeV 2 (which coincides with the lowest Q 2 of the experimental data included in our best fit) to the Q 2 corresponding to each specific data point at which the asymmetry is evaluated. To implement the collinear twist-3 evolution we use the HOPPET code [40], appropriately modified to include the kernels corresponding to the Sivers function [41]. Notice that while we do not include off diagonal terms in the twist-3 evolution case, this approximation is enough for our purposes: we will test whether the existing data can distinguish between an approach with no evolution (reference fit) and another where some scale dependence appears in the first moment of the Sivers function.
In Fig. 11 the u and d Sivers first moments extracted in the reference best fit with no Q 2 evolution (solid, black lines) are compared to those obtained by applying the collinear, twist-3 evolution described above (blue lines). Three values of Q 2 are shown: Q 2 = 1. The results corresponding to the reference fit with no Q 2 evolution (solid, black line) are compared to those obtained by applying a collinear twist-3 evolution (blue lines), as described in the text, for three values of Q 2 : 1.2 GeV 2 (long-dashed), 3.5 GeV 2 (short-dashed) and 40 GeV 2 (dotted). The bands correspond to the reference fit with no Q 2 evolution.
(dotted), the lowest and largest Q 2 values of the COMPASS measurements and Q 2 = 3.5 GeV 2 (shortdashed), which is approximately the mean value of the full data sample. The corresponding values of the χ 2 and best fit parameters are presented in Table IV.
The Sivers functions extracted in the reference fit are very similar to those obtained using the twist-3 evolution scheme at the experiment average value, Q 2 = 3.5 GeV 2 ; in fact, they are very similar in the region where data constraints are stronger, 0.03 < x < 0.3. Instead, they grow progressively apart when Q 2 is varied to reach its lowest and largest limits. For the u flavour, the Q 2 variation of the first moment due to the collinear twist-3 evolution is actually larger than the uncertainty band corresponding to the no-evolution case. This gives a positive message about the precision of the data, in particular about the new binning of COMPASS asymmetries: signals of collinear evolution could possibly be observed in the experimental data. This does not happen for the d flavour which, as we have already pointed out, is affected by a larger uncertainty. The right panel of Fig. 11 clearly shows that the error band corresponding to ∆ N f (1) d/p ↑ is of the same size (or even larger at small x) of the variation induced by Q 2 evolution. Notice that in this case the whole Q 2 scaling occurs only through x, leaving the k ⊥ part of the Sivers function unchanged.
Finally we turn to the discussion on TMD evolution effects. To extract the Sivers function within a full TMD-scheme, one needs to exploit an "input function", i.e. the value of the Sivers function at the initial Q 2 scale. Then, a TMD factorization scheme as that discussed in Ref. [10], and successively implemented in Refs. [17,42], can be applied to compute the Sivers function at any larger value of Q 2 . One should keep in mind that, within this approach, TMD parton densities change their shape in k ⊥ as Q 2 varies: in particular, their k ⊥ -distributions broaden and dilute as Q 2 increases. While the TMDs themselves (and their first moments) experience variations in their x-distributions too as Q 2 increases, in the azimuthal asymmetries these effects are expected to roughly cancel in the ratio. One complication of this type of analysis is the limited knowledge of the unpolarized functions at the kinematics of the available Sivers asymmetries. In fact, recent studies have suggested that the errors of factorization, at these kinematics, may not be under control [14,15,32,43]. This may affect all of the measurements in SIDIS.
While waiting for further studies to clarify this situation, one may ask questions regarding the data and assess to what extent TMD-evolution effects are visible. To address this, we consider a modification of our model, where we allow for the width of the Sivers function , k 2 ⊥ S , to become a function of Q 2 , according to where g 1 and g 2 are two free parameters to be determined by a best fit, and Q 0 = 1 GeV. The particular choice of Eq. (13), is intended to mimic the main feature of the scale dependence of TMDs, the broadening of the k ⊥ -distribution with variations of Q 2 . In the full TMD definition, this is partly regulated by the non-perturbative, universal function g K (see, for instance, Eq. (44) in [10]).  While a one to one correspondence between g K and our parameter g 2 cannot be made, it serves as a proxy to study the sensitivity of the data to TMD effects. The values of the χ 2 and of the best fit parameters obtained within this model are presented in Table V. Notice that there is no reduction in the value of χ 2 tot w.r.t. that of the reference fit in Table II, although one extra parameter is added to the fit. Fig. 12 shows the correlation between g 1 and g 2 resulting from our analysis; the χ 2 is colour-coded: yellow corresponds to its lowest values while red and purple to the highest accepted values. As it is clearly indicated by this plot, both the parameters g 1 and g 2 are affected by a rather large uncertainty. In particular, the central value of g 1 remains quite close to that extracted in the reference fit while its error increases significantly. The central value of g 2 , instead, turns out to be extremely small, but again affected by a very large uncertainty. Also in this case the two parameters are strongly correlated, and equally good description of the data can be obtained by using rather large and positive values of g 2 provided g 1 is sufficiently small. Paradoxically, even negative values of g 2 are acceptable if g 1 is allowed to grow large and positive, in such a way that the combination (g 1 + g 2 ln Q 2 Q 2 0 ) remains overall positive. Fig. 13 shows the effect of TMD evolution on the k ⊥ -distributions of the Sivers functions for u and d flavours, where blue lines represent the Sivers function at a fixed value of x = 0.1, for three different scales of the data, and the light-blue shaded bands correspond to the uncertainties for the reference fit (no evolution). As expected, the small best fit value of g 2 renders virtually no visible effect in the Sivers function.
Notice that, as in all of our analysis, the widths of the unpolarized TMDs are allowed to be different for each experiment, in accordance with the best description achieved within the approach of Ref. [11]. While more refinements are possible within the same gaussian model, for instance, to include Q 2 dependence as that of Eq. (13), this is unlikely to change the main result of this section, due to the large uncertainties on our g 2 parameter. We remark that while the differences on the unpolarized widths may be attributed to TMD-evolution, this remains as of today an open question. It is quite possible for other effects to play a role (see for instance Refs. [15,44,45]). Regardless of the poor knowledge on the unpolarized TMDs at the kinematics of the Sivers asymmetries, the central point is that additional Q 2 dependence, introduced via g 2 , does not render a result significantly different from that of the reference fit.
These results suggest that in a full TMD analysis, the current Sivers asymmetries will probably not constrain strongly the function g K . On the other hand, due the large uncertainties on g 2 , good compatibility between the Sivers asymmetries and extracted values of g K from other observables are likely to be achieved.

VI. CONCLUSIONS
In this paper we have performed a novel extraction of the Sivers function from SIDIS asymmetry measurements. We have exploited all available SIDIS data from HERMES [27], JLab [29] and COMPASSdeuteron [28], including the new re-analysis of the 2010 run of the COMPASS-proton experiment [9].
The increased statistics and precision of these new sets of data, together with a finer binning in Q 2 as well as in x, has allowed a critical re-analysis of the extraction procedure and its uncertainties. To do so we have adopted a simple and transparent parametric form of the Sivers function, as given in Eq. (1). The aim of this new approach is to attempt an extraction of the Sivers TMD based, as much as possible, on the sole information provided by experimental data. In this framework, it has also been possible to perform a very detailed and accurate study of the parameter space, to provide a reliable estimate of the uncertainties which affect the extracted functions, shedding light on the subtle interplay among experimental errors, theoretical uncertainties and model-dependent constraints.
With our particular choice of parametrization, see Eq.
(1), we started by assessing how much the measured SIDIS asymmetries could tell us about the flavour content of the Sivers function, and on its separation into valence and sea contributions. We found that the existing data can resolve unambiguously the total u-flavour (valence + sea) contribution, while leaving all other flavours largely undetermined (see Table I). We associate to the total d-flavour the additional contribution needed to describe the data, but further investigations possibly with more precise data are necessary. From the statistical point of view, we found that a good configuration was given by a parametric form that considered the contribution of total u and d flavours of the Sivers function, as in Eqs. (10) and (11). Any attempt to separate valence from sea contributions, namely u fromū and d fromd, resulted in a decrease on the quality of the fit, due to a lack of information in the experimental data presently available.
For this analysis we have performed two best fits: the first one was a very basic fit, which we referred to as the "reference fit", based on the most simple parametric form which could reproduce the main features of the Sivers function; the second fit included two extra free parameters, to make the parametrization more flexible in the small-x region, in such a way that possible sea contributions to the u and d flavours could be accounted for, at least partially. Although we could not separate sea from valence contributions within the Sivers first moments, this approach allowed us to obtain a much more realistic estimate of the uncertainties affecting the extracted functions at small values of x. Issues related to the large-x regime, where uncertainties becomes extremely large due to the absence of experimental data, were also discussed. Drawing well-founded conclusions on the low-x and large-x kinematic regimes will only become possible when new experimental information will become available from dedicated experiments which are presently being planned, like the EIC [34], or have just started to run, like the newly upgraded JLab12 [35].
A considerable part of our work was devoted to the study of scale dependence effects. We considered and compared 3 different scenarios: no-evolution, collinear twist-3 evolution and TMD-evolution.
Collinear twist-3 evolution, which proceeds only through x while not affecting the k ⊥ dependence of the Sivers function, was found to be quite fast. In fact, when spanning the range of Q 2 values covered by the experimental data (1.2 GeV 2 < Q 2 < 40 GeV 2 ) the extracted Sivers function shows variations that are larger than the error band for the reference fit. This suggests that the data can help to determine some scale dependence on the first moment of the Sivers function. These results justify a cautious optimism in the possibility of observing this kind of scale dependence in the SIDIS asymmetry experimental measurements.
Signals of TMD evolution, which instead affects mostly the k ⊥ dependence of the Sivers function (effects involving the x-dependence are expected to roughly cancel in the asymmetry ratios) turned out to be more elusive. Our attempts to estimate them resulted in a rather poor determination of the g 2 parameter, which regulates the logarithmic variation of the k ⊥ width with Q 2 , and is intended to mimic the behaviour of the non-perturbative function g K defined in the full TMD approach of [10]. Our best fit delivered a very small value of g 2 , with a large uncertainty. This does not mean that TMD evolution is slow. In fact, within the large uncertainty bands corresponding to this extraction, there is room for quite a large variety of different Q 2 behaviours. Unfortunately, the available experimental information is presently too limited to determine g 2 with a satisfactory precision. While further constraints from a full TMD analysis may help to ease this uncertainty, it is unlikely that, for instance, g K can be constrained via the Sivers asymmetries. However, compatibility with information on g K from some other data sets, such as SIDIS multiplicities, can probably be easily accomplished, as evidenced by our large uncertainties in g 2 .
Finally, we comment on the role of the unpolarized TMDs in the extraction of the Sivers function. As shown in Fig. 1, different assumptions about these functions can alter results significantly. Differently from previous analyses, our choice for the unpolarized functions is based on an approach that better describes most of the available unpolarized SIDIS data. This consideration actually releases the tension on the model of the Sivers function, otherwise encountered when trying to simultaneously fit COMPASS and HERMES data. Among other complications this issue may raise, we realized that this kind of tension can reduce the statistical significance of the analysis, since it increases the minimal χ 2 values. This can, for instance, lead to inadvertently over-fit the data by adding more parameters in order to reduce an artificially large χ 2 . This type of complications make it evident how critical it is to obtain a better knowledge of the unpolarized functions, not only for this but also for any other SIDIS asymmetries.
In conclusion, this type of analysis is an essential step which, after the first decade of pioneering studies, may lead us toward a new phase of high precision TMD physics. Our bottom-up approach to extract the Sivers function, while carefully keeping track of error estimation and of the sensitivity of the data to different TMD effects, may assume a relevant role as new measurements, with ever increasing statistics and precision, are becoming available from dedicated experiments.