Vector boson production at hadron colliders: transverse-momentum resummation and leptonic decay

We consider the transverse-momentum ($q_T$) distribution of Drell-Yan lepton pairs produced, via $W$ and $Z/\gamma^*$ decay, in hadronic collisions. At small values of $q_T$, we resum the logarithmically-enhanced perturbative QCD contributions up to next-to-next-to-leading logarithmic accuracy. Resummed results are consistently combined with the known ${\mathcal O}(\alpha_S^2)$ fixed-order results at intermediate and large values of $q_T$. Our calculation includes the leptonic decay of the vector boson with the corresponding spin correlations, the finite-width effects and the full dependence on the final-state lepton(s) kinematics. The computation is encoded in the numerical program DYRes, which allows the user to apply arbitrary kinematical cuts on the final-state leptons and to compute the corresponding distributions in the form of bin histograms. We present a comparison of our results with some of the available LHC data. The inclusion of the leptonic decay in the resummed calculation requires a theoretical discussion on the $q_T$ recoil due to the transverse momentum of the produced vector boson. We present a $q_T$ recoil procedure that is directly applicable to $q_T$ resummed calculations for generic production processes of high-mass systems in hadron collisions.


Introduction
The production of high-mass lepton pairs through the Drell-Yan (DY) mechanism [1] is a benchmark hard-scattering process at hadron colliders. It provides important tests of the Standard Model (SM) with both precise measurements of its fundamental parameters and, at the same time, stringent constraints on new physics.
It is thus a major task to achieve accurate theoretical predictions for the DY production cross section and related kinematical distributions. This requires, in particular, the evaluation of QCD radiative corrections, which can be perturbatively computed as power series expansion in the strong coupling α S . The total cross section [2] and the rapidity distribution [3] of the vector boson are known up to the next-to-next-to-leading order (NNLO) in perturbative QCD. Two independent fully exclusive NNLO calculations, which include the leptonic decay of the vector boson, have been performed [4,5,6]. Electroweak (EW) radiative corrections are also available for both W [7] and Z/γ * [8] production. Mixed QCD-EW corrections have been considered in Refs. [9,10,11,12].
A particularly relevant observable is the transverse-momentum (q T ) distribution of the vector boson. To obtain a precise measurement of the W mass it is important to have accurate theoretical calculations of the W and Z bosons q T spectra. 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 , QCD corrections are known up to O(α 2 S ) [13,14,15], and these results were extended in Refs. [16,17] with the inclusion of the dependence on the leptonic decay variables. Very recently the fully exclusive O(α 3 S ) computation of vector boson production in association with a jet has been performed in Ref. [18] (in the case of W production) and Ref. [19] (in the case of Z/γ * production).
The bulk of the vector boson cross section is produced in the small-q T region (q T ≪ m V ), where the reliability of the fixed-order expansion is spoiled by the presence of large logarithmic corrections, α n S (m 2 V /q 2 T ) ln m (m 2 V /q 2 T ) (with 0 ≤ m ≤ 2n−1), of soft and collinear origin. To obtain reliable predictions, these logarithmically-enhanced terms have to be evaluated and systematically resummed to all orders in perturbation theory [20]- [29]. In recent years, the resummation of smallq T logarithms has been reformulated [30]- [38] by using Soft Collinear Effective Theory (SCET) methods and transverse-momentum dependent (TMD) factorization.
The resummed and fixed-order calculations, which are valid at small and large values of q T , respectively, can be consistently matched at intermediate values of q T to achieve a uniform theoretical accuracy for the entire range of transverse momenta.
In this paper we compute the vector boson transverse-momentum distribution [39,40] by using the resummation formalism proposed in Refs. [25,26,27], which can be applied to a generic process in which a high-mass system of non strongly-interacting particles is produced in hadronic collisions [41,26,27], [42]- [52]. Other phenomenological studies of the DY q T distribution, which combine resummed and fixed-order perturbative results at different levels of theoretical accuracy, can be found in Refs. [53]- [69]. Within the studies in Refs. [53]- [69], the kinematical dependence on the momenta of the final-state leptons is considered only in the RESBOS calculation [55,56,68] and in the calculations of Refs. [57] and [67].
Hadron collider experiments can directly measure only the decay products of vector bosons in finite kinematical regions. Therefore, it is important to include the vector boson leptonic decay in the theoretical calculations, by retaining the kinematics of the final-state leptons. In this way it is possible to obtain predictions for the transverse-momentum distribution of the measured leptons. This is specially relevant in the case of W production where, because of the final-state neutrino, the transverse momentum of the vector boson can only be reconstructed through a measure of the hadronic recoil. Moreover, in both cases of W and Z production, the inclusion of the leptonic decay allows one to apply kinematical selection cuts, thus providing a more realistic simulation of the actual experimental analysis.
In Ref. [39,40] we have presented a resummed computation of the transverse-momentum spectrum for Z/γ * production at Tevatron energies. We have combined resummation at the nextto-next-to-leading logarithmic (NNLL) accuracy in the small-q T region with the fixed-order results at O(α 2 S ) in the large-q T region. This leads to a calculation with uniform theoretical accuracy from small to intermediate values of q T . In particular, the integral over the range 0 ≤ q T ≤ q T max (q T max is a generic upper limit in the small-q T region) of the q T distribution includes the complete perturbative terms up to NNLO. Moreover, at large values of q T the calculation implements a unitarity constraint that guarantees to exactly reproduce the NNLO value of the total cross section after integration over q T . In this paper we extend the NNLL+NNLO calculation of Ref. [40] to W boson production, and we include the leptonic decay of the vector boson with the corresponding spin correlations. The spin of the vector boson dynamically correlates the decaying lepton momenta with the transverse momentum acquired by the vector boson through its production mechanism. Therefore, the inclusion of the full dependence on the lepton decay variables in the resummed calculation requires a theoretical discussion on the treatment of the q T recoil due to the transverse momentum of the vector boson. We treat the q T recoil by introducing a general procedure that is directly applicable to q T resummed calculations for generic production processes of high-mass systems in hadron collisions. The calculation presented in this paper parallels the one performed in Ref. [43] for the case of SM Higgs boson production, with the non-trivial additional complication of dealing with the spin correlations that are absent in the Higgs boson case. Our vector boson computation is implemented in the numerical code DYRes, which allows the user to apply arbitrary kinematical cuts on the final-state leptons and to compute the corresponding relevant distributions in form of bin histograms. The code DYRes is publicly available and it can be downloaded from the URL address http://pcteserver.mi.infn.it/~ferrera/dyres.html .
The paper is organized as follows. In Sect. 2 we briefly review the resummation formalism of Refs. [25,26,27], and we discuss the main features of q T resummation for the DY process with full dependence on final-state lepton variables. In Sect. 3 we present our quantitative results for vector boson production at LHC energies. Section 3.1 is devoted to the q T spectrum of the vector boson after integration over the final-state leptons. We present results at different orders of logarithmic accuracy, we study the corresponding dependence on scale variations, and we briefly comment on uncertainties due to parton densities and on non-perturbative effects. In Sect. 3.2 we compare our numerical results for Z/γ * and W production with some of the available LHC data, and we also study the impact of transverse-momentum resummation on lepton kinematical variables. In Sect. 4 we summarize our results. The Appendix presents a detailed discussion of q T recoil and of its implementation.

