1 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^{\pm }\)- and Z-boson 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^{\pm }\)- 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^{\pm }\)- and Z-boson hadro-production 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 pseudo-rapidity of the decay leptons and typically have a precision of \({{\mathcal {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 \({{\mathcal {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 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 Sect. 2 we present the benchmark numbers for two representative sets of \(W^{\pm }\)- 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 Sect. 3 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 Sect. 4.

2 Benchmark computations

2.1 Set-up and validation

The set-up for benchmarking available QCD predictions for \(W^{\pm }\)- 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^{\pm }\)- 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^\pm \)- and \(Z/\gamma ^*\)-production cross sections [6] measured at a center-of-mass energy of \(\sqrt{s}=7\) TeV at the LHC. These data are given in form of pseudo-rapidity distributions for the decay electron or muon (\(W^{\pm }\)-production) and the decay lepton-pair (\(Z/\gamma ^*\)-production), respectively. The transverse momenta \(p_T^l\) and the pseudo-rapidities \(\eta _l\) of the decay leptons are subject to fiducial cuts. The cross sections for \(Z/\gamma ^*\)-production are measured at central as well as at forward pseudo-rapidities.

  • The data obtained by DØ on \(W^\pm \)-boson production at \(\sqrt{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_T^{l,\nu }\) 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^{\pm }\)-boson production at \(\sqrt{s}=8\) TeV at the LHC [9] has similar experimental precision as the chosen data set [6] collected at \(\sqrt{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_\mu \) scheme with input values \(G_F\), \(M_Z\), \(M_W\) and with \(\sin ^2(\theta _W)\) and the QED coupling \(\alpha (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]

$$\begin{aligned} \begin{array}{ll} ~~G_{\mu } = 1.16637\times 10^{-5} \; \mathrm{GeV}^{-2} \, , \\ ~M_Z = 91.1876 \; \mathrm{GeV}\, , \quad \quad &{}~\Gamma _{Z} = 2.4952 \; \mathrm{GeV} \, , \\ M_W = 80.379\; \mathrm{GeV}\, , \quad \quad &{}\Gamma _W = 2.085\; \mathrm{GeV} \, , \end{array} \end{aligned}$$
(1)

and for the relevant CKM parameters

$$\begin{aligned}&|V_{ud}| \,=\, 0.97401\, , \qquad |V_{us}| \,=\, 0.2265 \, , \nonumber \\&|V_{cd}| \,=\, 0.2265 \, , \qquad |V_{cs}| \,=\, 0.97320 \, , \nonumber \\&|V_{ub}| \,=\, 0.00361\, , \qquad |V_{cb}| \,=\, 0.04053 \, . \end{aligned}$$
(2)

The computations are performed in the \(\overline{\mathrm{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, \(\alpha _s^{(n_f=5)}(M_Z) = 0.1147\). These choices do not bias any of the predictions we discuss below. The renormalization and factorization scales \(\mu _R\) and \(\mu _F\) are set to \(\mu _R = \mu _F = M_V\), where \(M_V\) is the mass of the gauge boson V.

Fig. 1
figure 1

The NLO QCD cross sections for inclusive \(pp \rightarrow W^\pm +X \rightarrow l^\pm \nu + X\) as function of pseudo-rapidity \(\eta _l\), computed with DYNNLO, MATRIX and MCFM relative to FEWZ and using the ABMP16 PDFs. Cuts of \(p_T^{l,\nu } \ge 25\) GeV and \(M_T \ge 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 \(\eta _l\)-bin, while DYNNLO and MCFM results are shifted slightly to the left and the right

We run the following publicly available codes designed for the computation of the fully differential NNLO QCD predictions for the lepton rapidity distributions.

  • DYNNLO (version 1.5) [2, 3]Footnote 1 (No built-in computation of PDF uncertainties available.)

  • FEWZ (version 3.1) [16, 17]Footnote 2

  • MATRIX (version 1.0.4) [18]Footnote 3 (No built-in computation of PDF uncertainties available.)

    MATRIX uses the scattering amplitudes from OpenLoops [19].

  • MCFM (version 9.0) [20]Footnote 4

    MCFM uses the implementation of the NNLO computation of Ref. [5].

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_{\mathrm{cut}}\) for MATRIX as a cut on \(q_T\) and \(\tau _{\mathrm{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].Footnote 5DYTurbo 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].Footnote 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^\pm \)- and \(Z/\gamma ^*\)-production cross sections at \(\sqrt{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 \(\alpha _s^{(n_f=5)}(M_Z) = 0.1147\), independent of the order of perturbation theory for the cross sections \(\sigma _{\mathrm{LO}}\), \(\sigma _{\mathrm{NLO}}\) and \(\sigma _{\mathrm{NNLO}}\). We have found excellent agreement for all cross sections computed at the leading order (LO, not shown here) at \({{\mathcal {O}}}(10^{-5})\). At NLO we can directly compare and quantify the accuracy of the numerical integrations, as all codes employ fully local subtraction schemes except for DYNNLO, which uses \(q_T\)-subtraction also at NLO.Footnote 7FEWZ applies sector decomposition while MATRIX and MCFM all use by default the dipole subtraction [26].

Fig. 2
figure 2

Same as Fig. 1 for the NLO QCD cross sections for \(pp \rightarrow Z/\gamma ^*+X \rightarrow l^+l^- + X\) production at \(\sqrt{s}=7\) TeV as function of the pseudo-rapidity \(\eta _{ll}\). Cuts of \(p_{T1,2} \ge 25\) GeV, \(66 \ge M_{ll} \ge 116\) GeV and \(|\eta _{l_i}| \le 2.5\), \(i=1,2\) for the lepton pseudo-rapidities are applied

Fig. 3
figure 3

Same as Fig. 2 with one lepton at forward pseudo-rapidities and cuts of \(|\eta _{l_1}| \le 2.5\) and \(2.5 \le |\eta _{l_2}| \le 4.9\). Predictions by DYNNLO, MATRIX and MCFM relative to FEWZ (left) and zoom on MATRIX and MCFM results (right)

In Fig. 1 we plot the results for the \(W^\pm \)-production cross sections at \(\sqrt{s}=7\) TeV at NLO corresponding to the kinematics of the ATLAS data set [6] and we find agreement at the level of between FEWZ, MATRIX and MCFM in the entire \(\eta _l\) range, while for DYNNLO we find the values to be systematically enhanced by for both, \(W^+\)- and \(W^-\)-production. In Figs. 2 and 3 we show the NLO cross sections for \(Z/\gamma ^*\)-production at \(\sqrt{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 \(|\eta _{l_i}| \le 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 for \(\eta _{ll} \lesssim 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, \(|\eta _{l_1}| \le 2.5 \le |\eta _{l_2}| \le 4.9\). We observe agreement at the level of 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 \(\eta _{ll}\) bins. Finally, in Fig. 4 we plot the electron charge asymmetry distribution \(A_e\) in \(W^{\pm }\)-boson production at \(\sqrt{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_T^{e,\nu } \ge 25\) GeV, and on the right staggered cuts, \(p_T^{e} \ge 35\) GeV and \(p_T^{\nu } \ge 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 \({{\mathcal {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 \(\eta _{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 \(\eta _{e} \gtrsim 2.5\). The individual cross sections for \(W^\pm \)-production are computed at an accuracy of , 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 \(\eta _l (\eta _{ll})\), which will be addressed in Sect. 3 below in the light of power corrections in the slicing parameter.

Fig. 4
figure 4

The ratio of the NLO QCD corrections to the electron charge asymmetry distribution \(A_e\) in \(W^{\pm }\)-boson production at \(\sqrt{s}=1.96\) TeV computed with DYNNLO, MATRIX and MCFM to ones by FEWZ. Cuts are applied as in the DØ data selection [8]: \(p_T^{e,\nu } \ge 25\) GeV symmetric (left) and \(p_T^{e} \ge 35\) GeV and \(p_T^{\nu } \ge 25\) GeV staggered (right)

2.2 NNLO benchmark predictions

We now proceed to the QCD predictions at NNLO accuracy, again starting with \(W^\pm \)- and \(Z/\gamma ^*\)-production cross sections at \(\sqrt{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.

Fig. 5
figure 5

The pulls for the ATLAS data measured in inclusive \(pp \rightarrow W^\pm +X \rightarrow l^\pm \nu + X\) and \(pp \rightarrow Z/\gamma ^*+X \rightarrow l^+l^- + X\) production at \(\sqrt{s}=7\) TeV [6] with the statistical (inner bar) and the total uncertainties, including the systematic ones. The fiducial cuts on the decay leptons in the final state are indicated in the figure. The ABMP16 central predictions at NNLO are obtained with FEWZ and the deviations of the predictions from DYNNLO are shown (dashed) for comparison

Fig. 6
figure 6

Same as Fig. 5 using predictions by the MATRIX code with different values for the \(q_T\)-slicing cut: \(r^{\mathrm{min}}_{\mathrm{cut}}=0.15\%\) (dashed) and \(r^{\mathrm{min}}_{\mathrm{cut}}=0.05\%\) (dashed-dotted)

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-pseudo-rapidity (\(\eta _l\)) distribution in \(W^+\)-production in Fig. 5, the DYNNLO results are below the FEWZ ones by up to \({{\mathcal {O}}}(1\%)\) and by a few per mill for \(W^-\)-production, in both cases over the whole range in \(\eta _l\). For the di-lepton pseudo-rapidity distribution in the central \(Z/\gamma ^*\)-production the DYNNLO predictions are below the FEWZ ones by several per mill up to \(\eta _{ll} \le 1.5\) and tend to agree better with the FEWZ ones for larger \(\eta _{ll}\). For forward \(Z/\gamma ^*\)-production instead, we see the DYNNLO results being below the FEWZ ones by up to \({{\mathcal {O}}}(1-2\%)\) in the entire \(\eta _{ll}\) range, with significant deviations up to \({{\mathcal {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^{\mathrm{min}}_{\mathrm{cut}}=q_T^{\mathrm{min}}/M_V=0.8\%\) for the slicing cut on \(q_T\). The MATRIX code uses \(r^{\mathrm{min}}_{\mathrm{cut}}=0.15\%\) as the default and offers also the choice \(r^{\mathrm{min}}_{\mathrm{cut}}=0.05\%\) for the \(Z/\gamma ^*\)-production process. In detail we find MATRIX results being below the FEWZ ones by up to \({{\mathcal {O}}}(1\%)\) for the \(\eta _l\)-distribution in \(W^+\)-production, and by a few per mill for \(W^-\)-production. For the \(\eta _{ll}\)-distribution in \(Z/\gamma ^*\)-production at central rapidities we see the MATRIX numbers with the default value \(r^{\mathrm{min}}_{\mathrm{cut}}=0.15\%\) below the FEWZ ones by a few per mill for \(\eta _{ll} \le 1.0\), above the FEWZ ones in the bins around \(\eta _{ll} \simeq 1.5\) and in agreement with FEWZ for larger rapidities. On the other hand, the MATRIX results with \(r^{\mathrm{min}}_{\mathrm{cut}}=0.15\%\) for forward \(Z/\gamma ^*\)-production are below the FEWZ ones by up to \({{\mathcal {O}}}(1-2\%)\) in the entire \(\eta _{ll}\) range and show a significant deviation of \({{\mathcal {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^{\mathrm{min}}_{\mathrm{cut}} \rightarrow 0\) and suggests a uniform rescaling of each bin by the ratio \(\sigma ^{\mathrm{extrapolated}}_{\mathrm{NNLO}}/\sigma ^{r_{\mathrm{cut}}}_{\mathrm{NNLO}}\), see also Sect. 3. 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 \pm 2\) for \(W^+\)- and by \(4 \pm 2\) for \(W^-\)-production. Central Z-boson predictions would move upwards by \(2 \pm 1\) and the ones for forward Z-bosons by \(7 \pm 3\) . The uncertainty in those rescaling factors comes from the extrapolation uncertainty in \(\sigma ^{\mathrm{extrapolated}}_{\mathrm{NNLO}}\). Such shifts decrease, but do not eradicate the differences. The MATRIX results with the smaller value \(r^{\mathrm{min}}_{\mathrm{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 \(\eta _{ll} \le 1.0\) for \(Z/\gamma ^*\)-production at central rapidities and up to a few per cent for forward \(Z/\gamma ^*\)-production in the bins with \(\eta _{ll} \lesssim 2.0\). The computational demands for these MATRIX runs were huge.Footnote 8 The suggested rescaling factor \(\sigma ^{\mathrm{extrapolated}}_{\mathrm{NNLO}}/\sigma ^{r_{\mathrm{cut}}}_{\mathrm{NNLO}}\) turns out to be unity within the numerical accuracy of our computation for central \(Z/\gamma ^*\)-production. Predictions for forward Z-bosons would be shifted upwards uniformly by \(3 \pm 2\) and the observed differences, especially in the first \(\eta _{ll}\) bins, would still persist.

Fig. 7
figure 7

Same as Fig. 5 using predictions by the MCFM code and with different values for the jettiness slicing parameter: \(\tau _{\mathrm{cut}}=6\cdot 10^{-3}\) (dashed), \(\tau _{\mathrm{cut}}=1\cdot 10^{-3}\) (dotted), and \(\tau _{\mathrm{cut}}=4\cdot 10^{-4}\) (dashed-dotted)

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 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 \({{\mathcal {O}}}(3\%)\) for the \(\eta _l\)-distribution in \(W^+\)-production and up to \({{\mathcal {O}}}(2\%)\) for \(W^-\)-production, respectively. For central \(Z/\gamma ^*\)-production the \(\eta _{ll}\)-distribution is also \({{\mathcal {O}}}(2\%)\) below the FEWZ results for \(\eta _{ll} \le 1.5\) and up to \({{\mathcal {O}}}(2-3\%)\) for forward \(Z/\gamma ^*\)-production. In the first bins of the latter the deviations grow up to \({{\mathcal {O}}}(20\%)\). As discussed, MCFM uses N-jettiness subtraction and allows for different \(\tau _{\mathrm{cut}}\) choices for the jettiness slicing parameter. We use the default value, \(\tau _{\mathrm{cut}}=6\cdot 10^{-3}\) and two smaller ones, \(\tau _{\mathrm{cut}}=1\cdot 10^{-3}\) and \(\tau _{\mathrm{cut}}=4\cdot 10^{-4}\), the limitation being here the goal to reach an integration accuracy of a few units in \(10^{-4}\) in reasonable timeFootnote 9 with given computational resources. The decreasing values of \(\tau _{\mathrm{cut}}\) display the expected trend clearly in Fig. 7, namely, the smaller the choice of \(\tau _{\mathrm{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.

Fig. 8
figure 8

Compilation of the NNLO theory predictions of Figs. 5, 6 and 7. Only the results with the smallest slicing cuts are plotted: MATRIX with \(r_{\mathrm{cut}}=0.15\%\) for \(W^\pm \rightarrow l^\pm \nu \) and \(r_{\mathrm{cut}}=0.05\%\) for \(Z\rightarrow l^+l^-\) production; MCFM with \(\tau _{\mathrm{cut}}=4\cdot 10^{-4}\)

Given the level of agreement among the predictions at NLO accuracy, the deviations observed in Figs. 5, 6 and 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, \(\delta _{\mathrm{NNLO}} = (\sigma _{\mathrm{NNLO}}/\sigma _{\mathrm{NLO}}-1)\). Typically pure NNLO corrections \(\delta _{\mathrm{NNLO}}\) are rather small, and we illustrate those only in the case of largest corrections. For \(W^+\)-production \(\delta _{\mathrm{NNLO}}\) amounts to a few per mill for \(\eta _{l} \lesssim 1\) and grows to \({{\mathcal {O}}}(1-2\%)\) for larger rapidities \(\eta _{l} \gtrsim 1\), while instead for \(W^-\)-production \(\delta _{\mathrm{NNLO}}\) is of the size \({{\mathcal {O}}}(1\%)\) for \(\eta _{l} \lesssim 1\) and increases to a few per cent for larger rapidities. For the central \(Z/\gamma ^*\)-production the NNLO corrections \(\delta _{\mathrm{NNLO}}\) are only a few per mill for \(\eta _{ll} \lesssim 1.5\) and grow to \({{\mathcal {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/\gamma ^*\)-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 signals that the numerical precision of the studied computer programs does not match the formal accuracy of predictions at NNLO.

Fig. 9
figure 9

The pulls for the ATLAS data for \(pp \rightarrow Z/\gamma ^*+X \rightarrow l^+l^- + X\) production at forward rapidities measured at \(\sqrt{s}=7\) TeV [6], normalized to the ABMP16 predictions at NNLO obtained with FEWZ (version 3.1) compared to predictions by the DYNNLO (left) and the FEWZ codes (right). Shown are the LO (dotted), NLO (dashed) and NNLO (solid) predictions for each code

Fig. 10
figure 10

Same as Fig. 9 using predictions by the MATRIX (left) and the MCFM codes (right)

In Fig. 9 we focus on the \(\eta _{ll}\)-distribution for the forward \(Z/\gamma ^*\)-production obtained by ATLAS [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 \(\sigma _{\mathrm{NNLO}}\) with the smallest slicing cuts, \(r^{\mathrm{min}}_{\mathrm{cut}}=0.05\%\) and \(\tau _{\mathrm{cut}}=4\cdot 10^{-4}\). As already remarked above, we use ABMP16 PDFs at NNLO in all cases, independent of the perturbative order. Figures 9 and 10 clearly illustrate the significant corrections up to \({{\mathcal {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 \({{\mathcal {O}}}(5-10\%)\) in some \(\eta _{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.

Fig. 11
figure 11

The DØ data on the electron charge asymmetry distribution \(A_e\) in \(W^{\pm }\)-boson production at \(\sqrt{s}=1.96\) TeV with the statistical (inner bar) and the total uncertainties, including the systematic ones. Shown is the difference of \(A_e\) to the ABMP16 central predictions at NNLO obtained with FEWZ. The symmetric \(p_T^{e,\nu }\)-cuts of the decay leptons are indicated in the figure. The LO (dotted), NLO (dashed-dotted) and NNLO (dashed) predictions by the DYNNLO code (left) and by the MATRIX code (right) are displayed for comparison

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^{\mathrm{min}}_{\mathrm{cut}}=0.05\%\) compared to the FEWZ ones is small in the entire \(\eta _{ll}\) range except for the first \(\eta _{ll}\) bins and originates from the different NNLO cross sections as illustrated in Fig. 6. In contrast, the cross sections \(\sigma _{\mathrm{LO}}\), \(\sigma _{\mathrm{NLO}}\) and \(\sigma _{\mathrm{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^\pm \)-boson production at \(\sqrt{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 \(\alpha _s^{(n_f=5)}(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^\pm \)-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^{\mathrm{min}}_{\mathrm{cut}}=0.15\%\)) being mostly above the FEWZ numbers. The deviations increase with increasing electron pseudo-rapidity \(\eta _e\) and become significant for \(\eta _e \gtrsim 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 Sect. 3 to account for \(r^{\mathrm{min}}_{\mathrm{cut}}\) dependence has no effect.

Fig. 12
figure 12

Same as Fig. 11 using predictions by the MCFM code

Fig. 13
figure 13

Same as Fig. 11, now with staggered cuts \(p_T^{e} > p_T^{\nu }\) on the decay leptons as indicated in the figure for the DYNNLO code (left) and the MATRIX code (right)

Fig. 14
figure 14

Same as Fig. 13 using predictions by the MCFM code

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 \(\tau _{\mathrm{cut}}\) value, \(\tau _{\mathrm{cut}}=6\cdot 10^{-3}\), and the numerical integration accuracy of the individual \(W^\pm \)-boson cross sections is typically . 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.Footnote 10 The parameter cutoff provides the minimum value on any dimensionless variables, for instance, any invariant mass squared \(s_{ij}\) of any pair of partons scaled by their energies, such that \(s_{ij}/(E_i E_j)\) can never be less than cutoff. Cross sections \(\sigma _{\mathrm{NNLO}}\) for \(W^\pm \)-boson production in the DØ kinematics computed with cutoff \(=10^{-9}\) showed severe numerical instabilities. For the agreement of the NNLO MCFM results with FEWZ, we find a similar pattern of increasing deviations with increasing electron pseudo-rapidity \(\eta _e\), which become significant for \(\eta _e \gtrsim 1.0\). The two choices of smaller \(\tau _{\mathrm{cut}}\) values in MCFM, \(1\cdot 10^{-3}\) and \(4\cdot 10^{-4}\), lead to the same NNLO predictions for \(A_e\), within the numerical uncertainites.

Finally, Figs. 13 and 14 show the DØ data and the predictions of DYNNLO, MATRIX and MCFM relative to FEWZ for the electron charge asymmetry distribution \(A_e\) in the case of staggered cuts with \(p_T^{e} \ge 35\) GeV and \(p_T^{\nu } \ge 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.

3 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 \({{\mathcal {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.

3.1 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 \(\tau \), 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

$$\begin{aligned} a(p_a)+b(p_b)&\rightarrow V(q)+X(k_i)&\rightarrow L(q)+X(k_i) \, , \end{aligned}$$
(3)

where ab 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^\mu =(1,0,0,1)\) and \({{\bar{n}}}^\mu =(1,0,0,-1)\), the initial hadron momenta \(p_{a,b}\) can be parametrized in Born kinematics in terms of \(Q=\sqrt{q^2}\) and the rapidity of the gauge boson Y as

$$\begin{aligned} p_a^\mu \,=\, \frac{Q}{2}\,\exp (+Y)\,n^\mu \, , \qquad \qquad p_b^\mu \,=\, \frac{Q}{2}\,\exp (-Y)\,{{\bar{n}}}^\mu \, . \end{aligned}$$
(4)

The slicing parameter \(\tau \) for \(q_T\)-subtraction is defined by

$$\begin{aligned} \tau&=\, q_T^2/Q^2&=\, \left( \sum _i \vec {k}_{T,i}\right) ^2/Q^2 \, , \end{aligned}$$
(5)

in terms of the transverse momenta \(\vec {k}_{T,i}\) of the hadronic real emission radiation. Likewise, one has

$$\begin{aligned} \tau&=\, {{\mathcal {T}}}_0/Q&=\, \sum _i \mathrm{min} \left\{ 2 p_a \cdot k_i, 2 p_b \cdot k_i \right\} /Q^2 \, , \end{aligned}$$
(6)

for the leptonic 0-jettiness \({{\mathcal {T}}}_0\) (see [28, 31] for other variants of 0-jettiness definitions). The leptonic \({{\mathcal {T}}}_0\) is preferred due to smaller intrinsic power corrections, cf. [32] and used in MCFM [20]. The definitions of \(\tau \) 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

$$\begin{aligned} \sigma= & {} \int d\tau \, \frac{d\sigma }{d\tau } \,=\, \int ^{\tau _{\mathrm{cut}}}\!\! d\tau \, \frac{d\sigma }{d\tau } + \int _{\tau _{\mathrm{cut}}}\! d\tau \, \frac{d\sigma }{d\tau } \,=\, \sigma (\tau _{\mathrm{cut}}) \nonumber \\&+ \int _{\tau _{\mathrm{cut}}}\! d\tau \, \frac{d\sigma }{d\tau } \, , \end{aligned}$$
(7)

where \(\tau _{\mathrm{cut}}\) is the cut for the slicing of the phase space. The dependence of \(d\sigma /d\tau \) on \(\tau \) can be predicted from the universal factorization of QCD in soft and collinear limits. It has the structure

$$\begin{aligned} \frac{d\sigma }{d\tau }\sim & {} \delta (\tau ) + \sum _i \left[ \frac{\ln ^i \tau }{\tau }\right] _+ + \sum _j \tau ^{p-1} \ln ^j \tau + {{\mathcal {O}}}(\tau ^{p}) \, , \end{aligned}$$
(8)

where the \(+\)-distributions are the well-known leading threshold logarithms and the terms proportional to \(\tau ^{p-1}\) with \(p>0\) are integrable and denote power corrections in the soft and collinear limit. From the analytical integration one obtains for \(\sigma (\tau _{\mathrm{cut}})\) schematically

$$\begin{aligned} \sigma (\tau _{\mathrm{cut}})\sim & {} 1 + \sum _i \ln ^{i+1} \tau _{\mathrm{cut}} + \sum _j \tau ^{p}_{\mathrm{cut}} \ln ^{j} \tau _{\mathrm{cut}} + {{\mathcal {O}}}(\tau _{\mathrm{cut}}^{p+1}) \, .\nonumber \\ \end{aligned}$$
(9)

The crucial point to stress here is the scaling behavior of the power corrections \(\tau ^{p}_{\mathrm{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 \(\sigma ^{\mathrm{sub}}(\tau _{\mathrm{cut}})\) as

$$\begin{aligned} \sigma= & {} \sigma ^{\mathrm{sub}}(\tau _{\mathrm{cut}}) + \int _{\tau _{\mathrm{cut}}}\! d\tau \, \frac{d\sigma }{d\tau } + \Delta \sigma ^{\mathrm{sub}}(\tau _{\mathrm{cut}}) \, . \end{aligned}$$
(10)

Here the term \(\Delta \sigma ^{\mathrm{sub}}(\tau _{\mathrm{cut}}) = \sigma (\tau _{\mathrm{cut}}) - \sigma ^{\mathrm{sub}}(\tau _{\mathrm{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 \(\sigma ^{\mathrm{sub}}(\tau _{\mathrm{cut}})\) cancels only the leading soft and collinear singularities in \(\sigma (\tau _{\mathrm{cut}})\) in Eq. (9), then the residual power corrections in the presence of cuts on the decay leptons scale as \(\sqrt{\tau _{\mathrm{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 \(\sqrt{{{\mathcal {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_{\mathrm{cut}}=q_T/M_V \rightarrow 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_{\mathrm{cut}} \in [0.15,1]\%\) in steps of \(0.01\%\). It is then recommended to correct the kinematic distributions by rescaling uniformly with the ratio \(\sigma ^{\mathrm{extrapolated}}_{\mathrm{NNLO}}/\sigma ^{r_{\mathrm{cut}}}_{\mathrm{NNLO}}\). MCFM has improved the \(\tau _{\mathrm{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 \(\tau _{\mathrm{cut}}\), cf. Eq. (9). In addition, MCFM computes the cross section for an array of different \(\tau _{\mathrm{cut}}\) values and performs an automated fitting of the \(\tau _{\mathrm{cut}}\) dependence at NNLO with the following ansatz

$$\begin{aligned} \sigma (\tau _{\mathrm{cut}})^{\mathrm{NNLO}}= & {} \sigma _0 + c_1\, \tau _{\mathrm{cut}} \ln ^3(\tau _{\mathrm{cut}}/M_V) \nonumber \\&+ c_2\, \tau _{\mathrm{cut}} \ln ^2(\tau _{\mathrm{cut}}/M_V) + c_3\, \tau _{\mathrm{cut}} \, , \end{aligned}$$
(11)

where \(c_i\) are the fit parameters and the result is then extrapolated to \(\tau _{\mathrm{cut}} \rightarrow 0\). Note, that the functional form in Eq. (11) does not capture well the scaling of the leading power corrections proportional to \(\sqrt{\tau _{\mathrm{cut}}}\) in the case of gauge boson decays [31].

3.2 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

$$\begin{aligned} Z/\gamma ^* \rightarrow L(q)= & {} l_1(p_1)+l_2(p_2)\, ,\nonumber \\ W^\pm \rightarrow L(q)= & {} l^\pm (p_1) + \nu (p_2) \, , \end{aligned}$$
(12)

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-\sum _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

$$\begin{aligned} q^\mu= & {} \left( m_T\,\cosh (Y), q_T, 0, m_T\,\sinh (Y) \right) \, , \nonumber \\ p_1^\mu= & {} p_{T1} \left( \cosh (Y+\Delta y), \cos \phi , \sin \phi , \sinh (Y+\Delta y) \right) \, , \nonumber \\ p_2^\mu= & {} q^\mu - p_1^\mu \, , \end{aligned}$$
(13)

where \(m_T = \sqrt{Q^2+q_T^2}\) and the azimuthal angle \(\phi \) in the transverse plane is given by \(\vec {p}_{T1} \cdot \vec {q}_T = p_{T1}\, q_T\, \cos \phi \). Momentum conservation yields for the transverse momenta and rapidities of the leptons

$$\begin{aligned} p_{T1}= & {} \frac{Q^2/2}{m_T\,\cosh (\Delta y) - q_T \cos \phi } \, , \nonumber \\ p_{T2}= & {} \sqrt{\left( p_{T1}\right) ^2 - 2 p_{T1} q_T\cos \phi + q_T^2} \, , \end{aligned}$$
(14)

and

$$\begin{aligned} \eta _{1}= & {} Y+\Delta y \, , \nonumber \\ \eta _{2}= & {} Y+\frac{1}{2} \ln \left( \frac{m_T - p_{T1} \exp (+\Delta y)}{m_T - p_{T1} \exp (-\Delta y)}\right) \, . \end{aligned}$$
(15)

The leptonic final state phase space \(\Phi _{L}\) (neglecting lepton masses) reads in terms of the variables \(\phi \) and \(\Delta y\) in Eq. (13),

$$\begin{aligned} \Phi _{L}(q_T)= & {} \left( \int \prod _{i=1}^2\, \frac{d^4p_i}{(2\,\pi )^{3}}\, \delta ^+(p_i^2)\,\right) \, (2\,\pi )^4\, \delta ^{(4)}(q-p_1-p_2)\nonumber \\= & {} \frac{1}{4\,\pi ^2}\, \int _0^\pi d\phi \, \int _{-\infty }^{\infty } d\Delta y\, \frac{p_{T1}^2}{Q^2}\, \, , \end{aligned}$$
(16)

and \(p_{T1}\) implicitly depends on \(\phi \) and \(\Delta y\) through Eq. (14). The fiducial cuts on the decay leptons applied to the data discussed in Sect. 2 modify Eq. (16) by constraining the integration range. With the typical cuts on the transverse momenta and rapidities of the leptons, \(p_{T1,2} \ge p_T^{\mathrm{min}}\) and \(\eta _{1,2}^{\mathrm{min}} \le \eta _{1,2} \le \eta _{1,2}^{\mathrm{max}}\), the phase space \(\Phi _{L}\) becomes

$$\begin{aligned} \Phi _{L}(q_T)= & {} \frac{1}{4\,\pi ^2}\, \int _0^\pi d\phi \, \int _{-\infty }^{\infty } d\Delta y\, \frac{p_{T1}^2}{Q^2}\, \nonumber \\&\times \left( \prod _{i=1}^2\, \theta (p_{Ti}-p_T^{\mathrm{min}})\, \theta (\eta _i - \eta _i^{\mathrm{min}} ) \theta (\eta _i^{\mathrm{max}} - \eta _i) \right) \, .\nonumber \\ \end{aligned}$$
(17)

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

$$\begin{aligned} p_{T1}= & {} \frac{Q}{2\cosh (\Delta y)} + q_T \frac{\cos \phi }{2 \cosh ^2(\Delta y)} + {{\mathcal {O}}}(q_T^2/Q) \, , \nonumber \\ p_{T2}= & {} p_{T1} - q_T\cos \phi + {{\mathcal {O}}}(q_T^2/Q) \, , \end{aligned}$$
(18)

and

$$\begin{aligned} \eta _{1}= & {} Y+\Delta y \, , \nonumber \\ \eta _{2}= & {} Y - \Delta y - 2 \frac{q_T}{Q} \cos \phi \sinh (\Delta y) + {{\mathcal {O}}}(q_T^2/Q^2) \, . \end{aligned}$$
(19)

The symmetric cut on the transverse momenta \(p_{T1,2} \ge p_T^{\mathrm{min}}\) used for the ATLAS or DØ data considered above thus splits the \(\phi \) integration range depending on the sign of \(\cos \phi \) as

$$\begin{aligned} {\mathrm{min}}\left\{ p_{T1}, p_{T2}\right\}= & {} \left\{ \begin{array}{lcl} p_{T1}\, , &{} \qquad &{} \cos \phi < 0 \\ \\ p_{T1} - q_T \cos \phi \, , &{} \qquad &{} \cos \phi \ge 0 \end{array} \right. \, , \end{aligned}$$
(20)

such that the \(\theta \)-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_{T1,2} \ge p_T^{\mathrm{min}}\) on the transverse momenta, the phase space \(\Phi _{L}\) becomes

$$\begin{aligned} \Phi _{L}(q_T)= & {} \frac{1}{8\,\pi }\, \sqrt{1-(2 p_T^{\mathrm{min}})^2/Q^2} \nonumber \\&- \frac{1}{2\,\pi ^2}\,\frac{q_T}{Q}\, \frac{p_T^{\mathrm{min}} }{Q\,\sqrt{1-(2 p_T^{\mathrm{min}})^2/Q^2}} + {{\mathcal {O}}}(q_T^2/Q^2) \, .\nonumber \\ \end{aligned}$$
(21)

In case of rapidity cuts, the constraints are slightly more involved, see [34]. Considering the cuts \(|\eta _{1,2}| \le \eta ^{\mathrm{max}}\) for the selection of both leptons at central pseudo-rapidity, as applied to the ATLAS data for \(Z/\gamma ^*\)-boson production, the phase space \(\Phi _{L}\) at Born level becomes

$$\begin{aligned} \Phi _{L}(q_T)= & {} \frac{1}{8\,\pi }\, \tanh (\eta ^{\mathrm{max}}-Y) + {{\mathcal {O}}}(q_T/Q) \, . \end{aligned}$$
(22)

For small Y this constraint is less tight compared to the \(p_{T1,2}\) cuts in Eq. (21) for the typical values of \(p_T^{\mathrm{min}}\) used. To assess where the \(\theta \)-functions in Eq. (17) involving the pseudo-rapidities affect the phase space \(\Phi _{L}\) and, in particular, break azimuthal symmetry, one has to examine the regions where they overlap. The boundary of this region is given by \(|\eta _1| = |\eta _2|\), which determines a value \(\phi ^*\) through the condition

$$\begin{aligned} \cos \phi ^*= & {} \frac{Q}{2 q_T}\, \frac{\sinh (2Y)}{\sinh (2Y+\Delta y)} + {{\mathcal {O}}}(q_T/Q) \, . \end{aligned}$$
(23)

The constraint \(|\cos \phi ^*| \le 1\) restricts either the region of \(|\eta _1|\) or \(|\eta _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 \phi ^*| < 1\) imply the following scaling for the values of \(q_T/Q\)

$$\begin{aligned} \frac{q_T}{Q}> & {} \frac{1}{2}\, \left| \frac{\sinh (2Y)}{\sinh (2Y+\Delta y)} \right| \, , \end{aligned}$$
(24)

which can be approximated for small Y (and using \(|\eta _{1}| = |Y+\Delta y| \le \eta ^{\mathrm{max}}\))

$$\begin{aligned} \frac{q_T}{Q}\gtrsim & {} \frac{q^*_T}{Q} \,=\, \frac{|Y|}{\sinh (\eta ^{\mathrm{max}})} + {{\mathcal {O}}}(Y) \, . \end{aligned}$$
(25)

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.

Fig. 15
figure 15

The difference between the Born and real emission phase spaces \(\Phi _L(0)-\Phi _L(q_T)\) of the decay leptons relative to the Born one at fiducial cuts applied to ATLAS data set [6] for \(Z/\gamma ^*\)-boson production (\(Q=M_Z\)) for different values of the gauge boson pseudo-rapidity \(\eta _{ll}\). For the lepton momenta \(p_T^l \ge 20\) GeV are required. Left: Cuts selecting central pseudo-rapidities \(|\eta _{l_i}| \le 2.5\) for \(i=1,2\). Right: Cuts selecting one lepton at central pseudo-rapidity \(|\eta _{l_1}| \le 2.5\) and the other at forward pseudo-rapidity, \(2.5 \le |\eta _{l_2}| \le 4.9\). The vertical dashed line indicates the minimum value \(r^{\mathrm{min}}_{\mathrm{cut}}=0.15\%\) used in MATRIX as a slicing cut

The appearance of linear power corrections in the phase space \(\Phi _{L}\) can be illustrated by considering the deviations \(\left| 1 - \Phi _{L}(q_T)/\Phi _{L}(0)\right| \) 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/\gamma ^*\)-boson production. On the left in Fig. 15, the leptons are selected at central pseudo-rapidities \(|\eta _{l_i}| \le 2.5\) for \(i=1,2\) and we observe the presence of linear power corrections in \(q_T\) for central gauge boson pseudo-rapidities \(\eta _{ll} \lesssim 1\) due to the \(p_T\) constraint in Eq. (20). For \(r_{\mathrm{cut}}=q_T/Q=0.15\%\), which is the default value for \(r_{\mathrm{cut}}\) used in MATRIX as a slicing cut and indicated by the vertical dashed line in Fig. 15, their size amounts to . In contrast, for larger \(\eta _{ll}\) the pseudo-rapidity constraints dominate the phase space \(\Phi _{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 \(\eta _{ll} = 1.2\) and 1.8, and the corrections to \(\Phi _{L}\) for \(r_{\mathrm{cut}}=0.15\%\) are smaller by more than two orders of magnitude. It is interesting to compare these findings with the \(\eta _{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_{\mathrm{cut}}=0.15\%\)) and MCFM from FEWZ in Figs. 5, 6 and 7 for central \(Z/\gamma ^*\)-boson production. For \(\eta _{ll} \lesssim 1\) all slicing codes undershoot FEWZ, while they tend to agree well for \(\eta _{ll} \gtrsim 1.5\). The transition around \(\eta _{ll} \gtrsim 1.2\) when the linear power corrections in \(q_T\) in \(\Phi _L\) vanish, is most pronounced in the case of MCFM in Fig. 7.

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 \(|\eta _{l_1}| \le 2.5\) and \(2.5 \le |\eta _{l_2}| \le 4.9\) azimuthal symmetry is always broken by the \(p_T\) constraint in Eq. (20) and we observe sizable linear power corrections. For the chosen values of \(\eta _{ll}=1.3\), 2.6 and 3.4 they amount to , and at the value of \(r_{\mathrm{cut}}=0.15\%\) in MATRIX, and the relatively large size of these corrections is remarkable. Moreover, their \(\eta _{ll}\) dependence matches well with the pattern of the observed deviations of DYNNLO from codes with dipole subtraction at NLO in Fig. 3 on the left and of DYNNLO, MATRIX and MCFM from FEWZ at NNLO in Figs. 5, 6 and 7, where the DYNNLO, MATRIX and MCFM deliver significantly smaller results in the first and the last \(\eta _{ll}\) bins.

In Fig. 16 we plot the phase space \(\Phi _{L}\) for the fiducial cuts applied to the ATLAS data set [6] for \(W^\pm \)-boson production. In this case, the binning in the lepton pseudo-rapidity fixes \(\eta _l\) and trivially fulfills one of the \(\theta \)-functions in the integral for \(\Phi _{L}\) in Eq. (17). The other \(\theta \)-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_T^{l,\nu } \ge 25\) GeV for the lepton momenta, which do break azimuthal symmetry as shown in Eq. (20). They amount to corrections for the value \(r_{\mathrm{cut}}=0.15\%\), depending on the value of the lepton pseudo-rapidity \(\eta _l\), with larger corrections observed for central \(\eta _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 \(\eta _l\).

Fig. 16
figure 16

Same as Fig. 15 for the fiducial cuts applied to ATLAS data set [6] for \(W^\pm \)-boson production (\(Q=M_W\)) and different values of the lepton pseudo-rapidity \(\eta _l\). For the lepton momenta \(p_T^{l,\nu } \ge 25\) GeV are required

Fig. 17
figure 17

Same 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_{T1} \ge p_T^{\mathrm{min}}+\delta p_{T}\), \(p_{T2} \ge p_T^{\mathrm{min}}\) and \(Q=M_W\) for the electron pseudo-rapidity fixed at \(\eta _{e}=0\)

Finally, we briefly discuss staggered cuts on the transverse momenta, when \(p_{T1} \ge p_T^{\mathrm{min}}+\delta p_{T}\) and \(p_{T2} \ge p_T^{\mathrm{min}}\) for some \(\delta 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 \(\delta p_{T}\) the phase space \(\Phi _{L}\) evaluates at Born level then as

$$\begin{aligned} \Phi _{L}(q_T)= & {} \frac{1}{8\,\pi }\, \sqrt{1- (2 (p_T^{\mathrm{min}}+\delta p_T))^2/Q^2} + {{\mathcal {O}}}(q_T/Q) \, .\nonumber \\ \end{aligned}$$
(26)

In addition, with staggered cuts the \(\theta \)-functions in Eq. (17) which constrain the fiducial phase space give rise to the following condition for the azimuthal integration

$$\begin{aligned} {\mathrm{min}}\left\{ p_{T1} - \delta p_{T}, p_{T2}\right\}= & {} \left\{ \begin{array}{lcl} p_{T1}\, - \delta p_{T}, &{} \qquad &{} \cos \phi < \delta p_{T}/q_T \\ \\ p_{T1} - q_T \cos \phi \, , &{} \qquad &{} \cos \phi \ge \delta p_{T}/q_T \end{array} \right. \, , \end{aligned}$$
(27)

which reproduces the case of symmetric cuts in Eq. (20) for \(\delta p_{T} \rightarrow 0\). On the other hand, it is obvious that for \(\delta p_{T} > q_T\) the \(\theta \)-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_{\mathrm{cut}}\) or \(\tau _{\mathrm{cut}}\) in the slicing codes imply that \(q_T^{\mathrm{min}} \ll 1\) GeV, while measurements with staggered cuts in the experiments use values of \(\delta 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 \(\delta p_{T} = 0\), 0.1 and 10 GeV for an electron pseudo-rapidity of \(\eta _{e}=0\). Other choices of \(\eta _{e}\) give qualitatively the same results, see Fig. 16. As the value of \(\delta p_{T}\) increases, the deviations from the linear power corrections for the case of \(\delta p_{T} = 0\) GeV become apparent. For \(\delta p_{T}=10\) GeV, which corresponds to the choice in the DØ data selection of \(p_T^{e} \ge 35\) GeV and \(p_T^{\nu } \ge 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_{\mathrm{cut}}\) or \(\tau _{\mathrm{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 \(\Phi _{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 \(\eta _{ll}\) in the case of \(Z/\gamma ^*\)-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_T\)-subtraction schemes and has been applied to provide cross section predictions for Higgs-boson production with fiducial cuts [35].

4 Conclusions

We have investigated the accuracy of available NNLO QCD predictions for the hadroproduction of \(W^{\pm }\)- 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 pseudo-rapidities of the decay leptons lead to linear power corrections in the slicing parameter, i.e. \(q_T/Q\) and \(\sqrt{{{\mathcal {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 \(\Phi _{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 \({{\mathcal {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.