The small k T region in Drell-Yan production at next-to-leading order with the Parton Branching Method

The Parton Branching (PB) method describes the evolution of transverse momentum dependent (TMD) parton distributions, covering all kinematic regions from small to large transverse momenta k T . The small k T -region is very sensitive both to the contribution of the intrinsic motion of partons (intrinsic k T ) and to the resummation of soft gluons taken into account by the PB TMD evolution equations. We study the role of soft-gluon emissions in TMD as well as integrated parton distributions. We perform a detailed investigation of the PB TMD methodology at next-to-leading order (NLO) in Drell-Yan (DY) production for low transverse momenta. We present the extraction of the nonperturbative “intrinsic-k T ” distribution from recent measurements of DY transverse momentum distributions at the LHC across a wide range in DY masses, including a detailed treatment of statistical, correlated and uncorrelated uncertainties. We comment on the (in)dependence of intrinsic transverse momentum on DY mass and center-of-mass energy, and on the comparison with other approaches.


Introduction
The measurement of the vector boson transverse momentum, p T , in Drell-Yan (DY) production [1] allows one to investigate in detail many different aspects of the strong interaction sector of the Standard Model, and their impact on precision electroweak measurements.The very low p T region of the DY cross section is sensitive to the contribution from the nonperturbative transverse motion of partons inside the hadrons; additionally at low transverse momentum multiple soft gluon emissions have to be resummed; at larger transverse momenta perturbative higher-order contributions become dominant.The precise description of the Z/γ boson transverse momentum distribution has been investigated since the 1980's, and approaches like CSS [2] analytic resummation and parton-shower [3][4][5][6] numerical algorithms have been applied with different success.
In this work we explore the approach [7,8] to DY p T spectra based on the parton branching (PB) TMD methodology in momentum space proposed in [9,10], and we perform a detailed analysis of the small-p T region for wide ranges in center of mass energies and in DY masses.Though fitted only on deep-inelastic scattering (DIS) data from HERA experiments, the PB-TMD methodology has been shown to be capable of describing DY p T spectra at LHC energies [7] and at low energies [8] without any need for adjustment of parameters.This approach takes into account simultaneously soft gluon radiations and the transverse momentum recoils in the parton branchings along the QCD cascade.It provides a successful natural treatment of the multiple-scale problem of the DY transverse momentum for transverse momenta much smaller than DY masses but also of the DY with hard jet production [11].It also confirms the universality of the TMDs being able to describe both DIS and DY cross sections at all available center of mass energies [12].Alternative approaches based on parton showers in standard Monte Carlo event generators like PYTHIA8 [3] can also describe multi-differential DY cross section but it has been observed that they require intrinsic transverse momentum distributions strongly dependent on √ s [13,14].In order to describe the measurements at LHC energies, a Gaussian width exceeding the Fermi motion kinematics is needed.Approaches based on CSS [15] provide very precise analytic predictions for inclusive enough observables like the Drell-Yan cross section transverse momentum.In this paper we study in detail the low k T behavior of the PB-TMD parton distributions where both very soft gluon emission and intrinsic-k T contribute significantly and interplay.The results presented here provide a multi-scale economical and coherent approach demonstrating the sensitivity to nonperturbative TMD contributions and first steps in disentangling the intrinsic-k T contribution from the nonperturbative Sudakov one [16].We compare DY theoretical predictions with experimental measurements in wide ranges in center-of-mass energies, √ s and in DY masses, m DY , to extract the intrinsic-k T parameter from the transverse momentum distributions.We are carefully taking into account systematic and statistical uncertainties using the breakdown of experimental uncertainties provided by the full set of covariance matrices available in the recent Drell-Yan differential cross section measurement at 13 TeV [17] and we treat for the first time the scale uncertainties in the theoretical predictions as correlated uncertainties within a given mass bin.
The results for TMD parameters such as intrinsic-k T obtained from the DY analysis in this paper can be compared with analogous results obtained from TMD fits in the CSS coordinatespace framework, see e.g. the recent studies [18,19].A significant difference between these approaches and the approach of this paper concerns the treatment of collinear parton distribution functions (PDFs).As shown in Refs.[7,8], in the approach of this paper the inclusive DGLAP limit is recovered and fits of collinear distributions are made, e.g. from inclusive DIS structure functions, along with TMD distributions [20][21][22].In contrast, CSS approaches do not recover inclusive DGLAP and rather use an ansatz based on the operator product expansion of TMD distributions in terms of collinear PDFs, assuming collinear PDFs to be given by standard PDF sets.The PDF bias effect [19] which results from this has been shown to influence significantly the central values of the extracted distributions and dominate the systematic uncertainties in all the existing TMD determinations based on CSS approaches.The possibility to treat collinear and TMD distributions on the same footing and determine them without having to rely on existing PDF fits is a distinctive feature of the PB TMD approach.We believe that in the long run this could bring significant advantages in pursuing TMD phenomenology.
On the other hand, the results of this paper for intrinsic-k T can also be compared with the case of parton shower Monte Carlo event generators, such as PYTHIA [3] and HER-WIG [4].Monte Carlo tuning to experimental data shows that parton shower approaches require intrinsic-k T distributions dependent on the center-of-mass energy √ s [13,14], and a Gaussian width exceeding the Fermi motion kinematics.In contrast, in the approach of this paper we find that the width of the intrinsic-k T distribution has a much milder center-ofmass energy √ s dependence.We obtain more natural Gaussian width σ, σ = q s / √ 2, with q s close to 1 GeV resulting from fits to DY measurements from fixed-target to LHC energies.We propose in this paper that the different behavior, concerning intrinsic-k T distributions, between PB TMD and parton-shower approaches can be ascribed to the different treatment of the contributions to parton evolution from the nonperturbative Sudakov region, near the soft-gluon resolution boundary.See also [23] for a discussion of this and comparison of PB TMD and parton-shower results.
The paper is organized as follows.In Sec. 2 we briefly recall the basic elements of the calculational framework [7,9,10,[21][22][23]: we start with the PB TMD approach; next we give a few comments on the treatment of the small transverse momentum region in this approach; then we discuss the Monte Carlo computation of DY differential distributions.Sec. 3 is the central section of the paper, in which we perform fits to DY data and present results for the intrinsic-k T TMD parameter.We give conclusions in Sec. 4.

PB TMDs and DY production
To study the different contributions to the low-p T spectrum, at different m DY and different √ s, we calculate DY production cross section in the PB TMD method, which proceeds as described in Refs.[7,23].NLO hard-scattering matrix elements are obtained from the MADGRAPH5_AMC@NLO [24] at next-to-leading (NLO) event generator and matched with TMD parton distributions and showers obtained from PB evolution [9,10,20], using the subtractive matching procedure proposed in [7] and further analyzed in [25].
We will show that the application of PB TMD distributions leads to a non negligible contribution of pure intrinsic-k T , even if most of the small-k T contribution comes from the PB-evolution.We also show that the proper treatment of photon radiation from the DY decay leptons is rather important, especially in the DY mass region below the Z boson peak.The contribution of intrinsic-k T of heavy flavor partons is found to be negligible over the whole range since heavy quarks are not present in the initial configuration of the proton.

TMD distributions from the PB method
The PB evolution equations for TMD parton distributions A a (x, k, µ 2 ) of flavor a are given by [9] where k and q are 2-dimensional momentum vectors, z M is the soft resolution scale [10], z is the longitudinal momentum transferred at the branching, P (R) ab (α s , z) are the resolvable splitting functions * (whose explicit expressions for all flavor channels are given in [9]), and ∆ a are the Sudakov form factors The branching evolution (1) fulfills soft-gluon angular ordering [27][28][29], with the branching variable q ′2 being related to the transverse momentum q T of the parton emitted at the branching by It is shown in [10] that angular ordering is essential for the TMD distribution arising from the solution of Eq. ( 1) to be well-defined and independent of the choice of the soft-gluon resolution scale z M = 1 − ε for ε → 0. In contrast, p T ordering leads, for instance, to ambiguities in the definition of the TMD from the z → 1 region.Analogously to the case of ordinary (collinear) parton distribution functions, the distribution A a (x, k, µ 2 0 ) at the starting scale µ 0 of the evolution, in the first term on the right hand side of Eq. ( 1), is a nonperturbative boundary condition to the evolution equation, and is to be determined from experimental data.For simplicity we parameterize with the width of the Gaussian distribution given by σ = q s / √ 2, independent of parton flavor and x, where q s is the intrinsic-k T parameter.
On the other hand, it is found that PB-NLO-2018 Set2 provides a much better description, compared to PB-NLO-2018 Set1, of measured Z/γ transverse momentum spectra at the LHC [7], in low-energy experiments [23], and of di-jet azimuthal correlations near the back-to-back region at the LHC [40].This underlines the relevance of the angular-ordered coupling α s (q 2 T ) in regions dominated by soft-gluon emissions.
Based on this observation, in the following we will focus on the PB-NLO-2018 Set2 approach and perform fits to DY transverse momentum measurements to investigate the sensitivity of these measurements to the nonperturbative TMD intrinsic-k T parameter q s , and perform determinations of its value.
As discussed in [7,20], in order to complete the definition of the PB-NLO-2018 Set2 scenario the treatment of the coupling α s needs to be specified in the region of small transverse momenta q T < ∼ q 0 , where q 0 is a semi-hard scale on the order of a GeV.As in [7,20], we take setting q 0 = 1 GeV, which may be regarded as similar in spirit to the "pre-confinement" proposal in the context of infrared-sensitive QCD processes [41,42]  † .In the present study, we will perform a determination of the nonperturbative TMD parameter q s from DY transverse spectra by assuming the above behavior for α s .
To better illustrate the underlying physical picture, we give next a few further comments on nonperturbative contributions and the treatment of the small transverse momentum region in the PB TMD approach.
As implied by Eqs. ( 1) and ( 2), the PB TMD method incorporates Sudakov evolution via phase space integrations of appropriate kernels over the resolvable region, i.e. over momentum transfers z up to the soft-gluon resolution scale z M .For each branching evolution scale q ′2 , it is instructive to examine separately parton emissions with transverse momenta above the semi-hard scale q 0 , q T > q 0 , and below q 0 , Λ QCD < q T ∼ < q 0 .Using the angular ordering relation (3), these emissions are mapped respectively on the regions where z dyn = 1 − q 0 /|q ′ | is the dynamical resolution scale associated with the angular ordering [27][28][29]34].In region (a), the strong coupling ( 6) is evaluated at the scale of the emitted transverse momentum, α s (q 2 T ); the contribution from region (a) to the evolution in Eqs. ( 1), (2) corresponds to the perturbative Sudakov resummation (see e.g.[43,44]).In region (b), the strong coupling (6) freezes around the semi-hard scale q 0 ; the contribution from region (b) to the evolution is the nonperturbative Sudakov form factor in the PB TMD approach.
It is worth noting that the PB-NLO-2018 Set 2 framework provides a very natural and economical description of nonperturbative Sudakov effects, based on perturbative modeling of the Sudakov form factor (2) combined with the infrared α s behavior (6): it does not contain any additional nonperturbative functions and parameters, besides the scale q 0 .
In Fig. 1 we show parton distributions obtained with the PB approach using the starting distributions from PB-NLO-2018 Set2.We show distributions for the gluon and down quark parton densities for different values of z M : z M → 1 (default -regions (a + b) -red curve [20])  and z M = z dyn = 1 − q 0 /q ′ (region (a) only -blue curve obtained with the same parameters as PB-NLO-2018 Set2 except z M using UPDFEVOLV [45]).The distributions obtained from PB-NLO-2018 set2 with z M → 1 are significantly different from those applying z M = z dyn , illustrating the importance of soft contributions even for collinear distributions.In Ref. [46] it was found that limiting the z-integration leads to inconsistencies.In Ref. [47] a procedure to correct the z limitation is discussed.A detailed discussion on the role of soft gluons and the nonperturbative Sudakov form factor is given in Ref. [48].Please note that the intrinsick T distribution, since not part of the collinear calculation, does not affect the collinear parton densities.
In the transverse momentum distributions obtained with the PB-approach, the effect of the z M cut-off is even more visible.In Fig. 2 the transverse momentum distributions obtained for down and charm quarks are shown for PB-NLO-2018 Set2, with z M → 1, i.e. regions (a + b), with (red curve) and without intrinsic-k T distribution applied (blue curve -a Gauss distribution with q s = 0.00001 GeV).We also show the transverse momentum distribution contribution from region (a) alone, i.e. for z M = z dyn = 1 − q 0 /µ ′ , without intrinsic-k T (corresponding to the magenta curve of Fig. 1).The importance of the large z-region on the transverse momentum distributions is seen in the comparison with the predictions without intrinsic-k T distribution (blue and magenta curves).
The transverse momentum distributions show very clearly the large effect of the choice of z M for the soft region, while in the perturbative region k T > q 0 the effect becomes smaller with increasing k T .Applying such a scale, z M = z dyn = 1 − q 0 /µ ′ , removes emissions with q T < q 0 (there are still low-k T contributions, which come from adding vectorially all intermediate emissions).However, very soft emissions are automatically included with z M → 1.
As shown in Fig. 2, the effect of the intrinsic-k T distribution is much reduced at large scales, but the contribution of the region z dyn < z < 1 stays important for small k T .
It is interesting to observe that the charm density shows essentially no effect of an intrinsick T distribution: this is because charm is generated dynamically from gluons only, and there is no intrinsic charm density.

