Unpolarized Transverse Momentum Distributions from a global fit of Drell-Yan and Semi-Inclusive Deep-Inelastic Scattering data

We present an extraction of unpolarized transverse-momentum-dependent parton distribution and fragmentation functions based on more than two thousands data points from several experiments for two different processes: semi-inclusive deep-inelastic scattering and Drell-Yan production of lepton pairs. The baseline analysis is performed with a Monte Carlo replica method resumming large logarithms at N3LL accuracy. The resulting description of the data is very good ($\chi^2/N_{\rm dat} = 1.06$). For semi-inclusive deep-inelastic scattering, predictions for multiplicities are normalized by factors that cure the discrepancy with data introduced by higher-order perturbative corrections.

In this paper, we present an extraction of Transverse Momentum Distributions (TMDs) [1-4] using more than two thousand data points from several experiments for two different kinds of processes: Semi-Inclusive Deep Inelastic Scattering (SIDIS) and production of Drell-Yan (DY) lepton pairs, significantly improving our previous analysis [5].
Building maps of the internal partonic structure of nucleons is a crucial step towards understanding the interactions between quarks and gluons and the phenomenon of confinement. A steady progress in the last decades has led to more and more refined versions of such maps. The TMDs extracted in this article encode information about the three-dimensional distributions of quarks in momentum space. The level of sophistication of a TMD extraction essentially depends on two ingredients: the amount of analyzed data from different processes, namely how "global" the experimental information from which TMDs are extracted, and the perturbative accuracy reached in the theoretical formalism.
The extraction of TMDs is based on TMD factorization theorems, which provide a precise definition of the objects to be extracted, establishing their universality and evolution equations. In this case, the accuracy of the calculation is defined by the amount of large logarithms being resummed, that in turn defines the powers of the strong coupling α s to be included in the perturbative quantities [6][7][8][9].
Of particular relevance is also the combination of analyzed data from different processes. Similarly to the extraction of collinear Parton Distribution Functions (PDFs), we can talk about a global fit, i.e., a fit that leverages the universality of the parton distributions and uses different processes to constrain them. However, in SIDIS two types of TMDs enter the cross section: the TMD PDFs, describing how partons are arranged in the nucleon, and TMD Fragmentation Functions (FFs), describing how a parton produces a final-state hadron. The knowledge of TMDs, in particular of TMD FFs, would be greatly improved by using data from electron-positron annihilations into two almost back-to-back hadrons [10]. Unfortunately, this data is presently not available. There are measurements for the inclusive production of single hadrons [11] but, in this case, transverse momenta need to be defined with respect to the thrust axis: a careful description of the latter is non trivial, and a rigorous factorization theorem for this process has been discussed only recently (see, e.g., Refs. [12][13][14] and references therein). Therefore, for TMD extractions we currently talk about a global fit when data from SIDIS and DY processes are included.
At present, only two works have reached the stage of combining SIDIS and DY data in a full-fledged global TMD fit: the above mentioned extraction of Ref. [5], henceforth named PV17, and the extraction of Ref. [22], henceforth named SV19.
The PV17 extraction reached the NLL accuracy, was based on the calculation of observables at mean kinematics in each bin, and did not manage to describe the normalization of all datasets. In this work, we push the accuracy of the analysis to what we will refer to as N 3 LL − (only NNLO collinear FFs are currently missing in order to reach full N 3 LL accuracy). 1 Apart from the increase in perturbative accuracy, in this work we also include many measurements published after 2017.
The SV19 extraction reached the same perturbative accuracy and included essentially the same datasets. Our work has crucial differences in the selection of specific data points (we include a much larger number of points), the implementation of TMD evolution, the choice of nonperturbative components, and the handling of normalization for SIDIS data.
As we will describe in detail, we found it particularly difficult to describe in a satisfactory way the normalization of SIDIS data obtained in fixed-target experiments at moderate to low scales. When the analysis is performed at NLL accuracy, we can describe well shape and normalization of SIDIS data, but the description of high-energy DY data is very poor. Going to N 2 LL and N 3 LL − , the description of DY data significantly improves, but we fail to reproduce the normalization of SIDIS data, mainly due to the O(α s ) corrections to the hard factor. As a possible solution, we choose to adjust the normalization of the TMD predictions by comparing their integral upon transverse momentum and the corresponding collinear formula. With this procedure, we fix the normalization a priori, in a way that is independent of the results of the fit.
Our baseline fit is performed at N 3 LL − , using 2031 data points and obtaining χ 2 /N dat = 1.06. We also discuss variations of this baseline fit by changing the theoretical accuracy and the selected data.
The paper is organized as follows. In Sec. II we describe the theoretical framework used in the analysis. In Sec. III we explain how experimental data have been selected. Sec. IV presents the results of our extraction. Finally, in Sec. V we draw our conclusions. Two partons from two hadrons collide. They have transverse momenta k ⊥A and k ⊥B (not measured). They produce a virtual photon with (measured) transverse momentum qT = k ⊥A + k ⊥B with respect to the hadron collision axis.

