Precision studies for Drell-Yan processes at NNLO

We present a detailed comparison of the fixed-order predictions computed by four publicly available computer codes for Drell-Yan processes at the LHC and Tevatron colliders. We point out that while there is agreement among the predictions at the next-to-leading order accuracy, the predictions at the next-to-next-to-leading order (NNLO) differ, whose extent depends on the observable. The sizes of the differences in general are at least similar, sometimes larger than the sizes of the NNLO corrections themselves. We demonstrate that the neglected power corrections by the codes that use global slicing methods for the regularization of double real emissions can be the source of the differences. Depending on the fiducial cuts, those power corrections become linear, hence enhanced as compared to quadratic ones that are considered standard.


I. INTRODUCTION
The very high quality of data for many Standard Model (SM) scattering processes collected at the Large Hadron Collider in recent years makes it mandatory to use high precision theoretical predictions for physics analyses of these data. This is true especially for the vector-boson hadroproduction that is the prime benchmark process at hadron colliders. The theory computations need to account for radiative corrections, the dominating ones being due to quantum chromodynamics (QCD), where presently the next-to-next-to-leading order (NNLO) is the state-of-the-art [1]. In addition, the experimental cuts such as limits on the transverse momenta or the rapidities of the observed final state particles applied during the data selection have to be included in the theory calculations. Cross section predictions have to be provided for the respective fiducial regions.
For colorless final states such as the Drell-Yan process or more specifically, for W ± -and Zboson production including their decay, the QCD predictions are known to NNLO for fully exclusive kinematics [2][3][4][5]. These computations require the combination of squared matrix elements with three different multiplicities of partons in the final state and the subsequent cancellation of the soft and collinear singularities upon integration over their phase space to arrive at infrared finite results, a step which is performed with the help of dedicated subtraction schemes. In summary, the problem is considered being solved and the NNLO QCD predictions for W ± -and Z-boson production including the leptonic decay have been made available in several computer programs which perform the integration of the parton level predictions numerically by Monte Carlo methods. It has become apparent, however, that the level of accuracy of the published codes is not sufficient for the needs in analyses of experimental data. Significant deviations between the predictions of the different codes have been documented, for example, in an ATLAS analysis [6].
The motivation for the present study comes from the use of data on W ± -and Z-boson hadroproduction collected at the LHC and the Tevatron in the determination of parton distribution functions (PDFs) at NNLO accuracy. Those data consist of differential distributions in the pseudorapidity of the decay leptons and typically have a precision of O(1 − 2%), which is mainly dominated by the experimental systematics. This fact and the generally rather small size of the pure NNLO QCD corrections in the range of O(1%) or even less relative to the fiducial cross sections at next-to-leading order (NLO) lead us to investigate the precision of available QCD predictions. To that end, we focus on two aspects in this work. For one, we provide benchmark numbers for NNLO QCD predictions in kinematics which are representative for the bulk of the available experimental data. As a target we aim at predictions with a residual uncertainty of less than O(0.1 o / oo ) in each bin from the Monte Carlo integration for the cross section integrated over the fiducial region, whenever possible in order to have accurate results when comparing the different codes. Previous comparisons [7] of some of the published codes have limitations in the precision of the predictions. A second point concerns the investigation of the observed deviations among the codes in the light of the subtraction schemes used, which are either local or global, depending on whether the cancellations of the infrared singularities are performed locally in the integrand at each point in phase space or whether they are accomplished globally after integrating over a slice of the phase space. In particular, we illustrate the impact of fiducial cuts on the decay leptons for subtraction schemes with slicing.
The paper is organized as follows. In Sec. II we present the benchmark numbers for two representative sets of W ± -and Z-boson data from the LHC and the Tevatron. We first validate the different codes at NLO in QCD and then quantify the deviations at NNLO. In Sec. III we provide a brief review of global slicing methods and the power counting in the slicing parameter. Then we compute the effect and the size of power corrections on the example of the lepton decay phase space for the fiducial cuts of the LHC and Tevatron data considered in the previous section. We conclude in Sec. IV.

A. Set-up and validation
The set-up for benchmarking available QCD predictions for W ± -and Z-boson hadro-production cross sections up to NNLO in QCD in the fiducial phase space of the experimental measurement contains three aspects: the choice of the data sets, the list of input parameters and the selection of the NNLO QCD codes for the comparison of the theoretical predictions.
We choose two sets of data on W ± -and Z-boson production collected by the ATLAS experiment at the LHC and the DØ experiment at the Tevatron, respectively, which are statistically significant in current fits of PDFs.
• The ATLAS data set for the W ± -and Z/γ * -production cross sections [6] measured at a center-of-mass energy of √ s = 7 TeV at the LHC. These data are given in form of pseudorapidity distributions for the decay electron or muon (W ± -production) and the decay leptonpair (Z/γ * -production), respectively. The transverse momenta p l T and the pseudo-rapidities η l of the decay leptons are subject to fiducial cuts. The cross sections for Z/γ * -production are measured at central as well as at forward pseudo-rapidities.
• The data obtained by DØ on W ± -boson production at √ s = 1.96 TeV at the Tevatron [8] measures the electron charge asymmetry distributions and their dependence on the electron pseudo-rapidity. These data also probe forward kinematics. Also, the DØ data apply fiducial cuts, both symmetric as well as staggered, on the transverse momenta p l,ν T of the electron and the neutrino and on their pseudo-rapidities.
Another data set by ATLAS, the measurement of the muon charge asymmetry in W ± -boson production at √ s = 8 TeV at the LHC [9] has similar experimental precision as the chosen data set [6] collected at √ s = 7 TeV and also largely overlaps in kinematics. Similar considerations apply to data from CMS and LHCb, e.g., [10,11] Hence, we do not include these data in the benchmark comparison.
We use the G µ scheme with input values G F , M Z , M W and with sin 2 (θ W ) and the QED coupling α(M Z ) as output values. This scheme minimizes the impact of NLO electroweak corrections, see e.g. [12]. In detail, our SM input parameters are [13] and for the relevant CKM parameters The computations are performed in the MS factorization scheme with n f = 5 light flavors. Therefore we take the n f = 5 flavor PDFs of ABMP16 [14,15] as an input together with the value of the strong coupling, α (n f =5) s (M Z ) = 0.1147. These choices do not bias any of the predictions we discuss below. The renormalization and factorization scales µ R and µ F are set to µ R = µ F = M V , where M V is the mass of the gauge boson V.
We run the following publicly available codes designed for the computation of the fully differential NNLO QCD predictions for the lepton rapidity distributions.
The codes differ by the subtraction schemes used. FEWZ uses sector decomposition and employs a fully local subtraction scheme. DYNNLO and MATRIX both use q T -subtraction [2] at NNLO, which is a global phase space slicing method. MCFM uses N-jettiness subtraction [21,22] at NNLO, which is also a global phase space slicing method, see [23] for a recent review. Codes with global slicing do require a slicing parameter. These are r cut for MATRIX as a cut on q T and τ cut for MCFM as the jettiness slicing parameter.
The DYNNLO program is a legacy code, now superseded by MATRIX. It has been the first publicly available program containing the NNLO QCD predictions for fully exclusive kinematics and it is included in this list because of the ATLAS study [6] and its continued use in the analyses of experimental data. An improved reimplementation of the DYNNLO code is also part of program DYTurbo for fast predictions for Drell-Yan processes [24] 5 . DYTurbo includes the resummation of large logarithmic corrections, too. In the fixed-order mode, software profiling was employed to achieve code optimization, but it reproduces exactly the results of DYNNLO within numerical uncertainties due to the different integration method, hence its predictions are not included in our benchmark comparisons. Another code which we mention for reference is SHERPA-NNLO-FO [25] 6 . We do not include this code in the benchmark exercise and instead refer to the previous study [7].
We start with the validation up to NLO of the codes selected for the comparison. This serves a two-fold purpose. First, it provides a check on the input settings for the codes and second, it demonstrates the level of agreement for the benchmark results. We show first the results for the predictions for the W ± -and Z/γ * -production cross sections at √ s = 7 TeV corresponding to the kinematics of the ATLAS data set [6]. In all cases, here and below we use the ABMP16 PDFs at NNLO and the value α  FIG. 1: The NLO QCD cross sections for inclusive pp → W ± + X → l ± ν + X as function of pseudo-rapidity η l , computed with DYNNLO, MATRIX and MCFM relative to FEWZ and using the ABMP16 PDFs. Cuts of p l,ν T ≥ 25 GeV and M T ≥ 40 GeV for the transverse momenta and mass are applied as in the ATLAS data selection [6]. The error bars indicate the accuracy of the numerical integration, shown by the horizontal dashed lines for the FEWZ result. The MATRIX result is plotted in the center of each η l -bin, while DYNNLO and MCFM results are shifted slightly to the left and the right.  Fig. 1 for the NLO QCD cross sections for pp → Z/γ * + X → l + l − + X production at √ s = 7 TeV as function of the pseudo-rapidity η ll . Cuts of p T 1,2 ≥ 25 GeV, 66 ≥ M ll ≥ 116 GeV and |η l i | ≤ 2.5, i = 1, 2 for the lepton pseudo-rapidities are applied. subtraction schemes except for DYNNLO, which uses q T -subtraction also at NLO 7 . FEWZ applies sector decomposition while MATRIX and MCFM all use by default the dipole subtraction [26].
In Fig. 1 we plot the results for the W ± -production cross sections at √ s = 7 TeV at NLO corresponding to the kinematics of the ATLAS data set [6] and we find agreement at the level of O(1 o / oo ) between FEWZ, MATRIX and MCFM in the entire η l range, while for DYNNLO we find the values to be systematically enhanced by O(5 o / oo ) for both, W + -and W − -production. In Figs 3 we show the NLO cross sections for Z/γ * -production at √ s = 7 TeV in the ATLAS kinematics with the different selection cuts on the lepton pseudo-rapidities. In the case of both leptons at central pseudo-rapidities |η l i | ≤ 2.5 for i = 1, 2 in Fig. 2 we find agreement among all codes, except for a slight systematic off-set of the DYNNLO result by O(3 o / oo ) for η ll 1.0. In Fig. 3 we display the case when one lepton is required at central and the other one instead at forward pseudo-rapidity, |η l 1 | ≤ 2.5 ≤ |η l 2 | ≤ 4.9. We observe agreement at the level of O(1 − 2 o / oo ) between FEWZ, MATRIX and MCFM as shown in Fig. 3 on the right, while the DYNNLO results turn out to be larger by up to a few per cent in the first η ll bins. Finally, in Fig. 4 we plot the electron charge asymmetry distribution A e in W ± -boson production at √ s = 1.96 TeV for two choices of cuts applied in the selection of the DØ data [8]: on the left we have symmetric cuts, p e,ν T ≥ 25 GeV, and on the right staggered cuts, p e T ≥ 35 GeV and p ν T ≥ 25 GeV. We show the ratio of the DYNNLO, MATRIX and MCFM results with FEWZ. The level of agreement for A e is better than O(1%), except in regions, where NLO predictions are very small. For symmetric cuts in Fig. 4 on the left, this happens in the first bin, where the relative agreement deteriorates to a few per cent, and in the bin η e = 2.1, where A e vanishes and the corresponding ratios have not been plotted. For staggered cuts in Fig. 4 on the right, this is observed also in the first bin and in the region η e 2.5. The individual cross sections for W ± -production are computed at an accuracy of O(0.1 o / oo ), but in the asymmetry, one looses more than one order of magnitude in precision.
In summary, the validation shows that the predictions by MATRIX, MCFM and FEWZ are all in very good agreement at NLO. DYNNLO using q T -subtraction also at NLO delivers results, which are typically accurate up to a few per mill and deviate in particular for distributions with challenging kinematics like in Fig. 3, where agreement can only be reached at the level of a few per cent. Moreover, the deviations of the DYNNLO results from the ones by MATRIX, MCFM and FEWZ display a particular pattern as a function of the (di-)lepton pseudo-rapidities η l (η ll ), which will be addressed in Sec. III below in the light of power corrections in the slicing parameter.

