DYTurbo: fast predictions for Drell–Yan processes

Drell–Yan lepton pair production processes are extremely important for standard model (SM) precision tests and for beyond the SM searches at hadron colliders. Fast and accurate predictions are essential to enable the best use of the precision measurements of these processes; they are used for parton density fits, for the extraction of fundamental parameters of the SM, and for the estimation of background processes in searches. This paper describes a new numerical program, DYTurbo, for the calculation of the QCD transverse-momentum resummation of Drell–Yan cross sections up to next-to-next-to-leading logarithmic accuracy combined with the fixed-order results at next-to-next-to-leading order (O(αS2)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {O}(\alpha _{\mathrm {S}}^2)$$\end{document}), including the full kinematical dependence of the decaying lepton pair with the corresponding spin correlations and the finite-width effects. The DYTurbo program is an improved reimplementation of the DYqT, DYRes and DYNNLO programs, which provides fast and numerically precise predictions through the factorisation of the cross section into production and decay variables, and the usage of quadrature rules based on interpolating functions for the integration over kinematic variables.


Introduction
The Drell-Yan process denotes massive lepton-pair production in hadron-hadron collisions at high energies, as a e-mail: stefano.camarda@cern.ch (corresponding author) proposed by Sidney D. Drell and Tung-Mow Yan in 1970 [1], and first observed at the Alternating Gradient Synchrotron [2]. At the Large Hadron Collider (LHC) [3], the Drell-Yan process continues to play a fundamental role in probing the proton parton distribution functions (PDF), thereby providing valuable information on the u-and d-quark valence PDFs [4] and insight into the light-quark sea decomposition, in particular on the s-overd-quark ratio [5]. This process is also used to measure fundamental electroweak parameters such as the mass of the W boson [6], the weakmixing angle [7,8], and the W -boson width [9]. An accurate modelling of the Drell-Yan process is of paramount importance for searches of new physics phenomena beyond the Standard Model (SM) in final states with high dilepton invariant mass [10][11][12][13]. These experimental measurements need to be compared to accurate predictions based on high-order perturbative QCD and electroweak corrections. The Drell-Yan production total cross section and the vector boson rapidity distribution have been analytically computed up to the next-to-next-to-leading order (NNLO) in powers of the QCD coupling α S in Refs. [14,15] and [16], respectively. Fully exclusive parton-level NNLO calculations, which include the leptonic decay of the vector boson, have been implemented in publicly available Monte Carlo codes [17][18][19][20]. The transverse-momentum (q T ) distribution of the lepton pair at large (formally, non-vanishing) values of q T can be evaluated at O(α 3 S ) from the parton-level calculations of W /Z /γ * + jet production that have been performed in Refs. [21][22][23][24][25]. Various calculations that combine the QCD resummation formal-ism of logarithmically enhanced contributions at small-q T [26][27][28][29] with fixed-order perturbative results at different levels of theoretical accuracy have been performed in Refs. [30][31][32][33][34][35][36][37][38]. Analogous resummed calculations have been performed by applying Soft Collinear Effective Theory methods [39][40][41][42][43] and transverse-momentum dependent factorisation [44][45][46][47][48][49][50][51][52]. Electroweak (EW) [53][54][55][56][57][58][59][60][61] and mixed QCD-EW [62][63][64][65][66][67][68] radiative corrections have also been considered. A reliable estimate of the theoretical uncertainties requires various procedures, which also include variations of PDFs, renormalisation and factorisation scales, and SM parameters. It is thus necessary to rely on computing codes that allow fast calculations of these variations with small numerical uncertainties. The DYTurbo program, which is presented in this paper, aims at providing fast and numerically precise predictions of the Drell-Yan production cross sections, for phenomenological applications such as QCD analyses and extraction of fundamental parameters of the SM. The enhancement in performance over original programs is achieved by overhauling pre-existing codes, by factorising the differential cross section into production and decay variables, and by introducing the usage of one-dimensional and multi-dimensional numerical integration based on interpolating functions. The DYTurbo program is a reimplementation of the DYRes [36] and DYqT [33] programs for q T resummation, and of the DYNNLO [19] program for the finite-order perturbative QCD calculation up to NNLO. The DYRes [36] and DYqT [33] programs encode the q T resummed cross sections up to nextto-next-to-leading-logarithmic (NNLL) accuracy by using the resummation formalism proposed in Refs. [69][70][71]. The W +jet and Z /γ * +jet predictions at O(α S ) and O(α 2 S ) are reimplemented from the analytical calculations of Refs. [72][73][74], as encoded in DYqT, for the case of the triple-differential production cross sections as a function of rapidity y, invariant mass m, and transverse momentum q T of the lepton pair, and from the MCFM program [75], as encoded in DYRes and DYNNLO, for the full kinematical dependence of the decaying leptons. Software profiling was employed to achieve code optimisation. The most successful optimisation strategies leading to significant performance improvement were hoisting loop-invariant expressions out of loops, removing conditional statements from loops to allow the compiler performing automatic loop vectorisation, and manual loop unrolling. The DYTurbo software is based on a modular C++ structure, with a few Fortran functions wrapped and interfaced to C++. Multi-threading is implemented with OpenMP, and through the Cuba library by means of fork/wait system calls [76]. A flexible user interface allows setting the parameters of the calculation through input files and command line options. The results are provided in the form of text files and ROOT histograms [77]. Preliminary versions of the DYTurbo program were used by the ATLAS Collaboration in Refs. [6,78,79]. The DYTurbo program is publicly available [80].

