Extraction of partonic transverse momentum distributions from semi-inclusive deep-inelastic scattering, Drell-Yan and Z-boson production

We present an extraction of unpolarized partonic transverse momentum distributions (TMDs) from a simultaneous fit of available data measured in semi-inclusive deep-inelastic scattering, Drell-Yan and Z-boson production. To connect data at different scales, we use TMD evolution at next-to-leading logarithmic accuracy. The analysis is restricted to the low-transverse-momentum region, with no matching to fixed-order calculations at high transverse momentum. We introduce specific choices to deal with TMD evolution at low scales, of the order of 1 GeV$^2$. This could be considered as a first attempt at a global fit of TMDs.


I. INTRODUCTION
Parton distribution functions describe the internal structure of the nucleon in terms of its elementary constituents (quarks and gluons). They cannot be easily computed from first principles, because they require the ability to carry out Quantum Chromodynamics (QCD) calculations in its nonperturbative regime. Many experimental observables in hard scattering experiments involving hadrons are related to parton distribution functions (PDFs) and fragmentation functions (FFs), in a way that is specified by factorization theorems (see, e.g., Refs. [1,2]). These theorems also elucidate the universality properties of PDFs and FFs (i.e., the fact that they are the same in different processes) and their evolution equations (i.e., how they get modified by the change in the hard scale of the process). Availability of measurements of different processes in different experiments makes it possible to test factorization theorems and extract PDFs and FFs through so-called global fits. On the other side, the knowledge of PDFs and FFs allows us to make predictions for other hard hadronic processes. These general statements apply equally well to standard collinear PDFs and FFs and to transverse-momentum-dependent parton distribution functions (TMD PDFs) and fragmentation functions (TMD FFs). Collinear PDFs describe the distribution of partons integrated over all components of partonic momentum except the one collinear to the parent hadron; hence, collinear PDFs are functions of the parton longitudinal momentum fraction x. TMD PDFs (or TMDs for short) include also the dependence on the transverse momentum k ⊥ . They can be interpreted as three-dimensional generalizations of collinear PDFs. Similar arguments apply to collinear FFs and TMD FFs [3].
There are several differences between collinear and TMD distributions. From the formal point of view, factorization theorems for the two types of functions are different, implying also different universality properties and evolution equations [4]. From the experimental point of view, observables related to TMDs require the measurement of some transverse momentum component much smaller than the hard scale of the process [5,6]. For instance, Deep-Inelastic Scattering (DIS) is characterized by a hard scale represented by the 4-momentum squared of the virtual photon (−Q 2 ). In inclusive DIS this is the only scale of the process, and only collinear PDFs and FFs can be accessed. In semi-inclusive DIS (SIDIS) also the transverse momentum of the outgoing hadron (P hT ) can be measured [7,8]. If P 2 hT Q 2 , TMD factorization can be applied and the process is sensitive to TMDs [2]. If polarization is taken into account, several TMDs can be introduced [7,[9][10][11][12]. Attempts to extract some of them have already been presented in the past [13][14][15][16][17][18][19][20][21]. In this work, we focus on the simplest ones, i.e., the unpolarized TMD PDF f q 1 (x, k 2 ⊥ ) and the unpolarized TMD FF D q→h 1 (z, P 2 ⊥ ), where z is the fractional energy carried by the detected hadron h, k ⊥ is the transverse momentum of the parton with respect to the parent hadron, and P ⊥ is the transverse momentum of the produced hadron with respect to the parent parton. Despite their simplicity, the phenomenology of these unpolarized TMDs presents several challenges [22]: the choice of a functional form for the nonperturbative components of TMDs, the inclusion of a possible dependence on partonic flavor [23], the implementation of TMD evolution [4,24], the matching to fixed-order calculations in collinear factorization [25].
We take into consideration three kinds of processes: SIDIS, Drell-Yan processes (DY) and the production of Z bosons. To date, they represent all possible processes where experimental information is available for unpolarized TMD extractions. The only important process currently missing is electron-positron annihilation, which is particularly important for the determination of TMD FFs [24]. This work can therefore be considered as the first attempt at a global fit of TMDs.
The paper is organized as follows. In Sec. II, the general formalism for TMDs in SIDIS, DY processes, and Z production is briefly outlined, including a description of the assumptions and approximations in the phenomenological implementation of TMD evolution equations. In Sec. III, the criteria for selecting the data analyzed in the fit are summarized and commented. In Sec. IV, the results of our global fit are presented and discussed. In Sec. V, we summarize the results and present an outlook for future improvements.

II. FORMALISM
A. Semi-inclusive DIS In one-particle SIDIS, a lepton with momentum l scatters off a hadron target N with mass M and momentum P . In the final state, the scattered lepton momentum l is measured together with one hadron h with mass M h and momentum P h . The corresponding reaction formula is (l) + N (P ) → (l ) + h(P h ) + X . (1) The space-like momentum transfer is q = l − l , with Q 2 = −q 2 . We introduce the usual invariants x = Q 2 2 P · q , y = P · q P · l , z = P · P h P · q , γ = 2M x Q . (2) The available data refer to SIDIS hadron multiplicities, namely to the differential number of hadrons produced per corresponding inclusive DIS event. In terms of cross sections, we define the multiplicities as where dσ h N is the differential cross section for the SIDIS process and dσ DIS is the corresponding inclusive one, and where P hT is the component of P h transverse to q (we follow here the notation suggested in Ref. [26]). In the single-photon-exchange approximation, the multiplicities can be written as ratios of structure functions (see Ref. [8] for details): where In the numerator of Eq. (4) the structure function F XY,Z corresponds to a lepton with polarization X scattering on a target with polarization Y by exchanging a virtual photon in a polarization state Z. In the denominator, only the photon polarization is explicitly written (T , L), as usually done in the literature. The semi-inclusive cross section can be expressed in a factorized form in terms of TMDs only in the kinematic limits M 2 Q 2 and P 2 hT Q 2 . In these limits, the structure function F U U,L of Eq. (4) can be neglected [27]. The structure function F L in the denominator contains contributions involving powers of the strong coupling constant α S at an order that goes beyond the level reached in this analysis; hence, it will be consistently neglected (for measurements and estimates of the F L structure function see, e.g., Refs. [28,29] and references therein).
To express the structure functions in terms of TMD PDFs and FFs, we rely on the factorized formula for SIDIS [2,[30][31][32][33][34][35][36][37] (see Fig. 1 for a graphical representation of the involved transverse momenta): Here, H U U,T is the hard scattering part; f a 1 (x, k 2 ⊥ ; Q 2 ) is the TMD PDF of unpolarized partons with flavor a in an unpolarized proton, carrying longitudinal momentum fraction x and transverse momentum k ⊥ . The D a h 1 (z, P 2 ⊥ ; Q 2 ) is the TMD FF describing the fragmentation of an unpolarized parton with flavor a into an unpolarized hadron h carrying longitudinal momentum fraction z and transverse momentum P ⊥ (see Fig. 1). TMDs generally depend on two energy scales [2], which enter via the renormalization of ultraviolet and rapidity divergencies. In this work we choose them to be equal and set them to Q 2 . The term Y U U,T is introduced to ensure a matching to the perturbative fixed-order calculations at higher transverse momenta.
In our analysis, we neglect any correction of the order of M 2 /Q 2 or higher to Eq. (6). At large Q 2 this is well justified. However, fixed-target DIS experiments typically collect a large amount of data at relatively low Q 2 values, where these assumptions should be all tested in future studies. The reliability of the theoretical description of SIDIS at low Q 2 has been recently discussed in Refs. [39,40].
Eq. (6) can be expanded in powers of α S . In the present analysis, we will consider only the terms at order α 0 S . In this case H a U U,T (Q 2 ) ≈ e 2 a and Y U U,T ≈ 0. However, perturbative corrections include large logarithms L ≡ log z 2 Q 2 /P 2 hT , so that α S L ≈ 1. In the present analysis, we will take into account all leading and Next-to-Leading Logarithms (NLL  [38]): a virtual photon (defining the reference axis) strikes a parton inside a proton. The parton has a transverse momentum k ⊥ (not measured). The struck parton fragments into a hadron, which acquires a further transverse momentum P ⊥ (not measured). The total measured transverse-momentum of the final hadron is P hT . When Q 2 is very large, the longitudinal components are all much larger than the transverse components. In this regime, P hT ≈ zk ⊥ + P ⊥ .
In these approximations (α 0 S and NLL), only the first term in Eq. (6) is relevant (often in the literature this has been called W term). We expect this term to provide a good description of the structure function only in the region where P 2 hT Q 2 . It can happen that Y U U,T , defined in the standard way (see, e.g., Ref. [31]), gives large contributions also in this region, but it is admissible to redefine it in order to avoid this problem [25]. We leave a detailed treatment of the matching to the high P 2 hT ≈ Q 2 region to future investigations. To the purpose of applying TMD evolution equations, we need to calculate the Fourier transform of the part of Eq. (6) involving TMDs. The structure function thus reduces to where we introduced the Fourier transforms of the TMD PDF and FF according tõ B. Drell-Yan and Z production In a Drell-Yan process, two hadrons A and B with momenta P A and P B collide at a center-of-mass energy squared s = (P A + P B ) 2 and produce a virtual photon or a Z boson plus hadrons. The boson decays into a lepton-antilepton pair. The reaction formula is The invariant mass of the virtual photon is Q 2 = q 2 with q = l + l . We introduce the rapidity of the virtual photon/Z boson where the z direction is defined along the momentum of hadron A (see Fig. 2). The cross section can be written in terms of structure functions [42,43]. For our purposes, we need the unpolarized cross section integrated over dΩ and over the azimuthal angle of the virtual photon, The elementary cross sections are where θ W is Weinberg's angle, M Z is the mass of the Z boson, and B R (Z → + − ) is the branching ratio for the Z boson decay in two leptons. We adopted the narrow-width approximation, i.e., we neglect contributions for Q 2 = M 2 Z . We used the values sin 2 θ W = 0.2313, M Z = 91.18 GeV, and B R (Z → + − ) = 3.366 [44]. Similarly to the SIDIS case, in the kinematic limit q 2 T Q 2 the structure function F 2 U U can be neglected (for measurement and estimates of this structure function see, e.g., Ref. [45] and references therein).
The longitudinal momentum fractions of the annihilating quarks can be written in terms of rapidity in the following way Some experiments use the variable x F , which is connected to the other variables by the following relations The structure function F 1 U U can be written as (see Fig. 2 for a graphical representation of the involved transverse momenta) As in the SIDIS case, in our analysis we neglect the Y U U term and we consider the hard coefficients only up to leading order in the couplings, i.e., The structure function can be conveniently expressed as a Fourier transform of the right-hand side of Eq. (16) as C. TMDs and their evolution Evolution equations quantitatively describe the connection between different values for the energy scales. In the following we will set their initial values to µ 2 b and their final values as Q 2 , so that only Q 2 and µ 2 b need to be specified in a TMD distribution. Following the formalism of Refs. [2,34], the unpolarized TMD distribution and fragmentation functions in configuration space for a parton with flavor a at a certain scale Q 2 can be written as We choose the scale µ b to be where γ E is the Euler constant and This variable replaces the simple dependence upon ξ T in the perturbative parts of the TMD definitions of Eqs. (20), (21). In fact, at large ξ T these parts are no longer reliable. Therefore, theξ * is chosen to saturate on the maximum value ξ max , as suggested by the CSS formalism [2,34]. On the other hand, at small ξ T the TMD formalism is not valid and should be matched to the fixed-order collinear calculations. The way the matching is implemented is not unique. In any case, the TMD contribution can be arbitrarily modified at small ξ T . In our approach, we choose to saturateξ * at the minimum value ξ min ∝ 1/Q. With the appropriate choices, for ξ T = 0 the Sudakov exponent vanishes, as it should [46,47]. Our choice partially corresponds to modifying the resummed logarithms as in Ref. [48] and to other similar modifications proposed in the literature [25,49]. One advantage of these kind of prescriptions is that by integrating over the impact parameter ξ T , the collinear expression for the cross section in terms of collinear PDFs is recovered, at least at leading order [25]. We remind the reader that there are different schemes available to deal with the high-ξ T region, such as the the so-called "complex-ξ prescription" [50] or an extrapolation of the perturbative small-ξ T calculation to the large ξ T region based on dynamical power corrections [51]. The values of ξ max and ξ min could be regarded as arbitrary scales separating perturbative from nonperturbative regimes. We choose to fix them to the values The motivations are the following: • with the above choices, the scale µ b is constrained between 1 GeV and Q, so that the collinear PDFs are never computed at a scale lower than 1 GeV and the lower limit of the integrals contained in the definition of the perturbative Sudakov factor (see Eq. (30)) can never become larger than the upper limit; • at Q = Q 0 = 1 GeV, ξ max = ξ min and there are no evolution effects; the TMD is simply given by the corresponding collinear function multiplied by a nonperturbative contribution depending on k ⊥ (plus possible corrections of order α S from the Wilson coefficients).
At NLL accuracy, for our choice of scales K(ξ * , µ b ) = 0. Similarly, the C andĈ are perturbatively calculable Wilson coefficients for the TMD distribution and fragmentation functions, respectively. They are convoluted with the corresponding collinear functions according to In the present analysis, we consider only the leading-order term in the α S expansion for C andĈ , i.e., As a consequence of the choices we made, the expression for the evolved TMD functions reduces to The Sudakov exponent S can be written as where the functions A and B have a perturbative expansions of the form To NLL accuracy, we need the following terms [31, 52] We use the approximate analytic expression for α S at NLO with the Λ QCD = 340 MeV, 296 MeV, 214 MeV for three, four, five flavors, respectively, corresponding to a value of α S (M Z ) = 0.117. We fix the flavor thresholds at m c = 1.5 GeV and m b = 4.7 GeV. The integration of the Sudakov exponent in Eq. (30) can be done analytically (for the complete expressions see, e.g., Refs. [36,53,54]). Following Refs. [55][56][57], for the nonperturbative Sudakov factor we make the traditional choice with g 2 a free parameter. Recently, several alternative forms have been proposed [58,59]. Also, recent theoretical studies aimed at calculating this term using nonperturbative methods [60]. All these choices should be tested in future studies. In Ref. [61], a good agreement with data was achieved even without this term, but this is not possible when including data at low Q 2 .
In this analysis, for the collinear PDFs f a 1 we adopt the GJR08FFnloE set [62] through the LHAPDF library [63], and for the collinear fragmentation functions the DSS14 NLO set for pions [64] and the DSS07 NLO set for kaons [65]. 3 We will comment on the use of other PDF sets in Sec. IV C.
We parametrize the intrinsic nonperturbative parts of the TMDs in the following ways After performing the anti-Fourier transform, the f 1NP and D 1NP in momentum space correspond to The TMD PDF at the starting scale is therefore a normalized sum of a Gaussian with variance g 1 and the same Gaussian weighted by a factor λk 2 ⊥ . The TMD FF at the starting scale is a normalized sum of a Gaussian with variance g 3 and a second Gaussian with variance g 4 weighted by a factor λ F P 2 ⊥ /z 2 . The choice of this particular functional forms is motivated by model calculations: the weighted Gaussian in the TMD PDF could arise from the presence of components of the quark wave function with angular momentum L = 1 [67][68][69][70][71]. Similar features occur in models of fragmentation functions [38,67,72].
The Gaussian width of the TMD distributions may depend on the parton flavor a [23,38,73]. In the present analysis, however, we assume they are flavor independent. The justification for this choice is that most of the data we are considering are not sufficiently sensitive to flavor differences, leading to unclear results. We will devote attention to this issue in further studies.
Finally, we assume that the Gaussian width of the TMD depends on the fractional longitudinal momentum x according to where α, σ, and N 1 ≡ g 1 (x) withx = 0.1, are free parameters. Similarly, for fragmentation functions we have where β, γ, δ, and N 3,4 ≡ g 3,4 (ẑ) withẑ = 0.5 are free parameters. The average transverse momentum squared for the distributions in Eq. (36) and (37) can be computed analytically:

III. DATA ANALYSIS
The main goals of our work are to extract information about intrinsic transverse momenta, to study the evolution of TMD parton distributions and fragmentation functions over a large enough range of energy, and to test their universality among different processes. To achieve this we included measurements taken from SIDIS, Drell-Yan and Z boson production from different experimental collaborations at different energy scales. In this chapter we describe the data sets considered for each process and the applied kinematic cuts.
Tab. I refers to the data sets for SIDIS off proton target (Hermes experiment) and presents their kinematic ranges. The same holds for Tab. II, Tab. III, Tab. IV for SIDIS off deuteron (Hermes and Compass experiments), Drell-Yan events at low energy and Z boson production respectively. If not specified otherwise, the theoretical formulas are computed at the average values of the kinematic variables in each bin.

A. Semi-inclusive DIS data
The SIDIS data are taken from Hermes [74] and Compass [75] experiments. Both data sets have already been analyzed in previous works, e.g., Refs. [23,76], however they have never been fitted together, including also the contributions deriving from TMD evolution.
The application of the TMD formalism to SIDIS depends on the capability of identifying the current fragmentation region. This task has been recently discussed in Ref. [39], where the authors point out a possible overlap among different fragmentation regions when the hard scale Q is sufficiently low. In this paper we do not tackle this problem and we leave it to future studies. As described in Tabs. I and II, we identify the current fragmentation region operating a cut on z only, namely 0.2 < z < 0.74.
Another requirement for the applicability of TMD factorization is the presence of two separate scales in the process. In SIDIS, those are the Q 2 and P 2 hT , which should satisfy the condition P 2 hT Q 2 , or more precisely P 2 hT /z 2 Q 2 . We implement this condition by imposing P hT < min[0.2 Q, 0.7 Qz] + 0.5 GeV. With this choice, P 2 hT is always smaller than Q 2 /3, but in a few bins (at low Q 2 and z) P 2 hT /z 2 may become larger than Q 2 . The applicability of TMD factorization in this case could be questioned. However, as we will explain further in Sec. IV C, we can obtain a fit that can describe a wide region of P hT and can also perform very well in a restricted region, where TMD factorization certainly holds.
All these choices are summarized in Tabs. I and II.

Hermes data
Hermes hadron multiplicities are measured in a fixed target experiment, colliding a 27.6 GeV lepton beam on a hydrogen (p) or deuterium (D) gas target, for a total of 2688 points. These are grouped in bins of (x, z, Q 2 , P hT ) with the average values of (x, Q 2 ) ranging from about (0.04, 1.25 GeV 2 ) to (0.4, 9.2 GeV 2 ). The collinear energy fraction z in Eq. (2) ranges in 0.1 ≤ z ≤ 0.9. The transverse momentum of the detected hadron satisfies 0.1 GeV ≤ |P hT | ≤ 1.3 GeV. The peculiarity of Hermes SIDIS experiment lies in the ability of its detector to distinguish between pions and kaons in the final state, in addition to determining their momenta and charges. We consider eight different combinations of target (p, D) and detected charged hadron (π ± , K ± ). The Hermes collaboration published two distinct sets, characterized by the inclusion or subtraction of the vector meson contribution. In our work we considered only the data set where this contribution has been subtracted.

Compass data
The Compass collaboration extracted multiplicities for charge-separated but unidentified hadrons produced in SIDIS off a deuteron ( 6 LiD) target [75]. The number of data points is an order of magnitude higher compared to the Hermes experiment. The data are organized in multidimensional bins of (x, z, Q 2 , P hT ), they cover a range in (x, Q 2 ) from about (0.005, 1.11 GeV 2 ) to (0.09, 7.57 GeV 2 ) and the interval 0.2 ≤ z ≤ 0.8. The multiplicities published by Compass are affected by normalization errors (see the erratum to Ref. [75]). In order to avoid this issue, we divide the data in each bin in (x, z, Q 2 ) by the data point with the lowest P 2 hT in the bin. As a result, we define the normalized multiplicity as where the multiplicity m h N is defined in Eq. (3). When fitting normalized multiplicities, the first data point of each bin is considered as a fixed constraint and excluded from the degrees of freedom.

B. Low-energy Drell-Yan data
We analyze Drell-Yan events collected by fixed-target experiments at low-energy. These data sets have been considered also in previous works, e.g., in Ref. [56,57,77,78]. We used data sets from the E288 experiment [79], which measured the invariant dimuon cross section Ed 3 σ/dq 3 for the production of µ + µ − pairs from the collision of a proton beam with a fixed target, either composed of Cu or Pt. The measurements were performed using proton incident energies of 200, 300 and 400 GeV, producing three different data sets. Their respective center of mass energies are √ s = 19.4, 23.8, 27.4 GeV. We also included the set of measurements Ed 3 σ/dq 3 from E605 [80], extracted from the collision of a proton beam with an energy of 800 GeV ( √ s = 38.8 GeV) on a copper fixed target . The explored Q values are higher compared to the SIDIS case, as can be seen in Tab. III. E288 provides data at fixed rapidity, whereas E605 provides data at fixed x F = 0.1. We can apply TMD factorization if q 2 T Q 2 , where q T is the transverse momentum of the intermediate electroweak boson, reconstructed from the kinematics of the final state leptons. We choose q T < 0.2 Q + 0.5 GeV. As suggested in Ref. [79], we consider the target nuclei as an incoherent ensemble composed 40% by protons and 60% by neutrons.
As we already observed, results from E288 and E605 experiments are reported as Ed 3 σ d 3 q ; this variable is related to the differential cross section of Eq. (12) in the following way: where φ is the polar angle of q T and the third term is the average over φ. Therefore, the invariant dimuon cross section can be obtained from Eq. (12) integrating over Q 2 and adding a factor 1/π to the result Numerically we checked that integrating in Q 2 only the prefactor σ γ q (see Eq. (13)) introduces only a negligible error in the theoretical estimates. We also assume that α em does not change within the experimental bin. Therefore, for Drell-Yan we obtain where Q i,f are the lower and upper values in the experimental bin.

C. Z-boson production data
In order to reach higher Q and q T values, we also consider Z boson production in collider experiments at Tevatron. We analyze data from CDF and D0, collected during Tevatron Run I [81,82] at √ s = 1.8 TeV and Run II [83,84] at √ s = 1.96 TeV. CDF and D0 collaborations studied the differential cross section for the production of an e + e − pair from pp collision through an intermediate Z vector boson, namely pp → Z → e + e − + X.
The invariant mass distribution peaks at the Z-pole, Q ≈ M Z , while the transverse momentum of the exchanged Z ranges in 0 < q T < 20 GeV. We use the same kinematic cut applied to Drell-Yan events: q T < 0.2 Q + 0.5 GeV = 18.7 GeV, since Q is fixed to M Z .
The observable measured in CDF and D0 is apart from the case of D0 Run II, for which the published data refer to 1/σ × dσ/dq T . In order to work with the same observable, we multiply the D0 Run II data by the total cross section of the process σ exp = 255.8 ± 16 pb [85].
In this case, we add in quadrature the uncertainties of the total cross section and of the published data. We normalize our functional form with the factors listed in Tab. IV. These are the same normalization factors used in Ref. [78], computed by comparing the experimental total cross section with the theoretical results based on the code of Ref. [86]. These factors are not precisely consistent with our formulas. In fact, as we will discuss in Sec. IV C a 5% increase in these factors would improve the agreement with data, without affecting the TMD parameters.

D. The replica method
Our fit is based on the replica method. In this section we describe it and we give a definition of the χ 2 function minimized by the fit procedure. The fit and the error analysis are carried out using a similar Monte Carlo approach as in Refs. [23,87,88] and taking inspiration from the work of the Neural-Network PDF (NNPDF) collaboration (see, e.g., Refs. [89][90][91]). The approach consists in creating M replicas of the data points. In each replica (denoted x range 0.04 < x < 0.4     by the index r), each data point i is shifted by a Gaussian noise with the same variance as the measurement. Each replica, therefore, represents a possible outcome of an independent experimental measurement, which we denote by m h N,r (x, z, P 2 hT , Q 2 ). The number of replicas is chosen so that the mean and standard deviation of the set of replicas accurately reproduces the original data points. In this case 200 replicas are sufficient for the purpose. The error for each replica is taken to be equal to the error on the original data points. This is consistent with the fact that the variance of the M replicas should reproduce the variance of the original data points.

Hermes Hermes Hermes Hermes Compass Compass
A minimization procedure is applied to each replica separately, by minimizing the following error function: The sum runs over the i experimental points, including all species of targets N and final-state hadrons h. In each z bin for each replica the values of the collinear fragmentation functions D a h 1 are independently modified with a Gaussian noise with standard deviation equal to the theoretical error ∆D a h 1 . In this work we rely on different parametrizations for D a h 1 : for pions we use the DSEHS analysis [64] at NLO in α S ; for kaons we use the DSS parametrization [65] at LO in α S . The uncertainties ∆D a h 1 are estimated from the plots in Ref. [92]; they represents the only source of uncertainty in ∆m h N,theo . Statistical and systematic experimental uncertainties ∆m h N,stat and ∆m h N,sys are taken from the experimental collaborations. We do not take into account the covariance among different kinematic bins.
We minimize the error function in Eq. (46) with Minuit [93]. In each replica we randomize the starting point of the minimization, to better sample the space of fit parameters. The final outcome is a set of M different vectors of best-fit parameters, {p 0r }, r = 1, . . . M, with which we can calculate any observable, its mean, and its standard deviation. The distribution of these values needs not to be necessarily Gaussian. In fact, in this case the 1σ confidence interval is different from the 68% interval. The latter can simply be computed for each experimental point by rejecting the largest and the lowest 16% of the M values.
Although the minimization is performed on the function defined in Eq. (46), the agreement of the M replicas with the original data is expressed in terms of a χ 2 function defined as in Eq. (46) but with the replacement m h N,r → m h N , i.e., with respect to the original data set. If the model is able to give a good description of the data, the distribution of the M values of χ 2 /d.o.f. should be peaked around one.

IV. RESULTS
Our work aims at simultaneously fitting for the first time data sets related to different experiments. In the past, only fits related either to SIDIS or hadronic collisions have been presented. Here we mention a selection of recent existing analyses.
In Ref. [23], the authors fitted Hermes multiplicities only (taking into account a total of 1538 points) without taking into account QCD evolution. In that work, a flavor decomposition in transverse momentum of the unpolarized TMDs and an analysis of the kinematic dependence of the intrinsic average square transverse momenta were presented. In Ref. [76] the authors fitted Hermes and Compass multiplicities separately (576 and 6284 points respectively), without TMD evolution and introducing an ad-hoc normalization for Compass data. A fit of SIDIS data including TMD evolution was performed on measurements by the H1 collaboration of the so-called transverse energy flow [55,94].
Looking at data from hadronic collisions, Konychev and Nadolsky [57] fitted data of low-energy Drell-Yan events and Z-boson production at Tevatron, taking into account TMD evolution at NLL accuracy (this is the most recent of a series of important papers on the subject [56,77,95]). They fitted in total 98 points. Contrary to our approach, Konychev and Nadolsky studied the quality of the fit as a function of ξ max . They found that the best value for ξ max is 1.5 GeV −1 (to be compared to our choice ξ max ≈ 1.123 GeV −1 , see Sec. II C). Comparisons of best-fit values in the nonperturbative Sudakov form factors are delicate, since the functional form is different from ours. In 2014 D'Alesio, Echevarria, Melis, Scimemi performed a fit [78] of Drell-Yan data and Z-boson production data at Tevatron, focusing in particular on the role of the nonperturbative contribution to the kernel of TMD evolution. This is the fit with the highest accuracy in TMD evolution performed up to date (NNLL in the Sudakov exponent and O(α S ) in the Wilson coefficients). In the same year Echevarria, Idilbi, Kang and Vitev [15] presented a parametrization of the unpolarized TMD that described qualitatively well some bins of Hermes and Compass data, together with Drell-Yan and Z-production data. A similar result was presented by Sun, Isaacson, Yuan and Yuan [96].
In the following, we detail the results of a fit to the data sets described in Sec. III with a flavor-independent configuration for the transverse momentum dependence of unpolarized TMDs. In Tab. V we present the total χ 2 . The number of degrees of freedom (d.o.f.) is given by the number of data points analyzed reduced by the number of free parameters in the error function. The overall quality of the fit is good, with a global χ 2 /d.o.f. = 1.55 ± 0.05. Uncertainties are computed as the 68% confidence level (C.L.) from the replica methodology.

A. Agreement between data and theory
The partition of the global χ 2 among SIDIS off a proton, SIDIS off a deuteron, Drell-Yan and Z production events is given in Tab. VI, VII, VIII, IX respectively.

Semi-inclusive DIS
For SIDIS at Hermes off a proton, most of the contribution to the χ 2 comes from events with a π + in the final state. In Ref. [23] the high χ 2 was attributed to the poor agreement between experiment and theory at the level of the collinear multiplicities. In this work we use a newer parametrization of the collinear FFs (DSEHS [64]), based on a fit which includes Hermes collinear pion multiplicities. In spite of this improvement, the contribution to χ 2 from Hermes data is higher then in Ref. [23], because the present fit includes data from other experiments (Hermes represents less than 20% of the whole data set). The bins with the worst agreement are at low Q 2 . As we will discuss in Sec. IV C, we think that the main reason for the large χ 2 at Hermes is a normalization difference. This may also be due to the fact that we are computing our theoretical estimates at the average values of the kinematic variables, instead of integrating the multiplicities in each bin. Kaon multiplicities have in general a lower χ 2 , due to the bigger statistical errors and the large uncertainties for the kaon FFs.

Hermes
Hermes Hermes Hermes  For pion production off a deuteron at Hermes the χ 2 is lower with respect to the production off a proton, but still compatible within uncertainties. For kaon production off a deuteron the χ 2 is higher with respect to the scattering off a proton. The difference is especially large for K − . SIDIS at Compass involves scattering off deuteron only, D → h ± , and we identify h ≡ π. The quality of the agreement between theory and Compass data is better than in the case of pion production at Hermes. This depends on at least two factors: first, our fit is essentially driven by the Compass data, which represent about 75% of the whole data set; second, the observable that we fit in this case is the normalized multiplicity, defined in Eq. (41). This automatically eliminates most of the discrepancy between theory and data due to normalization.

Hermes
Hermes Hermes Hermes   Fig. 3 presents the agreement between the theoretical formula in (3) and the Hermes multiplicities for production of pions off a proton and a deuteron. Different x , z and Q 2 bins are displayed as a function of the transverse momentum of the detected hadron P hT . The grey bands are an envelope of the 200 replica of best-fit curves. For every point in P hT we apply a 68% C.L. selection criterion. Points marked with different symbols and colors correspond to different z values. There is a strong correlation between x and Q 2 that does not allow us to explore the x and Q 2 dependence of the TMDs separately. Studying the contributions to the χ 2 /points as a function of the kinematics, we notice that the χ 2 (Q 2 ) tends to improve as we move to higher Q 2 values, where the kinematic approximations of factorization are more reliable. Moreover, usually the χ 2 (z) increases at lower z values. Fig. 4 has same contents and notation as in Fig. 3 but for kaons in the final state. In this case, the trend of the agreement as a function of Q 2 is not as clear as for the case of pions: good agreement is found also at low Q 2 .
In Fig. 5 we present Compass normalized multiplicities (see Eq. (41)) for production of π − off a deuteron for different x , z , and Q 2 bins as a function of the transverse momentum of the detected hadron P hT . The open marker around the first P hT point in each panel indicates that the first value is fixed and not fitted. The correlation between x and Q 2 is less strong than at Hermes and this allows us to study different x bins at fixed Q 2 . For the highest Q 2 bins, the agreement is good for all x , z and P 2 hT . In bins at lower Q 2 , the descriptions gets worse, especially at low and high z. For fixed Q 2 and high z , a good agreement is recovered moving to higher x bins (see, e.g., the third line from the top in Fig. 5). Fig. 6 has same content and notation as in Fig. 5, but for h + ≡ π + . The same comments on the agreement between theory and the data apply.

Drell-Yan and Z production
The low energy Drell-Yan data collected by the E288 and E605 experiments at Fermilab have large error bands (see Fig. 7). This is why the χ 2 values in Tab. VIII are rather low compared to the other data sets.
The agreement is also good for Z boson production, see Tab. IX. The statistics from Run-II is higher, which generates smaller experimental uncertainties and higher χ 2 , especially for the CDF experiment.    Fig. 7 displays the cross section for DY events differential with respect to the transverse momentum q T of the virtual photon, its invariant mass Q 2 and rapidity y. As for the case of SIDIS, the grey bands are the 68% C.L. envelope of the 200 replicas of the fit function. The four panels represents different values for the rapidity y or x F (see Eq. (15)). In each panel, we have plots for different Q 2 values. The lower is Q, the less points in q T we fit (see also Sec. III B). The hard scale lies in the region 4.5 < Q < 13.5 GeV. This region is of particular importance, since these "moderate" Q values should be high enough to safely apply factorization and, at the same time, low enough in order for the nonperturbative effects to not be shaded by transverse momentum resummation.
In Fig. 8 we compare the cross section differential with respect to the transverse momentum q T of the virtual Z (namely Eq. (12) integrated over η) with data from CDF and D0 at Tevatron Run I and II. Due to the higher Q = M Z , the range explored in q T is much larger compared to all the other observables considered. The tails of the distributions deviate from a Gaussian behavior, as it is also evident in the bins at higher Q 2 in Fig. 7. The band from the replica methodology in this case is much narrower, due to the reduced sensitivity to the intrinsic transverse momenta at Q = M Z and to the limited range of best-fit values for the parameter g 2 , which controls soft-gluon emission. As an effect of TMD evolution, the peak shifts from ∼ 1 GeV for Drell-Yan events in Fig. 7 to ∼ 5 GeV in Fig. 8. The position of the peak is affected both by the perturbative and the nonperturbative part of the Sudakov exponent (see Sec. II C and [22]). Most of the contributions to the χ 2 comes from normalization effects and not from the shape in q T (see Sec. IV C).

FIG. 3:
Hermes multiplicities for production of pions off a proton and a deuteron for different x , z , and Q 2 bins as a function of the transverse momentum of the detected hadron P hT . For clarity, each z bin has been shifted by an offset indicated in the legend.

FIG. 4:
Hermes multiplicities for production of kaons off a proton and a deuteron for different x , z , and Q 2 bins as a function of the transverse momentum of the detected hadron P hT . For clarity, each z bin has been shifted by an offset indicated in the legend.
FIG. 5: Compass multiplicities for production of negative hadrons (π − ) off a deuteron for different x , z , and Q 2 bins as a function of the transverse momentum of the detected hadron P hT . Multiplicities are normalized to the first bin in P hT for each z value (see (41)). For clarity, each z bin has been shifted by an offset indicated in the legend.

B. Transverse momentum dependence at 1 GeV
The variables ξ min and ξ max delimit the range in ξ T where transverse momentum resummation is computed perturbatively. The g 2 parameter enters the nonperturbative Sudakov exponent and quantifies the amount of transverse FIG. 6: Compass multiplicities for production of positive hadrons (π + ) off a deuteron for different x , z , and Q 2 bins as a function of the transverse momentum of the detected hadron P hT . Multiplicities are normalized to the first bin in P hT for each z value (see (41)). For clarity, each z bin has been shifted by an offset indicated in the legend.
momentum due to soft gluon radiation that is not included in the perturbative part of the Sudakov form factor. As already explained in Sec. II C, in this work we fix the value for ξ min and ξ max in such a way that at Q = 1 GeV the unpolarized TMDs coincide with their nonperturbative input. We leave g 2 as a fit parameter.
Tab. X summarizes the chosen values of ξ min , ξ max and the best-fit value for g 2 . The latter is given as an average with 68% C.L. uncertainty computed over the set of 200 replicas. We also quote the results obtained from replica 105, since its parameters are very close to the mean values of all replicas. We obtain a value g 2 = 0.13 ± 0.01, smaller than the value (g 2 = 0.184 ± 0.018) obtained in Ref. [57], where however no SIDIS data was taken into consideration, and smaller than the value (g 2 = 0.16) chosen in Ref. [15]. We stress however that our prescriptions involving both ξ min and ξ max are different from previous works.
Tab. XI collects the best-fit values of parameters in the nonperturbative part of the TMDs at Q = 1 GeV (see Eqs. (34) and (35)); as for g 2 , we give the average value over the full set of replicas and the standard deviation based on a 68% C.L. (see Sec. III D), and we also quote the value of replica 105.
In Fig. 9 we compare different extractions of partonic transverse momenta. The horizontal axis shows the value of the average transverse momentum squared for the incoming parton, k 2 ⊥ (x = 0.1) (see Eq. (40)). The vertical axis shows the value of P 2 ⊥ (z = 0.5), the average transverse momentum squared acquired during the fragmentation process (see Eq. (40)). The white square (label 1) indicates the average values of the two quantities obtained in the present analysis at Q 2 = 1 GeV 2 . Each black dot around the white square is an outcome of one replica. The red region around the white square contains the 68% of the replicas that are closest to the average value. The same applies to the white circle and the orange region around it (label 2), related to the flavor-independent version of the analysis in Ref. [23], obtained by fitting only Hermes SIDIS data at an average Q 2 = 2.4 GeV 2 and neglecting QCD evolution. A strong anticorrelation between the transverse momenta is evident in this older analysis. In our new analysis, the inclusion of Drell-Yan and Z production data adds physical information about TMD PDFs, free from the influence of TMD FFs. This reduces significantly the correlation between k 2 ⊥ (x = 0.1) and P 2 ⊥ (z = 0.5). The 68% confidence region is smaller than in the older analysis. The average values of k 2 ⊥ (x = 0.1) are similar and compatible within error bands. The values of P 2 ⊥ (z = 0.5) in the present analysis turn out to be larger than in the older analysis, an effect that is due mainly to Compass data. It must be kept in mind that the two analyses lead also to differences in the x and z dependence of the transverse momentum squared. This dependence is shown in Fig. 10 (a) for k 2 ⊥ (x) and Fig. 10 1): average values (white square) obtained in the present analysis, values obtained from each replica (black dots) and 68% C.L. area (red); (2) results from Ref. [23], (3) results from Ref. [97], (4) results from Ref. [76] for Hermes data, (5) results from Ref. [76] for Hermes data at high z, (6) results from Ref. [76] for normalized Compass data, (7) results from Ref. [76] for normalized Compass data at high z, (8) results from Ref. [15].

C. Stability of our results
In this subsection we discuss the effect of modifying some of the choices we made in our default fit. Instead of repeating the fitting procedure with different choices, we limit ourselves to checking how the χ 2 of a single replica is affected by the modifications.
As starting point we choose replica 105, which, as discussed above, is one of the most representative among the whole replica set. The global χ 2 /d.o.f. of replica 105 is 1.51. We keep all parameters fixed, without performing any new minimization, and we compute the χ 2 /d.o.f. after the modifications described in the following.
First of all, we analyze Hermes data with the same strategy as Compass, i.e., we normalize Hermes data to the value of the first bin in P hT . In this case, the global χ 2 /d.o.f. reduces sharply to 1.27. The partial χ 2 for the different SIDIS processes measured at Hermes are shown in Table XII. This confirms that normalization effects are the main contribution to the χ 2 of SIDIS data and have minor effects on TMD-related parameters. In fact, even if we perform a new fit with this modification, the χ 2 does not improve significantly and parameters do not change much.  We consider the effect of changing the normalization of the Z-boson data: if we increase the normalization factors quoted in the last row of Tab. IV by 5%, the χ 2 quoted in the last row of Tab. IX drops to 0.66, 0.52, 0.65, 0.68. This effect is also already visible by eye in Fig. 8: the theoretical curves are systematically below the experimental data points, but the shape is reproduced very well.
We consider the sensitivity of our results to the parameterizations adopted for the collinear quark PDFs. The χ 2 /d.o.f. varies from its original value 1.51, obtained with the NLO GJR 2008 parametrization [62], to 1.84 using NLO MSTW 2008 [98], and 1.85 using NLO CJ12 [99]. In both cases, the agreement with Hermes and Z boson data is not affected significanlty, the agreement with Compass data becomes slightly worse, and the agreement with DY data becomes clearly worse.
An extremely important point is the choice of kinematic cuts. Our default choices are listed in Tabs. I-IV. We consider also more stringent kinematic cuts on SIDIS data: Q 2 > 1.5 GeV 2 and 0.25 < z < 0.6 instead of Q 2 > 1.4 GeV 2 and 0.2 < z < 0.7, leaving the other ones unchanged. The number of bins with these cuts reduces from 8059 to 5679 and the χ 2 /d.o.f. decreases to the value 1.23. In addition, if we replace the constraint P hT < Min[0.2 Q, 0.7 Qz] + 0.5 GeV with P hT < Min[0.2 Q, 0.5 Qz] + 0.3 GeV, the number of bins reduces to 3380 and the χ 2 /d.o.f. decreases further to 0.96. By adopting the even stricter cut P hT < 0.2 Qz, the number of bins drops to only 477, with a χ 2 /d.o.f. =1.02. We can conclude that our fit, obtained by fitting data in an extended kinematic region, where TMD factorization may be questioned, works extremely well also in a narrower region, where TMD factorization is expected to be under control.

V. CONCLUSIONS
In this work we demonstrated for the first time that it is possible to perform a simultaneous fit of unpolarized TMD PDFs and FFs to data of SIDIS, Drell-Yan and Z boson production at small transverse momentum collected by different experiments. This constitutes the first attempt towards a global fit of f a 1 (x, k 2 ⊥ ) and D a→h 1 (z, P 2 ⊥ ) in the context of TMD factorization and with the implementation of TMD evolution at NLL accuracy.
We extracted unpolarized TMDs using 8059 data points with 11 free parameters using a replica methodology. We selected data with Q 2 > 1.4 GeV 2 and 0.2 < z < 0.7. We restricted our fit to the small transverse momentum region, selecting the maximum value of transverse momentum on the basis of phenomenological considerations (see Sec. III). With these choices, we included regions where TMD factorization could be questioned, but we checked that our results describe very well the regions where TMD factorization is supposed to hold. The average χ 2 /d.o.f. is 1.55 ± 0.05 and can be improved up to 1.02 restricting the kinematic cuts, without changing the parameters (see Sec. IV C). Most of the discrepancies between experimental data and theory comes from the normalization and not from the transverse momentum shape.
Our fit is performed assuming that the intrinsic transverse momentum dependence of TMD PDFs and FFs can be parametrized by a normalized linear combination of a Gaussian and a weighted Gaussian. We considered that the widths of the Gaussians depend on the longitudinal momenta. We neglected a possible flavor dependence. For the nonperturbative component of TMD evolution, we adopted the choice most often used in the literature (see Sec. II C).
We plan to release grids of the parametrizations studied in this work via TMDlib [100] to facilitate phenomenological studies for present and future experiments.
In future studies, different functional forms for all the nonperturbative ingredients should be explored, including also a possible flavor dependence of the intrinsic transverse momenta. A more precise analysis from the perturbative point of view is also needed, which should in principle make it possible to relax the tension in the normalization and to describe data at higher transverse momenta. Moreover, the description at low transverse momentum should be properly matched to the collinear fixed-order calculations at high transverse momentum.
Together with an improved theoretical framework, in order to better understand the formalism more experimental data is needed. It would be particularly useful to extend the coverage in x, z, rapidity, and Q 2 . The 12 GeV physics program at Jefferson Lab [101] will be very important to constrain TMD distributions at large x. Additional data from SIDIS (at Compass, at a future Electron-Ion Collider), Drell-Yan (at Compass, at Fermilab), Z/W production (at LHC, RHIC, and at A Fixed-Target Experiment at the LHC [102]) will be very important. Measurements related to unpolarized TMD FFs at e + e − colliders (at Belle-II, BES-III, at a future International Linear Collider) will be invaluable, since they are presently missing.
Our work focused on quark TMDs. We remark that at present almost nothing is known experimentally about gluon TMDs [11,103], because they typically require higher-energy scattering processes and they are harder to isolate as compared to quark distributions. Several promising measurements have been proposed in order to extract both the unpolarized and linearly polarized gluon TMDs inside an unpolarized proton. The cleanest possibility would be to look at dijet and heavy quark pair production in electron-proton collisions at a future EIC [104,105]. Other proposals include isolated photon-pair production at RHIC [106] and quarkonium production at the LHC [107][108][109][110].
Testing the formalism of TMD factorization and understanding the structure of unpolarized TMDs is only the first crucial step in the exploration of the 3D proton structure in momentum space and this work opens the way to global determinations of TMDs. Building on this, we can proceed to deepen our understanding of hadron structure via polarized structure function and asymmetries (see, e.g., Refs. [111,112] and references therein) and, at the same time, to test the impact of hadron structure in precision measurements at high-energies, such as at the LHC. A detailed mapping of hadron structure is essential to interpret data from hadronic collisions, which are among the most powerful tools to look for footprints of new physics.