A. Drell-Yan
In the Drell-Yan (DY) process two hadrons h A and h B with four-momenta P A and P B , respectively, collide with center-of-mass energy squared s = (P A + P B ) 2 , producing a neutral vector boson γ * /Z with four-momentum q and large invariant mass Q = q 2 . The vector boson eventually decays into a lepton and an antilepton with four-momenta constrained by momentum conservation, q = l + l . The involved momenta and the respective transverse components are summarized in Fig. 1. The hadronic four momenta P A and P B can be chosen to identify the longitudinal direction z and define the transverse momentum q T of the γ * /Z. The rapidity of the neutral boson (or, equivalently, of the lepton pair) is defined as For our purposes, we need the cross section for this process initiated by unpolarized hadrons and integrated over the azimuthal angle of the exchanged boson. That cross section can be written in terms of two structure functions, F 1 U U and F 2 U U [30,31]. In the limit M 2 Q 2 (with M the mass of the incoming hadrons) and q 2 T Q 2 , the F 2 U U is suppressed. Accordingly, the cross section reads where α is the electromagnetic coupling, P is the phase space factor to account for potential cuts on the lepton kinematics, which turns out to have a relevant impact when high-precision data are taken into account (see, e.g., a recent analysis in Ref. [32]). 2 At low transverse momentum q 2 T Q 2 , the structure function can be expressed as a convolution over the partonic transverse momenta of two TMD PDFs: In the above equation, H DY is the hard factor, which can be computed order by order in the strong coupling α s and is equal to 1 at leading order. 3 This function encodes the virtual part of the hard scattering and depends on the hard scale Q and on the renormalisation scale µ. The unpolarized TMDs are denoted by f 1 . They depend on the renormalization scale µ and the rapidity scale ζ. The rapidity scales must obey the relation ζ A ζ B = Q 4 . Throughout the paper, we will set The following definition of the Fourier transform of the TMD PDFs has been used: 4 The structure of the TMD PDFs will be addressed in details in Sec. II C. The transverse momentum of the active quark and antiquark are denoted as k ⊥A,B . At low transverse momenta, the two variables x A,B take the values: The summation over a in Eq. (4) runs over the active quarks and antiquarks at the scale Q, and c a (Q 2 ) are the respective electroweak charges given by with where e a , V a , and A a are the electric, vector, and axial charges of the flavor a, respectively; V and A are the vector and axial charges of the lepton ; sin θ W is the weak mixing angle; M Z and Γ Z are mass and width of the Z boson. As discussed in Sec. III and summarized in Tab. II, for DY production the observable provided by the experimental collaborations is the (normalized) cross section differential with respect to |q T |. For each bin delimited by the initial (i) and final (f ) values of kinematical variables, the experimental values are compared with the following theoretical quantity: where the ffl symbol represents the integral divided by the width of the integration range. Hence, Eq. (10) corresponds to the cross section in Eq. (3) averaged over the transverse momentum and integrated over rapidity and invariant mass of the exchanged boson. The normalized cross section is obtained by dividing both sides of Eq. (10) by the appropriate fiducial cross section, which is computed by employing the DYNNLO code [35,36]. 5 The low-energy fixed-target experiments included in this analysis (E288, E605, E772, see Tab. II) measure the following cross section where E and q are the energy and the three-momentum of the photon, respectively. Given Eq. (11), in principle the experimental value in a given bin needs to be compared against the following theoretical quantity: However, since all the considered fixed-target experiments do not provide bins of |q T | but just the average transverse momentum values |q T |, the integration over |q T | is not considered. Moreover, the E288 provides only the average value y for the rapidity. Accordingly, the theoretical quantity considered for that experiment reads The E605 and E772 low-energy fixed-target experiments (see Tab. II) use, in place of the rapidity y, the variable x F , which is connected to the other kinematic variables as follows: Using Eq. (14), one obtains The E772 collaboration provides bins in x F and average transverse momentum values |q T |. Accordingly, in that case the experimental values are compared against the following theoretical quantity: For the sake of simplicity, we replaced y and Q with y and Q in the prefactor in front of the cross section and pull it out of the integral. The E605 experiment, instead, provides average values for both transverse momentum and x F and its data are compared against the following theoretical quantity: where, in this case, In SIDIS, a lepton with momentum l scatters off a hadron target N with mass M and four momentum P . In the final state, the scattered lepton momentum l is measured together with one hadron h with mass M h and four momentum P h . The other products of the scattering are undetected. Thus the reaction reads The (space-like) four-momentum transfer is q = l − l , with Q 2 ≡ −q 2 > 0. We use the standard SIDIS variables [37,38]: For transverse momenta, we will follow the definitions and notations discussed in Ref. [39,40] (see also Fig. 2). In particular, we define P hT as the hadron transverse momentum in the Breit frame, where P and q form a light-cone basis; as a consequence, P h has only transverse components and |P hT | 2 = −P 2 hT , which is frame independent. Similarly, we define q T as the photon transverse momentum in the hadron frame, where P and P h form a light-cone basis; as a consequence, q T has only transverse components and |q T | 2 = −q 2 T , which is also frame independent. The two momenta are related by [41,42] Diagram describing the relevant momenta involved in a SIDIS event in the Breit (nucleon-photon) frame. A virtual photon with momentum q (defining the reference axis) strikes a parton with momentum k inside a nucleon with momentum P . The parton has a transverse momentum k ⊥ (not measured). The struck parton with momentum p = k +q fragments into a hadron with momentum P h , which acquires a further transverse momentum P ⊥ (not measured) with respect to the fragmenting quark axis. 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, In the following, we will always work assuming that the invariant mass of the photon is large compared to the target and hadron masses (M 2 , M 2 h Q 2 ) and to the transverse momenta q T and P hT (q 2 T , P 2 hT Q 2 ). We neglect any power corrections that vanish in this limit, both kinematic and dynamical (higher twist), apart from some modifications to the normalization of the SIDIS observables (that could be seen as the effect of power corrections, see Sec. II D for more details). In this limit, Eq. (21) reduces to Several studies have been made concerning higher-twist corrections of various origin (see, e.g., Refs. [38,43,44] for recent works). A careful study of the impact of power corrections to the case of unpolarized TMDs has been discussed in Ref. [22]. Lately, important advances in the study of higher-twist TMD factorization have been published [42,[45][46][47], but they do not directly affect the observables we consider here.
The differential cross section for SIDIS can be witten in terms of two structure functions, F U U,T and F U U,L [37]. The subscripts refer to the lepton, the target, and the photon polarization, respectively. The second structure function is formally a twist four contribution and is suppressed in the limit considered here, thus we neglect it. The differential cross section at small transverse momentum [5,37] reads where α is the QED coupling constant and is the photon flux factor [37]. Neglecting target mass corrections O(γ) and O(γ 2 ), the prefactors can be approximated as Since we are interested only in the small-transverse-momentum limit, in Eq. (23) we have neglected the contributions from fixed-order calculations at high |q T | [48,49] and the matching of these on TMD factorization [50][51][52].
The unpolarized SIDIS structure function F U U,T is defined as [37] F U U,T x, z, |q T |, Q where the sum runs over quarks and antiquarks a. The hard factor H SIDIS can be computed order by order in the strong coupling α s and is equal to 1 at leading order. 6 The variable k ⊥ is the transverse momentum of the struck quark with respect to the nucleon axis, whereas P ⊥ is the transverse momentum of the produced hadron h with respect to the fragmenting quark axis (see Fig. 2). The variable b T is conjugated via Fourier transform to the transverse momentum q T . f a 1 (x, k 2 ⊥ ; µ, ζ A ) and D a→h 1 (z, P 2 ⊥ ; µ, ζ B ) are the unpolarized TMD PDF for a quark a in a nucleon and the unpolarized TMD FF for a quark with flavor a fragmenting into a hadron with flavor h, respectively;f a 1 (x, b 2 T ; µ, ζ A ) andD a→h 1 (z, b 2 T ; µ, ζ B ) are their Fourier transforms. The former is defined in Eq.(5), the latter is defined aŝ Their structure will be discussed in details in Sec. II C. The observable provided by the HERMES and COMPASS collaborations is the multiplicity, namely the ratio of the one-hadron inclusive cross section as a function of the transverse momentum of the hadron |P hT | over the fully inclusive one: The cross section for unpolarized DIS in the denominator of the multiplicities reads where the approximation is justified by neglecting the target mass corrections. At the perturbative order considered in this analysis, the longitudinal DIS structure function F L cannot be neglected, at variance with, e.g., Refs. [5,15].
The experimental values in each bin are compared against the quantity built by separately averaging the numerator and denominator of the multiplicity in Eq. (27) over the respective kinematics: The HERMES collaboration provides multiplicities in bins of |P hT |, whereas the COMPASS collaboration in bins of P 2 hT (see also Tab. III). In both cases, the observable can be calculated as in Eq. (29), but in the COMPASS case the average is on P 2 hT . Moreover, both collaborations introduce a cut on the invariant mass of the hadronic final states W 2 = (P +q) 2 (see Tab. III), which makes the upper integration limit x f a Q-dependent quantity.

C. Transverse Momentum Distributions (TMDs)
As a consequence of the renormalization of ultraviolet and rapidity divergences [6,53,54], TMD PDFs and FFs acquire a dependence on the renormalization scale µ and on the rapidity scale ζ. The evolution of TMDs from some initial values of the scales µ i , ζ i , to some final values µ f , ζ f , is given bŷ The anomalous dimension γ for the renormalization-group evolution in µ reads: where γ K is the cusp anomalous dimension and γ F α s (µ) = γ µ, µ 2 is the boundary condition [7]. The Collins-Soper kernel K, instead, is the anomalous dimension for the evolution in ζ [6]. The same structure holds for the TMD FF. In order to avoid the insurgence of large logarithms, the scales µ i and ζ i are conveniently fixed Since the coupling α s is computed at this scale (see Eq. (32)) the evolution of the TMD is perturbatively meaningful only at low values of |b T | such that the scale µ b is sufficiently larger than the Landau pole Λ QCD . This condition can be implemented by replacing the scale As suggested by the CSS formalism [6], b * saturates to b max at large |b T | guaranteeing that µ b * never enters the nonperturbative regime. However, this has also the effect of introducing power corrections scaling like (Λ QCD /|q T |) k [55], with k > 0, that in the region |q T | Λ QCD need to be accounted for by a nonperturbative function. At small |b T |, b * saturates to b min . Since µ f is of the order of the boson virtuality Q, this introduces subleading power corrections scaling like (|q T |/Q) n , with n > 0. Such a procedure has the advantage of facilitating a possible matching of the TMD formula, valid for |q T | Q, onto the fixed-order calculation valid for |q T | Q [8,56,57]. Accordingly, in the limit |b T | → 0 the Sudakov exponent vanishes.
Performing at the input scales the Operator Product Expansion (OPE) of the TMD PDFs (TMD FFs) around |b T | = 0 one gets: where the sum runs over quarks, antiquarks, and the gluon. The matching coefficients C are calculated as a perturbative expansion in powers of α s . In view of the power corrections introduced by the b * prescription, both the Collins-Soper kernel K [54] and the OPE in Eq. (35) need to be modified to account for nonperturbative effects. For the Collins-Soper kernel K, this results in a nonperturbative correction term, g K (b 2 T ), for which we choose a specific functional form: This correction gives rise in the evolution to a nonperturbative factor that goes like (ζ f /Q 2 0 ) g K /2 where Q 0 is an arbitrary scale at which this correction is parameterised; we set Q 0 = 1 GeV. In order not to affect the perturbative calculation at small |b T |, the term g K needs to vanish in the limit |b T | → 0. The nonperturbative corrections to the OPE can also be parameterised by a multiplicative function that generally depends on x or z and b T . The net result of the inclusion of the nonperturbative corrections into the evolved TMD PDF reads: and the same holds for the TMD FF where one introduces D 1 N P (z, b 2 T ; ζ, Q 0 ). Note that the number of active flavors n f in the perturbative quantities γ, C, and the hard function H of Eqs. (4) and (25), is separately determined by the scales µ, µ b * and Q, respectively. To be more precise, given a set of quark thresholds {m 1 , m 2 , m 3 , . . . }, the n f associated to each of the three scales above is computed by requiring that the scale lies between m n f and m n f +1 . Analogously, since the collinear distributions involved in the matching formula in Eq. (35) are computed at the scale µ b * , the value of n f for PDFs and FFs is chosen accordingly, i.e. it is the same used for the matching functions C.
The f 1 N P and D 1 N P factors (which we assume to be flavor-independent) incorporate both the correction to the evolution associated to the g K function and the correction to the respective OPE. For the TMD PDF we define and for the TMD FF the form is The nonperturbative factors f 1 N P , D 1 N P → 1 for b T → 0. The g i functions describe the dependence of the widths of the distributions on x and z: In total, the default configuration for the fit involves 21 free parameters: one associated to the nonperturbative part of the Collins-Soper kernel (Eq. (36)), 11 related to the nonperturbative part of the TMD PDF (Eqs. (38), (40)), and 9 for the nonperturbative part of the TMD FF (Eqs. (39), (41)).
The functional forms in Eqs. (38)-(41) are largely arbitrary. However, an important feature is that they are the Fourier transforms of the sum of a Gaussian, a weighted Gaussian (multiplied by k 2 ⊥ ) and, in the case of the TMD PDFs, a third Gaussian. They are therefore positive definite for all values of k 2 ⊥ . 7 The parameters λ and λ 2 in Eq. (38) are squared in order to avoid negative contributions (there is no need to square the parameter λ F in Eq. (39) because the fit always selects positive values for this parameter). The widths of the Gaussians, expressed by Eqs. (40), (41), are x (or z) dependent and vanish as x (or z) approaches 1. Our choice of the functional form is also inspired by model calculations of TMD PDFs (see, e.g., [58][59][60][61][62][63][64][65]) and TMD FFs (see, e.g, [66,67]). Many of these models predict the existence of terms that behave similarly to Gaussians and weighted Gaussians. The details of the functional dependence predicted by the models are related to the correlation between the spin of the quarks and their transverse momentum. In the case of fragmentaton functions, a different role can be played by different producton channels (e.g., direct production vs. production through the decay of hadronic resonances).
Finally, for the logarithmic ordering we use the same convention adopted in Ref. [7]. In particular, the orders of truncation of the perturbative ingredients relevant to the present analysis are summarized in Tab. I. At the time of this analysis, the full N 3 LL accuracy could not be achieved because NNLO collinear FFs were not available. Very recently, two analyses of collinear FFs were presented in Refs. [28,29] making an extraction of TMDs at full N 3 LL possible. We leave this study to a future pubblication.
Accuracy H and C K and γF γK PDF and αs evolution FF evolution NLL 0

D. Normalization factors for SIDIS
In Ref. [5] it was demonstrated that TMD factorization at NLL accuracy is able to successfully reproduce the normalization and shape of HERMES SIDIS multiplicities and the shape of the available COMPASS multiplicities. More recently, the COMPASS collaboration published a reanalysis of their data [68]. The NLL TMD predictions correctly reproduce normalization and shape of the new data [69]. However, when increasing the accuracy to N 2 LL or higher, the TMD formula severely underestimates the measurements [69,70] by nearly constant factors in each bin. Note that tensions between the TMD cross sections and the associated measurements exist also at large transverse momentum in SIDIS [48], DY [71], and electron-positron annihilation into two hadrons [72].
In the present study, as will be shown in Sec. IV, we confirm that we obtain an excellent description of both normalization and shape of the SIDIS multiplicities at NLL and that the N 2 LL results are much smaller. At N 3 LL the results increase slightly, but they are still far from the NLL ones and, therefore, from data. At average kinematics of the COMPASS measurements, we obtain the following ratios of multiplicities: The reason for the difference between the logarithmic orders is almost entirely due to the hard factor in Eq. (25) [34]. If we look at the explicit expression for the hard factor 8 with the standard choice µ = Q we can immediately see that just by introducing O(α s ) corrections, at Q = 2 GeV with α s ≈ 0.3, we reduce the structure function to about 60% of its original value. This change is not compensated by a similar reduction in the denominator of Eq. (27): the differences between the LO and NLO expressions of the inclusive DIS cross section are typically below 5% and the NLO results are actually larger than the LO ones.
If the NLL expression is much larger than the N 2 LL and N 3 LL ones, we may suspect that it should overshoot the data by a factor 1.5 at least. However, several works before the present one have shown a good agreement with data using a parton-model approach [15,16] or at NLL [17,18]. Moreover, the integral of the structure function over q T is equal to the value of its Fourier transform at b T = 0. Using any b * prescription with b min = 2e −γ E /Q, the integral of the NLL expression by construction corresponds to the LO expression of the collinear SIDIS structure function, independent of the TMD nonperturbative parameters: 9 The LO collinear SIDIS predictions are known to describe the data reasonably well and, if anything, they seem to be lower than the data [73,74]. Therefore, the integral of our NLL expression is in good agreement with data, which also indicates the absence of a large normalization error.
If the NLL predictions describe the data well and are a factor 2 or 1.5 above the N 2 LL and N 3 LL predictions, we propose to modify the normalization of the latter to recover a good agreement with data. An extended discussion of this issue can be found in Ref. [69]. We observe that the integral of the TMD formula, valid at low q T , should reproduce only part of the full collinear cross section. The only exception is the order O(α 0 s ) case, as we have seen above, since at that order there is no contribution from gluon radiation at high transverse momentum, beyond the TMD region.
However, at N 2 LL or higher, in the kinematics of fixed-target SIDIS experiments, the integral of the TMD region (i.e., the integral of the so-called W term in the language of Ref. [6]) is much smaller than the corresponding collinear cross section. The missing contribution to the integral should be recovered by the terms in the fixed-order calculation that are not included in the TMD resummed expression (the so-called Y term). Ideally, the Y term should be negligible in the low-q T region. This is not the case in the experimental regions under consideration: the Y term is finite but relatively large, even at q T = 0.
If we consider the contribution to the integral of the W term (i.e., the integral of Eq. (23)), with our b min prescription, at order α s we obtain, schematicallŷ where The double convolution over both x and z is defined as and the C TMD coefficients can be found in App. A. 8 There are different definitions for the hard factor in the literature, which are compensated by different definitions of the matching coefficients C in Eq. (35). Here we follow the definition of Ref. [34]. 9 Note that in the absence of a b min prescription, the integral would vanish.
The integral in Eq. (45) should be compared to the collinear expression at the same order (see, e.g., Ref. [75]) The QCD coefficients C ab 1 and C ab L are calculated in perturbation theory. The former can be written as The coefficient C nomix is the sum of all those terms that contain either a δ(1 − x) or a δ(1 − z) (or both). Some of these terms are present in C TMD , but not all. This definition holds at all orders in α s . The C L matching coefficients, instead, only contain "mixed" contributions. For convenience, we reproduce all coefficients at order α s in App. A. In order to increase the size of the TMD component, we consider the contribution of all the "nonmixed" terms C nomix . The reason behind this choice is that it might be possible to include such terms into a redefinition of the individual TMDs. Hence, we define and similarly for higher orders, and we introduce the following normalization factor: We stress that these normalization factors depend only on the collinear PDFs and FFs, are independent of the parametrization of the nonperturbative part of the TMDs, and can be computed before performing a fit of the latter. At NLL, ω(x, z, Q) = 1. Beyond NLL, the prefactor becomes larger than one and guarantees that the integral of the TMD part of the cross section reproduces most of the collinear cross section, as suggested by the data. On the contrary, without the enhancement due to the normalization factor, the integral of the TMD part of the cross section would be too small, requiring a large role of the high-transverse-momentum tail, which is not observed in the data. The impact of the normalization factor defined in Eq. (51) will be addressed in detail in Sec. IV.
As a consequence of our procedure, the theoretical expression for the SIDIS cross section in Eq. (23) becomes

III. DATA SELECTION
In this Section we describe the experimental data included in our global analysis. We consider a large number of datasets related to DY lepton pair production and SIDIS, for the observables discussed in Sec. II A and Sec. II B. The coverage in the x-Q 2 plane spanned by these datasets is illustrated in Fig. 3.
The majority of datasets analyzed in the present work was already included in the global analysis of SIDIS and DY data in Ref. [5] and in the fit of DY data discussed in Ref. [7]. For more details, we refer the reader to those references. The new datasets included in the present analysis are: • DY di-muon production from the collision of a proton beam with an energy of 800 GeV on a 2 H fixed target from E772 ( √ s = 38.8 GeV) [76]; • DY di-muon production from the PHENIX Collaboration [77]; • DY data at 13 TeV from the CMS Collaboration [78] and the ATLAS Collaboration [79].

A. Drell-Yan
Our analysis is based on TMD factorization, which is applicable only in the region |q T | Q. Therefore, in agreement with the choices of Refs. [7,22] we impose the following cut Table II summarizes all the DY datasets included in our analysis. For some DY datasets the experimental observable is given within a fiducial region. This means that kinematic cuts on transverse momentum p T and pseudo-rapidity η of the single final-state leptons are enforced (values reported in the next-to-last column of Tab. II). For more details we refer the reader to Ref. [7]. The second column of Tab. II reports, for each experiment, the number of data points (N dat ) that survive the kinematic cuts. The total number of DY data points considered in this work is 484. Note that for E605 and E288 at 400 GeV we have excluded the bin in Q containing the Υ resonance (Q 9.5 GeV).
As can be seen in Tab. II, the cross sections are released in different forms: some of them are normalized to the total (fiducial) cross section while others are not. When necessary, the required total cross section σ is computed using the code DYNNLO [35,36] with the MMHT14 collinear PDF set, consistently with the perturbative order of the differential cross section (see also Tab. I). More precisely, the total cross section is computed at NLO for NNLL accuracy, and NNLO for N 3 LL − accuracy. The values of the total cross sections at different orders can be found in Table 3 of Ref. [7]. For the ATLAS dataset at 13 TeV, the value of the fiducial cross section is 694.3 pb at NLO and 707.3 pb at NNLO.

B. SIDIS
The identification of the TMD region in SIDIS is not a trivial task and may be subject to revision as new data appears and the theoretical description is improved, as discussed in dedicated studies [38,94,95].
First of all, a cut in the virtuality Q of the exchanged photon is necessary to respect the condition Q Λ QCD needed for perturbation theory to be applicable. In this way also mass corrections and higher twist corrections can be neglected. In this work, we require that Q > 1.4 GeV. Studies of SIDIS in collinear kinematics employ similar cuts [29,96].
In order to restrict ourselves to the SIDIS current fragmentation region and interpret the observables in terms of parton distribution and fragmentation functions, we apply a cut in the kinematic variable z by requiring 0.2 < z < 0.7. The lower limit is the same used in the study of collinear fragmentation functions [29,96]. We used a slightly more restrictive upper limit, to avoid contributions from exclusive channels and to focus on a region where the collinear fragmentation functions have small relative uncertainties.
For what concerns the cut on transverse momentum, our baseline choice is with fixed parameters c 1 = 0.2, c 2 = 0.5 and c 3 = 0.3. This choice is more restrictive than a similar one made in Ref. [5], but less restrictive than the one made in Ref. [22]. It allows for many data points with |P hT | Q but also with 0.2 Q < |q T | < Q. In Sec. IV, we will discuss variations of the baseline SIDIS cut in Eq. (54) that give phenomenological support to our choice.
As for the datasets included in the present analysis, the main difference with Ref. [5] is that we include the new release of COMPASS data [68]. In this dataset, the vector-boson contributions have been subtracted. For   (53)), the observable delivered, the center-of-mass energy √ s, the range(s) in invariant mass Q, the angular variable (either y or xF ), possible cuts on the single final-state leptons, and the published reference (when available). The total number of DY data points amounts to 484.
the HERMES dataset we consistently select the vector-meson-subtracted dataset (.vmsub set). Moreover, we select the zxpt-3D-binning for HERMES multiplicities, since it presents a finer binning in |P hT |. The breakdown of the entire SIDIS dataset included in the present analysis is reported in Tab. III.
The second column of Tab. III shows the number of data points (N dat ) that respect the kinematic cuts for each dataset, with a total number of 1547 data points.
In conclusion, the total number of DY and SIDIS data points surviving our kinematic cuts is 2031.

