Non-perturbative QCD effects in $q_T$ spectra of Drell-Yan and $Z$-boson production

The factorization theorems for transverse momentum distributions of dilepton/boson production, recently formulated by Collins and Echevarria-Idilbi-Scimemi in terms of well-defined transverse momentum dependent distributions (TMDs), allows for a systematic and quantitative analysis of non-perturbative QCD effects of the cross sections involving these quantities. In this paper we perform a global fit using all current available data for Drell-Yan and $Z$-boson production at hadron colliders within this framework. The perturbative calculable pieces of our estimates are included using a complete resummation at next-to-next-to-leading-logarithmic accuracy. Performing the matching of transverse momentum distributions onto the standard collinear parton distribution functions and recalling that the corresponding matching coefficient can be partially exponentiated, we find that this exponentiated part is spin-independent and resummable. We argue that the inclusion of higher order perturbative pieces is necessary when data from lower energy scales are analyzed. We consider non-perturbative corrections both to the intrinsic nucleon structure and to the evolution kernel and find that the non-perturbative part of the TMDs could be parametrized in terms of a minimal set of parameters (namely 2-3). When all corrections are included the global fit so performed gives a $\chi^2/{\rm d.o.f.} \lesssim 1$ and a very precise prediction for vector boson production at the Large Hadron Collider (LHC).


I. INTRODUCTION
The Drell-Yan (DY) and vector boson (Z) production are fundamental processes for the Large Hadron Collider (LHC) physics as they represent a basic test for strong and electroweak interactions. In this paper we focus on the transverse momentum distributions of the leptons coming from Drell-Yan and Z-boson production when the sum of the transverse momentum of the dilepton pair is small compared to the center of mass energy and/or the invariant mass of the neutral boson. This particular regime is sensitive to non-perturbative QCD effects and has received a particular attention for a long time [1][2][3][4]. The new developments of the theory of transverse momentum dependent distributions (shortly referred to as TMDs) [5][6][7] forces us to a re-thinking of the analysis of non-perturbative QCD effects. In fact, the proper definition of the TMDs allows a systematic analysis of non-perturbative QCD in the observed transverse momentum distributions. These effects are fundamental to fix the precision of current and future experiments, including those at the LHC and at an Electron Ion Collider (EIC). The cross sections that we consider in this work can be formulated schematically according to the factorization formula [5][6][7] where F A,B are the transverse momentum dependent parton distribution functions (TMDPDFs). They depend on the dilepton invariant mass trough the scales ζ A and ζ B 1 , being ζ A ζ B = Q 4 , the intrinsic parton transverse momenta, the factorization scale µ and the Bjorken momentum fractions. Finally, H is the hard factor, which is spin independent and can be calculated adopting the standard perturbation theory. Notice that this formula can be easily extended to the case of semi-inclusive deep inelastic scattering (SIDIS) by replacing the two TMDPDFs by one TMDPDF and one TMD fragmentation function (TMDFF); analogously for electron-positron annihilation we would have the convolution of two TMDFFs 2 .
The TMDPDFs, like the ordinary collinear or integrated parton distribution functions (PDFs), are quantities which can be evolved between two different scales, so that all the non-perturbative QCD information encoded in a TMDPDF at one scale can be directly used at a different one by means of the appropriate evolutor. In impact parameter space we haveF whereR is the evolution kernel. Notice that here below we use µ 2 i,f = ζ i,f = Q 2 i,f for simplicity. In order to make the inverse Fourier transform of the TMDPDF in momentum space one has to define some prescriptions to treat the non-perturbative regime in impact parameter space. To fix the ideas we can refer to some aspects of the well-known Collins-Soper-Sterman (CSS) resummation framework [1]. In this approach the initial scale for the evolution is chosen directly in impact parameter space, namely Q i = C/b T where C is a numerical constant. This choice induces a Landau pole singularity, which is cured introducing a cutoff in impact parameter space, the well-known b * -prescription. Thus, the Landau pole is avoided by freezing the strong coupling constant in impact parameter space. However, this procedure, which is an important aspect of this approach, modifies the TMDs also in the region where their perturbative expression is supposed to hold.
In this paper we take a different path, fixing the initial scale in Eq. (2) directly as Q i ∼ Q 0 + q T , where Q 0 is a minimal energy scale at the border between perturbative and non-perturbative regimes (Q 0 ∼ 2 GeV). We find that the effect of the Landau pole is actually minimized by this choice. This confirms the result of Ref. [9], where it was argued that the completely resummed evolution kernel could be implemented in a large region of the impact parameter space. We also study the impact of a Q-dependent non-perturbative part using a specific model. We anticipate here that this correction plays a role only in the fit at next-to-next-to-leading-logarithmic (NNLL) accuracy.
In Ref. [9] the authors have deduced also the anomalous dimension of the TMDs up to next-to-next-to-leading order (NNLO), which is at the same level as the usual analysis of perturbative QCD in terms of collinear parton distribution functions. An expression of the evolutor suitable for higher order analysis is also provided in the same paper. This information is used now in this work to perform a global fit of data at NNLL, which, to our knowledge, represents a novelty in the field.
The functional form of the TMDPDF, viable both for low-energy Drell-Yan processes (say, processes with invariant masses Q ∼ a few GeV) and for high energies (say Q ∼ M Z or bigger) is still under debate. In this work we discuss and show how to construct a functional form of the TMDPDF which incorporates the full resummations with the highest available perturbative ingredients (for the moment a complete N 3 LL analysis is not possible as one misses, for instance, the four-loop cusp anomalous dimension).
A fundamental piece of information for TMDs comes from their limit when the transverse momentum is much larger than the hadronization scale (q T Λ QCD ). In this tail, the TMDPDF matches onto a Wilson coefficient and a PDF and this fact, to the best of our knowledge, is used in all models. It is well known from the literature that this matching coefficient of the TMDPDF onto the PDF can be partially exponentiated (see [10] and, more recently, [4]). Such exponentiated part is spin independent (i.e. the same for all leading-twist TMDs) and obeys a differential equation very similar to the one found for the evolution kernel. We provide an analytic solution of this equation in such a way that the exponentiated part of the matching coefficient is resummed consistently with the evolution kernel.
We find that this splitting of a TMDPDF into a Wilson coefficient and a PDF needs important non-perturbative corrections that we parametrize through a model, with a minimal set of parameters (namely two, called λ 1,2 in the text). These corrections are independent of the dilepton invariant mass Q. We have performed a study of several model forms and found that the (corrected) exponential form outlined in Section II is the best suited for the minimization of the χ 2 of our global fit. On top of this we have also explored the case of a Q-dependent non-perturbative contribution, which requires the use of a third parameter (λ 3 in the following). When this correction is included we find a general improvement of the fit, especially for the low-energy data. For this reason we think that this extra piece could play a relevant role also when studying the data from fixed-target SIDIS experiments, that are mostly collected at low Q 2 values. Notice that this kind of correction is universal, in the sense that it is the same for all TMDs.
The data that we have studied in this work come from the following experiments: E288 [11], E605 [12], R209 [13], CDF Run I [14], D0 Run I [15,16], CDF Run II [17], D0 Run II [18]. These data cover a large interval in center of mass energies, dilepton invariant masses and rapidities. The specific features (and role) of each experiment are further described below in the text.
The paper is organized as follows: in Section II we present in some detail the theoretical approach, discussing the construction of the TMDPDFs and the main aspects of the cross sections that we are going to study. We also comment on some features of the data under consideration. In Section III we show our phenomenological analysis with a selection of the most interesting results, in Section IV we comment of previous works on the subject and then, in Section V, we draw our conclusions.