Transverse momentum distributions of PB-NLO-2018
After having discussed the importance of the soft nonperturbative region to the transverse momentum distribution, we turn now to a discussion of the transverse component of the PB parton distributions of Ref. [20], which are used for comparison with measurements.
In previous investigations on Z-boson production at the LHC [7], as well as for low DY mass, m DY , and at low √ s [8], it was found that PB-NLO-2018 Set2 describes the measurements much better, while PB-NLO-2018 Set 1 gives too large a cross section at small DY lepton pair transverse momenta, p T (ℓℓ).
The difference between PB-NLO-2018 Set1 and Set2, which comes from the choice of renormalization scale (argument in α s ), is seen essentially in the low-k T region, where the nonperturbative Sudakov form factor (region (b)), with the integral z M → 1, plays an important role.In Fig. 3 (upper row) the distributions for up and charm quarks are shown when no intrinsic-k T distribution is included, the lower row shows distributions including the default intrinsic Gauss k T distributions of widths q s = 0.5 GeV.It is very interesting to observe that and µ = 100 GeV (right column) obtained from the PB approach based on PB-NLO-2018 Set2.Two distributions do not include intrinsic-k T : the blue curve corresponds to z M → 1 (regions a+b in text) and the magenta curve to z M = z dyn = 1 − q 0 /q (region a only).The red curve is the published one PB-NLO-2018 Set2 [20] and including intrinsic-k T and z M → 1.
the differences between the sets setting q s = 0 or not are very much reduced for heavy flavors since they are only generated dynamically (since heavy flavors are not present at the starting scale in the VFNS which is applied here).In principle an intrinsic charm contribution can be included in PB densities, however, this is not required from inclusive DIS data [35] used in  [20] as a function of k T at µ = 100 GeV and x = 0.01.In the upper row are shown distributions when no intrinsic-k T distribution is included (q s = 0.00001 GeV), and the lower row shows the default distributions with q s = 0.5 GeV.
the fit of PB-NLO-2018 Set2.In Fig. 4 the transverse momentum distribution for down quarks, with and without an intrinsic-k T distribution, is shown at different scales µ.While at low scales µ ∼ 50 GeV a significant effect of the intrinsic-k T distribution is observed for very small k T , at large scales µ ∼ 350 GeV this effect is much reduced.This scale dependence will result in a much smaller sensitivity to the intrinsic-k T distribution at high m DY .The lower panels show the full uncertainty of the TMD PDFs, as obtained from the fits [20].Shown is the ratio to each central value.The red band shows the uncertainty of PB-NLO-2018 Set2, the blue line has no uncertainty band.