Transverse-momentum resummation
In this Section we briefly recall the main features of the transverse-momentum resummation formalism that we use in this paper. A more detailed discussion of the resummation formalism can be found in Refs. [25,26,27,29]. In Ref. [40] we have considered NNLL resummation for the q T distribution of the vector boson after integration over the kinematical variables of the decaying leptons and the rapidity of the vector boson. In this paper we extend the results of Ref. [40] to include the entire kinematical dependence on the final-state leptons. The presentation in this Section parallels that of Sect. 2 in Ref. [40] and, in particular, we highlight the main differences that arise in the treatment of the rapidity of the vector boson and, especially, of the lepton kinematics.
We consider the inclusive hard-scattering process where the collision of the two hadrons h 1 and h 2 with momenta P µ 1 and P µ 2 produces the vector boson V (V = W + , W − , Z and/or γ * ) with total momentum q µ , which subsequently decays in the lepton pair l 3 l 4 , and X denotes the accompanying final-state radiation. We consider high values of the invariant mass M of the lepton pair (in general, M differs from the on-shell mass m V of the vector boson V ), and we treat the colliding hadrons and the leptons in the massless approximation (P 2 1 = P 2 2 = p 2 3 = p 2 4 = 0) throughout the paper. In a reference frame where the colliding hadrons are back-to-back, the momentum q µ is fully specified by the invariant mass M (M 2 = q 2 ), the two-dimensional transverse-momentum vector q T (with magnitude q T = q 2 T and azimuthal angle φ q T ) and the rapidity y (y = 1 2 ln q·P 2 q·P 1 ) of the vector boson. Analogously the momentum p µ j of the lepton l j (j = 3, 4) is specified by the lepton rapidity y j and transverse momentum p T j .
The kinematics of the lepton pair is completely specified by six independent variables (e.g., the three-momenta of the two leptons). For our purposes, it is convenient to use the vector boson momentum q µ to select four independent variables. Therefore, the final-state lepton kinematics is fully determined by the vector boson momentum q µ = p µ 3 + p µ 4 and by two additional and independent variables that specify the angular distribution of the leptons with respect to the vector boson momentum q µ . We generically denote these two additional kinematical variables as Ω = {Ω A , Ω B }. These two independent variables can be chosen in different ways. For instance, we can use longitudinally boost invariant variables such as the rapidity difference y 3 − y and the azimuthal angle φ 3 (or the azimuthal angle difference φ 3 − φ q T ) of the lepton l 3 and the vector boson in the hadronic back-to-back reference frame. Alternatively, we can use the polar and azimuthal angles {θ ′ , φ ′ } of one lepton in a properly specified rest frame of the vector boson (such as, for instance, the Collins-Soper rest frame [70]). Independently of the actual specification of the variables Ω, the most general fully-differential hadronic cross section is expressed in terms of the sixfold differential distribution where s = (P 1 + P 2 ) 2 = 2P 1 · P 2 is the square of the hadronic centre-of-mass energy. Obviously, the differential distribution also depends on the EW parameters (including the mass m V of the vector boson V ): unless otherwise specified, this dependence is not explicitly denoted throughout the paper.
The differential hadronic cross section can be written as where f a/h (x, µ 2 F ) (a = q f ,q f , g) are the parton densities of the colliding hadron h at the factorization scale µ F , dσ a 1 a 2 →l 3 l 4 are the differential partonic cross sections,ŝ = x 1 x 2 s is the square of the partonic centre-of-mass energy,ŷ = y − ln x 1 /x 2 is the vector boson rapidity with respect to the colliding partons, and µ R is the renormalization scale. Note that the partonic cross sections do not have any explicit dependence on hadronic kinematical variables, since the leptonic variables Ω are specified with respect to q µ . The partonic cross section dσ a 1 a 2 →l 3 l 4 is computable in QCD perturbation theory as a power series expansion in the QCD coupling α S .
In the region where q T ∼ M, the perturbative expansion of the partonic cross section starts at O(α S ). In this region the value of the auxiliary scales µ F and µ R can be chosen to be of the order of M, and the QCD perturbative series is controlled by a small expansion parameter α S (M 2 ). Therefore, fixed-order calculations of the partonic cross section are theoretically justified. The QCD radiative corrections are known analytically up to O(α 2 S ) after integration over the lepton angular variables [14,15] and with the inclusion of the full dependence on these angular variables [16,17]. The numerical results at O(α 2 S ) can be obtained also from the fully-exclusive calculations of Refs. [4,5,6]. Results at O(α 3 S ) can be derived from the recent numerical computations of W + jet production [18] and Z/γ * + jet production [19].
In the small q T region (q T ≪ M), the perturbative computation of the partonic cross section starts at O(α 0 S ) through the leading-order (LO) EW process q fqf ′ → V of quark-antiquark annihilation. In this region, the QCD radiative corrections are known up to NNLO (i.e., O(α 2 S )) in analytic form [74] by neglecting corrections of O(q T /M) (these corrections can directly be extracted from Refs. [14,15,16,17]). The complete (i.e., by including corrections of O(q T /M)) NNLO result can be obtained from the numerical computations of Refs. [4,5]. However, in the small q T region the convergence of the fixed-order perturbative expansion is spoiled by the presence of powers of large logarithmic terms, α n S (M 2 /q 2 T ) ln m (M 2 /q 2 T ) (with 0 ≤ m ≤ 2n − 1). In particular, these terms become singular in the limit q T → 0. To obtain reliable predictions these terms have to be resummed to all orders.
Within our formalism, the resummation is performed at the level of the partonic cross section, which is decomposed as follows: Here we have introduced a shorthand notation: the symbol dσ a 1 a 2 →l 3 l 4 denotes the multidifferential partonic cross section that appears as the last factor in the right-hand side of Eq. (3). The first term, dσ (res.) , on the right-hand side of Eq. (4) is the resummed component. It contains all the logarithmically-enhanced contributions (at small q T ) that have to be resummed to all orders in α S . The second term, the finite component dσ (fin.) , is free of such contributions and thus it can be evaluated at fixed order in perturbation theory. Note that part of the non-singular (i.e., not logarithmically-enhanced) contributions can also be included in dσ (res.) , and we comment later about this point.
The resummation of the logarithmic contributions has to be carried out in the impact parameter (b) space [20,21,22,53,23] to fulfil the important constraint of transverse-momentum conservation for inclusive multiparton radiation. The impact parameter b is the conjugate variable to q T through a Fourier transformation. The small-q T region (q T ≪ M) corresponds to the large-b region (bM ≫ 1) and the logarithmic terms ln(M 2 /q 2 T ) become large logarithmic contributions ln(M 2 b 2 ) in b space. The resummed component of the cross section is then obtained by performing the inverse Fourier transformation (or the Bessel transformation in Eq. (6)) from b space to q T space. The resummed component of the partonic cross section in Eq. (4) can be expressed as wherê and J 0 (x) is the 0th-order Bessel function. The factor dσ (0) b 1 b 2 →l 3 l 4 in the right-hand side of Eq. (5) is the Born level differential cross section for the partonic subprocess q fqf ′ → V → l 3 l 4 of quarkantiquark annihilation, where the quark flavours f and f ′ can be either different (if V = W ± ) or equal (if V = Z, γ * ). This factor is of purely EW origin, and it completely encodes the dependence on the lepton kinematical variables Ω. We postpone more detailed comments on dσ (0) (see Eq. (12) and the discussion therein). The QCD radiative corrections and their associated dependence on ln(M 2 b 2 ) are embodied in the resummed factor W a 1 a 2 ,b 1 b 2 →V , which depends on the produced vector boson V but it is independent of the decay leptons (in particular, it does not depend on Ω). The integrand W in Eq. (6) depends on b 2 = b 2 and the inverse Fourier transformation is recast in terms of the Bessel transformation through the integration over the azimuthal angle of b. Note that the resummation factorŴ a 1 a 2 ,b 1 b 2 →V depends on q 2 T and it does not contain any dependence on the azimuthal angle φ q T of q T . This azimuthal independence is a feature of transverse-momentum resummation [23] for the production processes of colourless systems (such as vector bosons) through quark-antiquark annihilation. In contrast, logarithmically-enhanced azimuthal correlations enter transverse-momentum resummation for processes initiated by gluongluon fusion [28] (such as Higgs boson production) and for production of systems that carry colour charges (such as heavy quarks) [71] through either quark-antiquark annihilation or gluon-gluon fusion.
The all-order resummation structure of W a 1 a 2 ,b 1 b 2 →V in Eq. (6) can be organized in exponential form [26,27]. The exponentiated structure is directly evident by considering the 'double' (N 1 , N 2 ) Mellin moments W † For the sake of simplicity, in this presentation we omit the explicit dependence on the parton indices {a 1 a 2 , b 1 b 2 }. This simplified notation applies to the case of a sole parton species or, more precisely, to flavour non-singlet partonic channels (see Refs. [26,27] for the general case).
where the dependence on b (and on the large logarithm ln(M 2 b 2 )) is denoted by defining and introducing the logarithmic expansion parameter L ≡ ln(Q 2 b 2 /b 2 0 + 1) with b 0 = 2e −γ E (γ E = 0.5772... is the Euler number). The scale Q ∼ M, named resummation scale [41], which appears in the right-hand side of Eq. (7), parametrizes the arbitrariness in the resummation procedure. Although W does not depend on Q when evaluated to all perturbative orders, its explicit dependence on Q occurs when it is computed by truncation of the resummed expression at some level of logarithmic accuracy (see Eq. (8)). Variations of Q around M can thus be used to estimate the size of yet uncalculated higher-order logarithmic contributions.
The contribution exp{G (N 1 ,N 2 ) } in the right-hand side of Eq. (7) includes the Sudakov form factor and collinear-evolution terms. This contribution (which does not depend on the factorization scale µ F ) is universal (i.e. process independent), namely, it is independent on the produced vector boson V and, more generally, it occurs in transverse-momentum resummation for all the processes that are initiated by quark-antiquark annihilation at the LO level. The generalized form factor exp{G (N 1 ,N 2 ) } contains all the terms that order-by-order in α S are logarithmically divergent as b → ∞ (or, equivalently, as q T → 0). The all-order expression of the form factor can be systematically expanded in terms of functions g (k) (α S L) of the resummation parameter α S (µ 2 R ) L (each function g (k) (α S L) resums terms α n S L n and it is defined such that g (k) (0) = 0). The resummed logarithmic expansion of G (N 1 ,N 2 ) in powers of α S (µ 2 R ) reads where the term L g (1) collects the leading logarithmic (LL) contributions, the function g (2) includes the next-to-leading logarithmic (NLL) contributions [24], g (3) controls the NNLL terms [72,32] and so forth. The function H depends on the specific process of vector boson production and it is due to hard-virtual and collinear contributions. This function does not depend on the impact parameter b (it includes all the perturbative contributions to W (N 1 ,N 2 ) V that behave as constants in the limit b → ∞) and, therefore, it can be expanded in powers of α S = α S (µ 2 R ) as The next-to-leading order (NLO) term H is known since a long time [73], and the NNLO term H has been obtained more recently by two independent calculations in Refs. [74] and [75]. The explicit form of the functions G (N 1 ,N 2 ) and H (1) (N 1 ,N 2 ) V and, in particular, their dependence on the Mellin moment indices (N 1 , N 2 ) can be found in Ref. [26] and in Appendix A of Ref. [27].
Incidentally, we recall that the generalized form factor exp{G} is known up to NNLL accuracy also for processes initiated by the gluon fusion mechanism [76,77,28,32], and that the O(α 2 S ) collinear coefficients (which contribute to the NNLO term in Eq. (9)) are also known for all possible partonic channels [78,74,75,79,29]. Owing to the universality structure of transverse-momentum resummation, these results and those for the qq annihilation channel (which contribute to vector boson production) can be directly implemented in resummed calculations for production processes of generic high-mass systems.
The finite component dσ (fin.) in Eq. (4) has to be evaluated starting from the usual fixedorder perturbative truncation of the partonic cross section and subtracting the expansion of the resummed part at the same perturbative order. We have where the subscript f.o. denotes the perturbative truncation at the order f.o. (NLO, NNLO and so forth). The customary fixed-order component [dσ a 1 a 2 ] f.o. (and consequently also the finite component) definitely contains azimuthal correlations with respect to q T , although these are not logarithmically-enhanced in the small-q T region.
To obtain NLL+NLO accuracy we have to include the functions g (1) and g (2)(N 1 ,N 2 ) in the generalized form factor G (N 1 ,N 2 ) of Eq. (8), the function H S ) ‡ . This matching procedure between resummed and finite contributions guarantees to achieve uniform theoretical accuracy over the entire range of transverse momenta. In particular, we remark that the inclusion of H in the resummed component at the NNLL+NNLO level is essential to achieve NNLO accuracy in the small-q T region (considering a generic upper limit value q T max , the integral over the range 0 ≤ q T ≤ q T max of the q T distribution at the NNLL+NNLO level includes the complete perturbative terms up to NNLO). An analogous remark applies to the inclusion of H We have so far illustrated the resummation formalism for the most general sixfold differential partonic cross section dσ a 1 a 2 →l 3 l 4 (and for the corresponding hadronic cross section in Eq. (3)).
Starting from dσ a 1 a 2 →l 3 l 4 and performing integrations over some kinematical variables, we can obtain resummed results for more inclusive q T -dependent distributions. For instance, integrating over the lepton kinematical variables Ω, we obtain the q T cross section dσ/(d 2 q T dM 2 dy) at fixed invariant mass and rapidity of the lepton pair. The corresponding resummed component of the partonic cross section, as obtained from Eq. (5), is is the Born level (EW) total cross section for the partonic subprocess q fqf ′ → V → l 3 l 4 . By performing an additional integration over the rapidity y of the vector boson (lepton pair), we obtain dσ/(d 2 q T dM 2 ) and the corresponding resummed component of the partonic cross section simply involves the integration overŷ of the resummed factorŴ (q 2 T , M,ŷ,ŝ) in Eqs. (5) and (11) (or, equivalently, the factor W(b, M,ŷ,ŝ) in Eq. (6)). After integration overŷ, the ensuing ‡ This classification of the resummed+matched expansion exactly coincides with that of Refs. [26,40]. We simply note that we are using labels that differ from those used in Refs. [26,40]. The various terms of the expansion are denoted here (analogously to Ref. [43]) with the labels NLL+NLO and NNLL+NNLO, whereas they were denoted in Refs. [26,40] with the corresponding labels NLL+LO and NNLL+NLO. The fixed-order labels NLO and NNLO used here directly refer to the perturbative accuracy in the small-q T region (which corresponds to the perturbative accuracy of the total cross section), whereas the labels LO and NLO used in Refs. [26,40] were directly referring to the perturbative accuracy in the large-q T region. resummed factor depends on M andŝ, and it can be conveniently expressed in exponentiated form [26] by considering 'single' N Mellin moments with respect to the variable z = M 2 /ŝ at fixed M. The resummed expression for these 'single' N moments is exactly obtained by simply setting N 1 = N 2 = N in Eqs. (7)- (9). Our resummed calculation of dσ/(d 2 q T dM 2 ) was discussed in Ref. [39,40], and it is implemented in the numerical code DYqT. In Refs. [39,40] we presented detailed quantitative results for vector boson production at Tevatron energies. Results from DYqT at LHC energies are presented in the following Sect. 3.1.
Within our formalism the resummation of the large terms ln(M 2 /q 2 T ) at small values of q T is achieved by first performing the Fourier transformation of the q T cross section (or, more precisely, of its singular behaviour in the small-q T region) from q T space to b space (incidentally, the renormalization scale µ R and the others auxiliary scales Q and µ F are kept fixed and, especially, independent of q T in the integration over q T of the Fourier transformation). In b space, the large logarithmic variable (whose dependence has to be resummed) is L, at large values of b. Note that in the context of the resummation approach, the parameter α S (µ 2 R ) L is formally considered to be of order unity. Therefore, the ratio of two successive terms in the expansion (8) is formally of O(α S (µ 2 R )) (with no L enhancement). In this respect the resummed logarithmic expansion in Eq. (8) is as systematic as any customary fixed-order expansion in powers of α S (µ 2 R ). Analogously to any perturbative expansions, the perturbative terms g (k) (N 1 ,N 2 ) (α S (µ 2 R ) L; M/µ R , M/Q) in Eq. (8) have an explicit logarithmic dependence on ln(M/µ R ) or ln(M/Q) (see, e.g., Eqs. (22) and (23) in Ref. [26]). Therefore, to avoid additional large logarithmic enhancements that would spoil the formal behaviour of the expansion in Eq. (8), the renormalization scale µ R has to be set at a value of the order of M ∼ Q. A completely analogous reasoning applies to the µ F dependence of H and, therefore, we should set µ F ∼ M. In other words, once the enhanced perturbative dependence on b 2 M 2 (i.e., on the two different scales M and 1/b) is explicitly resummed (albeit at a definite logarithmic accuracy), we are effectively dealing with a single-scale observable at the hard scale M and we can set µ R ∼ µ F ∼ M in both the resummed and finite components of the q T cross section in Eq. (4).
We remark that setting µ = O(M) (here µ generically denotes the auxiliary scales µ R , µ F , Q) does not mean that the q T cross section is physically controlled by parton radiation with intensity that is proportional to α S (M 2 ). The resummed form factor exp{G (N 1 ,N 2 ) } in Eq. (7) (and the ensuing logarithmic expansion in Eq. (8)) is produced by multiparton radiation with intensity that is proportional to α S (k 2 ) and k 2 is a dynamical scale that varies in the range M 2 > k 2 > 1/b 2 (see, for instance, Eq. (19) in Ref. [26]), where 1/b 2 can be physically identified with q 2 T at small values of q T . Setting µ ∼ M in Eqs. (7) and (8) corresponds, roughly speaking, to consider the scale range We recall [26] a feature of our resummation formalism. The small-q T singular contributions that are resummed in Eqs. (5) (or Eq. (11)) are controlled by the large logarithmic parameter ln(M 2 /q 2 T ), which corresponds to L = ln( In our resummation formula (7), we actually use the logarithmic parameter L = ln(Q 2 b 2 /b 2 0 + 1) [41]. The motivations to use the logarithmic parameter L are detailed in Ref. [26] (see, in particular, the Appendix B and the comments that accompany Eqs. (16)- (18) and Eqs. (74)- (75) in Ref. [26]), and here we simply limit ourselves to recalling some aspects. In the relevant resummation region bQ ≫ 1, we have L = L + O(1/(Q 2 b 2 )) and, therefore, L and L are fully equivalent to arbitrary logarithmic accuracy (in other words, the replacement L ↔ L simply modifies the partition of small-q T non-singular contributions between the two components in the right-hand side of Eq. (4)). However, L and L have a very different behaviour as b → 0 (and, thus, they differently affect the q T cross section in the large-q T region § ). When bQ ≪ 1, we have L ≫ 1 and, therefore, the replacement L → L in Eq. (7) would produce the resummation of large and unjustified perturbative contributions in the large-q T region (strictly speaking, the replacement L → L leads to a q T cross section that is even not integrable over q T when q T → ∞: see, in particular, Eqs. (131) and (132) of the arXiv version of Ref. [26] and related accompanying comments). In contrast, when bQ ≪ 1 we have L → 0 and G (N 1 ,N 2 ) → 0. Therefore, the use of L reduces the impact of unjustified large contributions that can be introduced in the small-b region through the resummation procedure. Moreover, the behaviour of the form factor exp{G (N 1 ,N 2 ) } at b = 0 is related to the integral over q T of the q T -dependent cross section and, since we have exp{G (N 1 ,N 2 ) } = 1 at b = 0, our resummation formalism fulfils a perturbative unitarity constraint [26]: after inclusion of the finite component as in Eq. (10), the integration over q T of our resummed q T cross sections recovers the fixed-order predictions for the total cross sections. Specifically, the integral over q T of dσ/(d 2 q T dM 2 dy) and dσ/(d 2 q T dM 2 ) at the NNLL+NNLO (NLL+NLO) accuracy completely and exactly (i.e., with no additional higher-order contributions) agrees with the rapidity distribution dσ/(dM 2 dy) and the total cross section dσ/dM 2 at NNLO (NLO) accuracy, respectively. In summary, the expressions (7) and (8) . After having combined the resummed calculation at N k LL accuracy with the complete N k LO calculation, as in Eqs. (4) and (10), these parametrically-small corrections produce residual terms that start to contribute at the N k+1 LO level. Therefore, the use of L has the purpose of reducing the impact of unjustified and large higher-order (i.e., beyond the N k LO level) contributions that can be possibly introduced at intermediate and large values of q T through the resummation of the logarithmic perturbative behaviour at small values of q T . In particular, no residual higher-order contributions are introduced in the case of the total (integrated over q T ) cross section (which is the most basic quantity that is not affected by logarithmically-enhanced perturbative corrections).
We add some relevant comments about the dependence of the resummed cross section on the kinematical variables Ω that specify the angular distribution of the leptons with respect to the vector boson. By direct inspection of Eqs. (5) and (11) we see that they involve exactly the same resummation factorŴ . The only difference between the right-hand side of these equations arises form the Born level factors dσ (0) /dΩ andσ (0) , which are related as follows through the integration over Ω: dσ with the normalization condition Although both factors depend on EW parameters (EW couplings, mass and width of the vector boson), they have a different dependence on the relevant kinematical variables. The vector boson distribution dσ (res.) /(d 2 q T dM 2 dy) (and, analogously, dσ (res.) /(d 2 q T dM 2 )) involves the Born level total cross sectionσ (0) (M 2 ), which depends on M 2 , whereas the less inclusive leptonic distribution dσ (res.) /(d 2 q T dM 2 dy dΩ) involves the Born level differential cross section dσ (0) /dΩ that additionally depends on Ω and also on the transverse momentum q T of the lepton pair (see the function F in the right-hand side of Eq. (12)). To our knowledge the q T dependence of dσ (0) /dΩ has not received much attention in the previous literature on transverse-momentum resummation and, therefore, we discuss this issue with some details in Appendix A. Physically, this q T dependence is a necessary consequence of transverse-momentum conservation and it arises as a q T -recoil effect in transverse-momentum resummation. At the LO in perturbation theory the lepton angular distribution is determined by the Born level production and decay process of the vector boson, which carries a vanishing transverse momentum. Through the resummation procedure at fixed lepton momenta, higher-order contributions due to soft and collinear multiparton radiation dynamically produce a finite value of the transverse momentum q T of the lepton pair, and this finite value of q T has to be distributed between the two lepton momenta by affecting the lepton angular distribution. This q T -recoil effect on the Born level angular distribution is a non-singular contribution to the q T cross section at small values of q T and, therefore, it is not directly and unambiguously computable through transverse-momentum resummation. In other words, the Born level function F in Eq. (12) has the form where F (0; M 2 , Ω) is uniquely determined, whereas the small-q T corrections of O(q T /M) has to be properly specified. In any physical computations of lepton observables (i.e., in any computations that avoid possible unphysical behaviour due to violation of momentum conservation for the decay process q = p 3 (l 3 ) + p 4 (l 4 )) through transverse-momentum resummation, a consistent q Trecoil prescription has to be actually (either explicitly or implicitly) implemented ¶ . Note that, after having combined the resummed and finite components as in Eqs. (4) and (10), the O(q T /M) recoil effects lead to contributions that start at O(α 3 S ) (i.e., N 3 LO) in the case of our resummed calculation at NNLL+NNLO accuracy (correspondingly, these contributions start at O(α 2 S ) in the case of NLL+NLO accuracy). Obviously there are infinite ways of implementing the q T -recoil effect, and in Appendix A we explicitly describe a very general and consistent procedure . Note that the q T -recoil effect completely cancels after integration over the leptonic variables Ω (see Eq. (13)).
Our resummed calculation of the sixfold differential distribution in Eq. (2) is implemented in the numerical code DYRes, which allows the user to apply arbitrary kinematical cuts on the momenta of the final-state leptons and to compute the corresponding relevant distributions in form of bin histograms. We add some comments on the numerical implementation of our calculation. In Eqs. (7)-(9) we have illustrated the structure of the resummed component in the double (N 1 , N 2 ) ¶ The dynamical treatment of the q T recoil is embedded in the formulation of transverse-momentum (k T ) factorization [80,81] of hard-scattering processes at high energy (at small x).
The q T -recoil issue is not a specific issue of leptonic decay in vector boson production. The issue is completely general (see Appendix A), and it arises in any q T resummed calculation for the production of a set of particles with measured momenta at fixed total transverse momentum q T (e.g., diphoton, diboson, or heavy-quark pair production). A noticeable exception (as discussed in Appendix A) is the production of a SM Higgs boson and its subsequent decay. In this case, due to the spin-0 nature of the Higgs boson, the q T dependence of the corresponding Born level function F (q T /M ; M 2 , Ω) can be entirely determined by kinematics [43], without the necessity of specifying q T -recoil effects of dynamical origin.
Mellin space. Through the inverse Mellin transformation, this structure can equivalently be expressed in terms of convolutions with respect to longitudinal momentum fractions x 1 and x 2 (see Eq. (3)). In the DYRes code, the Mellin inversion is carried out numerically. The results for the NNLO term H (2) V in Eq. (9) are presented in Ref. [74] in analytic form directly in (x 1 , x 2 ) space. These results have to be transformed in Mellin space. Then, the Mellin inversion requires the numerical evaluation of some basic N-moment functions that appear in the expression of H this evaluation has to be performed for complex values of N, and we use the numerical results of Ref. [82]. This implementation of the resummed component is completely analogous to that of the DYqT code [39,40] and of other previous computations [27]. Nonetheless, the efficient generation of 'vector boson events' according to the multidifferential distribution of Eq. (4) and the inclusion of the leptonic decay are technically non trivial, and this requires substantial improvements in the computational speed of the numerical code that evaluates the resummed component of the cross section. The fixed-order (NLO and NNLO) cross section in Eq. (10) and then the finite component of the cross section in Eq. (4) are evaluated through an appropriate modification of the DYNNLO code [5]: DYNNLO is particularly suitable to this purpose, since it is based on the q T subtraction formalism [83], which uses the transverse-momentum resummation formalism to construct the subtraction counterterms.
Using the resummation expansion parameter L in Eq. (7) and the matching procedure (which implements the perturbative unitarity constraint on the total cross section) with the complete fixed-order calculation, our resummation formalism [26] formally achieves a uniform theoretical accuracy in the region of small and intermediate values of q T , and it avoids the introduction of large unjustified higher-order contributions in the large-q T region. In the large-q T region, the results of the resummed calculation are consistent with the customary fixed-order results and, typically [26,40], show larger theoretical uncertainties (e.g., larger dependence with respect to auxiliary-scale variations) with respect to the corresponding fixed-order results. This feature is not unexpected, since the theoretical knowledge (and the ensuing resummation) of large logarithmic contributions at small q T cannot improve the theoretical predictions at large values of q T . In the large-q T region, where the resummed calculation shows 'unjustified' large uncertainties and ensuing loss of predictivity with respect to the fixed-order calculation, the reliability of the resummed calculation is superseded by that of the fixed-order calculation. In this large-q T region, we can simply use the theoretical results of the fixed-order calculation. In the computation of quantities that directly and explicitly depend on q T (e.g., the transverse-momentum spectrum of the vector boson), it is relatively straightforward to identify and select 'a posteriori' the large-q T region where the resummed calculation is superseded by the fixed-order calculation. In the present work, however, we are also interested in studying kinematical distributions of the vector boson decay products: our goal is thus to generate the full kinematics of the vector boson and its (leptonic) decay, to apply the required acceptance cuts, and to compute the relevant distributions of the lepton kinematical variables. In this framework, the actual results can become sensitive to the large-q T region in which the resummed calculation cannot improve the accuracy of the fixed-order calculation. To reduce this sensitivity, in the DYRes implementation of the resummed calculation we thus introduce a smooth switching procedure at large value of q T by replacing the resummed cross section in Eq. (4) as follows: where the function w(q T ) is defined as T (16) and the function f (q T ) is chosen as We have quantitatively checked that the value of the parameter q sw.
T can be 'suitably' chosen in the large-q T region, and that both parameters q sw.
T and ∆q T can be consistently chosen so as not to spoil our unitarity constraint (in Sect. 3.1 we show that the integral over q T of our NLL+NLO and NNLL+NNLO resummed results still reproduces well the NLO and NNLO total cross sections). We note that we do not introduce any switching procedure in the DYqT calculation (though, its introduction is feasible) since, as previously mentioned, the identification of the large-q T region is straightforward in the computation of dσ/(d 2 q T dM 2 ).
We recall [26] that the resummed form factor exp{G(α S , L)} of Eq. (7) is singular at very large values of b. The singularity occurs in the region b ∼ > 1/Λ QCD , where Λ QCD is the momentum scale of the Landau pole of the perturbative running coupling α S (µ 2 ). This singularity is the 'perturbative' signal of the onset of non-perturbative (NP) phenomena at very large values of b (which practically affect the region of very small transverse momenta). In this region NP effects cannot any longer be regarded as small quantitative corrections and they have to be taken into account in QCD calculations. A simple and customary procedure to include NP effects is as follows. The singular behaviour of the perturbative form factor exp{G(α S , L)} is removed by using a regularization procedure * * and the resummed expression in Eq. (7) is then multiplied by a NP form factor and it is inserted as integrand of the b space integral in Eq. (6). The regularization procedure that was used in the DYqT calculation [40] is the 'minimal prescription' of Ref. [84,61], which basically amounts to avoid the singularity of exp{G(α S , L)} by deforming the integration contour of Eq. (6) in the complex b plane. In the DYRes calculation of the present work, we use a different regularization procedure by freezing the b dependence of exp{G(α S , L)} before reaching its singular point. The freezing procedure follows the 'b * prescription' of Refs. [22,23] and it is obtained by performing the replacement in the b dependence of G(α S , L). The value of the parameter b lim has to be large (b lim M ∼ b lim Q ≫ 1) but smaller than the value of b at which the singularity of exp{G(α S , L)} takes place (note that the replacement in Eq. (18) has a negligible effect at small and intermediate values . The use of the b * freezing procedure improves the (numerical) performances of the DYRes code. Additional comments on NP effects are presented in Sect. 3.1. * * We recall that the resummed form factor exp{G(α S , L)} produces a strong suppression (G(α S , L) ∝ −α S L 2 ) in the large-b region where α S L 2 ∼ > O(1). Therefore, the choice of different regularization procedures mildly affects [20,22,23] the results since its effects are relevant only in the region b ∼ O(1/Λ QCD ) where the b integral is strongly damped by the form factor.

Numerical results at the LHC
In this Section we consider the processes pp → Z/γ * → l + l − and pp → W ± → lν l at LHC energies. We present our resummed results at NNLL+NNLO and NLL+NLO accuracy, and we compare them with some of the available LHC data. We compute the hadronic cross sections at NNLL+NNLO (NLL+NLO) accuracy by using the NNPDF3.0 NNLO (NLO) [85] parton densities functions (PDFs), with α S (m 2 Z ) = 0.118 and with α S (µ 2 R ) evaluated at 3-loop (2-loop) order. As in customary fixed-order calculations at high invariant mass (M = O(m Z )), we consider N f = 5 flavours of light quarks and we treat them in the massless approximation.
As for the EW 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 2014 [86] values G F = 1.1663787 × 10 −5 GeV −2 , m Z = 91.1876 GeV, Γ Z = 2.4952 GeV, m W = 80.385 GeV and Γ W = 2.085 GeV and in the case of W ± production, we use the (unitarity constrained) CKM matrix elements V ud = 0.97427, V us = 0.22536, V ub = 0.00355, V cd = 0.22522, V cs = 0.97343, V cb = 0.0414. Our calculation implements the leptonic decays Z/γ * → l + l − and W → lν l (we include the effects of the Z/γ * interference and of the finite width of the W and Z bosons) with the corresponding spin correlations and the full dependence on the kinematical variables of final state leptons. This allows us to take into account the typical kinematical cuts on final state leptons that are considered in the experimental analysis. As discussed in Sect. 2, the resummed calculation at fixed lepton momenta requires a q T -recoil procedure. We implement a procedure that is described in Appendix A, and that is practically equivalent to compute the Born level distribution dσ (0) /dΩ of Eqs. (5) and (12) in the Collins-Soper rest frame [70] (this is exactly the same procedure as used in other resummed calculations [55,56,57,52]). As explained in Sect. 2, the DYRes resummed calculation uses a smooth switching procedure (see Eq. (15)) in the large-q T region. In our numerical implementation the parameters in Eq. (17) are chosen to be ∆q T = M/(2 √ 2) and q sw. T = 3M/4.

Inclusive results at fixed q T
We start the presentation of our results by discussing some general features of the q T spectrum that can be addressed at the inclusive level, i.e. after integration over the lepton angular variables Ω and over the rapidity y of the lepton pair. Unless otherwise specified, the numerical results of this Subsection are obtained by using the code DYqT [39,40]. The code DYqT is publicly available and it can be downloaded from the URL address http://pcteserver.mi.infn.it/~ferrera/dyqt.html .
We first consider the dependence on the auxiliary scales µ F , µ R and Q. These scales have to be set at values of the order of the invariant mass M of the produced system, with no definite preference for specific values. Then, scale variations around the chosen central value can be used to estimate the size of yet uncalculated higher-order terms and the ensuing perturbative uncertainties. In the NNLL+NNLO studies of Refs. [26,42,43] on Higgs boson production, and in our previous work on vector boson production [40] the central reference values of the scales were chosen as µ F = µ R = 2Q = m F , where m F is the mass of the produced boson (the Higgs boson mass in the case of Higgs boson production, and the vector boson mass in the case of vector boson production). In the case of Higgs boson production, this choice gives consistent NLL+NLO and NNLL+NNLO results with a reduced scale dependence at NNLL+NNLO level and, in particular, with a nice overlap of the NLL+NLO and NNLL+NNLO uncertainty bands (see, e.g., Fig. 2 of Ref. [42]). In the case of vector boson production, our previous studies were focused on Tevatron energies, and a similar pattern was observed [40]. When moving to vector boson production at LHC energies, we notice that the factorization-scale dependence exhibits a (slightly) different behaviour. , we see that the factorization-scale bands at NLL+NLO and NNLL+NNLO accuracy never overlap, except for the tiny region around q T ∼ 7 GeV where they cross each other. The lack of overlap is particularly evident in the peak region, where the bulk of the events is placed: here the central NLL+NLO and NNLL+NNLO results differ by about 30%. We also notice that throughout the region of small and intermediate values of q T (q T ∼ < 30 GeV) the size of the NLL+NLO band is rather small and it is always (with the exception of a small region around q T ∼ 8 GeV) smaller than the size of the NNLL+NNLO band, and this suggests that an accidental cancellation of the µ F dependence may occur at the NLL+NLO level with this choice of central scale. In the right panel we observe a µ F -dependence behaviour that is qualitatively similar but quantitatively different from that in the left panel. If µ F = µ R = Q = m Z /2 is the central scale choice (right panel), the NLL+NLO and NNLL+NNLO bands are closer and they overlap at small transverse momenta. The overlap occurs in a limited region of q T that, nevertheless, includes the peak region. The shape of the spectra appears closer when going from NLL+NLO to NNLL+NNLO accuracy, and the NLL+NLO band is wider than the NNLL+NNLO one in the small and intermediate region of q T . Note that the central values of µ R are µ R = m Z and µ R = m Z /2 in the left and right plot, respectively, but we have checked that this difference has little effect: the observed different behaviour is mainly due to the different central value of µ F . In summary, the µ F dependence observed in the left panel of Fig. 1 suggests that the corresponding scale variation bands (and especially the NLL+NLO band) are likely to underestimate the perturbative uncertainties of the calculation. Based on these observations, in the rest of the paper, we will adopt µ F = µ R = Q = m Z /2 as reference values of the central scales.
The differences between the NLL+NLO and NNLL+NNLO results have a twofold origin. Part of the differences is due to the next-order radiative corrections in the partonic cross sections, and the remaining part is due to the increased order of the PDFs. To quantify the impact of these two different contributions, we have considered the result that is obtained by convoluting the NLL+NLO partonic cross sections with the NNLO PDFs. This result, at the central scales that are considered in Figs. 1(a) and 1(b), is reported (see the black dotted line) in the corresponding lower panel, and we can see that it is quite close to the NLL+NLO result with NLO PDFs. In other words, a large part of the quantitative differences between the NLL+NLO and NNLL+NNLO results is due to the corresponding differences at the level of the partonic cross sections.  The NLL+NLO and NNLL+NNLO results for the q T spectrum of on-shell Z boson produced at the LHC with different collision energies are presented in Fig. 2. We consider two centre-of-mass energies: √ s =8 TeV (Fig. 2 left) and √ s =14 TeV (Fig. 2 right). At each logarithmic accuracy we present the result at the central value µ F = µ R = Q = m Z /2 of the scales and a corresponding band. The bands provide an estimate of the perturbative uncertainties of the calculations due to missing higher-order contributions. The bands are obtained through independent variations of µ F , µ R and Q by following the procedure of Ref. [40]: we independently vary µ F , µ R and Q in the range m Z /4 ≤ {µ F , µ R , Q} ≤ m Z with the constraints 0.5 ≤ µ F /µ R ≤ 2 and 0.5 ≤ Q/µ R ≤ 2.
We remind the reader that the constraint on µ F /µ R is introduced to avoid large logarithmic contributions (ln(µ F /µ R ) terms from the evolution of the parton densities) in the perturbative expansion of the hard/collinear factor H V of Eq. (7). Analogously, the constraint on Q/µ R avoids large logarithmic terms (ln(Q/µ R )) in the resummed expansion of the form factor exp{G} of Eq. (7). The lower panels in Fig. 2 present the ratio of the scale-dependent NLL+NLO and NNLL+NNLO results with respect to the NNLL+NNLO result at the central value µ F = µ R = Q = m Z /2 of the scales.
The region of small and intermediate values of q T is shown in the main panels of Fig. 2. At fixed centre-of-mass energy the NNLL+NNLO q T spectrum is harder than the spectrum at NLL+NLO accuracy. At fixed value of q T the cross section sizeably increases by increasing the centre-of-mass energy from 8 TeV to 14 TeV. The shape of the NNLL+NNLO q T spectrum is slightly harder at the higher energy. The NLL+NLO scale-variation band is wider than the NNLL+NNLO band. The NLL+NLO and NNLL+NNLO bands overlap at small transverse momenta and remain very close by increasing q T (the differences with respect to the plot on the right-hand side of Fig. 1 are due to the additional dependence on µ R and Q). The NNLL+NNLO (NLL+NLO) scale dependence is about ±10% (±20%) at the peak, it decreases to about ±2% (±7%) at q T ≃ 10 GeV and increases to about ±6% (±10%) at q T ∼ 25 GeV. Since the NNLL+NNLO and NLL+NLO bands do not exactly touch each other in the region where q T ∼ > 8 GeV, one may argue that the 'true' perturbative uncertainty of the NNLL+NNLO result in this region is slightly larger than the size of the NNLL+NNLO scale dependence band (for instance, one may use [39] the difference between the NNLL+NNLO central scale result and the upper line of the NLL+NLO band in Fig. 2 to estimate the uncertainty of the NNLL+NNLO result).
The inset plots show the cross section in the large-q T region. The resummation results obtained with DYqT and reported in the inset plots are presented for completeness and mainly for illustrative purposes. At large values of q T (q T ∼ > m Z ) the resummed result looses predictivity, and its perturbative uncertainty becomes large. In this region of transverse momenta we see that the uncertainty band increases in going from the NLL+NLO to the NNLL+NNLO level. However, as already mentioned in Sect. 2, at high q T the resummation cannot improve the predictivity of fixed-order calculations and the DYqT result in Fig. 2 cannot be regarded as reference theoretical result. The resummed result has to be replaced by the standard fixed-order prediction. The NNLO (NLO) result (which is not shown in Fig. 2) lies inside the NNLL+NNLO (NLL+NLO) band and the former has a smaller scale dependence than the latter. We also note that, at high q T , the preferred reference central scales µ R and µ F of the fixed-order prediction should be of the order of m 2 Z + q 2 T (rather than of the order of m Z ).
We also recall that, increasing q T throughout the high-q T region, fixed-order QCD calculations are affected by additional and potentially-large logarithmic terms. These are collinear (fragmentation) contributions [87,88], which become more relevant by increasing the ratio q T /M, and soft (threshold) contributions [89,90], which become more relevant by increasing the ratio q T / √ s (or q T / √ŝ ).
We have so far discussed only uncertainties from missing higher-order contributions. Before moving to consider the case in which cuts on the final-state leptons are applied, we briefly discuss two additional sources of QCD uncertainties on the q T spectrum: the uncertainty from PDFs and that from NP effects. We consider these effects in turn.
Modern sets of PDFs include an estimate of the errors (mainly experimental errors) in their determination from global data fits, and this estimate can then be used to compute the ensuing PDF uncertainty on the QCD calculation of hadron collider observables. In Fig. 3 we consider Z boson production at NNLL+NNLO accuracy. In Fig. 3 (a) we report the NNLL+NNLO results of Fig. 2 (a) ( √ s =14 TeV) and the effect of the PDF uncertainty at 68% CL on the NNLL+NNLO calculation at the central scale value µ F = µ R = Q = m Z /2. In Fig. 3  NP effects are known to increasingly affect the transverse-momentum spectrum as q T decreases towards q T → 0. A detailed study of these effects is beyond the scope of the present work. We limit ourselves to roughly estimate the possible impact of such effects, and we use a very simple model in which the perturbative form factor exp{G(α S , L) in Eq. (7) is multiplied by a NP form factor S N P (b) = exp{−g N P b 2 }, which produces a Gaussian smearing of the q T distribution at small-q T values. We vary the value of the parameter g N P in a quite wide ('conservative') range, 0 ≤ g N P ≤ 1.2 GeV 2 , and in Fig. 3(a) (black band) we show the ensuing quantitative effects on the q T spectrum. In Fig. 3(b) the NP effects are normalized with respect to the perturbative NNLL+NNLO result at central value of the scales.
Comparing the lower panels of Fig. 2 with Fig. 3, we can first make an overall qualitative comment. Perturbative corrections make the q T spectrum harder in going from NLL+NLO to NNLL+NNLO accuracy, and this occurs at both small and intermediate values of q T . NP effects increase the hardness of the q T spectrum at small values of q T and they are negligible at intermediate values of q T . Therefore, we note a non trivial interplay of perturbative and NP effects. In particular, at small values of q T higher-order perturbative contributions can be mimicked by NP effects.
At the quantitative level, in Fig. 3 we see that the NNLL+NNLO result supplemented with NP effects is very close to the perturbative result except in the very low q T region (q T ∼ < 3 GeV), i.e. below the peak of the q T distribution. In the region 3 GeV ∼ < q T ∼ < 10 GeV, the size of the NP band is similar to that of the PDF uncertainty band. At larger values of q T , the NP effects vanish (the size of the NP band is smaller than about 2% starting from q T ∼ 15 GeV).
We note that our simple model treats the regularization of the perturbative form factor (through the 'minimal prescription', see Sect. 2) and the NP form factor in an uncorrelated way, and this produces a conservative estimate of NP uncertainties. In other words, the model underestimates the potential of the resummed calculation at very small values of q T . For instance, the NP model can be improved by correlating the interplay between the perturbative form factor (and, e.g., its scale variation dependence) and the NP form factor (and the value of g N P ), and further constraints on the NP model can be possibly obtained by inputs from comparisons with experimental data.
In summary, from our brief discussion on the possible impact of NP effects for vector boson production at the LHC, we conclude that our conservative estimate leads to quantitative effects that are small and well within the scale variation dependence, still in the very low q T region. A quantitatively similar conclusion applies to the effect of PDF uncertainties. Based on these observations (and for practical purposes), in the presentation of our results of Sect. 3.2 we limit ourselves to considering only the perturbative calculation and the corresponding scale variation uncertainties. We conclude this subsection by presenting a comparison between the DYqT results and the results of the 'multidifferential' program DYRes. When no cuts are applied on the final-state leptons, the q T spectrum of the on-shell vector boson obtained with DYRes has to be in agreement with the one obtained with the numerical program DYqT. We have numerically checked that this is indeed the case. For illustrative purposes, we show the results of a comparison in Fig. 4. Here we consider the q T spectrum for on-shell Z boson production at the LHC with √ s = 7 TeV. The more clearly visible in the lower panel, which presents the result of the calculation of the binned ratio between the DYRes and DYqT results at both NLL+NLO and NNLL+NNLO accuracy. The ratio is everywhere consistent with unity within the numerical accuracy of its computation (the numerical errors in the computation of the binned ratio are below about 1% at small values of q T , and they are still below about 2% in the region 30 GeV ∼ < q T ∼ < 50 GeV where the value of the cross section sizeably decreases). We recall (see Sect. 2) that, at the inclusive level, the DYqT and DYRes calculations involve differences in the numerical implementation and two additional differences related to the treatment of the very low q T and high-q T regions. At very low values of q T , the difference is due to the regularization procedure of the perturbative form factor for very large values of the impact parameter b: the DYqT calculation uses the 'minimal prescription', while the DYRes calculation uses the b * freezing procedure. In our actual calculation with DYRes the value of b lim in Eq. (18) is set to b lim = b max , where b max is the maximum value of b that can be reached before encountering the singularity of the perturbative form factor (setting b lim = b max we do not introduce any additional regularization parameter, analogously to the case of the 'minimal prescription'). The value of b max depends on the renormalization and resummation scales µ R and Q and, in the case of Z and W production around the central value of the scales, the typical value is b max Q ∼ 1.2 · 10 3 µ R /m Z . We have checked that the 'minimal prescription' and the choice b lim = b max give basically the same numerical results, also at very small values of q T (q T ∼ 1 GeV). This numerical agreement is also visible (lower panel in Fig. 4) from the ratio between the DYRes and DYqT results at low values of q T .
At large values of q T , the DYRes calculation implements the smooth switching procedure of Eqs. (15)- (17). The large-q T region is shown in the inset plot of Fig. 4, and here the differences between the DYqT and DYRes calculations are due to the smooth switching procedure. In the high-q T region the DYRes result at NNLL+NNLO (NLL+NLO) accuracy basically agrees with the customary fixed-order result at O(α 2 S ) (O(α S )). The differences between the DYRes and DYqT results (consistently) decrease in going from NLL+NLO to NNLL+NNLO accuracy, and they are small at the NNLL+NNLO level. At both level of logarithmic accuracy, the DYRes and DYqT results agree within their corresponding scale variation uncertainties (which are not shown in the inset plot), and the DYRes result has a reduced scale dependence (it matches the scale dependence of the corresponding fixed-order result). The introduction of the smooth switching procedure in the DYRes calculation has practically a negligible quantitative effect on the unitarity constraint that is fulfilled by the DYqT calculation. In Table 1 we report the total cross sections for both Z and W production at √ s = 7 TeV, and we compare the resummed DYRes results with the corresponding fixed-order results obtained with the DYNNLO code. We see that the NLL+NLO (NNLL+NNLO) total cross section agrees with the NLO (NNLO) result to better than 1% accuracy.

Cross section [pb]
NLO NLL+NLO NNLO NNLL+NNLO pp → Z → l + l − 904.3 ± 0.2 904.6 ± 0.4 949.1 ± 0.7 947.3 ± 0.9 pp → W (±) → l (±) ν 9819 ± 2 9813 ± 4 10337 ± 6 10328 ± 9 In the main plot of Fig. 4, we also present a complementary information on the results of the fixed-order calculations (dashed lines) at O(α S ) (red dashed) and O(α 2 S ) (blue dashed). At intermediate values of q T the differences between the resummed results at two subsequent orders are smaller than the differences between the corresponding fixed-order results at two subsequent orders. The differences between the resummed results and the corresponding fixed-

Vector boson production at the LHC
In this Section we consider (q T related) physical observables that depend on the individual lepton momenta and on the kinematics of the lepton pair. The dependence can be indirect, through the application of acceptance cuts, and direct, through the definition of the observable. Therefore, the resummed calculation presented in this Section are performed by using the numerical program DYRes.
We start our presentation by considering the measurements of the q T spectrum of dilepton pairs at the LHC with √ s = 7 TeV, as reported by the CMS [91] and ATLAS [92] Collaborations with an integrated luminosity of 36 pb −1 and 4.7 fb −1 , respectively. The cuts that define the fiducial region in which the measurements are performed (our corresponding resummed calculation of the Z/γ * spectrum is carried out in the same region) are as follows. In the case of the CMS analysis the invariant mass m ll of the lepton pair is required to be in the range 60 GeV < m ll < 120 GeV, and the leptons must be in the central rapidity region, with pseudorapidity |η l | < 2.1, and they have a transverse momentum p l T > 20 GeV. In the case of the ATLAS analysis the fiducial region is defined by: 66 GeV < m ll < 116 GeV, |η l | < 2.4 and p l T > 20 GeV.
The results of our resummed calculation are shown in Fig. 5 (a) and (b). The blue-solid (red-dashed) histogram is the NNLL+NNLO (NLL+NLO) prediction for the q T spectrum, which is normalized to the cross section in the fiducial region, and the points are the data with the corresponding experimental errors. The inset plot shows the high-q T region. To facilitate the comparison between the data and the perturbative calculation we consider their ratio with respect to a reference theoretical result. We choose the NNLL+NNLO result at central values of the scales (µ F = µ R = Q = m Z /2) as reference result. The lower panel shows the data and the scale dependent NNLL+NNLO prediction normalized to this reference theoretical prediction. The scale dependence band of the perturbative calculation is computed by varying µ F , µ R and Q as previously discussed in Sect. 3.1: we vary µ F , µ R and Q in the range m Z /4 ≤ {µ F , µ R , Q} ≤ m Z , with the constraints 0.5 ≤ {µ F /µ R , Q/µ R } ≤ 2. We see that our perturbative calculation is consistent with the data within the uncertainties. The scale variation bands at NLL+NLO and NNLL+NNLO accuracy overlap. Moreover, in going from NLL+NLO to NNLL+NNLO accuracy the perturbative uncertainty is reduced and the agreement between experimental data and theory prediction is improved. The perturbative uncertainty at NNLL+NNLO accuracy is about ±10% at the peak, it decreases to about ±4% at q T ∼ 10 GeV, and it increases again to about ±10% at q T = 40 GeV. The comparison between our theoretical prediction and the CMS and ATLAS data is qualitatively similar, the main difference being that, due to the larger data sample, the experimental errors in the ATLAS analysis are significantly smaller.
We add a comment on the large-q T region (see inset plots of Fig. 5), where the cross section is dominated by the fixed-order contribution. For very large q T , i.e. q T ≫ m Z , the physical hard scale of the process is of the order of q T and not of the order of m Z , and a sensible scale choice is µ F ∼ µ R ∼ q T . Therefore, it is not unexpected that our NNLL+NNLO calculations, which use µ F ∼ µ R ∼ m Z /2, slightly overshoot the CMS and ATLAS data in the last few high-q T bins. The size of the QCD corrections evaluated with µ F ∼ µ R ∼ q T would be smaller. Moreover, in the extreme region q T ≫ m Z a resummation of enhanced large-q T perturbative terms is in principle required [87,88].
In Fig. 6 we consider the q T spectrum of W ± bosons. We present a comparison of our resummed results with the pp → W → lν data collected by the ATLAS Collaboration [93] with an integrated luminosity of 31 pb −1 at √ s = 7 TeV. The fiducial region is defined as follows: the charged lepton has transverse momentum p l T > 20 GeV and pseudorapidity |η l | < 2.4, the missing transverse energy is p ν T > 25 GeV, and transverse mass m T = 2p l T p ν T (1 − cos(φ l − φ ν )) is constrained in the region m T > 40 GeV. We recall that, because of the presence of the neutrino in the final state, the q T of the W has to be reconstructed through the transverse energy of the hadronic recoil, which has a poorer experimental resolution than that of the lepton momentum. In the small q T region, the bin sizes of the experimental data are rather large, with only four bins in the region with q T < 55 GeV. For this reason in Fig. 6 we focus on the large q T region 55 GeV < q T < 300 GeV, while the small q T region is shown in the inset plot. The lower panel of Fig. 6, which covers the entire q T region of the data, presents the ratio of both data and theoretical results with respect to the reference theoretical result. This ratio and the scale variation bands are computed exactly in the same manner as in the case of Fig. 5. Looking at the ratio plot in the lower panel, we see that our NNLL+NNLO calculation describes the W production data within the perturbative uncertainties. The NNLL+NNLO perturbative uncertainty is about ±8% at the peak, it decreases to about ±4% at q T ∼ 15 GeV, and it increases again to about ±15% at q T = 50 GeV.
In Sect. 3.1 and in the first part of this Section, we have examined vector boson q T distributions (without and with the application of acceptance cuts) and we have computed and studied the effects that are produced by the all-order resummation of large logarithmically-enhanced terms Figure 6: Vector boson production at the LHC with lepton selection cuts. The NLL+NLO (red) and NNLL+NNLO (blue) normalized q T spectra for W ± production are compared with the ATLAS data of Ref. [93]. The ratio in the lower panel and the scale variation bands are obtained as in Fig. 5.
at small values of q T . Our related calculations are performed at complete NNLL+NNLO (and NLL+NLO) accuracy. In the following part of this Section, we consider other observables that are related to the q T distributions but in which fixed values of q T are not directly measured. These observables are inclusive over q T within certain q T ranges. Since the bulk of the vector boson cross section is produced at small values of q T , if the observable (indirectly) probes the detailed shape of the production cross section in the small-q T region, the observable itself can be very sensitive to high-order radiative corrections and to the q T resummation effects that we can explicitly compute. This reasoning illustrates and justifies the physical (and quantitative) relevance of q T resummation for other q T -related observables. In the second part of this Section we study the quantitative impact of q T resummation on some observables.
At the formal level, our study of other observables implies that we are resumming high-order logarithmic corrections (in case they are present) that appear in the computation of those observables. Strictly speaking, this resummation has to be performed on an observable-dependent basis (see, e.g., Ref. [96]). Therefore, our observable-independent treatment (based on transversemomentum resummation) cannot guarantee that we formally achieve exact NNLL+NNLO accuracy for all these observables. Nonetheless we are able to correctly take into account all the leading-logarithmic contributions, all the complete (with and without logarithmic enhancement) perturbative terms up to the NNLO level ‡ ‡ , and a substantial part of subleading logarithmic terms beyond the NNLO accuracy.
This statement about resummation is a consequence of the following discussion. The observabledependent logarithmic terms (in case they are present) are due to multiple radiation of soft and collinear partons in the inclusive final-state: these logarithmic corrections are computed by approximating the QCD scattering amplitudes in the soft and collinear limits and, then, by integrating the final-state QCD radiation over the corresponding phase space with appropriate (observabledependent) kinematical approximations. In our transverse-momentum resummation procedure we correctly take into account the NNLL dynamics (the behaviour of the QCD scattering amplitudes in the soft and collinear limits) of soft and collinear radiation, and we treat the phase space of the final-state QCD radiation with consistent kinematical approximations that are specific of the q T spectrum. However, the observable-dependent kinematical approximations can only differ beyond the leading-logarithmic level (to leading-logarithmic level, a strong-ordering approximation in the energy/angle of the emitted partons is sufficient), and these differences do not spoil the leadinglogarithmic accuracy of our resummed calculation. In this respect, it is important to remark the role of the q T recoil (see Sect. 2 and Appendix A) on the kinematics of the produced (observed) lepton pair. We treat the q T recoil in a kinematically consistent way (though it necessarily involves non logarithmic approximations that are uniformly of O(q T /M) throughout the small-q T region), and such a treatment is necessary to correctly correlate the dynamical q T resummation effect with the ensuing q T dependence of the measured (computed) observable.
In summary, the application of our q T resummed calculations to the computation of other observables is physically (and, thus, quantitatively) and formally (as we have just discussed) justified. A detailed specification of the subleading-logarithmic accuracy of the q T resummed calculation at the formal (analytical) level requires (and deserves) observable-dependent investigations, which can be performed in future studies.
Among other observables, we first consider the measurement † of the φ * distribution from pp → Z/γ * → l + l − data at √ s = 7 TeV as reported by the ATLAS Collaboration [94] with an integrated luminosity of 4.6 fb −1 . The φ * observable is defined as φ * = tan(π/2 − ∆φ/2) sin(θ * ), where ∆φ is the azimuthal angle between the leptons and the angle θ * is defined by cos θ * = tanh((η l + − η l − )/2) where η l + (η l − ) is the rapidity of the positively (negatively) charged lepton. The cuts that define the fiducial region are those of the ATLAS analysis of the q T spectrum: 66 GeV < m ll < 116 GeV, p l T > 20 GeV and |η l | < 2.4.
The φ * variable at small values of φ * is correlated to q T and, therefore, it is strongly sensitive to q T resummation effects. A detailed discussion on the relation between φ * and q T is ‡ ‡ For observables that are inclusive over the region that includes q T = 0, the NNLO accuracy is achieved through our detailed matching procedure (see Sect. 2) with the fixed-order calculation.
† An analogous measurement of the φ * distribution at the LHC was reported by the LHCb Collaboration [95] with an integrated luminosity of 0.94 fb −1 . In the small-φ * region, the bin sizes of the LHCb measurement are rather large (with respect to those of the ATLAS measurement [94]), with only two (four) bins in the region φ * < 0.1 (φ * < 0.2).

Figure 7:
The NLL+NLO (red) and NNLL+NNLO (blue) normalized φ * distribution for Z/γ * production at the LHC is compared with the ATLAS data of Ref. [94]. The NLL+NLO and NNLL+NNLO central results are computed at the scales µ R = µ F = Q = m Z /2. The ratio in the lower panel and the scale variation bands are obtained as in Fig. 5.
presented in Ref. [96], where the resummation of the ln φ * terms is carried out in analytic form up to NNLL+O(α 2 S ) accuracy ‡ , and it turns out to be strictly related and very similar to q T resummation. Ensuing phenomenological studies are presented in Refs. [67,97,68].
In Fig. 7 we report the ATLAS data of the φ * distribution (normalized to the measured cross section in the fiducial region) and the comparison with the results of our resummed calculation. The NLL+NLO and NNLL+NNLO central results are computed at the scales µ R = µ F = Q = m Z /2. The scale variation bands at NLL+NLO (red) and NNLL+NNLO (blue) accuracy and ‡ The analytical treatment of Ref. [96] does not reach complete NNLO accuracy at small values of φ * since the analogue of the vector boson coefficient H the reference NNLL+NNLO result for the ratio in the lower panel of Fig. 7 are computed as in Figs. 5 and 6. We observe that the scale variation bands at the two subsequent orders overlap, and that the NNLL+NNLO perturbative uncertainty is substantially smaller than the NLL+NLO one. The NNLL+NNLO result is consistent with the data within the uncertainties in both the small-φ * and large-φ * regions (the large-φ * region is shown in the inset plot). The NNLL+NNLO perturbative uncertainty is about ±10% for φ * < 0.01, it decreases to about ±5% at φ * ∼ 0.05, and it increases again to about ±10% at φ * ∼ 0.2.
We add a comment on the results that we have shown in Figs. 5-7. We recall that all the results presented in this Section are obtained in a purely perturbative framework. In Sect. 3.1 we have discussed the possible impact of the inclusion of a NP form factor, and we have seen (Fig. 3) that NP effects should lead to a deformation of the perturbative result that is well within the scale variation uncertainties of the NNLL+NNLO calculation. In Figs. 5-7 we observe that all the resummed perturbative predictions are consistent with the data within our estimation of perturbative uncertainties. Owing to the agreement between the theoretical NNLL+NNLO predictions and the experimental data in the very small q T /φ * region, we cannot draw any precise quantitative conclusion about the definite size of NP effects in the Z/γ * , W ± and φ * distributions that we have considered. We can only conclude that NP effects have to be small in order not to spoil the agreement between the data and the corresponding NNLL+NNLO results in Figs. 5-7. We conclude this Section by considering other observables. We study the impact of q T resummation on the kinematical distributions that are relevant for the measurement of the W mass. We consider pp → W − → l −ν l with √ s = 7 TeV and we apply the following selection cuts: the charged lepton has transverse momentum p l T > 30 GeV and rapidity |η l | < 2.4, the missing transverse momentum is p ν T > 30 GeV, and the transverse mass m T has m T > 60 GeV. We also apply a cut, p W T < 30 GeV, on the transverse momentum p W T of the W boson (lepton pair). The results of our calculation of the m T distribution and of the lepton momentum distributions are presented in Fig. 8 and 9, respectively. The reference scale choice of the calculation is µ F = µ R = Q = m W /2. In both figures we present the results of the fixed-order calculation at LO (cyan dotted), NLO (green solid) and NNLO (black dot-dashed) accuracy and we compare them with the results of the q T resummed calculation at NLL+NLO (red dashed) and NNLL+NNLO (blue solid) accuracy. The lower panels show the ratio between the various results and the NNLL+NNLO result (the ratio LO/(NNLL+NNLO) is not reported in the lower panels).
The m T distribution in the range m T < 90 GeV is presented in Fig. 8. We can consider two regions: the large-m T region, around m T ∼ m W (we recall that we use m W = 80.385 GeV), and the small-m T region. In the large-m T region, m T ∼ > 70 GeV, we see that the perturbative prediction is extremely stable against radiative corrections, and the stability is present both in going from NLO to NNLO accuracy and with inclusion of resummation. This is a consequence of the well known fact that the transverse mass is weakly sensitive to the transverse momentum of the W boson. Formally, the m T distribution has no logarithmic corrections of the type ln(|m T − m W |/m W ), and our q T resummed calculation does not spoil the stability of the fixed-order expansion. On the contrary, in the small-m T region, we observe that the fixed-order predictions become unreliable. The LO distribution is large at m T = 60 GeV, and both the NLO and NNLO distributions become negative at m T ∼ 60 GeV. This (mis-)behaviour is due to the fact that the constraints p l T > 30 GeV and p ν T > 30 GeV produce an unphysical boundary (and a stepwise behaviour) of the m T distribution at m T = m T step = 60 GeV in the LO calculation. The boundary is due to the LO kinematics p l T + p ν T = q T = 0, and it disappears at higher orders since q T = 0. The LO boundary induces (integrable) logarithmic singularities of the type ln(1 − m T step /m T ) 2 at NLO and beyond [98]. These logarithmic terms are resummed to all order by q T resummation, and the singularities are absent in the resummed prediction [98], which is well behaved at the LO boundary m T = m T step . We also note that the differences between the NLL+NLO and NNLL+NNLO results are small at m T ∼ 60 GeV.
In Figs. 9 (a) and (b) we present the p l T and p ν T distributions, respectively. In the limit in which the W boson is produced on shell, these distributions have an LO kinematical boundary at m W /2. The finite width of the W boson (partially) smears this effect: at LO both the p l T and p ν T (a) (b) Figure 9: Effect of q T resummation for pp → W − → l −ν l production at the LHC: (a) lepton p T distribution and (b) missing p T distribution. The fixed-order and resummed results are denoted as in Fig. 8.
distributions are strongly peaked at m W /2 (Jacobian peak) and quickly drop for p T ∼ > m W /2. The almost stepwise behaviour of the LO distribution produces large radiative corrections at NLO and beyond (in the limit in which the W boson is produced on shell, these large corrections would be integrable logarithmic singularities at each perturbative order [98]). The NLO and NNLO distributions indeed display an unphysical peak at p T ∼ 42 GeV, which is an artifact of such large corrections (singularities in the on-shell limit). The resummed predictions at NLL+NLO and NNLL+NNLO accuracy are free of such instabilities and display a smooth shoulder behaviour around the LO boundary for on-shell production. The perturbative instabilities of the fixed-order calculation at small values of p T (p l T ∼ 30 GeV and p ν T ∼ 30 GeV) are analogous to those that we have previously discussed in the case of the m T distribution in the region m T ∼ 60 GeV (see Fig. 8). In the case of the p T distributions, it is the constraint m T > 60 GeV that produces the LO boundaries at p l T = p ν T = 30 GeV, an LO stepwise behaviour and ensuing instabilities at each subsequent perturbative order. The resummed calculation is perturbatively stable in the small-p T region, and the differences between the NLL+NLO and NNLL+NNLO results are small throughout the entire region with p T ∼ < 45 GeV. In the large-p T region (p T ∼ > 45 GeV) both the p l T and p ν T distributions display radiative corrections that are relatively large. This is not unexpected since in this region of transverse momenta the NLO calculation is essentially the first perturbative order at which both the p l T and the p ν T distributions are non vanishing (in the on-shell limit, the O(α S ) and O(α 2 S ) result would be an LO and an NLO prediction, respectively).

Summary
In this paper we have considered the transverse-momentum (q T ) distribution of DY high-mass lepton pairs produced, via Z/γ * and W bosons decay, in hadronic collisions. We have presented a perturbative QCD study based on transverse-momentum resummation up to NNLL accuracy.
We have combined small-q T resummation with the known O(α 2 S ) fixed-order result at small, intermediate and large values of q T .
We have followed the resummation formalism developed in Refs. [25,26,27] to implement transverse-momentum resummation and the matching with the result at O(α 2 S ). In particular, our calculation includes the complete NNLO contributions at small values of q T (i.e., in any regions that include q T = 0) and it exactly reproduces the complete NNLO total cross section after integration over q T . This leads to theoretical 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 predictivity of small-q T resummation is superseded by that of the customary fixed-order expansion, and our resummed calculation can be smoothly joined onto the O(α 2 S ) calculation. The resummed calculation can be systematically expanded at various orders of logarithmic accuracy (e.g., NLL+NLO and NNLL+NNLO accuracy), and its theoretical uncertainties due to uncalculated higher-order QCD corrections can be studied by comparing the results at two subsequent orders and by performing systematic studies on factorization, renormalization and resummation scale dependence.
We have performed such a study for the case of vector boson production at LHC energies, and we have briefly illustrated the uncertainties due to parton densities and the possible impact of non-perturbative effects.
In the present paper we have extended the resummed calculation presented in Ref. [40] for Z/γ * production by considering also W ± production and by including the leptonic decay of the vector boson with the corresponding spin correlations, the finite-width effects and the full dependence on the final-state leptonic variables. We have compared our resummed results for Z/γ * and W production with some of the available data of the ATLAS and CMS experiments at the LHC, applying the same kinematical cuts on final state leptons that are considered in the experimental analyses. We find that the data are well described by our predictions within the perturbative uncertainties. We have also considered the impact of transverse-momentum resummation on observables, which are different from the vector boson q T , that depend on the lepton kinematical variables. In particular, we have studied the φ * distribution in Z/γ * production and the leptonic transverse-momentum, the missing transverse-momentum and the transverse-mass distributions in W production.
Our calculation is implemented in the parton-level Monte Carlo numerical code DYRes which allows the user to apply arbitrary kinematical cuts on the vector boson and the final-state leptons, and to compute the corresponding relevant distributions in the form of bin histograms. These features make our program a useful tool for DY studies at the Tevatron and the LHC. A version of the DYRes code is publicly available.
The production and decay mechanisms of the vector boson are dynamically correlated by the non vanishing spin of the vector boson. The inclusion of the lepton decay (with the spin correlations and the full dependence on the kinematical variables of the two leptons) in the resummed calculation requires a general theoretical discussion on the q T recoil due to the transverse mo-mentum of the vector boson. This discussion is not limited to the specific case of vector boson production. We have presented a general and explicit procedure to treat the q T recoil. The procedure is directly applicable to q T resummed calculations for production processes of generic high-mass systems in hadron collisions.
Acknowledgements. We would like to thank Stefano Camarda, Luca Perrozzi and Jan Stark for useful discussions. This research was supported in part by the Swiss National Science Foundation ( A Appendix: Lepton angular distribution and q T recoil in transverse-momentum resummation This Appendix is devoted to the q T -recoil issue that we have introduced and illustrated in Sect. 2 (see Eqs. (12)- (14) and accompanying comments). To our knowledge the issue has not received much attention in the previous literature on transverse-momentum resummation. We present a detailed discussion of the issue and a general, explicit and consistent procedure to implement the q T recoil in transverse-momentum resummation. Our procedure explicitly exhibits the degree of freedom involved in the implementation of the q T recoil and, moreover, it gives an explicit formal parametrization of the ensuing ambiguities. The procedure is straightforwardly applicable to implement the q T recoil (and, possibly, estimate related uncertainties) in calculations based on transverse-momentum resummation.
The q T -recoil issue is not specific of the lepton angular distribution for vector boson decay, but it regards transverse-momentum resummation for generic production processes. For simplicity of presentation, in the following we consider in detail vector boson production and the DY process. Then we discuss the generalization to generic processes.
We begin our discussion by considering the computation of the DY multidifferential cross section in Eq. (3) at the LO in perturbative QCD. At this order the hadronic cross section (and the corresponding partonic cross section) is directly and exactly (i.e., with no small-q T approximation) proportional to the Born level angular distribution dσ (0) /dΩ in Eq. (5). We have For the purpose of our general discussion of q T recoil, we write the lepton angular distribution in the following form: Note that, following the general notation of Eq. (3), we have not specified the actual definition of the angular variables Ω, and we have written the left-hand side of Eq. (20) in a Lorentz invariant form. The relation (20) is written in the form of a proportionality relation: the additional proportionality factors that are not explicitly denoted in the right-hand side are not relevant for our following discussion (in particular, they are independent of the lepton momenta {p 3 , p 4 } and, thus, of Ω). The factor |M for the partonic process where k i (i = 1, 2) is the momentum of the colliding parton a i from the initial-state hadron h i (P i ) (see Eq. (1)), with the kinematics In our specific case of vector boson production, the Born level partonic process is the qq annihilation process q fqf ′ → V → l 3 l 4 (i.e., {a 1 , a 2 } = {q f ,q f ′ }). All the other factors in the right-hand side of the relation (20) are related to the kinematical phase space of the final-state leptons and, in particular, they enforce the kinematical constraint q = p 3 + p 4 .
The LO calculation of the cross section kinematically relates the parton and hadron momenta k i and P i . In particular, at the LO we have q T = 0 and, specifically, the LO value with where (see Sect. Higher-order perturbative contributions produce logarithmically-enhanced ('singular') terms at small q T that can be resummed to all orders, leading to the resummation factorŴ in Eq. (5). These logarithmic terms are due to multiple radiation of soft and collinear partons, and this soft and collinear radiation is factorized [29] with respect to the Born level amplitude M (0) a 1 a 2 →l 3 l 4 of Eq. (20). As a consequence, after q T resummation the angular distribution of the decaying leptons is still given by the Born level function dσ (0) /dΩ in Eq. (20), and this function thus appears as a multiplicative factor in front of the resummation factorŴ of the resummed component of the vector boson q T cross section (see Eq. (5)). Strictly speaking [29], in the limit q T ≪ M that is relevant for resummation, the angular distribution can be expressed in terms of the LO distribution dσ (0) /dΩ LO in Eq. (19), namely the expression (20) with the LO kinematics of Eqs. (23) and (24), which in particular has q T = 0. Indeed, after soft/collinear factorization and resummation, any residual dynamical effect on the process in Eq. (21) (and on M (0) and dσ (0) /dΩ) is due to hard-parton radiation. Hard radiation produces O(q T /M) corrections that lead to non-singular contributions if q T ≪ M: these corrections can be formally approximated by their limiting behaviour as q T → 0 and, thus, neglected in the computation of the resummed component (see Eqs. (5) and (11)) and included in the finite component (see Eq. (4)).
Neglecting these O(q T /M) corrections is a perfectly suitable procedure for the resummed calculation of the vector boson q T cross section (see Eq. (11)). However, performing the resummation at fixed lepton momenta, the momentum of the vector boson must be fully specified by the lepton momenta and, in particular, q T = p T 3 + p T 4 is not vanishing. The resummation factorŴ (see Eq. (5)) produces a smearing of the LO distribution δ (2) (q T ) of Eq. (19) and finite values of q T : to avoid unphysical results (e.g., events with q T = 0 and p T3 + p T4 = 0) the factor dσ (0) /dΩ in Eq. (5) cannot be the LO angular distribution dσ (0) /dΩ LO (which has p T 3 +p T 4 = 0) in Eq. (19). In other words, the non-vanishing value of q T has to be distributed between the two lepton momenta and this leads to the q T -recoil issue that we have illustrated in Sect. 2 (see Eqs. (12)- (14) and accompanying comments). The resummed calculation requires the specification of a q T -recoil prescription that has to be consistent (and physically sensible), although this can be done in many (infinitely many) different ways.
Actual resummed calculations performed in the literature do not mention the q T -recoil issue. The calculations of Refs. [55,56,57] directly refer to the use of the Collins-Soper (CS) rest frame [70]. The procedure to compute the factor dσ (0) /dΩ in Eq. (5) is as follows. The lepton angular variables Ω are specified to be the polar and azimuthal angles {θ ′ CS , φ ′ CS } of one of the leptons in the CS rest frame. The LO distribution dσ (0) /dΩ LO in Eq. (19) is then expressed in terms of {θ ′ CS , φ ′ CS } (since the LO distribution has q T = 0, in this case {θ ′ CS , φ ′ CS } exactly coincide with the lepton scattering angles in the centre-of-mass frame of the LO colliding parton momenta in Eq. (23)) and this leads to an unambiguously defined angular function F q fqf ′ →l 3 l 4 (θ ′ CS , φ ′ CS ) (this function is actually independent of φ ′ CS ) that is used to define (see Eq. (12)) the angular distribution dσ (0) /dΩ of the resummed component of the cross section (see Eq. (5)). This is a perfectly defined and consistent procedure, but it hides the actual implementation of O(q T /M) corrections through an implicit prescription for the q T recoil: the definition of the CS rest frame is q T dependent and a q T dependence is introduced by identifying/equating the angles {θ ′ CS , φ ′ CS } of the LO and resummed calculations (additional comment on this are presented in a paragraph after Eq. (32)).
Here we explicitly present a consistent q T -recoil procedure and an entire class of q T -recoil prescriptions. Our viewpoint is as follows: the non-vanishing value of q T of dynamical origin that is produced by resummation leads to a q T -recoil that can be 'kinematically absorbed' ¶ by the momenta k 1 and k 2 of the colliding partons of the underlying hard-scattering process (see Eq. (21)). As specified below, there are infinitely-many ways of implementing this kinematical recoil on the colliding partons in a consistent manner (i.e., without modifying the logarithmically-enhanced perturbative terms at small q T ): they differ by corrections that are of O(q T /M) order-by-order in the perturbative expansion (after having matched the resummed calculation with the complete N k LO calculation, as in Eqs. (4) and (10), these corrections start to contribute at the N k+1 LO level). ¶ The q T recoil issue does not arise in the context of transverse-momentum (k T ) factorization [80] for high-energy (small-x) hard-scattering processes. In this formulation the q T recoil is dynamically (and uniquely) embedded in the factorization formula. The parton densities of the colliding hadrons are k T dependent and the hard-scattering colliding partons have ensuing non-vanishing transverse momenta k iT that enter as integration variables in the factorization formula.
According to our procedure, the lepton angular distribution dσ (0) /dΩ to be used in the resummed calculation (see Eq. (5)) is exactly given by the expression in Eq. (20). The phase space factor in the right-hand side of Eq. (20) is directly given in terms of the physical (measured) lepton momenta p 3 and p 4 (with p 3 + p 4 = q). The momentum k 1 (then, k 2 = q − k 1 ) to be used to compute the Born level scattering amplitude in Eq. (20) is given by the following parametrization: where and k µ 1T is a two-dimensional vector that is transverse to both P µ 1 and P µ 2 (i.e., k 1T lies in the q T plane) and that fulfils the following constraints: We note that, following the definition in Eqs. (25) and (26), k µ 1 and k µ 2 are well defined 'physical' parton momenta: they fulfil the kinematics in Eq. (22) and they have positive definite energies, k 0 1 > 0 and k 0 2 > 0 (the constraint in Eq. (28) guarantees that the four-momentum k µ 1 has positive definite energy and, then, k 0 2 > 0 follows from q 0 > 0). Therefore the scattering amplitude M We also note that k µ 1 is invariant under longitudinal boosts of the hadronic centre-of-mass frame, provided k µ 1T is boost invariant.
At fixed values of q µ , P µ 1 and P µ 2 , Eqs. (25)- (28) give the most general expression of k 1 that respects the Born level kinematics in Eqs. (21) and (22) and the LO kinematics in Eqs. (23) and (24). This expression is parametrized by the arbitrary (though constrained) transverse-momentum vector k 1T . By choosing different values of k 1T , we can obtain an entire class of consistent q T -recoil prescriptions. For example, two 'obvious' possible choices are as follows: A) set k 1T = q T /2 (and thus k 2T = q T /2): from Eq. (26) we obtain and we have B) set k 1T = 0 (and thus k 2T = q T ): from Eq. (26) we obtain z 1 = 1 and we have We also note that, after integration over the lepton angular variables Ω, we consistently obtain the Born level total cross sectionσ (0) a 1 a 2 →l 3 l 4 (M 2 ) of the resummation formula (11). Indeed, after the Ω integration of Eq. (20), the result does no longer depend on the lepton momenta and, since it is a Lorentz invariant quantity, the result can only depends on the invariant (k 1 + k 2 ) 2 = 2k 1 · k 2 = q 2 = M 2 , which is independent of k 1T . In other words, the dependence on the arbitrary parameter k 1T completely cancels in lepton-inclusive observables.
Using our q T -recoil procedure, we can compute the corresponding lepton angular function F q fqf ′ →l 3 l 4 of Eq. (12). This function is the product of two factors. One factor is a purely kinematical origin (it derives from the phase space factor in the right-hand side of Eq. (20)), and it depends on the specification of the angular variables Ω. The other factor, denoted as F (D) in the following (for simplicity we omit the subscript q fqf ′ → l 3 l 4 ), has a dynamical origin and it depends on the Born level factor |M (0) (k 1 , k 2 ; p 3 , p 4 )| 2 in the right-hand side of Eq. (20). Since |M (0) (k 1 , k 2 ; p 3 , p 4 )| 2 is a Lorentz invariant scalar quantity and the momenta {k 1 , k 2 ; p 3 , p 4 } are constrained by momentum conservation, F (D) can only depend on the dimensionless variable 4k 1 · p 3 /(2k 1 · k 2 ) = 4k 1 · p 3 /M 2 . Considering the centre-of-mass frame of k 1 and k 2 , we have 4k 1 · p 3 /M 2 = 1 − cos θ ′ 13 , where θ ′ 13 is the scattering angle between k 1 and p 3 . In other words, F (D) = F (D) (θ ′ 13 ) and θ ′ 13 is the lepton scattering angle in a particular rest frame of the vector boson momentum q µ (the centre-of-mass frame of k 1 and k 2 ). Our q T -recoil procedure can thus be reinterpreted in terms of generation of lepton-pair events. Considering a definite (with respect to the hadronic collision frame) rest frame of the vector boson momentum, the lepton-pair event and the individual lepton momenta are generated in that frame according to the corresponding Born level angular distribution; then the lepton pair distribution is boosted to the hadronic collision frame through the corresponding Lorentz transformation. Since there is an infinite numbers of vector boson rest frames, this event-generation procedure has an infinite degree of arbitrariness. Applying a three-dimensional rotation to a vector boson rest frame, we obtain another vector boson rest frame and, thus, the infinite numbers of vector boson rest frames depends on the two scalar parameters of the three-dimensional rotation. Accordingly, our q T -recoil procedure depends on the arbitrary two-dimensional vector k 1T , namely on two parameters (the magnitude |k 1T | and the azimuthal angle of k 1T ). In other words, the relation between the LO momenta k µ i (LO) in Eq. (23) and the q T -recoiled momenta k i obtained through Eqs. (25)- (28) can be reinterpreted as a Lorentz transformation of the colliding parton momenta from the hadronic collision frame to a specified vector boson rest frame. This interpretation directly relates our q T -recoil procedure with the specific CS frame procedure (as already mentioned and described in the initial part of this Appendix) that is directly used in other resummed calculations [55,56,57]. It can be explicitly checked that the CS frame procedure used to specify the lepton angular distribution dσ (0) /dΩ in the resummed calculation of Refs. [55,56,57] corresponds to the choice k 1T = k 2T = q T /2 (see Eqs. (30) and (31)) within our class of q T -recoil prescriptions.
Owing to our explicit parametrization and implementation of the q T -recoil procedure, the quantitative effects of various q T -recoil prescriptions can be directly investigated in applications of the numerical program DYRes (and of other resummed calculations). Obviously (as discussed in Sect. 2), different q T -recoil prescriptions have no effects on quantities that are fully inclusive over the leptonic variables Ω. In general, our expectations are as follows. We expect that the quantitative differences produced by various q T -recoil prescriptions are small for lepton non-inclusive observables that are mostly sensitive to either the small-q T region (in this region the q T -recoil effects are non-singular and thus subdominant with respect to the singular logarithmic contributions) or the high-q T region (in this region the q T -recoil effects are suppressed by the smooth switching procedure of Eqs. (15)-(17)), while relatively larger differences can appear in case of sensitivity to the region of intermediate values of q T . Moreover, these quantitative differences are expected to decrease in going from NLL+NLO to NNLL+NNLO accuracy, since the non-vanishing q T -recoil effects start to formally contribute at the N k+1 LO level in the N k LL+N k LO calculation. In Sect. 3.2 we have presented our quantitative results obtained with the DYRes code. As stated at the beginning of Sect. 3, these results are obtained (analogously to those in Refs. [55,56,57]) by computing the Born level angular distribution dσ (0) /dΩ in the CS rest frame, i.e. by setting k 1T = k 2T = q T /2 in the actual implementation of our q T -recoil procedure (this corresponds to use the prescription A in Eqs. (29)-(31)). Setting k 1T = q T /2, we have also considered other variants of the q T -recoil prescriptions and we have examined the quantitative differences that are produced on the observables that are examined in Sect. 3.2 (i.e., the observables in Figs. 5-9). We have found that various q T -recoil prescriptions produce differences that are in agreement with our general expectations and, in particular, at NNLL+NNLO accuracy these differences lead to small quantitative effects: typically, the effects are much smaller than the scale-variation uncertainties (estimated as in Sect. 3.2). For instance, comparing the q T -recoil prescriptions A (see Eqs. (29)-(31)) and B (see Eq. (32)), we obtain quantitative differences that are at most at the percent level (e.g., in the case of the φ * distribution of Fig. 7 in the region 0.3 ∼ < φ * ∼ < 1, and in the case of the lepton-p T and missing-p T distributions of Fig. 9 in the region 45 GeV ∼ < p T ∼ < 50 GeV): these differences are definitely smaller than the scale uncertainty (at both the NLL+NLO and NNLL+NNLO levels) and, especially, they are of the same size as (and, hence, hardly distinguishable from) the pure numerical errors of the DYRes calculation in the NNLL+NNLO mode.
The q T -recoil issue that we have discussed in this Appendix is not specific of vector boson production and the ensuing leptonic decay. The issue affects q T resummed calculations for any process of the type h 1 + h 2 → F (p 3 , p 4 , p 5 , . . . ) + X (we use the same notation as in Eq. (1)) where the final-state high-mass system F has total transverse momentum q T and the momenta p 3 , p 4 , p 5 , . . . of its 'decay products' are directly measured. Owing to the universality (processindependent) structure of transverse-momentum resummation [29], the q T -recoil procedure that we have introduced in this Appendix is directly applicable to all these processes. The only key difference with respect to vector boson production is that the Born level scattering amplitude M (0) (k 1 , k 2 ; p 3 , p 4 ) in Eq. (20) is replaced by a properly computable (all-loop) hard-virtual amplitude M(k 1 , k 2 ; p 3 , p 4 , p 5 , . . . ; α S (M 2 )) (see Sect. 4 in Ref. [29]), which embodies QCD virtual radiative corrections ( M is computable as power series in α S (M 2 )). Strictly speaking [29], the q T resummed cross section at small values of q T is proportional to the angular dependent distribution of the momenta {p 3 , p 4 , p 5 , . . . } as computed from M at q T = 0 (i.e., with the LO momenta k µ i (LO) of Eq. (23)). The ensuing q T -recoil issue can be directly solved by our q T -recoil procedure. Indeed, the hard-virtual amplitude M (k 1 , k 2 ; p 3 , p 4 , p 5 , . . . ; α S (M 2 )) has the same kinematical properties as its Born level counterpart M (0) = M (0) : therefore, the q T recoil can be directly implemented by simply evaluating M (k 1 , k 2 ; p 3 , p 4 , p 5 , . . . ; α S (M 2 )) with the q T -recoiled momenta k i of Eqs. (25)- (28).
We add some final comments on spin correlations and on the specific process of SM Higgs boson production and its decay in colourless particles (e.g., H → γγ, H → W W → ℓνℓν, H → ZZ → 4ℓ) [43]. The q T resummed Higgs boson cross section at fixed momenta of the decay products is proportional to the angular distribution as obtained (analogously to Eq. (20)) by the corresponding Born level scattering amplitude M H→l 3 l 4 l 5 ... only depends on observable momenta, while |M (0) g 1 g 2 →H (k 1 , k 2 ; q)| 2 only depends on (k 1 + k 2 ) 2 = q 2 because of Lorentz invariance. As a consequence, our q T -recoil procedure (and its dependence on the definition of k 1 and k 2 ) has no effect on the angular distribution of the Higgs boson decay products. The angular distribution of the resummed calculation can be directly obtained [43] by supplementing the Born level total cross sectionσ (0) gg→H (M 2 ) with the (kinematical and dynamical) Higgs boson decay factor. This specific example also clearly illustrates that the q T -recoil issue that we have introduced in Sect. 2 and discussed in this Appendix is directly related and due to the vector boson spin and the spin correlations between the production and decay subprocesses of the vector boson. This kind of relation between q T recoil and spin correlations is valid for generic production processes.