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

It has been observed in the literature that measurements of low-mass Drell–Yan (DY) transverse momentum spectra at low center-of-mass energies s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s}$$\end{document} are not well described by perturbative QCD calculations in collinear factorization in the region where transverse momenta are comparable with the DY mass. We examine this issue from the standpoint of the Parton Branching (PB) method, combining next-to-leading-order (NLO) calculations of the hard process with the evolution of transverse momentum dependent (TMD) parton distributions. 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 s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s}$$\end{document} to determine the width qs\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_s$$\end{document} 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 s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s}$$\end{document} even in the region of pT/mDY∼1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_\mathrm{T}/m_\mathrm{DY}\sim 1$$\end{document} the contribution of multiple soft gluon emissions (included in the PB-TMDs) is essential to describe the measurements, while at larger masses (mDY∼mZ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_\mathrm{DY}\sim m_{{\mathrm{Z}}}$$\end{document}) and LHC energies the contribution from soft gluons in the region of pT/mDY∼1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_\mathrm{T}/m_\mathrm{DY}\sim 1$$\end{document} 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 1GeV) is ∼ 46 GeV [3], while at lower center-of-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. In a study [45] based on the Monte Carlo event generator Herwig, good agreement with measurements at low √ s was found after changing parameters for the parton shower and intrinsic transverse momentum. The description of these measurements is discussed in terms of TMDs in Refs. [17,[46][47][48]. In Ref. [17] these measurements were compared with collinear NLO predictions and significant discrepancies were observed.
In this paper we apply the TMD parton densities obtained using the PB method (fitted [49] 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 and the matching of the PB-TMDs with the NLO calculation (Sect. 2). 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 examine the role of both the perturbative evolution and the non-perturbative (intrinsic-k T ) distribution (Sect. 3). We provide a discussion to put these results in a broader context (Sect. 4), and finally give conclusions (Sect. 5).

PB-TMDs and DY production at NLO
In this section we recall the basic elements of the PB approach, and illustrate the main features of applying it to DY production at different center-of-mass energies, from fixedtarget experiments to the LHC.

Collinear and TMD densities from the PB method
The approach proposed in [37] allows evolution equations for both collinear and TMD parton distributions to be solved numerically with the PB method. In this approach, the concept of resolvable and non-resolvable branchings is applied by using Sudakov form factors. A soft-gluon resolution scale z M is introduced to separate resolvable and non-resolvable branchings. The Sudakov form factors, which describe the evolution without resolvable branching from one scale μ 0 to another scale μ, are given in terms of the resolvable splitting probabilities P (R) ba (α s , z) as follows, where 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. A detailed description of the PB method is given in Refs. [38,49]. The TMD parton density distributions are obtained from the non-perturbative starting distributions A 0,b (x , k 2 T,0 , μ 2 0 ) after convoluting with a perturbative evolution kernel K ba . As described in [49], we have We use a factorized form for the starting distribution A 0 , for simplicity, (in general, the k T,0 distribution can be also flavor and x-dependent), with σ 2 = q 2 s /2 independent of the parton flavor and x, with a constant value q s = 0.5 GeV. Also, the evolution kernels K ba in Eq. (2) 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 [51,52]. 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 Deep-inelastic scattering measurements from HERA are used in Ref. [49] to determine the free parameters of the               GeV and x = 0.01. In the lower panels the full uncertainty of the TMDs is shown, as obtained from the fits [49] starting distributions at scale μ 0 ∼ 1 GeV. The fits were performed using the open-source fitting platform xFitter [53] and a new development described in Ref. [38,49] of the numerical techniques [54]. Collinear and TMD distributions were extracted including the determination of experimental and theoretical uncertainties. In Ref. [49] two sets of parton distributions are described: Set 1, which uses the evolution scale as argument in the running coupling α s , similar to what is used in HERAPDF 2.0 NLO [55], and Set 2, which uses the transverse momentum in the evolution of α s .
For soft gluon resolution z M → 1 and strong coupling α s → α s (μ 2 ) it was verified [38,49] 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 [56][57][58][59] 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. In Fig. 2 we show the TMD distributions for up-quarks at x = 0.01 and μ = 10 and 100 GeV. The plots in Figs. 1 and 2 are made using the TMDplotter tool [60,61]. Collinear densities are given in a format compatible with LHAPDF [62].

Matching PB-TMDs with NLO calculations
We employ the approach proposed in Ref. [36] to perform the matching of PB-TMDs with the NLO calculation of DY production. In this subsection we briefly describe a few technical aspects of the computation and analyze numerically the contributions of PB-TMDs and NLO in the matching procedure.
Following [36], MadGraph5_aMC@NLO (version 2.6.4, hereafter labelled MC@NLO) [34] is used to calculate the Drell-Yan process at NLO, i.e., including O(α s ) corrections to the hard-scattering matrix element, together with the NLO PB parton distributions (Set 2) of Ref. [49]. As in [36], motivated by the angular ordering in PB evolution, we use Herwig6 [63,64] subtraction terms in MC@NLO. A similar method to describe DY production at leading order is proposed in Ref. [65].
A matching scale μ m (parameter SCALUP) separates the contribution of the real emission treated by the matrix element calculation and the contribution from the PB-TMD.
The hard process is calculated at a scale μ and the longitu- , with the sum running over the decay products and the final jet. The same scales are used in the PB-TMD. The scale μ m =SCALUP is also used as an upper limit for the transverse momentum (the calculation are performed with the Cascade3 package [66,67] (version 3.0.X)). We employ Rivet [68] to analyze output files.
In Fig. 3 we show results, at different center-of-mass energies √ s, for the DY lepton-pair 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 inclusion of PB-TMDs, is shown for LHC energies of √ s = 13 TeV and for DY masses around the Z-mass. At high √ s = 13 TeV and at sufficiently large DY mass, the predictions with and without PB-TMDs become similar at large transverse momenta, supporting the simple expectation that for p T /m DY 1 the trans-

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 Z-production 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 [49]. This will be further discussed in Sect. 3.2.
The mass spectra in Fig. 5 are generally well described by the PB-TMD + NLO calculation. For the region of highest masses at lowest √ s (NuSea experiment), we see in the top panel of Fig. 5 that the description of experimental data by the PB-TMD + NLO calculation deteriorates. This is because we enter the large-x region where the parton densities [49] used in the calculation, which are determined from fits to HERA data [55], 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 [50], 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 experimen- 5 Drell-Yan mass distribution production measured by NuSea [42,43], R209 [41] and PHENIX [40] compared to predictions at NLO using PB-TMDs. For NuSea also the prediction using NNPDF3.0 [50] is shown Fig. 6 Transverse momentum spectrum of Drell-Yan production measured by NuSea [42,43], R209 [41], PHENIX [40] compared to predictions at NLO using PB-TMDs 123 Fig. 7 Transverse momentum spectrum of Z production measured by CMS [7] compared to predictions at NLO using PB-TMDs tal measurements 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. (3)). The quality of the description of the measurements (including independent variations of the factorization and renormalization scales by a factor of two up and down) is good with χ 2 /nd f = 1.08, 1.27, 1.04 for NuSea, R209 and PHENIX, respectively. The χ 2 values are calculated using the full p T range. In the above discussion we have shown results for NuSea, R209 and PHENIX as representative of a broad range of different center-of-mass energies. We have obtained similar results for other data sets in this energy range such as E605 [44].
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 /nd f = 0.8 for p T < 80 GeV). As discussed in Ref. [36], the drop in the prediction at large transverse momenta comes from missing NLO contributions to Z+ jet production, i.e., O(α 2 s ) terms in the hard process calculation.

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 /nd f as a function of q s obtained Fig. 8 The χ 2 /nd f as a function of the width of the intrinsic transverse momentum distribution, 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 from the transverse momentum distributions of NuSea [42,43], R209 [41], PHENIX [40] (as shown in Fig. 6). For the calculation of χ 2 /nd f 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 /nd f 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 [49], determined from fits to inclusive DIS data from HERA which are not sensitive to intrinsic-k T .

