1 Introduction

The measurement of the vector boson transverse momentum, \(p_{\textrm{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/\gamma \) 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_{\textrm{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_{\textrm{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_{\textrm{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 \(\sqrt{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_{\textrm{T}}\) behavior of the PB-TMD parton distributions where both very soft gluon emission and intrinsic-\(k_{\textrm{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_{\textrm{T}}\) contribution from the nonperturbative Sudakov one [16]. We compare DY theoretical predictions with experimental measurements in wide ranges in center-of-mass energies, \(\sqrt{s}\) and in DY masses, \(m_{\text {DY}}\) , to extract the intrinsic-\(k_{\textrm{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_{\textrm{T}}\) obtained from the DY analysis in this paper can be compared with analogous results obtained from TMD fits in the CSS coordinate-space 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_{\textrm{T}}\) can also be compared with the case of parton shower Monte Carlo event generators, such as Pythia [3] and Herwig [4]. Monte Carlo tuning to experimental data shows that parton shower approaches require intrinsic-\(k_{\textrm{T}}\) distributions dependent on the center-of-mass energy \(\sqrt{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_{\textrm{T}}\) distribution has a much milder center-of-mass energy \(\sqrt{s}\) dependence. We obtain more natural Gaussian width \(\sigma \), \( \sigma = q_s / \sqrt{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_{\textrm{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 Sect. 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. Section 3 is the central section of the paper, in which we perform fits to DY data and present results for the intrinsic-\(k_{\textrm{T}}\) TMD parameter. We give conclusions in Sect. 4.

2 PB TMDs and DY production

To study the different contributions to the low-\(p_{\textrm{T}}\) spectrum, at different \(m_{\text {DY}}\)  and different \(\sqrt{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] 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_{\textrm{T}}\), even if most of the small-\(k_{\textrm{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_{\textrm{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.

2.1 TMD distributions from the PB method

The PB evolution equations for TMD parton distributions \( \mathcal{A}_a ( x, \textbf{k}, \mu ^2) \) of flavor a are given by [9]

$$\begin{aligned}{} & {} { \mathcal{A}}_a(x,\textbf{k}, \mu ^2)= \Delta _a ( \mu ^2 ) \ { \mathcal{A}}_a(x,\textbf{k},\mu ^2_0) \nonumber \\{} & {} \qquad + \sum _b \int {{d^2 \textbf{q}^{\prime } } \over {\pi \textbf{q}^{\prime 2} } } \ { {\Delta _a ( \mu ^2 )} \over {\Delta _a ( \textbf{q}^{\prime 2} ) } } \ \Theta (\mu ^2-\textbf{q}^{\prime 2}) \ \Theta (\textbf{q}^{\prime 2} - \mu ^2_0) \nonumber \\{} & {} \qquad \times \int _x^{z_M} {{dz}\over z} \; P_{ab}^{(R)} (\alpha _s,z) \;{ \mathcal{A}}_b\left( {x \over z}, \textbf{k}+(1-z) \textbf{q}^\prime , \textbf{q}^{\prime 2}\right) ,\nonumber \\ \end{aligned}$$
(1)

where \(\textbf{k}\) and \(\textbf{q}\) are 2-dimensional momentum vectors, \(z_M\) is the soft resolution scale [10], z is the longitudinal momentum transferred at the branching, \(P_{ab}^{(R)} (\alpha _s, z ) \) are the resolvable splitting functionsFootnote 1 (whose explicit expressions for all flavor channels are given in [9]), and \( \Delta _a\) are the Sudakov form factors

$$\begin{aligned}{} & {} \Delta _a ( z_M, \mu ^2, \mu ^2_0 ) \nonumber \\{} & {} \quad =\exp \left( - \sum _b \int ^{\mu ^2}_{\mu ^2_0} {{d \textbf{q}^{\prime 2} } \over \textbf{q}^{\prime 2} } \int _0^{z_M} dz \ z \ P_{ab}^{(R)}\left( \alpha _s, z \right) \right) . \end{aligned}$$
(2)

The branching evolution (1) fulfills soft-gluon angular ordering [27,28,29], with the branching variable \(\textbf{q}^{\prime 2}\) being related to the transverse momentum \(q_{\textrm{T}}\) of the parton emitted at the branching by

$$\begin{aligned} q_{\textrm{T}}= (1-z) \, | \textbf{q}^{\prime } |. \end{aligned}$$
(3)

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 - \varepsilon \) for \(\varepsilon \rightarrow 0\). In contrast, \(p_T\) ordering leads, for instance, to ambiguities in the definition of the TMD from the \(z \rightarrow 1\) region.

Analogously to the case of ordinary (collinear) parton distribution functions, the distribution \( { \mathcal{A}}_a(x,\textbf{k},\mu ^2_0) \) at the starting scale \(\mu _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 \( { \mathcal{A}}_a(x,\textbf{k},\mu ^2_0) \) in the form

$$\begin{aligned} \mathcal{A}_{0,a} (x, \textbf{k},\mu _0^2) = f_{0,a} (x,\mu _0^2) \cdot \exp \left( -| \textbf{k}|^2 / 2 \sigma ^2\right) / ( 2 \pi \sigma ^2),\quad \end{aligned}$$
(4)

with the width of the Gaussian distribution given by \( \sigma = q_s / \sqrt{2} \), independent of parton flavor and x, where \(q_s\) is the intrinsic-\(k_{\textrm{T}}\) parameter.

The scale at which the strong coupling \(\alpha _s\) is to be evaluated in Eqs. (1) and (2) is a function of the branching variables. Two scenarios are studied in Refs. [9, 20]:

$$\begin{aligned}{} & {} \textrm{i}): \;\;\; \alpha _s = \alpha _s (\textbf{q}^{\prime 2}) \nonumber \\{} & {} \textrm{ii}): \;\;\; \alpha _s = \alpha _s (\textbf{q}^{\prime 2} (1-z)^2) = \alpha _s (q_{\textrm{T}}^2) \end{aligned}$$
(5)

In scenario i), it is shown in [9] that Eq. (1), in the collinear case, i.e. once it is integrated over all transverse momenta, reproduces exactly the DGLAP evolution [30,31,32,33] of parton densities. In scenario ii), it is discussed in [34] how, upon integration over transverse momenta and suitable treatment of the resolution scale, Eq. (1) returns the CMW coherent branching evolution [29].

In Ref. [20], fits to precision DIS HERA measurements [35] based on Eqs. (1) and (4), combined with NLO DIS matrix elements, are performed for both scenarios i) and ii), using the fitting platform xFitter [36, 37]. It is found that fits to DIS measurements with good \(\chi ^2\) values can be achieved in either case. Correspondingly, PB-NLO-HERAI+II-2018 set 1 (abbreviated as PB-NLO-2018 Set1) (with the DGLAP-type \( \alpha _s (\textbf{q}^{\prime 2}) \)) and PB-NLO-HERAI+II-2018 set2 (abbreviated as PB-NLO-2018 Set2) (with the angular-ordered CMW-type \( \alpha _s (q_{\textrm{T}}^2)\)) are obtained, both having intrinsic-\(k_{\textrm{T}}\) parameter in Eq. (4) set to \(q_s = \) 0.5 GeV [20]. All PB TMD parton distributions (and many others) are accessible in TMDlib and via the graphical interface TMDplotter [38, 39].

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/\gamma \) 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 \( \alpha _s (q_{\textrm{T}}^2)\) 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_{\textrm{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 \(\alpha _s\) needs to be specified in the region of small transverse momenta \(q_{\textrm{T}}\lesssim q_0\), where \(q_0 \) is a semi-hard scale on the order of a GeV. As in [7, 20], we take

$$\begin{aligned} \alpha _s = \alpha _s(\max (q^2_{0},q_{\textrm{T}}^2)), \end{aligned}$$
(6)

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].Footnote 2 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 \(\alpha _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 \(\textbf{q}^{\prime 2}\), it is instructive to examine separately parton emissions with transverse momenta above the semi-hard scale \(q_0\), \(q_{\textrm{T}}> q_0\), and below \(q_0\), . Using the angular ordering relation (3), these emissions are mapped respectively on the regions

(7)

where \( z_{\mathrm{{dyn}}} = 1 - q_0 / | \textbf{q}^{\prime } |\) 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, \( \alpha _s (q_{\textrm{T}}^2)\); 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 \(\alpha _s \) behavior (6): it does not contain any additional nonperturbative functions and parameters, besides the scale \(q_0\).

Fig. 1
figure 1

Integrated gluon and down-quark distributions at \(\mu = 4\) GeV (left column) and \(\mu = 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 \rightarrow 1\) (regions \((a+b)\) in text). The blue curve corresponds to \(z_M=z_{\textrm{dyn}}\) with \(q_0=1\) GeV (region (a) only). The ratio plots show the ratios to the one for \(z_M \rightarrow 1\)

Full size image

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 \rightarrow 1\) (default - regions \((a+b)\) - red curve [20]) and \(z_M=z_{\textrm{dyn}}= 1 - q_0/\textbf{q}^{\prime }\) (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 \rightarrow 1 \) are significantly different from those applying \(z_M=z_{\textrm{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 intrinsic-\(k_{\textrm{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 \rightarrow 1\), i.e. regions \((a+b)\), with (red curve) and without intrinsic-\(k_{\textrm{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_{\textrm{dyn}}= 1 - q_0/\mu ^\prime \), without intrinsic-\(k_{\textrm{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_{\textrm{T}}\) distribution (blue and magenta curves).

Fig. 2
figure 2

Transverse momentum distributions of down and charm quarks at \(\mu = 4\) GeV (left column) and \(\mu = 100\) GeV (right column) obtained from the PB approach based on PB-NLO-2018 Set2. Two distributions do not include intrinsic-\(k_{\textrm{T}}\): the blue curve corresponds to \(z_M \rightarrow 1\) (regions a+b in text) and the magenta curve to \(z_M=z_{\textrm{dyn}}= 1 - q_0/q\) (region a only). The red curve is the published one PB-NLO-2018 Set2 [20] and including intrinsic-\(k_{\textrm{T}}\) and \(z_M \rightarrow 1\)

Full size image

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_{\textrm{T}}> q_0\) the effect becomes smaller with increasing \(k_{\textrm{T}}\). Applying such a scale, \(z_M=z_{\textrm{dyn}}= 1 - q_0/\mu ^\prime \), removes emissions with \(q_{\textrm{T}}< q_0\) (there are still low-\(k_{\textrm{T}}\) contributions, which come from adding vectorially all intermediate emissions). However, very soft emissions are automatically included with \(z_M \rightarrow 1\).

As shown in Fig. 2, the effect of the intrinsic-\(k_{\textrm{T}}\) distribution is much reduced at large scales, but the contribution of the region \(z_{\textrm{dyn}}< z < 1\) stays important for small \(k_{\textrm{T}}\).

It is interesting to observe that the charm density shows essentially no effect of an intrinsic-\(k_{\textrm{T}}\) distribution: this is because charm is generated dynamically from gluons only, and there is no intrinsic charm density.

2.2 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_{\text {DY}}\) , and at low \(\sqrt{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_{\textrm{T}}(\ell \ell )\).

Fig. 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_{\textrm{T}}\) at \(\mu =100\) GeV and \(x=0.01\). In the upper row are shown distributions when no intrinsic-\(k_{\textrm{T}}\) distribution is included (\(q_s=0.00001\) GeV), and the lower row shows the default distributions with \(q_s=0.5\) GeV

Full size image

The difference between PB-NLO-2018 Set1 and Set2, which comes from the choice of renormalization scale (argument in \(\alpha _\textrm{s}\)), is seen essentially in the low-\(k_{\textrm{T}}\) region, where the nonperturbative Sudakov form factor (region (b)), with the integral \(z_M \rightarrow 1\), plays an important role. In Fig. 3 (upper row) the distributions for up and charm quarks are shown when no intrinsic-\(k_{\textrm{T}}\) distribution is included, the lower row shows distributions including the default intrinsic Gauss \(k_{\textrm{T}}\) distributions of widths \(q_s = 0.5\) GeV. It is very interesting to observe that 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 the fit of PB-NLO-2018 Set2.

Fig. 4
figure 4

TMD parton density distributions for down quarks of PB-NLO-2018 Set2 with (red curve) and without (blue curve) intrinsic-\(k_{\textrm{T}}\) distribution as a function of \(k_{\textrm{T}}\) at different scales \(\mu \) 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

Full size image

In Fig. 4 the transverse momentum distribution for down quarks, with and without an intrinsic-\(k_{\textrm{T}}\) distribution, is shown at different scales \(\mu \). While at low scales \(\mu \sim 50 \) GeV a significant effect of the intrinsic-\(k_{\textrm{T}}\) distribution is observed for very small \(k_{\textrm{T}}\), at large scales \(\mu \sim 350\) GeV this effect is much reduced. This scale dependence will result in a much smaller sensitivity to the intrinsic-\(k_{\textrm{T}}\) distribution at high \(m_{\text {DY}}\) .

Fig. 5
figure 5

Left: The mass distribution of DY lepton pairs at 13 TeV [51] compared to predictions of MCatNLO+CAS3 with PB-NLO-2018 Set1 (red curve), PB-NLO-2018 Set2 (blue curve) and without QED corrections (green curve). Right: The spectrum of photons transverse momentum in \(Z\rightarrow \mu ^+ \mu ^- \gamma \) at 7 TeV [55] compared to MCatNLO+CAS3 PB-NLO-2018 Set2 including QED radiation for a transverse momentum of the DY pair \(p_{\textrm{T}}(\ell \ell )< 10\) GeV. The bands show the scale uncertainty

Full size image

2.3 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 \(\alpha _\textrm{s}(m_{Z})=0.118\) at NLO.

The factorization scale \(\mu \), used in the calculation of the hard process is set to \(\mu = \frac{1}{2} \sum _i \sqrt{m^2_i +p^2 _{t,i}}\), 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 \(\mu \) in the hard process is set to \(\mu = m_{\text {DY}}\,\), in the case of a real emission it is set to \(\mu = \frac{1}{2} \sum _i \sqrt{m^2_i +p^2 _{t,i}}\). The generated transverse momentum is limited by the matching scale \(\mu _{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. Figure 5 (left) shows the DY mass distribution as measured by CMS [51] at 13 TeV together with predictions of MCatNLO+CAS3.Footnote 3 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 ”noQED”). A rather good description of the DY mass spectrum over a large range on \(m_{\text {DY}}\)  is obtained both with PB-NLO-2018 Set1 and Set2. Only at \(m_{\text {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_{\text {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_{\textrm{T}}\) region.

Fig. 6
figure 6

The \(p_{\textrm{T}}(\ell \ell )\) dependent DY cross section for different \(m_{\text {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

Full size image

3 The transverse momentum spectrum of DY lepton pairs

The transverse momentum spectrum of DY pairs at \(\sqrt{s}=13\) TeV has been measured for a wide \(m_{\text {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_{\textrm{T}}(\ell \ell )\), while PB-NLO-2018 Set2 describes the measurements rather well, without any further adjustment of parameters,Footnote 4 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_{\text {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_{\textrm{T}}(\ell \ell )\). All calculations predict too low a cross section at large transverse momentum 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 multi-jet merging technique gives an excellent description even for largest \(p_{\textrm{T}}(\ell \ell )\). For all further distributions, we restrict the investigations to \(p_{\textrm{T}}(\ell \ell )\) below the peak region (i.e. \(p_{\textrm{T}}(\ell \ell )\) \(\lesssim 8\) GeV).

3.1 Influence of the intrinsic-\(k_{\textrm{T}}\) distribution on DY transverse momentum distributions

Given the rather successful description of the DY \(p_{\textrm{T}}(\ell \ell )\)-spectrum with MCatNLO+CAS3 using PB-NLO-2018 Set 2 in the low \(p_{\textrm{T}}(\ell \ell )\)-region, we investigate below the importance of the intrinsic-\(k_{\textrm{T}}\) distribution. In PB-NLO-2018 the intrinsic-\(k_{\textrm{T}}\) distribution is parameterized as a Gauss distribution with zero mean and a width \(\sigma ^2 =q_s^2/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_{\textrm{T}}\) distribution, we show in Fig. 7 the MCatNLO+CAS3 predictions for the low \(p_{\textrm{T}}(\ell \ell )\)-spectrum of DY production at different DY masses \(m_{\text {DY}}\)  for different intrinsic-\(k_{\textrm{T}}\) distribution (with different \(q_s\) parameter values) compared to the CMS measurement  [17]. We observe that sensitivity to intrinsic-\(k_{\textrm{T}}\) is more pronounced at small \(p_{\textrm{T}}(\ell \ell )\) values. This sensitivity decreases with increasing mass, as expected from Fig. 4

Fig. 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_{\textrm{T}}(\ell \ell )\) for different \(m_{\text {DY}}\)  regions. Only the lowest \(p_{\textrm{T}}(\ell \ell )\) 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_{\textrm{T}}(\ell \ell )\) values show the maximal values included in the fit described in Sect. 3.2

Full size image

In the following we describe a determination of the Gaussian width \(q_s\) for different DY masses, \(m_{\text {DY}}\) , at different \(\sqrt{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_{\text {DY}}\) -bin \(q_s\) in steps of 0.1 to 0.3 GeV in the range \(q_s = 0.1, \dots , 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]).

3.2 Fit of the Gauss width \(q_s\) in pp at \(\sqrt{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_{\textrm{T}}\) parameter \(q_s\), since it covers a wide \(m_{\text {DY}}\) -range with high precision and that a detailed uncertainty breakdown, discussed in Sect. 3.2.1, is provided. The measurements of \(Z\)-production obtained from LHCb [58] are discussed in Sect. 3.2.2, while measurements at lower center-of-mass energies are shown in Sect. 3.3.

3.2.1 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_{\text {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.Footnote 5

In order to determine the intrinsic-\(k_{\textrm{T}}\) we vary the \(q_s\) parameter and calculate a \(\chi ^2\) to quantify the model agreement with the measurement. We evaluate the following expression,Footnote 6

$$\begin{aligned} \chi ^2 = \sum _{i,k} (m_i-\mu _i) C_{ik}^{-1}(m_k-\mu _k), \end{aligned}$$
(8)

with \(m_i\) being the measurement and \(\mu _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_{ik}^\text {meas.}\), and the statistical and scale uncertainties in the prediction,

$$\begin{aligned} C_{ik} = C_{ik}^\text {meas.} + C_{ik}^\text {model-stat.} + C_{ik}^\text {scale}. \end{aligned}$$
(9)

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,

$$\begin{aligned} C_{ik}^\text {model-stat.} = \sigma ^2_{i,\text {stat.}} \; \delta _{ik}, \end{aligned}$$
(10)

where \(\sigma _{i,\text {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_{\text {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,

$$\begin{aligned} C_{ik}^\text {scale} = \sigma _{i,\text {scale}} \; \sigma _{k,\text {scale}}, \end{aligned}$$
(11)

where \(\sigma _{i,\text {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_{\textrm{T}}(\ell \ell )<8\,\text {GeV}\,\). We reduce this range further in the first two regions to \(p_{\textrm{T}}(\ell \ell )<6\,\text {GeV}\,\) for \(50<m_{\text {DY}}\,<76\,\text {GeV}\,\) and \(p_{\textrm{T}}(\ell \ell )<7\,\text {GeV}\,\) for \(76<m_{\text {DY}}\,<106\,\text {GeV}\,\) to stay in the region of sensitivity and not be biased by missing higher orders in the predictions affecting high \(p_{\textrm{T}}(\ell \ell )\) shape. The obtained \(\chi ^2/\text {n.d.f}\) (reduced \(\chi ^2)\)) values are shown in Fig. 8 as a function of \(q_s\).

Fig. 8
figure 8

The reduced \(\chi ^2/\text {n.d.f}\) distribution as a function of \(q_s\) for different \(m_{\text {DY}}\) regions obtained from a comparison of the MCatNLO+CAS3 prediction with the measurement by CMS  [17]. The points represent the obtained \(\chi ^2\) values. The lines represent the curves used for the uncertainty estimate (see text), which is materialized by the shaded areas

Full size image
Table 1 Results of the fit on individual \(m_{\text {DY}}\) intervals for the CMS measurement  [17]. The “data” uncertainty is the one estimated using \(\min (\chi ^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 \(\text {n.d.f}\)
Full size table

Within each region, we consider the value of \(q_s\) for which the smallest \(\chi ^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 \(\chi ^2(q_s) < \min (\chi ^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_{\text {DY}}\)  interval are compatible with each other. The most precise determination is obtained from the \(Z\) peak region, \(76<m_{\text {DY}}\,<106\,\text {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_{\textrm{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 \(\chi ^2\) including the considered bins in all mass ranges. For this, we construct a new covariance matrix \(C_{ik}^\text {comb.}\) as a sum over the 650 uncertainty sources included in the detailed breakdown. We consider that each systematic uncertainty is fully correlated between \(m_{\text {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_{\text {DY}}\) region, and therefore we construct a block-diagonal matrix from the covariance matrices in each \(m_{\text {DY}}\) region. The statistical uncertainty in the prediction is diagonal. We consider that the uncertainties in the QCD scales are not correlated between \(m_{\text {DY}}\)  regions and use a block-diagonal matrix.

The \(\chi ^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,

$$\begin{aligned} q_s = 1.04\pm 0.03(\text {data})\pm 0.05(\text {scan})\pm 0.05(\text {binning})\,\text {GeV}\,. \end{aligned}$$

This value and its uncertainty are shown as a black line and shaded area on Fig. 9 for comparison with the individual \(m_{\text {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 \(\chi ^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].

Fig. 9
figure 9

The values of \(q_s\) in each \(m_{\text {DY}}\) -bin as obtained from Ref.  [17]. Indicated is also the combined fit value of \(q_s\)

Full size image
Fig. 10
figure 10

The reduced \(\chi ^2/\text {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 \(\chi ^2+1\). The best fit value is \(q_s = 0.74 \pm 0.15\) GeV . The value of \(q_s = 1.04 \pm 0.08\) GeV  as obtained from the measurements in Ref.  [17] is indicated by a black vertical line

Full size image
Fig. 11
figure 11

The reduced \(\chi ^2/\text {n.d.f.}\) distribution as a function of \(q_s\) obtained from a comparison of the MCatNLO+CAS3 PB-NLO-2018 Set2 prediction with the measurements at lower center-of-mass energies. The colored shaded band shows the \(\chi ^2\) variation of one unit for each data set. The value of \(q_s =1.04 \pm 0.08\) GeV  is shown as the grey band. Left-top: ATLAS measurement in 2 mass bins that we analysed at \(\sqrt{s}= 8\) TeV  (\(\text {n.d.f.} = 4\) for each mass bin) [61] and CMS in pPb at \(\sqrt{s}=~8.1\) TeV  [67]. Right-top: Tevatron measurements - D0 at \(\sqrt{s}= 1.8\) TeV  (\(\text {n.d.f.} = 4\)) [62], CDF at \(\sqrt{s}= 1.8\) TeV  (\(\text {n.d.f.} = 5\)) [63] and \(\sqrt{s}= 1.96\) TeV  (\(\text {n.d.f.} = 6\)) [64] Bottom: Measurements at lower energies - PHENIX at \(\sqrt{s}= 200\) GeV  (\(\text {n.d.f.} = 12\)) [65] and E605 at \(\sqrt{s}= 38.8\) GeV  (\(\text {n.d.f.} = 11\)) [66]

Full size image

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.

3.2.2 Z production at high rapidities at 13 TeV 

The LHCb collaboration [58] has measured \(Z\)-production at \(\sqrt{s} = 13\) TeV in the forward region, covering a rapdity range of \(2< |y| < 4.5 \).

The \(\chi ^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 \pm 0.15 \) GeV , where the uncertainty comes from a variation of \(\chi ^2\) by one unit and from the step size of the \(q_s\) scan.

3.3 The Gauss widths \(q_s\) from lower center of mass energies

The ATLAS collaboration has measured the production of DY from collisions at \(\sqrt{s} = ~8\) TeV  in several DY mass bins, out of which only the two with \(44< m_{\text {DY}}\,< 66\) GeV  and \(66< m_{\text {DY}}\,< 116\) GeV  are relevant for \(p_{\textrm{T}}(\ell \ell )< 10 \) GeV  [61]. In Fig. 11 we show the \(\chi ^2/\text {n.d.f}\) as a function of \(q_s\) obtained from these two mass bins (\(\text {n.d.f}=8\)).

The Tevatron experiments D0 [62] and CDF have measured transverse momenta of DY lepton pairs created in \(p{\bar{p}}\) collisions at lower center-of-mass energies (1.8 TeV [63] and 1.96 TeV [64]). The PHENIX collaboration measured DY production at \(\sqrt{s}=200\) GeV [65], and E605 [66] at \(\sqrt{s}=38.8\) GeV . The Drell–Yan differential cross section in \(p_{\textrm{T}}(\ell \ell )\) has also been measured in pPb data at \(\sqrt{s} = ~8.1\) TeV  by CMS  [67]. Figure 11 shows the impact that the \(q_s\) choice has on \(\chi ^2/\text {n.d.f}\) for these different measurements.

3.4 Consistency between determinations of intrinsic \(k_{\textrm{T}}\) width

A global fit of \(q_s\) is obtained by calculating \(\chi ^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_{\textrm{T}}\) distribution at lower collision energies has been analyzed using the entire range of \(p_{\textrm{T}}(\ell \ell )\), while at higher center-of-mass energies we investigate up to the peak region in the transverse momentum distribution.

Table 2 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\)
Full size table
Fig. 12
figure 12

The reduced \(\chi ^2/\text {n.d.f.}\) distribution \((\text {n.d.f.} = 81\)) as a function of \(q_s\) obtained from a comparison of the MCatNLO+CAS3 PB-NLO-2018 Set2 prediction with the measurement of Refs.  [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

Full size image

The \(\chi ^2/\text {n.d.f}\) distribution as a function of \(q_s\), for all the data together, is shown in Fig. 12. The \(\chi ^2\) distribution exhibits a minimum at around \(q_s = 1.0 \) GeV, which is consistent with the value obtained as described above.

Figure 13 displays the value of \(q_s\) as a function of \(m_{\text {DY}}\)  and \(\sqrt{s}\) obtained from the different 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 \(\chi ^2\) variation of one unit and step size of the \(q_s\) scan only. The value of \(q_s=1.04 \pm 0.08\) GeV   as derived from the measurements in Ref. [17], is compatible for all ranges of \(m_{\text {DY}}\) , and also holds true for various values of \(\sqrt{s}\). The obtained value is also found to be compatible for pPb data.

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 \pm 0.08\) GeV . This value, in contrast to standard Monte Carlo event generators, has no strong dependence on the center-of-mass energy as well as on the mass of the produced Drell–Yan lepton pair \(m_{\text {DY}}\) . The results of this section indicate that the treatment of soft emissions in the region with the strong coupling of Eq. (6) applied in PB-NLO-2018 Set2 leads to intrinsic-\(k_{\textrm{T}}\) distributions with width parameter \(q_s\) consistent with Fermi motion kinematics, and mildly varying with energy.

4 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_{\textrm{T}}\) differential cross sections performed recently at the LHC at \(\sqrt{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.

Fig. 13
figure 13

Left: the value of \(q_s\) as a function of the DY-mass as obtained from the measurements in Refs.  [17, 58, 61,62,63,64,65,66,67]. Right: same as a function of \(\sqrt{s}\). The value of \(q_s = 1.04 \pm 0.08\) GeV  as obtained from the measurements in Ref.  [17] is indicated

Full size image

As shown previously, the measured DY cross section at low \(p_{\textrm{T}}\) favours a choice of the strong coupling \(\alpha _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 emissions with \(q_T > q_0\), populating the phase space region \( z < z_{\mathrm{{dyn}}}\) (where \(z_{\mathrm{{dyn}}} = 1 - q_0 / | \textbf{q}^{\prime } |\), with \(| \textbf{q}^{\prime } |\) being the scale of the branching), while it is evaluated at the semi-hard scale \(q_0\) for emissions with . The contribution to the Sudakov evolution from the parton branching in the phase space region \( z < z_{\mathrm{{dyn}}}\) gives the perturbative resummation of Sudakov logarithms, while the contribution from the parton branching in the phase space region gives the nonperturbative Sudakov form factor. Therefore the PB-NLO-2018 Set2 contains two sources of nonperturbative effects: i) the nonperturbative TMD distribution at low evolution scale \(\mu _0\), and ii) the nonperturbative Sudakov form factor, specified by the “pre-confinement” scale prescription to continue the branching evolution to the infrared region . The former includes the intrinsic-\(k_T\) width parameter \(q_s\), corresponding to Fermi motion in the hadron beam, while the latter is characterized by the semi-hard scale parameter \(q_0\). At low \(k_{\textrm{T}}\), the contribution of nonperturbative Sudakov form factor interplays with the contribution of the intrinsic transverse momentum.

The main result of the present work is the extraction, within the PB-NLO-2018 Set2 framework, of the intrinsic-\(k_T\) Gauss distribution with zero mean and width parameter \(q_s= \sqrt{2} \sigma \) from the measured \(p_{\textrm{T}}\) dependence of the DY cross sections obtained recently at the LHC at \(\sqrt{s}= 13\) TeV [17], for different DY masses \(m_{\text {DY}}\) , between 50 GeV and 1 TeV. These measurements provide a complete decomposition of the different systematic uncertainties and their covariance matrices. To compare to the data, we have used DY production at NLO obtained with the MadGraph5_aMC@NLO event generator matched with the PB TMD distributions PB-NLO-2018 Set2, with a given parameter \(q_s\) value. We performed a scan over a large range of values \(q_s\) on the transverse momentum spectrum below the peak, i.e. the sensitive part to intrinsic-\(k_{\textrm{T}}\), and considering separately each experimental source of uncertainty and their correlations. The theory scale uncertainties have been considered to be fully correlated inside each \(m_{\text {DY}}\) bin and uncorrelated between \(m_{\text {DY}}\) bins. We found the value \(q_s = 1.04 \pm 0.08\) GeV , consistent for the different \(m_{\text {DY}}\) . The obtained value is in agreement to the expected value from Fermi-motion in protons. It has been cross checked that this value is compatible with \(q_s\) values obtained from other DY measurements at different center-of-mass energies \(\sqrt{s}\) and for a variety of DY masses. The global picture shows no strong dependence of the intrinsic-\(k_{\textrm{T}}\) on the center-of-mass energy or on the DY mass, which contrasts with tuned standard Monte Carlo event generators that need a strongly increasing intrinsic Gauss width with \(\sqrt{s}\) and with \(m_{\text {DY}}\) .

We suggest that the remarkably stable value of \(q_s\) that we obtain in our study can be attributed to the contribution of the nonperturbative Sudakov form factor and the treatment of the region near the soft-gluon resolution boundary.