Calculation of the DY cross section
The cross section of DY production is calculated at NLO with MADGRAPH5_AMC@NLO [24].In the MCatNLO method, the collinear and soft contributions of the NLO cross section are subtracted, as they will be later included when parton shower, or as in our case, TMD parton densities are applied.As in earlier studies, we use CASCADE3 [49] to include TMD parton distributions and parton shower to the MCatNLO calculation (a detailed investigation of the effect of TMD parton distributions and parton showers applied in the CASCADE3 Monte Carlo generator is given in Ref. [25]).We use the HERWIG6 subtraction terms in MCatNLO, since they are based on the same angular ordering conditions as the PBTMD parton distribution sets, PB-NLO-2018 Set1 and PB-NLO-2018 Set2, described in the previous section.The validity and consistency using HERWIG6 subtraction terms in MCatNLO together with PB TMD distributions has been studied in detail in the appendix of Ref. [25].The predicted cross sections (labeled as MCatNLO+CAS3 in the following) are calculated using the integrated versions of the NLO parton densities PB-NLO-2018 Set1 and PB-NLO-2018 Set2 together with α s (m Z ) = 0.118 at NLO.
The factorization scale µ, used in the calculation of the hard process is set to µ = , with the sum running over all final state particles, in case of DY production over all decay leptons and the final jet.For the generation of transverse momentum according to the PB-TMD distributions, the factorisation scale µ in the hard process is set to µ = m DY , in the case of a real emission it is set to µ = 1 2 i m 2 i + p 2 t,i .The generated transverse momentum is limited by the matching scale µ m =SCALUP [49].Since there are no PB-fragmentation functions available yet, the final state parton shower in CASCADE3 is generated from PYTHIA [50], including photon radiation of the lepton pair.
A good description of the final state QED corrections, and in particular the kinematic effect of the real photon radiations, is essential in order to achieve a precise description of the DY transverse momentum.Fig. 5 (left) shows the DY mass distribution as measured by CMS [51] at 13 TeV together with predictions of MCatNLO+CAS3 ‡ .The bands show the scale uncertainty coming from a variation of the renormalization and factorization scale by a factor of two up and down, avoiding the extreme values (7-point variation).The DY mass is calculated from the so-called dressed-leptons (see for example [53,54]), where photons radiated within a cone of radius of R < 0.1 are merged to the lepton before the momenta are calculated.We show predictions based on PB-NLO-2018 Set1 and Set2, and also, for illustration, when photon radiation is turned off in the final-state shower (labeled as "no-QED").A rather good description of the DY mass spectrum over a large range on m DY is obtained both with PB-NLO-2018 Set1 and Set2.Only at m DY greater than a few hundred GeV the predictions tend to become smaller than the measurement (while still within the uncertainties).However, this is the region where the partonic x becomes large and not well constrained by the fit to HERA data [35] used for the PB-NLO-2018 TMD extraction [20].In the region of m DY below the Z-pole, one can observe the importance of QED corrections.In Fig. 5 (right) we show the photon transverse momentum spectrum in Z-production as measured by CMS [55] at 7 TeV in comparison with MCatNLO+CAS3 including QED radiation.The photon spectrum is well described at low E T < 40 GeV, while the high E T spectrum predicted by the parton shower falls below the measurement, since the precision of parton showers are limited for the high p T region.