II. THEORETICAL FRAMEWORK
As stated above, the processes which are relevant in this work are DY and Z-boson production, A B → γ * → + − , A B → Z → + − , where A and B are the incoming hadrons whose center of mass energy is denoted as √ s in the following. In the case of DY, say γ * production, one can write the differential cross section as where, adopting the standard notation, For the case of Z-boson production the analogous observable is where The approximate value of τ γ,Z holds only when Q 2 q 2 T , which is a good approximation for all sets of data we have used. The term Y γ,Z is the usual remnant piece used to recover the complete perturbative limit of the cross section at high transverse momenta. For the hard coefficient we use H(Q 2 , µ 2 ) = |C V (Q 2 , µ 2 )| 2 , where C V is the Wilson coefficient that can be extracted from the finite terms of the calculation of the (full QCD) quark form factor in pure dimensional regularization and is known up to three loops [19,20]. At one-loop we have: Since, as shown in the equation above, the hard coefficient has an explicit dependence on −Q 2 , in this work we have used the π-resummed coefficient as suggested in Ref. [21]. This resummation ensures a slightly better convergence of the perturbative expansion in the hard coefficient, although the effects are not numerically relevant as for the Higgs boson production. In the following the factorization scale µ is fixed at the same value as the invariant mass of the exchanged vector boson, Q or M Z , so that the TMDPDFs are completely evolved up to this energy scale and the logarithms in the hard part are minimized. For some sets of data one has to integrate on the rapidity y within the range |y| ≤ y i max = − 1 2 lnτ i , i = γ, Z. In order to finally build the TMDPDFs we need to include the completely resummed evolution kernel and discuss the operator product expansion of the TMDPDFs onto the PDFs. This is done in the following sections.

A. Evolution kernel
The basic quantities which encode all QCD non-perturbative information are the evolved TMDPDFs. The evolution of a generic quark-TMDPDF is given by Eq. (2), where the evolution kernelR is [9] In this equation γ F is the anomalous dimension of the TMDPDF and is related to the anomalous dimension of the hard factor as explained in Ref. [6]. The function D can be resummed to obtain the D R function, that will replace D in Eq. (8), via the renormalization group equation where Γ cusp is the cusp anomalous dimension. Solving this equation one finds and re-writing α s (µ b ) as a series in α s (µ), where µ b = 2e −γ E /b T and β(α s ) is the QCD β-function [9], one finds where we have used the notation The final form of the full resummed kernel is The resummation of the evolution kernel is valid only up to a certain maximum value of the impact parameter [8,9], beyond which the evolution kernel becomes completely non-perturbative. Given that the evolution kernel can be derived from the soft matrix element [8], which is process and spin independent, the non-perturbative correction to the evolution kernel is also process and spin independent. This correction to the D R term provides then a correction to the evolution kernel depending on ln(ζ f /ζ i ), as can be deduced from Eq. (8). In this work we have tried to fix this non-perturbative part phenomenologically, once the complete perturbative resummation is performed. The complete form of the evolution kernel is theñ where D NP is the non-perturbative piece of the D function.