Predictions with DYTurbo
The DYTurbo program provides predictions for W and Z /γ * -boson (collectively denoted as V -boson) production cross sections, fully differential in the four momenta of the decay leptons, and inclusive over final-state QCD radiation. The cross sections can be computed by performing the resummation of logarithmically-enhanced contributions in the small-q T region of the leptons pairs at leading-logarithmic (LL), next-to-leading-logarithmic (NLL), and NNLL accuracy, and also including the corresponding finite-order QCD contributions at next-to-leading order (NLO) and NNLO. The logarithmically-enhanced terms are resummed by using the resummation formalism of Ref. [70] in impact-parameter space. The structure of the cross section calculations is summarised in Eqs. (1) and (4), and we refer the reader to the discussion in Refs. [19,33,36] for details on the theoretical formulation. Upon integration of final-state QCD radiation, the fully-differential Drell-Yan cross section is described by six kinematic variables corresponding to the momenta of the two leptons. To the purpose of reducing the complexity of the calculation, it is useful to reorganise the fully-differential Drell-Yan cross section by factorising the dynamics of the boson production, and the kinematics of the boson decay. The cross section is therefore expressed as a function of the transverse momentum q T , the rapidity y and the invariant mass m of the lepton pair, and three angular variables corresponding to the polar angle θ and azimuth φ of the lepton decay in a given boson rest frame and to the azimuth φ V of the boson in the laboratory frame. However, the cross section does not depend on φ V , since in unpolarised hadron collisions the initial-state hadrons, i.e. the incoming beams, are to very good approximation azimuthally symmetric. Therefore the dependence of the cross section on φ V is not considered further. In the following a distinction will be made between fiducial cross sections, where kinematic requirements are applied on the final state leptons, and total or full-lepton phase space cross sections. The former requires the evaluation of the fivefold differential cross sections, the latter are (q T ,m,y)-dependent triple-differential cross sections integrated over cos θ and φ . At NLL+NLO and NNLL+NNLO, the q T -resummed cross section for V -boson production can be written as where dσ res is the resummed component of the cross-section, dσ asy is the asymptotic term that represents the fixed-order expansion of dσ res , and dσ f.o. is the V +jet finite-order cross section integrated over final-state QCD radiation. All the cross sections are differential in q 2 T . The resummed component dσ res is the most important term at small q T . The finite-order term dσ f.o. gives the larger net contribution at large q T . The fixed-order expansion of the resummed component dσ asy embodies the singular behaviour of the finite-order term, providing a smooth behaviour of Eq. (1) as q T approaches zero. The two finite-order terms of Eq. (1) and the finite-order factor H V (N)NLO in dσ res (see Eq. (2)) are calculated up to the same power in α S . The resummed component and its fixed-order expansion are given by 1 where Q denotes the auxiliary resummation scale [70] that is introduced in dσ res and, consistently, in dσ asy . The term dσ V LO (q T ) is the leading-order (LO) cross section evaluated for non-vanishing values of q T according to a given q Trecoil prescription [36], namely, with values of θ and φ that correspond to a chosen dilepton rest frame. The factor H V is the hard-collinear coefficient function. The term G is the exponent of the Sudakov form factor and it is originally expressed as a function of the impact parameter b, which is the Fourier-conjugate variable to q T . This term embodies the resummation of the logarithmically-enhanced contributions at LL, NLL or NNLL accuracy in b space. In order to parameterise non-perturbative QCD effects, the Sudakov form factor includes a non-perturbative contribution, whose simplest form is a Gaussian form factor. The b space expression of the Sudakov form factor is then evaluated in q T space by numerically performing the (inverse) Fourier transformation. The function Σ V (q T /Q) arises from the finite-order expansion of H V × exp{G}, and it matches the singular behaviour of dσ f.o. in the region q T → 0. An additional feature of the DYTurbo program is the possibility of computing finiteorder cross sections at LO, NLO and NNLO without the resummation of logarithmically-enhanced contributions. At NLO and NNLO, the finite-order cross section for V -boson production is computed by using the q T -subtraction formalism [81], and it is expressed as the sum of three components: with dσ CT (N)LO given by The LO cross-section term dσ V LO = dσ V LO (q T )δ(q 2 T ) is evaluated at q T = 0, and dσ V+jet is the V +jet cross section. 2 A uni- 1 The convolution with PDFs and the sum over different initial-state partonic contributions are implied in the shorthand notation of Eqs. (2), (3) and (5). Analogously, the inverse Fourier transformation from b space to q T space is implied in Eq. (2). 2 More precisely the term dσ V+jet (N)LO in Eq. (4) has to be evaluated with q T > q Tcut , the lower integration limit on q T in Eq. (5) has to be understood to be q Tcut and the square bracket term in Eq. (4) has to be evaluated in the limit q Tcut → 0. tarity constraint is implemented in the resummation formalism [70] so as to recover exactly the finite-order result upon integration over q T of the full-lepton phase space resummed cross section. The unitarity constraint leads to the following relation: The terms dσ res (N)NLL and dσ asy (N)LO can be, in general, multiplied by a switching function w(q T , m) above a given q T threshold, to the purpose of reducing the contribution of the resummed calculation in the large-q T region, where small-q T resummation cannot improve the accuracy of the finite-order calculation. The switching function can spoil the unitarity constraint of Eq. (6) by an amount which is smaller when the chosen q T threshold is larger. The default choice in DYTurbo is a Gaussian switching function, as used in DYRes. The Drell-Yan cross section predictions are obtained by integrating over the kinematic variables of the two leptons, and over additional variables related to QCD radiation, convolutions and integral transforms, as described in the following Sections. The integral transformations are evaluated by means of one-dimensional quadrature rules based on interpolating functions. The numerical integration over the other variables is performed with two different methods. The first method is based on the Vegas algorithm [82] as implemented in the Cuba library [83]. The second method employs a combination of one-dimensional and multi-dimensional numerical integrations based on interpolating functions. The one-dimensional integrations are performed by means of Gauss-Legendre quadrature rules, with nodes and weights evaluated with the Elhay-Kautsky method [84,85]. The multi-dimensional integrations are evaluated with the Cuhre algorithm [86,87] as implemented in the Cuba library [83] and in the Cubature package [88], and with a tensor product of Clenshaw-Curtis quadrature rules as implemented in the Cubature package. The Vegas integration method is available for all terms of the resummed and fixed-order calculations, and allows evaluating predictions for any arbitrary observable, for total and fiducial cross sections. The numerical integration based on interpolating functions is available for all the terms in the case of total cross sections, and for all the terms except the finite-order term at O(α 2 S ) in the case of fiducial cross sections. This integration method allows calculating only the cross sections as functions of q T , m, and y. Of these two methods, the former is the most versatile, whereas the latter allows reaching relative uncertainties in the predicted cross sections well below 10 −3 in a time frame that is significantly shorter than that required by the DYNNLO and DYRes programs. The EW parameters G F , α(m Z ), m W , m Z and sin 2 θ W of the Drell-Yan LO cross section are set by choosing three parameters as input, and calculating the others according to tree-level rela- at NNLL + NNLO accuracy for the Z -boson total cross section at √ s = 8 TeV as a function of the boson transverse momentum and the separated contribution of various terms: resummed component, asymptotic term, finite-order term, sum of asymptotic and finite-order terms tions. In the following the G μ scheme is used, in which 1876 GeV, and sin 2 θ W and α are calculated at tree level. The default values of the renormalisation (μ R ), factorisation (μ F ) and resummation scales are fixed to μ R = μ F = 2Q = m. The prescriptions necessary to obtain the resummed results (i.e. the q T -recoil prescription, the switching function w(q T , m) and the prescription to avoid the Landau singularity) have been chosen following Ref. [36]. Figure 1a shows results for Zboson production in proton-proton collisions at √ s = 8 TeV with the CT10nnlo set of parton density functions [89], and default choices of QCD scales and EW parameters. The relative contributions of the various terms to the Z -boson total cross section are illustrated in Figure 1b. The evaluation of each term is described in the following subsections.

