Drell-Yan lepton-pair production: $q_T$ resummation at N$^3$LL accuracy and fiducial cross sections at N$^3$LO

We present high-accuracy QCD predictions for the transverse-momentum ($q_T$) distribution and fiducial cross sections of Drell-Yan lepton pairs produced in hadronic collisions. At small value of $q_T$ we resum to all perturbative orders the logarithmically enhanced contributions up to next-to-next-to-next-to-leading logarithmic (N$^3$LL) accuracy, including all the next-to-next-to-next-to-leading order (N$^3$LO) (i.e. $\mathcal{O}(\alpha_S^3)$) terms. Our resummed calculation has been implemented in the public numerical program DYTurbo, which produces fast and precise predictions with the full dependence on the final-state leptons kinematics. We consistently combine our resummed results with the known $\mathcal{O}(\alpha_S^3)$ fixed-order predictions at large values of $q_T$ thus obtaining full N$^3$LO accuracy also for fiducial cross sections. We show numerical results at LHC energies discussing the reduction of the perturbative uncertainty with respect to lower-order calculations.

After the successful operation of the first two runs of the Large Hadron Collider (LHC) at CERN and the discovery of the long sought Higgs boson, a major task of the high-energy physics community has become a direct investigation of the electroweak symmetry breaking mechanism. In the absence of clear direct signals of new physics phenomena, precision studies give us a unique opportunity to search for possible deviations from Standard Model (SM) predictions. In this scenario it is clear that theoretical predictions for SM cross sections and associated distributions at an unprecedented level of accuracy are indispensable to fully exploit the discovery potential provided by the collected and forthcoming collider data.
The electroweak (EW) vector boson production, through the Drell-Yan (DY) mechanism [1,2], is the most "classical" hard-scattering process in hadronic collisions. The large production rates and clear experimental signatures make this processes important for detector calibration, as luminosity monitor and to probe underlying event. Moreover it plays a fundamental role in the contest of SM precision studies [3, 4, 5, 6, 7] and for searches of physics signals beyond the SM [8,9,10,11]. It is thus essential to provide accurate theoretical predictions, through detailed computations of the higher-order radiative corrections in QCD and in the EW theory, for vector boson production cross sections and related kinematical distributions. Among the various kinematical distributions the vector boson transverse-momentum (q T ) spectrum plays a special role. Precise knowledge of the Z boson q T distribution gives important information on the W boson spectrum which in turn directly affects the measurement of the W boson mass [12,13,14].
The next-to-next-to-leading order (NNLO) corrections in the QCD coupling α S have been computed for the total cross section [15,16], the rapidity distribution [17] and at fully differential level including the leptonic decay of the vector boson [18,19,20,21]. More recently next-to-nextto-next-to-leading order (N 3 LO) QCD calculations of the total cross section have been performed in Refs. [22,23]. The next-to-leading order (NLO) EW corrections, the mixed QCD-EW and QCD-QED corrections have also been computed [24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43]. In the large-q T region, where q T is of the order of the invariant mass of the lepton pair M , the fixed-order QCD corrections for the q T distribution are known up to O(α 2 S ) in analytic form [44,45,46,47,48] and up to O(α 3 S ) numerically through the fully exclusive NNLO calculation of vector boson production in association with jets [49,50,51,52,53]. However the bulk of the vector boson cross section lies in the small-q T region (q T M ) where the reliability of the fixed-order expansion is spoiled by the presence of large logarithmic corrections of the type ln(M 2 /q 2 T ) due to the initial-state radiation of soft and/or collinear partons. In order to obtain reliable perturbative QCD predictions, the enhanced-logarithmic terms have to be evaluated and systematically resummed to all orders in perturbation theory [54,55,56,57,58,59]. Resummed calculations at different levels of theoretical accuracy have been performed in Refs. [60,61,62,63,64,65,66,67,68] also applying methods from Soft Collinear Effective Theory [69,70,71,72,73,74,75] and transverse-momentum dependent factorisation [76,77,78,79,80,81].
In this Letter we apply the QCD transverse-momentum resummation formalism of Refs. [57,61,64] for the case of Z/γ * boson production up to N 3 LL accuracy. We analytically include all the N 3 LO terms at small-q T reaching full N 3 LL+N 3 LO accuracy in the small-q T region * . We implement our resummed calculation in the public numerical program DYTurbo [82] which provides fast and numerically precise predictions both for resummed and fixed-order QCD calculations including the full kinematical dependence of the decaying lepton pair with the corresponding spin correlations and the finite value of the Z boson width. We consistently match our resummed predictions with the NNLO numerical results at large-q T calculated in Refs. [50,53] and reported in Ref. [67] thus including the O(α 3 S ) corrections for the entire spectrum of q T . By using the connection between the q T resummation and the q T subtraction formalism [83] for fixed-order calculations we analytically [84] expanded the resummed results thus providing predictions for fiducial cross section of the Drell-Yan process both at N 3 LL+N 3 LO and at N 3 LO which, to our knowledge, have never appeared in the literature. Higher-order calculations beyond NLO QCD are definitely an hard task and are based on forefront and highly-specialized computations and numerical codes. The calculation presented in this paper is released as public software [85] with the aim of facilitating an efficient and wide spread of the results to the theoretical and experimental communities.
We briefly review the resummation formalism developed in Refs. [57,61,64] highlighting the main aspect relevant for our calculation. We consider the process where V denotes the vector boson † produced by the colliding hadrons h 1 and h 2 with a centreof-mass energy s, while l 3 and l 4 are the final state leptons produced by the V decay.
The hadronic cross section, fully differential in the lepton kinematics, is completely specified in terms of the transverse-momentum q T (with q T = q T 2 ), the rapidity y and the invariant mass M of the lepton pair, and by two additional variables Ω that specify the angular distribution of the leptons with respect to the vector boson momentum. 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,ŝ = 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, µ R and µ F are the renormalization and factorization scales. The last factor in the right-hand side of Eq. (2) is multi-differential partonic cross sections, computable in perturbative QCD as a series expansion in the strong coupling α S = α S (µ R ), which will be denoted in the following by the shorthand notation [dσ a 1 a 2 →l 3 l 4 ].
The partonic cross section can be decomposed as where the first term on the right-hand side of Eq. (3) is the resummed component which contains all the logarithmically-enhanced contributions of the type α n S M 2 /q 2 T ln m (M 2 /q 2 T ) (with 0 ≤ m ≤ 2n − 1) that have to be resummed to all orders, while the second term is the finite component which can be evaluated at fixed order in perturbation theory.
We perform the resummation in the impact-parameter space b [55]. The resummed component can then be written as dσ (res.) where J 0 (x) is the 0th-order Bessel function, the factor dσ is the Born level differential cross section for the partonic subprocess qq → V → l 3 l 4 .
The function W V (b, M,ŷ,ŝ) can be expressed in an exponential form by considering the 'double' (N 1 , N 2 ) Mellin moments with respect to the variables z 1 = e +ŷ M/ √ŝ and where we have introduced the logarithmic expansion parameterL ≡ ln( . is the Euler number). The scale Q ∼ M is the resummation scale [87], which parameterizes the arbitrariness in the resummation procedure.
The process dependent function H V [88,89] includes the hard-collinear contributions and it can be expanded in powers of α S as The universal (process independent) form factor exp{G} in the right-hand side of Eq. (5) contains all the terms that order-by-order in α S are logarithmically divergent as b → ∞ (i.e. q T → 0). The resummed logarithmic expansion of G reads where the functions g (n) control and resum the α k SL k (with k ≥ 1) logarithmic terms in the exponent of Eq. (5) due to soft and collinear radiation. At NLL+NLO we include the functions g (1) , g (2) and H (1) V , at NNLL+NNLO we also include the functions g (3) and H (2) V [90,91]. In order to reach full N 3 LL+N 3 LO accuracy in the small-q T region (i.e. including all the O(α 3 S ) terms) we have included the functions g (4) [92,93,94] and H V has been determined by exploiting its relation with the matching coefficients of the transverse-momentum dependent parton densities (TMD) calculated in Refs. [95,96] (see also Refs. [97,98,99]). The Mellin moments of the function H V have been calculated using the method of Ref. [100], and the FORM [101] packages summer [102] and harmpol [103]. The evolution of parton densities in Mellin space, and the Mellin moments of the splitting functions are calculated with the package QCD-PEGASUS [104], the Mellin inversion and the Fourier-Bessel inverse transform from the impact-parameter space are performed numerically as discussed in Ref. [82].
The function G is singular when α SL = π/β 0 (where β 0 is the one-loop coefficient of the QCD β function) which corresponds to the region of transverse-momenta of the order of the scale of ‡ For the sake of simplicity the explicit dependence on parton indices (which are relevant for the exponentiation in the multiflavour space) and the Mellin indices are understood. The interested reader can find the details in Ref. [57] (in particular Appendix A) and Ref. [86].
the Landau pole of the QCD coupling or b −1 ∼ Λ QCD . This signals that a truly non-perturbative (NP) region is approached and perturbative results (including resummed ones) are not reliable. In this region a model for NP QCD effects, which has to include a regularization of the singularity of the function G, is necessary. In our calculation we explicitly implemented the so-called Minimal Prescription [105,106,107] which has the advantage to regularize the Landau singularity in resummed calculations without introducing higher-twist power-suppressed contributions of the type O(Λ QCD /Q). Power-suppressed contributions can certainly be relevant at very small transverse-momentum (q T ∼ Λ QCD ) and should be eventually included, taking into account the delicate interplay with the leading-twist term, in order to correctly describe the experimental data in that region. In this Letter we have included the NP contribution in the form of a NP form factor S N P = exp{−g N P b 2 } with g N P = 0.6 GeV 2 § which multiplies the perturbative form factor exp{G(α S ,L)}, leaving a more detailed analysis of the inclusion of the power-suppressed NP contribution to future work.
Finally the finite component dσ .
We have performed the analytic expansion of the resummed component Eq. (4) up to O(α 3 S ) while the fixed-order cross section at large q T (formally at q T > 0) can be obtained from the the fullyexclusive computation of vector boson production in association with a jet at LO, NLO [109] and NNLO [49,50,51,52,53]. We observe that both the fixed-order cross section and the expansion of the resummed part are separately divergent with the same small-q T limit and the finite component formally satisfies the equation [82] lim We have checked that our analytic expression for the expansion of the resummed part agrees in the small-q T limit with the NNLO fixed-order results reported in Ref. [67] at permille level down to q T ∼ 4 GeV ¶ .
In the following we consider Z/γ * production and leptonic decay at the LHC. We present resummed predictions at NLL+NLO, NNLL+NNLO and N 3 LL+N 3 LO accuracy, matching our computation with the fixed-order results at large-q T respectively at LO, NLO and NNLO. The hadronic cross section is obtained by convoluting the partonic cross section in Eq. (3) with the parton densities functions (PDFs) from the NNPDF3.1 set [110] at NNLO with α S (m 2 Z ) = 0.118 where we have evaluated α S (µ 2 R ) at (n + 1)-loop order at N n LL+N n LO accuracy. In the case of Z production, because of the axial coupling, additional Feynman diagrams with quark loops contribute to the cross-section at O(α 2 S ) and O(α 3 S ). Their contribution cancels out for each isospin multiplet when massless quarks are considered. The effect of a finite top-quark mass in the third generation has been considered and found extremely small at O(α 2 S ) [111,46] while the finite mass top-quark contribution at O(α 3 S ) remains to be derived [112]. Therefore these contributions have currently been neglected in our calculation. We use the so called G µ scheme for EW couplings § This value is of the same order of the ones typically fitted in the literature, see e.g. Refs. [107,108,63].  with input parameters G F = 1.1663787 × 10 −5 GeV −2 , m Z = 91.1876 GeV, Γ Z = 2.4952 GeV, m W = 80.379 GeV. Our calculation implements the leptonic decays Z/γ * → l + l − and we include the effects of the Z/γ * interference and of the finite width Γ Z of the Z boson 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. The resummed calculation at fixed lepton momenta requires a q Trecoil procedure. We implement the general procedure described in Ref. [64] which is equivalent to compute the Born level distribution dσ (0) of Eq. (4) in the Collins-Soper rest frame [113].
We have applied the resummation formalism to the production of l + l − pairs from Z/γ * decay at the LHC ( √ s = 13 TeV) with the following fiducial cuts: the leptons are required to have transverse momentum p T > 25 GeV, pseudo-rapidity |η| < 2.5 while the lepton pair system, is required to have invariant mass 66 < M < 116 GeV and transverse momentum q T < 100 GeV ‖ .
In Fig. 1 we show the resummed component (see Eq. (9)) of the transverse-momentum distribution in the small-q T region. In order to estimate the size of yet uncalculated higher-order terms and the ensuing perturbative uncertainties we present the dependence of the resummed component on the auxiliary scales µ F , µ R and Q. The scale dependence band is obtained through independent variations of µ F , µ R and Q in the range M/2 ≤ {µ F , µ R , 2Q} ≤ 2M with the constraints 0.5 ≤ {µ F /µ R , Q/µ R , Q/µ F } ≤ 2 * * . The lower panel shows the ratio of the distribution with ‖ In order to match with the NNLO numerical results at large-q T we follow the kinematical selection cuts applied in Ref. [67]. * * In order to estimate the Q scale dependence of the resummed component we set the logarithmic expansion parameter to be L = ln(Q 2 b 2 /b 2 0 ) which is equivalent toL in the small-q T region.  respect to the N 3 LL+N 3 LO prediction at the central value of the scales µ F = µ R = 2Q = M . We observe that the NLL+NLO and NNLL+NNLO scale dependence bands do not overlap thus showing that the NLL+NLO scale variation underestimates the true perturbative uncertainty. This feature was observed and discussed in Ref. [64]. Conversely the NNLL+NNLO and N 3 LL+N 3 LO scale variation bands do overlap in the entire region q T < 30 GeV (except that they nearly overlap in the window 1 < q T < 4 GeV) thus suggesting that, from NNLL+NLO, missing higher order corrections are correctly estimated by scale variations. We also observe that the scale dependence is reduced by a factor of 2 (or more) going from NNLL+NNLO to N 3 LL+N 3 LO: the scale variation at N 3 LL+N 3 LO (NNLL+NNLO) is around ±0.8% (±2.5%) at the peak (q T ∼ 4 GeV), then it reduces at ±0.3% (±0.8%) level at q T ∼ 12 GeV and increase up to ±0.4% (±1.4%) level at q T ∼ 25 GeV. Finally, we note that in the low q T region non-perturbative effects are expected to become important. In particular by considering variations of the NP parameter in the range 0.3 ≤ g N P ≤ 0.9 GeV 2 we obtain the following additional uncertainties for the N 3 LL+N 3 LO resummed prediction in Fig. 1: ±0.3% at q T ∼ 25 GeV, ±0.6% at q T ∼ 12 GeV and ±0.7% at q T ∼ 4 GeV. For q T 4 GeV the NP uncertainties rapidly increase at few percent level. These NP uncertainties have a delicate interplay with the non-perturbative parton densities uncertainties which deserve a careful analysis.
In Fig. 2   lower panel shows the K-factors K N n LO defined as the ratio of between the N n LL+N n LO and the N n−1 LL+N n−1 LO predictions (with n = 2, 3). By looking at the K-factors we observe that the impact of the N 3 LL+N 3 LO (NNLL+NNLO) corrections with respect to the previous order is around −4% (−19%) at the peak, then it becomes −0.1% (+22%) at q T ∼ 30 GeV and increase to +3% (+55%) at q T ∼ 90 GeV.
773.7 ± 1 759.8 ± 2 749.6 ±2.5 Table 1: Fiducial cross sections at the LHC ( √ s = 13 TeV): fixed-order results and corresponding resummation results obtained with the DYTurbo numerical program. The uncertainties on the values of the cross sections plots refer to an estimate of the numerical uncertainties in the integration.
By exploiting the connection between the q T resummation and the q T subtraction formalism [83] we are able to provide fixed-order results for fiducial cross sections up to N 3 LO ‡ ‡ . In Table 1 we report the predictions for the cross section in the fiducial region at NLO, NNLO and N 3 LO, ‡ ‡ A fully consistent N 3 LO calculation for hadronic cross sections would require PDFs at the corresponding order which are currently not available. Uncertainties from missing higher order PDFs have been studied in Refs. [114,115,116]. NLL+NLO, NNLL+NNLO, N 3 LL+N 3 LO fixing the auxiliary scales to their central values. The N 3 LO corrections decrease the NNLO cross-section at −1.5% level and the resummation effects further enhance the N 3 LO result by +0.5%. We observe that the K-factor between the N 3 LO and NNLO results is 0.985 which is comparable with results reported in Table I of Refs. [22,23].
Generally speaking scale dependence cannot be regarded as a consistent estimate of the perturbative uncertainty because the effect due to uncalculated higher-order terms is typically larger than conventional scale dependence (see e.g. the comments on the scale variation bands of Fig. 1). A more realistic uncertainty estimate of the "true" perturbative uncertainty can be obtained, for instance, by comparing two subsequent orders of the expansion at central values of the scales and using half of the difference between them to assign the perturbative uncertainty [117]. This procedure leads (see Table 1) to an uncertainty of about ±2% (±4%) at NLO (NLL+NLO), ±0.7% (±1%) at NNLO (NNLL+NNLO) and ±0.8% (±0.7%) at N 3 LO (N 3 LL+N 3 LO). As expected (see comments below on fiducial cuts) this procedure shows a better convergence on the uncertainties of the resummed perturbative expansion.
In order to judge the numerical stability of the matching procedure and the effects of the power-suppressed terms we consider the contrubution of the finite part of the cross section. We show in Fig. 3 the finite part of the cross section for central values of the scales at O(α S ), O(α 2 S ) and O(α 3 S ). The effect of the finite component smoothly vanish as q T → 0 (see Eq. 9) and gives a small contribution to the matched result in the small q T region: the integral over the ranges 4 < q T < 20 GeV and 1 < q T < 4 GeV of the LO finite component represents respectively the 1.5% and 0.12% of the NLL+NLO fiducial cross section in Table 1, the O(α 2 S ) correction in the same ranges is respectively the 0.10% and −0.04% of the NNLL+NNLO result while the O(α 3 S ) correction in the range 4 < q T < 20 GeV the 0.16% of the N 3 LL+N 3 LO. As previously observed below q T ∼ 4 GeV the agreement with the O(α 3 S ) results of Ref. [67] (see Fig. 2) deteriorates. However the finite component gives a tiny contribution in the small q T region and we thus avoided to include in our results the finite part at O(α 3 S ) for q T < 4 GeV.
The results in Table 1 have been obtained applying the symmetric lepton p T cuts previously defined. It is well known [118,117] that in the case of symmetric cuts fixed-order calculations are affected by perturbative (soft-gluon) instabilities at higher orders. The results in Table 1 are obtained with a lower integration limit for the finite part of the cross section fixed to q T cut = 0.5 GeV and the numerical uncertainties include an estimate of the corresponding systematic uncertainty. More accurate fixed-order results and an estimate of such uncertainty can be obtained by evaluating the q T cut → 0 extrapolation or by a direct calculation of perturbative power corrections of the type O((q T cut /M ) p ) with p > 0 which are neglected for q T cut > 0 [119,120,121,122,73,123]. We stress however that the inclusion of such contributions cannot improve the physical predictivity of the fixed-order results in case of symmetric cuts which are affected by sizable theoretical instabilities produced by the soft-gluon effects.
We have performed the implementation of both the q T resummation and q T subtraction formalism for Drell-Yan processes up to N 3 LL+N 3 LO and N 3 LO in the DYTurbo numerical program. In this Letter we have illustrated the first numerical results for the case of Z/γ * production and leptonic decay at the LHC. [4] ATLAS Collaboration, "Measurement of the forward-backward asymmetry of electron and muon pair-production in pp collisions at √ s = 7 TeV with the ATLAS detector," JHEP 09 (2015) 049, arXiv:1503.03709 [hep-ex].