$W^+W^-$ production at NNLO+PS with MiNNLO$_{\rm PS}$

We consider $W^+W^-$ production in hadronic collisions and present the computation of next-to-next-to-leading order accurate predictions consistently matched to parton showers (NNLO+PS) using the MiNNLO$_{\rm PS}$ method. Spin correlations, interferences and off-shell effects are included by calculating the full process $pp \to e^+\nu_e \mu^-\bar{\nu}_\mu$. This is the first NNLO+PS calculation for $W^+W^-$ production that does not require an a-posteriori multi-differential reweighting. The evaluation time of the two-loop contribution has been reduced by more than one order of magnitude through a four-dimensional cubic spline interpolation. We find good agreement with the inclusive and fiducial cross sections measured by ATLAS and CMS. Both NNLO corrections and matching to parton showers are important for an accurate simulation of the $W^+W^-$ signal, and their matching provides the best description of fully exclusive $W^+W^-$ events to date.


Introduction
Precision phenomenology has evolved to one of the cornerstones of todays physics programme at the Large Hadron Collider (LHC). Without clear hints for new physics, the precise measurement of production rates and distributions of Standard Model (SM) processes provides a valuable path towards the observation of deviations from the SM picture. The production of vector-boson pairs is among the most important LHC signatures in that respect. Those processes are crucial to constrain or measure anomalous interactions among SM particles, such as anomalous couplings among three vector bosons (triple-gauge couplings), as any small deviation from the expected rates or shapes of distributions could be a signal of new physics. W + W − production has the largest cross section among the massive diboson processes and it provides direct access to triple-gauge couplings, which appear already in the leading perturbative contribution to the cross section. The measurement of this process at the LHC is a direct probe of the gauge symmetry structure of electroweak (EW) interactions and of the mechanism of EW symmetry breaking in the SM. Moreover, W + W − final states are an irreducible background to Higgs measurements in the H → W + W − decay channel and to direct searches for BSM particles decaying into two leptons, missing energy, and/or jets. The W + W − cross section has been measured at both the Tevatron [1][2][3] and the LHC (at 7 TeV [4-7], 8 TeV [8][9][10][11] and 13 TeV [12][13][14][15]). The high sensitivity to anomalous triplegauge couplings has been exploited in various indirect BSM searches [1, 2, 4-8, 10, 11, 16-19] and the irreducible W + W − background has been extensively studied in the context of H → W + W − decays in refs. [20][21][22][23][24][25][26][27][28][29][30][31][32].
The theoretical description of fiducial cross sections and kinematic distributions has been greatly improved by the calculation of next-to-next-to-leading order (NNLO) corrections in QCD perturbation theory, which have become the standard for 2 → 1 and 2 → 2 colour-singlet production . With γγγ production even the first 2 → 3 LHC process was recently pushed to NNLO accuracy [66,67]. In comparison to LHC measurements NNLO corrections are crucial for a more accurate and precise description of data. On the other hand, the validity of fixed-order calculations is challenged in kinematical regimes sensitive to soft and collinear radiation through the appearance of large logarithmic contributions. In such regimes an all-order description is mandatory to obtain physically meaningful predictions. The analytic resummation of large logarithmic contributions is usually restricted to a single observable or at most two observables, see e.g. ref. [68] for the recent next-to-next-to-next-to-logarithmic (N 3 LL) result of the W + W − transversemomentum (p T,W W ) spectrum and the joint resummation of logarithms in p T,W W and in the transverse momentum of the leading jet (p T,j 1 ) at next-to-next-to-logarithmic (NNLL) accuracy. By contrast, parton showers are based on a numerical resummation approach with limited logarithmic accuracy, but they include all-order effects in all regions of phase space at the same time. Moreover, the fully exclusive description of the final state enables a full-fledged hadron-level simulation that is indispensable for experimental analyses.
In order to meet the experimental demands for having both high-precision predictions and exclusive hadron-level events, an enormous effort is made by the theory community to include higher-order corrections in parton showers. Almost two decades ago the matching of next-to-leading order (NLO) QCD predictions and parton showers (NLO+PS) was formulated in seminal publications [69][70][71]. More recently, the first NNLO+PS approaches have been developed for colour-singlet processes [72][73][74][75][76], and the MiNNLO PS approach of refs. [75,76] was very recently extended to heavy-quark pair production [77]. The methods of refs. [72,75,76] originate from the MiNLO procedure [72,78], which upgrades a NLO calculation for colour singlet plus jet production to become NLO accurate for both zero-jet and one-jet observables by exploiting features of the all-order structure of the transverse momentum resummation formula. In ref. [79], a numerical extension of the MiNLO procedure to higher jet multiplicities was presented and applied to Higgs production in association with up to two jets. Most NNLO+PS applications have been done for simple 2 → 1 LHC processes or 1 → 2 decays so far, such as Higgs-boson production [75,76,80,81], Drell-Yan (DY) production [74-76, 82, 83], Higgsstrahlung [84][85][86], which is still a 2 → 1 process with respect to QCD corrections, and the H → bb decay [87,88]. There are a few notable exceptions where NNLO+PS matching was achieved for more involved colour-singlet processes, namely W + W − [89], Zγ [90], γγ [91] and ZZ [92] production. Moreover, with top-quark pair production the very first NNLO+PS calculation for a coloured initial and final state has been presented in ref. [77].
In the case of W + W − production at the LHC, substantial advancements have been made in the theoretical description of the process in terms of both fixed-order and allorder calculations. W -boson pairs are produced in quark annihilation at LO, which was calculated several decades ago for on-shell W bosons [93]. NLO QCD corrections were obtained in the on-shell approximation first [94,95], and in refs. [96][97][98][99] the leptonic W decays with off-shell effects and spin correlations were accounted for. Also, NLO EW corrections are known both for on-shell W bosons [100][101][102] and including their off-shell treatment [103][104][105]. The simplest O(α 2 s ) contribution is the loop-induced gluon fusion channel. Being separately finite and enhanced by the large gluon luminosities, its LO cross section is known already for a long time [99,[106][107][108][109][110][111][112][113][114][115]. The full NNLO QCD corrections were first obtained for the inclusive cross section in the on-shell approximation [59], while the fully differential NNLO calculation for off-shell W bosons was presented in ref. [60], using the qq → V V two-loop helicity amplitudes [116][117][118]. Recently, NNLO corrections were studied for polarized W + W − production [119]. Also NLO QCD corrections to the loop-induced gluon fusion contribution, which are formally of O(α 3 s ), were evaluated using the gg → V V two-loop helicity amplitudes of refs. [120,121]: first in an approximation without quark initial states [122] and later including all relevant contributions [123]. To date the most advanced fixed-order prediction for W + W − production combines all of those contributions and is available in the Matrix framework [124]: the combination of NNLO QCD [60] and NLO EW predictions has been achieved in ref. [105] using Matrix and OpenLoops [125][126][127]. Approximate N 3 LO predictions (labelled as nNNLO) have been calculated by combining the NNLO quark-initiated cross section with the NLO gluoninitiated cross section in ref. [123], where the nNNLO cross section has also been combined with NLO EW corrections.
All-order predictions for the W + W − process have been obtained for various observables using state-of-the-art resummation techniques: threshold resummation at NLO+NNLL was presented in ref. [128], b-space resummation was used to obtain the NNLO+NNLL transverse momentum spectrum of the W + W − pair [129] and the NNLO+NNLL jet-vetoed cross section was computed in ref. [130]. More recently, the Matrix+RadISH framework was introduced [68,131,132], which combines NNLO calculations in Matrix with highaccuracy resummation through the RadISH formalism [133][134][135]. For all 2 → 1 and 2 → 2 colour-singlet processes the Matrix+RadISH code makes NNLO+N 3 LL predictions for the transverse momentum of the colour singlet, NNLO+NNLL predictions for the transverse momentum of the leading jet, as well as their joint resummation at NNLO+NNLL publicly available. In particular, ref. [68] has applied this resummation framework as an example to W + W − production, presenting state-of-the-art predictions for the p T,W W spectrum, the p T,j 1 spectrum, the jet-vetoed cross section and the p T,W W spectrum with a jet veto. Indeed, one important aspect of the theoretical description of W + W − production is the correct modelling of the jet veto (see refs. [68,130,[136][137][138][139] for example), which is applied by the experimental analyses to suppress backgrounds involving top-quarks (tt and tW ). A strict veto against jets in the final state increases the sensitivity to higher-order QCD effects due to potentially large logarithms of the ratio of the small jet-veto scale over the large invariant mass of the system. Such terms challenge the reliability of fixed-order predictions and induce large uncertainties in theory predictions that are typically not covered by scale-variation procedures, especially when extrapolating cross-sections measured in the fiducial region to the total phase space. In particular, the tension with NLO+PS predictions observed in earlier W + W − measurements [9, 140] challenged the validity of lower-order Monte Carlo predictions for W + W − production [139]. Only through the calculation of NNLO corrections [59,60] this tension could be released, and their combination with all-order resummation confirmed that the jet-vetoed W + W − cross section is under good theoretical control [68,130]. Moreover, it was shown that resummation effects are eventually required to obtain reliable predictions in the tails of some kinematical distributions, for instance in the invariant mass distribution of the W + W − pair [141] when a jet-veto is imposed. These issues show the relevance of fully flexible, hadron-level Monte Carlo predictions with state-of-the-art perturbative precision for the W + W − production process, which is achieved by the the combination of NNLO corrections with parton-shower simulations.
In this paper, we obtain NNLO+PS predictions for W + W − production using the MiNNLO PS method. For the first time NNLO QCD corrections are directly included during the generation of W + W − events, without any post-processing or reweighting being required. In fact, this is also the first time a NNLO W + W − calculation independent of a slicing cutoff is performed (cf. refs. [59,60]). To this end, we have applied the recently developed MiNNLO PS method [75,76] and its extension to 2 → 2 reactions presented in ref. [90]. At variance with the NNLOPS calculation of ref. [89], our new MiNNLO PS generator does not include any of the approximations or limitations related to the reweighting approach used in ref. [89]. In particular, ref. [89] had to resort to a number of features of the W -boson decays, such as the fact that the full angular dependence of each vector-boson decay can be parametrized through eight spherical harmonic functions [156] and the fact that QCD corrections are largely independent of the off-shellness of the vector bosons, in order to simplify the parametrization of the nine dimensional W + W − → e + ν e µ −ν µ Born phase space. Moreover, the discretization of the residual variables in the parametrization of the Born phase space for the reweighting limits the numerical accuracy in phase-space regions sensitive to coarse bins. Not rarely, such regions can be relevant for BSM searches, especially when situated in the tails of kinematic distributions. Without those limitations, our new MiNNLO PS calculation provides the most flexible and most general simulation of W + W − signal events with NNLO accuracy at the LHC. For the two-loop contribution, we use the helicity amplitudes for the production of a pair of off-shell vector bosons [118] from the public code VVAMP [157] and exploit their implementation for all qq → 4 leptons processes in the Matrix framework [124,158]. The evaluation of these two-loop amplitudes turns out to be the major bottleneck in our calculation. In order to deal with this we substantially speed up the evaluation time by using a four-dimensional cubic spline interpolation procedure of the two-loop coefficients entering the helicity amplitudes.
In the present calculation we consider all topologies that lead to two opposite-charge leptons and two neutrinos in the final state ( + ν −ν ) with off-shell effects, interferences, and spin correlations. As a basis we exploit the W + W − +jet generator of ref. [153] and include NNLO QCD corrections to W + W − production through the MiNNLO PS method. The ensuing MiNNLO PS generator is implemented and will be made publicly available within the Powheg-Box-Res framework [70,71,152,159], which provides a general interface to parton showers. This is necessary for a complete and realistic event simulation. Especially, non-perturbative QCD effects using hadronization and underlying event models, as well as multiple photon emissions through a QED shower can be included. Those can induce sizable corrections in jet-binned cross sections, on the lepton momenta (especially invariant mass distributions/line shapes), and other more exclusive observables measured at the LHC. In our calculation and throughout this paper we omit the loop-induced gluonfusion contribution, as it is already known to higher-order in QCD [122,123] and can be evaluated with known tools at LO+PS, such as the gg2ww generator [160,161] used by ATLAS and CMS. In fact, also a NLO+PS generator was presented for this process recently [162] in the Powheg-Box-Res framework. Finally, we define W + W − signal events free of top-quark contamination by exploiting the four-flavour scheme with massive bottom quarks and drop all contributions with final-state bottom quarks. Refs. [59,60] have shown for both total and fiducial rates at NNLO that this approach agrees within ∼1-2% with an alternative procedure to obtain top-free W + W − predictions. The latter one is defined in the five-flavour scheme and exploits the resonance structure of top-quark contributions to extract the part of the cross section independent of the top-quark width.
This manuscript is organized as follows: in section 2 we provide all details about our calculation and implementation. In particular, we introduce the process and its resonance structures (section 2.1), describe the MiNNLO PS formulae (section 2.2) and the practical implementation in Powheg-Box-Res+Matrix (section 2.3). We also discuss in detail how we obtain the full two-loop contributions by interpolating the basic two-loop coefficients entering the helicity-amplitudes and how we validated this procedure (section 2.4). In section 3, after describing the setup and the set of fiducial cuts used in the analysis (section 3.1), we present phenomenological results for MiNNLO PS and compare them against MiNLO , NNLOPS, NNLO, analytic resummation, and data for both integrated cross sections (section 3.2) and differential observables (section 3.3). We conclude and summarize in section 4. Figure 1: Sample LO diagrams in the different-flavour channel + ν −ν for (a) t-channel W + W − production, (b) s-channel Z/γ → W + W − production, and (c) DY-type production.
2 Outline of the calculation

Description of the process
We study the process for any combination of massless leptons , ∈ {e, µ, τ } with different flavours = . For simplicity and without loss of generality we consider only the process pp → e + ν e µ −ν µ + X here, which we will refer to as W + W − production in the following. By including all resonant and non-resonant topologies leading to this process, off-shell effects, interferences and spin correlations are taken into account. Sample LO diagrams are shown in figure 1, including (a) double-resonant t-channel W + W − production, The corresponding production of opposite-charge same-flavour leptons pp → + ν −ν + X involves the same type of W + W − diagrams as shown in figure 1, but also additional ZZ diagrams as shown in figure 2. By focusing on the different-flavour case ( = ) we avoid the complications originating from the mixing of the W + W − and ZZ topologies. In fact, as shown in refs. [58,104,147], W + W − and ZZ interference effects can be largely neglected and, to a very good approximation, predictions for the two processes can be added incoherently.
An important aspect of W + W − production is that its cross section is subject to a severe contamination from top-quark contributions. Not only does this affect W + W − measurements at the LHC, which usually employ a jet veto, a b-jet veto, or both to suppress top-quark backgrounds, it also renders the theoretical definition of the W + W − cross section cumbersome. Indeed, resonant top-quark contributions enter radiative corrections to W + W − production through interference with real-emission diagrams involving two bottom quarks in the final state. Those interference terms are numerically so large that they easily provide the dominant contribution to the cross section. Specifically, in the inclusive phase space genuine W + W − contributions are more than one order of magnitude smaller. Therefore, the consistent removal of the top-quark contamination is mandatory to define a top-free W + W − cross section. To this end, we exploit the four-flavour scheme (4FS), where bottom quarks are treated as being massive, do not enter in the initial state and diagrams with real bottom-quark radiation are separately finite. This allows us to drop all contributions with final-state bottom quarks, thereby cancelling the top-quark contamination and obtaining top-free W + W − results. We note that there exists an alternative approach to define a top-free W + W − cross section that can be used in the five-flavour scheme (5FS). However, this approach is much less practical as it requires the repeated evaluation of the cross section (and distributions) with increasingly small values of the top-quark width Γ t to extract the top-free W + W − cross section as the contribution that is not enhanced by 1/Γ t . Indeed, it was shown in ref. [59] at the inclusive level and in ref. [60] for the fully-differential case that the 4FS and the 5FS definition of the W + W − cross section agree at the level of ∼1-2%. For the sake of simplicity, the easier 4FS approach is employed throughout this paper.

The MiNNLO PS method
We employ the MiNNLO PS method to build a NNLO+PS generator for W + W − production. The method was introduced in ref. [75], optimized for 2 → 1 scattering processes in ref. [76], and generalized to 2 → 2 colour-singlet scattering processes in ref. [90]. In the following we recall the basic ideas and essential ingredients of MiNNLO PS , adapting the notation of refs. [75,76,90].
MiNNLO PS formulates a NNLO calculation fully differential in the phase space Φ F of a produced colour singlet F with invariant mass Q, in such a way that it can be subsequently matched to a parton shower. It starts from a differential description of colour singlet plus jet (FJ) production in the Powheg approach [70,71,152] and it achieves NNLO accuracy for F production by modifying the content of theB(Φ FJ ) function. With Φ FJ we have denoted the FJ phase space, ∆ pwg is the Powheg Sudakov form factor, Φ rad (p T,rad ) is the phase space (transverse momentum) of the second-hardest radiation, and B and R denote the squared tree-level matrix elements for FJ and FJJ production, respectively. The central ingredient of the MiNNLO PS method is the modified B(Φ FJ ) function, which describes the F process at NNLO and the FJ process at NLO, including both zero and one QCD emissions, respectively. The content of the curly brackets generates the second QCD emission according to the Powheg mechanism, with a default Powheg cutoff of Λ pwg = 0.89 GeV. Additional radiation that contributes at O(α 3 s (Q)) and beyond to all orders in perturbation theory is added by the parton shower.
The MiNNLO PSB (Φ FJ ) function can be expressed as follows [75,76,90] where p T refers to the transverse momentum of the color singlet. The overall sum runs over all flavour structures FJ of FJ production, while F denotes the flavour structures of the Born process pp → F. With F ← FJ we denote a projection of the flavour structures, which is trivial in the case of W + W − production, since the Born is always qq initiated. All quantities with index F have to be evaluated in the Born kinematics Φ F , which requires a suitable projection Φ FJ → Φ F as introduced in appendix A of ref. [75]. The notation [X] (i) is used for the i-th term in the perturbative expansion of a quantity X.S F (p T ) represents the Sudakov form factor and dσ FJ is the differential fixed-order cross section, as defined in eqs. (2.9) and (2.11) of ref. [75], respectively. The last term in eq. (2.3) is the central contribution added by the MiNNLO PS method to achieve NNLO accuracy. The precise definition and derivation of D F (p T ) is discussed below. The factor F corr (Φ FJ ) encodes the dependence of the Born-like NNLO corrections upon the full Φ FJ phase space, as discussed in detail in section 3 of ref. [75] and section 3.3 of ref. [90]. A few comments are in order: a crucial feature of the MiNNLO PS method is that the renormalisation and factorisation scales are evaluated as µ R ∼ µ F ∼ p T . As a consequence, each term contributes to the total cross section with scales µ R ∼ µ F ∼ Q according to the following power counting formula: This implies that, when including terms up to second order in α s (p T ) in eq. (2.3), upon integration over p T , the cross section is NLO accurate, as observed first in ref. [72]. By deriving also all (singular) contributions in eq. (2.3) at third order in α s (p T ), NNLO accuracy is achieved after integration over p T [75]. Indeed, D F (p T ) consistently adds the relevant singular α 3 s (p T ) corrections, while regular contributions at this order can be safely omitted as a consequence of the counting in eq. (2.4). In fact, two results for D F (p T ) have been derived [75,76] that differ only by terms of O(α 4 s ) and higher. Their derivation stems from the analytic formulation of the NNLO cross section differential in p T and Φ F : where R f includes only non-singular contributions at small p T , and The luminosity factor L F (p T ) contains the parton densities, the squared hard-virtual matrix elements for F production up to two loops as well as the NNLO collinear coefficient functions, and its expression is given in eq. (3.5) of ref. [90]. As discussed in detail in ref. [75], by choosing a suitable resummation scheme (µ R ∼ µ F ∼ p T ) and matching scheme (factoring outS F (p T ) from R f as well), and by making eq. (2.5) accurate to third order in α s (p T ), the relevant corrections to achieve NNLO accuracy upon integration over p T are derived. In the original MiNNLO PS formulation of ref. [75] the expansion was truncated beyond third order in α s (p T ), so that D F (p T ) would be derived as which breaks the total derivative of the starting formula in eq. (2.5). Instead, ref. [76] suggested a new prescription that preserves the total derivative by keeping into account additional terms beyond accuracy, so that we use as our default choice throughout this paper. The relevant expressions for its evaluation, including the ones of the [D F (p T )] (i) coefficients, are reported in appendix C and D of ref. [75] and in appendix A of ref. [76], where the flavour dependence can be simply included through the replacements F . We further note that the flavour dependence ofS F (p T ) and L F originates entirely from the hard-virtual coefficient function H F , which for a general 2 → 2 hadronic process depends on both the flavour and the Born phase space. This dependence propagates to the Sudakov form factor through the replacementB F and H (2) F are unambiguously defined in section (3.3) of ref. [90].

Practical implementation in Powheg-Box-Res+Matrix
As a starting point, we exploit the W + W − +jet generator developed in ref. [153] for Powheg-Box-V2 [152] and integrated it into the Powheg-Box-Res framework [159]. To this end, we had to adapt the Powheg-Box-Res code to automatically find all relevant resonance histories for W + W − +jet production. This was required, because the automatic generation of resonance histories is not fully functional for processes with a jet in the final state. As described in detail in ref. [159] and recalled in section 2.2 of ref. [90], the correct implementation of all resonance histories is necessary to take advantage of the efficient phase-space sampling within Powheg-Box-Res. We have then upgraded the W − W + +jet generator to include NNLO accuracy for W + W − production by means of the MiNNLO PS method. This has been achieved by making use of the general MiNNLO PS implementation for colour singlet production developed in ref. [90] and adapting it consistently to the 4FS.
As far as the physical amplitudes are concerned, all tree-level real and double-real matrix elements (i.e. for + ν −ν +1,2-jet production) are evaluated through the Powheg-Box interface to Madgraph 4 [163] developed in ref. [164]. The + ν −ν +jet one-loop amplitude is obtained from GoSam 2.0 [165], neglecting one-loop fermion box diagrams, which have been shown to give a negligibly contribution, but slow down the code substantially (cf. ref. [153]). 1 The Born-level and one-loop amplitudes for + ν −ν production have been extracted from MCFM [166]. The (one-loop and) two-loop qq → + ν −ν helicity amplitudes that were derived in ref. [118] are obtained through their implementation in Matrix by suitably adapting the interface created in ref. [90]. Those amplitudes are known only in the massless approximation, but the effect of including massive quark loops is expected to be negligible because of the smallness of closed fermion-loop contributions. For a fast evaluation of the two-loop amplitudes, we have generated interpolation grids, as discussed in detail in the next section.
The calculation of D F (p T ) in eq. (2.8) involves the evaluation of several convolutions with the parton distribution functions (PDFs), which are performed through hoppet [167]. Moreover, the collinear coefficient functions require the computation of polylogarithms, for which we employ the hplog package [168].
Finally, we report some of the most relevant (non-standard) settings we have used to produce W + W − events. We refer the reader to ref. [76] for a detailed discussion on these settings. In particular, to avoid spurious contributions from higher-order logarithmic terms at large p T we consistently introduce modified logarithms with the choice of p = 6, as defined in eq. (10) of ref. [76]. At small p T , we use the standard MiNNLO PS scale setting in eq. (14) of ref. [76], while we activate the option largeptscales 1 to set the scales entering the NLO W + W − +jet cross section at large p T as in eq. (19) of ref. [76]. We use those scale settings with the parameter Q 0 = 0 GeV, and instead regularize the Landau singularity by freezing the strong coupling and the PDFs for scales below 0.8 GeV. We turn on the Powheg-Box option doublefsr 1, which was introduced and discussed in detail in ref. [169]. As far as the parton-shower settings are concerned, we have used the standard ones (also for the recoil scheme).

Fast evaluation of the two-loop amplitude
As discussed before, the two-loop helicity amplitudes for the production of a pair of offshell vector bosons were computed in ref. [118] and the relevant coefficients functions to construct the amplitudes can be obtained from the publicly available code VVAMP [157]. Using those results all qq → 4 leptons amplitudes have been implemented in the Matrix framework [124,158]. To exploit this implementation for our calculation, we have compiled Matrix as a C++ library and linked it to our MiNNLO PS generator using the interface created in ref. [90].
The evaluation of these two-loop amplitudes turns out to be the bottleneck of the calculation. In fact, it takes on averaget VVAMP ≈ 1.9 s to evaluate a single phase-space point, while the evaluation of the tree-and one-loop amplitudes are orders of magnitude faster. Therefore, even though we provide the option to run the code using the exact twoloop amplitudes, all of the results of this paper have been obtained using a four-dimensional cubic spline interpolation procedure for the set of independent two-loop coefficient functions that are required for the evaluation of the two-loop helicity amplitudes. In the following, we present this procedure in detail.

Coefficient functions of the qq → + ν −ν helicity amplitudes
We start by recalling some relevant formulae in ref. [118] for the helicity amplitudes. Specifically, the physical process is denoted by: where p i are the momenta of the corresponding particles and each of the two off-shell W bosons decays into a neutrino-lepton pair, such that p 3 = p 5 + p 6 and p 4 = p 7 + p 8 . We denote by M λλ 1 λ 2 the bare helicity amplitudes of a general vector-boson pair production process, where λ represents the handedness of the partonic current, while λ 1 and λ 2 stand for the helicities of the two leptonic currents. There are just two independent helicity amplitudes M LLL and M RLL , since all the other helicity configurations can be recovered by permutations of external legs [118]. The bare helicity amplitudes are the building blocks of the dressed helicity amplitudes M λLL , which are process specific and for W + W − production read where λ = L, R. In the previous expression, α EW refers to the EW coupling constant, θ W to the mixing angle, and m W and Γ W to the W -boson mass and decay width, respectively. Since a W boson can just couple to left-handed lepton currents, it is clear that M λRL = M λLR = M λRR = 0. As shown in ref. [118], for four-dimensional external states the expression of the bare helicity amplitudes can be written in a compact form using the spinor-helicity formalism: where the two indices i and j are determined by the handedness of the partonic current: Eq. (2.11) depends on nine complex scalar coefficients E j , which are functions of the invariant masses p 2 3 and p 2 4 of the two vector bosons and of the two Mandelstam invariantss andt, defined as (2.12) Each coefficient E j receives a contribution from four different classes of diagrams C where i 1 , i 2 represent the colours of the incoming quark and anti-quark, respectively, and Q λ,W + W − ,[C] q q denotes a coupling factor, which is the only process specific ingredient entering eq. (2.13). Following the labeling introduced in ref. [118] for the diagram classes, we have for W + W − production: • class A and B, including all diagrams where the two vector bosons are attached to the fermion line, with the W + boson adjacent to the incoming quark or antiquark, respectively, whose coupling factors read which is identical to zero for λ = R; • class C, containing diagrams where both vector bosons are attached to a fermion loop, where with n g being the number of massless quark generations; • class F V , collecting form-factor diagrams where the production of the two W bosons is mediated either by a virtual photon (V = γ * ) or a Z boson (V = Z), as shown in figure 1b. 2 In that case we have where e q and I 3 q are the electric charge and isospin number of the incoming quark q, and m Z and Γ Z the Z-boson mass and decay width, respectively.
Since the functions E j admit a perturbative expansion as . Moreover, any loop correction to the corresponding form-factor diagrams just amounts to a function F(s) which multiplies the tree level structure, so that at two loops (2.18) The tree-level coefficients evaluate to constants: The dependence ons in F (2) (s) just enters through the ratio ofs with the squared of the renormalization scale µ R . By setting µ R = √s , the non-vanishing E [F ],(2) j coefficients also become constants, which are known [170][171][172][173]. Note that the correct renormalization scale dependence will be recovered through the MiNNLO PS formulae (cf. appendix D of ref. [75]). Therefore, only the coefficient functions belonging to families C = {A, B, C} need to be interpolated. This reduces the number of real-valued functions that need to be interpolated from 90 to 54.

Generation of interpolation grids
As a first step, we have generated rectilinear grids (i.e comprised of congruent parallelotopes) for each of the 54 non-constant two-loop coefficient functions E [C],(2) j defined in eqs. (2.13) and (2.17), whose exact values have been computed through VVAMP on a set of given phasespace points (p 2 3 , p 2 4 ,s,t) and stored. All results have been obtained by fixing the centreof-mass energy to √ s = 13 TeV. As it turns out, a suitable parametrization of the four-dimensional phase-space points (p 2 3 , p 2 4 ,s,t) is crucial to obtain a good interpolation performance. Moreover, a finer binning is required in those phase-space regions that have a large contribution to the overall integral of the multi-differential cross section, such as resonance-enhanced regions in p 2 3 and p 2 4 around the two W -boson masses. To this end, our grids are defined on a four-dimensional unit hypercube [0, 1] 4 with fifty equally spaced bins, where each element ( is uniquely mapped to a physical phase-space point. The first two axes x 1 and x 2 are related to the invariant masses p 2 3 and p 2 4 through the transformations where z 1 and z 2 are continuous functions of x 1 and x 2 , respectively. The lower and upper bounds on z 1 and z 2 have been chosen in such a way that the physical range of invariant mass values is covered. Specifically, z 1, min and z 1, max are fixed by the choice 40 GeV 2 < p 2 3 < s/100. Through energy conservation z 2, min and z 2, max depend directly on p 2 3 , as it has been made explicit in eq. (2.20). However, their exact expressions, which we omit here, have been tuned such that the physical mass range of p 2 4 is covered efficiently. The two functions z 1 (x 1 ) and z 2 (x 2 ) are defined piecewise on three subranges of the two intervals in eq. (2.20). In the central subrange, z 1 and z 2 correspond to a linear mapping, which guarantees that p 2 3 and p 2 4 follow a Breit-Wigner distribution. For the other two subranges of z 1 and z 2 polynomial functions are used such that the off-shell regions are covered by a sufficient number of grid points.
The variables x 3 and x 4 also have a physical interpretation, since they are related to the relativistic velocity β W + and the cosine of the scattering angle cos θ W + of one of the vector bosons in the center of mass frame. In particular, we define where a s/t and b s/t determine the range of values allowed for the two physical quantities. Instead of setting a s/t = 0 and b s/t = 1, we use small cutoffs to avoid numerical instabilities at the kinematic edges. β W + and cos θ W + can be expressed in terms ofs andt as we can expresss andt in terms of the hypercube variables x 3 and x 4 . As illustrated in ref. [118], the behaviour of the coefficients E is not always smooth over the two-dimensional phase space (β W + , cos θ W + ) and it can even be divergent close to the highly relativistic (β W + → 1) or highly collinear (| cos θ W + | → 1) regions. One possibility to improve the description of this rapidly changing functional behaviour is to combine different grids to cover the whole phase space (p 2 3 , p 2 3 ,s,t), instead of simply increasing the number of bins for selected axes. For the case at hand, using four precomputed grids for each of the 54 real-valued functions has proven to significantly improve the performance of the interpolator in some phase-space regions. Even though the definition of the grids is unchanged for x 1 and x 2 , by adjusting the values of a s/t and b s/t in eq. (2.21) we defined four slightly overlapping grids in the (β W + , cos θ W + ) phase space, to properly cover the above mentioned singular regions.

Interpolation and validation
At the beginning of each run the grids just need to be read and loaded into memory. Then, for each E [C],(2) j any value can be computed by properly interpolating between the values stored in the precomputed grids. To perform this task, we make use of the N -dimensional interpolation library Btwxt [174], which just requires the input grids to be rectilinear. The interpolation is achieved through N -dimensional cubic splines [175], which are multivariate piecewise polynomials of degree three. Specifically, Btwxt employs cubic Hermite splines, where each polynomial in a given N -dimensional interval is specified by its values and its first derivatives at the corners of the interval itself. The values of the first derivatives are computed according to the Catmull-Rom implementation [176].
In order to quantify the accuracy of our four-dimensional interpolation strategy we use an adimensional parameter , which describes the deviation of the interpolated result for the two-loop contribution from its exact expression on a given phase-space point and is defined as where the four independent Born flavour configurations F = {qq,qq; for q = u-type or q = d-type} have been considered separately. In eq. (2.25), H qT (2) (ex) F refers to the hard function in the q T -scheme [177] returned by Matrix using the evaluation of the exact two-loop coefficients through VVAMP, while H qT(2) (int) F stands for its value obtained by Matrix using the interpolation of the two-loop coefficients from the precomputed grid results. Note that the conversion between the q T -scheme H qT(2) (int) F and the MiNNLO PS H (2) F in section 2.2 has been given in eq. (3.22) of ref. [90].
In figure 3, the distribution of the values of F is displayed for a selected flavour channel (namely theūu one). All the other flavour channels have the same qualitative behaviour. The figure shows the impact of increasing the number of precomputed grids on the F parameter from one (orange histogram) to four (blue histogram). Besides being essential for a simultaneously accurate description of physical observables over a wide phase-space region, our choice of covering the phase space with four separate grids improves the accuracy of the predictions for the two-loop contribution and yields a more symmetric F distribution. We further notice that the bulk of the interpolator predictions (roughly 80%) has an accuracy greater than 5% (i.e | F | ≤ 5%), while almost 95% lie inside the interval | F | ≤ 100%. The remaining fraction of F values consists of phase-space points where the interpolator poorly reproduces the correct two-loop result. In figure 3, where the F distribution is reported both in linear and in logarithmic scale, this fraction is clearly visible in the overflow bins at the edges of the histograms. In most cases, these poorly predicted values are associated to phase-space points falling outside the grid boundaries and thus requiring extrapolation of the two-loop coefficient functions E [C],(2) j outside the grid edges. However, this also means that most of these points lie in kinematical regions where the cross section is strongly suppressed.
To deal with instabilities, a basic rescue-system is introduced. This mechanism takes care of computing the exact E (ex) F distribution is concentrated. In figure 4 this distribution is shown together with the median and the value of the first three moments of the distribution. As it can be inferred from the positive skewness value (or from the fact that the mean and median do not coincide), the distribution is asymmetric, which is why our acceptance interval for H qT (2) (ex) F is not centered around the mean, but rather it extends to higher values to partially account for the long distribution tail on the right of the peak. This simple  criterium suffices to catch the small fraction of F outliers that would compromise the stability of the results. Some phase-space points remain that elude the rescue-system and where the two loop coefficients are not computed accurately, but we have checked that they have a negligible impact on the physical results, as it will be discussed below (see figure 6). Clearly, the advantage of using the interpolation approach compared to the full evaluation of the two-loop coefficients is the time performance. Indeed, the average time required by VVAMP (t VVAMP ) and the interpolator (t int ) to evaluate the two-loop contribution turns out to differ by three orders of magnitude, while the improvement is still roughly a factor forty if one uses the rescue-system (t int ): As complementary information, figure 5 shows the time distributions of t VVAMP (blue histogram) and t int (orange histogram). The bulk of the VVAMP evaluation times (roughly 80%) is located in the time interval 0.5s < t VVAMP < 2.0 s, with a small, but not negligible fraction of phase-space points requiring a CPU time between 5.0s < t VVAMP < 7.5 s, and about 0.1% exceeding 20 s (visible from the overflow bin). When using the interpolator, more than 99% of the evaluations just require some hundredths of a second. The small number of phase-space points with a CPU time t int > 0.5 s are associated to the values catched by the rescue-system. Those timings have been obtained on machines with Intel Haswell Xeon E5-2698 processors with 2.3 GHz per core.
Our implementation of a faster evaluation of the two-loop amplitudes through interpolation is tested by looking at its impact on physical predictions, especially on some relevant differential distributions in the inclusive phase space. All results have been obtained at the level of the Monte Carlo integration of the cross section (i.e Powheg stage 2), so that no parton shower radiation or hadronization effects, which would not be relevant for the validation, have been included.
First, it is worth mentioning that the code with the interpolation of the two-loop amplitudes reproduces accurately the exact inclusive cross section, with discrepancies of the order of about 0.4 permille, which are well within statistical uncertainties. Then, in figure 6 we show representative plots that compare the exact VVAMP predictions (blue, solid line) against the results with interpolation (black, dotted line) for the rapidity of the W + W − pair (y W W ), the rapidity difference between the two W -bosons (∆y W − , W + ), the rapidity of positively-charged W boson (y W + ), the experimental definition of the transverse mass of the W + W − pair the transverse momentum of the negatively-charged W boson (p T,W − ) and the missing transverse momentum (p T,miss ). We stress that many more distributions than those included in this manuscript have been carefully examined and verified to show a very good agreement between the analytic and interpolated results. Moreover, in order to highlight the phasespace regions where the two-loop contribution gives a large contribution to the cross section, a third curve (red, dotted line) has been included in all plots, obtained by setting the The lower panel of the plots displays the bin-by-bin ratio using the VVAMP curve as a reference.
From all plots, it is evident that the interpolator reproduces correctly the differential distributions in all kinematical regimes, with only small fluctuations at very high values of ∆y W − , W + (at most of the order of 2%) or high values of m T,W W or p T,W − (where differences are well below one percent). These are the regions where the two-loop contribution has the largest impact on the cross section. Indeed, from the y W W distribution it is evident that H (2) F has a 5-6% effect on the cross section, and contributes uniformly to this observable, while for instance for |∆y W − , W + | > 4.5 and |y W + | > 4 it induces a positive correction that reaches more than 30% and 10%, respectively. For transverse-momentum distributions, such as p T,W − or p T,miss , the two-loop contribution has a positive impact for relatively low transverse momenta (roughly for p T < 100 GeV) of at most 10%, while it yields an increasingly negative correction for very large transverse momenta.

Input parameters and settings
We consider √ s = 13 TeV proton-proton collisions at the LHC and present predictions for pp → + ν −ν + X production with = e and = µ. The EW parameters are determined in the G µ scheme, therefore computing the EW coupling as α Gµ = √ 2 G µ m 2 W (1 − m 2 W /m 2 Z )/π and the mixing angle as cos θ 2 W = m 2 W /m 2 Z . We use the following PDG [178] values as inputs: G F = 1.16639 × 10 −5 GeV −2 , m W = 80.385 GeV, Γ W = 2.0854 GeV, m Z = 91.1876 GeV, and Γ Z = 2.4952 GeV. We set the CKM matrix to unity, which, because of unitarity and the fact that we consider only massless external quarks is a very good approximation, as explained in ref. [89]. As described in section 2.1, the four-flavour scheme with N f = 4 massless quark flavours and massive bottom and top quarks is used to define a top-free W + W − cross section by removing all contributions with final-state bottom quarks. Accordingly, we use the N f = 4 NNLO set of the NNPDF3.0 [179] parton densities. More precisely, in case of MiNLO and MiNNLO PS the PDF grids are read from the lhapdf interface [180], copied into hoppet grids [167] and evaluated by hoppet for scales below the internal PDF infrared cutoff through DGLAP evolution with the number of active flavours kept fixed to the one at the internal PDF infrared cutoff, as described in ref. [76]. The central renormalization and factorization scales are set following the usual setting for MiNNLO PS and MiNLO discussed in section 2.3. Perturbative uncertainties are estimated from customary 7-point variations, i.e. by varying µ R and µ F around the central scale by a factor of two while respecting the constraint 0.5 ≤ µ R /µ F ≤ 2.
We compare our MiNNLO PS (and MiNLO ) predictions to the NNLOPS results presented in ref. [89]. Those are based on a MiNLO calculation with µ R = µ F = p T,W W , but use NNLO predictions for the reweighting with where m ν and p T, ν (m ν and p T, ν ) are the invariant masses and the transverse momenta of the reconstructed W bosons. The setting in eq. (3.1) is therefore the effective scale used in the NNLOPS calculation of ref. [89], where the perturbative uncertainties are obtained from 7-point scale variations that are assumed correlated in the reweighting. For the p T,W W spectrum and the jet-vetoed cross section we also compare against more accurate analytically resummed predictions obtained with Matrix+RadISH [68,131,132], where we have chosen  Q res by a factor of two around its central value, while keeping µ R and µ F fixed to µ 0 . For some observables it is instructive to also compare to fixed-order NNLO predictions with both the scale settings in eq. (3.1) and the ones in eq. (3.2), which we have obtained with Matrix [124,158]. In this case, perturbative uncertainties are again estimated from 7-point scale variations. As pointed out before, we do not include the loop-induced gluon-fusion contribution in all NNLO results throughout this paper and study the genuine NNLO corrections to the qq initiated process. The leading-order gluon-gluon initiated contribution enters at NNLO and NLO QCD corrections to it are known and can be incoherently added to the predictions presented here through a dedicated calculation, which is beyond the scope of this paper. Finally, for the matching to the parton shower we employ Pythia8 [181] with a A14 tune [182] (specifically py8tune 21). Since in this study our focus is on the comparison with other theory predictions, we do not include any effect from hadronization, underlying event models, or a QED shower. Such effects can, however, be directly included and studied by any user of our program through the Pythia8 interface of Powheg-Box-Res.
Since the W + W − cross section is finite at the LO without any cuts, we present results in the fully inclusive W + W − phase space, referred to as setup-inclusive. Additionally, we consider two sets of fiducial cuts, which are summarized in table 1. The first one corresponds to an earlier ATLAS 13 TeV analysis [13] and it is identical to that used in the NNLOPS calculation of ref. [89], which allows us to compare directly our MiNNLO PS predictions with the fiducial NNLOPS results of ref. [89]. We refer to it as fiducial-1-JV in the following. We note that fiducial-1-JV involves a two-fold jet-veto requirement, vetoing all jets in the rapidity region |η j | < 2.5 and separated from the leptons by ∆R ej > 0.3 with p T,j > 25 GeV and all jets in the rapidity region |η j | < 4.5 and separated from the leptons by ∆R ej > 0.3 with p T,j > 30 GeV. The second setup instead corresponds to the most recent ATLAS 13 TeV measurement of ref.
[12], and it was used to study high-accuracy resummed predictions for W + W − production in ref. [68]. This setup, referred to as fiducial-2-JV in the following, is useful for two reasons. First, at variance with fiducial-1-JV, it includes a single jet-veto cut for jets with p T,j > 35 GeV. This allows us to directly compare against the NNLO+NNLL resummed predictions for the p T,W W spectrum with a jet veto [68]. Note that to facilitate this comparison, we have removed the jet rapidity (η j ) requirement from fiducial-2-JV [12], which has a numerically tiny effect. Second, fiducial-2-JV is used to compare against data, since ref.
[12] measured the fiducial cross section as a function of the jet-veto cut to validate theory predictions for the jet-vetoed W + W − cross section. Let us recall that jet-veto requirements are crucial for experimental W + W − analyses in order to suppress the large top-quark backgrounds. In addition, we introduce fiducial-1-noJV and fiducial-2-noJV for the same fiducial setups as given in table 1, but without any restriction on the jet activity. Those are relevant to study the p T,W W distribution inclusive over jet radiation as well as the cross section as a function of the jet-veto cut. Besides jet-veto requirements, the two setups in table 1 involve standard cuts on the transverse momentum (p T, ) and pseudorapidity (η ) of the charged leptons as well as a lower cut on the invariantmass of the dilepton pair (m ) and on the missing transverse momentum (p T,miss ). Setup fiducial-2-JV includes also a lower cut on the transverse momentum of the dilepton pair (p T, ), while setup fiducial-1-JV cuts on the so-called relative missing transverse (p T,miss,rel ). The latter denotes the component of the missing transverse momentum vector perpendicular to the direction of the closest lepton in the azimuthal plane, and is defined as p T,miss,rel = p T,miss · sin |∆φ| for ∆φ < π/2 , p T,miss for ∆φ > π/2 , where ∆φ denotes the azimuthal angle between the missing transverse momentum vector vector and the nearest lepton.

Integrated cross sections
We start the presentation of phenomenological results by discussing integrated cross sections in table 2. In particular, we report predictions in the fully inclusive and the two fiducial phase spaces introduced in section 3.1 for MiNLO , MiNNLO PS , NNLOPS [89] as well as two fixed-order NNLO predictions obtained with Matrix [124,158] using the scale settings of eq. (3.1) and eq. (3.2). We summarize our main observations in the following: • It is clear that NNLO accuracy is crucial for an accurate prediction of the W + W − cross section, since the MiNLO result is about 12% lower than the MiNNLO PS one not only for setup-inclusive, but also after including the fiducial-1-JV and fiducial-2-JV cuts. In fact, in all cases the MiNNLO PS prediction is outside the upper uncertainty boundary of the MiNLO one. This is not surprising since for W + W − production also at fixed order the NLO uncertainty band does typically not include the central value of the NNLO prediction [60,89]. Additionally, the precision at NNLO is substantially improved, with scale uncertainties reduced by almost an order of magnitude.
• The NNLO-accurate predictions compare well against one another. They are compatible within their respective scale uncertainties and the central predictions are all setup-inclusive fiducial-1-JV fiducial-2-JV   [89] and to the NNLO cross section obtained with different settings of µ R and µ F . All NNLO corrections to qq-induced W + W − production are taken into account, while the loop-induced gg contribution is excluded. In the last rows, the comparison to CMS and ATLAS data is shown. For the measured inclusive cross sections we have assumed a branching fraction of BR(W ± → ± ν ) = 0.108987, consistently evaluated with our inputs, and applied one factor for each of the two W bosons. The measured fiducial cross sections have been divided by a factor two so that they correspond to pp → + ν − ν production with = e and = µ. In addition, we have subtracted the loop-induced gluonfusion contribution from the central value of the data. For the inclusive results and the fiducial-1-JV result we used the prediction for gg (non-resonant) cross section quoted in table 5 of the ATLAS analysis in ref. [13]. For the fiducial-2-JV result we used the ggLO result in table 2 of ref. [123]. The ATLAS measurement of ref.
[13] includes resonant Higgs contributions, which have been subtracted from that data as well, using the corresponding prediction quoted in table 5 of that paper.
within less than 2%. Indeed, small differences are expected from the fact that those predictions differ by terms beyond NNLO accuracy. Note that the NNLOPS and the NNLO calculations with µ 0 = 1 2 (m T,W + + m T,W − ) are very close, both in terms of central values and uncertainties, since the former is actually reweighted to the latter prediction in the inclusive phase space. The fact that the inclusive MiNNLO PS result is about 1.2% below the NNLOPS one is due to the different scale choice and treatment of terms beyond accuracy. Indeed, the second NNLO prediction with a scale choice of µ 0 = m T,W W is even slightly lower than the MiNNLO PS one.
• The agreement among predictions with NNLO accuracy gets even better as far as fiducial cross sections are concerned. This could be an indication that the jet vetos that are applied in the fiducial phase spaces reduce the impact of terms beyond NNLO accuracy.
• Some caution is advised regarding the quoted scale uncertainties. First of all, the quoted uncertainties generally appear to be quite small, and potentially at the edge of providing a realistic estimate of the true uncertainty. Clearly, the MiNLO uncertainty does not cover the inclusion of NNLO corrections through MiNNLO PS . As far as the NNLO-accurate results are concerned, the MiNNLO PS uncertainties in the inclusive phase space are even a factor of two smaller than the ones of the other predictions. Moreover, the fixed-order NNLO uncertainties further decrease when the jet-veto requirements are imposed. Such behaviour is not new [183], but at least the showered results show an increased uncertainty when imposing fiducial cuts. Still, especially for the jet-vetoed predictions one may consider more conservative approaches to estimate the perturbative uncertaintes, see for instance refs. [183,184].
• Finally, there is a good agreement of MiNNLO PS results with data from ATLAS and CMS in both inclusive and fiducial phase-space regions. The measured cross sections agree mostly within one and at most within two standard deviations.

Differential distributions
We now turn to discussing differential distributions. We start in section 3.3.1 and section 3.3.2 with comparing our MiNNLO PS to MiNLO and NNLOPS predictions in the inclusive and the fiducial phase space, respectively. This allows us, on the one hand, to study the effect of NNLO corrections through MiNNLO PS with respect to MiNLO and, on the other hand, to assess the compatibility of the MiNNLO PS predictions with the known NNLOPS results. Then in section 3.3.3 we move to distributions sensitive to soft-gluon radiation that require the inclusion of large logarithmic corrections to all orders in QCD perturbation theory either through a parton shower or through analytic resummation. Unless indicated otherwise, the plots are organized as follows: there is a main frame, which shows differential distributions for the MiNNLO PS (blue, solid), MiNLO (black, dotted), and NNLOPS (magenta, dash-dotted) predictions. In a lower frame we show binby-bin ratios of all curves to the central MiNNLO PS result. In some cases, where it is instructive to compare to the fixed-order results, we show fixed-order NNLO distributions for µ 0 = 1 2 (m T,W + + m T,W − ) (green, long-dashed) and/or for µ 0 = m T,W W (red, dashed). We note that we refrain from showing fixed-order NNLO predictions for most observables as the NNLOPS results correspond to a scale setting of µ 0 = 1 2 (m T,W + + m T,W − ) and are, in general, numerically very close to the respective fixed-order NNLO cross section. Additionally, for the p T,W W distribution we show NNLO+N 3 LL (green, double-dash dotted) and NNLO+NNLL (brown, dash-double dotted) predictions, and for the p T,W W spectrum with a jet veto as well as the jet-vetoed cross section we show NNLO+NNLL (green, doubledash dotted) and NLO+NLL (brown, dash-double dotted) results.

Inclusive phase space
We start by discussing distributions in the inclusive phase space. We have considered a large number of relevant distributions of both the leptonic final states and of the reconstructed W -bosons. A selection of those, which reflect some general features, is presented in figure 7. Since experimentally the W bosons can not be directly reconstructed and the fully inclusive phase space can not be covered by the detectors in any case, we follow here a more theoretical motivation and study observables related to the reconstructed W bosons rather than their decay products. In particular, figure 7 shows the transverse-momentum spectrum of the W + boson (p T,W + ), the rapidity distribution of the W -boson pair (y W W ), the rapidity difference between the two W bosons (∆y W − , W + ), and the invariant-mass distribution of the W -boson pair (m WW ).
For the p T,W + spectrum, the MiNNLO PS prediction is in full agreement with the NNLOPS result, which is particularly striking in the low-p T,W + region since scale uncertainties are only at the level of ±1%. At larger values of p T,W + , the uncertainty bands of the NNLO+PS accurate predictions widen and reach about ±5%. This indicates that this region is predominantly filled by higher-order (real) radiative corrections with at least one jet, and that the formal accuracy is somewhat decreased by one order. Indeed, in the region p T,W + 100 GeV the NNLO+PS predictions become fully compatible with the MiNLO result, also in terms of the size of the uncertainty bands. By contrast, for smaller p T,W + we observe large NNLO corrections with respect to MiNLO that reach almost 20% and substantially reduced scale uncertainties.
Also for the y W W and ∆y W − , W + distributions we find fully compatible results with overlapping uncertainty bands when comparing MiNNLO PS and NNLOPS predictions. While the NNLO corrections compared to MiNLO are relatively flat for y W W , we find that the corrections increase substantially at larger values of ∆y W − , W + , reaching ∼ +30% for y W W 3. This behaviour was observed already in ref. [89] and it is reassuring to see that this large effect is not a remnant of the scale setting in the NNLOPS calculation, but a genuine NNLO correction.
Similarly sizeable NNLO corrections are observed also at large values of m WW . This is also one of the few regions of phase space that we found where MiNNLO PS and NNLOPS predictions do not agree at the level of a few percent. While up to m WW 500 GeV the MiNNLO PS and NNLOPS results are fully compatible, they start deviating at larger invariant masses, reaching differences of about 20% at m WW = 1.8 TeV. Those differences can be traced back to the different scale settings in the MiNNLO PS and NNLOPS calculations. Indeed, comparing the additional NNLO results shown for the m WW distribution, we notice a relatively large spread between the green long-dashed curve with scale setting µ 0 = 1 2 (m T,W + + m T,W − ) and the red dashed curve with µ 0 = m T,W W , which is of the same size as (or even slightly larger than) the observed differences between MiNNLO PS and NNLOPS. As expected, the NNLOPS result is close to the NNLO one with µ 0 = 1 2 (m T,W + + m T,W − ), while the MiNNLO PS prediction is somewhat in-between the two NNLO results, but slightly closer to the one with µ 0 = m T,W W . Thus, the origin of the observed differences are terms beyond NNLO accuracy. Although the uncertainty bands increase to about 10% towards large m WW , the two NNLO+PS accurate predictions do not (or only barely) overlap for m WW 1 TeV, indicating that plain 7-point scale variations do not represent a realistic estimate of the actual size of uncertainties in that region of phase space. One may ask the question whether one of the two scale choices can be preferred. Although one might assume that m WW would be the natural scale of the m WW distribution, the situation is actually not that clear. This was discussed in some detail in ref. [89], and it boils down to the fact that for s-channel topologies m WW would be the more natural scale, while for t-channel topologies the transverse masses of the W bosons reflect better the natural scale of the process. Since both topologies appear in W + W − production already at the LO and they interfere, it is hard to argue in favour of any of the two scale choices. As a result, and since now there are two NNLO+PS accurate predictions available, the difference between the two should be regarded as an uncertainty induced by terms be-yond NNLO accuracy. Moreover, one could introduce a different setting of the hard scales at high transverse momenta in the MiNNLO PS calculation (i.e. in the W + W − +jet part) as a further probe of missing higher-order terms.
In summary, we find that MiNNLO PS and NNLOPS predictions are in excellent agreement for essentially all observables we considered in the inclusive phase space that are genuinely NNLO accurate. This indicates the robustness of NNLO+PS predictions for such observables. For the few exceptions, like large m WW , we could trace back the origin of the differences to terms beyond accuracy that are induced by the different scale settings. Moreover, in all cases the NNLO corrections substantially reduce the scale uncertainties with respect to the MiNLO prediction. Notice, however, that in the bulk region of the inclusive phase space the MiNNLO PS uncertainty bands are about a factor of two smaller than the NNLOPS ones, as already observed for the fully inclusive cross section.

Fiducial phase space
We continue our comparison by considering differential distributions in the fiducial-1-JV phase space. Here, we have selected a set of leptonic observables that are directly measured by the experimental analyses, cf. refs. [12,13], and which represent well the general features of all observables we considered. To this end, figure 8 shows the distributions in the transverse momentum of the leading lepton (p T, 1 ), in the invariant mass (m ), transverse momentum (p T, ), rapidity (y ) and azimuthal difference (∆φ ) of the dilepton pair, and in an observable particularly sensitive to new physics effects defined through the separation in η of the two leptons: |cos(θ )| = |tanh (∆η /2)| . (3.4) As for the setup-inclusive in the previous section, we find full compatibility between MiNNLO PS and NNLOPS predictions. With fiducial cuts, even the differences induced by terms beyond accuracy are reduced and the scale-uncertainty bands of the two calculations are of similar size. Also in this case, an important observation is that the inclusion of NNLO corrections on top of the MiNLO is crucial not only for the correct normalization, but for many observables also to capture relevant shape effects. Moreover, the NNLO-accurate predictions are substantially more precise due to their strongly reduced uncertainty bands with respect to MiNLO . We further notice that the impact of parton-shower emissions on observables with NNLO accuracy is quite moderate. Nevertheless, at phase-space boundaries where the fixed-order accuracy is reduced and the perturbative expansion breaks down due to effects from soft QCD radiation, the parton shower is absolutely crucial for a physical description. For instance, this can observed in the p T, distribution, where we have added the fixed-order NNLO prediction for comparison. Since the p T,miss > 20 GeV requirement in fiducial-1-JV setup translates directly into a p T, > 20 GeV cut at LO, where the two leptons are back-to-back with the two neutrinos, the region p T, ≤ 20 GeV is filled only upon inclusion of higher-order corrections and is effectively only NLO accurate. As a result, the boundary region p T, ∼ 20 GeV becomes sensitive to soft-gluon effects that induce large logarithmic corrections and a perturbative instability [185] at p T, = 20 GeV in the dσ/dp T,ℓ1 [fb/GeV] fixed-order NNLO prediction. This unphysical behaviour is cured through the matching to the parton shower in the MiNNLO PS and NNLOPS calculations.
It is clear that our new MiNNLO PS predictions compare very well with the previous NNLOPS results, and that the two tools can be used equivalently to produce W + W − cross sections and distributions at NNLO accuracy matched to parton showers. This is also an indication of the robustness of NNLO+PS predictions for observables that are genuinely NNLO accurate. Given the limitation of the NNLOPS calculation regarding the necessity of multi-dimensional reweighting, the advantage of the MiNNLO PS generator is that those results can be obtained directly at the level of the event generation. However, in the few phase-space regions where differences between the two calculations can be observed, those differences indicate relevant corrections beyond NNLO accuracy. Since plain 7-point scale variations do not always cover those discrepancies, they should be regarded as a perturbative uncertainty.
In the next section we will move to observables that are subject to large logarithmic corrections and where differences between the MiNNLO PS and NNLOPS generator are larger. Thus, their assessment as an uncertainty becomes particularly important.

Observables sensitive to soft-gluon effects
In figure 9, we study the transverse-momentum spectrum of the W + W − pair (p T,W W ) in the fiducial-1-noJV phase space. We refrain from showing the corresponding distribution in setup-inclusive and within the fiducial-2-noJV phase space, as we found them to be almost identical concerning the relative behaviour of the various predictions. At small values of p T,W W , large logarithmic contributions break the validity of the expansion in the strong coupling constant at a given fixed order, which requires their inclusion all orders in perturbation theory either through a parton shower or through an analytic resummation. The left figure displays the region 0 ≤ p T,W W ≤ 50 GeV and, indeed, the NNLO prediction, which is shown in the main frame only, becomes unphysical for small values of p T,W W . If we compare MiNNLO PS and NNLOPS results in that region, we observe differences of about −10% to +5%. By and large, those are covered by the respective uncertainty bands. However, it is clear (and expected) that for such an observable, which is sensitive to infrared physics, the differences between the two calculations become larger. In particular, both predictions are only NLO accurate in the tail of the p T,W W distribution and at small transverse momenta the parton shower limits the accuracy of the calculation effectively to leading-logarithmic (LL) or partial (i.e. at leading colour) next-to-LL (NLL) accuracy. Therefore, differences of the order of those that we observe between MiNNLO PS and NNLOPS are understood. Also the comparison against the high-accuracy analytic resummation results at NNLO+N 3 LL and NNLO+NNLL is quite good, which also agree within −10% to +5% with the MiNNLO PS prediction for p T,W W < 20 GeV and are even fully compatible in the intermediate region up to 50 GeV. The resummed predictions do not favour either MiNNLO PS or NNLOPS results, but rather show similar differences to the two. On the other hand, the agreement is actually quite remarkable considering the fact that the region p T,W W < 20 GeV is entirely described by the substantially less accurate parton shower. Given the fact that for some bins the NNLO+N 3 LL and NNLO+NNLL predictions are outside the uncertainty bands dσ/dp T,WW [fb/GeV] pp→ℓ + ν ℓ ℓ' − ν ℓ' @LHC 13 TeV of the NNLO+PS accurate predictions though, the estimated uncertainties from µ R and µ F variations appear insufficient to reflect the actual size of uncertainties and one should consider additional handles to better assess the uncertainties of the parton shower at small p T,W W . Indeed, the NNLL prediction has a much larger uncertainty band in this region (induced by the variation of Q res ) even though it is more accurate. In the right plot of figure 9, we show the range 0 ≤ p T,W W ≤ 250 GeV. In the tail of the distribution, MiNNLO PS and NNLOPS (as well as MiNLO ) predictions are in perfect agreement with fully overlapping uncertainty bands. In the lower frame we show an additional curve that is ratio of the central fixed-order NNLO prediction with µ 0 = 1 2 (m T,W + + m T,W − ) to the one with µ 0 = m T,W W . It is very interesting to observe that the ratio corresponds almost exactly to the NNLOPS/MiNNLO PS ratio at smaller p T,W W . We recall that µ 0 = 1 2 (m T,W + + m T,W − ) is the scale used in the reweighting of the NNLOPS prediction, while µ 0 = m T,W W is somewhat more similar to the one within the MiNNLO PS approach. This suggests that the differences originating from terms beyond accuracy at small p T,W W between the MiNNLO PS and NNLOPS are predominantly induced by the different scale settings in the two calculations. In fact, for any distribution (of the various ones we considered) where the NNLOPS/MiNNLO PS ratio becomes larger than a couple of percent, we observe that the corresponding ratio of fixed-order NNLO predictions is either very similar or even larger.
In figure 10 we consider the W + W − transverse momentum spectrum in the presence of a jet veto of p veto T,j 1 = 35 GeV using the fiducial-2-JV setup. The relative behaviour between the MiNNLO PS , NNLO+PS, NNLO+NNLL and NLO+NLL results at small transverse momenta is relatively similar to the one observed for the p T,W W distribution without jet dσ/dp T,WW [fb/GeV] pp→ℓ + ν ℓ ℓ' − ν ℓ' @LHC 13 TeV veto in setup fiducial-1-noJV. One main difference is that for this observable, which requires double differential resummation in p T,W W and p T,j 1 , the analytically resummed results are less accurate and therefore feature larger uncertainty bands, rendering them more compatible with the showered results. Indeed, the NLL uncertainty band is strongly increased at small p T,W W and much larger than the NNLO+PS one, which, as argued before, also points to the fact that the scale uncertainties of the latter are somewhat underestimated, given that the parton shower is less accurate than the NLL calculation in that region. Another interesting region for this observable is around p T,W W values of 35 GeV, i.e. of the order of the jet-veto cut. The region p T,W W ≥ p veto T,j 1 is filled for the first time at NNLO, which is effectively only LO accurate, since at LO it is p T,W W = 0 and at NLO p T,W W = p T,j 1 . Therefore, large logarithmic contributions challenge the perturbative expansion around p T,W W = p veto T,j 1 and the fixed-order NNLO prediction develops a perturbative instability, as visible in the main frame of the left plot in figure 10. This instability is partially cured by the analytic resummation approach, which resums Sudakov logarithms in the limit where p T,W W and p T,j 1 are much smaller than the hard scale, but not all logarithmic contributions of the form log(p T,W W −p veto T,j 1 ), which would require additional resummation when one or more hard jets are present. By contrast, the NNLO+PS calculations cure this instability entirely as they resum all relevant classes of logarithms (although with limited accuracy). Therefore, the MiNNLO PS and NNLOPS calculations provide a more physical prediction at and above threshold, while below the threshold they are in good agreement with the analytically resummed predictions. If we look at region above threshold in the right plot of figure 10, we σ(p T,j1 < p T,j1 veto ) [fb] pp→ℓ + ν ℓ ℓ' − ν ℓ' @LHC 13 TeV notice that the NNLO result drops substantially for p T,W W values above p veto T,j 1 , and also the NNLO+NNLL prediction is only slightly larger. Hence, this region of phase space is almost exclusively filled by the parton shower. Consequently, the transverse-momentum spectrum of a colour singlet in presence of a jet veto could be a good observable to tune the parton shower in experimental analyses.
In figure 11 we study the jet-vetoed cross section as a function of the jet-veto cut p veto T,j 1 , which is defined as σ(p T,j 1 < p veto T,j 1 ) = p veto T,j 1 0 dp T,j 1 dσ dp T,j 1 , (3.5) and the jet-veto efficiency given by where σ int is the integrated cross section in the fiducial-1-noJV phase space. Again the results for setup-inclusive and fiducial-2-noJV are very similar and are not shown. The interesting region is at small jet-veto cuts, where the validity of the perturbative expansion is broken by large logarithmic contributions in p veto T,j 1 , while for larger values the results tend towards their respective integrated cross sections. As it can be seen from the main frame, in the low p veto T,j 1 region the pure fixed-order result at NNLO becomes indeed unphysical and turns actually negative. When comparing MiNNLO PS and NNLOPS predictions, we find them to be in reasonable agreement within their respective uncertainties, with the NNLOPS one tending a bit faster towards zero for p veto T,j 1 20 GeV. In that region, the resummed NNLO+NNLL and NLO+NLL results tend even faster towards zero, with the NNLO+NNLL curve being about 20% below the MiNNLO PS one at p veto T,j 1 = 5 GeV.  Figure 12: Jet-vetoed cross section in the fiducial-2-noJV phase space compared to data. As described in the caption of table 2 the data has been adjusted by subtracting the ggLO contribution quoted in table 2 of ref. [123] and by dividing out a factor of two.
This region is dominated by the parton shower, which resums only the LL (partial NLL) contributions. Clearly, the actual uncertainties in the NNLO+PS calculations are not covered by plain µ R and µ F variations. As argued for the p T,W W spectrum, additional handles would need to be considered to better assess the parton-shower uncertainties for very small p veto T,j 1 cuts. Indeed, the NLL result features much wider uncertainties, despite being more (similarly) accurate in that region of phase space. However, we stress that such small p veto T,j 1 cuts are usually not relevant for experimental W + W − analyses. Moreover, as pointed out before, there have been suggestions to include more conservative uncertainty estimates for jet-vetoed predictions [183,184]. We leave their proper assessment to future work, as those effects are currently not accessed by any W + W − measurement. For instance, looking at the fiducial phase-space definitions of refs. [12,13] that are considered in this paper, jetveto cuts of p veto T,j 1 = 25 GeV, 30 GeV and 35 GeV are used. For those values, MiNNLO PS predictions are in perfect agreement with the NNLO+NNLL resummation, and even down to p veto T,j 1 ∼ 15 GeV they differ by less than 5% with overlapping uncertainties. When comparing the predicted jet-vetoed cross section as a function of p veto T,j 1 in the fiducial-2-noJV setup against data in figure 12, it is clear that the MiNNLO PS and the NNLO+NNLL prediction are fully compatible in the relevant region. The agreement with data is good in either case, with the data points either marginally overlapping within one standard deviation or being just outside this range. One should bear in mind however that the normalization of the theory prediction can be increased by ∼ 5% just by using a different PDF set, which yields even better agreement with data, as shown in ref. [68]. Apart from that, it is clear that the inclusion of NNLO corrections brings the theory predictions closer to data.

Conclusions
In this paper we have presented the matching of NNLO-accurate predictions with parton showers for W + W − production at the LHC using the MiNNLO PS method. We have performed the calculation consistently in the four-flavour scheme with massive bottom quarks. By dropping contributions with final-state bottom quarks, which are regulated by the finite bottom mass, we generate top-free W + W − events. We have presented an extensive comparison of our MiNNLO PS predictions against the NNLOPS results of ref. [89]. We find excellent agreement with the latter results, with only minor differences in phase-space regions where they are expected from the different treatment of terms beyond accuracy, most importantly the ones related to different scale settings. Especially for genuine NNLO observables MiNNLO PS and NNLOPS are compatible within less than a few percent. Larger differences (but still within uncertainties) can be observed for observables that are sensitive to the limited accuracy of the parton shower, in particular for small values of the W + W − transverse momentum or for very small jet-veto cuts. For those observables we also compared to high-accurate analytically resummed predictions. The agreement is very reasonable considering the limited (leading logarithmic) accuracy of the parton shower. In particular for phase-space regions relevant for experimental W + W − analyses, i.e for jet-veto cuts of 25 GeV to 35 GeV (and higher), the MiNNLO PS prediction essentially coincides with the NNLO+NNLL result. Indeed, we find good agreement both for inclusive and fiducial cross sections with experimental data, as well as for the cross section as a function of the jet-veto cut.
We found that the major bottleneck in the computation is the evaluation of the twoloop amplitude. To improve the speed, we have constructed four-dimensional grids of the coefficients that encode all the information required to reproduce the full two-loop contribution. We then use those grids to obtain the coefficients at any given phase-space configuration through four-dimensional cubic spline interpolation and reconstruct the twoloop amplitude on-the-fly. As a result, the evaluation time of the two-loop contribution has been reduced by a factor of forty, becoming subleading with respect to the other parts of the calculation. We have performed a thorough and extensive validation of the results with interpolation against the ones without. The MiNNLO PS code can be used either with or without interpolator, the latter option being about five times slower.
The MiNNLO PS approach has various positive features. First, NNLO corrections are calculated on-the-fly during the generation of the events, with no need for further reprocessing or reweighting of the events. Second, no unphysical merging scale needs to be introduced to separate event samples of different multiplicities, a concept already introduced and discussed in detail in ref. [72]. Third, when combined with transverse-momentum ordered parton showers, the matching guarantees that the logarithmic accuracy of the parton shower simulation is preserved. We note that, because the logarithmic accuracy of parton showers is only leading logarithmic, it is often taken for granted that this accuracy is preserved. On the contrary, this is a subtle, but crucial point in any approach that combines NNLO and parton showers. In the MiNNLO PS case this requirement is immediately fulfilled for transverse-momentum ordered showers, since zero, one, and two emissions are included in the NNLO weight, where the second-hardest is generated following the Powheg matching procedure, while additional emissions are generated by the parton shower.
We expect that the MiNNLO PS code associated to the work presented in this paper, which enables an accurate fully-exclusive hadron-level generation of W + W − events, will be highly valuable for experimental measurements, which require an accurate simulation of W + W − production either as signal or as background to other processes. The code will be publicly released within Powheg-Box-Res and will supersede the previous NNLOPS approach based on multi-differential reweighting in Powheg-Box-V2. Nevertheless, we stress that it may be useful to take differences between MiNNLO PS and NNLOPS predictions (or fixed-order NNLO with different scale settings) to assess residual uncertainties in certain phase-space regimes. Alternatively, different scale settings in the MiNNLO PS calculation and handles in the parton shower could be used to probe the size of terms beyond accuracy.
Finally, any meaningful comparison to data of differential distributions in chargeneutral vector-boson pair production processes should include not only NNLO QCD accuracy for the qq channel, but also NLO QCD corrections to the loop-induced gg process, and possibly NLO EW corrections in the high-energy tails. We leave such studies to future work. [17] CMS collaboration, Measurement of the W + W − cross section in pp collisions at sqrt(s) = 8 TeV and limits on anomalous gauge couplings, CMS-PAS-SMP-14-016.