B. Operator Product Expansion of the TMDPDF onto a PDF, spin-independent exponentiation and resummation
In this section we recall the perturbative limit of the TMDPDFs for high transverse momenta, building the TMD-PDF consistently.
The main constraint on the structure of TMDPDFs comes from their limit at high transverse momentum, q T Λ QCD . In this case the TMDPDFs split into a Wilson coefficient and a PDF, see the specific literature for its proof [5,6]. Taking into account this perturbative tail we start writing the TMDPDF for a generic nucleon N as The coefficientsC Q / for the quark-quark and the quark-gluon channel at one loop are [6] C Q / and the one-loop splitting kernels appearing above are given by The important point is that under renormalization group equation one has where the double logarithms can be partially exponentiated in a way similar to what presented in Ref. [10] (see also Ref. [22]):C where The coefficients of the perturbative expansions of Γ cusp and γ V can be found in the Appendix C of Ref. [9]. Choosing h Γ,γ (b T ; µ b ) = 0 we have, at fixed order, and correspondinglyĈ Thus, the partially resummed TMDPDF can be written as The expression above still contains large logarithms L T that need to be resummed when α s L T is of order 1. This resummation is also necessary to have an expression consistent with the fully resummed evolution kernel derived in the previous Section. Solving the evolution equations for D, Eq. (9), and for h Γ and h γ , Eq. (21), we find respectively Eq. (10) and The method of integration is the same used in Ref. [9] for the evaluation of D R , and it can be immediately applied to find h R γ as the equation for h γ (second line of Eq. (25)) has the same functional form as Eq. (10), so that In the case of h Γ we note that L T = αs(µ) αs(µ b ) dα β(α) , and in this way Eq. (25) (first line) can be solved as One can perform the above integration first expanding the β-function and finally re-writing α s (µ b ) in terms of α s (µ) at the correct order, as shown in Ref. [9]. The result is Joining all pieces together we get the fully resummed perturbative part of the TMDPDF: An important point to underline in this formula is that the whole factor in front of the integral is universal among all leading-twist TMDPDFs. In fact, D R is a piece of the universal evolution kernel [8], and h R Γ,γ are deduced respectively from the cusp and the TMDPDF anomalous dimensions, which are independent of the specific TMDPDF. This fact shows that the treatment of this factor is not restricted to the unpolarized TMDPDF, but can be applied as well to all (leading-twist) polarized TMDPDFs. To the best of our knowledge this observation has never been used for the explicit construction of transverse momentum dependent distributions.
An explicit representation of Eq. (29) is given in Fig. 1 for a fixed scale µ 2 = ζ = Q 2 i . We notice that even for such a low-energy scale as Q i ∼ 2 GeV, the region of high values of the impact parameter is suppressed. A further suppression is found when the non-perturbative input for the TMDPDF is considered, as detailed in the next Section. We notice that the change between the next-to-leading-log (NLL) and NNLL curves is driven by the matching coefficientĈ q←j , whose perturbative order in each curve is detailed in Table I. In this table we report also the orders in α s included in each step of the logarithmic resummation. Comparing the two plots in Fig. 1 we can also observe that at high values of the impact parameter, where we expect bigger non-perturbative effects, the values of the TMDPDF are basically independent of Q i . On the other hand for higher values of Q i and low values of b T the TMDPDF is particularly enhanced.

