Reweighting the Sivers function with jet data from STAR

The reweighting procedure that using Bayesian statistics incorporates the information contained in a new data set, without the need of re-fitting, is applied to the quark Sivers function extracted from Semi-Inclusive Deep Inelastic Scattering (SIDIS) data. We exploit the recently published single spin asymmetry data for the inclusive jet production in polarized $pp$ collisions from the STAR Collaboration at RHIC, which cover a much wider $x$ region compared to SIDIS measurements. The reweighting method is extended to the case of asymmetric errors and the results show a remarkable improvement of the knowledge of the quark Sivers function.


Introduction
Single spin asymmetries (SSAs) in inclusive and semiinclusive processes are an invaluable tool to deepen our knowledge of the internal structure of nucleons as well the hadronization mechanism. Despite the wealth and richness of available data and an extensive theoretical effort carried out in the last decades, the understanding of their origin still represents a formidable challenge from the phenomenological point of view. Related to this issue, there is nowadays a general consensus that the threedimensional (3D) picture of hadrons and the corresponding multi-parton correlations would lead to a better comprehension of the hadronic structure.
The 3D hadron structure in momentum space is encoded in a new class of parton distribution and fragmentation functions, the so-called Transverse Momentum Dependent distributions and fragmentation functions (TMDs), which depend on the collinear momentum fraction carried by the parton and its intrinsic transverse momentum. Among the eight leading-twist nucleon TMD distributions, the Sivers function [1,2] is one of the most studied and plays a seminal role. In fact, the Sivers distribution has several distinctive features: it is naively time-reversal odd, its existence is related to a nonvanishing parton orbital angular momentum and, more importantly, it is expected to be process dependent, having opposite signs in SIDIS and Drell-Yan (DY) processes [3,4].
For two-scale processes, such as SIDIS and DY, where Q 2 Q 1 ∼ Λ QCD , TMD factorization [5][6][7] is proven and SSAs are described in terms of convolutions of TMDs. These two scales are the virtuality of the exchanged boson and the transverse momentum of the observed hadron in SIDIS, and the invariant mass and transverse momentum of the lepton pair in DY. Therefore, such measurements are sensitive to the non perturbative transverse motion of bound partons in the nucleon encoded in TMDs. For processes with one characteristic perturbative scale, or Q 1 Q 2 Λ QCD , the collinear factorization at twist-3 level [8,9] plays the central role and SSAs are generated by the correlations of multi-parton densities in the nucleon, the so-called collinear twist-3 functions [8][9][10][11]. It was theoretically proven and demonstrated phenomenologically that in the intermediate region, Q 2 Q 1 Λ QCD , the two formalisms are related [12][13][14][15][16][17]. TMDs can be expressed in terms of collinear and twist-3 functions via Operator Product Expansion [18][19][20]. The integral relations [21][22][23] based on operator definitions also show the close relation between TMD and twist-3 distributions.
Inclusive jet production in pp collisions is one example of single-scale processes that can be described within the twist-3 approach, see for instance Ref. [24]. Nevertheless, other approaches, such as the so-called generalized parton model (GPM) [25][26][27], where a factorized formulation in terms of TMDs is assumed as the starting point, have been successfully applied in the analysis of SSAs for inclusive particle production in pp collisions, as well as in inclusive jet production [28]. The GPM is indeed a very powerful method to study processes where factorization is not established and can, hopefully, shed light on factorization breaking effects.
In the last decade, a color gauge invariant formulation of the GPM, named CGI-GPM, has been proposed [29] and then extensively developed in Refs. [30][31][32]. Its main feature is the inclusion of initial-(ISI) and final-(FSI) state interactions within a one-gluon exchange approximation. As a result, the Sivers function becomes non-universal, and its calculable process dependence can be absorbed into the partonic cross sections. Hence, in the evaluation of physical observables, one can still use the quark Sivers function obtained from SIDIS fits, but now convoluted with modified partonic cross sections, such that the expected sign change from SIDIS to DY is restored. Moreover, this modified GPM formalism has a very close connection [29] with the collinear twist-3 approach. In all these respects, it is extremely interesting to explore its potentiality and its implications.
In the spirit of testing the compatibility of the extraction of the Sivers function, as obtained by best-fitting the SIDIS azimuthal asymmetries, we analyze the recent SSA data for inclusive jet production in pp collisions from the STAR Collaboration at RHIC [33], within the GPM and the CGI-GPM approaches. These data cover a wide region of x, expanding the range explored in SIDIS measurements, and will allow us to improve and extend our current knowledge on the quark Sivers function.
As a global fit would at present require prohibitive machine power, here we employ an equivalent, but less numerically costly procedure, the so-called reweighting method [34][35][36][37][38][39] within Bayesian statistics. Reweighting allows us to properly include the information from the new set of data in the phenomenological analysis, estimate their impact on the extraction of the quark Sivers distributions, and determine via Bayesian statistics the fitted parameters and their errors. At the same time, this will also allow us to test the relevance of the expected process dependence of the Sivers function.
The paper is organized as follows. In Section 2 we recall the formalism applied to compute the SIDIS azimuthal Sivers asymmetry and the corresponding single spin asymmetries for inclusive jet production in pp collisions. In Section 3, we discuss the reweighting procedure and extend it to the case of asymmetric errors. In Section 4 we apply reweighting to obtain the statistical impact of the new STAR data sets to the extraction of the Sivers function. Finally, in Section 5, we draw our conclusions.