The transverse momentum spectrum of DY lepton pairs
The transverse momentum spectrum of DY pairs at √ s = 13 TeV has been measured for a wide m DY range by CMS [17].We use this measurement for comparison with predictions of MCatNLO+CAS3 based on PB-NLO-2018 Set1 and PB-NLO-2018 Set2, as shown in Fig. 6.As already observed in previous investigations [7,23,25,40], the PB-NLO-2018 Set1 gives too high a contribution at small transverse momenta p T (ℓℓ), while PB-NLO-2018 Set2 describes the measurements rather well, without any further adjustment of parameters § , underlining the role of evaluating the strong coupling at the transverse momentum scale.In order to illustrate the importance of QED corrections, we show in addition a prediction based on PB-NLO-2018 Set2 without including QED final state radiation (labeled noQED).Especially in the low m DY region, the inclusion of QED radiation is essential, not only changing the total cross section but rather strongly modifying the shape of the transverse momentum distribution p T (ℓℓ).All calculations predict too low a cross section at large transverse momentum ‡ We use the Rivet package [52] for the calculation of the final distributions.§ The predictions shown here are slightly different compared to the predictions in [17] because we use here a lower minimum k T cut and because of a bug in the treatment of QED radiation in Rivet, corrected in version 3.  due to missing higher-order contributions in the matrix element.In Refs.[11,56,57] it is shown explicitly that including higher orders in the matrix element through the TMD multijet merging technique gives an excellent description even for largest p T (ℓℓ).For all further distributions, we restrict the investigations to p T (ℓℓ) below the peak region (i.e.p T (ℓℓ) < ∼ 8 GeV).