B. NNLO benchmark predictions
We now proceed to the QCD predictions at NNLO accuracy, again starting with W ± -and Z/γ *production cross sections at √ s = 7 TeV as measured by ATLAS data set [6]. As a baseline for the comparison, we have computed the NNLO QCD predictions with the ABMP16 PDFs [14] and FEWZ as this is the only available code which implements a fully local subtraction scheme at NNLO. We emphasize that our conclusions do not depend on the choice for the default theoretical prediction. Note, that the ATLAS data have been released after completion of the ABMP16 PDFs and were not included in the fit. However, these data are in a good agreement with the predictions.
In Fig. 5 we show the pulls of data and the differences of NNLO predictions obtained from DYNNLO normalized to the predictions computed with the FEWZ code for the same distributions as studied in the previous section at NLO. The nominal relative accuracy of the numerical integration is in all cases at the level of a few units in 10 −4 , thus negligible in the plots. For the lepton-pseudorapidity (η l ) distribution in W + -production in Fig. 5, the DYNNLO results are below the FEWZ ones by up to O(1%) and by a few per mill for W − -production, in both cases over the whole range in η l . For the di-lepton pseudo-rapidity distribution in the central Z/γ * -production the DYNNLO predictions are below the FEWZ ones by several per mill up to η ll ≤ 1.5 and tend to agree better with the FEWZ ones for larger η ll . For forward Z/γ * -production instead, we see the DYNNLO results being below the FEWZ ones by up to O(1 − 2%) in the entire η ll range, with significant deviations up to O(7%) in first bins.
In Fig. 6 we show the same comparison, now with the NNLO predictions obtained with the MATRIX code. Our findings regarding the deviations from the FEWZ results are qualitatively similar, quantitatively slightly smaller to those from DYNNLO, with both codes being based on the q T -slicing method. All DYNNLO results have been obtained with the default minimum value r min cut = q min T /M V = 0.8% for the slicing cut on q T . The MATRIX code uses r min cut = 0.15% as the default and offers also the choice r min cut = 0.05% for the Z/γ * -production process. In detail we find MATRIX results being below the FEWZ ones by up to O(1%) for the η l -distribution in W + -production, and by a few per mill for W − -production. For the η ll -distribution in Z/γ * -production at central rapidities we see the MATRIX numbers with the default value r min cut = 0.15% below the FEWZ ones by a few per mill for η ll ≤ 1.0, above the FEWZ ones in the bins around η ll 1.5 and in agreement with FEWZ for larger rapidities. On the other hand, the MATRIX results with r min cut = 0.15% for forward Z/γ * -production are below the FEWZ ones by up to O(1 − 2%) in the entire η ll range and show a significant deviation of O(5%) in the first bin. In order to deal with the residual dependence on the slicing cut in q T MATRIX extrapolates the total rates for r min cut → 0 and suggests a uniform rescaling of each bin by the ratio σ extrapolated NNLO /σ r cut NNLO , see also Sec. III. In Fig. 6 this rescaling has not been applied. If done, it would lead to upward shifts of the central values obtained with MATRIX by 5 ± 2 o / oo for W + -and by 4 ± 2 o / oo for W − -production. Central Z-boson predictions would move upwards by 2 ± 1 o / oo and the ones for forward Z-bosons by 7 ± 3 o / oo . The uncertainty in those rescaling factors comes from the extrapolation uncertainty in σ extrapolated NNLO . Such shifts decrease, but do not eradicate the differences. The MATRIX results with the smaller value r min cut = 0.05% lead to better agreement with the FEWZ results, i.e., there are systematic upward shifts in Fig. 6. In detail, these amount to a few per mill for η ll ≤ 1.0 for Z/γ * -production at central rapidities and up to a few per cent for forward Z/γ * -production in the bins with η ll 2.0. The computational demands for these MATRIX runs were huge. 8 The suggested rescaling factor σ extrapolated NNLO /σ r cut NNLO turns out to be unity within the numerical accuracy of our computation for central Z/γ * -production. Predictions for forward Z-bosons would be shifted upwards uniformly by 3 ± 2 o / oo and the observed differences, especially in the first η ll bins, would still persist.
Finally, in Fig. 7 we repeat the benchmark study with NNLO predictions obtained with the MCFM code, in which case the numerical integration accuracy is typically O(1 o / oo ) and negligible in the plots. We do find substantial deviations of the MCFM results at NNLO, being below the FEWZ ones for all distributions considered. Differences amount to O(3%) for the η l -distribution in W + -production and up to O(2%) for W − -production, respectively. For central Z/γ * -production the η ll -distribution is also O(2%) below the FEWZ results for η ll ≤ 1.  Only the results with the smallest slicing cuts are plotted: MATRIX with r cut = 0.15% for W ± → l ± ν and r cut = 0.05% for Z → l + l − production; MCFM with τ cut = 4 · 10 −4 .
forward Z/γ * -production. In the first bins of the latter the deviations grow up to O(20%). As discussed, MCFM uses N-jettiness subtraction and allows for different τ cut choices for the jettiness slicing parameter. We use the default value, τ cut = 6 · 10 −3 and two smaller ones, τ cut = 1 · 10 −3 and τ cut = 4 · 10 −4 , the limitation being here the goal to reach an integration accuracy of a few units in 10 −4 in reasonable time 9 with given computational resources. The decreasing values of τ cut display the expected trend clearly in Fig. 7, namely, the smaller the choice of τ cut , the closer the MCFM result to that by FEWZ. Nevertheless, the differences remain. In order to compare those differences easier, we collect the best prediction for each code at NNLO in a single figure in Fig. 8.
Given the level of agreement among the predictions at NLO accuracy, the deviations observed in Figs. 5-7 need to be put into perspective by looking at the size of the pure NNLO corrections alone, which we define bin-by-bin through the deviation of the NNLO K-factor from one, δ NNLO = (σ NNLO /σ NLO − 1). Typically pure NNLO corrections δ NNLO are rather small, and we illustrate those only in the case of largest corrections. For W + -production δ NNLO amounts to a few per mill for η l 1 and grows to O(1 − 2%) for larger rapidities η l 1, while instead for W − -production δ NNLO is of the size O(1%) for η l 1 and increases to a few per cent for larger rapidities. For the central Z/γ * -production the NNLO corrections δ NNLO are only a few per mill for η ll 1.5 and grow to O(2 − 3%) for larger di-lepton rapidities. Thus, the observed differences between considered codes are actually similar in size to that of the pure NNLO corrections, even exceeding them at times. The case of forward Z/γ * -production features larger higher order corrections and will be discussed in detail next. The comparable size of the NNLO corrections and differences among the predictions signal that the numerical precision of the studied computer programs does not match the formal accuracy of predictions at NNLO.
In Fig. 9 we focus on the η ll -distribution for the forward Z/γ * -production obtained by AT-LAS [6]. The particular fiducial cuts on the decay leptons for these data lead to sizeable QCD corrections at higher orders, which we illustrate in Fig. 9, where we display at LO, NLO and NNLO accuracies obtained with DYNNLO (left) and with FEWZ (right). The same comparison is performed in Fig. 10 for the results of the MATRIX and the MCFM codes, where we display σ NNLO with the smallest slicing cuts, r min cut = 0.05% and τ cut = 4 · 10 −4 . As already remarked above, we use ABMP16 PDFs at NNLO in all cases, independent of the perturbative order. Figs. 9 and 10 clearly illustrate the significant corrections up to O(50%) in first bins, when increasing the perturbative order from LO to NLO, while the change from NLO to the NNLO QCD predictions still amounts to corrections of O(5 − 10%) in some η ll bins. The LO results in Figs. 9 and 10 are all in perfect agreement and the deviations in the NLO predictions by DYNNLO have already been discussed above.
The observed pattern of the higher order corrections for the predictions with FEWZ in Fig. 9 (right) and with MATRIX in Fig. 10 (left) is very similar. The overall offset of the pulls for the MATRIX results with r min cut = 0.05% compared to the FEWZ ones is small in the entire η ll range except for the first η ll bins and originates from the different NNLO cross sections as illustrated in Fig. 6. In contrast, the cross sections σ LO , σ NLO and σ NNLO from DYNNLO in Fig. 9 (left) and from MCFM in Fig. 10 (right) show a different trend. The pulls for the ATLAS data follow rather closely the respective NLO predictions across the entire range in rapidities. The NNLO predictions from those codes do undershoot the data by several per cent, which causes the significant deviations displayed in Figs. 5 and 7. Next we continue the benchmark studies with DØ data on the electron charge asymmetry distribution A e , which has been obtained as a function of the electron pseudo-rapidity from W ± -boson production at √ s = 1.96 TeV at the Tevatron [8]. This observable is also subject to larger higher order corrections so that we illustrate again the size of the LO, NLO and NNLO predictions obtained, as before, in all cases with the NNLO ABMP16 PDFs and α (n f =5) s (M Z ) = 0.1147 and we plot the difference to the NNLO predictions computed with the FEWZ code. The DØ data had already been included in the fit of the ABMP16 PDFs and a good description of those data in the fit had been reached.
In Fig. 11 we plot in addition to the DØ data on A e the LO, NLO and NNLO predictions by the DYNNLO code (left) and the MATRIX code (right), keeping again a relative numerical integration accuracy of a few units in 10 −4 for the respective W ± -boson cross sections. The LO and NLO curves illustrate the sizable higher order corrections and those predictions agree among these codes. With the given accuracy of the DØ data on A e , also the NNLO corrections are relevant, but we see both, the DYNNLO and the MATRIX results (here with r min cut = 0.15%) being mostly above the FEWZ numbers. The deviations increase with increasing electron pseudo-rapidity η e and become significant for η e 1.0, where the size of the difference exceeds the size of the pure NNLO corrections. For the asymmetry A e any overall rescaling of cross sections as suggested for the MATRIX code and described in Sec. III to account for r min cut dependence has no effect. In Fig. 12 we show the same study, now comparing to the results obtained with the MCFM code. The NNLO MCFM result has been computed with the default τ cut value, τ cut = 6 · 10 −3 , and the numerical integration accuracy of the individual W ± -boson cross sections is typically O(1 o / oo ). In addition, deviating from the default settings of MCFM, the parameter cutoff has been changed to 10 −6 from 10 −9 , which is its default value 10 . The parameter cutoff provides the minimum value on any dimensionless variables, for instance, any invariant mass squared s i j of any pair   relative to FEWZ for the electron charge asymmetry distribution A e in the case of staggered cuts with p e T ≥ 35 GeV and p ν T ≥ 25 GeV. The NNLO results for MCFM were obtained with cutoff = 10 −6 as in the case of symmetric cuts to avoid numerical instabilities. Interestingly, the very good agreement among all codes already observed at NLO in Fig. 4, is found now in all cases also at NNLO. This is in complete contrast to the case of symmetric cuts in Figs. 11 and 12.

III. POWER CORRECTIONS
Our comparison in the previous section was based on computer codes implementing different approaches to the regularization of double real singular emissions. The global slicing methods neglect power corrections, hence one may assume that those are at least partly responsible for the observed differences, which we explore in this section. There are two sources of power corrections, namely those intrinsic to the q T and T N factorization on the one hand, which have been studied extensively in the literature before, see e.g. [27][28][29][30], and fiducial power corrections on the other, which have been studied in detail more recently [31]. The latter ones arise whenever one considers fiducial cuts or leptonic observables and formally dominate.

A. Review of global slicing methods
We briefly review global slicing methods applied to the Drell-Yan process and with emphasis on the power counting in the slicing parameter τ, see e.g., [28,31]. Specifically, we focus on the q T and N-jettiness subtractions [2,21,22] for the hadro-production of a gauge boson V, which decays to a leptonic final state L with the following kinematics where a, b are the initial hadrons with momenta p a,b , q is the gauge boson momentum and X(k i ) denotes hadronic final states with momenta k i . In the laboratory frame, spanned by light-like vectors n µ = (1, 0, 0, 1) andn µ = (1, 0, 0, −1), the initial hadron momenta p a,b can be parametrized in Born kinematics in terms of Q = q 2 and the rapidity of the gauge boson Y as The slicing parameter τ for q T -subtraction is defined by in terms of the transverse momenta k T,i of the hadronic real emission radiation. Likewise, one has for the leptonic 0-jettiness T 0 (see [28,31] for other variants of 0-jettiness definitions). The leptonic T 0 is preferred due to smaller intrinsic power corrections, cf. [32] and used in MCFM [20].
The definitions of τ in eqs. (5) and (6) vanish at Born level and resolve additional radiation in an infrared-safe manner, so that the phase space integration for the cross section can be written as where τ cut is the cut for the slicing of the phase space. The dependence of dσ/dτ on τ can be predicted from the universal factorization of QCD in soft and collinear limits. It has the structure where the +-distributions are the well-known leading threshold logarithms and the terms proportional to τ p−1 with p > 0 are integrable and denote power corrections in the soft and collinear limit. From the analytical integration one obtains for σ(τ cut ) schematically The crucial point to stress here is the scaling behavior of the power corrections τ p cut , i.e. the value of the exponent p. For the production of a stable gauge boson V, p takes positive integer values, while the subsequent decay with cuts on the leptonic final state changes the scaling of the power corrections [31], such that p rises in steps of half-integers, i.e., p = 1/2, 1, 3/2 and so on. This will be discussed in more detail below.
The scaling of the power corrections has consequences for the particular subtraction scheme, which is then implemented via a global subtraction term σ sub (τ cut ) as σ = σ sub (τ cut ) + τ cut dτ dσ dτ + ∆σ sub (τ cut ) .
Here the term ∆σ sub (τ cut ) = σ(τ cut ) − σ sub (τ cut ) parametrizes the residual power corrections. It is neglected in slicing methods, giving rise to an intrinsic error of these methods. If the global subtraction term σ sub (τ cut ) cancels only the leading soft and collinear singularities in σ(τ cut ) in eq. (9), then the residual power corrections in the presence of cuts on the decay leptons scale as √ τ cut . This implies enhanced corrections of the order q T /Q for the q T subtraction, as will be explained below, or of the order of √ T 0 /Q for the N-jettiness subtraction, as detailed in [31] with a power counting argument.
The phase space slicing codes under consideration employ different strategies for dealing with power corrections. MATRIX performs an extrapolation of r cut = q T /M V → 0 for the total rate of the process computed with q T -subtraction by evaluating the cross section at fixed values in the interval r cut ∈ [0.15, 1]% in steps of 0.01%. It is then recommended to correct the kinematic distributions by rescaling uniformly with the ratio σ extrapolated NNLO /σ r cut NNLO . MCFM has improved the τ cut dependence by implementing the leading power corrections of [32,33] (see also [28]), which are derived for the production of stable gauge bosons V and scale as τ cut , cf. eq. (9). In addition, MCFM computes the cross section for an array of different τ cut values and performs an automated fitting of the τ cut dependence at NNLO with the following ansatz σ(τ cut ) NNLO = σ 0 + c 1 τ cut ln 3 (τ cut /M V ) + c 2 τ cut ln 2 (τ cut /M V ) + c 3 τ cut , where c i are the fit parameters and the result is then extrapolated to τ cut → 0. Note, that the functional form in eq. (11) does not capture well the scaling of the leading power corrections proportional to √ τ cut in the case of gauge boson decays [31].