C. Scale choices and non-perturbative inputs for the TMDPDF
Despite the resummations performed above the formulas so obtained are still of perturbative origin, and as such they have to be corrected by non-perturbative contributions. The full resummations allow nevertheless an improvement of the convergence of the perturbative series on a large portion of the impact parameter space. Below we definitely switch to the notation: Parametrizing the non-perturbative large-b T region of the TMDPDF, we write it at some initial scale Q i as Notice that in the equation above we have parametrized the non-perturbative contribution in the same way as we did for the evolution kernel in Eq. (14). We observe thatF pert q/N (x, b T ; Q 2 i , µ i ) goes to zero for high values of b T at fixed Q i , in agreement with our expectations, see Fig. 1, and the non-perturbative contributions are not expected to alter this behavior.
Another relevant issue in Eq. (30) concerns the choice of the scales Q i and µ i . In order to minimize the value of the logarithms we choose µ i = Q i . Next we notice that the splitting into a coefficient and a PDF is valid only at high transverse momentum, so that we expect that the choice Q i = Q 0 + q T (where Q 0 is a fixed low scale) minimizes the logarithms generated by this splitting. The scale Q 0 works as a minimum matching scale between the TMDPDF and the PDF, such that it sits at the border between the perturbative and non-perturbative regimes; in particular we choose Q 0 ∼ 2 GeV.
As explained below Eq. (7), in the cross section we fix the factorization scale µ = Q, so that collecting all results in Eqs. (2), (14) and (30) we can write the resummed TMDPDF that enters into the factorization theorem as We point out that in the last factor in Eq. (32),F NP q/N (x, b T ; Q), we have included the non-perturbative Q dependence coming from the evolution kernel (Eq. (14)), More explicitly, the TMDPDF is implemented as where, as we already mentioned, In order to fix the arguments of the non-perturbative partF NP , we need to consider the following constraints: • It must correct the behavior ofF pert q/N at large values of b T , where the perturbative expansion looses its convergence properties and the Landau pole singularity shows up, both in the evolution kernel and in the matching coefficient of the TMDPDF onto the PDF.
• It has to be such that in order to guarantee that the perturbative series is not altered where its convergence properties are sound.
We have not included a dependence on x, as data eventually do not need such correction and to keep the model simple enough. In Eq. (35) we are assuming that the values of x are not extremely small (say x > 10 −3 ), in which case the whole TMD formalism should be re-considered.
We have studied several parametrizations of the non-perturbative part (Gaussian, polynomial, etc.) and the final one which better provides a good fit of the data, with the minimum set of parameters and D NP = 0, is As discussed below in the text the data for Z-boson production are basically sensitive just to the parameter λ 1 , that is to the exponential factor and not to the second power-like term that, controlling the large b T region, is more sensitive to small q T data. The global fit so performed allows to fix, to a certain precision, the value of this non-perturbative constant. In other words, this fit can be used to fix the amount of non-perturbative QCD corrections in the transverse momentum spectra. As commented above, the parameter λ 2 corrects the behavior of the TMDPDF at high values of b T and results necessary to describe the data at low dilepton invariant mass and low q T . Considering now a nonzero D NP , this results in a Q-dependent factor in the non-perturbative model (see the studies of Refs. [23,24] and more recently Refs. [5,8]). Thus, from Eqs. (33) and (36), by setting D NP = λ 3 b 2 T /2, we havẽ We anticipate here that the sensitivity of the data to this extra factor with λ 3 is not very strong, although we observe an improvement in the χ 2 . This is a consequence of the fact that the fully resummed D function is actually valid on a region of impact parameter space which is broad enough for the analysis of the sets of available data (notice that we have in all cases a dilepton invariant mass Q > 4 GeV). It might be that at lower values of Q such corrections could be more significant. On the other hand one expects that also the factorization theorem should be revised when the values of Q become of the order of the hadronization scale. It is then possible that the non-perturbative corrections to the evolution kernel happen there where the basic hypothesis of the factorization theorem (Q q T ∼ Λ QCD ∼ O(1 GeV)) become weaker and so are more difficult to extract. A more detailed study in this direction is beyond the scope of this paper.
Finally in Fig. 2 we show the effect of the model of Eq. (36) on the TMDPDF at low scale, Q = 2 GeV, where we expect that its impact is more substantial. We see that the non-perturbative correction affects mainly the intermediate b T region, while the high-b T region keeps naturally suppressed.
As a general remark one has to keep in mind that in practical calculation we finally Fourier transform the product of two TMDPDFs. The integration in impact parameter space is done numerically over a suitable b T range. We have checked that the region outside the endpoints of this integration does not affect the final result. In fact, the points for very small b T are relevant only for extremely high transverse momenta, which is not the case of our work. At very high b T the TMD are completely negligible (see Fig. 2).
The results for the global fit are detailed in the next sections.

III. PHENOMENOLOGICAL ANALYSIS
We now move to a comparison of our theoretical estimates based on the above approach to available measurements on unpolarized cross sections. We perform a fit on different sets of data at various energies and covering a large q T range. This will allow us to fix the parameters entering the non-perturbative model, test the scale evolution of TMDs and the validity of the whole approach and highlight its main features. We will also discuss and quantify the uncertainties on the parameters so extracted, coming both from the statistical analysis in the fitting procedure as well as from some theoretical assumptions of our calculations. We will start with the selection of data discussing their role in our fit and then present our results.

A. Data
The selection of experimental data used in our fit is an important issue. Definitely we need data at moderate center of mass energies covering the small-q T region (up to 1-2 GeV) and intermediate dilepton invariant mass values (below 10 GeV). These come mainly from fixed-target experiments. On the other hand to access larger q T values and even larger scales we have to include also high-energy collider experimental data, like those from Tevatron at the Z 0 pick. In both cases we will keep fulfilling the requirement q T Q, region of application of our approach. Notice that to conform with the standard notation adopted in experimental analysis in the following we use M = Q for the dilepton invariant mass.
These two classes of data are indeed complementary and essential to test the scale evolution of TMDs over a suitable range of scale values and to quantify the role of the non-perturbative part entering these distributions.
More precisely we consider the following data sets (see also Tabs. II and III, where we collect further details): • moderate energy p-Cu and pp data, 0 < q T < 1-2 GeV, 4 GeV < M < 25 GeV: -R209 data [13] at √ s = 62 GeV; • high-energy pp data, 0 < q T < 20 GeV, M = M Z : -CDF Run I [14] at √ s = 1.8 TeV; -D0 Run I [15,16]   While for the low-energy data we consider the invariant differential cross section in the virtual boson momentum, for the high-energy data sets we use their ratio of the q T dilepton distribution normalized to the experimental total cross section. In such a case, we compute this numerator following our approach and use the normalization factor as obtained with the DYNNLO code of Catani et al. [27,28]. The use of this ratio avoids the problem of the discrepancy between D0 and CDF experiments (see Tab. III) that could cause a source of systematics and/or tension between data sets.
As a first attempt of our global fit we include all data shown in Tables II and III. From this we realize that the E605 data set gives a very large, and anomalous, contribution to the total χ 2 . We then decide to study separately and in some detail the role played by these data. This is the outcome of this dedicated analysis: • a separate fit of these data alone gives a large χ 2 by itself (χ 2 points 10). The bad overall χ 2 is not simply a matter of tension between these data and the others; • the different invariant mass bins seem not to be consistent among them and not to follow a proper scale dependence (even considering a Q-dependent non-perturbative part in the TMD); • the parameters resulting from the global fit are almost unaffected by the inclusion of these data; • in the literature the inclusion of this data set seems also problematic. In Ref. [2] only the two bins at the lowest dilepton invariant masses are included, without any comment on the rest of data.
For these reasons we prefer to exclude the E605 data set in the following global analysis. We will further comment on the role of these data when comparing our study with other fits, Section IV.