Resummed component
The resummed component of the q T -resummed cross section (see Eq. (2)) can be factorised as the product of the LO cross section dσ V LO and the term W = H V × exp{G}. In these two terms, only the LO cross section depends on the lepton angular variables, and their integration is factorised as follows. The dependence of the cross section on cos θ is dσ (cos θ ) ∝ (1 + cos 2 θ ) + a cos θ , whereas an explicit dependence on φ enters only in the case of fiducial cross sections, due to the kinematic requirements on the final-state leptons. In the case of full-lepton phase space cross sections, the integration over the angular variables is obtained through the following substitutions: cos θ → d cos θ = 0, where d = d cos θ dφ . In the more general case of fiducial cross sections, the integrals in Eqs. (7) and (8) are as follows where K is the acceptance function of the kinematic requirements. The integrals in Eq. (9) are evaluated by first searching all the values of cos θ corresponding to the extremities of the region defined by the kinematic requirements at fixed values of φ . For each pair of cos θ extreme values, the integrals are evaluated analytically in d cos θ . In a second step, the integration along dφ is performed by means of Gauss-Legendre quadrature. In the case of the full-lepton phase space cross sections, the expressions in Eqs. (7) and (8) do not depend on q T and y, which allows further simplifications. In contrast, for the fiducial cross sections, the K acceptance function in Eq. (9), and so the integrals, depend in general on q T and y. Such a dependence varies between different q T -recoil prescriptions and it is of O(q T /m) at small q T . The W term is expressed through the Sudakov form factor exp{G} in b space. The q T -dependent cross section is obtained by means of a two-dimensional inverse Fourier transformation, which is expressed as a zeroth-order inverse Hankel transformation by exploiting the azimuthal symmetry of the W function in the transverse plane: (10) whereW is the expression of W in b space, J 0 (x) is the zeroth-order Bessel function and s is the centre-of-mass energy. The integral transformation of Eq. (10) is computed by means of a double-exponential formula for numerical integration [90][91][92]. The convolution with PDFs is more efficiently performed by considering double Mellin moments of the partonic functionsŴ ab , defined aŝ where where x 1,2 = m/ √ s e ±y , c is a real number which lies at the right of all the poles of the integrand, and F N i , with i = a, b, are Mellin moments of PDFs, f i (x), defined as: The integral transformation of Eq. (12) is computed by means of Gauss-Legendre quadrature, and the PDFs are evolved [33,36] from the factorisation scale μ F to the scale b 0 /b (b 0 = 2e −γ E , and γ E is the Euler number) by using the program Pegasus QCD for the evolution of PDFs in Mellin space [94]. To perform the Mellin inversion, it is necessary to calculate the Mellin moments F i (N ) at values of N along a contour of integration in the complex plane. Parameterising the PDFs in a simple form such as where α, β are constants and P(x) is a polynomial, Mellin moments for arbitrary complex N can be calculated through a simple formula involving the Γ function: Thanks to the analytic continuation of Eq. (15) in the region of the complex plane with (N ) < 0, when PDFs are expressed with this form, the integration contour in Eq. (12) can be optimised by bending towards negative values of (N ), as depicted schematically in Fig. 2, allowing for a faster convergence of the Mellin inversion integral. Such a strategy is adopted in DYRes and in Refs. [94,95]. As a drawback, PDFs need to be parameterised as in Eq. (14), or an approximation of the PDFs that follows this form has to be evaluated, which is significantly time consuming. In DYTurbo, the Mellin moments of PDFs are evaluated numerically, by using Gauss-Legendre quadrature to calculate the integrals of Eq. (13). However these integrals can be evaluated numerically only for (N ) > 0. As a consequence the integration contour of the inverse Mellin transform cannot be bent towards negative values of (N ), and a standard contour along the straight line [c − i∞, c + i∞] is used (see Fig. 2). This procedure results in a slower convergence of the integration in Eq. (12), for which about twice as many function evaluations are required, but it has the great advantage of allowing usage of PDFs with arbitrary parameterisation, without requiring knowledge of their functional form, and without requiring any time consuming evaluation of an approximation of PDFs in the form of Eq. (14). The integration over the V -boson rapidity, y, is factorised as follows. In the case of total cross sections, the values of the angular integrals in Eqs. (7) and (8) do not depend on y. The only dependence on the rapidity in Eq. (12) is in the expression and the integrals of Eq. (16) are evaluated analytically using the following relation: where y 0 and y 1 are the lower and upper y-bin boundaries. In the case N 1 = N 2 , Eq. (17) further simplifies to y 1 −y 0 . When the y-bin boundaries are larger than the allowed kinematic range |y| ≤ y max , with y max = ln( √ s/m), Eq. (17) simplifies to 2πi δ(N 1 −N 2 ), and the double Mellin inversion is reduced to a single Mellin inversion [71] by setting N 1 = N 2 = N . In the case of fiducial cross sections, the values of θ 0 and θ 1 in Eq. (9) also depend on y, and the integrals are evaluated numerically for all pairs of N 1 and N 2 by means of Gauss-Legendre quadrature. The integration over the Vboson transverse momentum, q T , can be performed analytically in the case of the full-lepton phase space cross sections, since the expressions in Eqs. (7) and (8) do not depend on q T , and the only term that depends on q T is J 0 (bq T ) of Eq. (10). By using the relation dx x J 0 (x) = x J 1 (x), the integration over q T in a bin of boundaries q 0 T and q 1 T can be evaluated as Similarly to Eq. (10), the integral of Eq. (19) is computed by means of a double-exponential formula for numerical integration, and by performing two separate integrations corresponding to the terms J 1 (bq 1 T ) and J 1 (bq 0 T ). The information of the one-loop (two-loop) virtual correction to the LO subprocess is contained in the H V function. In the computation of the fixed-order cross section of Eq. (4), the H V function is evaluated in x-space, i.e. without performing a Mellin transformation, and the convolution with PDFs is performed by integrating over the variables z 1,2 = e ±ŷ m/ √ŝ . The corresponding integrals are calculated with Gauss-Legendre quadrature.