C. Error treatment
The considered experimental datasets are released with a set of systematic and statistical uncertainties. As already pointed out in Ref. [7], a proper treatment of the experimental uncertainties is extremely important in order to obtain a reliable extraction of TMDs. Thus, we choose to treat systematic uncertainties as fully correlated only if, in the corresponding publication, it is explicitly specified that they are correlated. The statistical uncertainties, instead, are always considered as uncorrelated.
At variance with Ref. [7], in this analysis we do not make use of the iterative t 0 -prescription [97] for the  treatment of correlated normalization uncertainties. This prescription is usually introduced to avoid the underestimation of the predictions caused by the so-called D'Agostini bias [98]. After performing the fit with and without the t 0 -prescription, we found that our analysis is not affected by the D'Agostini bias and therefore we saw no reason to introduce the t 0 -prescription in the computation of the χ 2 .
On top of experimental systematic uncertainties, there can be several sources of systematic theoretical errors. The first one comes from the choice of the underlying collinear parton distribution and fragmentation functions. In our case we choose the MMHT2014 [99] collinear PDFs and the DSS collinear FFs. Since the HERMES collaboration provides multiplicities for pions and kaons separately (see Tab. III), we use DSS14 [100] for π ± and DSS17 [96] for K ± . Given the nature of the PDF and FF set used in this analysis, their uncertainties are computed using the Hessian method [100][101][102]. We observed that PDF and FF uncertainties are significantly correlated across bins. In order to account for this correlation, we decomposed the corresponding uncertainties into a fully correlated part that amounts to 80% of the total, while we treated the remaining 60% as uncorrelated. 10 Moreover, for the observables measured by COMPASS one needs a collinear set of FFs for unidentified charged hadrons. Since the dedicated DSS07 FF set for unidentified charged hadrons [103] does not provide an estimate of the uncertainties, we computed the COMPASS multiplicities by using the sum of the DSS14 and DSS17 sets for pions and kaons. 11 The associated uncertainty is calculated by propagating to the multiplicity the Hessian errors associated to each of the two hadronic component. As pointed out in Ref. [22], the choice of specific sets for the collinear distributions may have a sizeable impact on the final result. In our analysis, we did not consider alternative collinear sets and postpone this study to a future publication. Likewise, we leave for future work the study of other sources of theoretical uncertainties such as higher-twist corrections, TMD flavor dependence and the choice of the perturbative scales.