B. Fitting procedure and results
We perform different fits to test various aspects of our approach. We start by considering only the high-energy data set from Tevatron [14][15][16][17][18], adopting for the non-perturbative factor,F NP in Eq. (36), the exponential part alone without the correcting power-like behaviour. This results in a very good fit (χ 2 d.o.f. = 0.44) showing how the model is well sound and that a one-parameter fit is enough to describe the Z-boson spectrum. We also tried a Gaussian functional form -that implies a Gaussian form also in momentum space (as commonly adopted in many parametrizations of TMDs) -but the resulting χ 2 is worse.
The use of the exponential form, instead of the more usual Gaussian-like form, has been suggested in Refs. [29] and [30] based on the fact that Euclidean correlation functions in quantum field theory are usually exponential and not Gaussian at large distances. To our knowledge it is the first time that this hypothesis has been tested explicitly on a phenomenological analysis and within the TMD formalism.
This preliminary study allows us to draw some interesting conclusions: • the large-q T spectrum (q T above the pick, say q T > ∼ 5 GeV) at Tevatron is very well reproduced by the perturbative, and therefore calculable, part of the TMDs. This is a non trivial result showing the consistency of the approach as well as the resummation procedure.
• The non-perturbative piece is necessary to describe the data points below and around the pick, with a strong preference in favor of an exponential damping instead of a Gaussian shape in b T space. As the data are at fixed dilepton invariant mass (M Z ) one cannot check the Q-dependence of the non-perturbative model considering this set of data alone.
• Being the low-q T region populated by very few data points, to which the non-perturbative part should be more sensitive, our fit at NNLL appears somehow over-determined with a χ 2 d.o.f.

