Non-perturbative contributions to vector-boson transverse momentum spectra in hadronic collisions

Experimental measurements of Drell-Yan (DY) vector-boson production are available from the Large Hadron Collider (LHC) and from lower-energy collider and fixed-target experiments. In the region of low vector-boson transverse momenta $q_T$, which is important for the extraction of the $W$-boson mass at the LHC, QCD contributions from non-perturbative Sudakov form factors and intrinsic transverse momentum distributions become relevant. We study the potential for determining such contributions from fits to LHC and lower-energy experimental data, using the framework of low-$q_T$ factorization for DY differential cross sections in terms of transverse momentum dependent (TMD) distribution functions. We investigate correlations between different sources of TMD non-perturbative effects, and correlations with collinear parton distributions. We stress the relevance of accurate DY measurements at low masses and with fine binning in transverse momentum for improved determinations of long-distance contributions to Sudakov evolution processes and TMDs.


Introduction
The production of photons, weak bosons and leptons at large momentum transfer Q Λ QCD in high-energy hadronic collisions is described successfully by factorization [1] of short-distance hard-scattering cross sections, computable at finite order in QCD perturbation theory as power series expansions in the strong coupling α s , and non-perturbative long-distance parton distribution functions (PDFs), determined from fits to experiment. It was realized long ago, however, that even for Q Λ QCD additional dynamical effects need to be taken into account to describe physical spectra in the vector-boson transverse momentum q T when the multiple-scale region q T Q is reached [2][3][4][5]. These amount to i) perturbative logarithmically-enhanced corrections in α k s ln m Q/q T (m ≤ 2k), which go beyond finite-order perturbation theory and call for summation to all orders in α s , and ii) non-perturbative contributions besides PDFs, which correspond to the intrinsic transverse momentum distributions in the initial states of the collision and to non-perturbative components of Sudakov form factors.
The summation of the logarithmically-enhanced corrections to Drell-Yan (DY) lepton pair hadroproduction has since been accomplished systematically by methods based on the CSS formalism [6]. It has been fully computed through next-to-next-to-leading-logarithmic (NNLL) accuracy, which requires calculations up to two-loop level, and partial results at three and four loops are already available for some of the coefficients needed for higher logarithmic accuracy [7,8]. On the other hand, nonperturbative effects besides PDFs in DY production are included in the formalism of transverse momentum dependent (TMD) parton distribution functions [9]. Intrinsic transverse momentum distributions enter as boundary conditions to the renormalization group evolution equations for TMDs, while nonperturbative Sudakov effects are taken into account via non-perturbative contributions to the kernel of the evolution equations associated with TMD rapidity divergences [10][11][12][13][14].
The purpose of this work is to examine the combined determination of the nonperturbative rapidity-evolution kernel and intrinsic transverse momentum k T distribution from fits to measurements of transverse momentum spectra in DY lepton-pair production at the Large Hadron Collider (LHC) and in lower-energy experiments, including Tevatron, RHIC and fixed-target experiments. To this end, we employ the calculational framework developed in [15][16][17][18][19][20]. We investigate to what extent the two sources of non-perturbative effects are correlated, and study the role of different data sets, from the high-precision DY LHC data to the lower-energy DY data, in disentangling them. We also analyze how these two non-perturbative contributions are correlated with non-perturbative contributions encoded in PDF sets. Quantifying these effects will be important both for strong interaction investigations of hadron structure and for determinations of precision electroweak parameters, as the low-q T DY region is relevant for the extraction of the W-boson mass at the LHC.
The paper is organized as follows. In Sec. 2 we briefly describe the factorization formula, evolution equations and perturbative coefficients which constitute the theoretical inputs to our analysis. In Sec. 3 we present the results of the numerical studies and fits to experimental data. We give conclusions in Sec. 4.