Asymptotic term and counter-term
The asymptotic term of Eq. (3) and the counter-term of Eq. (5) are computed using the function Σ V (q T /Q), which embodies the singular behaviour of dσ f.o. in the limit q T → 0. In the finite-order case the counter-term contributes at q T = 0. Accordingly, the LO cross section is evaluated at q T = 0 and the function Σ V (q T /Q) is integrated over the auxiliary variable q T . At variance, in the resummed case the asymptotic term is a function of q T , and the LO cross section is evaluated for nonzero values of q T according to a given q Trecoil prescription. As for the resummed term, the integration over the angular variables is factorised in the LO cross section by using Eqs. (7) and (8) or Eq. (9). The function Σ V (q T /Q) is evaluated in x-space, i.e. without performing a Mellin transformation, and the convolution with PDFs is performed by integrating over the variables z 1,2 with Gauss-Legendre quadrature. In the case of full-lepton phase space cross sections, the q T dependence of the asymptotic term and of the function Σ V (q T /Q) is fully embodied in a set of four functionsĨ n (q T /Q) with n = 1, . . . , 4 [70]. The integration over q T of the asymptotic term is performed by integrating theĨ n (q T /Q) functions with Gauss-Legendre quadrature. In the case of fiducial cross sections, the values of θ 0 and θ 1 of Eq. (9) also depend on q T , and the integrals where q 0 T and q 1 T are the lower and upper q T -bin boundaries, are evaluated numerically by means of Gauss-Legendre quadrature.

