Difficulties in the description of Drell-Yan processes at moderate invariant mass and high transverse momentum

We study the Drell-Yan cross section differential with respect to the transverse momentum of the produced lepton pair. We consider data with moderate invariant mass Q of the lepton pair, between 4.5 GeV and 13.5 GeV, and similar (although slightly smaller) values of the transverse momentum q_T. We approach the problem by deriving predictions based on standard collinear factorization, which are expected to be valid toward the high-q_T end of the spectrum and to which any description of the spectrum at lower q_T using transverse-momentum dependent parton distributions ultimately needs to be matched. We find that the collinear framework predicts cross sections that in most cases are significantly below available data at high q_T. We discuss additional perturbative and possible non-perturbative effects that increase the predicted cross section, but not by a sufficient amount.

energies for moderate values of the invariant mass Q and in the region q T Q. We focus on the predictions based on collinear factorization and examine their ability to describe the experimental data in this regime. We find in fact that the predicted cross sections fall significantly short of the available data even at the highest accessible values of q T . We investigate possible sources of uncertainty in the predictions based on collinear factorization, and two extensions of the collinear framework: the resummation of high-q T threshold logarithms, and transverse-momentum smearing.
None of these appears to lead to a satisfactory agreement with the data. We argue that these findings also imply that the Drell-Yan cross section in the "matching regime" q T Q is presently not fully understood at fixed-target energies.
We note that a similar problem has been reported in [31] for semi-inclusive deep inelastic scattering (SIDIS) processes in the region of large transverse momenta, where large disagreements have been observed also in this case between fixed-order calculations and experimental data. The discrepancies we report here arguably appear more serious since the calculation of the Drell-Yan cross section relies on the very well constrained PDFs, while SIDIS is also sensitive to the comparably more poorly known fragmentation functions.

II. MOTIVATION: FROM TMDS TO THE MATCHING REGIME
As described in the Introduction, the regime q T Q may be addressed in terms of TMD factorization, and numerous studies using fixed-target Drell-Yan data have been carried out [6][7][8][9][10][11][12][13][14][15], which however only address the region q T 1.5 GeV and do not make any attempt to perform a matching to a fixed-order calculation at higher q T . Indeed, extending the description to the whole q T -spectrum is a delicate task. To understand the related issues, it is worth to summarize here the basic ideas behind the most common matching procedures (for detailed expositions, we refer the reader to dedicated studies, e.g., [26,32]).
The low-q T formula for the cross section, which embodies TMD physics, has the following expression: where W pert contains soft gluon resummation and W NP the non-perturbative terms. The observed transverse momentum distribution is thus a convolution of the two contributions. Since perturbative calculations would hit the Landau pole at large values of b, one common solution is to freeze the impact parameter b beyond a threshold b max , by introducing the function b * (b), constructed in such a way that b * b when b b max , and b * = b max when b > b max . With increasing q T , one expects a smooth transition from TMD physics to collinear factorization. A common way to describe this transition is the following: a correction term (so called "Y-term") is added to Eq. (1), in order to approximate the sub-leading (in powers of q T /Q) contributions that are not present in the resummed formula. It is given by the difference between the fixed-order and asymptotic cross sections: where the asymptotic piece is obtained by isolating the terms in the fixed-order expression that are most divergent for q T /Q → 0. In an ideal situation, at some point as q T increases towards Q, the asymptotic term in Eq. (2) cancels with the W term, so that the sum W + Y approaches the fixed-order cross section (see, e.g., Sec. 1.4 of [32]). The matching procedure can pose serious problems when Q is not very high, as was shown in [33] for the case of SIDIS. 2 The observed problems can be summarized as follows: the high-q T tail of the TMD formula shows sensitivity to the non-perturbative parameters and to the details of the b * function, preventing a proper cancellation with the Y term. It is straightforward to check that this behavior is also present in Drell-Yan at fixed-target kinematics. To give an example, in Fig. 1 we show the effect of extrapolating the TMD fitted in [13] to high q T . Although the asymptotic curve drops very rapidly at some point, signaling that O (q T /Q) corrections should become dominant, the TMD extends far beyond, owing to the non-perturbative Sudakov contribution. We note that also the behavior of the W pert and W NP for b → 0 is expected to play a role here. All in all, while on the one side the shape of the data seems to suggest that TMD physics is indeed involved in some form up to transverse momenta as high as 2.5 GeV, one has to admit at the same time that presently there is not a good understanding of the TMD formalism in this region. The matching procedure is afflicted by large uncertainties, and the TMD tail is largely affected by nonperturbative elements, such as the functional form of b * (see figure for details).  [36]). bmax is kept fixed at 1.123 GeV −1 . The asymptotic curve is also plotted (at LO, to be consistent with the fit). Right: matched curve obtained from the same TMD, with the procedure described in [34]. The damping functions are taken as in Sec. IX of the same article. Data are taken from [37].