Theoretical inputs
We start from the TMD factorization formula for the differential cross section for DY lepton pair production where Q 2 , q T and y are the invariant mass, transverse momentum and rapidity of the lepton pair, and the TMD distributions F f ←h fulfill evolution equations in rapidity and in mass We further perform the small-b b b operator product expansion of the TMD F f ←h as follows, where f f ←h are the PDFs, C f ← f are the matching Wilson coefficients, and f NP are functions 1 to be fitted to data, encoding non-perturbative information about the intrinsic transverse momentum distributions. The non-perturbative component of the rapidity-evolution kernel D f and the distribution f NP are the main focus of this paper. The TMD distributions in Eq. (1) depend on the scales µ, ζ. To set these scales, we will use the method of the ζ-prescription proposed in [15]. (See e.g. [21] for recent examples of alternative scale-setting.) The summation of the logarithmically-enhanced corrections at low q T is achieved through Eqs. (1)-(4) by computing perturbatively the quantities H, C, γ F and Γ cusp as series expansions in powers of α s . In Table 1 we summarize the perturbative orders used for each of these quantities in the calculations that follow. We refer to the logarithmic accuracy specified by these orders as NNLL. 2 Table 1: Summary of perturbative orders used for each part of the DY cross section.
The rapidity evolution kernel D contains perturbative and nonperturbative components. The perturbative expansion for D is currently known up to three-loops [22][23][24][25]. Using the b * prescription [6], we model D as where D f res [26] is the resummed perturbative part of D f , g is an even function of b vanishing as b → 0, and with the parameter B NP to be fitted to experimental data. For the function g(b) we will use the models and fitting respectively the parameters g K , c 0 and g * K to experimental data. The quadratic model in Eq. (7) has traditionally been used since the pioneering works [27][28][29][30]. The model in Eq. (9) contains the perturbative quadratic behavior at small |b| but it goes to a constant behavior at large |b|, fulfilling the asymptotic condition ∂D/∂ ln b 2 = 0, in a similar spirit to parton saturation in the s-channel picture [31] for parton distribution functions. The model in Eq. (8) is an intermediate model between the previous two, being characterized by a linear rise at large |b|. In the following we will refer to the non-perturbative component of the rapidity-evolution kernel, modeled according to Eqs. (7)-(9), as D NP .
The nonperturbative contribution to D f in Eq. (5) also influences the rapidity scale fixing with the ζ-prescription [18]. In fact, once the nonperturbative correction is included in D f , one is to use ζ NP given by [18] Only the perturbative part ζ pert , computed in [16], was used in the fits [17]. The expression in Eq. (10) converges to ζ pert in the limit b → 0. We will use this expression in the fits of the next section.
The modeling of the TMD through the function f NP allows one to fit data at different energies. In particular it allows the nonperturbative behavior of the TMD to be described for large values of b. In [15,17,32] it has been observed that a modulation between Gaussian and exponential models is necessary. This can be provided by the following model, where the interpolation of Gaussian/exponential regimes is dependent on the Bjorken x-variable, and λ 1,..,5 > 0.