IV. RESULTS
In this section, we present the results obtained for the extraction of unpolarized quark TMDs from a global analysis including both DY and SIDIS data (see Sec. III). This work represents an important upgrade with respect to Ref. [7]. In fact, in the SIDIS process a single hadron is measured in the final state, allowing us to extract information about fragmentation functions. Therefore, the final result of this work is the simultaneous extraction of both TMD PDFs and FFs at N 3 LL − . At the moment, this is the most precise extraction of TMDs that has been achieved on more than two thousand data points. In Sec. IV A we present the quality of the fit, in Sec. IV B we discuss the extracted TMD distributions, and in Sec. IV C we investigate the effect of variations on the baseline fit configuration.

A. Fit quality
In this section, we discuss the quality of the baseline fit performed at N 3 LL − imposing the cuts on |q T |/Q as discussed in Sec. III. The error analysis is performed with the so-called bootstrap method, which consists in fitting an ensemble of 250 Monte Carlo replicas of the experimental data (see Ref. [5] for more details).
The most complete statistical information about the TMDs is given by the full ensemble of 250 replicas. For some purposes, however, it is useful to define a single, representative result instead of the full replica ensemble. In order to estimate the quality of our fit, the most appropriate indicator is the χ 2 value of the best fit to the central (not fluctuated) experimental data (χ 2 0 ). We refer to this fit as the "central replica." It is possible to analyze also the average of the χ 2 over all replicas ( χ 2 ) as well as the χ 2 of the mean replica (χ 2 m ), constructed as the average of all replicas [7]. These values should be very close to each other. Tab. IV reports the breakdown of the χ 2 0 values normalized to the number of data points (N dat ) for DY and SIDIS datasets. As already discussed in Ref. [7], in the presence of bin-by-bin correlated uncertainties the total χ 2 can be expressed as the sum of two contributions where χ 2 D is given by the standard formula for N experimental data points exp i and statistical and uncorrelated uncertainties σ 2 i = σ 2 i,stat + σ 2 i,uncor , but involving theoretical predictions shifted by the correlated uncertainties, where σ