�=��� ���
To avoid these (and other) problems, the authors in [34] proposed a modified matching procedure. Without entering into details, we only underline that this procedure forces the use of pure fixed-order calculation at intermediate values of q T , by suppressing the tail of the TMD cross section with a damping function. Just to give a qualitative example, in Fig. 1 we show the effect of using the same damping function in our case. An alternative approach to suppress the TMD contribution at high transverse momentum was proposed also in [35].
In conclusion, the Drell-Yan q T spectra at low invariant mass are presently not understood beyond the region q T Q typical of TMD fits. In the following, we will approach the problem from high q T ∼ Q, where collinear factorization is expected to offer a suitable framework for describing the cross section. Undoubtedly the collinearfactorized cross section will be an important ingredient for a better understanding of the regime q T Q, where it will be especially important for carrying out the proper matching of the resummed cross section.

III. COLLINEAR FACTORIZATION AND COMPARISON TO FIXED-TARGET DATA
In this section we show the comparison of fixed-order perturbative QCD calculations to Drell-Yan data from Fermilab, CERN and RHIC experiments, mainly for proton-proton collisions. The center-of-mass energies of the experiments taken into account lie in the range 20 GeV ≤ √ s ≤ 60 GeV (except for RHIC, where √ s=200 GeV), while the invariant mass of the Drell-Yan lepton pair lies in the range 4.5 ≤ Q ≤ 13.5 GeV. For all our theoretical predictions, we use the DYqT [23,24] and CuTe [25] codes, obtaining completely equivalent results for the fixed-order differential cross sections, at both LO QCD (O (α s )) and NLO QCD O α 2 s . These codes also provide an all-order resummation of logarithms in q T /Q in the cross section, which become relevant toward low q T . This enables us to study the asymptotic expansion of the resummed result, which we will make use of below. We note that we have also performed cross-checks using the numerical codes of Refs. [21] and [22]. Throughout this paper, the CT14 PDF set [38] will be our default choice.