Finite-order term
The real-emission corrections are embodied in the (N)LO finite-order term of Eq. (1) and in the V +jet term of Eq. (4) for the resummed and fixed-order predictions, respectively. Since DYTurbo provides results that are inclusive over finalstate QCD radiation, the two terms are fully equivalent. 3 Two independent calculations of this term are implemented. The first calculation, which is based on the code MCFM [75], is fully differential with respect to the lepton angular variables and the final-state QCD radiation. The second calculation, which is inclusive over the lepton angles and the QCD radiation, implements the analytic results of Refs. [72][73][74], and it relies in part on the code taken from DYqT [33]. The MCFM implementation of the lowest-order term dσ V+jet LO can be evaluated by using either the Vegas integration method or the numerical integration based on interpolating functions. The MCFM implementation of the next-order term dσ V+jet NLO is the most complex part of the calculation, and it can be evaluated only with the Vegas algorithm. The reason is that this NLO calculation is based on the Catani-Seymour dipole subtraction scheme [96], in which for each point in the phase space where the real radiation is evaluated, a set of counterterm dipoles are computed corresponding to various different phase-space points. As in any local subtraction procedure, the resulting integrand presents discontinuities and it cannot be efficiently approximated by interpolating functions. The implementation of the analytic calculation of Refs. [72][73][74] yields the triple-differential production cross sections as a