Determination of f NP and D NP from fits to experiment
We next present results of performing TMD fits to experimental data for DY differential cross sections, by employing the theoretical framework described in the previous section. We consider DY measurements both at the  Figure 1: Results of the TMD global fit to DY measurements from LHC and lower-energy experiments. PDF 1.14 HERAPDF2.0 [62] 0.97 CT14 [63] 1.59 MMHT14 [64] 1.34 PDF4LHC [65] 1.53  3 We restrict the fit to data in the low transverse momentum region by applying the cut q T /Q < 0.2 to the data sets. 4 The values of the fitted TMD parameters in Eqs. (6),(8) (for D NP ) and in Eq. (11) (for f NP ) and their associated uncertainties are shown in Fig. 1. Since PDFs enter the TMD fit through Eq. (4), the results in Fig. 1 are presented for different PDF sets. The corresponding χ 2 values are given in Table 2. We observe that the values of the fitted parameters λ i (see Eq. (11)) in Fig. 1 vary more significantly among different PDF sets than the values of the fitted parameters B NP and c 0 (see Eqs. (6), (8)), corresponding to the fact that the λ i parameters in f NP are related to the x-dependence of the distributions, while the rapidity evolution kernel is x-independent.
The correlations among TMD parameters for different PDF sets are illustrated in Fig. 2. Light colors in the pictures of Fig. 2 indicate low correlations; dark colors indicate high correlations. Shades of blue denote negative correlations; shades of brown denote positive correlations. In particular, the correlation between the parameters c 0 (controlling the long-distance behavior of the rapidity evolution kernel in Eq. (8)) and λ 1 (controlling the intrinsic transverse momentum distribution in Eq. (11)) is fairly low in the case of the HERAPDF set, but it increases in the NNPDF3.1 case, and is higher still in the CT14 and MMHT14 cases. We note that the latter two PDF sets do not include LHC data in the fits, while the NNPDF3.1 does. The χ 2 values in Table 2 are lowest for the HERAPDF and NNPDF3.1 cases. Next, we wish to focus on the role of present (and future) LHC measurements to investigate the sensitivity to the nonperturbative contributions in D NP and f NP . To this end we will perform fits to LHC data only, using a smaller number of parameters. That is, we model D NP as in Eqs. (5)-(9), depending on two parameters, B NP and either g K or c 0 or g * K , and we take a form for f NP which is simplified with respect to Eq. (11), namely, we take an x-independent simple gaussian depending on one parameter λ 1 only, which provides a measure of the intrinsic transverse momentum in terms of a gaussian width. We then perform 3-parameter fits to LHC DY data [33][34][35][36][37][38][39], fitting λ 1 , B NP and either g K or c 0 or g * K , as well as 2-parameter fits to the same data, fitting only B NP and either g K or c 0 or g * K , and fixing λ 1 to λ 1 = 0.001 GeV 2 to simulate the cases of nearly zero intrinsic transverse momentum (as in purely collinear approaches). The results from the 3-parameter and 2-parameter fits, using the PDF set NNPDF3.1, are summarized in Table 3.  We see that the 3-parameter fits (cases 2, 4 and 6 in Table 3) yield results, both for the χ 2 values and for the values of the fitted TMD parameters, which are not dissimilar from the global fit results given earlier, supporting the overall consistency of the TMD picture of low-energy and high-energy DY data. These three cases correspond to the three different long-distance behaviors of the rapidity-evolution kernel D(µ, b) in Eqs. (7)- (9). Case 2 and case 6, in particular, while giving fits of comparable quality, correspond to very different physical pictures of the nonperturbative component of D. Case 2 extends the quadratic behavior to large distance scales (see Eq. (7)). In contrast, case 6 fulfills the saturating condition ∂D/∂ ln b 2 = 0 at large |b| (see Eq. (9)). This is, to our knowledge, the first time that a full fit to low-q T DY data is performed in the hypothesis of long-distance saturating behavior of the rapidity-evolution kernel. The 2-parameter fits (cases 1, 3 and 5 in Table 3), on the other hand, show significantly different behaviors, characterized by somewhat higher χ 2 values and especially by significantly different values of the D NP fitted parameters. This indicates that, although most of the sensitivity to the intrinsic transverse momentum distribution comes from the lower-energy measurements, non-negligible f NP effects are present at the LHC too. In particular, Table 3 suggests that without any intrinsic transverse momentum distribution it may be possible to describe DY data at the LHC but this would lead to a different determination for B NP and the rapidity evolution kernel. That is, intrinsic transverse momentum effects may be reabsorbed by changes in the D NP fit.
To further analyze the sensitivity of LHC DY measurements to f NP and gain insight into the results of Table 3, we next consider the ratio where dσ T MD is the DY differential cross section computed from the full TMD fit, and dσ test is the DY differential cross section computed by setting f NP = 1 in the full fit. In Fig. 3 we plot the numerical results for the ratio (12) versus the DY lepton-pair transverse momentum q T for different values of the DY lepton-pair invariant mass Q. For reference, in Fig. 3 we also plot the theoretical uncertainty band on the full TMD result which comes from scale variation, taken according to the ζ prescription of Sec. 2. We see that in the lowest q T bins the nonperturbative effects, evaluated according to the ratio in Eq. (12), exceed the perturbative uncertainty, evaluated from scale variation in the ζ prescription. The comparison of Table 3 and Fig. 3 confirms that sensitivity to f NP is present in LHC data but may be reabsorbed by varying D NP . We explore the above point, associated with correlations between D NP and f NP , by analyzing the b dependence of the rapidity evolution kernel D(µ, b) in Fig. 4. We plot results for D from the different cases in Table 3, at µ = M Z and µ = 5 GeV. Consider first the upper right panel (µ = M Z ). The two red curves correspond to the nonperturbative quadratic model in Eq. (7). The solid red curve is the result of the 3-parameter fit in Table 3 (case 2), while the dashed red curve is the result of the 2-parameter fit in Table 3 (case 1). Similarly, the two yellow curves correspond to the saturating model in Eq. (9) (solid yellow is the 3-parameter fit, while dashed yellow is 2-parameter), and the two blue curves correspond to the linear model in Eq. (8) (solid blue is the 3-parameter fit, while dashed blue is 2-parameter).
For each of the three modeled large-distance behaviors of D(µ, b), the difference between the solid and dashed curves in the upper right panel of Fig. 4 measures the correlation between the D NP and f NP nonperturbative effects, namely, it measures the impact of the intrinsic k T on the determination of the rapidity evolution kernel. We see that in each case this impact is non-negligible. If we look at the analogous results for lower masses in the upper left panel  Table 3. In the lower panels the result for the global DY+SIDIS fit [20] is also plotted.
(µ = 5 GeV), we see that for the quadratic model particularly (red curves) the impact of intrinsic k T increases (even exceeding the uncertainty bands). That is, although the quality of the fit from the quadratic model is shown in Table 3 (case 2) to be comparable to that of the saturating and linear models, the quadratic model requires a much more pronounced dependence than the others on the intrinsic k T distribution, which is revealed especially at low masses. Apart from the intrinsic k T correlations, the differences among the three solid curves in the upper panels of Fig. 4 illustrate the current status in the determination of the large-|b| behavior of the non-perturbative rapidity evolution kernel from fits to experimental data. As expected, the sensitivity of current LHC measurements to the long-distance region is limited, which results into sizeable uncertainty bands at large |b|. This sensitivity could be enhanced by precision measurements of the low-q T DY spectrum at the LHC, with fine binning in q T , for low masses µ M Z (see e.g. first results from LHCb [66]).
For comparison, in the lower panels of Fig. 4 we also report the result for D which is obtained from the global fit to Drell-Yan and semi-inclusive deep inelastic scattering data [20] (grey curves in the two lower panels of Fig. 4). The global fit includes, besides LHC data, also data from low-energy experiments. This fit is performed assuming the linear model in Eq. (8). It is interesting to observe that the grey curves at µ = M Z and µ = 5 GeV in the lower panels, compared to the blue curves obtained from the same linear model, are lower and closer to the yellow curves (saturating behavior), reflecting the role of low-energy data in determining long-distance features of D.