Influence of the intrinsic-k T distribution on DY transverse momentum distributions
Given the rather successful description of the DY p T (ℓℓ)-spectrum with MCatNLO+CAS3 using PB-NLO-2018 Set 2 in the low p T (ℓℓ)-region, we investigate below the importance of the intrinsic-k T distribution.In PB-NLO-2018 the intrinsic-k T distribution is parameterized as a Gauss distribution with zero mean and a width σ 2 = q 2 s /2 [20] (see Eq. ( 4)), where q s was fixed by default at q s = 0.5 GeV.
In order to illustrate the sensitivity range of the intrinsic-k T distribution, we show in Fig. 7 the MCatNLO+CAS3 predictions for the low p T (ℓℓ)-spectrum of DY production at different DY masses m DY for different intrinsic-k T distribution (with different q s parameter values) compared to the CMS measurement [17].We observe that sensitivity to intrinsic-k T is more pronounced at small p T (ℓℓ) values.This sensitivity decreases with increasing mass, as expected from Fig.In the following we describe a determination of the Gaussian width q s for different DY masses, m DY , at different √ s.The prediction is obtained from a calculation of MCatNLO+CAS3 using TMD distributions obtained with the PB-NLO-2018 Set2 parameters for the collinear distribution, but with different q s values.We scan for each m DY -bin q s in steps of 0.1 to 0.3 GeV in the range q s = 0.1, . . ., 2.0 GeV.At higher DY transverse momenta, higher order contributions have to be taken into account (a study using multijet merging is given in Refs.[11,56,57]).

