The transverse momentum spectrum of low mass Drell-Yan production at next-to-leading order in the parton branching method

The transverse momentum spectrum of low mass Drell-Yan (DY) production at low center-of-mass energies $\sqrt{s}$ is calculated by applying transverse momentum dependent (TMD) parton distributions obtained from the Parton Branching (PB) evolution method, combined with the next-to-leading-order (NLO) calculation of the hard process in the MCatNLO method. We compare our predictions with experimental measurements at low DY mass, and find very good agreement. In addition we use the low mass DY measurements at low $\sqrt{s}$ to determine the width $q_s$ of the intrinsic Gauss distribution of the PB-TMDs at low evolution scales. We find values close to what has earlier been used in applications of PB-TMDs to high-energy processes at the Large Hadron Collider (LHC) and HERA. We find that at low DY mass and low $\sqrt{s}$ even in the region of $p_t/m_{DY} \sim 1$ the contribution of multiple soft gluon emissions (included in the PB-TMDs) is essential to describe the measurements, while at larger masses ($m_{DY} \sim m_{Z}$) and LHC energies the contribution from soft gluons in the region of $p_t/m_{DY}\sim 1$ is small.

The transverse momentum spectrum of DY production at lower mass m DY allows one to study in more detail the non-perturbative contribution, as the phase space for perturbative evolution is reduced. However, the measurement of the transverse momentum at low mass of the DY pair is experimentally very challenging, since one has to measure down to low transverse momenta of the decay leptons, where experimental background and misidentification of the DY lepton pairs can be significant. At the LHC the lowest DY mass used for the low transverse momentum spectra (p T > ∼ 1 GeV) is ∼ 46 GeV [3], while at lower centerof-mass energies DY measurements covering the low p T region for lower masses exist from PHENIX [40] at √ s = 200 GeV, from R209 [41] at √ s = 62 GeV, and from NuSea [42,43] and E605 [44] at √ s = 38.8 GeV. The description of these measurements is discussed in terms of TMDs in Refs. [17,45,46]. In Ref. [17] these measurements are compared with collinear NLO predictions and significant differences were observed. In this paper we apply the TMD parton densities obtained using the PB method (fitted [47] to inclusive deep-inelastic scattering (DIS) precision data from HERA) together with an NLO calculation of DY production [34] precisely in the same manner as in Ref. [36], but now to treat low-mass DY production. We first briefly review the main elements of the PB approach (Sec. 2) and the matching of the PB-TMDs with the NLO calculation (Sec. 3). Then we show that these low energy measurements are very well described with the PB-MCatNLO approach in the whole region of p T /m DY (in contrast to the observation in Ref. [17]) and discuss the role of the non-perturbative (intrinsic-k T ) distribution (Sec. 4). We finally give conclusions (Sec. 5).

