The small $k_{\rm T}$ region in Drell-Yan production at next-to-leading order with the Parton Branching Method

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

where $\bf k$ and $\bf 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 functions^***Using transverse momentum dependent splitting functions as described in Ref.[26] would require using off-shell matrix elements and a completely new fit to inclusive structure functions. (whose explicit expressions for all flavor channels are given in [9]), and $\Delta_{a}$ are the Sudakov form factors

The branching evolution (1) fulfills soft-gluon angular ordering [27, 28, 29], with the branching variable ${\bf q}^{\prime 2}$ being related to the transverse momentum $q_{\rm 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-\varepsilon$ for $\varepsilon\to 0$. In contrast, $p_{T}$ ordering leads, for instance, to ambiguities in the definition of the TMD from the $z\to 1$ region.

Analogously to the case of ordinary (collinear) parton distribution functions, the distribution ${{\cal A}}_{a}(x,{\bf 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 ${{\cal A}}_{a}(x,{\bf k},\mu^{2}_{0})$ in the form

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_{\rm 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. [20, 9]:

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}({\bf 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_{\rm T}^{2})$) are obtained, both having intrinsic-$k_{\rm 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_{\rm 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_{\rm T}$ parameter $q_{s}$, and perform determinations of its value.

As discussed in  [20, 7], 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_{\rm T}\ \hbox{\raise 2.0pt\hbox{$<$} \kern-13.0pt\lower 3.0pt\hbox{$\sim$}}\ q_{0}$, where $q_{0}$ is a semi-hard scale on the order of a GeV. As in [20, 7], 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] ^†††Different forms of the extension to small $q_{\rm T}$ could be considered. However, this will entail new fits both to precision DIS data and to DY data.. 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 ${\bf q}^{\prime 2}$, it is instructive to examine separately parton emissions with transverse momenta above the semi-hard scale $q_{0}$, $q_{\rm T}>q_{0}$, and below $q_{0}$, $\Lambda_{\rm{QCD}}<q_{\rm T}\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\thinspace$\sim$}\hss}\raise 1.0pt\hbox{$<$}}q_{0}$. Using the angular ordering relation (3), these emissions are mapped respectively on the regions

where $z_{\rm{dyn}}=1-q_{0}/|{\bf 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_{\rm 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}$.

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}\to 1$ (default - regions $(a+b)$ - red curve [20]) and $z_{M}=z_{\rm dyn}=1-q_{0}/{\bf 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}\to 1$ are significantly different from those applying $z_{M}=z_{\rm 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_{\rm 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}\to 1$, i.e. regions $(a+b)$, with (red curve) and without intrinsic-$k_{\rm 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_{\rm dyn}=1-q_{0}/\mu^{\prime}$, without intrinsic-$k_{\rm 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_{\rm 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_{\rm T}>q_{0}$ the effect becomes smaller with increasing $k_{\rm T}$. Applying such a scale, $z_{M}=z_{\rm dyn}=1-q_{0}/\mu^{\prime}$, removes emissions with $q_{\rm T}<q_{0}$ (there are still low-$k_{\rm T}$ contributions, which come from adding vectorially all intermediate emissions). However, very soft emissions are automatically included with $z_{M}\to 1$.

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

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

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_{\small\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_{\rm T}(\ell\ell)$.

The difference between PB-NLO-2018 Set1 and Set2, which comes from the choice of renormalization scale (argument in $\alpha_{\mathrm{s}}$), is seen essentially in the low-$k_{\rm T}$ region, where the nonperturbative Sudakov form factor (region $(b)$), with the integral $z_{M}\to 1$, plays an important role. In Fig. 3 (upper row) the distributions for up and charm quarks are shown when no intrinsic-$k_{\rm T}$ distribution is included, the lower row shows distributions including the default intrinsic Gauss $k_{\rm 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.

In Fig. 4 the transverse momentum distribution for down quarks, with and without an intrinsic-$k_{\rm T}$ distribution, is shown at different scales $\mu$. While at low scales $\mu\sim 50$ GeV a significant effect of the intrinsic-$k_{\rm T}$ distribution is observed for very small $k_{\rm 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_{\rm T}$ distribution at high $m_{\small\text{DY}}$.

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_{\mathrm{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_{\small\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. Fig. 5 (left) shows the DY mass distribution as measured by CMS [51] at 13 TeV together with predictions of MCatNLO+CAS3 ^‡‡‡We use the Rivet package [52] for the calculation of the final distributions.. 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_{\small\text{DY}}$ is obtained both with PB-NLO-2018 Set1 and Set2. Only at $m_{\small\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_{\small\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_{\rm T}$ region.

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

In the following we describe a determination of the Gaussian width $q_{s}$ for different DY masses, $m_{\small\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_{\small\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. [56, 11, 57]).

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_{\rm T}$ parameter $q_{s}$, since it covers a wide $m_{\small\text{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.

The ATLAS collaboration has measured the production of DY from $\mathup{{{p}}}$ ​​$\mathup{{{p}}}$ collisions at $\sqrt{s}=~8$ TeV in several DY mass bins, out of which only the two with $44<m_{\small\text{DY}}<66$ GeV and $66<m_{\small\text{DY}}<116$ GeV are relevant for $p_{\rm 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_{\rm 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.

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_{\rm T}$ distribution at lower collision energies has been analyzed using the entire range of $p_{\rm T}(\ell\ell)$, while at higher center-of-mass energies we investigate up to the peak region in the transverse momentum distribution.

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_{\small\text{DY}}$ and $\sqrt{s}$ obtained from the different measurements in Refs. [17, 58, 67, 61, 63, 62, 64, 65, 66]. 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_{\small\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_{\small\text{DY}}$. The results of this section indicate that the treatment of soft emissions in the region $z_{\rm dyn}\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\thinspace$\sim$}\hss}\raise 1.0pt\hbox{$<$}}z<z_{M}$ with the strong coupling of Eq. (6) applied in PB-NLO-2018 Set2 leads to intrinsic-$k_{\rm T}$ distributions with width parameter $q_{s}$ consistent with Fermi motion kinematics, and mildly varying with energy.