Conclusion
Transverse momentum spectra in DY lepton pair production have been measured at the LHC and at lower-energy collider and fixed-target experiments. The low-q T end of DY spectra is important for the extraction of the W-boson mass and for hadron structure investigations.
In this paper we have carried out a study of low-q T DY spectra based on the TMD factorization approach in Eq. (1), using the ζ prescription (10) to treat the double scale evolution in Eqs. (2), (3). This approach contains the perturbative TMD resummation through the coefficients in Table 1 and the non-perturbative TMD contributions through f NP (intrinsic k T ) and D NP (non-perturbative Sudakov) in Eqs. (4) and (5) (besides the non-perturbative collinear PDFs in Eq. (4)). As such, it can be contrasted with other approaches in the literature: on one hand, low-energy approaches based on fixed-scale parton model which include non-perturbative TMD contributions but do not include any perturbative resummation and/or evolution of TMDs; on the other hand, high-energy approaches based on purely perturbative resummation and non-perturbative collinear PDFs, which do not include any non-perturbative TMD contributions.
We have limited ourselves to considering the low-q T region q T Q, and not addressed issues of matching with finite-order perturbative corrections which are essential to treat the region q T ∼ Q (see e.g. [52,54,56,58]).
Using this theoretical framework, we have performed fits to low-q T DY measurements from the LHC and from lower-energy experiments. The ultimate goal of these fits is to extract universal (non-perturbative) TMD distributions to be used in factorization formulas of the type (1), much in the spirit of the approaches discussed in [67][68][69]. This will be essential to bring the use of TMDs for phenomenological analyses on a similar level as that of ordinary parton distributions. The determination of non-perturbative TMDs from fits to experimental measurements is complementary to determinations from lattice QCD -see e.g. ongoing lattice studies of D NP [70,71]. In this work we have focused on studying the sensitivity of LHC and lower-energy DY experiments to non-perturbative f NP and D NP contributions, and examining their correlations with different extractions of collinear PDFs. To this end, we have defined model scenarios for D NP in Eqs. (7)-(9) and f NP in Eq. (11).
We have presented results from global DY fits (Figs. 1,2 and Table 2) and from LHC fits ( Table 3 and Fig. 3). These results indicate that, while the strongest sensitivity to the intrinsic k T is provided by the low-energy data, neglecting any intrinsic k T at the LHC worsens the description of the lowest q T bins in the DY spectrum, giving higher χ 2 values in the fit (see differences between cases 1 and 2, between cases 3 and 4, and between cases 5 and 6 in Table 3), and causes a potential bias in the determination of the rapidity evolution kernel D(µ, b) (see differences between cases 1 and 2, between cases 3 and 4, and between cases 5 and 6 in Fig. 4). A quantitative measure of the size of nonperturbative TMD effects is provided in Fig. 3 and compared with perturbative theoretical uncertainties estimated from scale variations. Given the strong reduction of these uncertainties achieved through the high logarithmic accuracy of perturbative resummations and the use of the ζ prescription for scale-setting, the residual uncertainty due to nonperturbative TMD effects is found to play a non-negligible role for the DY spectrum at the LHC in the low-q T region, which increases with decreasing DY masses.
On the other hand, we see from the comparison of cases 2, 4 and 6 in Fig. 4 that the large-|b| behavior of D is not yet constrained at present by available data both at low energy and at the LHC. We have investigated and contrasted the hypotheses of quadratic behavior, which has traditionally been considered by extrapolation from the perturbative result, and saturating behavior at long distances. We have observed in particular that the latter, besides being consistent with current LHC fits, is also compatible with the result of a global fit based on an intermediate linear model, but including low-energy DY and SIDIS data. Given the extraordinarily high experimental accuracy achieved in DY processes at the LHC, this opens new opportunities for future LHC analyses. Specifically, extending measurements of the DY transverse momentum q T , for low q T Q and with fine binning ≤ 1 GeV, into the so far unexplored region of low masses Q < 40 GeV will provide valuable new information on D at large |b|, and thus enable improved determinations of TMDs.