E866
The E866/NuSea experiment [39] was a fixed-target Drell-Yan experiment designed to measure the internal structure of the nucleon, in particular the asymmetry of down and up antiquarks in the sea, using di-muon events originating from the collision of an 800-GeV proton beam with hydrogen and deuterium targets ( √ s = 38.8 GeV). The measurement of the q T -distribution of the muon pair is presented in [40], a Fermilab PhD thesis, and results are given in terms of the differential cross section: (3) Data are reported for different bins in x F = 2p L / √ s, ranging from −0.05 to 0.8, and are integrated over different ranges in the invariant mass Q of the muon pair.
The comparison of our LO and NLO theoretical calculations with the experimental data is shown in Fig. 2 for the bin 0.15 ≤ x F ≤ 0.35 and for the invariant mass range 4.2 GeV ≤ Q ≤ 5.2 GeV. The lower part of the plot shows the ratio (data-theory)/theory. The error margins of the data points correspond to the sum in quadrature of statistical and systematic uncertainties, including also an overall normalization uncertainty of 6.5%, as indicated in [40]. Our theoretical predictions are computed at the average Q value and x F of each bin (Q = 4.7 GeV and x F = 0.25 in the case of Fig. 2). The left plot of Fig. 2 shows the comparison of the experimental data with NLO QCD O α 2 s predictions for central values of the factorization and renormalization scales, µ R = µ F = Q. The 90% confidence interval of the CT14 PDF set [38] is included in the plot, but the corresponding variation is barely visible.
An immediate observation from Fig. 2 is that the NLO cross section is below the E866 data at high transverse momenta, q T 3 GeV, even within the relatively large uncertainties that the data have here. The NLO cross section falls below the data even much more severely at lower q T closer to the "matching regime" with TMD physics, where the experimental uncertainties are much smaller. This provides further evidence to our observation above that this regime is presently not well understood theoretically. At the same time we emphasize that data from [40], integrated over q T , are in good agreement with theoretical predictions and are commonly used in global PDF fits [41,42] (see, for instance, Section 5.1 of [40], where the only relevant discrepancy concerns the lowest mass point ( Q 4.4 GeV) for 0.05 < x F < 0. ). This suggests that TMD physics may be the main player for the cross section up to relatively high q T , since the tail at very large q T makes only a small contribution to the cross section.
The right plot of Fig. 2 shows the effect of varying the renormalization and factorization scales independently in the range Q/2 < µ R , µ F < 2Q, both for the LO QCD (O (α s )) and the NLO QCD O α 2 s calculation. The fact that, for q T 2.5 GeV, the NLO uncertainty band overlaps with (and is eventually included in) the LO uncertainty band provides some indication that perturbation theory is well-behaved for this process. On the other hand, we also observe that the NLO scale uncertainty band is only marginally more narrow than the LO one.
We have also considered different PDF choices (CTEQ 10 [43], NNPDF 2.3 [44] and MSTW2008 [41]), obtaining very similar results: the different curves lie within the uncertainty bands shown in the right plot of Fig. 2. Such a mild PDF dependence was expected, since the PDFs are well constrained and have small uncertainties in the xrange probed in this process. We conclude that PDF uncertainties (unless they are grossly underestimated by the parameterizations) cannot explain the discrepancy between theory and data at high q T .
The comparison between data and theory for other x F bins ( Fig. 3) and for a different invariant mass range (Fig. 4) gives the same qualitative results. The upper part of each plot contains the NLO QCD O α 2 s prediction (blue) with its uncertainty band obtained through the customary scale variation (Q/2 < µ R , µ F < 2Q) around the central value Q of the invariant mass range. The lower part of each plot again shows the ratio (data-theory)/theory. We also plot the asymptotic expansion of the resummed calculation (red lines). The asymptotic result coincides with the fixed order prediction in the region of very low transverse momenta, but it becomes very small (and eventually negative) with increasing q T . We show the asymptotic piece in order to obtain a rough guide concerning the region where the fixed-order calculation may start to become reliable [32]: ideally, when q T is large enough that the difference between the fixed-order and asymptotic calculations (the so-called "Y term") exceeds the full ("W + Y ") cross section, one should switch from W + Y to the fixed-order result to obtain more reliable predictions. This occurs for q T values around 1-2 GeV in the present case. Figures 3 and 4 show the same qualitative features seen above: the overall agreement between theory and high-q T data is poor. In general, the disagreement between data and theoretical predictions seems to become worse with increasing Feynman-x F and to be only mildly dependent on the invariant mass Q of the lepton pair.

R209
The R209 experiment [45,46] (two proton beams colliding at a center-of-mass energy of √ s = 62 GeV) was carried out at the CERN ISR (Intersecting Storage Rings) to search for new particles and test scaling models. The differential cross section dσ/dq 2 T for the production of a muon pair with transverse-momentum q T is reported in [47] for the invariant mass range 5 GeV < Q < 8 GeV. The low transverse momentum part of these data has been included in extractions of TMDs [7,9]. Studies of the whole q T spectrum can be found in [48,49].
Comparisons of our NLO results to the R209 data are shown in Fig. 5. Again NLO is below the data at high q T , although the discrepancy is not as statistically significant in this case as for the E866 data. We note that a GeV, including a 90% confidence interval from the CT14 PDF set [38]. Right: LO QCD and NLO QCD theoretical uncertainty bands obtained by varying the renormalization and factorization scales independently in the range Q/2 < µR, µF < 2Q.
FIG. 3: E866: comparison between experimental data and NLO QCD predictions for different xF bins. We also show the low-qT asymptotic part of the cross section. For details, see text. similar gap between data and theory was reported in [49] in the context of a LO calculation. There, the so-called "k T -factorization" formalism was claimed to account for the discrepancy. In contrast, in [48] the W + Y formalism was reported to match the data over the whole q T range.