Tests of numerical precision
In order to validate the numerical precision of the resummed calculation, three closure tests are performed: the comparison of the fixed-order expansion of the resummed component (asymptotic term) and the finite-order term at small q T , the comparison of the term H V × dσ V LO and the resummed component upon q T integration, and comparisons of the integration methods available in DYTurbo, namely the Vegas algorithm and the multi-dimensional numerical integration based on interpolating functions, referred to as Quadrature integration in the plots. The numerical tests of this section are performed in full-lepton phase space, using the CT10nnlo set of parton density functions and with default values of the QCD scales and EW parameters. As discussed in Sect. 2, the function dσ asy embodies the singular behaviour of dσ f.o. when q T → 0, yielding the relation (see Eq. (4) of Ref. [70]) or, equivalently, We note that dσ f.o. and dσ asy in Eq. (22) separately diverge proportionally to q −2 T (modulo powers of log q T ) as q T → 0. Computing such a relation at small values of Q T provides a stringent test of the numerical precision of the asymptotic and finite-order terms. The triple-differential cross sections dσ asy and dσ f.o. , as functions of q T , m and y, are evaluated at the fixed values y = 0 and m = m V , with V = W, Z , for proton-proton collisions at √ s = 13 TeV. The result of the test is shown in Figure 3 for the NLL+NLO and NNLL+NNLO calculations. In all the cases, the relation of Eq. (22) is shown at values of q T as low as q T = 0.01 GeV. As a second closure test, the unitarity constraint of Eq. (6), which relates the H V × dσ V LO and dσ res terms, is tested. Computing such a relation provides a stringent test of the numerical precision of the procedure. The triple-differential cross sections dσ res as a function of q T , m, and y are integrated in q T from zero to infinity, and in the range of rapidity |y| ≤ y max , and they are compared with the H V × dσ V LO (0) double differential cross sections as a function of m and y, integrated in the same range of rapidity. The switching function w(q T , m), which reduces the contribution of the resummed calculation in the large-q T region, is not used in this test. Figure 4 shows the result of such a comparison at LO, NLO and NNLO for Z /γ * -boson production in proton-proton collisions at √ s = 13 TeV, for 180 equally-spaced bins of m in the range [20,200] GeV. In all cases the relation ∞ 0 dq 2 T dσ res / H V × dσ V LO = 1 is verified, with deviations from unity that are smaller than 10 −6 . The terms H V × dσ V LO and dσ res are evaluated in xspace and Mellin-space, respectively. Therefore, computing such a relation also provides a test of the numerical precision of the Mellin inverse transformation in Eq. (12). Similar level of agreement is observed by performing this closure test as a function of the rapidity. The results at NNLL+NNLO from DYTurbo using the numerical integration based on interpolating functions and the Vegas algorithm are compared for the Z -boson differential cross section in proton-proton collisions at √ s = 8 TeV. The Z -boson invariant mass range is required to be 80 GeV < m < 100 GeV. Ratios of these results are shown in Fig. 5. The scatter in the central values of the points is at the permille level, except for the points at the very edges of the kinematic phase space, and for the finite-order term at high q T where deviations of two permille are observed. From left to right: resummed component, asymptotic term, finite-order term, and total differential cross section