B. Fiducial cuts
For the discussion of the fiducial cuts on the decay leptons, we follow the presentation in [31,34]. First, we need to further specify the leptonic final state L in eq. (3), which reads Z/γ * → L(q) = l 1 (p 1 ) + l 2 (p 2 ) , and p 1,2 are the lepton momenta, q = p 1 + p 2 . In the presence of hadronic final states X(k i ) from additional real emission radiation in eq. (3), using q = p a + p b − i k i , the gauge boson momentum can be expressed through Q, Y and a non-vanishing transverse momentum q T , such that in components where m T = Q 2 + q 2 T and the azimuthal angle φ in the transverse plane is given by p T 1 · q T = p T 1 q T cos φ. Momentum conservation yields for the transverse momenta and rapidities of the leptons and η 1 = Y + ∆y , The leptonic final state phase space Φ L (neglecting lepton masses) reads in terms of the variables φ and ∆y in eq. (13), and p T 1 implicitly depends on φ and ∆y through eq. (14). The fiducial cuts on the decay leptons applied to the data discussed in Sec. II modify eq. (16) by constraining the integration range.
With the typical cuts on the transverse momenta and rapidities of the leptons, p T 1,2 ≥ p min T and η min 1,2 ≤ η 1,2 ≤ η max 1,2 , the phase space Φ L becomes It has been pointed out in [31] that the presence of cuts on the leptons' transverse momenta breaks azimuthal symmetry and leads to linear power corrections in q T . The expansion of eqs. (14) and (15) for small q T up to quadratic corrections in q T gives and The symmetric cut on the transverse momenta p T 1,2 ≥ p min T used for the ATLAS or DØ data considered above thus splits the φ integration range depending on the sign of cos φ as such that the θ-functions in eq. (17) give rise to different integrands in the respective regions and the linear power corrections in q T do not vanish in the phase space integral. Considering only the cut p T 1,2 ≥ p min T on the transverse momenta, the phase space Φ L becomes In case of rapidity cuts, the constraints are slightly more involved, see [34]. Considering the cuts |η 1,2 | ≤ η max for the selection of both leptons at central pseudo-rapidity, as applied to the ATLAS data for Z/γ * -boson production, the phase space Φ L at Born level becomes For small Y this constraint is less tight compared to the p T 1,2 cuts in eq. (21) for the typical values of p min T used. To assess where the θ-functions in eq. (17) involving the pseudo-rapidities affect the phase space Φ L and, in particular, break azimuthal symmetry, one has to examine the regions where they overlap. The boundary of this region is given by |η 1 | = |η 2 |, which determines a value φ * through the condition The constraint | cos φ * | ≤ 1 restricts either the region of |η 1 | or |η 2 | in the phase space (17), but not both. Thus, azimuthal symmetry is not broken by the pseudo-rapidity cuts, but may still be affected by p T cuts through eq. (20). Solutions in the physical range for | cos φ * | < 1 imply the following scaling for the values of q T /Q which can be approximated for small Y (and using This defines a lower bound on q T for linear power corrections to appear as a result of broken azimuthal symmetry due to the pseudo-rapidity cuts. For values q T < q * T azimuthal symmetry is restored and only quadratic power corrections arise. As eqs. (23)- (25) indicate, the transition between these two regions of q T is sharp up to corrections.
The appearance of linear power corrections in the phase space Φ L can be illustrated by considering the deviations |1 − Φ L (q T )/Φ L (0)| from the Born level leading power results for q T = 0. In Fig. 15 we show them for the fiducial cuts applied to the ATLAS data in case of Z/γ * -boson production. On the left in Fig. 15, the leptons are selected at central pseudo-rapidities |η l i | ≤ 2.5 for i = 1, 2 and we observe the presence of linear power corrections in q T for central gauge boson pseudo-rapidities η ll 1 due to the p T constraint in eq. (20). For r cut = q T /Q = 0.15%, which is the default value for r cut used in MATRIX as a slicing cut and indicated by the vertical dashed line in Fig. 15, their size amounts to O(0.5 o / oo ). In contrast, for larger η ll the pseudo-rapidity constraints dominate the phase space Φ L and azimuthal symmetry is restored, resulting in quadratic power corrections in q T for small enough q T , see eqs. (23)- (25). In Fig. 15 on the left this feature is illustrated for η ll = 1.2 and 1.8, and the corrections to Φ L for r cut = 0.15% are smaller by more than two orders of magnitude. It is interesting to compare these findings with the η ll dependence of the differences at NLO of DYNNLO from codes using local subtraction in Fig. 2 and with the deviations at NNLO of DYNNLO, MATRIX (with r cut = 0.15%) and MCFM from FEWZ in Figs. 5, 6 and 7 for central Z/γ * -boson production. For η ll 1 all slicing codes undershoot FEWZ, while they tend to agree well for η ll 1.5. The transition around η ll 1.2 when the linear power corrections in q T in Φ L vanish, is most pronounced in the case of MCFM in Fig. 7.  Fig. 15 for the fiducial cuts applied to ATLAS data set [6] for W ± -boson production (Q = M W ) and different values of the lepton pseudo-rapidity η l . For the lepton momenta p l,ν T ≥ 25 GeV are required.
In Fig. 15 on the right we plot the same study for the ATLAS cuts with one lepton at central and the other at forward pseudo-rapidity. In this case, due to the non-overlapping regions |η l 1 | ≤ 2.5 and 2.5 ≤ |η l 2 | ≤ 4.9 azimuthal symmetry is always broken by the p T constraint in eq. In Fig. 16 we plot the phase space Φ L for the fiducial cuts applied to the ATLAS data set [6] for W ± -boson production. In this case, the binning in the lepton pseudo-rapidity fixes η l and trivially fulfills one of the θ-functions in the integral for Φ L in eq. (17). The other θ-function constrains the neutrino's pseudo-rapidity in the whole integration range, so that the pseudo-rapidity cuts do not restore the azimuthal symmetry for small q T . The linear power corrections in q T , which we observe in Fig. 16 originate from the constraints p l,ν T ≥ 25 GeV for the lepton momenta, which do break azimuthal symmetry as shown in eq. (20). They amount to corrections O(0.4 − 0.8 o / oo ) for the value r cut = 0.15%, depending on the value of the lepton pseudo-rapidity η l , with larger corrections observed for central η l . This pattern is in line with the observed deviations in Fig. 1 at NLO between DYNNLO and local subtraction results, and in Figs. 5, 6 and 7 at NNLO between DYNNLO, MATRIX and MCFM on the one and FEWZ on the other side, where all slicing codes give consistently lower results than FEWZ and the deviations display little dependence on the lepton pseudo-rapidity η l .
Finally, we briefly discuss staggered cuts on the transverse momenta, when p T 1 ≥ p min T + δp T and p T 2 ≥ p min T for some δp T > 0, as realized for instance in the DØ measurement of the electron charge asymmetry distribution A e discussed above. For staggered cuts with δp T the phase space Φ L evaluates at Born level then as  Fig. 15 for the fiducial cuts applied to the DØ data [8] for the electron charge asymmetry distribution A e , using staggered cuts p T 1 ≥ p min T + δp T , p T 2 ≥ p min T and Q = M W for the electron pseudorapidity fixed at η e = 0.
In addition, with staggered cuts the θ-functions in eq. (17) which constrain the fiducial phase space give rise to the following condition for the azimuthal integration which reproduces the case of symmetric cuts in eq. (20) for δp T → 0. On the other hand, it is obvious that for δp T > q T the θ-functions for the cuts on the transverse momenta in eq. (17) have no regions of common overlap, and therefore, do not affect the azimuthal symmetry, which is unbroken in this case. The typical scales for the slicing cuts r cut or τ cut in the slicing codes imply that q min T 1 GeV, while measurements with staggered cuts in the experiments use values of δp T of the order of a few GeV. Therefore, due to unbroken azimuthal symmetry, linear power corrections are largely absent.
We illustrate the effect of the staggered cuts in Fig. 17 for the fiducial cuts applied to the DØ data [8] on the electron charge asymmetry distribution A e . We apply a series of cuts δp T = 0, 0.1 and 10 GeV for an electron pseudo-rapidity of η e = 0. Other choices of η e give qualitatively the same results, see Fig. 16. As the value of δp T increases, the deviations from the linear power corrections for the case of δp T = 0 GeV become apparent. For δp T = 10 GeV, which corresponds to the choice in the DØ data selection of p e T ≥ 35 GeV and p ν T ≥ 25 GeV, we observe in Fig. 16 only quadratic power corrections for small q T . This is in line with the general good agreement between the results of DYNNLO, MATRIX, MCFM and FEWZ observed in Figs. 13 and 14. In addition, the observed pattern for the convergence of slicing codes in the limit of vanishing slicing cuts r cut or τ cut also agrees with studies performed in the presence of staggered cuts with MATRIX in [18] and with MCFM in [20].
In summary, the changes in the relative size of the power corrections in Φ L are correlated with better or worse agreement between the cross sections generated with phase space slicing codes (DYNNLO, MATRIX and MCFM) and the one with a local subtraction (FEWZ). This holds in particular for the dependence on the gauge boson pseudo-rapidity η ll in the case of Z/γ * -boson production and the effect of staggered cuts. This observation does not prove that the neglected power corrections are the pure source of the differences among the predictions of the various codes. Nevertheless, it is at least a warning sign that the computation of the NNLO corrections can only be considered a solved problem if the numerical precision of the Monte Carlo integration is under better control than the size of the correction itself, which calls for an improvement of the presently available codes and/or subtraction methods at NNLO accuracy. The linear power corrections can be uniquely predicted by factorization. This fact has been used in [34] to propose a method for the inclusion of all fiducial-cut induced power corrections in the framework of q Tsubtraction schemes and has been applied to provide cross section predictions for Higgs-boson production with fiducial cuts [35].