E288
The E288 experiment [37] measured the invariant cross section Ed 3 σ/d 3 q, at fixed photon rapidity, for the production of µ + µ − pairs in the collision of a proton beam with a fixed target composed of either Cu or Pt. The measurements were performed using proton incident energies of 200, 300 and 400 GeV, producing three different data sets. The respective center of mass energies are √ s = 19.4, 23.8, 27.4 GeV. Our results are shown in Figs. 6,7,8,9. The comparison to data shows the same features as before. We have tested the importance of nuclear effects by computing the cross sections also with the nCTEQ15 [50] and CT14 [38] nuclear PDFs. These turn out to lead to almost indistinguishable results. We note that the low transverse momentum part of the E288 data has been used for extractions of TMDs [7,9,10,13,14]. E605 We also consider the set of measurements of Ed 3 σ/d 3 q in the E605 [51] experiment, extracted from an 800-GeV proton beam incident on a copper fixed target ( √ s = 38.8 GeV). Results at fixed x F = 0.1 are shown in Fig. 10. The low transverse momentum part of these data has also been included in extractions of TMDs [7,9,13].

PHENIX
Finally, we also compare to the recent measurement [52] performed by the PHENIX collaboration at the Relativistic Heavy Ion Collider in pp collisions at √ s=200 GeV. The experimental points are taken from Fig. 33 of [52] and compared to LO QCD and NLO QCD, including theoretical uncertainties, in Fig. 11. The asymptotic expansion of the W term to NLO is also shown. Evidently, the comparison between NLO and the data is overall satisfactory in this case. It thus appears that there is a qualitative difference between the fixed-target and collider regimes.

IV. THRESHOLD RESUMMATION
As we have seen in Fig. 2, the NLO corrections to the q T -differential cross sections are quite sizable. It is therefore important to investigate in how far beyond-NLO perturbative corrections might be relevant for obtaining a better agreement with the data. For the kinematics relevant for the Fermilab and CERN experiments, the invariant mass and transverse momentum of the Drell-Yan pair are such that the production is relatively close to partonic threshold, where a new class of logarithms (separate from that mentioned above at low q T ) arises. The summation of these logarithms to all orders is known as threshold resummation. We note that large corrections from threshold resummation have been found previously in purely hadronic single-inclusive processes such as pp → πX [53,54], which motivates a corresponding study for the high-q T Drell-Yan cross section pp → γ * X → + − X that will be carried out in this section. The relevant formalism has been developed in Refs. [55][56][57][58][59], although in most of these papers only fixed-order (NNLO) expansions of the resummed cross sections have been considered, and in [60] for the closely related high-q T Higgs production cross section. We follow here the approach taken in the latter reference.

Factorized cross section and Mellin moments
For simplicity we will focus here just on the transverse momentum distribution of the lepton pair and integrate over the full range of allowed rapidities of the virtual photon. We therefore consider where FIG. 9: Additional plots for E288: experimental data vs. NLO QCD predictions for η=0.03 and different invariant mass bins.  with m T ≡ Q 2 + q 2 T the transverse mass. The integrated differential cross section for h 1 h 2 → − + X may be written in factorized form as