Benchmark results
This section provides benchmark results of DYTurbo to DYRes at NNLL+NNLO for cross sections differential in q T , and benchmark results of DYTurbo to DYRes and DYNNLO for fully-integrated fiducial cross sections at NNLL+NNLO and NNLO. Z -boson and positively-charged W -boson production cross sections. The Z -boson fiducial phase space is defined by the lepton transverse momentum p T > 20 GeV, the lepton pseudorapidity |η | < 2.4, and invariant mass of the lepton pair in the range 66 GeV < m < 116 GeV. All comparisons of predictions for the resummed term are validated at the better than 1% level while the comparison of the sum of the asymptotic and finite-order terms are validated at the ∼2% level for q T above 40 GeV. In particular, the positively-charged W -boson predictions show well that the sum of the asymptotic and finite-order terms converges to zero at low q T , as expected (it is also consistent with zero for the Z -boson predictions, within the Vegas uncertainties, which are highly correlated bin-to-bin). The DYTurbo and DYRes results are compared in Fig. 8 after summing all terms. Also in this case good agreement is observed between the two codes, within the numerical uncertainty of the Vegas integration.

Benchmark of fully-integrated cross-section results
Benchmark results for fully-integrated fiducial cross section at NNLL+NNLO from DYRes [36] and at NNLO from DYNNLO [19] are shown in Table 1 and compared with the corresponding results calculated with DYTurbo. 4 The predictions are evaluated for proton-proton collisions at the centre-of-mass energy √ s = 8 TeV, and according to the fiducial definition and QCD and EW settings of Ref. [99]. The Z -and W -boson fiducial phase space is defined by the charged lepton and neutrino transverse momentum p ,ν T > 25 GeV, the charged lepton pseudorapidity |η | < 2.5, and invariant mass of the lepton pair larger than 50 GeV for Zboson production and larger than 1 GeV for W -boson production. The results for DYNNLO shown in the table are taken from Table 12 of Ref. [99]. The results of DYTurbo are in agreement with the results of the other programs considered in Ref. [99]. Differences as large as 1% are observed between the NNLO and the NNLL + NNLO results, which are mostly due to recoil effects in the lepton kinematics (the unitarity constraint of Eq. (6) between fixed-order and resummed cal- 4 The NNLO results in Table 1 are obtained with a minimum value of r = q T /m fixed to r cut = 0.002 and their corresponding numerical uncertainties do not include the systematic uncertainty from the r cut → 0 extrapolation. A more accurate NNLO result and an estimate of such uncertainty can be obtained by evaluating the cross section at different values of r cut and carrying out the limit r cut → 0 [98]. culations does not apply in the presence of lepton kinematic cuts). An additional source of difference between the NNLO and the NNLL + NNLO results is the inclusion of the switching function w(q T , m), which affects the cross sections at the per mille level.