Collinear and TMD densities from the PB method
The PB method [37,38] allows evolution equations for collinear and TMD parton distributions to be solved numerically in an iterative procedure, by making use of the concept of resolvable and non-resolvable branchings and by applying Sudakov form factors to describe the evolution from one scale to another without resolvable branching. The method is based on introducing the soft-gluon resolution scale z M into the QCD evolution equations to separate resolvable and non-resolvable emissions, and treating them via, respectively, the resolvable splitting probabilities P (R) ba (α s , z) and the Sudakov form factors Here a, b are flavor indices, α s is the strong coupling, z is the longitudinal momentum splitting variable, and z M < 1 is the soft-gluon resolution parameter. The PB method is described in detail in Refs. [38,47].            The TMD parton density distributions are obtained from convolving the evolution kernels K ba , given in terms of the splitting functions P (R) ba , Sudakov form factors ∆ a and phase space constraints, with the non-perturbative starting distributions A 0,b (x , k 2 T,0 , µ 2 0 ). As described in [47], we have In general, the starting distribution A 0 can have flavor-dependent and x-dependent k T,0 distributions. However, for maximal simplicity we use here a factorized form, in which the intrinsic k T,0 distribution is given by a Gauss distribution with σ 2 = q 2 s /2 for all parton flavors and all x, with a constant value q s = 0.5 GeV. Also, the evolution kernels K ba in Eq. (1) do not include any non-perturbative component. In principle, non-perturbative contributions to Sudakov form factors could be introduced in the K ba kernels of the PB method, and parameterized in terms of non-perturbative functions to be fitted to experimental data (similarly to what is done in other approaches, e.g. [19][20][21][22]). For simplicity, however, at present we take the kernels K ba to be purely perturbative. The PB method enables the explicit calculation of the kinematics at every branching vertex, once the evolution scale is specified in terms of kinematic variables. In Ref. [37] it was pointed out that angular ordering gives transverse momentum distributions which are stable with respect to variations of the resolution parameter z M . In angular ordering, the angles of the emitted partons increase from the hadron side towards the hard scattering [49,50]. The transverse momentum of the i's emitted parton q t i can be calculated in terms of the angle Θ i of the emitted parton with respect to the beam directions from q t,i = (1 − z i )E i sin Θ i . Associating the "angle" E i sin Θ i with µ i gives Collinear and TMD parton distributions were obtained in [47] from fits of the parameters of the starting distributions at scale µ 0 ∼ 1 GeV to the inclusive-DIS precision measurements from HERA, after PB evolution and convolution with the DIS hard-scattering coefficient functions at NLO. The fits were performed using the open-source fitting platform xFitter [51] and a new development described in Ref. [38,47] of the numerical techniques [52]. Collinear and TMD distributions were extracted including the determination of experimental and theoretical uncertainties. Two different sets of parton distributions were obtained: Set 1, which corresponds at collinear level to HERAPDF 2.0 NLO [53], and Set 2, which differs by the choice of the scale used in the running coupling α s , namely, it uses the transverse momentum (instead of the evolution scale), corresponding to the angularordering approach. For soft gluon resolution z M → 1 and strong coupling α s → α s (µ 2 ) it was verified [38,47] numerically, with a numerical accuracy of better than 1 % over a range of five orders of magnitude both in x and in µ, that DGLAP evolution equations [54][55][56][57] are recovered from PB evolution.
In Fig. 1 the Set 1 and Set 2 collinear densities are shown for up-quark and down-quark at evolution scales of µ = 10 and 100 GeV. The plots in Figs. 1,2 are made using the TMDplotter tool [58,59]. Collinear densities are available in a format compatible with the one employed in LHAPDF [60], and can be used in calculations of physical processes at NLO.
In Fig. 2 we show the TMD distributions for up-quarks at x = 0.01 and µ = 10 and 100 GeV.

PB-TMDs and DY production at NLO
In this section we employ the approach proposed in Ref. [36] to perform the matching of PB-TMDs with the NLO calculation of DY production. We briefly describe a few technical aspects of the computation and analyze numerically the contributions of PB-TMDs and NLO in the matching procedure.
In order to avoid double counting between the contribution of the real emission treated by the matrix element calculation and the contribution from the PB-TMD (which in this respect plays a role analogous to that of parton showers), a matching scale µ m needs to be defined. This scale is determined by the NLO calculation and is transferred to the user via the parameter SCALUP (included in the LHE file).
The PB-TMDs depend, as indicated in Eq. the decay products and the final jet. In order to allow the full phase space to be covered for the transverse momentum in the PB-TMD, the factorization scale µ is set to the invariant mass of the hard process µ = √ŝ for the underlying Born configuration. For the real emission configuration the scale is changed to µ = 1 2 i m 2 i + p 2 t,i , as in the MC@NLO calculation. Finally, the transverse momentum is constrained to be smaller than the matching scale µ m =SCALUP. The calculation are performed with the CASCADE3 package [63,64] (version 3.0.X), which allows us to read LHE files and to produce output files to be analyzed with Rivet [65].
In Fig. 3 we show results, at different center-of-mass energies √ s, for the DY leptonpair transverse momentum distribution, obtained from the MC@NLO calculation at a purely partonic level (LHE level) using HERWIG6 subtraction terms (red solid curves in the plots), and from the MC@NLO calculation after inclusion of PB-TMDs (blue solid curves). It is interesting to observe that the contribution coming from the real hard partonic emission is small at low center-of-mass energies and at low p T , but increases with increasing √ s, thus allowing one to study the contribution of multiple soft emissions in detail. In Fig. 4 the distribution in transverse momentum, with subtraction terms and after in-