The formalism
Here we briefly recall the main aspects of the theoretical formalism employed in this study. All details can be found in the papers quoted below.
The expression for the azimuthal Sivers asymmetry in SIDIS, p ↑ → hX, is given by [40,41] where is the azimuthal modulation triggered by the correlation between the nucleon spin and the quark intrinsic transverse momentum. This effect is embodied in the Sivers function [42] which appears in the number density of unpolarized quarks, q, with intrinsic transverse momentum k ⊥ = k ⊥ (cos ϕ, sin ϕ), inside a transversely polarized proton p ↑ , with polarization vector S T = S T (cos φ S , sin φ S ) and mass M p , moving along the z direction. The dependence of structure functions and TMDs on the hard scale Q 2 is omitted for brevity. For the inclusive jet production in pp collisions, the SSA is defined as where dσ ↑(↓) denotes the single-polarized cross section, in which one of the protons is polarized along the transverse direction ↑(↓) with respect to the production plane. For this process, within the GPM as well as the CGI-GPM approach, only the Sivers effect, from quarks and gluons, can be at work. 1 The gluon contribution is negligible in the region of moderate and forward rapidities, as well as at small x F (x F = 2P jL / √ s, with P jL the longitudinal jet momentum) as shown, for both approaches, in Refs. [32,43]. The gluon Sivers effect can therefore be safely ignored in this study. Within the framework of the CGI-GPM, the numerator of the asymmetry is given by [29] where α s is the strong coupling constant, s is the centerof-mass (c.m.) energy, P j is the jet momentum, andŝ,t,û are the usual Mandelstam variables for the partonic subprocess ab → cd. Furthermore, f b/p (x b , k ⊥b ) is the TMD distribution for an unpolarized parton b inside the unpolarized proton. Notice that in a leading-order (LO) approach the jet is identified with the final parton c. Moreover, H Inc ab→cd are the perturbatively calculable hard scattering functions. In particular, the H Inc ab→cd functions where a is a quark or an antiquark can be found in Ref. [29]. The GPM results can be obtained from Eq. (4) by simply replacing H Inc ab→cd with the standard unpolarized partonic cross sections, H U ab→cd . Finally, the unpolarized cross section, dσ, at denomina- In the computation of the jet SSA we will adopt the jet transverse momentum P jT as the factorization scale.