(α)
i,corr is the α-th (100%) correlated uncertainty associated with the i-th experimental data point and λ α is the so-called nuisance parameter. In Eq. (55), the χ 2 λ is a penalty term due to the presence of correlated uncertainties and is completely determined by the nuisance parameters, The optimal value of the nuisance parameters is obtained by minimizing the total χ 2 in Eq. (55) with respect to them. Because the shifted predictions in Eq. (56) provide a better visual assessment of the fit quality, we will consistently display them for all observables used in the fit. From Tab. IV, the global χ 2 0 value is 1.06, indicating that the description of the whole dataset is very good. 12 This means that the fit is able to simultaneously describe experimental data coming from two different processes over a wide kinematic range. As can be seen in Tabs. II and III, the low-energy dataset comes from fixed-target DY experiments and SIDIS observables, while the high-energy dataset comes from collider experiments at the LHC and Tevatron at energies higher by more than two orders of magnitude.
It is important to notice that the correlated penalty term χ 2 λ gives a significant contribution to the total χ 2 0 value. This means that the shifts induced by correlated uncertainties are often large. The χ 2 λ /N dat = 0.29 obtained in this analysis is larger than the one in Ref. [7], mainly because of the different treatment of the theoretical uncertainties related to collinear PDFs and FFs, which are considered here as 80% correlated.

SIDIS
For HERMES and COMPASS multiplicities, which represent about 75% of the total number of data points considered in this work, the description obtained by our global analysis is very good. From Tab. IV, the values of χ 2 0 /N dat are almost always smaller than 1. In the case of HERMES multiplicities, we note that the largest contribution to the χ 2 0 comes from the π + channel, for both proton and deuteron targets. This result is consistent with the findings in both Refs. [5,15] and [22]. Since kaon multiplicities are affected by larger statistical errors, and the collinear FFs for kaons display large uncertainties, the corresponding χ 2 0 value is lower. The comparison between theoretical results for the SIDIS multiplicities of Eq. (30) and HERMES data for the production of charged pions and kaons off a deuteron target is shown in Fig. 4. Each column corresponds to a specific x bin. Each row corresponds to a specific final-state channel. The results are displayed as functions   of the transverse momentum |P hT | of the measured final-state hadron. Points with different markers and colors correspond to different representative z bins, and are offset for a better visualization as indicated in the plot legend. The light blue rectangles are the theoretical results and correspond to the 68% Confidence-Level (CL) band (namely excluding the largest and smallest 16% of the replicas).
We note that for K − production the central column, corresponding to the 0.2 < x < 0.35 bin, does not include the magenta points for the highest 0.6 < z < 0.8 bin because of the kinematic cut in z. Similarly, in all panels there are only three green points (for the lowest z bin) because of the z-dependent cut in Eq. (54), which leads to the exclusion of a larger number of |q T | bins at lower values of z. The theoretical results display larger uncertainties for K − production (second row) because of the combined effect of larger experimental errors and larger uncertainties in the kaon collinear FFs.     5 refers to the same HERMES multiplicities with same conventions and notation as in Fig. 4 but off a proton target. We remark that for K − production the kinematic cuts have a more drastic effect because the magenta points for the |P hT | distributions at the largest z bin are excluded for the two largest x bins considered (central and rightmost panels of second row from top).