IV. CONCLUSIONS
We have investigated the accuracy of available NNLO QCD predictions for the hadroproduction of W ± -and Z-bosons, including their leptonic decays and keeping fully exclusive kinematics. Such predictions are available from several publicly available codes and we have chosen DYNNLO, FEWZ, MATRIX and MCFM to compute benchmark values for kinematics which are representative for measurements of differential distributions in the pseudo-rapidities of the decay leptons from the LHC and Tevatron. The uncertainties in the cross sections from the numerical Monte Carlo integration have been limited to few units in 10 −4 and are negligible in all cases. At NLO there is perfect agreement among the results from FEWZ, MATRIX and MCFM and, partially with those of DYNNLO as well. However, at NNLO accuracy we found differences among the predictions by the same codes that are comparable to the NNLO correction itself. We demonstrated that the observed systematic differences can be understood in terms of the subtraction schemes employed, FEWZ using a fully local subtraction scheme and DYNNLO, MATRIX and MCFM applying phase space slicing schemes. We have illustrated, how the fiducial cuts on the transverse momenta and pseudorapidities of the decay leptons lead to linear power corrections in the slicing parameter, i.e. q T /Q and √ T 0 /Q for the q T and N-jettiness subtractions. Also the deviations share certain patterns across the range of pseudo-rapidities in the considered distributions, which have been correlated with the appearance of linear power corrections in the lepton decay phase space Φ L as a function of q T . The latter serves as a simple and efficient model to study power corrections in cross sections for the gauge boson production with their subsequent leptonic decay. For most of the distributions considered the pure NNLO QCD corrections on top of the NLO ones are rather small, often in the range of O(1%), while the deviations among the codes investigated in this study are not substantially smaller, often even of the same size or larger, hinting towards a significant intrinsic uncertainty in the computation of the NNLO QCD corrections for those observables.
In summary, with the continuous increase in the precision of the experimental measurements, the theory predictions are pressed to provide cross sections at NNLO (or beyond) where the systematic uncertainties due to choices of particular schemes or algorithms for the computation can be safely neglected in comparison to the experimental uncertainties. The results of our study call for further improvements in this direction.