The reweigthing method for TMDs
We now briefly illustrate the reweighting method that we will apply in our study. This is a well established technique in the context of collinear PDF extractions. Since the seminal work of Giele and Keller [34], the method has been used in both the Bayesian framework [35][36][37][38][39]44] and in the Hessian approach [45][46][47]. The reweighting procedure allows to assess the impact of new data sets on extractions of distributions describing hadron structure, avoiding a new global fit. It also indicates whether these additional data are consistent with the data sets used for the original extraction.
Let us consider a model for TMDs depending on a ndimensional set of parameters a = {a 1 , . . . , a n }. Traditionally, a χ 2 minimization procedure is employed in order to estimate the values of parameters that describe the experimental data. Let us suppose that a set of data y = y 1 , . . . , y N dat (with an associated covariance matrix C) is measured. Then, the χ 2 is defined as 2 where we have indicated by y i [a] the values computed using the theoretical model. The best-fit set a 0 , determined through the χ 2 minimization procedure, will have a corresponding minimum value χ 2 0 [a 0 ]. The uncertainty on the extracted TMDs can either be calculated in terms of Hessian eigensets or generated by applying Monte Carlo (MC) 2 If only uncorrelated uncertainties, σ i , are given, the new χ 2 re- procedures. In an ideal case, both methods yield the correct results and can be used in phenomenology. However, the former relies on the Gaussianity of the underlying distributions and does not necessarily account for the tails of distributions or for the potential presence of multiple solutions to the minimum of χ 2 in the parameter space.
In order to circumvent these complications, we will use a reweighting procedure based on the Bayesian inference. This will allow us to exploit the well known advantages of Bayesian inference, i.e. the ability to construct probability distributions for the parameters, and study the influence of the new data on the prior knowledge. Let us assume that some prior knowledge, say theoretical, exists on the parameters and that they are described by probability density functions π(a), i.e. the prior distributions. If no prior knowledge exists, π(a) are flat distributions.
We now generate the parameter sets, a k (k = 1, . . . , N set ), with corresponding χ 2 k ∈ [χ 2 0 , χ 2 0 + ∆χ 2 ], where ∆χ 2 is the desired tolerance corresponding to n parameters and at a certain confidence level (CL). For this purpose, we adopt a Markov Chain Monte Carlo procedure and produce N set parameter sets, employing a Metropolis-Hastings algorithm with an auto-regressive generating density [48]. Such approach allows for an efficient exploration and reconstruction of the parameter space, starting from the information contained in the error matrix obtained from the minimization procedure. The χ 2 k for every parameter set a k is calculated as: Using Bayes theorem, one can then calculate the posterior density given the data set as where L(y|a) is the likelihood, and Z = P(y) is the evidence, that ensures a normalized posterior density. P(a|y) will therefore incorporate the impact of the data on our knowledge of TMDs. Different choices for the likelihood have been discussed in the literature. Following Refs. [34,[36][37][38][39], here we adopt the likelihood definition as obtained by taking L(y|a) dy as the probability to find the new data confined in a differential volume dy around y. This results in defining the weights w k as follows: Note that Eqs. (10) and (11) are related to a symmetric error calculation. It is well known that when a parameter is symmetrically distributed (e.g. according to a Gaussian distribution) the mean value will coincide with the median value and the uncertainty at a certain confidence level, for instance at 68% (1σ) or 95% (2σ) will be symmetric. Still, this does not ensure that any distribution for any observable O(a k ), depending on Gaussian distributed parameters a k , is itself Gaussian.
In fact, asymmetric distributions may arise. Thus, mean and median for O(a k ) would, in general, be different: the more the distribution is asymmetric, the more the uncertainty on the observable is not properly described by a symmetrized error, as discussed for example in Ref. [49]. To overcome this potential issue, we extend the reweighting method, providing asymmetric uncertainties. Such errors are calculated using the rv discrete function of the SciPy Python library [50]. This function is able to build a discrete weighted distribution, providing automatically the mean value, the median and the asymmetric uncertainty interval at a specific confidence level. The interval is estimated considering equal areas around the median, and the endpoints of these areas are given automatically by rv discrete. For our analysis we will adopt the median as central value, and the uncertainty will be provided at a 2σ CL. We will perform this procedure with the data and parametrizations considered in Ref. [51] and explicitly demonstrate the equivalence of the reweighting procedure and the global fit.
The Bayesian inference is easily generalized to evaluate the impact of the new data y new , in our case the data from STAR Collaboration at RHIC [33]. As a matter of fact, the new evidence, i.e. the new data, will change the weights w k → w new k ≡ w k (χ 2 k + χ 2 new,k ) for each parameter set a k , where χ 2 new,k = χ 2 k [a k , y new ], Eq. (7). Therefore, the probability distribution for every parameter, as well as any other observable that depends on the TMDs, is expected to change. They can be calculated either with initial priors π(a) and weights w k (χ 2 k + χ 2 new,k ) or, equivalently, with posteriors P(a|y) after the global fit as the priors of reweighting and weights w k (χ 2 new,k ). Notice that, as expected, weights are multiplicative, w k (χ 2 k + χ 2 new,k ) ∝ w k (χ 2 k )w k (χ 2 new,k ), as χ 2 is additive. The resulting posteriors will contain the impact of the new data and the new values of the parameters, and the observables can be evaluated with Eqs. (10), (11) or with rv discrete function. As already mentioned, we will adopt the latter in our analysis.
In the following, we will apply the reweighting technique for the first time to a TMD function, and in particular to the quark Sivers functions.

Results
We will now present and discuss our results on the reweighting of the quark Sivers function. Before showing the outcome of the reweighting procedure using jet A N data from STAR, let us briefly illustrate the model parametrization chosen to describe the Sivers function.
Here we adopt the quark Sivers function extracted from SIDIS data in Ref. [51]. More precisely, we use the so called "reference fit", which results in a minimum χ 2 dof = 0.99 for N SIDIS dat = 220 data points (see Table 2 of Ref. [51]). Notice that this choice is motivated by the simplicity of this fit, which makes it particularly suitable for the purposes of this work. The relatively small number of parameters will allow us to highlight the effects of reweighting the Sivers function. Another important aspect is that this fit was performed paying special attention to the amount of information one can actually infer from data, reducing the assumptions which could bias the extraction. In this respect, we will be able to show and quantify the impact of a new set of data on our knowledge of the Sivers function using Bayesian inference.
The "reference fit" consists in a factorized x and k ⊥ dependence (the latter being Gaussian-like and flavor independent) for the up-and down-quark Sivers function. Within this parametrization, the Sivers function reads where q = u, d, and where ∆ N f (1) q/p ↑ (x) is the Sivers first k ⊥ -moment: Thus, the model depends on a total of five parameters: N u , N d , β u , β d , and k 2 ⊥ S . In the computation of any asymmetry, special care should be taken in the calculation of the corresponding unpolarized cross section. Here we follow the same approach adopted in Ref. [51], and compute the unpolarized SIDIS cross section, which appears at denominator in Eq. (1), by applying the k ⊥ -widths extracted in the multiplicity analysis of Ref. [52]. For the jet single spin asymmetry, we adopt the corresponding k ⊥ -width resulting from the fit of HERMES SIDIS data.
For the collinear parton densities, we use the CTEQ6L1 set of PDFs [53] and the DSS set of fragmentation functions [54]. Finally, for this fit, we generate N set = 2 · 10 5 parameter sets, adopting the corresponding ∆χ 2 at 2σ CL for N = 5 parameters, i.e. ∆χ 2 = 11.31, as tolerance.

Impact of jet data on the quark Sivers extraction
We now proceed by illustrating our final results. Very recently, the STAR Collaboration at RHIC has published new measurements of the single-spin asymmetries in inclusive jet production from polarized proton-proton collisions, at two different c.m. energies: √ s = 200 GeV and √ s = 500 GeV [33]. This new data set, amounting to N jet dat = 18 points, covers a wide range in x F (0.1 x F 0.6), and is useful to further constrain the quark Sivers function in the large-x region, where information from SIDIS data is either scarce or even absent. STAR measurements refer to electromagnetic jet production from polarized pp scattering. Two sets of data are collected at each energy: one fully inclusive and one imposing a cut on the photon multiplicity (n γ > 2). The latter is the most suitable for our analysis since it is not contaminated by single-photon and isolated π 0 production. Nevertheless, for completeness, we have also considered this data set and we will comment on the corresponding results.
In the following, apart from showing the impact of the A N jet data on the extraction of the Sivers function, we will also address whether any signal pointing towards a process dependence of the Sivers function itself can be observed.
As far as SIDIS is concerned, for the extraction of the Sivers function we will refer to Ref. [51], applying a LO approximation. Indeed, extractions of the Sivers function at higher perturbative orders exist, as in Refs. [55][56][57], but all the extracted Sivers functions are in good agreement, confirming the findings of Ref. [51] on the weak dependence of the asymmetries on the scale. We plan to perform similar TMD phenomenological analyses and reweighting of the Sivers functions to higher perturbative orders in the future.
In what follows, the predictions based on the Sivers functions as extracted by fitting only the SIDIS data are dubbed as "SIDIS", while the asymmetries computed after the reweighting procedure, by using the jet data as the new evidence as described in Section 3, are indicated by a "SIDIS+jet" label. The central value and the asymmetric uncertainty bands at 2σ CL are calculated using the procedure explained in Section 3, and adopting the corresponding weights for the "SIDIS" and the "SIDIS+jet" (GPM and CGI-GPM) cases. Figs. 1 and 2 show the results of the reweighting procedure for the Sivers function in the GPM and the CGI-GPM, respectively. On the upper (lower) panels, the comparison with data at √ s = 200 (500) GeV is shown.
The results in Figs. 1 and 2 clearly show the effect of the reweighting procedure. As we will see later, the jet SSA data measured by the STAR Collaboration allow to drastically reduce the uncertainty on the quark Sivers functions (especially for d quarks) and, in turn, on the corresponding estimates for the single-spin asymmetries in p ↑ p → jet X. It is important to emphasize the role of the STAR data, that extend up to x F 0.6 with remarkably small uncertainties, offering valuable information on the large x-region which in SIDIS remains largely uncovered. The reweighting procedure clearly indicates that jet data offer a crucial complementary information to SIDIS data. This analysis also points in favor of a good compatibility between the two data sets.
To better interpret the results shown in Figs. 1 and 2, in Fig. 3 we present the probability densities before and after reweighting for the five parameters of the "reference fit" and the χ 2 dof . The corresponding values and uncertainties, computed according to the procedure described at the end of Section 3, are gathered in Table 1. Notice that the results for the "SIDIS" case are fully compatible  Table 1: Summary of the results of the reweighting procedure for the fitted parameters. The quoted asymmetric uncertainties are at 2σ CL. with the findings of Ref. [51]. Some comments are in order: (i) The flavor independent Sivers Gaussian width (lower left panel in Fig. 3) does not vary significantly. This signals a mild role of TMD-evolution effects in the available jet SSA data. It is indeed worth to notice that at large x F we probe large P jT values (up to 4 ÷ 5 GeV depending on the data set). Since P jT represents the factorization scale adopted for this observable, we reach Q 2 scales significantly larger than those probed in SIDIS. However, asymmetries are ratios of cross sections where evolution and higher order effects tend to cancel out [58]. Although our parametrization does not have the complete features of TMD evolution, results of Refs. [51] are compatible with full TMD evolution at higher logarithmic accuracy [55][56][57]. (ii) The β q (q = u, d) parameters (mid panels in Fig. 3) that govern the large-x behavior of the Sivers function change, but in a different way when applying the GPM or the CGI-GPM formalisms: this is due to the fact that in the CGI-GPM approach color factors change the role of the u-and d-quark Sivers contributions. In particular, for the dominant channels in the forward rapidity region, like qg → qg, H inc in the CGI-GPM approach presents, roughly, a change of sign w.r.t. H U in the GPM (which are all positive). This implies that while the slightly positive A N at large x F in the GPM is driven by the positive sign of the up-quark Sivers function, in the CGI-GPM is given by the negative down-quark Sivers function. (iii) Concerning the normalization parameters (upper panels in Fig. 3) we see that while N u changes slightly, N d turns out to be smaller in size in the CGI-GPM approach. On the other hand, as mentioned above, the corresponding Sivers function for d-quarks is less suppressed in the large-x region. (iv) The χ 2 dof after the reweighting, calculated for N SIDIS+jet dat = 238 points, is different for the GPM and the CGI-GPM cases, slightly favoring the former approach.
To better visualize the effect of the reweighting procedure on the Sivers function, in Fig. 4 we show the comparison between the prior Sivers first moments and the corresponding moments after reweighting in the GPM (left panels) and in the CGI-GPM (right panels) frameworks. More specifically, in Fig. 4(a) we compare the first k ⊥ -moments before and after the reweighting, while in Fig. 4(b) we show the ratio of each Sivers first moment to its corresponding central value. As previously mentioned, the reweighting allows to constrain the Sivers function in the large-x region, not covered by the current SIDIS data. No significant difference is observed in the low-x region, as the model parametrization does not have any parameter controlling the low-x behavior. The reduction of the uncertainty is dramatic, especially for the d-quark Sivers function in the CGI-GPM approach, that SIDIS leaves largely unconstrained. This is confirmed by the much narrower reweighted probability density for the β d parameter (see mid-right panel of Fig. 3).
In order to quantify the uncertainty reduction, in Fig. 5 we show the ratio between the relative errors on the Sivers function before and after the reweighting procedure, both for the GPM (left panels) and CGI-GPM (right panels) formalisms. In the GPM approach, we see an uncertainty reduction of about 60% (80%) for the reweighted u(d)quark Sivers function at x > 0.2, while in the CGI-GPM case, and in the same kinematical region, we have ∼ 80% and ∼ 90% for the u and d Sivers function, respectively. Before concluding this section, let us comment on the use of the jet data set, where no cut on the photon multiplicity is imposed. The corresponding results of the reweighting are indeed very similar to those already shown in all respects, apart from the fact that, in this case, the resulting χ 2 dof slightly favors the CGI-GPM approach.

Conclusions and outlook
In the present analysis we applied, for the first time, a reweighting procedure to a TMD parton density, the quark Sivers distribution function extracted from SIDIS data. By exploiting the recently published single spin asymmetry data for inclusive jet production in polarized pp collisions from STAR [33], we showed the feasibility of such a procedure, which represents a valuable alternative to a full global fit. This allowed us, for the first time, to combine SIDIS azimuthal asymmetries data and the single spin asymmetries measured in p ↑ p → jet X processes in a global analysis. Moreover, by using two different approaches, the GPM and the CGI-GPM, we could also attempt to assess the degree of process dependence of the Sivers function, beside its sign change. Although the reweighted χ 2 dof slightly favors the GPM approach, further investigations are needed to have a clear discrimination between the two formalisms.
By including the jet SSA data from STAR we were able to significantly improve and extend our present knowledge of the quark Sivers function towards larger x values, not probed in current SIDIS data. In particular, their high precision allows to remarkably reduce the uncertainties of the Sivers function in such a kinematical region. We found that for the u-quark Sivers distribution, such reduction is about 60% in the GPM and 80% in the CGI-GPM frameworks, while for the d-quark case we observed a reduction of about 80% and 90% for GPM and CGI-GPM respectively.
This work has to be considered as an exploratory study to show, on one side, the potentiality of the reweighting procedure in the TMD sector and, on the other side, to refine the behavior of the Sivers function in a region not explored in SIDIS processes. The natural extension of this study will be a global analysis including also SSAs for inclusive pion production. This will allow us to simultaneously apply the reweighting procedure to the Collins fragmentation function, to transversity and to the Sivers function, as extracted by best-fitting the azimuthal asymmetries in SIDIS and e + e − processes. We also expect that forthcoming SIDIS measurements from COMPASS [59], JLab [60] and the future Electron Ion Collider [61,62] will play a crucial role in unravelling the nucleon structure in its full complexity. This rather ambitious program will provide important information on the impact of inclusive SSA data in the determination of these TMDs as well as a test of the compatibility of their extractions.