Multiplicity [GeV
In Fig. 6, we show the result of our fit for the COMPASS SIDIS multiplicities for the production of unidentified negatively charged hadrons off a deuteron target. For each Q and x bin, each panel displays the multiplicity on a logarithmic scale as a function of P 2 hT /Q 2 . Again, points with different markers and colors correspond to different representative z bins, as indicated in the plot legend. As before, the light-blue rectangles correspond to the 68% CL theoretical results. The results for unidentified negatively charged hadrons h − are obtained by simply adding the results for negatively charged pions and kaons, h − ∼ π − + K − .
We note that the agreement is good for almost all bins, which is reflected in small χ 2 0 values in Tab. IV. The    situation worsens for the lowest Q bin (1.3 < Q < 1.73 GeV), particularly for x 0.02. We also remark that for some Q and x bins the theoretical uncertainties for the largest z bin compatible with our kinematic cuts (0.6 < z < 0.8) are significantly larger than for other z bins because of much larger uncertainties in the collinear FFs. Finally, looking at the table of panels from top to bottom, we realize that the lowest z-bin distributions (red diamonds) are present only for the largest Q bins, and vice-versa, because of the kinematic cut in Eq. (54) Fig. 7 refers to the same COMPASS multiplicities with same conventions and notation as in Fig. 6 but for unidentified positively charged hadrons h + . Again, the light-blue rectangles correspond to the 68% CL theoretical results, and are obtained by adding the results for positively charged pions and kaons, h + ∼ π + +K + . Comments similar to Fig. 6 can be made about the agreement between data and theory.

Drell-Yan
DY data represents approximately 25% of the full set of analyzed data. From Tab. IV it is evident that most of low-energy DY data from fixed-target experiments (E605, E288, E772), but also from PHENIX and STAR, can be fitted with low χ 2 values, much lower than high-energy DY data from collider experiments like, e.g., those at the LHC. As already pointed out in Ref. [7], this most likely originates from the fact that low-energy DY data are affected by larger errors and collinear PDFs at these kinematics have larger uncertainties.
From Tab. IV, we also note that the quality of our fit for the ATLAS datasets is poor. In particular, the description worsens for the first two low-rapidity bins of both ATLAS 7 TeV and ATLAS 8 TeV datasets, the worst case being at |y| < 1 for ATLAS 7 TeV. Several effects might be responsible for this result. Since the experimental observable is a normalized cross section, systematic errors cancel in the ratio producing measurements with very small error bars. Fitting these data is very difficult, also because small theoretical effects can give significant contributions to the χ 2 . Moreover, different implementations of phase-space cuts on the final-state leptons could lead to modifications in both the shape and the normalization of the theoretical observable (see, e.g., Ref. [32,105,106]). We leave this issue for future studies. At variance with Ref. [22], we obtain our results without excluding any extra data points on top of the ones exceeding the maximum value of |q T |/Q in Eq. (53). It is interesting to comment the results of the fit for those datasets that were not included in the previous analysis of Ref. [7] (see Sec. III). For E772, we are able to obtain a good description only for data points above the peak of the Υ resonance. For Q < 9 GeV, the quality of the fit worsens. At variance with Ref. [22], we keep the Q < 9 GeV bins because there is no evident motivation to exclude them.
The new CMS dataset at √ s = 13 TeV is important because it extends the kinematic coverage considered in Ref. [7]. This dataset is very nicely described, even better than the ones at the smaller center-of-mass energies √ s = 7, 8 TeV. This is probably due to the fact that the CMS 13 TeV dataset is densely binned in rapidity. In order to visualize the quality of our fit, in Figs. 8-11 we present the comparison between experimental data and theoretical results for a representative selection of the DY dataset. In the upper panels of each plot, we display the cross section differential in |q T |, while in the lower panels we show the ratio of data to theory. As for the SIDIS case, the light-blue bands are the 68% CL theoretical results. Fig. 8 displays the DY cross section for the E288 dataset at beam energy E beam = 200 GeV for different bins in Q. For DY fixed-target observables, we calculate the cross section at mean values of |q T | (not integrating upon the bin limits). Hence, we display the 68% CL uncertainty as a band rather than a series of rectangles. We remark that the uncertainty band is larger for lower Q bins; this trend is induced by larger correlated uncertainties for smaller invariant masses of the lepton pair. In Fig. 9, we compare the theoretical results for the cross section for DY in pp collisions at the Tevatron. Black data points in the left panel correspond to the results of Run I of the CDF experiment, while in the right panel the results for Run II are reported. The lower panels show the corresponding ratio of experimental data to theoretical results. The latter are displayed as light-blue rectangles, each one corresponding to the integral of the cross section within the corresponding bin limits. The size of the rectangle is given by the 68% CL. The quality of the fit for CDF data is comparable to the one in Ref. [7]. In Fig. 10, we compare the theoretical results for the DY cross section in pp collisions at the LHC. From left to right, the black data points refer to the measurements by the CMS Collaboration at increasing √ s = 7, 8, 13 TeV. For √ s = 7, 8 TeV, the results are normalized to the fiducial cross section. As in previous figures, the lower panels display the ratio of experimental data to the 68% CL theoretical results. For √ s = 7, 8 TeV, the quality of the fit is comparable to that in Ref. [7]. For the new dataset at √ s = 13 TeV, the agreement in the displayed rapidity bin is excellent, but remains very good also for higher rapidities. In Fig. 11, we compare the theoretical results for the DY cross section in pp collisions normalized to the fiducial cross section for the ATLAS data at √ s = 7 TeV. From left to right, we consider three representative bins at increasing rapidity |y|. The leftmost one corresponds to the worst described bin in our global fit, with χ 2 /N data = 13.5. As in previous figures, the lower part of each panel displays the ratio of experimental data to the 68% CL theoretical results. The quality of the fit increases at more forward rapidities (from left to right). The same trend is observed at √ s = 8 TeV, but not for CMS at 13 TeV.

B. TMD distributions
We now discuss the TMD distributions extracted from our baseline fit with N 3 LL − accuracy. Tab. V displays the list of our 21 fitting parameters with their mean value and standard deviation. The majority of the parameters is well constrained. The only parameter that is compatible with zero is γ 2 .
Parameter Average over replicas 0.316 ± 0.025 α 1 1.29 ± 0.19 0.0055 ± 0.0006 β 1 10.23 ± 0.29 0.455 ± 0.050 σ 3 12.71 ± 0.21 β 2 4.17 ± 0.13 δ 2 0.167 ± 0.006 γ 2 0.0007 ± 0.0110 The λ parameter measures the relative weight of the first Gaussian and the weighted-Gaussian in the nonperturbative part of the TMD PDF in Eq. (38). The value of this parameter is close to 2, indicating that the contribution of the weighted-Gaussian component is important. In Eq. (38), the parameter λ 2 measures the relative weight of the first Gaussian and the third Gaussian; this parameter is small but not compatible with zero, which means that also this component of the TMD PDF is important to reach a good description of experimental data.
Our parametrization of the nonperturbative part of TMD FFs in Eq. (39) contains just the combination of a Gaussian and a weighted Gaussian: this is sufficient to describe the data in an accurate way.
The λ F parameter measures the relative weight of the two components; its value is close to 0.1, indicating that the contribution of the weighted Gaussian is small. Nevertheless, it has non-trivial consequences on the tail of the TMD FF, as we will show below.
The g 2 parameter is a key ingredient to the extraction of the Collins-Soper kernel, discussed in Sec. IV B 1. The same parameter was used in the analysis of Ref. [5]. It is interesting to observe that the value obtained in the present global fit is smaller by almost a factor of 4 with respect to Ref. [5]. This may be due to the higher theoretical accuracy of the present analysis and to the role of the very precise high-energy Drell-Yan measurements, which also determine the very small standard deviation of g 2 .
In Fig. 12, we show a graphical representation of the correlations among the 21 fitting parameters. Using the color code indicated in the legend, it is easy to realize that the nondiagonal elements are very small except for some (anti-)correlation among the β 1 , δ 1 and γ 1 parameters that control the z-dependent width of the Gaussians in the TMD FF (see Eqs. (39), (41)). The overall absence of large correlations suggests that the model parametrization of the non perturbative parts of TMDs is appropriate.
In Fig. 13, we show the unpolarized TMD PDF for the up quark in the proton at µ = √ ζ = Q = 2 GeV (left panel) and 10 GeV (right panel) as a function of the quark transverse momentum |k ⊥ | for three different values of x, namely x = 0.001, 0.01, and 0.1. The bands correspond to the 68% CL.
The TMD seems to be wider at intermediate x = 0.01, but has also a high tail at x = 0.001. As already mentioned, a significant role is played by the weighted Gaussian and by the second Gaussian in Eq. (38). This may be a sign of the presence of contributions from different quark flavors and/or from different spin configurations (see Sec. II C).
It is worth noticing that in both left and right panels the TMD PDF at x = 0.001 shows the largest error band, particularly at low |k ⊥ |. This is due to the lack of experimental points in that kinematic region (see

Collins-Soper kernel
It is interesting to study the Collins-Soper kernel [6,109] that drives the evolution of TMDs in terms of the rapidity scale ζ. Recent discussions of this crucial component of the TMD formalism have been presented in Refs. [110,111] and estimates based on lattice QCD have been proposed in Refs. [112][113][114].
The Collins-Soper kernel, as written in Eq. (36), is composed of two parts. The first part can be calculated perturbatively at N k LL accuracy, and is computed at b * : where K (n,0) and γ (n) K are coefficients of the perturbative expansion (see, e.g., Ref. [34]). Note that the integral on the r.h.s. is directly computed by means of numerical integration, thus providing a fully resummed result. The second part, denoted as g K , cannot be computed in perturbation theory and is one of the results of our fit. Only the full Collins-Soper kernel can be compared with other works.
In Fig. 15, we show the Collins-Soper kernel as a function of |b T | by conventionally keeping the scale µ fixed at 2 GeV, for our present analysis (MAPTMD22, green band) and for four other analyses in the literature [5,7,20,22]. The solid lines at low |b T | follow the perturbative result. For MAPTMD22, PV19 [7] and PV17 [5], they correspond to setting b min = 0 for sake of comparison with the other SV19 [22], SV17 [20] results. The slight differences between the curves are due to the different logarithmic accuracies of the perturbative calculations: the PV17 analysis was performed at NLL, the SV17 analysis at N 2 LL, the PV19, SV19 and MAPTMD22 at N 3 LL. The size of the bands around the solid lines corresponds to one standard deviation of the parameter g 2 around its best-fit value. The b * prescription modifies the curves starting from |b T | ≈ 1 GeV −1 . The behavior at high |b T | is driven by g K and is different for the various analyses.
The dashed curves show the effect of using our prescription b min = 2e −γ E /µ ≈ 1.123/µ in MAPTMD22, PV19 and PV17. This implies that at low |b T | the Collins-Soper kernel saturates to a finite value, as indicated by the dashed lines. As the scale increases, this modification occurs at lower and lower values of |b T | and becomes less relevant.
In the case of the TMD PDF for a quark q in the proton at µ = √ ζ = Q, one has [115,116]: where the Fourier transformf q 1 of the TMD PDF has been defined in Eq. , compared with the PV17 [5], SV17 [20], PV19 [7], and SV19 [22] analyses. For the MAPTMD22, PV17, and PV19 curves, the uncertainty bands represent the 68% CL. The corresponding dashed lines show the effect of including the bmin-prescription (see text).

the TMD PDFf
is defined as [115]: (60) In order to obtain meaningful values for the average squared transverse momenta, i.e., finite, positive, and not dominated by the perturbative tails of the TMDs, we shift the value of |b T | from 0 to a value well inside the nonperturbative region [116]. In this way the Bessel functions J 0,1 tame the power-law behavior of the TMD at large transverse momentum. The choice of the specific value for |b T | is of course arbitrary and has a significant effect on the associated numerical values. We choose |b T | = 1.5 b max that guarantees that the average squared transverse momenta are positive across the x, Q values considered in the fit. Accordingly, Eq. (60) becomes: where the subscript r stands for regularized. We have checked that the results are consistent when choosing either the integral or differential expressions in Eq. (60). The same arguments apply to the regularized average squared transverse momentum produced during the hadronization of the quark q into the final state hadron h [42,115,116]: where the Fourier transformD q→h 1 of the TMD FF is defined in Eq. (26) and the first Bessel moment of the TMD FFD q→h (1) 1 is defined as [42]: In Fig. 16, we show the scatter plot of k 2 ⊥ u r (up quark contribution) at x = 0.1 vs. P 2 ⊥ u→π + r ("favored" fragmentation) at z = 0.5. The blue circles (denoted by MAPTMD22) correspond to the single replicas while the red square is the average over all replicas for the N 3 LL − analysis. Three different values of Q = 1, 2, 4.75 GeV are included to show the evolution of the average transverse momenta with the scale. The orange circles indicate the PV17 replicas [5] at NLL and Q = 1 GeV with the black cross being the average. No regularization is needed for the values extracted in the PV17 analysis, since the involved TMDs at Q = 1 reduce entirely to their nonperturbative components. By comparing MAPTMD22 at Q = 1 GeV to PV17, we observe that the former produces much less anti-correlation between k 2 ⊥ and P 2 ⊥ than the latter, probably because of the inclusion of very precise DY data.

C. Variations on the fit configurations
In this subsection we discuss the results obtained by modifying some of the baseline settings. In Sects. IV C 1 and IV C 2 we present fits at NNLL and NLL accuracy, respectively, comparing them to the baseline N 3 LL − . Finally, in Sec. IV C 3 we study the impact of adopting different cuts in |q T | on the SIDIS dataset.
This is the case of the Sivers TMDs where the computation of the polarized cross section for the Sivers effect presently cannot go beyond the NNLL level [117][118][119], hence demanding unpolarised TMDs at the same level of accuracy.
To this aim, we perform a new global fit at NNLL. However, when lowering the perturbative accuracy, it is possible to obtain acceptably good fits only by excluding those datasets whose precision requires the highest theoretical accuracy. Specifically, we found that only by removing the ATLAS dataset we were able to achieve an acceptable global description at NNLL accuracy. As a matter of fact, in Tab. VI the value of χ 2 in this configuration, namely for fixed-target DY and SIDIS, is lower than at N 3 LL − where ATLAS data is included.
Because of the difference in the perturbative accuracy as well as in the dataset, we do not expect to get compatible values for the best fit parameters between the NNLL and N 3 LL − fits. For instance, we obtain λ = 12 ± 10 GeV −1 and λ F = 340 ± 280 GeV −2 at NNLL, to be compared to λ = 1.8 ± 0.3 GeV −1 and λ F = 0.08 ± 0.01 GeV −2 at N 3 LL − .
The λ and λ F parameters control the relative weight of the weighted Gaussian in the non perturbative part of the TMD PDF and FF, respectively, and control the size of the DY and SIDIS spectrum at middle to large values of |q T |. The large values obtained in the NNLL imply that the weighted Gaussian dominates for both TMD PDF and FF parametrizations. This behavior may be partially induced by the lack of perturbative corrections of the NNLL fit with respect to the N 3 LL − one, which are compensated by nonperturbative effects.

Global fit at NLL
We performed also a global analysis at NLL accuracy. Similarly to the NNLL fit discussed above, it might be useful to have also a global fit at NLL accuracy for contexts where it is not possible to reach the highest accuracy. However, in order to obtain acceptable χ 2 values at this order, we had to exclude all collider DY data and the E772 fixed-target DY dataset. Thus, we reduced the dataset to the SIDIS data and the remaining fixed-target DY data, i.e., E605 and E288.
We point out that in our approach the integral of the W -term at NLL is equal by construction to the SIDIS q T -integrated collinear cross section. As a consequence, the value of the prefactor ω in Eq. (52) is automatically equal to 1. Moreover, for this fit we consistently used MMHT2014 at LO for the collinear PDFs and we used DSS at NLO for collinear FFs. As shown in Tab. I, at NLL both collinear PDFs and FFs should be evaluated at LO, but we choose FFs at NLO because no recent extractions at LO are currently available.
As can be seen in Tab. VI, we obtain low χ 2 values for all included datasets. It is interesting to compare this result with that of the PV17 analysis [5]. In that work, the normalization of COMPASS data was fixed by the first bin in P 2 hT . Here we demonstrate that we can obtain an excellent description of the most recent COMPASS data without any adjustment of the normalization.

Cut in |qT |/Q
A crucial ingredient of any phenomenological analysis of TMDs is the introduction of the kinematic cut |q T |/Q 1 on the dataset to ensure that TMD factorization is valid. Our default choices are discussed in Sec. III. It is interesting to study how the global quality of our fit changes upon variations of this cut in order to get quantitative information on the range of validity of TMD factorization.
To this purpose, we changed the parameters c 1 , c 2 , and c 3 in Eq. (54) devised for the SIDIS data, while keeping the constraint |q T |/Q < 0.2 for the DY data. We refer the reader to Ref. [7] for an analogous study of the effect of the |q T |/Q cut on the DY data.
We consider five different configurations for the cut on SIDIS data: (e) A fifth cut inspired by the PV17 analysis [5], namely the same as in the previous case but with c 1 = 0.2, c 2 = 0.7 and c 3 = 0.5; this is the least conservative choice.
In Fig. 17, we show how the global χ 2 /N dat changes when considering the five configurations described above. By observing that more conservative choices do not necessarily correspond to better χ 2 values, we conclude that the TMD formalism is able to describe SIDIS data that fails to fulfil the formal requirement |q T |/Q 1. In this respect, we notice that the global χ 2 /N dat of cut (d) is smaller than the baseline fit, despite including a larger amount of data, some of which at |q T | Q, i.e., well outside the region where TMD factorization is valid.
In Fig. 18, we better illustrate the situation by showing the comparison between COMPASS data and theoretical predictions from our baseline fit for the SIDIS multiplicity for positively charged hadrons as a function of |P hT |/Q in the bin 1.3 < Q < 1.73 GeV, 0.02 < x < 0.032, 0.3 < z < 0.4. The upper panel of the plot displays the multiplicity while the lower panel shows the ratio of experimental data to the theoretical predictions. The solid circles indicate data points included in the fit while empty squares refer to those that do not survive the cut. Remarkably, the agreement remains very good up to |P hT |/Q 0.5, well beyond the largest |P hT | allowed by the cut. We also stress that this behavior is not specific of the considered bin but is a general feature also of other bins, as well as of the HERMES dataset.
In conclusion, from our analysis it emerges that the validity of the TMD formalism in the kinematic region covered by COMPASS and HERMES seems to extend well beyond the customary cut |q T |/Q 1. This evidence justifies in a quantitative way our choice for the cut |q T |/Q in Eq. (54) for the baseline fit, and explains why we obtain values of χ 2 /N dat close to one also with less conservative cuts. Moreover, it suggests that the applicability of TMD factorization in SIDIS might be defined in terms of |P hT | rather than |q T |, calling for more extensive studies in this direction.

V. CONCLUSIONS AND OUTLOOK
In this article, we presented an extraction of unpolarized Transverse-Momentum Dependent Parton Distribution Functions and Fragmentation Functions (TMD PDFs and TMD FFs, respectively), which we refer to as MAPTMD22.
We analyzed 2031 data points collected by several experiments: 251 data points from Drell-Yan (DY) production measured at Tevatron, LHC and RHIC, 233 points from fixed-target DY (see Tab. II) and 1547 data points from Semi-Inclusive Deep Inelastic Scattering (SIDIS) measured by the HERMES and COMPASS collaborations (see Tab. III).
Our description of the experimental observables is based on TMD factorization at a perturbative accuracy defined as N 3 LL − , in the sense that we use TMD evolution at the N 3 LL level, hard factor and matching coefficients at order α 2 s , collinear PDFs at NNLO, and collinear FFs at NLO (see Tab. I). We constructed the full TMDs by combining the perturbative components, regularized by means of the b * prescription defined in Eq. (33), and nonperturbative parts described as a sum of Gaussians and weighted Gaussians (see Eqs. (38) and (39)). We used 21 free parameters: 11 related to the TMD PDF, 9 for the TMD FF, and one for the Collins-Soper kernel driving the TMD evolution. We assumed these parameters to be the same for all quark flavors.
The TMD formalism is applicable only if observed transverse momenta are much smaller than the hard scale Q of the considered process. The details of our data selection are explained in Sec. III. In the DY case, for the final lepton pair transverse momentum q T we adopted a simple criterion, |q T | < 0.2 Q, which is common to other works in the literature. In the SIDIS case, our selection criterion is more restrictive than the one made in Ref. [5], but less restrictive than the one made in Ref. [22]. For SIDIS data, where the detected hadron momentum is P hT ≈ −zq T , our cut includes many data points with |P hT | Q but also with 0.2 Q < |q T | < Q. We also tested the quality of our fit by exploring other criteria for the |q T |/Q cut. Surprisingly, we find that fits of quality comparable to the baseline result can be obtained with less conservative cuts such that |q T | Q for several data points. We believe that the definition of the range of applicability of the TMD factorization formula deserves further studies.
We found it difficult to reproduce the normalization of SIDIS data measured in fixed-target experiments at moderate to low scales. The formalism works well at NLL order, but severely underestimates the data at the N 2 LL order and above. The reason can be traced back to the fact that the integral upon transverse momentum of the TMD cross section coincides with the known collinear result only at NLL, while it misses other contributions at higher orders. A rigorous solution of the problem would imply the calculation of all these contributions (including the so-called Y -term in the language of Ref. [6]), which are currently not under control, and are beyond the scope of this work. We decided to modify our predictions by including pre-computed normalization factors that are independent of the fit, and that are identified by comparing the integral upon transverse momenta of the TMD formula to the corresponding collinear calculation, as explained in Sec. II D.
The fit is performed with the bootstrap method, leading to an ensemble of 250 replicas. We reach a very good overall agreement with data, with χ 2 /N dat = 1.06 for the central replica (defined in Sec. IV A). The description of all individual datasets is satisfactory, except for a few cases, in particular for the ATLAS data (see Tab. IV and Fig. 11). The discrepancy with the ATLAS data may be due to corrections not fully included in our analysis (see, e.g., Ref. [32]) which, in spite of being small, can have a large impact in the comparison with very high precision data.
The resulting TMDs are shown in Figs. 13 and 14. It is interesting to note that they deviate from a simple Gaussian shape (especially for the FFs) and the TMD shape changes in a nontrivial way as a function of x or z. Tables with grids of the obtained TMDs will be made publicly available at the NangaParbat website 13 and the TMDlib website 14 [120].
The availability of TMD analyses at increasing precision will be essential in guiding detector design and feasibility studies for the future Electron Ion Collider, and will also have a strong impact on precision measurements of observables particularly sensitive to the hadron structure, such as W mass measurements at hadron colliders (see, e.g., Refs. [121][122][123]).
Our work can be extended by improving the perturbative accuracy of the analysis (see, e.g., Refs. [32,124]), including theoretical uncertainties in terms of suitable scale variations (along the lines of Ref. [125]), improving the treatment of correlated uncertainties, and investigating possible flavor dependence of TMDs [15,24].