Fit of the Gauss width q s in pp at √ s = 13 TeV
The transverse momentum distribution of DY leptons has been measured by the CMS collaboration [17].This is the basic measurement for the determination of the intrinsic-k T parameter q s , since it covers a wide m DY -range with high precision and that a detailed uncertainty breakdown, discussed in subsection 3.2.1, is provided.The measurements of Z-production obtained from LHCb [58] are discussed in subsection 3.2.2, while measurements at lower center-of-mass energies are shown in subsection 3.3.

DY production over a wide DY mass range
The CMS collaboration has measured Drell-Yan production at 13 TeV [17] covering a range of DY mass m DY = [50,76,106,170,350,1000] GeV.The measurement is provided with a detailed uncertainty breakdown, corresponding to a complete treatment of experimental uncertainties including correlations between bins of the measurement for each uncertainty source separately.Note that we use the fully detailed breakdown of the experimental uncertainties provided on the CMS public website ¶ .¶ The corresponding HEPdata records only contain summarised information  In order to determine the intrinsic-k T we vary the q s parameter and calculate a χ 2 to quantify the model agreement with the measurement.We evaluate the following expression || , with m i being the measurement and µ i being the prediction for data point i.The covariance matrix C ik is decomposed into a component describing the uncertainty in the measurement, C meas.

ik
, and the statistical and scale uncertainties in the prediction, || The code used with the full covariance matrix is available in Ref. [59], an earlier version to be used directly with Rivet is in Ref. [60].The covariance matrix of the measurement is taken directly from the supplementary material provided by CMS.The statistical uncertainty in the prediction, arising from the use of a Monte Carlo simulation, is accounted for as a small diagonal contribution without correlations between bins, where σ i,stat.is the bin-by-bin statistical uncertainty.We also treat, for the first time, the scale uncertainties of the theoretical prediction as a correlated uncertainty, for a given m DY range, allowing for a global shift of all bins together within the band defined by the symmetrized envelope of the scale uncertainties.This contribution to the covariance matrix is constructed as follows, where σ i,scale encodes the scale variation for each bin.We first extract independent values of q s for each invariant mass region considered in the measurement, considering only the region most sensitive to q s , p T (ℓℓ) < 8 GeV.We reduce this range further in the first two regions to p T (ℓℓ) < 6 GeV for 50 < m DY < 76 GeV and p T (ℓℓ) < 7 GeV for 76 < m DY < 106 GeV to stay in the region of sensitivity and not be biased by missing higher orders in the predictions affecting high p T (ℓℓ) shape.The obtained χ 2 /n.d.f (reduced χ 2 )) values are shown in Fig. 8 [17].The "data" uncertainty is the one estimated using min(χ 2 ) + 1, the "scan" uncertainty accounts for the step size of the q s scan, and the "bins" uncertainty is derived by varying the number of bins included in the fit.The number of bins used in the fit gives n.d.f.The values of q s in each m DY -bin as obtained from Ref. [17].Indicated is also the combined fit value of q s .Within each region, we consider the value of q s for which the smallest χ 2 is obtained as our "best fit" value.We construct a one-sigma confidence region as the set of all q s values for which χ 2 (q s ) < min(χ 2 ) + 1.When possible, this region is determined graphically using a linear interpolation between scan points.When the minimum is too narrow for a reliable determination of the uncertainty using this method, we use instead a quadratic interpolation between the lowest three points and add an uncertainty equal to one half-bin-width (0.05 GeV) in quadrature.In addition, we include an uncertainty derived by repeating the procedure with modified fit boundaries.The values obtained using this method are listed in Table 1 and a comparison is shown in Fig. 9.
The values derived from each m DY interval are compatible with each other.The most precise determination is obtained from the Z peak region, 76 < m DY < 106 GeV, followed by the regions around it.The sensitivity at high mass suffers from larger statistical uncertainties in the measurement.This independence of the intrinsic-k T with the DY mass contrasts with the need to tune the Parton Shower parameters for different masses in standard Monte Carlo events generators (see [17] -Fig.6 for a data comparison with MADGRAPH5_AMC@NLO interfaced with PYTHIA Parton Shower).
Having obtained compatible results, we proceed to deriving a combined fit by calculating a joint χ 2 including the considered bins in all mass ranges.For this, we construct a new covariance matrix C comb.ik as a sum over the 650 uncertainty sources included in the detailed breakdown.We consider that each systematic uncertainty is fully correlated between m DY regions and construct their covariance matrices in the same way as in Eq. (11).The statistical uncertainties (data and Monte Carlo) in the measurement feature nontrivial correlations due to the use of unfolding but are independent in each m DY region, and therefore we construct a block-diagonal matrix from the covariance matrices in each m DY region.The statistical uncertainty in the prediction is diagonal.We consider that the uncertainties in the QCD scales are not correlated between m DY regions and use a block-diagonal matrix.
The χ 2 values obtained using the combined covariance matrix are shown in Fig. 8.The best fit value, extracted in the same way as for separate regions, is, q s = 1.04 ± 0.03(data) ± 0.05(scan) ± 0.05(binning) GeV.This value and its uncertainty are shown as a black line and shaded area on Fig. 9 for comparison with the individual m DY bins.A cross-check has been performed by interpolating the prediction for each bin between q s values and searching for the minimum of the χ 2 distribution using a finer q s grid.It returned values within the uncertainties quoted above.The TMD distributions including the new q s value are available in TMDlib and TMDplotter [38,39].
To make consistency checks of the obtained value of q s and to examine possible trends of its dependence on DY mass and centre-of mass energy, the DY measurements at high rapidity and lower collision energies have been analysed.Since for these measurements no full error breakdown are available, we treat all uncertainties as being uncorrelated and do not include systematic uncertainty coming from the scale variation in the theoretical calculation.