�=���-��� ���
where σ 0 = 4πα 2 /(9Q 2 ), f a/h1 and f b/h2 are the PDFs, and whereŝ = sx 1 x 2 is the partonic center-of-mass energy squared. In the second line we have written out the variables that the dimensionless hard-scattering functions ω ab may depend on. It is convenient to write the kinematical arguments of ω ab aŝ Note that in terms of these we have q T / √ŝ =ŷ T r/(1+r) and Q 2 /ŝ =ŷ 2 T (1−r)/(1+r). In (6) we have also introduced the corresponding hadronic variable From Eq. (5) we see that y T ≤ 1. Likewise, we also haveŷ T ≤ 1, which immediately leads to the limits on x 1 and x 2 given in Eq. (6). The rapidity-integrated cross section thus takes the form of a convolution of the hard-scattering functions ω ab with the PDFs. The perturbative expansion of the ω ab reads There are two LO partonic channels, qq → γ * g and qg → γ * q. As before, µ F and µ R in (6) denote the factorization and renormalization scales. Forŷ T → 1 the partonic center-of-mass energy is just sufficient to produce the lepton pair with mass Q and transverse momentum q T . Therefore,ŷ T = 1 sets a threshold for the process. As is well known [61,62], the partonic cross sections receive large logarithmic corrections near this threshold. At the kth order of perturbation theory for the ω ab , there are logarithmically enhanced contributions of the form α k s ln m (1 −ŷ 2 T ), with m ≤ 2k. These logarithmic terms are due to soft and/or collinear gluon radiation and dominate the perturbative expansion when the process is kinematically close to the partonic threshold. We note thatŷ T becomes especially large when the partonic momentum fractions approach their lower integration limits. Since the PDFs rise steeply towards small argument, this enhances the relevance of the threshold regime, and the soft-gluon effects are relevant even when the hadronic center-of-mass energy is much larger than the produced transverse mass and transverse momentum of the final state. In the following, we discuss the resummation of the large logarithmic corrections to all orders in α s .
The resummation of the soft-gluon contributions is carried out in Mellin-N moment space, where the convolutions in Eq. (6) between parton distributions and subprocess cross sections factorize into ordinary products. We take Mellin moments of the hadronic cross section in y 2 T : where the corresponding moments of the partonic hard-scattering functions arẽ and where thef a,b (N + 1, µ 2 F ) are the moments of the parton distributions. The threshold limitŷ 2 T → 1 corresponds to N → ∞, and the leading soft-gluon corrections now arise as terms ∝ α k s ln m N , m ≤ 2k. The NLL resummation procedure we present here deals with the "towers" α k s ln m N for m = 2k, 2k − 1, 2k − 2.