Notice that at NLL accuracy this is not the case, showing once again the importance of a NNLL analysis.
As a second step we include also lower-energies DY data and perform a global fit. This fit, as already pointed out, covers the very important small-q T region and therefore is expected to strongly constrain the non-perturbative piece. This is what we will discuss in the following.
Differently from the Tevatron data, the low-energy data are affected by large uncertainties in the normalization of the cross section. For this reason we allow for two extra parameters in the fit, namely the normalization factors for E288 and R209 data: N E288 and N R209 . Notice that these extra two parameters are strongly related to the present accuracy of the available experimental data, while the two (three) parameters entering our TMDPDF, λ 1 , λ 2 (λ 3 ) are the main goal of our analysis.
The global fit of Tevatron and low-energy data adopting a simple exponential factor in the non-perturbative model is very bad and the description of the Tevatron data results somehow spoiled. This result does not change if we allow a Q dependence of the non-perturbative model, alike in Eq. (37). This is a signal that the low-q T region, so sensitive to the large-b T behaviour of the TMD, requires more care. For this reason we exploit various simple extensions of the functional form ofF NP and end up with the power-like piece (quadratic in b T ) entering Eqs. (36), (37): this is what better describes the data.
We perform a fit both at NLL and NNLL. When adopting the NLL approximation we use the next-to-leading-order (NLO) collinear parton distributions, while at NNLL we use the NNLO PDFs. In both cases we adopt Q i = Q 0 + q T and for the collinear PDFs we use the MSTW08 set [31]. We also check the role played by the collinear PDFs adopting another set, namely the CTEQ10 [32] and we find a complete consistency among the results.
One of the main goals of this work consists in the fits performed at NNLL accuracy with full resummation. The NLL fits are mainly used as a check of convergence of the theory and other phenomenological aspects. We have tested both a Q-independent and Q-dependent parametrization of the non-perturbative inputs as given respectively in Eq.  Tables IV and V. The NNLL-NNLO calculation gives a better overall χ 2 and definitely improves the description of Tevatron data. Concerning the low-energy data, while the E288 fit is practically unaffected, for R209 we have a strong correlation between the normalization factor and the relative χ 2 that eventually gives a good description in both approximations. By looking at Tab. IV we see that for each high-energy data set, that means also large-q T and large-Q values, the χ 2 for data points improves significantly going from NLL to NNLL accuracy. This is somehow expected since in this region the logarithmic terms are more important. On the other hand we slightly loose in two low-energy data sets, still achieving a better overall χ 2 .  [11,13] and Tevatron data [14][15][16][17][18], with D NP = 0 (Eq. (36)), Qi = Q0 + qT , at NNLL and NNL accuracies and with the collinear parton distributions from MSTW08 [31] at NNLO and NLO.
The normalization factor for E288 data is always within their experimental uncertainty, while for R209 it is a bit larger. One observes however that the central values for R209 data are not so accurate since there are no official tables from the Collaboration, so that one can expect a larger normalization error for these data. The parameter λ 2 is unaffected by the higher order contributions being the same in both approximation, while λ 1 presents some differences even if within the relative errors. The technical reason for this shift is the appearance at NNLL of the one-loop contribution of the coefficientsĈ as outlined in Tab. I and visible also in Fig. 1. We expect that higher order contributions on this coefficients would stabilize the result.
The statistical error is estimated requiring a 68% confidence level, corresponding to a ∆χ 2 = 4.72 for four parameters. Concerning the theoretical errors, we study the dependence on the initial scale Q i = Q 0 + q T (where Q 0 = 2 GeV) in two ways: i) we check the impact of a change in Q 0 allowing m charm ∼ 1.3 GeV ≤ Q 0 ≤ 2.7 GeV, where the lowest value of Q 0 is about the charm threshold and the highest value is limited by the energy of the lowest energy bin of data; ii) keeping Q 0 = 2 GeV, we vary Q 0 + q T /2 ≤ Q i ≤ min (Q 0 + 2q T , Q). In the first case the fit is practically unaffected concerning the values of the parameters λ 1 , λ 2 . For the second case the scale dependence instead has some impact on the these values. In particular from Tab. V one can see that at NLL the theoretical error is of the same order of the statistical one and that there is a clear reduction of the scale dependence at NNLL. At this order the main uncertainty on the fitted parameters comes from the statistical error. In Fig. 3 we show the comparison of our theoretical estimates with Tevatron at NLL and NNLL accuracies, while in Figs. 4 and 5 we show the corresponding comparison for low-energy data. In Fig. 6 we present the impact of the statistical error on the TMDPDF at NNLL at Q = M Z and Q = 5 GeV.  E288 pN E l = 400 GeV, √s=27.4 GeV, y=0.03 5<M<6 GeV 6<M<7 GeV 7<M<8 GeV 8<M<9 GeV 11<M<12 GeV 12<M<13 GeV 13<M<14 GeV   FIG. 4. Comparison of our theoretical estimates for d 3 σ/d 3 q with E288 at three different energies [11]. The results are obtained from the global fit with D NP = 0 (Eq. (36)), Qi = Q0 + qT , at NLL accuracy with collinear PDFs at NLO (upper panels) and at NNLL accuracy with NNLO PDFs (lower panels). For the collinear PDFs we use the MSTW08 set [31]. T with R209 data [13]. The results are obtained from the global fit with D NP = 0 (Eq. (36)), Qi = Q0 + qT , at NLL accuracy with collinear PDFs at NLO (left panel) and at NNLL accuracy with NNLO PDFs (right panel). For the collinear PDFs we use the MSTW08 set [31]. The study of the Q dependence of the fit within the model of Eq. (37) is summarized in Tables VI and VII. It is important to note that the Q-dependent term cannot substitute any other piece of the fit. In other words putting λ 1 = 0 and/or λ 2 = 0 and leaving only the term with λ 3 = 0 in the model of Eq. (37) does not produce a reasonable fit. The estimates of the error in the the fit and the parameters are done using the same criterium (68% confidence level) as in the case λ 3 = 0.
With a three-parameter fit we have an improvement of the total χ 2 , specially at NNLL. In fact, at NLL the uncertainties somehow mask the benefits of the introduction of a new correction and we do not obtain a significative change of the χ 2 . The core of the improvement in the χ 2 is that all low-energy experiments are much better described at NNLL, while in the case of λ 3 = 0 the change from NLL to NNLL is more controversial for this data subset. In other words, while the introduction of a new parameter in the fit is not conclusive from a NLL analysis, it seems instead more appropriate when all perturbative pieces are developed at NNLL. The final χ 2 /d.o.f. ∼ 0.8 at NNLL confirms that this kind of precision is necessary to extract important information on the non-perturbative structure of TMDs, that a NLL analysis may hide.
Comparing the values of λ 1,2 with and without the inclusion of the Q-dependent correction we observe that λ 1 is stable and practically not affected by the introduction of λ 3 . On the other hand the central value of λ 2 manifests a sensible change, which is evidently due to the fact that both λ 2 and λ 3 cooperate to improve the description of the low-energy data. Looking at Tab. VI one sees also that the Tevatron data (χ 2 /points 1) are probably overparametrized in a NNLL fit with λ 3 = 0. This is due to the fact that the non-perturbative effects for Z-boson production are probably important only for a small fraction of the data at low transverse momentum. Also this fact is not evident when looking just at the NLL fit. Notice that the range of q T explored in this work is basically the same as the one considered by different authors (see, for instance, Ref. [2] [11,13] and Tevatron data [14][15][16][17][18] with D NP = 0 (Eq. (37)), Qi = Q0 + qT , at NNLL and NNL accuracies and with the collinear parton distributions from MSTW08 [31] at NNLO and NLO.   [11,13] and Tevatron data [14][15][16][17][18], with D NP = 0 (Eq. (37)), Qi = Q0 + qT , at NNLL and NNL accuracies and with the collinear parton distributions from MSTW08 [31] at NNLO and NLO.
For completeness we show in Fig. 7 the best fit curves at NNLL for this case. In the case of Tevatron data we find that the peak region is slightly enhanced providing so a better description of the data.
We point out that the parameter λ 3 parametrizes a non-perturbative correction to the evolution kernel and as such, it is spin independent and the same for TMDPDFs and TMDFFs.. Thus we expect that this correction should be included also in the analysis of SIDIS data. Further study in this direction is left for future work.
As a first application and test of this phenomenological analysis, we use the present TMD formalism and our estimates at NNLL-NNLO accuracy to predict the q T dependence for Z-boson production at CMS (here we use the model with λ 3 = 0). Our results, shown in Fig. 8, are compared with the data from Ref. [33] and we observe the pick region is nicely described. The band which appears in the figure shows the error due to the statistical uncertainty on λ 1 , which is the main source of error according to our analysis. The extremely good agreement (χ 2 /points ∼ 0.9) supports once again the overall picture.

IV. COMPARISON WITH PREVIOUS WORKS
One of the goals of this paper is to provide a framework for the analysis of transverse momentum distributions in which all ingredients coming from perturbation theory are under control and used to their maximum extent. Only in such a case the parametrization of the non-perturbative inputs can be reliably treated. Although such an attempt is not included in most of the studies available in the literature, different intermediate steps have been discussed in several works. In the following we focus on the most relevant ones from our point of view.
A detailed phenomenological analysis of low-and moderate-energy DY data, aimed at extracting the transverse momentum dependence of TMDs at leading-order accuracy, appeared in Ref. [34]. The relevance of this study was the attempt to describe within a unified, even if simplified, picture the role of TMDPDFs and TMDFFs in different  36)), Qi = Q0 + qT , at NNLL-NNLO accuracy compared to CMS experimental data [33]. The band comes from the statistical error on the fitted parameter λ1.
hadronic processes. In Ref. [3] the authors present an analysis of transverse momentum distributions of vector boson production at NNLL. The needed Fourier transforms to momentum space are done deforming the integration contour in b-space and calculating moments. They consider only very high dilepton invariant mass and explicitly avoid a complete treatment of non-perturbative effects. In this respect, within the TMD formalism that we have developed, one can achieve a direct comparison of the non-perturbative inputs with low-energy data. The evolution of TMDs allows in this way a complete fixing of this non-perturbative part, and so of the precision of the theoretical prediction.
The authors of Ref. [4] perform also an analysis of Z-boson production. Although in their formalism they do not consider the theory for TMDs, the expression for the cross section agrees with ours when looking at the DY case or Z-boson production at high transverse momentum (of course they do not claim any universal structure which could eventually be used in SIDIS). The resummations provided in our work are different from the ones in Ref. [4] in the sense that in their "collinear anomaly" part they perform a sum of logarithms which is valid up to values of the strong coupling and impact parameter such that α s L 2 T ∼ 1. Notice that one can re-obtain the "collinear anomaly" factor re-expanding D R , h R Γ,γ in α s and counting α s L 2 T ∼ 1. Given that this is not the highest possible resummation that one can perform, the Landau pole does not appear explicitly in their resummed expression (although the perturbative series has intrinsically a Landau pole problem). The authors in any case realize that some non-perturbative input is necessary and they suggest some Gaussian non-perturbative (Q-independent) part in impact parameter space, without performing any fit of Z-boson production data. A non-perturbative correction to the "collinear anomaly" factor is suggested in [35]. In the present work we perform a complete resummation of the logs with the counting α s L T ∼ 1, which is more relevant when low-energy data are included, and we check that an exponential non-perturbative (Qindependent) correction in impact parameter space works better than the Gaussian one for the Z-production data. A more complex analysis is then done here to describe also the low-energy DY data. We nevertheless agree with the authors of Ref. [4] that non-perturbative parameters cannot be fixed looking just at the Z-production case and that resummations allow for a better control of the scale dependence of the final result.
The most comprehensive work aimed at considering both low-energy and high-energy data and including an explicit non-perturbative parametrization of the cross section is the one of Refs. [2] and [36] (the so called BLNY model) which implements a model within the standard CSS approach. Their analysis considers the cusp anomalous dimension at order α 2 s and the rest of the coefficients at order α 1 s , and was done before the recent developments of the TMD formalism appeared. They have no partial exponentiation of the matching coefficient between the TMDPDF and the PDF. In their approach they use the well known b * prescription, where b * = b T / 1 + b 2 T /b 2 max and b max is a parameter, modeling the Q-dependent non-perturbative corrections in terms of the parameter g 2 , and the Qindependent corrections in terms of g 1 , g 3 . Consequently, as shown in Ref. [36], all the parameters g 1,2,3 heavily depend on the choice of b max . In that work they find as a best value b max ∼ 1.5 GeV −1 , consistent with the low-b T behavior of D R as observed in Ref. [9]. Notice that the parameter g 2 is the same for DY, SIDIS and e + e − processes and for all leading twist TMDs, while g 1,3 are specific of the unpolarized TMDPDF. An important point is that the DY data from which these parameters are extracted have Q ≥ 4 GeV, while SIDIS data (say from HERMES or COMPASS experiments) on which these parameters are often used, cover lower Q 2 ranges. Notice that in BLNY the model of the evolution kernel (which is Q-dependent) plays a fundamental role and is crucial in the fit of data.
The point of view of our work is instead completely different. In fact we have shown that in order to analyze the Drell-Yan data (so for dilepton invariant mass Q ≥ 4 GeV) the Q-dependent non-perturbative corrections play a minor role and most of the non-perturbative effects are Q-independent. In other words, we have found that the extraction of the parameters b max and g 2 in the BLNY model is done using data for which the evolution kernel can instead be almost entirely predicted theoretically. The differences in our fit between the cases with λ 3 = 0 and λ 3 = 0 are important but not essential for the success of the fit. In any case, the final χ 2 that we get with λ 3 = 0 improves notably the results of BLNY. A second issue concerns the fact that the parameters (b max , g 2 ) extracted from DY data are used at scales different from those of their extraction. It is possible that using these values in low energy SIDIS data (which usually have 1 GeV < Q < 4 GeV) can cause some tension. While this problem could be avoided with our parametrization of the TMDs, it needs a dedicated study which is beyond the scope of this paper and we leave for future work. Finally for a closer comparison to the fit of Ref. [2], where they include also the two bins of E605 with the lowest dilepton mass, we have also checked the impact of these two bins on our global fit with λ 3 = 0 at NNLL. The number of points of the two bins is 14, and so the impact of these points on our global fit is relatively small. The central values and the errors of λ 1,2,3 do not change while the final χ 2 /d.o.f. = 0.95, confirming our considerations in Section III A.
Several modifications of the BLNY model have been proposed in the literature and sometimes used in actual fits. In Ref. [37] the authors propose a different parametrization for the g 2 term, based on phenomenological considerations of SIDIS data. Again, we can apply the same considerations as in the original BLNY model: while we do not exclude a non-perturbative structure of the evolution kernel, its relevance cannot be clearly stated doing a NLL fit given the current status of DY data.
In Ref. [38] a first attempt of a global fit of SIDIS and DY data is performed, by using the b * prescription and resumming large logarithms at NLL accuracy. The fit of data is only qualitatively correct and the final χ 2 is not declared. In the present work we have shown the importance of including higher perturbative orders and resummations in the parametrization of the TMDPDFs, together with other ingredients.
Recently a fit of DY data has been performed in Ref. [39]. The authors propose a modification of the BLNY model, changing both the Q-dependent part (i.e. the evolution kernel), and the Q-independent part (inserting new arbitrary parameters). Concerning the fit of the DY data they fit the three standard parameters of the BLNY model plus the normalization of all experimental data sets they use. The fit is done considering only a subset of low energy data. Moreover the partial χ 2 's of each experiment are not shown and the total χ 2 for the DY data declared is about 1.5. When they consider SIDIS data they increase the number of parameters, namely the ones for the TMD fragmentation function, and a normalization factor for each z-bin. Given the amount of variables in the fit, the almost arbitrary choice of data and the final χ 2 they get, their result cannot be considered conclusive.

V. CONCLUSIONS
The TMD formalism is a powerful tool to analyze perturbative and non-perturbative effects in q T spectra. In this work we have fitted the DY and Z-boson production data to fix the non-perturbative part of TMDPDFs: all results are collected in Tables V and VII. We have stressed the fact that in order to have a reliable fixing of the nonperturbative inputs one has to provide a fully resummed expression for the perturbative part. The fully resummed cross section in fact is less sensitive to the factorization scale dependence and this allows a more stable extraction of the non-perturbative pieces of the TMDs. We have provided a full study of the perturbative inputs. In particular we have used the TMD evolution kernel at NNLL [9] which, to our knowledge, was never used before in a global fit of this kind. We have also discussed the matching of TMDPDFs onto PDFs with the exponentiation and fully resummation of the corresponding coefficient. We have argued that the exponentiated part of this matching coefficient is spin independent and should be included in the analysis of other types of TMDPDFs. This part is fundamental to have a reliable description of the TMDs both at NLL and NNLL.
One of the important aspects of our perturbative analysis is that the factorization scale is fixed in momentum space instead of the more usual impact parameter space. This choice provides a good stability of the perturbative series and offers a new understanding of the data and of the model dependence of the TMDs. We find that the NNLL fit clarifies several issues about the non-perturbative nature of TMDs.
The model-dependent non-perturbative inputs for the TMDPDF are studied in order to minimize the number of non-perturbative parameters and to provide a good description of the data. We find that the Z-boson data are better described by an exponential damping factor in impact parameter space rather than a Gaussian one. The associated parameter, called λ 1 in the text, has a stable value within the errors, which are mainly of statistical origin. The lowenergy data explore values of the impact parameter higher than those covered in the case of Z-boson production. We find that a polynomial correction with a new parameter, called λ 2 in the text, plays a relevant role in this respect and both corrections, induced by λ 1,2 , do not depend on the dilepton invariant mass Q. The values of these parameters can be fixed by fitting data for DY and Z-boson production and the NNLL resummation greatly reduces the theoretical error on this determination. Particular attention has been paid to the study of the Q dependence of the non-perturbative model. The insertion of this contribution (parametrized by λ 3 ) provides only an improvement of the χ 2 at the price of adding a new parameter to the fit. Nevertheless, given the actual uncertainties on data and collinear PDFs, the need for this correction in the fits cannot be firmly established. Increasing the precision of the experimental data can be crucial to fix this issue.
One important outcome of this work is that a fully resummed evolution kernel, together with other exponentiations and resummations in the various matching coefficients that appear in the cross section, avoid an excessive use of a modelization of the cross section, making the predictions more stable.
We consider this work as a first step towards the proper understanding of non-perturbative effects in transverse momentum distributions. Several important perturbative pieces, recently calculated [40][41][42][43], can be used in an approximate N 3 LL analysis and will be included in a forthcoming publication.
We point out that fixing the non-perturbative part of transverse momentum distributions can improve substantially the theoretical precision needed for the current LHC experiments, as our prediction for Z-boson q T spectrum at CMS shows.
Finally, we comment on the use of this formalism for SIDIS processes. The parameters λ 1,2 are specific of the unpolarized TMDPDF and can be used also in the analysis of SIDIS data. A different value of these parameters is expected for the fragmentation function. The parameter λ 3 is a universal correction and, as such, it is the same in DY and SIDIS processes. One important aspect that must be pointed out is that most of SIDIS data at our disposal are given for 1 GeV < Q < 4 GeV, that is in an energy range lower than the one studied in this work. For this reason, it can be that the model with D NP = 0 (maybe different from the one discussed in this work) is better suited for the description of the SIDIS data in this regime and we expect that a NNLL study would be fundamental to understand them.