Time performance
In this section various tests of time performance are discussed. The computation time requested to calculate crosssection predictions for DYTurbo and DYRes is compared and used to assess the performance improvement of DYTurbo. The amount of time required to perform a calculation as a function of threads provides a test of the scaling behaviour of the multi-threading implementation. The timeperformance tests are run on a server mounting two AMD Opteron 6344 CPUs with 12 cores each. The fully-integrated fiducial cross section of Z -boson production, as defined for the results shown in Table 1, is computed with DYRes and DYTurbo at NNLL+NNLO. The DYRes calculation took 40 hours for an uncertainty of 0.4%, whereas DYTurbo took 8 hours, yielding a factor of 5 in the improvement of the time performance. Figure 9 shows the speedup factors for crosssection calculations as a function of the number of threads, where the speedup is defined as the ratio of elapsed time of the multi-threaded calculation divided by the reference (a) (b) Fig. 9 Computing time as a function of the number of threads for Z -boson production at √ s = 8 TeV with different target precision: a NNLO results with the Vegas integration method, and b NNLL + NNLO results with the Quadrature integration method elapsed time of the one-thread calculation. Assuming that the one-thread calculation has a parallelisable time fraction and a non-parallelisable time fraction, and the multithreading process has an overhead time proportional to the number of threads N , the speedup curve can be parameterised as where s is the ratio between non-parallelisable and parallelisable times, and o is the overhead time per thread. The measured speedup factors are well described by this model, with overhead times compatible with zero, and with fractions of non-parallelisable time which are smaller when the target precision is higher. Indeed most of the non-parallelisable time is spent in the program initialisation, which becomes negligible for long runs with high target precision. We conclude this section reporting typical running times for fast and numerically precise DYTurbo predictions with the numerical integration based on interpolating functions. Figure 10 shows NLL+NLO and NNLL+NNLO predictions for the Z -boson production at 13 TeV in fulllepton phase space, integrated in the range of invariant mass [66,116] GeV and in the range of rapidity |y| ≤ y max ∼ 5.3. The predictions are computed in 100 equally-spaced q T bins from zero to 25 GeV. The predictions are evaluated with a target in the relative numerical uncertainty of 10 −4 for each term, and using simultaneously 20 parallel threads. The computation time required to perform the full calculation is 4 min at NLL + NLO and 3.4 h at NNLL + NNLO. The computation of the resummed component required 6 s at NLL+NLO and 10 s at NNLL + NNLO, the computation of the asymptotic term required 0.2 s at NLL + NLO and 0.7 s at NNLL+NNLO, the computation of the finite-order term required 4 min at NLL + NLO and 3.4 h at NNLL+ NNLO. In these examples, as in all other time-performance tests, the great majority of the computation time is spent to evaluate the finite-order term. For applications as PDF fits, where very fast predictions are required, this part of the calculation could be computed by using APPLGRID [100].

Conclusions
The DYTurbo program provides fast and numerical precise predictions of Drell-Yan processes, through a new implementation of the DYqT, DYRes and DYNNLO numerical codes. The cross-section predictions include the calculation of the QCD transverse-momentum resummation up to nextto-next-to-leading logarithmic accuracy combined with the fixed-order results at next-to-next-to-leading order (O(α 2 S )). They also include the full kinematical dependence of the decaying lepton pair with the corresponding spin correlations and the finite-width effects. The enhancement in per-formance over previous programs is achieved by code optimisation, by factorising the cross section into production and decay variables, and with the usage of numerical integration based on interpolating functions. The resulting cross-section predictions are in agreement with the results of the original programs. The great reduction of computing time for performing cross-sections calculation opens new possibilities for the usage of Drell-Yan processes for PDF fits, for the extraction of fundamental parameters of the SM, such as the mass of the W boson and the weak-mixing angle, and for the estimation of background processes in searches for physics beyond the SM.

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