Resummation to NLL
In Mellin-moment space, threshold resummation results in exponentiation of the soft-gluon corrections. In case of lepton pair production at high q T , the resummed y-integrated cross section reads [60]: where for simplicity we have suppressed the arguments of all functions other than the Mellin variable N . Each of the "radiative factors" ∆ a,b N , J c N , ∆ (int)ab→γ * c N is an exponential. The factors ∆ a,b N represent the effects of soft-gluon radiation collinear to initial partons a and b. The function J c N embodies collinear, soft or hard, emission by the parton c that recoils against the lepton pair. Large-angle soft-gluon emission is accounted for by the factors ∆ (int)ab→γ * c N , which depend on the partonic process under consideration. Finally, the coefficients C ab→γ * c contain N -independent hard contributions arising from one-loop virtual corrections and non-logarithmic soft corrections. The structure of the resummed expression is similar to that for the large-q T W -boson production cross section [55][56][57][58] or, in the massless limit, to that for prompt-photon production in hadronic collisions [63,64].
The expressions for the radiative factors are where Q 2 0 ≡ q T (q T + m T ). Each of the coefficients C = A a , B a , D ab→γ * c is a power series in the coupling constant α s , C = ∞ i=1 (α s /π) i C (i) . The universal LL and NLL coefficients A (1) a are well known [65,66]: with where C g = C A = N c = 3, C q = C F = (N 2 c − 1)/2N c = 4/3, γ q = −3/2C F and γ g = −(11C A − 2N f )/6. The process-dependent coefficients D ab→γ * c may be obtained as for the Higgs production cross section considered in [60]: The coefficient is evidently just proportional to a combination of the color factors for each hard parton participating in the process. This simplicity is due to the fact that there is just one color structure for a process with only three external partons. The final ingredients for the resummed cross section in (12) are the lowest-order partonic cross sections in Mellinmoment space,ω (0) ab (N ), and the coefficients C ab→γ * c . The expressions for the former are presented in Appendix A. At NLL accuracy, we only need to know the first-order terms in the expansion C ab→γ * c = 1 + These coefficients may be obtained by comparison to the full NLO results given in Ref. [18] (see also [58]), after transforming to moment space. Our explicit results for the one-loop coefficients C (1) ab→γ * c are given in Appendix B. In order to organize the resummation according to the logarithmic accuracy of the Sudakov exponents we expand the latter to NLL as [67] ln ∆ a N (α s (µ 2 R ), where γ E is the Euler constant. The LL and NLL functions h (1,2) and f (1,2) are given in Appendix C.

Matching and inverse Mellin transform
When performing the resummation, one wants to make full use of the available fixed-order cross section, which in our case is NLO (O(α 2 s )). Therefore, one matches the resummed result to the fixed-order expression. This is achieved by expanding the resummed cross section to O(α 2 s ), subtracting the expanded result from the resummed one, and adding the full NLO cross section: whereω (res) ab is the resummed cross section for the partonic channel ab → γ * c as given in Eq. (12). In this way, NLO is fully taken into account, and the soft-gluon contributions beyond NLO are resummed to NLL. The procedure avoids any double-counting of perturbative orders.
Since the resummation is achieved in Mellin-moment space, one needs an inverse Mellin transform in order to obtain a resummed cross section in y T space. This requires a prescription for dealing with the Landau poles at λ = 1/2 and λ = 1 in Eqs. (C1)-(C2) arising from the singularity in the perturbative strong coupling constant at scale Λ QCD . We employ the "Minimal Prescription" developed in Ref. [67], for which one uses the NLL expanded forms Eq. (17) and Eqs.(C1)-(C2), and chooses a Mellin contour in complex-N space that lies to the left of the poles at λ = 1/2 and λ = 1 in the Mellin integrand.

Numerical results
Numerical results for the above formalism are shown in Figs. 12   resummed cross section (which above was determined after integration over all y) by the ratio of NLO cross sections integrated over the y-bin used in experiment and integrated over all y, respectively: This approximation assumes that the rapidity dependence is similar at NLO and in the resummed case, an expectation that was confirmed in Ref. [68] for the closely related prompt-photon cross section. We first notice from Figs. 12 and 13 that the NLO expansion of the resummed formula (black dashed curve) accurately reproduces the NLO result (blue solid curve, with uncertainty bands). This provides some confidence that threshold resummation correctly describes the dominant parts of the cross section to all orders, and that subleading contributions not addressed by resummation are reasonably small. In the left part of Fig. 12 we also show the scale uncertainty band for the NLL matched result (red dot-dashed curve), which is barely broad enough to be visible. Evidently, resummation leads to a strong reduction in scale dependence, as one would expect from a result that incorporates the dominant contributions to the cross section at all orders.
Overall, we find a further significant increase of the cross section due to NLL resummation, with respect to the NLO results shown in Sec. III. The enhancement is more pronounced for the case of E288 than for E866 since, for a given Q, at E288 energy one is closer to threshold because of the lower c.m.s. energy. However, despite the increase, the NLL result unfortunately still remains well below the E288 and E866 experimental data at high q T . We thus conclude that NLL high-q T threshold resummation is not able to lead to a satisfactory agreement with the data.

V. INTRINSIC-kT SMEARING AND POWER CORRECTIONS
The factorized cross section given in Eq. (6) receives corrections that are suppressed by inverse powers of Q ∼ q T . Little is known so far about the structure and size of such power corrections for the high-q T Drell-Yan cross section. It is an interesting question whether the discrepancies between perturbative predictions and the high-q T experimental data seen above might be explained by power corrections. We will try here to address this question from a phenomenological point of view.
As a simple way of modeling power corrections we estimate below the impact of a non-perturbative partonic "intrinsic" transverse momentum k T on the Drell-Yan q T spectrum. Such an "intrinsic-k T smearing" is a phenomenological model that has been invoked in the early literature in cases where collinear factorization was found to underestimate transverse momentum spectra, like for inclusive prompt photon and pion production in hadronic collisions (see for instance [8,69,70]). For inclusive processes such as these and the high-q T Drell-Yan process considered here, no general factorization theorem is known that would extend to arbitrary kinematics of the partonic process. For prompt photons, factorization has been established, however, for near-threshold kinematics and low k T in the framework of the "joint resummation" formalism [71][72][73], and for high-energy (small-x) dynamics [74]. A technical challenge for all these approaches is the potential for an artificial singularity when the total transverse momentum of the initial state partons is comparable to the observed transverse momentum. A method for dealing with this issue was proposed in Ref. [75] and found to give rise to power corrections to the cross section. A full treatment of the Drell-Yan cross section may require implementation of perturbative joint resummation along with a study of corrections in inverse powers of Q or q T . Rather than pursuing this elaborate framework, for the purpose of obtaining a simple estimate of the potential size of such higher-order perturbative and power-suppressed non-perturbative effects, we resort to an implementation of a simple model of intrinsic-k T smearing that will be described now.

Overview of the formalism
The collinear factorization formula for the process h 1 h 2 → γ * X may be adapted from Eq. (6) and reads at LO (O (α s )): where as before the f a/h (x a , Q 2 ) are the usual collinear PDFs for partons a = q,q, g in hadron h. If one allows the incoming partons to have a small transverse momentum k T , Eq. (20) becomes [8]: where the functions F a/h are a generalization of the PDFs, including a dependence on transverse momentum. Notice that the partonic Mandelstam invariants must be modified with the inclusion of k T , and consequently a factor s/(x a x b s) must be inserted to account for the modification of the partonic flux (see Appendix A of [8]). The modification of the partonic four-momenta is most often done according to two criteria: (1) the partons remain on-shell: p aµ p µ a = 0, and (2) the light-cone momentum fractions retain the usual meaning, e.g.: x a = p + a /P + a . This leads to the following choice, in terms of Minkowski components [8,76]: and likewise for the other parton's momentum. Note that we use LO cross sections in Eq. (21) since a higher-order formulation is not really warranted for our simple model. As mentioned above, the framework must become unreliable when k aT or k bT become of the order of the observed transverse momentum, and arguably well before. Large values of k aT can make the partonic Mandelstam in the denominators of the LO hard-scattering cross sections unphysically small. In [8], the following condition was chosen to limit the size of, for example, k aT : This ensures that each parton moves predominantly along the direction of its parent hadron, and that its energy does not exceed the hadron's energy. However, for √ s 40 GeV (E866 and E605 experiments), this condition implies that k aT may still reach values as high as 20 GeV. In our numerical analysis we therefore prefer to introduce an additional cutoff k T max on both k aT and k bT and will test the dependence of the results on this cutoff.
For the generalized PDFs in Eq. (21), the most common choice is where k 2 T is independent of flavor 3 and momentum fraction x a , but does depend logarithmically on Q 2 because of soft gluon radiation. Instead of Eq. (24), one could also consider using the transverse momentum dependent PDFs extracted from the low-q T spectra of Drell-Yan experiments (as given for instance in Refs. [7,9,13,14]). However, these functions show a non-negligible tail at large k T , where they lose physical meaning. Hence, if they are used inside a convolution such as Eq. (21), the result will strongly depend on the choice of the cutoff k T max , since the integrations (21) include contributions from this tail. This dependence will be mostly unphysical and is, in fact, precisely a manifestation of the artificial singularity arising in the partonic scattering functions at really large k aT and k bT . For this reason, we stick with Eq. (24); however, we tune k 2 T to the width of the TMD PDFs taken from [13], evolved to the given Q 2 . This is shown in Fig. 14 where the dashed lines show the evolved TMD of Ref. [13], evolved to Q= 4.7 GeV, normalized by dividing by its integral over d 2 k T . We compare it to a pure Gaussian with a width tuned in such a way that the two distributions become very similar, except for the high-k T tail. This "equivalent Gaussian" turns out to have a width of k 2 T = (0.95 GeV) 2 . It is this Gaussian that we use for our numerical studies presented below.
Our choice of an x-independent Gaussian width in Eq. (24) is motivated by the fact that the x-dependence of k 2 T is still not well constrained in the present TMD fits [13]. Different parametrizations have been proposed in the literature [7,9], including also x-independent choices [10,14,77]. A dependence of k 2 T on x is a natural feature in the joint resummation formalism [75]. In any case, for the mostly exploratory study presented here, an x-independent value of k 2 T appears adequate. Since our goal is to give an upper limit for the k T -smearing effects, we use the largest value of k 2 T found in [13] (see Fig. 10 there), which occurs at x = 0.06.

Numerical results
In Fig. 15 we show the effect of k T -smearing, Eq. (21), for E866 kinematics. The k 2 T of the Gaussian is taken as in Fig. 14. The impact of smearing on the cross section overall remains mild, as long as the cutoff k T max is chosen below 2 GeV. Especially the regime q T Q is only little affected by k T -smearing. We conclude that, although k T -smearing does somewhat improve the comparison with the data, its effects do not appear to be sufficiently large to lead to a satisfactory agreement. We note that at lower c.m.s. energies as relevant for E288, one is forced to choose smaller cutoffs since the reach in q T is more limited in these cases.

VI. CONCLUSIONS
We have shown that theoretical predictions based on fixed-order perturbation theory fail to describe Drell-Yan data from Fermilab and CERN ISR at large values q T ∼ Q of the transverse momentum of the lepton pair, the experimental cross sections being significantly larger than the theoretical ones. This is the region where collinearfactorized perturbation theory is expected to accurately describe the cross section. This disagreement is observed for several experiments, and across a range of different kinematics in x F , y and Q, although admittedly the experimental uncertainties are in some cases quite large.
We have on the other hand found an essentially satisfactory agreement between perturbative calculations and experimental points in the case of PHENIX data taken at √ s = 200 GeV, suggesting that the disagreement is present only in the fixed-target regime. Indeed, at yet higher energies, ATLAS Drell-Yan data ( √ s = 8 TeV) have been shown to be consistently described by NNLO QCD supplemented with NNNLL resummation (see, for instance, Fig. 10 and Fig. 11 in [78]), even though some tension is still present in the lowest invariant mass bins (see Fig. 18 in [79]).
Barring the possibility of sizable normalization uncertainties in the experiments, it is important to identify the theoretical origins of the discrepancies observed in the fixed-target regime. We have first implemented perturbative threshold resummation and found that it improves the situation somewhat; a significant discrepancy remains, however. This leaves the investigation of power-suppressed corrections, which we have modeled by implementing a simple Gaussian intrinsic-k T smearing into the LO cross section. We find that this again helps somewhat, but does not lead to a satisfactory description of the data. Ultimately, a more detailed study of power corrections may be required in this case. Generically, on the basis of resummed perturbation theory [75], one would expect even power corrections of the form λ 2 / Q 2 (1 − y 2 T ) 2 , possibly modified by logarithms, where λ is a hadronic mass scale. Given the kinematics of the experiments, it is hard to see how such corrections could become of the size needed for an adequate description of the data.
Our findings are in line with those reported for the SIDIS cross section in Ref. [31]. We close by stressing the importance of obtaining a thorough understanding of the full Drell-Yan and SIDIS q T -spectra in the fixed-target regime. Low-q T Drell-Yan and SIDIS cross sections measured at fixed-target experiments are a prime source of information on TMDs. At present, the theoretical description for the important "matching regime" around q T = 2 GeV is not robust, as we have argued. Given the shape of the experimental spectra, it appears that TMD physics may extend to such large q T and may well remain an important ingredient even beyond. This view is corroborated by the fact that the q T -integrated Drell-Yan cross section is well described by fixed-order perturbation theory at these energies. In any case, a reliable interpretation of data in terms of TMDs, including the matching to collinear physics, is only possible if the cross sections are theoretically understood over the full transverse-momentum range, which includes the regime of q T ∼ Q we have addressed here.
The explicit expressions for the LL and NLL functions in Eq. (17) are: and are the first two coefficients of the QCD β function.