Z production at high rapidities at 13 TeV
The LHCb collaboration [58] has measured Z-production at √ s = 13 TeV in the forward region, covering a rapdity range of 2 < |y| < 4.5.The χ 2 distribution is shown in Fig. 10 summed over the rapidity range of the DY lepton pair as a function of q s .A minimum is obtained for q s = 0.74 ± 0.15 GeV, where the uncertainty comes from a variation of χ 2 by one unit and from the step size of the q s scan.The reduced χ 2 /n.d.f distribution as a function of q s summed over all rapidity regions obtained from a comparison of the MCatNLO+CAS3 prediction with the measurement by LHCb [58].The shaded area corresponds to χ 2 + 1.The best fit value is q s = 0.74 ± 0.15 GeV.The value of q s = 1.04 ± 0.08 GeV as obtained from the measurements in Ref. [17] is indicated by a black vertical line.

The Gauss widths q s from lower center of mass energies
The ATLAS collaboration has measured the production of DY from collisions at √ s = 8 TeV in several DY mass bins, out of which only the two with 44 < m DY < 66 GeV and 66 < m DY < 116 GeV are relevant for p T (ℓℓ) < 10 GeV [61].In Fig. 11 we show the χ 2 /n.d.f as a function of q s obtained from these two mass bins (n.d.f = 8).
The Tevatron experiments D0 [62] and CDF have measured transverse momenta of DY lepton pairs created in pp collisions at lower center-of-mass energies (1.8 TeV [63] and 1.96 TeV [64]).The PHENIX collaboration measured DY production at √ s = 200 GeV [65], and E605 [66] at √ s = 38.8GeV.The Drell-Yan differential cross section in p T (ℓℓ) has also been measured in pPb data at √ s = 8.1 TeV by CMS [67]. Figure 11 shows the impact that the q s choice has on χ 2 /n.d.f for these different measurements.