Low mass DY production
We next apply the framework described in the previous section, based on the matching of PB-TMDs with NLO, to the evaluation of DY spectra at low DY masses.

Mass and transverse momentum spectra
We start with the DY mass spectrum at low masses and low √ s. In Fig. 5 we present theoretical predictions obtained from PB-TMDs and NLO matrix elements using MC@NLO matching, and compare them with experimental measurements for different center-of-mass energies from NuSea [42,43], R209 [41] and PHENIX [40]. We also show the theoretical uncertainties coming from the determination of the PB-TMDs as well as from the variation of the scale in the perturbative calculation. As already observed in Ref. [36] for the case of Zproduction at the LHC, the contribution to uncertainties from the parton density turns out to be small compared to the one from the scale uncertainty. Not included are the uncertainties coming from the variation of the intrinsic Gauss distribution (q s ), as this parameter was not constrained by the fits to HERA data [47]. This will be further discussed in Subsec  [42,43], R209 [41] and PHENIX [40] compared to predictions at NLO using PB-TMDs. For NuSea also the prediction using NNPDF3.0 [48] is shown.
deteriorates. This is because we enter the large-x region where the parton densities [47] used in the calculation, which are determined from fits to HERA data [53], are poorly constrained. The description in this region can be readily improved by using parton density sets from global fits. We show this in Fig. 5 by plotting the result from the set NNPDF3.0 [48], obtained from global fits that include NuSea data [42,43]. On the other hand, for the lowest mass region m DY < 6 GeV of NuSea the mass spectrum is well described. We use this region to investigate the transverse momentum spectrum. In Fig. 6 we present theoretical predictions from PB-TMDs and NLO matrix elements for transverse momentum spectra, and again we compare them with experimental measure-ments for different center-of-mass energies from NuSea [42,43], R209 [41] and PHENIX [40]. The PB-TMDs used in the calculation include an intrinsic (non-perturbative) transverse momentum spectrum parameterized as a Gauss distribution with width σ 2 = q 2 s /2 (see eq. (2)). The quality of the description of the measurements (including independent variations of R209 [41], PHENIX [40] compared to predictions at NLO using PB-TMDs. the factorization and renormalization scales by a factor of two up and down) is good with χ 2 /ndf = 1.08, 1.27, 1.04 for NuSea, R209 and PHENIX, respectively. In Fig. 7 we show the transverse momentum spectrum of Z-bosons at LHC energies of √ s = 13 TeV as measured by CMS [7] and compare it with predictions using the same method of the above low-energy predictions and of Ref. [36], with the PB-TMD Set 2. We observe a very good description of the measurement (with χ 2 /ndf = 0.8 for p T < 80 GeV).
The drop in the prediction at large transverse momenta comes from missing higher order contributions in the hard process calculation, as discussed in Ref. [36].

Determination of the non-perturbative (intrinsic) transverse momentum distribution
The low-mass DY measurements can be used to constrain the intrinsic transverse momentum distribution. In Fig. 8 we report the calculated χ 2 /ndf as a function of q s obtained from the transverse momentum distributions of NuSea [42,43], R209 [41], PHENIX [40] (as shown in Fig. 6). For the calculation of χ 2 /ndf we use the full experimental uncertainties (except an overall normalization uncertainty) and the central values for the theory predictions (without inclusion of pdf and scale uncertainties, leading to a larger χ 2 /ndf as the one reported in the previous subsection). A clear minimum is found for NuSea and R209 measurements, with values of q s ∼ 0.3 − 0.4 GeV. On the other hand, the PHENIX measurement shows little sensitivity to the choice of q s , which is understandable since only two values for p T < 1 GeV are measured, while the other experiments have a finer binning. It is interesting to note that the values of intrinsic transverse momentum determined from low-mass DY are rather close to the value of q s = 0.5 GeV that was assumed in PB-Set2 [47], determined from fits to inclusive DIS data from HERA which are not sensitive to intrinsic-k T . obtained from a comparison of the measurements (NuSea [42,43], R209 [41], PHENIX [40]) with a prediction at NLO using PB-TMDs. For the theory prediction only the central value is taken, but no uncertainty from scale variation is included.

