The three-dimensional distribution of quarks in momentum space

We present the distribution of unpolarized quarks in a transversely polarized proton in three-dimensional momentum space. Our results are based on consistent extractions of the unpolarized and Sivers transverse momentum dependent parton distributions (TMDs).

The antipode of taking a picture of a black hole is to take a picture of the inside of a proton, unveiling its internal constituents, confined in the most common element of the visible universe by the strong forces of Quantum Chromodynamics (QCD). Using data obtained from the scattering of a hard virtual photon off a proton, we map the density of quarks in three dimensions, i.e., as a function of their longitudinal momentum (along the photon's direction) and their transverse momentum (orthogonal to the photon). If the proton is unpolarized, the distribution is cylindrically symmetric: we determine it using recent results from our group [1]. If the proton is polarized in the transverse plane, the distributions of up and down quarks turn out to be distorted in opposite directions. This distortion, known as Sivers effect [2], is related to quark orbital angular momentum. We determine its details with the same formalism used for the unpolarized distribution. In this way, we obtain a consistent picture of the full 3-dimensional momentum distribution of quarks in a transversely polarized proton. Our study constitutes a benchmark for future determinations of multidimensional quark distributions, one of the main goals of existing and planned experimental facilities [3,4,5,6].
We consider a frame where the proton has momentum P with space component in the +ẑ direction, is polarized in the +ŷ direction, and is probed by a spacelike virtual photon with momentum q (with Q 2 = −q 2 ) in the −ẑ direction. We define the xy plane as transverse and we denote it with the subscript T . We consider the light-cone + direction (t +ẑ)/ √ 2 and we define it as longitudinal. If Q 2 is much larger than the proton's mass M 2 , the proton's momentum is approximately longitudinal (P + is the dominant component).
Our goal is to reconstruct the distribution of unpolarized quarks inside the nucleon as a function of three components of their momentum. In the frame we are considering, the distribution of a quark with flavor a in a transversely polarized nucleon N ↑ can be written in terms of two Transverse Momentum Distributions (TMDs) as [7] ρ a N ↑ (x, k x , k y ; where f a 1 is the unpolarized TMD and f ⊥a 1T is the Sivers TMD [2], k is the momentum of the quark, k T the modulus of its transverse component, and x = k + /P + is its longitudinal momentum fraction. Q 2 plays the role of a resolution scale. Recent extractions of f 1 have been published in Refs. [1,8,9,10]. Several parametrizations of f ⊥ 1T have been released up to now [11,12,13,14,15,16,17,18,19,20]. At variance with these works, in this paper we start from a recent determination of f 1 by our group [1] and we extract f ⊥ 1T using the same formalism, namely for the first time we reconstruct the 3-dimensional quark density of Eq. (1) in a fully consistent way within the TMD framework. Later publications have appeared [21,22,23] which adopt the same strategy; in the following, we will discuss a comparison with their results.
Both unpolarized and Sivers TMDs appear in the cross section of polarized Semi-Inclusive Deep-Inelastic Scattering (SIDIS) and vector-boson production processes. For SIDIS we consider the process (l)+N(P) → (l )+h(P h )+X, where a lepton with momentum l scatters off a nucleon target N with mass M and momentum P. In the final state, the scattered lepton with momentum l = l − q is detected, together with a hadron h with momentum P h and transverse momentum P hT . We define the usual SIDIS variables x Bj = Q 2 /(2P · q), y = P · q/(P · l), and z = P · P h /(P · q). In this study, we neglect power corrections of order M 2 /Q 2 and P 2 hT /Q 2 , which allow us also to identify x Bj = x. At leading twist and for a transversely polarized nucleon target N ↑ , the SIDIS cross section can be parametrized in terms of five structure functions [24]: where α is the fine structure constant, φ h and φ S indicate the azimuthal orientations of P hT and the target polarization S T in the transverse plane, respectively, the structure functions depend only on (x, z, P 2 hT , Q 2 ), and The structure function F UU,T can be obtained from the unpolarized cross section after integrating upon all azimuthal angles. The polarized structure function is experimentally measurable through the single spin asymmetry (SSA) Factorization theorems make it possible to write the structure functions at small transverse momentum (P 2 hT Q 2 ) in terms of TMDs and to derive their evolution equations. The latter ones are more involved than in the collinear framework because TMDs generally depend on two scales, µ 2 and ζ, that renormalize ultraviolet and rapidity divergences, respectively [25]. These two scales are usually chosen to be equal to the virtual photon mass: The unpolarized TMD f 1 enters the structure function F UU,T , while the Sivers TMD f ⊥ 1T enters the structure function F sin(φ h −φ S ) UT,T . Both structure functions can be defined as convolutions of TMDs upon quark transverse momenta [24], or as Fourier transforms of a product of functions in b T [26]. At leading order in the strong coupling α s (LO), they read where D a→h 1 is the Fourier-transformed expression of the corresponding TMD fragmentation function that describes how the parton a converts into a hadron h with transverse momentum P hT and carrying a fraction z of the parton energy. The Fourier transform of the unpolarized TMD is defined as where J l is the spherical Bessel functions of order l. Note that there is a factor 2π difference compared to the definition in the extraction of Ref. [1], denoted as Pavia17, which has been taken into account in the rest of the article. A similar definition holds for D a→h 1 . In Eq. (6), we have also introduced the first derivative of the Sivers function in Fourier space [26]: The limit of this formula for b T → 0 corresponds to the definition of the first k T -moment of the Sivers function: which is an x-dependent function and is related to the so-called Qiu-Sterman function [27,28]. The precise connection with the Qiu-Sterman function is nontrivial when considering higher-order corrections (see, e.g., [29,30,31]). However, these differences are relevant beyond the order we consider in our analysis.
The unpolarized TMD f 1 and the Sivers TMD f ⊥ 1T appear also in the process A ↑ (P A , S AT )+B(P B ) → γ * /W ± /Z 0 +X, where a hadron A with momentum P A and transverse polarization S AT scatters off an unpolarized hadron B with momentum P B , producing a vector boson with four-momentum q and rapidity y = 1 2 log[(q 0 + q z )/(q 0 − q z )], where P A points towards theẑ direction [1].
At leading twist and for q T q, the cross section can be parametrized in terms of five structure functions [32]. The relevant terms for the Sivers effect can be expressed as [1,21] dσ where Q 2 = q 2 is the invariant mass of the final state, φ q and φ S indicate the azimuthal orientations of q T and S AT in the transverse plane, respectively, and for V = γ * , W ± , Z 0 we have where s = (P A + P B ) 2 , N c is the number of colors, G F is the Fermi weak coupling constant, and B W/Z R is the branching ratio for the decay of vector bosons W ± and Z 0 with mass M W and M Z , respectively [33].
Again, the structure function F 1 T U for the Sivers effect is measurable through the SSA where x A = e y Q/ √ s, x B = e −y Q/ √ s, and at LO the structure functions read where the symbol implies adding the contribution of the flavor sum with A ↔ B. For V = W ± , the |V W aa | 2 are the elements of the CKM matrix and a, a run over light quark and antiquark flavors corresponding to W ± production: W + → ud, us, cd, cs , W − → du, dc, su, sc .
In this work, we take the unpolarized functions f 1 and D 1 from the Pavia17 extraction [1]. We extract the Sivers function using the very same approach: it is based on the TMD framework formulated in Ref. [25], which in turn elaborates on the original work of Collins, Soper, Sterman [34] (hence, in the following we refer to it as the CSS approach). The renormalization group evolution of TMDs is encoded in the so-called Sudakov form factor S , which contains the contribution of large logarithms. In this work, we perform the resummation of these logarithms at the next-to-leading-logarithmic (NLL) accuracy, as defined in detail in Ref. [10]. 1 The expression of S greatly simplifies if the starting scale of evolution is chosen as where γ E is the Euler constant. However, at large b T the TMD evolution runs into a nonperturbative region and becomes unreliable. In the CSS approach, this pathology is cured by the so-called b * -prescription, which amounts to replacing where b * is an arbitrary function of b T with appropriate asymptotic conditions [25]. In accordance with the extraction of the unpolarized TMD [1], in this analysis we adopt the following function With this choice, at large b T the function b * (b T ) saturates to b max , as already suggested by the CSS approach, and the scale µ b freezes at 1 GeV. In this way, the perturbative contributions to the TMD smoothly merge into the nonperturbative region, described by a parametric function (see below). At small b T (large k T ), the TMD formalism is not valid and must match onto the fixed-order formalism. The way the matching is implemented is not unique and the TMD contribution can be arbitrarily modified in this region. At variance with the standard CSS approach, in Eq. (17) we modify the high-transverse-momentum behavior of TMDs as −→ b min , which implies µ b → Q and preserves a meaningful definition of the integrals inside the Sudakov form factor S [1]. The latter prescription partially corresponds to modifying the resummed logarithms as in Ref. [35] (and similarly in Refs. [36,37]).
At NLL accuracy, the TMD evolution of the Sivers function from a starting scale Q 2 0 to a generic scale Q 2 is formally very similar to the unpolarized TMD f 1 [1]: The g K (b T ) is the above mentioned universal parametric function that describes the nonperturbative evolution. Together with the perturbative Sudakov form factor S , it is the same function that drives the evolution of the unpolarized TMD f 1 and is taken from the Pavia17 fit [1]. Without this information, it would not be possible to reliably calculate the Sivers function at the experimental scales. At the initial scale Q = Q 0 = 1 GeV, we have b max ≡ b min and µ b = Q 0 : the exponentials reduce to unity and evolution effects are switched off. The Sivers function has an intrinsic nonperturbative b T -distribution given by the function f ⊥ 1T NP , which needs to be determined from experimental data. For perturbative small values of b T , it can be shown that consistently f ⊥ 1T NP (x, b 2 T ) → 1 [10]. Hence, in the limit b T → 0 from Eq. (19) we recover Eq. 1T .
1T , we apply the same evolution as the collinear parton density f 1 using the HOPPET code [38]. This is an approximation of the full evolution [39,40,41,42,43,44,21]. In order to estimate the impact of the collinear evolution, we compared predictions obtained with our assumptions and predictions with no evolution. We found no significant difference in SIDIS kinematics (containing almost all data used in our fit) because of the limited range in Q 2 being spanned. For vector boson production, the difference becomes more relevant, but in both cases the theoretical predictions are small compared to the data, which are few and affected by large errors. In conclusion, the approximation in the implementation of collinear evolution does not affect the results of our present fit. The situation will certainly change when more and more precise data at high Q will be available. 1 At this accuracy, in the general formula of the Operator Product Expansion the hard functions and the matching coefficients can be neglected. 4 The Sivers function must satisfy the positivity bounds [45] for any value of x and k T . This constraint is essential to guarantee that the quark density distribution of Eq. (1) is positive everywhere. Therefore, it is convenient to parametrize the nonperturbative function f ⊥ 1T NP of Eq. (19) in momentum space. At the initial scale (Q 0 = 1 GeV), we write the Sivers function as The nonperturbative term f ⊥ 1T NP is given by where f 1NP is the corresponding nonperturbative term of the unpolarized TMD f 1 , and is consistently taken from the Pavia17 extraction [1]. More details about the explicit form of the involved functions can be found in Appendix A.
The M 1 , λ S are free parameters, and K(x) is a normalization factor to guarantee that the weighted integral of f ⊥ is parametrized as where T n (x) are Chebyshev polynomials of order n. The unpolarized collinear parton densities f a 1 are taken from the GJR parametrization [47], consistently with the Pavia17 fit. The flavor-dependent factor G a max , and the constraint |N a Siv | ≤ 1, are introduced to guarantee that the Sivers function of Eq. (21) satisfies the positivity condition of Eq. (20) (see Appendix A for more details). The free parameters N Siv , α, β, A, B are different for up, down, and sea quarks.
The actual total number of free parameters is 17. We fix them by fitting experimental data for the single transversespin asymmetries A sin(φ h −φ S ) UT of Eq. (4) for SIDIS measurements, and A V N of Eq. (12) for vector-boson-production measurements. In our fit, we include SIDIS measurements by the HERMES [48], COMPASS [49,50] and JLab collaborations [51], and W/Z-production measurements taken by the STAR collaboration [52]. Usually, the SIDIS asymmetries are presented as projections of the same dataset in x, z, and P hT . To avoid fully correlated measurements, we fit only the x projection because it has a direct impact on the x-dependence of the collinear function f ⊥(1) 1T . Similarly, for the STAR dataset we include only one of the projections of the measurements, specifically the data projected in rapidity. We select data by applying the same criteria used in the Pavia17 fit for unpolarized TMD, i.e., Q 2 > 1.4 GeV 2 , 0.20 < z < 0.74 and P hT < min[0.2Q, 0.7Qz] + 0.5 GeV [1]. With these kinematic cuts, we have a total of 125 data points: 30 from HERMES, 82 from COMPASS (32 from the 2009 analysis, and 50 from the 2017 analysis), 6 from JLab, and 7 from STAR.
Similarly to our previous Pavia17 extraction and to other studies of parton densities [53,54,55], we perform the fit using the bootstrap method. The method consists in creating M different replicas of the n original data by randomly shifting them with a Gaussian noise with the same variance as the experimental measurement. Each replica represents the possible outcome of an independent measurement. The number M is fixed by accurately reproducing the mean and standard deviation of the original data points. In our case, it turns out M = 200, which is also consistent with our Pavia17 fit [1]. We denote the replicated measurements as A Siv r , with the r index running from 1 to M, and with A Siv th {p r } the outcome of the calculated asymmetry using our functional form with the set of parameters {p r }. Once  replicas are generated, a minimization procedure is applied to each replica separately to search for the parameter values, {p r0 }, that minimize the error function where the covariance matrix is constructed as ∆A Siv corr i ∆A Siv corr j .
For each pair of experimental points (i, j), the covariance matrix contains the contributions of the statistical ∆A Siv stat and uncorrelated systematic ∆A Siv sys,unc experimental errors, the theoretical error ∆A Siv th due to the uncertainty in the unpolarized TMDs, as well as the correlated experimental uncertainties ∆A Siv corr (like, for example, a 7.3% target polarization correlated uncertainty for the HERMES data). Following the procedure outlined in Ref. [10], we apply the iterative t 0 -prescription [56] in order to avoid the D'Agostini bias that would lead to underestimate the predictions. The initial parameter values are chosen randomly within reasonable intervals. For each replica, the goodness of the fit is evaluated using the usual χ 2 test, which corresponds to the error function of Eq. (24), but with the original experimental data instead of the replicated ones. The maximal information about our results is given by the full ensemble of 200 replicas, combined with the corresponding unpolarized TMD replicas. To report our results in a concise way, we adopt the following choice: for any result (χ 2 values, parameter values, resulting distribution functions) we quote intervals containing 68% of the replicas, obtained by excluding the upper 16% and lower 16% values. These intervals correspond to the 1σ confidence level only if the observable's values follow a Gaussian distribution, which is not true in general. When it is not possible to draw uncertainty bands, we report the results obtained using replica 105, which was selected as a representative replica, since its parameters are the closest to the average ones both in the unpolarized and polarized case.
In Tab. 1 we give the value of the parameters obtained from our fit. For each one, we quote the central 68% of the 200 replica values (by quoting the average ± the semi-difference of the upper and lower limits). Parameters of replica 105, used for the multidimensional plots, are also given.
We obtain an excellent agreement between the experimental measurements and our theoretical prediction, with an overall value of χ 2 /d.o.f. = 1.12 ± 0.04 (total χ 2 = 121 ± 5). In Appendix B, we collected all figures that show the quality of our fit. Our parametrization is able to describe very well the COMPASS 2009 data set [49] [48] is somewhat worse (30 points with χ 2 = 46.2 ± 3.8; see Fig. B.7). We checked that the largest contribution to the χ 2 comes from the subset of data with K − in the final state [57]. Looking at the previous figures it is important Hatched bands from PV11 [15], EIKV [17], TC18 [18], JAM20 [20] parametrizations, and at different Q 2 as indicated in the figure. to notice, as a check of the results validity, that our predictions well describe also the z and P hT distributions, even if those projections of the data were not included in the fit (see Appendix B for more details).
The agreement with vector-boson-production STAR measurements [52] is worse than the SIDIS case, with a χ 2 = 13.97±0.6 for 7 points. However, the lower number of points (see Fig. B.8) indicates that STAR data have less influence on the global fit than the SIDIS data. In any case, we observe that our predictions follow the sign of the measurements, being negative for W + and positive for W − and Z 0 . The agreement is similar for the data points projected in p T not included in the fit (see Appendix B for more details).
In Fig. 1, we show the first transverse moment x f ⊥(1) 1T (Eq. (9), multiplied by x) as a function of x at Q 0 = 2 GeV for the up (left panel) and down quark (right panel). We compare our results (solid band) with other parametrizations available in the literature [15,17,18,20] (hatched bands, as indicated in the figure). In agreement with previous studies, the distribution for the up quark is negative, while for the down quark is positive and both have a similar magnitude. The Sivers function for sea quarks is very small and compatible with zero.
The authors of Ref. [21] also find results very similar to the ones in Fig. 1 when they fit the same SIDIS data and COMPASS Drell-Yan data with pion beams [58]. In this case, they also compute predictions for W ± and Z 0 production at STAR kinematics which are very close to our fitted bands displayed in Fig. B.8. Their strategy is very similar to the one adopted in this work but at higher perturbative accuracy, although their unpolarized TMDs are not obtained from an actual fit. However, when they include the STAR data in the global fit they artificially increase the statistical weight of those data by a factor ∼ 13. Their global χ 2 largely deteriorates and the uncertainty on the Sivers function significantly increases. Our finding is that because of large experimental errors STAR data does not affect much our final results when including them in the global fit, as discussed in detail in Appendix B.
The authors of Ref. [23] also perform a consistent extraction of both unpolarized and Sivers TMDs, and build contour plots of the density distribution in Eq. (1) similar to Fig. 2. A direct comparison is more difficult because the evolution of TMDs is achieved in a different framework, and the classification of the perturbative accuracy does not match the standard described in Ref. [10]. The displayed x-dependence of their Qiu-Sterman function (or related first k T -moment of the Sivers function as in Eq. (9)) is roughly similar, at least for up and down quarks. However, the sea-quark channel shows large oscillations at large x, which entail a strong breaking of the positivity constraint of Eq. (20).
In general, the result of a fit is biased whenever a specific fitting functional form is chosen at the initial scale. In our case, we tried to reduce this bias by adopting a flexible functional form, as it is evident particularly in Eq. (23). Nevertheless, we stress that our extraction is still affected by this bias and extrapolations outside the range where data exist (0.01 x 0.3) should be taken with due care. At variance with previous studies, in the denominator of the asymmetries in Eqs. (4) and (12) we are using unpolarized TMDs that were extracted from data in our previous Pavia17 fit, with their own uncertainties. Therefore, our uncertainty bands in Fig. 1 represent a realistic estimate of the statistical error of the Sivers function.  Figure 2: The density distribution ρ a p ↑ of an unpolarized quark with flavor a in a proton polarized along the +y direction and moving towards the reader, as a function of (k x , k y ) at Q 2 = 4 GeV 2 . Left panels for the up quark, right panels for the down quark. Upper panels for results at x = 0.1, lower panels at x = 0.01. For each panel, lower ancillary plots represent the 68% uncertainty band of the distribution at k y = 0 (where the effect of the distortion due to the Sivers function is maximal) while left ancillary plots at k x = 0 (where the distribution is the same as for an unpolarized proton). Results in the contour plots and the solid lines in the projections correspond to replica 105 (see text).
In Fig. 2, we show the density distribution ρ a p ↑ of an unpolarized quark a in a transversely polarized proton defined in Eq. (1), at x = 0.1 (upper panels) and x = 0.01 (lower panels) and at the scale Q 2 = 4 GeV 2 . The proton is moving towards the reader and is polarized along the +y direction. Since the up Sivers function is negative, the induced distortion is positive along the +x direction for the up quark (left panels), and opposite for the down quark (right panels). At x = 0.1 the distortion due to the Sivers effect is evident, since we are close to the maximum value of the function shown in Fig. 1. The distortion is more pronounced for down quarks, because the Sivers function is larger and at the same time the unpolarized TMD is smaller. The peak positions are approximately (k x ) max ≈ 0.1 GeV for up quarks and −0.15 GeV for down quarks. At lower values of x, the distortion disappears. These plots suggest that a virtual photon hitting a transversely polarized proton effectively "sees" more up quarks to its right and more down quarks to its left in momentum space.

8
The existence of this distortion requires two ingredients. First of all, the wavefunction describing quarks inside the proton must have a component with nonvanishing angular momentum. Secondly, effects due to final state interactions should be present [59], which in Feynman gauge can be described as the exchange of Coulomb gluons between the quark and the rest of the proton [60]. In simplified models [61], it is possible to separate these two ingredients and obtain an estimate of the angular momentum carried by each quark [62]. It turns out that up quarks give almost 50% contribution to the proton's spin, while all other quarks and antiquarks give less than 10% [15]. We will leave this model-dependent study to a future publication. A model-independent estimate of quark angular momentum requires the determination of parton distributions that depend simultaneously on momentum and position [63,64]. Nevertheless, the study of TMDs, and of the Sivers function in particular, can provide important constraints on models of the nucleon [65] and test lattice QCD computations [66].
In the near future, more data are expected from experiments at Jefferson Laboratory and CERN. Pioneering measurements in Drell-Yan processes with pion beams have been already reported [58], but they are not included in the present analysis because we do not have yet a consistent description of quark unpolarized TMDs in a pion. In the longer term, the recently approved Electron Ion Collider project [3,4] will provide a large amount of data in different kinematic regions compared to present experiments [67]. With this abundance of data, we will be able to reduce the error bands, extend the range of validity of the extractions to lower and higher values of x, and obtain a much more detailed knowledge of the 3-dimensional distribution of partons in momentum space.
The functional form we chose to parametrize the Sivers function is built in order to automatically satisfy the positivity bound It is given by 3) The f 1NP is the nonperturbative part of the unpolarized TMD f 1 and is taken from the Pavia17 extraction [1]: The distributions of values for the parameters N 1 , α, σ, λ are obtained from the Pavia17 fit and can be found in the NangaParbat repository 3 . For convenience, we report here the values for the relevant replica 105 In Eq. (A.3), the M 1 , λ S are free parameters. The K(x) is a normalization factor to guarantee that the weighted integral of f ⊥ 1T NP is 1 and the proper definition of the first k T -moment of the Sivers function is recovered in Eq. (A.2). It is given by The Fourier transform of f ⊥ 1T NP in Eq. (A. 3) reads The first transverse moment is parametrized as where T n (x) are Chebyshev polynomials of order n, and the unpolarized collinear parton densities f a 1 are taken from the GJR parametrization [47], consistently with the Pavia17 fit. The flavor-dependent factor G a max is defined as Together with the constraint |N a Siv | ≤ 1, it is introduced into Eq. (A.9) to guarantee that the Sivers function of Eq. (A.2) satisfies the positivity condition of Eq. (A.1). For the relevant replica 105 we have

Appendix B. Comparison with data
In this section, we present figures showing the quality of our fit. In all plots, our results are represented by solid bands whenever data points are actually included in the fit, while hatched bands are predictions for the same measurements projected over different kinematic variables. The predictions are obtained by integrating upon the variables over which they are not projected.
In HERMES Figure B.7: HERMES Sivers asymmetries from SIDIS off a proton target (H) with production of π + , π 0 , π − , K + , K − in the final state [48], presented as a function of x, z, P hT . Same notation as in Fig. B [52], presented as function of rapidity y and p W T . Same notation as in Fig. B.3.