Comments on the low-mass region
It has been observed in [17] that perturbative fixed-order calculations at O(α s ) and O(α 2 s ) 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 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 [69]), 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 middle 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]).
Our calculation thus indicates that at low energies QCD contributions beyond fixed order (O(α s ), O(α 2 s ), etc.) are important to describe the region p T /m DY ∼ 1, unlike the case of LHC energies where fixed order calculations are sufficient to describe the region p T /m Z ∼ 1. We have taken into account all-order contributions through the PB-TMD evolution formalism, and found that this allows one to describe well the transverse momentum spectra.
To sum up, the DY transverse momentum in the lowmass 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 [70], with PB-TMDs and MC@NLO (alternative methods of matching are e.g. those inspired by [12]).

Discussion
To put the results of this work in a broader context, one may start from a simple scenario in which one hopes to describe highp T dynamics by perturbative NLO calculations combined with collinear parton densities, and lowp T dynamics by non-perturbative TMDs based, in the simplest model, on intrinsic-k T Gauss distributions. One may wonder whether these two elements, NLO collinear calculations for perturbative highp T physics and intrinsic-k T TMD distributions for nonperturbative lowp T physics, are sufficient to provide a satisfactory description of the transverse momentum spectrum over all kinematic regions. The analysis of this paper illustrates that this simple approach cannot be guaranteed to give the correct physical picture in all phase space configurations. The key element which is missing in this simple approach is QCD multiple-parton radiation, and the analysis of this paper shows that (predominantly infrared) components of this radiation become essential in the region p T ∼ 1 − 10 GeV ∼ O(m DY ) of low-energy DY experiments. It also shows, more specifically, that such effects are essential for the transverse momentum spectrum, while they do not influence very much the mass spectrum integrated over transverse momenta.
If such contributions are to be included, one could imagine doing this in different manners. In this work we have done this by the PB method. This may be regarded as being well-suited to this problem, because (1) it includes multiple-parton radiation through the evolution of TMDs, (2) it incorporates the intrinsic-k T distribution as a nonperturbative boundary condition to a well-defined branching evolution equation in terms of perturbatively calculable kernels, and (3) it is matched through MC@NLO to NLO hard-scattering functions. It thus contains the three main inputs which are essential to the physical picture described above. In particular, the kernels describing multi-parton radiation through TMD evolution are given in terms of Sudakov form factors, real-emission splitting functions, and angular-ordering phase space constraints, which are important to correctly take into account infrared gluon emission.
The analysis performed in this paper leads to different conclusions from those which have appeared in the literature pointing to difficulties [17] in describing the lowmass and low-energy DY measurements and to the "q t crisis" scenario [71,72]. The analysis in this paper indicates that, provided infrared multi-parton radiation is included (e.g., through PB-TMD evolution), theoretical predictions describe low-mass and low-energy DY measurements well. It further shows that such measurements provide enhanced sensitivity to intrinsic k T compared to the case of high-energy experiments. They can thus be usefully exploited for determinations of nonperturbative TMDs. The analysis in this paper also stresses the difference between the behavior discussed above for the region p T ∼ m DY of low-energy DY experiments and the behavior in the region measured at the LHC with p T ∼ m DY ∼ 100 GeV. In the latter, no large correction is expected to purely-collinear finite-order perturbative calculations. This confirms that arguments purely based on scaling in the ratio m DY / p T are not sufficient, due to both the running of the strong coupling, and the role of infrared emission.
Other approaches would be possible as well. For instance, parton showers take into account multiple parton radiation in a manner alternative to the PB-TMD method. They can be matched to NLO matrix elements. Most parton shower Monte Carlo also model intrinsic-k T effects. In this respect, it is note-worthy that the Herwig study [45] found good agreement with DY measurements at low energy, provided parameters for the parton shower and intrinsic k T were suitably tuned, and it should be interesting to also reanalyze this in Pythia and other Monte Carlo generators. The agreement with DY data found in [45] underlines the relevance of infrared multiple emissions (taken into account, in this calculation, by showering) for the DY region of the low-energy experiments.
However, significant differences exist between the parton shower approach and the PB-TMD approach. One significant difference is that in the PB-TMD method nonperturbative TMD densities are defined and determined from fits to experimental data, which places constraints on fixed-scale inputs to evolution, while in parton showers the parton densities are not used to constrain evolution, and instead nonperturbative physics parameters are tuned. This may have an impact on the size of intrinsic-k T effects in the two approaches. On one hand, in the case of the PB method we have seen in this work that intrinsic k T q s / √ 2 with q s ∈ (300, 500) MeV provides predictions which describe well DY measurements across the energy range from NuSea √ s = 38.8 GeV to the LHC √ s = 13 TeV. On the other hand, to our knowledge it is not yet clear at present whether tuning of parton shower generators to LHC and low-energy data would result in similarly mild s-dependence of the intrinsic k T , or whether it would require a much stronger s-dependence.
A further significant difference between the shower and PB-TMD approaches is that in the shower calculation [45] the showering scale is lowered, with respect to the case of the LHC, to describe the low-energy region. In contrast, in the PB-TMD calculation of this paper the initial evolution scale is not changed, and the same starting scale μ 0 1 GeV is applied for the LHC and for the lower-energy NuSea, R209 and PHENIX experiments. We think that the investigation of these differences and their interpretation will be important questions to be examined, particularly to elucidate contributions from low-momentum regions.
Another possible approach is based on analytic CSS [12] resummation. In this formalism too the contributions from multiple soft-gluon emission, intrinsic k T and NLO hardscattering functions can be included. The formulation is however very different from that in the PB method. In particular, the matching procedure [12] (involving the so-called W and Y terms) differs from the matching used in this paper, which is of the type studied in [70]. Also the way to include intrinsic transverse momentum effects (in b or k T space) differs between CSS and PB. We expect the region p T ∼ 1-10 GeV of low-energy DY experiments to be particularly sensitive to such differences in the matching and intrinsic k T effects. We therefore think that much is to be learnt from a detailed comparison in this region.

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 Zproduction 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 /nd f ∼ 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 [49].
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.
The results which we have presented in Figs. 3 and 6, in particular, provide a new perspective on the "q t crisis" recently discussed in the literature (see e.g. contributions in Refs. [17,71,72]) with regard to measurements of transverse momentum spectra at low mass. Figure 3 illustrates that, in the kinematic region p T /m DY ∼ 1 of experiments at low center-of-mass energies √ s, hard real emission does not dominate the transverse momentum spectrum, in contrast to the case of the analogous kinematic region around the Z boson mass at the LHC. Correspondingly, NLO collinear calculations are not sufficient to describe the region of lowenergy DY measurements and multi-parton radiation contributions need to be taken into account. On the other hand, Fig. 6 shows that once the matching of NLO and multi-parton contributions is accomplished, as is done in the present study using the PB-TMD formalism, low-mass DY measurements can be well described over a broad range of center-of-mass energies √ s including the NuSea, R209 and PHENIX experiments. The matching in the present paper is carried out via PB-TMDs and MC@NLO with an approach similar to [70]. Low-mass DY data can thus be exploited to extract information on non-perturbative TMD dynamics. discussions on the PHENIX measurement. We thank A. Siodmok for pointing us to their study on Herwig. FH acknowledges the support and hospitality of DESY, Hamburg while part of this work was being done. STM thanks the Humboldt Foundation for the Georg Forster research fellowship and gratefully acknowledges support from IPM. QW and HY acknowledge the support by the Ministry of Science and Technology under Grant no. 2018YFA040390 and by the National Natural Science Foundation of China under Grant no. 11661141008.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: This is a theoretical study and has no associated experimental data.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .