Production of Drell--Yan lepton pairs in hadron collisions: transverse-momentum resummation at next-to-next-to-leading logarithmic accuracy

We consider the transverse-momentum (q_T) distribution of Drell--Yan lepton pairs produced in hadron collisions. At small values of q_T, we resum the logarithmically-enhanced perturbative QCD contributions up to next-to-next-to-leading logarithmic accuracy. At intermediate and large values of q_T, we consistently combine resummation with the known next-to-leading order perturbative result. All perturbative terms up to order alpha_S^2 are included in our computation which, after integration over q_T, reproduces the known next-to-next-to-leading order result for the Drell--Yan total cross section. We show and discuss the reduction in the scale dependence of the results with respect to lower-order calculations, estimating the corresponding perturbative uncertainty. We present a preliminary comparison with Tevatron Run II data.


Introduction
The hadroproduction of the vector bosons W and Z/γ * , also known as the Drell-Yan (DY) process [1], is the process to which parton model ideas (previously developed for deep inelastic lepton-hadron scattering) were first applied in the context of hard-scattering processes in hadronhadron collisions.
At high-energy hadron colliders, such as the Tevatron and the LHC, vector bosons are produced with large rates and with relatively-simple experimental signatures. The vector boson production process is thus relevant for various reasons. It is important for detector calibration; it provides us with strong tests of perturbative QCD and, in particular, it gives stringent information on the parton densities of the colliding hadrons; it represents an important background for new-physics searches. Owing to these reasons, it is essential to have accurate theoretical predictions for vector boson production cross sections and related kinematical distributions.
These predictions are based on perturbative QCD and are obtained as power series expansions in the strong coupling α S . The total cross section [2] and the rapidity distribution of the vector boson [3] are known up to the next-to-next-to-leading order (NNLO) in QCD perturbation theory. Fully exclusive NNLO calculations, including the leptonic decay of the vector boson, are also available [4,5,6]. Electroweak corrections up to O(α) have been computed both for W [7] and Z/γ * production [8].
In this paper we consider the transverse-momentum (q T ) spectrum of the vector boson. The q T spectra of the W and Z bosons are particularly important since the uncertainties in their shape directly affect the measurement of the W mass. In the large-q T region (q T ∼ m V ), where the transverse momentum is of the order of the vector boson mass m V , perturbative QCD calculations based on the truncation of the perturbative series at a fixed order in α S are theoretically justified. In this region, the QCD radiative corrections are known up to the next-to-leading order (NLO) [9,10,11]. Nonetheless the bulk of the vector boson events is produced in the small-q T region (q T ≪ m V ), where the convergence of the fixed-order expansion is spoiled by the presence of large logarithmic terms, α n S ln m (m 2 V /q 2 T ). To obtain reliable predictions, these logarithmically-enhanced terms have to be systematically resummed to all perturbative orders [12]- [20]. The resummed and fixed-order calculations at small and large values of q T can then be consistently matched at intermediate values of q T , to obtain QCD predictions for the entire range of transverse momenta.
We use the transverse-momentum resummation formalism proposed in Refs. [20,21]. The formalism is valid for a generic process in which a high-mass system of non strongly-interacting particles is produced in hadron-hadron collisions. The method has so far been applied to the production of the Standard Model (SM) Higgs boson [22,21,23], single vector bosons [24], W W [25] and ZZ [26] pairs, slepton pairs [27], and DY lepton pairs in polarized collisions [28]. The study of Ref. [24] is mainly based on next-to-leading logarithmic (NLL) resummation at small q T and on the leading-order (LO) calculation at large q T . In this paper we extend the analysis and the results of Ref. [24], combining the most advanced perturbative information that is available at present: next-to-next-to-leading logarithmic (NNLL) resummation at small q T and the NLO calculation at large q T . Other phenomenological studies of the vector boson q T distribution, which combine resummed and fixed-order perturbative results at various levels of theoretical accuracy, can be found in Refs. [29].
The paper is organized as follows. In Sect. 2 we briefly review the resummation formalism of Refs. [20,21] and its application to vector boson production. In Sect. 3 we present numerical results for Z/γ * production, and we comment on their comparison with the Tevatron Run II data [30,31]. We also study the scale dependence of our results to the purpose of estimating the corresponding perturbative uncertainty. In Sect. 4 we summarize our results.

Transverse-momentum resummation
We briefly recall some of the main points of the transverse-momentum resummation formalism of Refs. [20,21]. Here we consider the specific case of DY lepton pair production, i.e. the production of a vector boson V (V = W + , W − , Z/γ * ) that subsequently decays in a lepton pair.
The inclusive hard-scattering process is where h 1 and h 2 are the colliding hadrons with momenta p 1 and p 2 , V is the vector boson (which decays in the lepton pair l 1 , l 2 ) with invariant mass M and transverse momentum q T , and X is an arbitrary and undetected final state.
According to the QCD factorization theorem the q T differential cross section dσ V /dq 2 T can be written as (2) where f a/h (x, µ 2 F ) (a = q,q, g) are the parton densities of the colliding hadron h at the factorization scale µ F , dσ V ab /dq 2 T are the perturbative QCD partonic cross sections, s (ŝ = x 1 x 2 s) is the square of the hadronic (partonic) centre-of-mass energy, and µ R is the renormalization scale.
In the region where q T ∼ M (in practice, we always consider the case in which M is close to the mass m V of the vector boson), the QCD perturbative series is controlled by a small expansion parameter, α S (M), and fixed-order calculations are theoretically justified. In this region, the QCD radiative corrections are known up to next-to-leading order (NLO) [9].
In the small-q T region (q T ≪ M), the convergence of the fixed-order perturbative expansion is spoiled by the presence of powers of large logarithmic terms, α n S ln m (M 2 /q 2 T ). To obtain reliable predictions these terms have to be resummed to all orders.
We perform the resummation at the level of the partonic cross section, which is decomposed as The first term on the right-hand side contains all the logarithmically-enhanced contributions, which have to be resummed to all orders in α S , while the second term is free of such contributions and can thus be evaluated at fixed order in perturbation theory. Using the Bessel transformation between the conjugate variables q T and b (b is the impact parameter), the resummed component dσ (res.) V ab can be expressed as where J 0 (x) is the 0th-order Bessel function. Considering the Mellin N-moments W N of W with respect to the variable z = M 2 /ŝ at fixed M, the resummation structure of W V ab, N can be organized in exponential form † were we have defined the logarithmic expansion parameter L ≡ ln(Q 2 b 2 /b 2 0 ), and b 0 = 2e −γ E (γ E = 0.5772... is the Euler number). The scale Q (Q ∼ M), which appears on the right-hand side of Eq. (5), is the resummation scale [21]. Although W V N (i.e., the product H V N × exp{G N }) does not depend on Q when evaluated to all perturbative orders, its explicit dependence on Q appears when W V N is computed by truncation of the resummed expression at some level of logarithmic accuracy (see Eq. (6) below). Variations of Q around M can thus be used to estimate the size of yet uncalculated higher-order logarithmic contributions.
The universal ‡ form factor exp{G N } contains all the terms that order-by-order in α S are logarithmically divergent as b → ∞ (or, equivalently, q T → 0). The resummed logarithmic expansion of the exponent G N is defined as follows: where the term L g (1) collects the leading logarithmic (LL) contributions, the function g N includes the next-to-leading leading logarithmic (NLL) contributions [17], g N controls the NNLL terms [32,33,34] and so forth. The explicit form of the functions g (1) , g (2) N and g (3) N can be found in Ref. [21]. The process dependent function H V N does not depend on the impact parameter b and it includes all the perturbative terms that behave as constants as b → ∞. It can thus be expanded in powers of α S = α S (µ 2 R ): where σ (0) V is the partonic cross section at the Born level. The first-order coefficients H in Eq. (7) are known since a long time [32], while the second-order coefficients H qq←ab,N were computed only recently [5]. † For the sake of simplicity we consider here only the case of the diagonal terms in the flavour space of the partonic indices a, b. For the general case and a detailed discussion of the resummation formalism, we refer to Ref. [21]. ‡ The form factor does not depend on the type of produced vector boson. More generally, all the hard-scattering processes initiated by quark-antiquark annihilation have the same form factor.
Within a straightforward ('naive') implementation of Eq. (5), the resummation of the large logarithmic contributions would affect not only the small-q T region, but also the region of large values of q T . This can easily be understood by observing that the logarithmic expansion parameter L diverges also when b → 0. To reduce the impact of unjustified higher-order contributions in the large-q T region, the logarithmic variable L in Eq. (5) is actually replaced byL ≡ ln (Q 2 b 2 /b 2 0 + 1) [21,22]. This replacement has an additional and relevant consequence: after inclusion of the finite component (see Eq. (8)), we exactly recover the fixed-order perturbative value of the total cross section upon integration of the q T distribution over q T (i.e., the resummed terms give a vanishing contribution upon integration over q T ).
We now turn to consider the finite component of the transverse-momentum cross section (see Eq. (3)). Since dσ (fin.) V does not contain large logarithmic terms in the small-q T region, it can be evaluated by truncation of the perturbative series at a given fixed order. In practice, the finite component is computed starting from the usual fixed-order perturbative truncation of the partonic cross section and subtracting the expansion of the resummed part at the same perturbative order.
Introducing the subscript f.o. to denote the perturbative truncation of the various terms, we have: This matching procedure between resummed and finite contributions guarantees to achieve uniform theoretical accuracy over the region from small to intermediate values of transverse momenta.
At large values of q T , the resummation (and matching) procedure is eventually superseded by the customary fixed-order calculations (their theoretical accuracy in the large-q T region cannot be improved by resummation of the logarithmic terms that dominate in the small-q T region).
In summary, the inclusion of the functions g (1) , g in the resummed component, together with the evaluation of the finite component at LO (i.e. at O(α S )), allows us to perform the resummation at NLL+LO accuracy. This is the theoretical accuracy used in our previous study [24] of the DY q T distribution. Including also the functions g  N (see, e.g., Eq. (24) in Ref. [21]), is not yet known. In the following, we assume that the value of A (3) is the same as the one [35,36] that appears in resummed calculations of soft-gluon contributions near partonic threshold. Using the recently computed H V (2) N coefficient [5], we are thus able to present the complete result for the DY q T -distribution up to NNLL+NLO accuracy. We point out that the NNLL+NLO (NLL+LO) result includes the full NNLO (NLO) perturbative contribution in the small-q T region. In particular, the NNLO (NLO) result for the total cross section is exactly recovered upon integration over q T of the differential cross section dσ V /dq T at NNLL+NLO (NLL+LO) accuracy.
We conclude this section with some comments on the numerical implementation of our calculation. Within our formalism, the resummation factor W V N (b, M) is directly defined, at fixed M, in the space of the conjugate variables b and N. To obtain the hadronic cross section, we have to perform inverse integral transformations: the Bessel transformation in Eq. (4) and an inverse Mellin transformation. These integrals are carried out numerically. The Mellin inversion requires the numerical evaluation of some basic N-moment functions that appear in the expression of the the second-order coefficients H V (2) qq←ab,N [5]: this evaluation has to be performed for complex values of N, and we use the numerical results of Ref. [37]. We recall [21] that the resummed form factor exp{G N (α S (µ 2 R ), L)} is singular at the values of b where α S (µ 2 R ) L ≥ π/β 0 (β 0 is the first-order coefficient of the QCD β function). Performing the Bessel transformation with respect to the impact parameter b (see Eq. (4)), we deal with this singularity as we did in Ref. [21], by using the regularization prescription of Refs. [38,39]: the singularity is avoided by deforming the integration contour in the complex b space.
3 Numerical results for Z/γ * production at the Tevatron In this section we consider Z/γ * production in pp collisions at Tevatron energies. We present our resummed results at NNLL+NLO accuracy, we compare them with the NLL+LO results (the NLL+LO results in Ref. [24] were obtained by using the MRST2004 NLO parton densities [40]), and we comment on the comparison with Tevatron Run II data [30,31].
The hadronic q T cross section at NNLL+NLO (NLL+LO) accuracy is computed by using the MSTW2008 NNLO (NLO) parton densities [41], with α S (µ 2 R ) evaluated at 3-loop (2-loop) order. This choice of the order of the parton densities and α S is fully justified both in the small-q T region (where the calculation of the partonic cross section includes the complete NNLO (NLO) result and is controlled by NNLL (NLL) resummation) and in the intermediate-q T region (where the calculation is constrained by the value of the NNLO (NLO) total cross section).
As for the electroweak couplings, we use the so called G µ scheme, where the input parameters are G F , m Z , m W . In particular, we use the PDG 2008 [42] values G F = 1.16637 × 10 −5 GeV −2 , m Z = 91.1876 GeV, Γ Z = 2.4952 GeV, m W = 80.398 GeV. Our calculation implements the decays γ * → l + l − and Z * → l + l − at fixed value of the invariant mass of the l + l − pair. In particular, we include the effects of the γ * Z interference and of the finite width of the Z boson. Nonetheless, the numerical results presented below are obtained by simply using the narrow-width approximation and neglecting the photon contribution. We find that this approximation works to better than 1% accuracy in the inclusive region of lepton invariant mass that is covered by the D0 data § .
As discussed in Sect. 2, the resummed calculation depends on the factorization and renormalization scales and on the resummation scale Q. Our convention to compute factorization and renormalization scale uncertainties is to consider independent variations of µ F and µ R by a factor of two around the central values µ F = µ R = m Z (i.e. we consider the range m Z /2 ≤ {µ F , µ R } ≤ 2 m Z ), with the constraint 0.5 ≤ µ F /µ R ≤ 2. Similarly, we follow Ref. [24] and we choose Q = m Z /2 as central value of the resummation scale, considering scale variations in the range m Z /4 < Q < m Z .
In Fig. 1 (left panel) we present the NLL+LO q T spectrum at the Tevatron Run II ( (3)), multiplied by a factor of 10 to make it more visible, is also shown for comparison (dotted line). We see that the LO result diverges to +∞ as q T → 0. The resummation of the small-q T logarithms leads to a well-behaved distribution: it vanishes as q T → 0, has a § The measured q T spectra are inclusive over the following regions of lepton invariant mass: 70-110 GeV [30] and 65-115 GeV [31]. ¶ Analogous results at the Tevatron Run I ( √ s = 1.8 TeV), obtained by using the MRST2004 parton densities [40], were presented in the left panel of Fig. 6 of Ref. [24].

Figure 1: The q T spectrum of Z bosons at the Tevatron Run II: results at NLL+LO (left panel) and NNLL+NLO (right panel) accuracy. Each result is compared to the corresponding fixed-order result (dashed line) and to the finite component (dotted line) in Eq. (8).
kinematical peak at q T ∼ 2 GeV, and tends to the corresponding LO result at large values of q T . The finite component smoothly vanishes as q T → 0 and gives a small contribution to the NLL+LO result in the low-q T region.
The results in the right panel of Fig. 1 are analogous to those in the left panel, although systematically at one order higher. The q T spectrum at NNLL+NLO accuracy (solid line) is compared with the NLO result (dashed line) and with the NLO finite component of the spectrum (dotted line). The NLO result diverges to −∞ as q T → 0 and, at small values of q T , it has an unphysical peak (the top of the peak is above the vertical scale of the plot) that is produced by the numerical compensation of negative leading and positive subleading logarithmic contributions. The contribution of the NLO finite component to the NNLL+NLO result is smaller than 1% at the peak and becomes more important as q T increases: it is about 8% at q T ∼ 20 GeV, about 20% at q T ∼ 30 GeV and about 53% at q T ∼ 50 GeV. A similar quantitative behaviour is observed by considering the contribution of the NLO finite component to the NLO result; the contribution is about 10% at q T ∼ 20 GeV, about 22% at q T ∼ 30 GeV and about 60% at q T ∼ 50 GeV. In the region of intermediate values of q T (say, around 20 GeV), the difference between the NNLL+NLO and NLO results is larger than the size of the NLO finite component. This difference is produced by the logarithmic terms (at NNLO and beyond NNLO) that are included in the resummed calculation at NNLL accuracy. At large values of q T the contribution of the NLO finite component sizeably increases. This behaviour indicates that the logarithmic terms are no longer dominant and that the resummed calculation cannot improve upon the predictivity of the fixed-order expansion.
Comparing the left and right panels of Fig. 1, we see that the q T spectrum is slightly harder at NNLL+NLO accuracy than at NLL+LO accuracy. The height of the peak at NNLL+NLO is lower than at NLL+LO. This is mainly due to the fact that the NNLO total cross section, which fixes the value of the q T integral of our NNLL+NLO result, is only about 3% larger than the NLO total cross section, whereas in the region of intermediate values of q T the cross section at NLO (and, correspondingly, at NNLL+NLO) is definitely larger than at LO (and, correspondingly, at NLL+LO); this leads to a reduction of the cross section at small q T . In Fig. 2 we show the scale dependence of the NLL+LO (dashed lines) and NNLL+NLO (solid lines) results. In the left panel we consider variations of the renormalization and factorization scales. The bands are obtained by varying µ R and µ F as previously described in this section. We note that, in the region of small and intermediate transverse momenta (q T ∼ < 30 GeV), the NNLL+NLO and NLL+LO bands overlap. This feature, which is not present in the case of the fixed-order perturbative results at LO and NLO (see Figs. 2 and 3 in Ref. [24]), confirms the importance of resummation to achieve a stable perturbative prediction. In the region of small and intermediate values of q T , the main difference between the NNLL+NLO and NLL+LO predictions is in the size of the scale variation bands. Going from NLL+LO to NNLL+NLO accuracy, we observe a reduction of the scale dependence from ±4% to ±3% at the peak, from ±7% to ±3% at q T ∼ 20 GeV, and from ±7% to ±5% at q T ∼ 50 GeV. We point out that the q T region where resummed perturbative predictions are definitely significant is a wide region from intermediate to relatively-small (say, close to the peak of the distribution) values of q T . In fact, at very small values of q T (e.g. q T ∼ < 5 GeV) the size of non-perturbative effects is expected to be important, while in the high-q T region (e.g. q T ∼ > 60 GeV) the resummation of the logarithmic terms cannot improve the predictivity of the fixed-order perturbative expansion. The inset plot in the left panel of Fig. 2 shows the region from intermediate to large values of q T . At large q T , the NLL+LO and NNLL+NLO results deviate from each other, and the deviation increases as q T increases. As previously stated, this behaviour is not particularly worrying since, in the large-q T region, the resummed results loose their predictivity and can (should) be replaced by customary fixed-order results.
In the right panel of Fig. 2 we consider resummation scale variations. The bands are obtained by fixing µ R = µ F = m Z and varying Q between m Z /4 and m Z . Performing variations of the resummation scale, we can get further insight on the size of yet uncalculated higher-order logarithmic contributions at small and intermediate values of q T . We find that the scale dependence at NNLL+NLO (NLL+LO) is about ±5% (±12%) in the region of the peak, and about ±5% (±16%) in the region where q T ∼ 20 GeV. We note that in a wide region of q T values, 5 GeV ∼ < q T ∼ < 50 GeV, the resummation scale dependence is reduced by, roughly, a factor of 2 in going from the NLL+LO to the NNLL+NLO result. Comparing the left and right panels of Fig. 2, we see that, at NNLL+NLO accuracy, the resummation scale dependence is larger than (though, comparable to) the µ F and µ R dependence.
The integral over q T of the resummed NNLL+NLO (NLL+LO) spectrum is in agreement (for any values of µ R , µ F and Q) with the value of the corresponding NNLO (NLO) total cross section to better than 1%, thus checking the numerical accuracy of our code. We also note that the large-q T region gives a little contribution to the total cross section (see some numerical results in Sect. 3.2 of Ref. [24]); therefore, the total cross section constraint mainly acts as a perturbative constraint on the resummed spectrum in the region from intermediate to small values of q T . Figure 3: The normalized q T spectrum of Z bosons at the Tevatron Run II. The NNLL+NLO result is compared with the D0 data of Refs. [30,31]. The bands are obtained as described in the text.
The D0 Collaboration has measured the normalized q T distribution, 1 σ dσ dq T , from data at the Tevatron Run II in the e + e − [30] and µ + µ − [31] channels. In the left panel of Fig. 3 we report the D0 data and our corresponding results at NNLL+NLO accuracy. The NNLL+NLO band represents our estimate of the perturbative uncertainty, and it is obtained by performing scale variations as follows. We independently vary µ F , µ R and Q in the ranges m Z /2 ≤ {µ F , µ R } ≤ 2m Z and m Z /4 ≤ Q ≤ m Z , with the constraints 0.5 ≤ µ F /µ R ≤ 2 and 0.5 ≤ Q/µ R ≤ 2. The constraint on the ratio µ F /µ R is the same as used in the left panel of Fig. 2; it has the purpose of avoiding large logarithmic contributions (powers of ln(µ 2 F /µ 2 R )) that arise from the evolution of the parton densities. Analogously, the constraint on the ratio Q/µ R avoids large logarithmic contributions (powers of ln(Q 2 /µ 2 R )) in the perturbative expansion of the resummed form factor exp{G N } (see Eq. (6)). We recall (see e.g. Eq. (19) of Ref. [21]) that the exponent G N of the form factor is obtained by q 2 integration of perturbative functions of α S (q 2 ) over the range b 2 0 /b 2 ≤ q 2 ≤ Q 2 . To perform the integration with systematic logarithmic accuracy, the running coupling α S (q 2 ) is then expressed in terms of α S (µ R ) (and ln(q 2 /µ 2 R )). As a consequence, the renormalization scale We do not apply additional constraints on the ratio Q/µ F , since the form factor does not depend on µ F . µ R should not be too different from the resummation scale Q, which controls the upper bound of the q 2 integration.
The D0 data and the NNLL+NLO band are presented in the left panel of Fig. 3. The inset plot shows the region from q T = 30 GeV up to q T = 100 GeV. A quick inspection of the figure shows that the data are described quite well by the NNLL+NLO perturbative predictions.
Differences and similarities between theoretical calculations and the data are more clearly visible by considering their fractional difference with respect to a 'reference' theoretical result. We choose the NNLL+NLO result at central values of the scales (i.e. µ F = µ R = m Z , Q = m Z /2) as 'reference' theory, and we show the ratio (X−theory)/theory in the right panel of Fig. 3. The label X refers to either the experimental data or the NNLL+NLO (solid lines) and NLO (dashed lines and dotted line) results, including their scale dependence.
Considering the right panel of Fig. 3, we first comment on the scale uncertainty band of the NNLL+NLO result (solid lines). Such uncertainty is about ±6% at the peak, it decreases to about ±4-5% in the region up to q T = 10 GeV, and then it increases, reaching the size of about ±12% at q T = 50 GeV. In the region beyond q T ∼ 60 GeV the resummed result looses predictivity, and its perturbative uncertainty becomes large. The right panel of Fig. 3 also shows the scale variation band of the NLO result. The NLO band (dashed lines) is obtained by varying µ F and µ R (the NLO calculation does not depend on the resummation scale Q) as in the NNLL+NLO calculation. We comment on the comparison between the NLO and NNLL+NLO bands. At large values * * of q T , the NLO and NNLL+NLO bands overlap (the NLO and NNLL+NLO are certainly consistent), and the NLO result has a smaller uncertainty. At intermediate values of transverse momenta, the NLO result is lower than the NNLL+NLO result, and the corresponding scale variation bands do not overlap. We recall (see the discussion in Sect. 3.1 of Ref. [24]) that in this region the NLO band underestimates the true perturbative uncertainty of the NLO result; indeed, the NLO band and the corresponding LO band do not overlap at intermediate values of q T (see Figs. 3 and 5 in Ref. [24]). To get some quantitative insight into the 'true' perturbative uncertainty of the NLO calculation in this q T region, we can consider wider scale variations and, in particular, we can lower the values of µ F and µ R . In the right panel of Fig. 3, we show the NLO band that we obtain by independently varying µ F and µ R in the range m Z /4 ≤ {µ F , µ R } ≤ 2m Z , with the constraint 0.5 ≤ µ F /µ R ≤ 2: this band is delimited by the dotted line and the lower dashed line (the region between the dotted line and the central values of the dashed band roughly corresponds to scale variations in the range m Z /4 ≤ {µ F , µ R } ≤ m Z ). We note that lowering the scales at NLO improves the consistency between the NLO and NNLL+NLO results. We also note that we have considered similar enlarged scale variations at NNLL+NLO accuracy, and we have checked that they do not significantly modify the NNLL+NLO band in the intermediate-q T region. This confirms the better stability of the NNLL+NLO calculation with respect to scale variations. In the small-q T region, the NLO result is theoretically unreliable. The NLO band quickly deviates from the NNLL+NLO band as q T decreases.
The right panel of Fig. 3 shows that the NNLL+NLO result is consistent with the D0 data, even at very low values of q T . We note that the resummed result is obtained in a perturbative framework. At low values of q T , non-perturbative effects are important and are expected (see, e.g., the final part of Sect. 5 in Ref. [24]) to shift the resummed result such as to improve the agreement with the data. In the region where q T ∼ < 50 GeV, the experimental errors and the corresponding NNLL+NLO errors overlap, with the sole exception of a couple of data points at very low q T . In the same region, the perturbative uncertainty of the NNLL+NLO result turns out to be comparable with the size of the experimental errors. As pointed out by the D0 Collaboration [30,31], the NLO result tends to undershoot the data in the region of intermediate values of q T : NNLL resummation improves the agreement with the data in this q T region.

Summary
In this paper we have considered the q T spectrum of DY lepton pairs produced in hadron collisions, and we have presented a perturbative QCD study based on transverse-momentum resummation at the NNLL order.
We have followed the formalism developed in Refs. [20,21], which is valid for the production of a generic high-mass system of non strongly-interacting particles in hadron collisions. The formalism combines small-q T resummation at a given logarithmic accuracy with the fixed-order calculations. It implements a unitarity constraint that guarantees that the integral over q T of the differential cross section coincides with the total cross section at the corresponding fixed-order accuracy. This leads to QCD predictions with a controllable and uniform perturbative accuracy over the region from small up to large values of q T . At large values of q T , the resummation formalism is superseded by customary fixed-order calculations.
We have considered the explicit case of DY lepton pairs from the decay of a Z boson produced at the Tevatron Run II. Using the recently computed NNLL coefficient H V (2) N [5], we have extended the NLL+LO resummed calculation presented in Ref. [24] to the NNLL+NLO accuracy. The NNLL corrections are not large and make the q T spectrum slightly harder. We have performed a study of the scale dependence of the calculation to estimate the corresponding perturbative uncertainty. In a wide region of transverse momenta (5 GeV ∼ < q T ∼ < 50 GeV) the size of the scale uncertainties is considerably reduced in going from NLL+LO to NNLL+NLO accuracy.
We have compared the resummed calculation with the results of measurements [30,31] of the normalized q T spectrum at the Tevatron Run II. The perturbative uncertainty of the NNLL+NLO results turns out to be comparable with the experimental errors. The NNLL+NLO results (without the inclusion of any non-perturbative effects) are consistent with the experimental data in a wide region of transverse momenta. Comparing the NNLL+NLO and NLO results, we have also shown that NNLL resummation improves the agreement with the data at intermediate values of q T . As is well known (and theoretically expected), the NLO result fails to describe the data at small values of q T .
More detailed comparisons with available data on vector boson production and further studies of theoretical uncertainties, including the impact of non-perturbative effects, are left to future investigations.
Note added. After the completion of this paper, the value of the coefficient A (3) (see the related comment in the final part of Sect. 2) for q T resummation was derived in Ref. [43]. We have checked the quantitative effect of this value of A (3) on our results for the q T distribution of Z bosons. We find that the effect is generally very small. The largest effect is produced in the region of very low values of q T ; for instance, in the case of Z production at the Tevatron, the quantitative effect is at the level of about 2% (4%) at q T ≃ 2 GeV (1 GeV).