Consistency between determinations of intrinsic k T width
A global fit of q s is obtained by calculating χ 2 for different measurements, as shown in Table 2, including the corresponding center-of-mass energies, collision types and the number of fitted data points, resulting in a total of 81 data points.The impact of intrinsic-k T distribution at lower collision energies has been analyzed using the entire range of p T (ℓℓ), while at higher center-of-mass energies we investigate up to the peak region in the transverse momentum distribution.  .
The χ 2 /n.d.f distribution as a function of q s , for all the data together, is shown in Fig. 12.The χ 2 distribution exhibits a minimum at around q s = 1.0 GeV, which is consistent with the value obtained as described above.
Figure 13 [17,58,[61][62][63][64][65][66][67].The minimum of global DY data fit is close to q s = 1 GeV and consistent with the CMS measurement [17] shown separately by a black line.ferent measurements in Refs.[17,58,[61][62][63][64][65][66][67].For the data which do not provide detailed uncertainty breakdown and are mainly used for the cross checks and comparison purpose, the uncertainty bars of q s shown in the figures are obtained from the χ 2 variation of one unit and step size of the q s scan only.The value of q s = 1.04 ± 0.08 GeV , as derived from the measurements in Ref. [17], is compatible for all ranges of m DY , and also holds true for various values of √ s.The obtained value is also found to be compatible for pPb data.Refs.[17,58,[61][62][63][64][65][66][67].Right : same as a function of √ s.The value of q s = 1.04 ± 0.08 GeV as obtained from the measurements in Ref. [17] is indicated.
To summarize, we have obtained a value for the width of the Gauss distribution for modeling the intrinsic-k T distribution inside protons of q s = 1.04 ± 0.08 GeV.This value, in contrast to standard Monte Carlo event generators, has no strong dependence on the centerof-mass energy as well as on the mass of the produced Drell-Yan lepton pair m DY .The results of this section indicate that the treatment of soft emissions in the region z dyn ∼ < z < z M with the strong coupling of Eq. ( 6) applied in PB-NLO-2018 Set2 leads to intrinsic-k T distributions with width parameter q s consistent with Fermi motion kinematics, and mildly varying with energy.

Conclusion
In this paper we have carried out a detailed application of the PB-TMD methodology, which is reviewed in the first part of the paper, and used it to describe the DY low transverse momentum distributions across a wide range of DY masses.Within this methodology, we have presented the extraction of the intrinsic-k T nonperturbative TMD parameter from fits to the measurements of DY p T differential cross sections performed recently at the LHC at √ s = 13 TeV, for DY masses between 50 GeV and 1 TeV.We have compared this with extractions from other DY measurements at different center-of-mass energies and masses.As shown previously, the measured DY cross section at low p T favours a choice of the strong coupling α s scale to be taken as the transverse momentum of each parton emission, as in angular-ordered CMW parton cascades.This corresponds to the TMD parton distribution set PB-NLO-2018 Set2.In this paper we use PB-NLO-2018 Set2 with "pre-confinement" scale q 0 of 1 GeV.The strong coupling is evaluated at the emitted transverse momentum q T for

Figure 1 :
Figure 1: Integrated gluon and down-quark distributions at µ = 4 GeV (left column) and µ = 100 GeV (right column) obtained from the PB approach based on PB-NLO-2018 Set2.The red curve is the published PB-NLO-2018 Set2[20] and corresponds to z M → 1 (regions (a + b) in text)..The blue curve corresponds to z M = z dyn with q 0 = 1 GeV (region (a) only).The ratio plots show the ratios to the one for z M → 1.

TMDplotter 2 . 2 . 4 Figure 3 :
Figure 3: TMD parton density distributions for down and charm quarks of the published PB-NLO-2018 Set1 (red curve) and PB-NLO-2018 Set2 (blue curve)[20] as a function of k T at µ = 100 GeV and x = 0.01.In the upper row are shown distributions when no intrinsic-k T distribution is included (q s = 0.00001 GeV), and the lower row shows the default distributions with q s = 0.5 GeV.

Figure 4 :
Figure 4: TMD parton density distributions for down quarks of PB-NLO-2018 Set2 with (red curve) and without (blue curve) intrinsic-k T distribution as a function of k T at different scales µ and x = 0.01.The lower panels show the full uncertainty of the TMD PDFs, as obtained from the fits[20].Shown is the ratio to each central value.The red band shows the uncertainty of PB-NLO-2018 Set2, the blue line has no uncertainty band.

1 CMS, 13 CMS, 13 1 CMS, 13 Figure 6 :
Figure6: The p T (ℓℓ) dependent DY cross section for different m DY regions as measured by CMS[17] compared to MCatNLO+CAS3 predictions based on PB-NLO-2018 Set 1 (red curve) and Set 2 (blue curve).Also shown are predictions without the inclusion of final state QED radiation from the leptons (green curve).The band shows the 7-point variation of the renormalization and factorization scales.

Figure 7 :
Figure 7: Drell-Yan cross section ratios of MCatNLO+CAS3 predictions for different q s values over CMS measurement [17] as a function of p T (ℓℓ) for different m DY regions.Only the lowest p T (ℓℓ) values are shown.The points error bar show the statistical uncertainties and the gray bands the total experimental uncertainties.The gray area at the highest p T (ℓℓ) values show the maximal values included in the fit described in section 3.2.

Figure 8 :
Figure8: The reduced χ 2 /n.d.f distribution as a function of q s for different m DY regions obtained from a comparison of the MCatNLO+CAS3 prediction with the measurement by CMS[17].The points represent the obtained χ 2 values.The lines represent the curves used for the uncertainty estimate (see text), which is materialized by the shaded areas.

Figure 9 :
Figure9: The values of q s in each m DY -bin as obtained from Ref.[17].Indicated is also the combined fit value of q s .

Figure 10 :
Figure 10: The reduced χ 2 /n.d.f distribution as a function of q s summed over all rapidity regions

Figure 13 :
Figure13: Left: the value of q s as a function of the DY-mass as obtained from the measurements in

Table 1 :
as a function of q s .Results of the fit on individual m DY intervals for the CMS measurement

Table 2 :
displays the value of q s as a function of m DY and √ s obtained from the dif-All data sets with the corresponding center-of-mass energies, collision types and the number of degrees of freedom used for the global fit of q s .