Comments on the low-mass region
It has been observed in [17] that NLO calculations in collinear factorization are not able to describe the measurements of DY transverse momentum spectra at fixed-target experiments in the region p T /m DY ∼ 1. We remark that this is consistent with the observation which we have made in Fig. 3 that, in this kinematic region, the contribution from the real hard emission is small compared to the contribution from multiple parton radiation, embodied in the PB-TMD evolution. Indeed, Fig. 3 indicates that a purely collinear NLO calculation would not give a realistic description of the DY spectrum for p T /m DY ∼ 1 at low energies. On the other hand, Fig. 4 illustrates that the situation is very different at the LHC: in the region around the Z mass shown in Fig. 4, hard real emission dominates the transverse momentum spectrum for p T /m DY ∼ 1, so that a purely collinear NLO calculation gives a good approximation to the DY process for p T /m DY ∼ 1 at the LHC.
The comparison of theoretical predictions with transverse momentum measurements from NuSea [42,43] in the top left panel of Fig. 6 confirms that the inclusion of multiple parton emissions, taken into account by the PB-TMD evolution equation [38] (see also discussion in [66]), is essential to describe the region p T /m DY ∼ 1 at low energies. This physical picture is supported by the comparison of theoretical predictions with measurements at the increasingly high energies of R209 [41] and PHENIX [40] in the top right and bottom panels of Fig. 6. Going up to LHC energies in Fig. 7, we see that the PB-TMD + NLO calculation describes the spectrum well all the way up to transverse momenta p T ∼ m DY (while for even higher p T a deficit is observed due to the missing DY + jet NLO correction -see discussion in [36]).
To sum up, the DY transverse momentum in the low-mass region is sensitive to both finite-order QCD contributions and all-order QCD multi-parton radiation. Theoretical predictions depend on the matching procedure between these contributions. Once this is accomplished, low-mass DY measurements are well described and can provide a wealth of information on non-perturbative QCD dynamics. In this paper the matching is performed, in the spirit of [67], with PB-TMDs and MC@NLO (alternative methods of matching are e.g. those inspired by [12]).

Conclusion
We have investigated the transverse momentum spectra of DY lepton-pair production at small DY masses and low center-of-mass energies by matching PB-TMD distributions to NLO calculations via MC@NLO. We use the same PB-TMDs and MC@NLO calculations as we have used for Z-production at LHC energies in Ref. [36]. We observe a very good description of the measurements by the NuSea collaboration at √ s = 38 GeV, R209 at √ s = 62 GeV and PHENIX at √ s = 200 GeV, with values of χ 2 /ndf ∼ 1 for all measurements. We use the low-mass DY measurements to determine the best value for the width of the intrinsic Gauss distribution, and find a value of q s ∼ 0.3 − 0.4 GeV, slightly smaller than q s = 0.5 GeV used in the PB-TMD Set 2 distributions [47]. The very good description of low-mass DY measurements is achieved by a combination of a collinear NLO calculation (including the appropriate subtraction terms to avoid double counting) with the PB-TMDs. We find that, at low DY mass and low √ s, even in the region of p T /m DY ∼ 1 the contribution of QCD multi-parton radiation (included in the evolution of PB-TMDs in terms of Sudakov form factors, resolvable splitting functions and phase space constraints) is essential to describe the measurements, while at larger masses (m DY ∼ m Z ) and LHC energies this contribution is small in the region of p T /m DY ∼ 1.