Accurate single- and double-differential resummation of colour-singlet processes with MATRIX+RadISH: $W^+W^-$ production at the LHC

We present the combination of fully differential cross sections for colour-singlet production processes at next-to-next-to-leading order (NNLO) QCD obtained with MATRIX and all-order resummation through RadISH. This interface allows us to achieve unprecedented accuracy for various transverse observables in $2 \rightarrow 2$ production processes. As an important application we consider $W^+ W^-$ production at the LHC, more precisely the full leptonic process $pp\to \ell^+\ell^{\prime\, -}\nu_{\ell}{\bar\nu}_{\ell^\prime}+X$ with $\ell'\neq\ell$, and we present resummed predictions for differential distributions in presence of fiducial selection cuts. In particular, we resum the transverse-momentum spectrum of the $W^+ W^-$ pair at next-to-next-to-next-to-leading logarithmic (N$^3$LL) accuracy and match it to the integrated NNLO cross section. The transverse-momentum spectrum of the leading jet in $W^+ W^-$ production is calculated at NNLO+NNLL accuracy. Finally, the joint resummation for the transverse-momentum spectrum of the $W^+ W^-$ pair in the presence of a jet veto is performed at NNLO+NNLL. Our phenomenological study highlights the importance of higher-order perturbative and logarithmic corrections for precision phenomenology at the LHC.


Introduction
Precision phenomenology has become of major importance in the rich physics programme at the Large Hadron Collider (LHC). Lacking any direct observations of physics beyond the SM (BSM), precision measurements provide a valuable alternative in the discovery of newphysics phenomena through small deviations from the Standard Model (SM) picture. The vast amount of data collected at the LHC continuously decreases experimental uncertainties, thereby demanding accurate predictions for various observables in many relevant physics processes.
Differential distributions in colour-singlet production processes and associated QCD radiation play a special role in this context. Being measured by reconstructing the colour singlet from its leptonic decay products, these processes usually provide very clean experimental signatures. Therefore, such processes are often characterised by particularly small experimental uncertainties, which in the case of Drell-Yan production can reach subpercent precision and are at the few-percent level even for vector-boson pair production processes. Due to their precise measurement, leptonic processes induced by vector-boson decays provide prime signatures to extract SM parameters, to constrain parton densities, or to calibrate event-generation tools used in experimental analyses. Not least, vector-boson processes have a high sensitivity to new-physics phenomena, and they can be exploited to set stringent limits on BSM effects, which induce shape distortions in kinematic distributions.
It is well known that in kinematical regimes dominated by soft and collinear QCD radiation the perturbative expansion in the strong coupling constant α s does not provide a physical description of differential observables. An important class of distributions is that of transverse observables, which do not depend on the rapidity of the radiation. For colour-singlet processes prime examples of such observables are the transverse momentum of the colourless final state p T or of the leading accompanying jet p J T . In phase space regions dominated by soft and collinear radiation, the perturbative expansion of the cross section is marred by large logarithms L = ln(v), where the symbol v denotes a general transverse observable, such as p T or p J T . A resummation of these logarithmically-enhanced terms to all orders in α s is required to obtain physical results when v becomes small. The logarithmic accuracy is customarily defined in terms of the logarithm of the cumulative cross section ln σ(v). One refers to the dominant terms α n s L n+1 as leading logarithmic (LL), to terms α n s L n as next-to-leading logarithmic (NLL), to α n s L n−1 as next-to-next-to-leading logarithmic (NNLL), and so on.
A variety of formalisms to perform the resummation of large logarithmic contributions for transverse observables in colour-singlet processes has been developed over the last four decades [24][25][26][27][28][29][30][31][32][33][34][35][36][37]. The most accurate description of the p T spectrum has been achieved so far for Higgs boson production and for the neutral and charged Drell-Yan production in refs. [37][38][39][40], with N 3 LL resummation matched to NNLO QCD predictions for F +jet production of refs. [41][42][43][44], where F denotes the respective colour singlet 1 . Similarly, the resummation for p J T has been calculated at NNLL accuracy [46][47][48][49], and it has been used to produce predictions for jet-vetoed cross sections up to N 3 LO+NNLL accuracy in Higgs boson production [50]. Very recently, a formalism to simultaneously resum both classes of logarithms at NNLL accuracy has been presented in ref. [51], where it was applied to obtain the Higgs p T spectrum with a jet veto at NNLL matched to the NLO prediction for H+jet production.
Although these resummation formalisms can be generalised to the production of an arbitrary colourless final state, there is only a limited number of calculations that have included high-accuracy resummation in 2 → 2 processes due to their higher complexity.
In this paper, we present the Matrix+RadISH framework for high-accuracy resummed calculations at the multi-differential level. By developing a general interface between the Matrix and the RadISH codes, a substantial advancement over previous resummed predictions is achieved for a large number of non-trivial colour-singlet processes at the LHC. All the 2 → 1 and 2 → 2 processes available in the public release of Matrix are included, and essentially any colour-singlet process for which two-loop amplitudes become available can be added. Through its powerful and versatile parton-level Monte Carlo generator, Matrix+RadISH provides an accurate description of several transverse observables in colour-singlet processes. In particular, we perform transverse-momentum resummation of the colour-singlet final state at N 3 LL accuracy, φ * η [70] resummation for the Drell-Yan process at N 3 LL accuracy, transverse-momentum resummation of the leading jet (and equivalently jet-veto resummation) at NNLL accuracy, and double-differential resummation in the transverse momentum of the colour singlet and of the leading jet at NNLL accuracy. The latter allows for the consistent calculation of the p T (p J T ) spectrum with a veto cut on p J T (p T ). The matching is performed to the integrated cross section at NNLO QCD accuracy.
The resummation is formulated in the RadISH formalism of refs. [35,37], which allows us to resum transverse observables at high accuracy. The RadISH code also contains the implementation of the double-differential resummation on the basis of ref. [51]. The fixedorder component, the phase space and the relevant perturbative ingredients are evaluated through the computational framework Matrix [55,69].
As a first phenomenological application we study the production of W + W − pairs in 2 Here and in what follows the fixed-order accuracy refers to the one of the integrated cross section, not of the spectrum. Thus, matching to the (N)NLO cross section implies that accuracy only after integration over the resummed observable, and (N)LO accuracy in the spectrum (or accordingly for the F +jet process).
hadronic collisions. More precisely, we consider the full leptonic process with two charged leptons of different flavour and the two corresponding neutrinos in the final state. By evaluating all resonant and non-resonant contributions we include off-shell effects and spin correlations. In this paper, we advance the current state of the art of predictions in terms of accuracy for various observables at the multi-differential level in W + W − production.
In particular, we compute the fiducial cross section as a function of the jet-veto cut at NNLO+NNLL accuracy, and compare our predictions with recent ATLAS data [71]. Furthermore, we calculate fiducial predictions for the p J T spectrum at NNLO+NNLL accuracy and for the transverse momentum of the W + W − pair at NNLO+N 3 LL accuracy. Finally, we present double-differentially resummed predictions for the transverse-momentum spectrum of the W + W − pair at NNLO+NNLL accuracy in presence of a veto on p J T . The manuscript is organised as follows: In section 2 we give a general introduction to Matrix+RadISH. In particular, we review the computation of NNLO corrections for colour-singlet production within Matrix (section 2.1), provide the relevant formulae for the resummation of transverse observables in the RadISH formalism (section 2.2), discuss their matching (section 2.3), and give details on the practical implementation (section 2.4). In section 3 we discuss the case of W + W − production at the LHC, and we present results for fiducial predictions of transverse observables both at the single-differential and at the double-differential level. The main results are summarised in section 4, and we provide a practical description of how to use the Matrix+RadISH interface in appendix A.

Description of Matrix+RadISH
In this section, we discuss the calculation of all-order resummation matched to fixed-order predictions with Matrix+RadISH. Our implementation is completely general, and it can be applied to essentially any colour-singlet process. The combination of Matrix and RadISH facilitates the calculation of consistently resummed and matched predictions for several observables. The results are accurate up to NNLO in QCD perturbation theory for the fully differential cross section of the produced colour-singlet final state. In particular, the framework allows us to evaluate the following resummed predictions at unprecedented precision for 2 → 2 colour-singlet production processes: • single-differential resummation for the transverse-momentum spectrum of a coloursinglet (p T ) up to NNLO+N 3 LL accuracy, • single-differential resummation for the φ * η distribution in the Drell-Yan process up to NNLO+N 3 LL accuracy, • single-differential resummation for the transverse-momentum distribution of the leading jet (p J T ) (and equivalently for the jet-vetoed cross section) up to NNLO+NNLL accuracy, • double-differential resummation of p T and p J T logarithms up to NNLO+NNLL accuracy, allowing us to evaluate the p T (p J T ) spectrum with a veto on p J T (p T ) at NNLO+NNLL accuracy.
The calculations are fully differential in all Born level observables, and arbitrary fiducial cuts can be applied to the final state particles. 3

Higher-order corrections with Matrix
Matrix is a general framework for fixed-order calculations in QCD and EW perturbation theory, covering a large number of primary LHC scattering processes. The public release of Matrix [69,72] evaluates NNLO QCD predictions for 2 → 1 and 2 → 2 coloursinglet processes [13, 14, 16-19, 22, 23] 4 , including all possible leptonic decay channels of the massive vector bosons, while consistently accounting for resonant and non-resonant diagrams, off-shell effects and spin correlations. More recently, Matrix predictions have been further advanced by including important effects beyond NNLO QCD: The dominant next-to-NNLO (N 3 LO) QCD corrections have been implemented by calculating the loopinduced gluon fusion contribution at NLO QCD for ZZ [76] and W + W − [77] production, and the combination of NNLO QCD and NLO EW corrections has been achieved for all the leptonic final states of massive diboson processes [78]. 5 A new release of Matrix with these corrections is currently in preparation.
In this paper, we make use of the general implementation of fully differential NNLO cross sections in QCD perturbation theory for colour-singlet processes within Matrix. The computation of NNLO QCD corrections requires the evaluation of tree-level contributions with zero, one and two additional partons, of one-loop contributions with zero and one parton and of purely virtual two-loop contributions. Their combination in a fully differential (exclusive) calculation at NNLO QCD is highly non-trivial since infrared (IR) divergences affect real and virtual contributions in different ways, preventing a straightforward combination of these components. To overcome these issues, Matrix features a fully general implementation of the q T -subtraction method [80] at NNLO QCD, which is briefly described below. In this context, an automated extrapolation procedure has been implemented to calculate integrated cross sections in the limit in which the q T -subtraction cutoff parameter goes to zero [69]. The core of the Matrix framework [69] is the Monte Carlo program Munich 6 , which is capable of computing both NLO QCD and NLO EW [81,82] corrections to arbitrary SM processes. All tree-level and one-loop amplitudes are supplied by Open-Loops [83][84][85] through an automated interface. For validation and stability tests of the employed amplitudes a similar interface to the Recola amplitude generator [86,87] has been implemented. At two-loop level, for massive diboson production the public C++ library VVamp [88] is used that implements the qq → V V and gg → V V helicity amplitudes of refs. [89,90]. 7 For V γ [92] and γγ [93] production we rely on private implementations of the respective amplitudes.
The q T -subtraction formalism [80] exploits the fact that the behaviour of the cross section at small transverse momentum of a colour-singlet final-state system has a universal (process-independent) structure that is explicitly known up to NNLO QCD through the formalism of transverse-momentum resummation [25,27]. This knowledge is sufficient to construct a non-local, but process-independent IR subtraction counterterm for this entire class of processes. In the q T -subtraction method, the NNLO cross section for a general process pp → F + X, where F is a colourless system, is written as where dσ F +jet NLO is the cross section for the production of the system F and a jet at NLO accuracy, which can be evaluated by using one of the available NLO subtraction methods [94][95][96][97] to cancel the corresponding IR divergencies. In fact, unless the transverse-momentum of the colour singlet approaches zero, dσ F +jet NLO is finite. The process-independent counterterm dσ CT NNLO cancels the remaining divergence in the limit of vanishing transverse momentum, and it is constructed by expanding the transverse-momentum resummation formula [25,27] up to NNLO. The computation is completed by the last term on the right-hand side of eq. (2.1) that depends on the hard-collinear function H F NNLO up to NNLO [98][99][100][101]. The practical implementation of the q T -subtraction formalism in Matrix deserves some additional discussion. The contribution in the square bracket in eq. (2.1) is formally finite, but each individual term dσ F +jet NLO and dσ CT NNLO is separately divergent. Since the subtraction is not local, a technical cut-off r cut on the dimensionless quantity r = p T /M , where p T is the transverse momentum and M is the invariant mass of the colourless system, is introduced, rendering both terms separately finite. Below this cut-off dσ F +jet NLO and dσ CT NNLO are assumed to be identical, which is correct up to power-suppressed terms. These power-suppressed terms vanish only in the limit r cut → 0, but their impact is controlled by monitoring the dependence of the cross section on r cut . To this end, Matrix simultaneously computes the cross section at several r cut values without the need of repeated CPU-intensive runs, which is used to extrapolate the cross section in the r cut → 0 limit by fitting the results at finite r cut values. The extrapolated result and an estimate of the respective uncertainty are provided at the end of every Matrix run. We note that the q T -subtraction method works very similar to a phase space slicing method in this way, with r cut acting as a slicing parameter.
To perform the resummation of large logarithmic contributions, we have implemented a general interface to combine the Matrix framework with the RadISH code, which is introduced in the next section. In this context, Matrix provides all the fixed-order parts of the calculation as well as the Born level phase space points and the hard coefficients needed for the calculation of the resummed component. The latter are passed to RadISH which evaluates the resummation for the observable under consideration. In this paper we present a detailed study of the phenomenological implications for W + W − production only, but our implementation is completely general, and can directly be used for any of the other colour-singlet processes available in Matrix.

Resummation of large logarithmic contributions with RadISH
The RadISH approach for the resummation of transverse observables in colour-singlet processes has been presented in refs. [35,37] and will be summarised in the following. By exploiting the factorisation properties of squared QCD amplitudes and the recursive infrared collinear safety (rIRC) of the considered observables the resummation is formulated directly in momentum space, thereby obtaining a more differential description of the QCD radiation than that in customary conjugate-space formulations. The resummation is numerically evaluated via efficient Monte Carlo methods, yielding a powerful formalism similar in spirit to a semi-inclusive parton shower, but with the consistent inclusion of higher-order logarithmic contributions and full control over the formal accuracy. Thanks to its versatility, the approach can be exploited to formulate the resummation for the entire class of transverse observables, i.e. those which do not depend on the rapidity of the radiation, in a unique framework. This enabled the recent extension to double-differential resummation of the transverse-momentum spectrum of the colour singlet with a jet veto in ref. [51], as we will discuss below.
The RadISH formulae for the resummation of transverse observables are conveniently expressed in terms of the cumulative cross section where the transverse observable v = V (Φ B , k 1 , . . . , k n ) is a function of the Born phase space Φ B of the produced colour singlet and of the momenta k 1 , ..., k n of n real emissions. The all-order structure of the cumulative distribution σ(v), differential in the Born phase space, can be expressed as where M denotes the matrix element for n real emissions and V(Φ B ) is the resummed form factor that encodes the purely virtual corrections [102]. The phase spaces of the i-th emission k i and of the Born configuration are denoted by [dk i ] and dΦ B , respectively. For observables which fulfil rIRC safety [103] it is possible to establish a well-defined logarithmic counting in the squared amplitude, thereby providing a systematic way to identify the contributions that enter at a given logarithmic order [103,104]. This is achieved by decomposing the squared amplitude defined in eq. (2.3) in n-particle-correlated blocks containing the correlated portion of the squared n-emission soft amplitude and its virtual corrections [37] such that blocks with n particles start contributing one logarithmic order higher than blocks with (n − 1) particles.
The rIRC safety of the observables is further exploited to ensure that the divergences of virtual origin, contained in the V(Φ B ) factor of eq. (2.3), cancel those appearing at all perturbative orders in the real matrix elements. Indeed, the rIRC safety of the observable allows us to introduce a resolution scale q 0 on the transverse momentum of the radiation such that neglecting radiation softer than q 0 in the computation of V (Φ B , k 1 , . . . , k n ) only introduces terms suppressed by powers of q 0 . This unresolved radiation can thus be neglected when computing V (Φ B , k 1 , . . . , k n ), and it can be exponentiated to cancel the divergences of virtual origin at all orders. Resolved radiation, i.e. radiation harder than q 0 , must instead be generated exclusively since it is constrained by the measurement function Θ (v − V (Φ B , k 1 , . . . , k n )) in eq. (2.3). The rIRC safety of the observable also ensures that the limit q 0 → 0 can be taken safely since the dependence of the results upon the resolution scale is power-like.
The discussion above is completely general, and it can be applied to any transverse observable. We start by discussing a particular class of observables, that of inclusive observables, which depend solely on the total transverse momentum of QCD radiation. For clarity, we will consider the case of the transverse momentum of a general colour singlet. Nevertheless, the same formulae can be applied to any inclusive transverse observable such as the φ * η angle in Drell-Yan production [70], which was resummed at N 3 LL accuracy in ref. [39]. All the ingredients for the N 3 LL p T resummation have been computed in refs. [100,101,[105][106][107][108][109][110]. In the RadISH formalism, the resummation is numerically evaluated by setting the resolution scale q 0 to a small fraction δ > 0 of the transverse momentum of the block with largest k t , henceforth denoted by k t1 . As a result, the cumulative cross section in momentum space at N 3 LL accuracy for the production of a colour singlet of mass M , fully differential in the Born variables, reads [37] where the first line contains the full NLL correction, the first set of curly brackets (second to fifth line) starts contributing at NNLL accuracy, and the second set of curly brackets (from line six) is a pure N 3 LL correction. The luminosity factors L in eq. (2.4) are evaluated at different orders, and they involve the parton luminosities, the process-dependent squared Born amplitude and hard-virtual corrections H (n) , and the coefficient functions C (n) ci , which have been evaluated to second order for gg-and qq-initiated processes in refs. [100,101,105,106]. The factorsP (0) are the regularised splitting functions. We refer the reader to section 4 of ref. [37] for the definition of the luminosity factors and their ingredients. We further defined ζ si ≡ k tsi /k t1 and introduced the notation dZ[{R , k i }] to denote an ensemble that describes the emission of n identical independent blocks. In this notation, we define the average of a function (2.5) We stress that the ln 1/δ divergence that appears in the exponential prefactor of eq. (2.5) cancels exactly against that contained in the resolved real radiation, which is instead encoded in the nested sums of products on the right-hand side of the same equation.
To obtain eq. (2.4), we exploited the rIRC safety of the observable to expand all the ingredients in eq. (2.4) about k t1 since ζ i = k ti /k t1 ∼ O(1). Indeed, rIRC safety guarantees that blocks with k ti k t1 are fully cancelled by the term exp{−R (k t1 ) ln(1/δ)} of eq. (2.5). Such an expansion is not strictly necessary, but makes a numerical implementation much more efficient. Because eq. (2.4) was expanded about k t1 , it contains explicitly the derivatives of the radiator R, which is given by with α s = α s (µ R ) and µ R ∼ M being the renormalisation scale. The functions g i are reported in eqs. (B.8-B.11) of ref. [37]. With this choice, the logarithmic accuracy is effectively defined in terms of L = ln(Q/k t1 ), where Q ∼ M is the resummation scale, whose variation is used to probe the size of missing logarithmic higher-order corrections in eq. (2.4). We refer the reader to ref. [37] for further details. Though eq. (2.4) is valid for inclusive observables, the RadISH formalism can be systematically extended to any transverse observable in colour-singlet production, as stressed in ref. [37]. It is particularly instructive to consider the case of jet-veto resummation, which is currently available at NNLL [46-49, 111, 112]. At this logarithmic accuracy the resummation must include an additional clustering correction since the jet algorithm can cluster two independent emissions close in the pseudo-rapidity η or in the azimuthal angle φ. Moreover, the resummation must account for a correlated correction, which amends the inclusive treatment of the correlated squared amplitude for two emissions [113], accounting for configurations where the two correlated emissions are not clustered in the same jet [47]. The analytical result for the NNLL resummation reads [47] where the functions g 1 , g 2 and g 3 are the same appearing in the radiator of eq. (2.7) for the colour singlet p T , and now L = ln(Q/p J T ). For generalised k t algorithms [114][115][116][117][118] the clustering and the correlated corrections at NNLL accuracy in the limit of small jet radius R read 8 where C is C F = 4/3 or C A = 3 for incoming quarks or incoming gluons, respectively.
The inclusion of small-R resummation at LL R accuracy was found to have, for gluoninitiated processes, a very moderate effect on the result for values of R typically used in phenomenological applications [50]. In this work we neglect these effects, under the assumption that their impact is negligible. Further studies in this direction are beyond the scope of this paper.
We can now recast the resummation for p J T in the RadISH language as Here where ∆y ab and ∆φ ab are the difference in rapidity and in azimuthal angle between two emissions a and b, and C is the colour factor associated with the incoming hard leg . The function C is defined as the ratio of the correlated part of the double-soft squared amplitude and the product of the two single-soft squared amplitudes, cfr. eq. (30) of ref. [51]. At NNLL accuracy the formulations of eqs. (2.8) and (2.11) are equivalent, the only difference being in the treatment of subleading terms. The jet veto resummation of eq. (2.8) has the advantage that it is fully analytic, thereby allowing for a simple and fast implementation in a numerical code. The RadISH formulation, on the other hand, features less compact formulae which need to be evaluated numerically. For the single-differential jetveto resummation we employ the implementation of the analytic formulae in eq. (2.8) since their evaluation is faster. By casting the jet-veto resummation in the RadISH formulation, however, one gains a more differential description of the radiation than the one provided by eq. (2.8). This fact can be exploited to formulate the joint resummation of logarithms of p T and p J T within this framework by noticing that the two observables share the same Sudakov radiator R NNLL [47] at NNLL accuracy. The double-differential resummation is then achieved by supplementing the phase space constraint for the inclusive p T resummation of eq. (2.4) with a veto requirement, and by adding the clustering and correlated corrections that we described above. The interested reader can find the resulting formulae in the supplementary material of ref. [51].
The RadISH code implements the above formulae for the resummation of transverse observables using Monte Carlo methods. We refer the reader to section 4.3 of ref. [37] for details of the Monte Carlo evaluation and on event generation in the RadISH code.

Matching of resummation and fixed-order predictions
The calculation of a physical prediction for a general observable v across its entire differential spectrum requires a consistent matching between the fixed-order distribution valid in the hard region (large v) and the resummed result valid in the soft/collinear region (small v). This implies that, on the one hand, resummation effects have to vanish at large v, while, on the other hand, the fixed-order contribution should vanish at small v. In order to suppress resummation effects at large v, we map the limit k t1 → Q, where the logarithms vanish, onto k t1 → ∞ by introducing modified logarithms where p is a positive real parameter. Its value is chosen such that the resummed component decreases faster than the fixed-order spectrum for v 1. As a consequence, the logarithms in the Sudakov radiator (2.7), its derivatives and the luminosity factors have to be replaced byL. For consistency, eqs. (2.4) and (2.11) are supplemented by the following Jacobian in accordance with the replacement in eq. (2.15): This prescription leaves the Θ functions in eq. (2.4) and eq. (2.11) unchanged and modifies the final result by introducing power corrections in (Q/k t1 ) p beyond the nominal accuracy.
In order to perform the matching of the resummed calculation in the RadISH formalism to the fixed-order prediction, it is useful to introduce the cumulative fixed-order cross section, since the resummation is defined at this level, where dσ F NNLO is the fully differential NNLO cross section of eq. (2.1), and σ F NNLO is the NNLO cross section integrated over radiation. Note that in all expressions throughout this section we have dropped the explicit dependence on the Born phase space Φ B for the sake of brevity, which shall be understood implicitly for all cross sections. Thus, the cumulative NNLO cross section in eq. (2.17) is fully differential in the Born kinematics such that arbitrary fiducial cuts on the colour-singlet final state can be applied. We stress that the integral from 0 to v of dσ F NNLO is well-defined since all IR divergences have been canceled through the q T subtraction procedure. For brevity, we also drop the superscript F referring to a general colour-singlet final state from the NNLO cross section in the following equations.
We recall that the NNLO accuracy as denoted in eq. (2.17) applies to cumulative (integrated) cross sections, whereas the differential spectrum at large v is only NLO accurate. In the next section we will show results at the cumulative and at the differential level; in both cases, we shall label the fixed-order accuracy according to the accuracy defined at the cumulative level. 9 There is some freedom when defining a procedure to match resummation and fixed-order predictions: At a given perturbative order various schemes can be defined that differ from one another only by subleading terms, i.e. beyond the formal accuracy of the calculation. In the Matrix+RadISH interface we offer the possibility to choose between two different 9 Note that this convention is different from that used in ref. [37,39,40]. matching schemes to assess the associated uncertainties (see appendix A.4.1). The first scheme is a customary additive scheme, which at NNLO+NNLL is defined as Here [. . .] N k LO indicates the expansion of the expression inside the bracket truncated at N k LO, such that the second term is the expansion of the resummed cross section σ NNLL (v) up to NNLO (i.e. O(α 2 s )). It subtracts the logarithmically enhanced contributions at small v from the fixed-order component, thereby turning it finite and avoiding a double counting between the first and the third term. The second scheme is a multiplicative scheme, as formulated in refs. [39,119]. The quantity σ asym. NNLL is defined as the asymptotic (v 1) limit of the resummed cross section This prescription ensures that in the limit v → 0 eq. (2.19) reduces to the resummed prediction, and that for v 1 it reproduces the fixed-order result. The main difference between the two procedures is that the multiplicative approach is more robust against numerical instabilities at very low v, since a potential miscancellation of the logarithmic terms between the NNLO result and the expansion is suppressed by the resummation factor σ NNLL (v). Analogous formulae can be derived at NNLO+N 3 LL and at NLO+NLL accuracy. The detailed matching formulae for the multiplicative matching scheme are reported in appendix A of ref. [39]. We finally note that in both matching schemes the cumulative cross section at v → ∞ tends to σ NNLO and by construction the differential distribution fulfils the unitarity constraint, i.e. its integral yields the NNLO cross section.
The matching procedures immediately generalise to the double-differential case by rewriting the equations with the double-cumulative cross section σ(v 1 , v 2 ), obtained by integrating over the double-differential distribution dσ(v 1 , v 2 )/dv 1 dv 2 . In particular, the multiplicative scheme used for the double-differential case is We observe that by using the multiplicative scheme (2.21) one automatically recovers the Unless explicitly stated otherwise, multiplicative matching is used for all results presented in this paper. We recall that the matched cross sections in eqs. (2.18), (2.19) and (2.21) shall be understood as being fully differential in the Born level momenta, allowing us to apply arbitrary fiducial cuts on the colour-singlet final states.
The multiplicative scheme has one further advantage: when matching the resummation at NNLL to the fixed-order prediction at NNLO using eq. (2.19) or eq. (2.21), the terms that are constant in v at O(α 2 s ) are automatically included in the resummation, since the fixedorder contribution multiplies the resummed cross section. This procedure correctly resums the whole tower of α n s L 2n−4 contributions, which are formally part of the N 3 LL correction, in the matched cross section. This accuracy is often called 'primed' accuracy, and the predictions are sometimes denoted as NNLO+NNLL . Although all our NNLO+NNLL predictions have that accuracy, we refrain from explicitly adopting the primed notation in the remainder of this paper and assume it to be understood. In particular, the constant terms at O(α 2 s ) contain the two-loop virtual corrections and the second-order collinear coefficient functions. In an additive matching, on the other hand, the terms α n s L 2n−4 are included only up to n = 2 through the fixed-order contribution, so that they simply amount to a constant shift in the matched cross section at the cumulative level. Since the constant terms are still unknown analytically for jet-veto resummation and the joint resummation of p T and p J T , if we were to use the additive matching scheme, the additional logarithmic corrections would not be included.

Numerical implementation and validation
In the following, we briefly discuss some details of the practical implementation and the numerical validation of the Matrix+RadISH results. All fixed-order ingredients are evaluated by Matrix, which provides the cumulative distributions at NLO and at NNLO. Matrix also generates the Born level phase space used to integrate the resummed component, and for a given kinematics it computes the process-dependent Born matrix element as well as the hard-virtual corrections at NLO and at NNLO. For each Born event RadISH then produces the initial-state radiation using the numerical algorithm described in section 2.2 to perform the resummation of large logarithmic contributions. In parallel, RadISH also evaluates the expansion and the asymptotic limit of the resummed cross section entering eqs. (2.18) and (2.19). 10 Furthermore, we perform on-the-fly variations of the renormalisation, the factorisation and the resummation scale for each Born event. After all ingredients have been integrated to a sufficient numerical precision, the matching is performed according to either eq. (2.18) or eq. (2.19) as a post-processing step of the calculation.
Our calculations have been extensively validated by performing a variety of tests: for 2 → 1 processes, we have compared our resummed results against those obtained independently with the standalone version of the RadISH code where the required perturbative ingredients for Higgs boson and Drell-Yan production are implemented, finding full agreement up to numerical uncertainties for the resummed predictions. We have further checked for these processes that we find also full agreement at the level of matched cross sections calculated with Matrix+RadISH and those obtained by matching RadISH distributions with fixed-order predictions from MCFM [120,121].
A very powerful check of the robustness of our predictions for more complex coloursinglet processes can be achieved by comparing the expansion of the resummed cross section to the fixed-order result in the small-v limit, which we have performed for several 2 → 1 and MATRIX+RadISH (fiducial-noJV)  2 → 2 processes. At NLL (NNLL) accuracy, the resummation predicts the logarithmically enhanced contributions appearing in the differential fixed-order distribution at small v up to order α s (α 2 s ). The constant terms at the level of the NLO (NNLO) cumulative cross section, on the other hand, are not contained in the expansion up to α s (α 2 s ) of the NLL (NNLL) cross section, such that σ (N)NLO and σ (N)NLL (N)NLO differ by a constant in the v → 0 limit. The correct constant terms are included only in the resummed expression at the next logarithmic order (or when considering the respective primed accuracies), i.e. the difference σ NLO and [σ NNLL ] NLO as well as σ NNLO and [σ N 3 LL ] NNLO tend to zero in the small-v limit.
Thus, in order to perform a non-trivial check on the validity of our calculation we consider W + W − production and study the small-p W W T behaviour of the differential p W W T distribution in the left plot of figure 1, where p W W T is the transverse momentum of the W + W − pair. In the upper (lower) panel of the right plot of figure 1 we show the corresponding cumulative (differential) cross section as a function of ln(p W W T /GeV). These results are obtained with the settings described in section 3.2.1 in the fiducial phase space defined as "fiducial-noJV". 11 By using dedicated high-statistics runs we were able to push 11 Note that the dynamic scales and the fiducial cuts for the resummation and its expansion are based on the comparison down to remarkably low values of p W W T = 0.01 GeV with statistical uncertainties at the permille level. The right plot of figure 1 indicates an excellent agreement for the p W W T distribution between the fixed-order prediction and the respective logarithmic terms both at NLO and at NNLO. 12 The dip around 1-2 GeV is due to the fact that the NNLO spectrum and the NNLL expansion become negative. In particular the relative difference in the lower panel shows the striking cancellation between those terms at the few-permille level for very small p W W T . At the cumulative level, on the other hand, the upper frame of the right plot of figure 1 nicely confirms that the NLO−NLL exp difference yields a constant at small p W W T , while the fact that the NNLO−N 3 LL exp difference tends to zero validates our implementation of the constant terms at N 3 LL. Finally, in the lower panel of that figure we show the absolute differences of the differential distributions by taking the derivatives with respect to ln p W W T . Although the information might appear slightly redundant, by considering the absolute difference here we unambiguously prove that there are no potential residual differences also for the subleading logarithmic contributions.
We have performed similar checks for other 2 → 1 and 2 → 2 colour-singlet processes, with various choices of the renormalisation scale µ R , the factorisation scale µ F , and the resummation scale Q. In particular, we performed the exact same study for fully inclusive W + W − production, i.e. without fiducial cuts. In all cases we found an excellent agreement between the expansion of the resummed cross section and the fixed order result.
3 Resummed predictions for W + W − production at the LHC As a first application of the Matrix+RadISH framework introduced in section 2 we consider W + W − production at the LHC. This process plays a prominent role in the vast physics programme at the LHC since it has the largest cross section among all diboson production processes. Experimental measurements of the W + W − cross section at the Tevatron [122,123] and at the LHC [71,[124][125][126][127][128][129][130][131] provide crucial tests of the EW gauge sector of the SM and of the mechanism of EW symmetry breaking. Moreover, the production of W -boson pairs is an important probe of BSM physics. Since the dynamics of W + W − production is sensitive to the value of the trilinear gauge coupling already at the Born level, W + W − measurements put stringent bounds on the strength of anomalous trilinear gauge couplings in indirect searches for new physics [124,125,127,129]. W + W − production contributes also as irreducible background to direct searches for BSM particles decaying into leptons, missing energy and/or jets, and to Higgs boson measurements in the H → W + W − decay channel. Since the two neutrinos prevent the reconstruction of the full the Born level kinematics and thus correspond to those applied at fixed order only in the limit where extra QCD radiation is soft or collinear, and in particular in the p W W T → 0 limit. 12 We remind the reader that, when differential distributions are shown, the label NLO and NNLO refers to the O(αs) and the O(α 2 s ) prediction, respectively. event kinematics, an accurate theoretical description of the W + W − final state is essential to enhance the experimental sensitivity in BSM searches and Higgs boson measurements. Moreover, W + W − analyses generally apply a rather stringent cut on the jet activity in order to deplete the signal contamination due to top-quark backgrounds. Such a jet veto introduces large logarithmic contributions in the calculation of the fiducial cross section, which challenges the validity of fixed-order predictions and induces an additional uncertainty in the extrapolation to the inclusive W + W − cross section.
As a consequence, W + W − production has received much attention in the past years in a joint effort to reduce theoretical uncertainties in order to match the increasing precision of the experimental data. Next-to-leading order (NLO) QCD [132,133] predictions for on-shell W bosons have been known for years, and leptonic W boson decays were included in refs. [134][135][136][137]. Due to large O(α s ) corrections, a substantial effort has been put into calculating contributions at O(α 2 s ). The simplest O(α 2 s ) correction is the loopinduced gg → W + W − subprocess, which is enhanced by the gluon luminosities. Leading order (LO) predictions for the loop-induced gluon fusion channel were studied in refs. [137][138][139][140][141][142]. The full NNLO QCD corrections were calculated first for the inclusive W + W − cross section in the on-shell approximation in ref. [16], and were later advanced to the fully differential level including leptonic decays in ref. [17]. Contrary to the common lore, the quark-initiated contributions were found to be the dominant NNLO QCD corrections at the inclusive level, with the loop-induced gg-initiated channel contributing only ∼ 30% of the full O(α 2 s ) correction. Various efforts have been made to go beyond NNLO QCD accuracy: NLO EW corrections were calculated for stable W bosons [143][144][145], and with a consistent treatment of the leptonic decays [78,146]. The presumably dominant O(α 3 s ) corrections are known from the calculation of the loop-induced gluon fusion contribution at NLO QCD [77,147,148]. The first combination of NNLO QCD predictions with NLO EW corrections was presented very recently in ref. [78], and ref. [77] performed their further combination with NLO QCD corrections to the loop-induced gg channel, yielding the best fixed-order prediction for the W + W − process to date.
Analytic resummation approaches for different observables have been subject to various studies of W + W − production. The transverse-momentum spectrum of on-shell W + W − pairs was calculated at NNLO+NNLL accuracy in ref. [55], and the resummation of jet veto logarithms at NNLO+NNLL accuracy was performed in ref. [59]. In particular, the proper modelling of jet-vetoed W + W − cross sections has attracted much attention in the theory community [58,59,149,150], which has shown that higher-order corrections in both the perturbative and the logarithmic expansion are crucial to obtain reliable predictions and uncertainty estimates in presence of a jet veto cut in the fiducial phase space.

Outline of the calculation
We consider the process where the charged final-state leptons and the corresponding (anti-)neutrinos have different flavours ( = ). All resonant and non-resonant contributions to this process are accounted for, including off-shell effects and spin correlations, employing the complex-mass scheme [151] without any resonance approximation. We have full control on the momenta of the final-state leptons, which allows us to evaluate resummed cross sections in presence of fiducial selection criteria as employed by the experiments. Our calculation applies to an arbitrary combination of (massless) leptonic flavours, , ∈ {e, µ}, and in order to compare against experimental data, we evaluate the process pp → + − ν ν + X with , = e or µ and = . For the sake of simplicity, this process will be denoted as W + W − production in the following.
In figure 2 (a-c) we show representative Feynman diagrams at LO for W + W − production. They are driven by quark annihilation in the initial state and include t-channel W + W − topologies (panel a), s-channel W + W − topologies (panel b), and s-channel Drell-Yan-type topologies (panel c). Figure 2 (d), on the other hand, shows a loop-induced diagram that is driven by gluons in the initial state, which enters the cross section at O(α 2 s ) and is part of the NNLO corrections. In an NNLO calculation the loop-induced gluon fusion contribution is, however, effectively only LO accurate and has Born kinematics. Without associated QCD radiation at fixed order, it contributes trivially to the differential observables we are considering in this paper. At the resummed level, the LL resummation of this contribution starts at O(α 3 s L 2 ), and it is thus of the same size as N 3 LL corrections to the qq channel. In principle the resummation of the loop-induced gg channel can be included alongside the resummation of the qq channel. 13 However, in terms of their respective resummation both contributions can be treated completely independently, and we refrain from considering the loop-induced gluon fusion contribution in what follows, since we reckon that its proper treatment requires to go beyond an effective accuracy of LO+LL. In fact, NLO QCD corrections to the loop-induced gluon fusion contribution have already been evaluated within the Matrix framework for both ZZ [76] and W + W − [77] production. They are strongly enhanced by the large gluon luminosity, and contributing at O(α 3 s ) they constitute the dominant N 3 LO corrections to the cross sections of these two processes. These corrections involve diagrams with real QCD radiation that would directly contribute to the resummed observables under consideration. Therefore, a matching of the NLO QCD cross section with NNLL resummation for the loop-induced gluon fusion contribution would be a useful addition to the predictions presented below, which could be achieved within Matrix+RadISH, but it is beyond the scope of the present paper and left for future studies.
The calculation of higher-order corrections to the W + W − production process in QCD perturbation theory is affected by a subtle interplay with contributions stemming from the production of off-shell top quarks, which mix through t → W b decays [16,17,152,153] with real-radiation diagrams involving final-state bottom quarks. Due to the large cross section of top-quark processes at the LHC, these contributions induce a sizeable contamination of the W + W − cross section. In order to deal with this problem, we employ the four-flavour scheme (4FS), where bottom quarks are treated as massive and do not appear as initial-state particles. The bottom-quark mass renders partonic subprocesses with bottom quarks in the final state separately finite. We then evaluate top-free W + W − cross sections by omitting all subprocesses with real bottom-quark emissions, which in turn are considered to be part of the (off-shell) top-pair background. An alternative definition of the top-free W + W − cross section can be obtained in the five-flavour scheme (5FS), where bottom quarks are treated as massless, by exploiting the scaling behaviour of the cross section with the top-quark width, which has been explained in detail in refs. [16,17]. As it has been shown there, the top-free W + W − cross sections resulting from the 4FS and 5FS prescriptions agree within 1% − 2%, both at the inclusive level and with fiducial cuts. Consequently, we exploit the simpler 4FS prescription throughout this paper. We note that this approach requires the use of consistent sets of parton distribution functions (PDFs) with n f = 4 light parton flavours.
We recall that our implementation allows us to obtain resummed predictions for different transverse observables in W + W − production at the LHC. We calculate the transversemomentum spectrum of W + W − pairs at NNLO+N 3 LL accuracy as well as the transversemomentum spectrum of the leading jet (and equivalently the jet-vetoed cross section) at NNLO+NNLL accuracy. Even more notably, we also perform the simultaneous resummation of both observables at the double-differential level at NNLO+NNLL accuracy, which allows us to evaluate the spectrum of one of the two observables in presence of a veto on the other one. For all differential observables presented here it is the first time that such accuracies are achieved for a nontrivial process, i.e. beyond 2 → 1 scattering like Higgs boson production or the Drell-Yan process. In addition, arbitrary fiducial selection criteria can be applied in phase space of the leptonic final states.

Setup
We present resummed predictions for various observables in pp → + − ν ν +X production with , = e or µ and = at the LHC with √ s = 13 TeV. We employ the the G µ scheme to evaluate the EW coupling α = √ 2 G F m 2 W 1 − m 2 W /m 2 Z /π and the mixing angle cos    [155], respectively, by exploiting the Lhapdf interface [156] with the corresponding values of the strong coupling. Jets are defined according to the anti-k t algorithm [118] with a jet radius of R = 0.4 and no rapidity requirements. We employ the following dynamical setting for the central factorisation and renormalisation scales, and set the central resummation scale to where M W W is the invariant mass of the W + W − pair and p W W T is its transverse momentum. Perturbative uncertainties are estimated by performing seven-point scale variations of µ F and µ R around µ 0 by a factor of two in either direction, keeping 1/2 ≤ µ F /µ R ≤ 2 for Q = Q 0 , and by varying Q by a factor of two around Q 0 in either direction for µ F = µ R = µ 0 . The total scale uncertainty is evaluated as the envelope of the resulting nine variations. For the exponent of the modified logarithm in eq (2.15) we choose a value of p = 4 when evaluating the resummed p W W T spectra and p = 5 for the jet-vetoed cross section as well as the p J T distribution. We have explicilty checked that the dependence of the results on variations of p is negligible. Moreover, no non-perturbative corrections are included in our results.
In table 1 we summarise the set of fiducial cuts used in our study. Those cuts are the same as in the fiducial selection of the 13 TeV analysis of ATLAS in ref. [71], with the only exception that we do not apply any rapidity requirement on the jets since this would generate logarithmic corrections of non-global nature, which would spoil the formal accuracy  of our resummed results. 14 Furthermore, we introduce two notations: fiducial-JV will be used to refer to the full set of fiducial cuts in table 1, while fiducial-noJV denotes the same fiducial setup, but without any restriction on the jet activity.

Jet-vetoed cross section
We consider the cumulative cross section in the fiducial phase space as a function of a veto on the transverse momentum of the leading jet (p J,veto T ) as generically defined in eq. (2.2) where we identify the upper bound of the integral with the jet-veto cut p J,veto T , i.e. (3.4) The corresponding jet-veto efficiency is defined as where σ(p J,veto T → ∞) is the integrated cross section in the fiducial-noJV phase space. 14 We have checked that the impact of neglecting the rapidity cut on jets is at the permille level for the fiducial cross section.
In figure 3 we show results for the jet-vetoed cross section and the jet-veto efficiency as a function of p J,veto T . In the case of the cumulative cross section, we compare our predictions with the cross-section measurements from the ATLAS experiment at 13 TeV [71]. Since our resummed predictions do not include the loop-induced gluon fusion contribution, we have subtracted its LO prediction from the experimental results to facilitate a meaningful comparison. We observe that the NNLO result is in remarkable agreement with the NNLO+NNLL prediction down to jet-veto cuts of p J,veto T 15 GeV. Below ∼ 10 GeV the fixed-order result becomes unphysical, and the central prediction eventually turns negative. The fixed-order scale uncertainty band vastly underestimates missing higher-order corrections since this region is dominated by large logarithms. The uncertainties of the resummed results increase below 10 GeV, and by comparing NNLO+NNLL and NLO+NLL results it is clear that higher-order corrections in both the fixed-order and the logarithmic expansion are substantial and mandatory to provide reliable scale uncertainties, reduced to the fewpercent level. 15 Due to the relatively large experimental uncertainties, all predictions are in reasonably good agreement with the data. To better appreciate resummation effects in the comparison with data, it would be required to push the measurements to much lower jet-veto cuts.
Looking at the jet-veto efficiency in the right plot of figure 3, at 35 GeV, which is the value of the jet veto used in the fiducial phase space definition, the efficiency is 65%, and agreement between the NNLO and the NNLO+NNLL results is at the few permille level. Scale uncertainties for the efficiencies are calculated by considering fully correlated scale variations between the numerator and denominator of eq. (3.5). The scale uncertainty at NNLO+NNLL is about 5% at p T ∼ 5 GeV, and it decreases towards larger values of the jet veto, being about 2% at p T ∼ 35 GeV. The inclusion of the higher-order corrections reduces significantly the perturbative uncertainties, as illustrated by the NLO+NLL and the NNLO+NNLL bands.

Differential distributions
We now move to differential results in the fiducial-noJV phase space. In figure 4  20 GeV, where the NNLO distribution is marred by large logarithmic contributions and starts diverging. The NNLO+NNLL curve is instead well-behaved, and it has perturbative uncertainties of only a few percent down to p J T ∼ 5 GeV. Below that value, the uncertainty band becomes wider, reaching up to 20% in the first two bins. We note, however, that these regions are likely beyond experimental reach. In the tail of the NNLO+NNLL prediction the uncertainties gradually increase from about 5% at p J T = 50 GeV to about 10% at p J T = 400 GeV. The 15 We note that these results resemble closely the comparison between NNLO, NNLOPS and MiNLO in figure 6 of ref. [79].  Sudakov peak of the resummed p J T spectrum is around 4-5 GeV both at NLO+NLL and at NNLO+NNLL.
In figure 5 we show the transverse-momentum spectrum of the W + W − pair at NNLO as well as matched predictions at NLO+NLL and at NNLO+N 3 LL. For p W W T 20 GeV, the NNLO result starts diverging and thus loses predictivity. In that region only the matched results are well-behaved since the resummation of large logarithmic contributions becomes indispensable to obtain a physical prediction. The peak of the NNLO+N 3 LL spectrum is around 5-6 GeV, and it is shifted with respect to the NLO+NLL prediction by about 1 GeV, which peaks around 4-5 GeV. 16   To further study resummation effects in the region p W W T 50 GeV, we compare the result at NNLO+N 3 LL to the one at NNLO+NNLL in the left plot of figure 6. The effect of the N 3 LL corrections on the central prediction is quite sizeable, especially below 10 GeV, where it reaches almost 10%. The uncertainty bands, however, largely overlap, with the central prediction at NNLO+N 3 LL being fully contained within the NNLO+NNLL band. The inclusion of NNLO+N 3 LL corrections reduces the scale uncertainties by a bit less than a factor of two.
To illustrate uncertainties due to higher-order effects of the NNLO+N 3 LL prediction in the small-p W W T region, we compare the results for two different matching schemes, defined in eqs. (2.18) and (2.19), in the right plot of figure 6. At this accuracy, the two predictions contain the same ingredients and are compared on equal footing. We observe an excellent agreement between the two prescriptions, which indicates that our predictions exhibit a very mild dependence on the choice of the matching scheme. Only at very small transverse momenta (p W W T 2 GeV) we observe minor differences between the multiplicative result and the additive result due to the higher sensitivity of the additive matching to the exact cancellation between the fixed-order cross section and the expansion. The advantage of the multiplicative scheme is confirmed by the fact that the multiplicative matching is in  perfect agreement with the pure N 3 LL result, as it should be at very small transverse momenta, whereas the additive result is slightly different. Since the cancellation in the additive matching prescription is numerically challenging in this region, dedicated runs as those performed in section 2.4 would be required to achieve a more stable result.
We continue our studies by considering the transverse-momentum spectrum of the W + W − system in presence of a jet veto. To this end, we perform the joint resummation of large logarithmic terms in both p W W T and p J T . In the left plot of figure 7, we compare resummed predictions for the p W W T spectrum in the fiducial-JV phase space at NLO+NLL and NNLO+NNLL accuracies to the NNLO result. By exploiting the multiplicative matching scheme at the double-differential level, defined in eq. (2.21), we include the constant terms in the resummed prediction, and the integral of the NLO+NLL (NNLO+NNLL) p W W T spectrum yields the NLO+NLL (NNLO+NNLL) jet-vetoed cross section. The NNLO curve develops a perturbative instability [158] right at p W W T = 35 GeV, which corresponds to the value of the jet-veto cut. This instability is caused by an incomplete cancellation between soft contributions in the real and virtual amplitudes, and it leads to an integrable logarithmic divergence at the threshold. Since p W W T = p J T holds at LO, and consequently p W W T < 35 GeV, the p W W T region above the jet-veto cut is filled only by higher-order corrections, and the effective perturbative accuracy is reduced by one order. This is indicated  Figure 7. Transverse-momentum spectrum of the W + W − pair with a jet-veto requirement in the fiducial-JV phase space at NLO+NLL (blue, dotted), NNLO (green, dashed), and NNLO+NNLL (red, solid) accuracy (left plot), and comparing NNLO (green, dashed) and NNLO+NNLL (red, solid) predictions to the NNLOPS result (blue, dotted) of ref. [79] (right plot). To facilitate the comparison in the plot on the right, we reanalysed the NNLOPS events and we recalculated our predictions with the NNPDF3.0 NNLO PDF set [157].
by the widening of the uncertainty band of the NNLO prediction for p W W T > 35 GeV, which effectively becomes NLO accurate in that region.
The sensitivity to the perturbative instability at the threshold is largely cured by the NNLO+NNLL result, resumming a large part of the relevant logarithms. However, a slight sensitivity remains due to the fact that our approach resums Sudakov logarithms in the limit where p W W T and p J T are much smaller than the hard scale, while additional logarithmic terms contribute when hard jets are present. Nevertheless, the large differences between the NNLO+NNLL and NNLO results indicate the importance of resummation at the threshold. We observe even larger resummation effects for transverse momenta below 20 GeV, where the NNLO result becomes unphysical and resummation is mandatory to retain a reliable prediction. The resummed spectra at NLO+NLL and NNLO+NNLL are consistent with each other at small p W W T , with fully overlapping uncertainty bands up to p W W T ∼ 30 GeV. The theoretical uncertainty of the NNLO+NNLL prediction is at the few-percent level in that region and roughly a factor of two smaller than the NLO+NLL uncertainty band. It reaches O(10%) uncertainties only for p W W T 6 GeV, where the logarithmic corrections become larger. When approaching the perturbative instability, the differences between NLO+NLL and NNLO+NNLL progressively increase, with barely overlapping uncertainties below the threshold. At the threshold theory uncertainties of the NNLO+NNLL result reach 30%. The large differences between the resummed results and the widening of the uncertainty bands further indicate that additional logarithmic terms contribute at the jet-veto threshold, which are not resummed. Beyond threshold the NLO+NLL result becomes unreliable, being effectively only LO accurate, while the NNLO+NNLL prediction is effectively NLO accurate with an enlarged uncertainty band at the 20%-30% level.
In conclusion, our results indicate that the inclusion of higher-order corrections both in the fixed-order and the logarithmic expansion is particularly relevant for an accurate description of the p W W T distribution in presence of a jet veto. To further study resummation effects for this observable, we compare the results at NNLO+NNLL with the NNLOPS results of ref. [79] in figure 7. By and large, the two results are in reasonable agreement. At very small p W W T their uncertainty bands overlap. There is a gap between them for 5 GeV p W W T 15 GeV, because the NNLO+NNLL band almost vanishes at p W W T ∼ 10 GeV, and thus underestimates the actual uncertainties in that region, and also the NNLOPS band, which misses uncertainties related to the shower-starting scale, is somewhat underestimated. For p W W T 25 GeV, the NNLO+NNLL result provides the more accurate and therefore more reliable prediction. When approaching the threshold at 35 GeV, the two results agree within uncertainties due to the rather large NNLO+NNLL band. The NNLOPS result is smooth in that region and develops no perturbative instability. At higher p W W T values, the NNLOPS prediction becomes significantly larger than the NNLO+NNLL result since the shower generates additional QCD radiation, as discussed in ref. [79]. The comparison shows that in the low-p W W T region high-accuracy resummation is required for a precise prediction of the spectrum. For values of p W W T at and above the threshold, the NNLOPS result gives a more reliable description of the spectrum in the fiducial region, as it includes effects that are not included in other approaches, albeit with a limited formal accuracy.

Summary
In this paper we have introduced a general framework for the computation of accurate crosssection predictions including multi-differential resummation. By developing an interface between the codes Matrix and RadISH we have combined fully differential NNLO QCD predictions for 2 → 1 and 2 → 2 colour-singlet production processes with high-accuracy allorder resummation. In particular, Matrix+RadISH evaluates the transverse-momentum spectrum of colour singlets and the φ * η distribution for the Drell-Yan process up to NNLO+N 3 LL accuracy, the transverse-momentum spectrum of the leading jet and jet-veto resummation up to NNLO+NNLL accuracy, as well as the joint resummation of the colour-singlet and the leading-jet transverse momentum up to NNLO+NNLL. Thereby we provide a powerful and versatile parton-level Monte Carlo generator allowing for an accurate description of transverse observables in colour-singlet production processes with arbitrary cuts on the Born kinematics. This framework can be extended to describe any observable differential in the Born kinematics including resummation effects by exploiting a suitable scheme for the kinematic recoil within the resummation. Moreover, the framework facilitates the combination of future developments within Matrix and RadISH, such as advancements towards processes beyond colour-singlet production.
As a first phenomenological application we have studied the production of W + W − pairs at the LHC. More precisely, the full leptonic final state of two charged different-flavour leptons and two neutrinos has been considered, including spin correlations, interferences and off-shell effects. We have presented resummed results at 13 TeV for several kinematic distributions in presence of fiducial cuts on the leptonic final states and the associated QCD radiation. The inclusion of higher-order corrections in both the logarithmic and the fixedorder expansion turns out to be essential to achieve a precise description of differential observables at the few-percent level, as demanded by the experimental analyses. The accurate modelling of the cross section in presence of a veto against hard QCD radiation is particularly important in that respect since a jet veto is commonly employed in W + W − measurements to suppress top-quark backgrounds. Our NNLO+NNLL result yields residual uncertainties at the few-percent level, we have shown that NNLO predictions are reliable down to jet-veto cuts of roughly 15 GeV, and we find agreement within one standard deviation with the cross section measured by ATLAS as a function of the jet-veto cut between 30 GeV and 60 GeV.
For the differential spectra in the transverse momentum of the leading jet and the W + W − pair, which we have evaluated up to NNLO+NNLL and NNLO+N 3 LL accuracy, respectively, we obtain scale uncertainties that are generally below 5%. Furthermore, our results highlight the importance of high-accuracy resummation in the region of the Sudakov peak, and of the O(α 2 s ) corrections in the tail of the distribution. Indeed, scale variations at O(α s ) significantly underestimate the actual size of O(α 2 s ) corrections at large transverse momenta. At small transverse momenta, N 3 LL corrections to the W + W − transversemomentum distribution are still quite sizable, with differences of about 5% -10% when comparing NNLO+N 3 LL to NNLO+NNLL.
The matching of the resummed cross section, valid in the soft-collinear region, and the fixed-order predictions, valid in the hard region, can be achieved in different ways. We have shown that a multiplicative matching procedure, which at NNLL resums additional logarithmic contributions originating from the NNLO terms that are constant in the resummation variable, has the further advantage of being numerically more stable at small transverse momenta than the additive matching procedure. Furthermore, we found that at NNLO+N 3 LL there is essentially no dependence of the W + W − transverse-momentum spectrum on the choice of the matching scheme.
Finally, we have studied the transverse-momentum distribution of W + W − pairs in presence of a 35 GeV jet-veto requirement by simultaneously resumming both classes of logarithmic terms. At fixed order no realistic description of this observable can be obtained. Resummed results at NLO+NLL and NNLO+NNLL indicate good perturbative convergence at small transverse momenta (below ∼ 30 GeV), and the inclusion of the NNLO+NNLL corrections decreases scale uncertainties by roughly a factor of two in that regime. The most delicate region is the threshold where the W + W − transverse momentum is close to the value of the jet-veto cut since this region is subject to a perturbative instability at fixed order in α s . The joint resummation of jet-veto logarithms and logarithmic terms in the W + W − transverse momentum improves the stability of the spectrum in the threshold region by resumming part of the relevant logarithms. We have further compared our NNLO+NNLL results to NNLOPS predictions, which corroborates the importance of NNLL resummation for an accurate modeling at small transverse momenta. In the threshold region the NNLOPS result is smooth, and we found it to be in reasonable agreement with NNLO+NNLL predictions, well within the respective scale uncertainties.
In all presented results we refrained from including the loop-induced gluon fusion contribution, which is formally part of the NNLO corrections, because it is effectively only LO accurate and Born-like. Although at fixed order it contributes trivially to all observables we have considered, i.e. as a constant shift in the cross sections, at the resummed level it is of the same size as N 3 LL corrections to the quark channel. Since their resummation can be considered completely independently, we leave a proper treatment of the loop-induced gluon fusion contribution, and specifically the combination of its NLO corrections with NNLL resummation, for future work.
We reckon that the predictions presented in this paper for the specific case of W + W − production as well as the Matrix+RadISH framework in general will be a very useful addition to current fixed-order and parton-shower predictions and tools. 17 With this framework we hope to advance the sensitivity to transverse observables in colour-singlet processes for precision measurements and new-physics searches at the LHC and future pp colliders.

A.1 Compilation and setup of a process
We assume that the MATRIX+RadISH_v1.0.0.tar.gz package has been downloaded and extracted, and that a working installation of Lhapdf [156] is present on the machine such that the shell command lhapdf-config is recognised when it is run in a terminal. We also require a working Python 2.7 installation. 18 If these requirement have been met, the command $ ./matrix --radish executed inside the MATRIX+RadISH_v1.0.0 folder will launch the Matrix shell with the RadISH interface enabled. 19 An interactive Python session for the compilation will appear, with instructions printed on the screen. This is followed by the usual steps in the Matrix shell to select a process via its ${process_id}, 20 e.g. for Drell-Yan production |===>> ppz01 and to agree with the terms to use Matrix+RadISH by typing (a few times)

|===>> y
This requires you to acknowledge the work of various groups that went into the computation performed by Matrix+RadISH for the present process by citing the references collected in the file CITATION.bib, which is provided with the results in every run. The compilation shell will then execute the following steps: • linking to Lhapdf [156]; • download and installation of OpenLoops [83][84][85] (skipped if already installed); • installation of CLN [159] (skipped if already installed); • installation of GiNaC [160] (skipped if already installed); • installation of Hoppet [161] (skipped if already installed); • installation of Chaplin [162] (skipped if already installed); • installation of RadISH [35,37,51] (skipped if already installed); • compilation of Matrix process (asked for recompilation if executable exists); • download of the relevant tree-level and one-loop amplitudes through OpenLoops (skipped if they already exist); • setting up of the process folder under the path run/${process_id}_MATRIX . 18 Currently the scripts necessary to run the Matrix shell are not yet compliant with Python 3 standards. 19 For a comprehensive list of the additional options of the ./matrix command see section 3.2 of ref. [69]. 20 The resummation is available for all colour-singlet processes contained in the public release of Matrix.
Once completed, the shell will exit and the process is ready to be run from the created process folder. As instructed on screen, enter that folder (cd run/${process_id}_MATRIX) and start a run by following the instructions in the next section.

A.2 Running a process in the interactive shell
Once the process has been compiled, we launch the running shell with the usual command $ ./bin/run_process inside the folder (cd run/${process_id}_MATRIX). With this interactive steering interface all run-related settings, inputs, and options will be handled just as in a standard Matrix run. 21 After typing a name for the run, e.g.

|===>> run_name_resummation
a list of several commands is printed on the screen. Before starting the run, the only differences with respect to an ordinary Matrix run are the inputs of the files parameter.dat and distribution.dat of the current run inside the folder input. We open them by typing

|===>> parameter |===>> distribution
and we modify the settings to turn on the resummation as discussed in appendix A.4. Once the input cards have been adjusted to obtain the desired results, the run is started through the usual command:

|===>> run
From now on, no human intervention is needed. Once the run is finished, the results from the run are collected in the respective folder result/${run_name_resummation}, which is printed on screen at the end of the run. The only difference, compared to an ordinary fixed-order Matrix run, after a run with resummation has been completed, is that there will be additional results inside the folder result/${run_name_resummation}. In particular, there will be the following information: • The same results as in a fixed-order run are saved, and the general structure is identical.
• If resummation for the respective orders has been turned on in the file parameter.dat, as discussed in appendix A.4.1, inside the folders NLO-run and NNLO-run a folder MATRIX+RadISH is present.
variable description path_to_radish Path to the RadISH installation; not required in most cases. path_to_hoppet Path to hoppet-config; not required in most cases. path_to_chaplin Path to the Chaplin installation; not required in most cases. Table 2. Additional parameters available the file MATRIX_configuration.
• The folder MATRIX+RadISH contains the resummed differential distributions (ending with NLL.dat, NNLL.dat or N3LL.dat), the corresponding fixed-order differential distributions (ending with NLO_QCD.dat or NNLO_QCD.dat), the matched distributions at differential level (ending with NLO+NLL.dat, NNLO+NNLL.dat or NNLO+N3LL.dat), which correspond to the best and final prediction, and, for completeness, the expansion of the resummed cross section at the relevant order (ending with exp_NLO_QCD.dat or exp_NNLO_QCD.dat).
• There is a further subfolder cumulant inside MATRIX+RadISH that contains the respective cumulative cross sections obtained from the differential distributions.

A.3 Additional settings in the configuration file
Before turning to physics-related and technical settings relevant for a Matrix+RadISH run in the next section, the additional parameters in the file MATRIX_configuration are discussed. The main file MATRIX_configuration can be found in the folder config inside the Matrix+RadISH main folder. This file is linked during each setup of a process (see appendix A.1) into the folder input of the respective process folder. The additional options for a Matrix+RadISH run in the file MATRIX_configuration are listed in table 2.
A.4 Additional settings in parameter.dat and distribution.dat The various parameters in the input files are described in section 4.1 of ref. [69]. Here we will discuss only the additional settings in the files parameter.dat and distribution.dat that are relevant when performing a run in Matrix+RadISH.

A.4.1 Settings in parameter.dat
All main parameters, related to the run itself or the behaviour of the code, are specified in the file parameter.dat. Most of them should be completely self-explanatory, and we will focus our discussion on the essential ones with a direct connection to the resummation. The settings are organized into certain groups. Here, we will limit the discussion to the groups which contain additional information when the Matrix+RadISH interface is enabled, for the sample case of Zγ production (where applicable) as in ref. [69]. dynamic_scale This parameter is present already in a fixed-order run to choose between the specified fixed renormalisation and factorisation scales (scale_ren/scale_fact) and dynamic ones, and we refer the reader to the information given section 5.1.1.2 of ref. [69] for details. In a resummed calculation it is important to only use sensible choices for dynamical scales, i.e. any choice that depends on QCD radiation must smoothly approach the scale at Born level in the limit where the QCD radiation becomes soft or collinear. In particular, this allows for any dynamical scale choice that only involves the Born level momenta of the color singlets.
dynamic_scale_res This parameter allows the user to choose between the specified fixed resummation scale (scale_res) and a dynamic one. The dynamic scale, however, is limited to the invariant mass of the colour singlet, which is also the recommended setting.
factor_central_scale_res A relative factor that multiplies the central resummation scale; particularly useful for a dynamic resummation scale. It is recommended to set it to a fraction of the invariant mass of the color singlet (setting dynamic_scale_res = 1) such as 0.25 or 0.5.
scale_variation_res Switch to turn on and off scale variations of the resummation scale. If set to 1, a variation up and down will be done while keeping the renormalisation and factorisation at their central values. The total uncertainty is calculated as the envelope of the renormalisation and factorisation scale variations (1, 7, or 9 variations) and of the additional two variations of the resummation scale, for a total of 3, 9, or 11 variations.
variation_factor_res This (integer) value determines by which factor with respect to the central resummation scale the scale variation is performed. add_RadISH_NLL Switch to turn on and off the NLL resummation through RadISH in the NLO run and the matching to NLO QCD. Available for all RadISH observables (settings of RadISH_observable below).
add_RadISH_NNLL Switch to turn on and off the NNLL resummation through RadISH in the NNLO run and the matching to NNLO QCD. Available for all RadISH observables (settings of RadISH_observable below). Can only be turned on if add_RadISH_N3LL below is turned off.
add_RadISH_N3LL Switch to turn on and off the N 3 LL resummation through RadISH in the NNLO run and the matching to NNLO QCD. Available only for transverse-momentum resummation of the color singlet (RadISH_observable = 1 below). Can only be turned on if add_RadISH_NNLL above is turned off.
loop_induced This switch is present already in a fixed-order run for certain processes (such as ZZ, W + W − , . . . ) where a loop-induced gg contribution enters at the NNLO. In a resummation (for now) the switch has to be set to 0, since the effectively LO-accurate gg contribution would trivially enter the resummed observables under consideration.
4. transverse momentum of the colour singlet with a jet veto (double-differential resummation); requires to set a jet veto in the settings for fiducial cuts in parameter.dat with the following mandatory settings: # Jet algorithm ... jet_R_definition = 0 # (0) pseudo -rapidity (1)  where only the value of the jet veto (${veto-value}) may be chosen freely between 1 and 1.e99, a typical setting would be for instance define_pT jet = 30., and the value of the jet radius (${jet-radius-value}) may be chose either to be 0.4 or 0.8; 22 5. transverse momentum of the leading jet with a veto on the transverse momentum of the color singlet (double-differential resummation); requires to set a veto on the transverse momentum of the color singlet via RadISH_pTsystem_veto between 1 and 1.e99, a typical setting would be for instance RadISH_pTsystem_veto = 30. ; and to set the jet definition in parameter.dat with the following mandatory settings: # Jet algorithm ... jet_R_definition = 0 # (0) pseudo -rapidity (1) rapidity jet_R = $ { jet -radius -value } # DeltaR ... # Jet cuts ... define_eta jet = 1. e99 # requirement on jet pseudo -rapidity ( upper cut ) define_y jet = 1. e99 # requirement on jet rapidity ( upper cut ) n_observed_min jet = 0 # minimal number of observed jets ( with cuts above ) n_observed_max jet = 99 # maximal number of observed jets ( with cuts above ) modlog_p This (integer) value sets the exponent p in the modified logarithms defined in eq. (2.15). Has to be chosen between 1 and 5. We recommend a setting of modlog_p = 4 for the transverse momentum of the color singlet and φ * η , and a setting of modlog_p = 5 for the transverse momentum of the leading jet.
matching_scheme Switch to choose between the additive matching scheme in eq. (2.18) or the multiplicative matching in eq. (2.19). We recommend using the multiplicative matching for the various reasons discussed in section 2.3. In order to compare results obtained with both schemes for the same run, first complete a run with one matching scheme; then 22 Please contact the authors if other values for the jet radius are desired.
choose the other matching scheme by modifying the file parameter.dat of the respective run inside the folder input and restart the combination (and matching) of the results with the command $ ./bin/run_process run_name_resummation --run_mode run_matching If enabled (default) in the file parameter.dat the previous results will be saved in a subfolder saved_result_XX inside the result folder of the respective run, where XX is an increasing number starting at 01.
number_of_events This (integer) value sets the number of events used for each Born event in the RadISH Monte Carlo implementation. We recommend to use a value between 2000 and 5000, and not lower than 1000 to avoid numerical artefacts. A further increase of this number may be useful for dedicated runs with very fine binning at small transverse momenta, but one must bear in mind that this will slow down the computation.
min_born_events_radish This (integer) value sets the minimum number of Born events for each parallel run of the resummation with RadISH. We recommend to use a value between 2000 and 5000, and not lower than 1000 to avoid numerical artefacts. In heavily parallelised runs, this ensures that each parallel run of the resummation with RadISH has a minimum number of Born events, so that the numerical integration of each of them converges sufficiently before combining all parallel runs.

A.4.2 Settings in distribution.dat
In the file distribution.dat the user defines histograms for distributions which are filled during the run. Each distribution is represented by one block of certain parameters. We refer the reader to section 5.1.3 of ref. [69] for details.
In order to perform a resummation run it is mandatory that there is at least one block with distributiontype = RadISH_observable. This special distributiontype automatically selects the correct observable according to the setting of the RadISH_observable inside the file parameter.dat of the respective run, described in the previous section. In the first version of Matrix+RadISH all other distributions, i.e. those with distributiontype different from RadISH_observable, are ignored for the calculation of the resummation and the matching. Their fixed-order results will be calculated and provided in the result folder as in an ordinary fixed-order run though. Thus, the resummation and matching are evaluated only for distribution blocks with distributiontype = RadISH_observable, which corresponds to the observable selected for the resummation.
The starting point, selected either via startpoint or via the first entry of edges, of any distribution block using distributiontype = RadISH_observable must be set to 0 in order to warrant the correct matching of fixed-order prediction and resummation. The user is free to add an arbitrary number of distribution blocks with distributiontype = RadISH_observable, provided that the distributionname of each distribution block is different. Below, we provide an example for the case of a regular and an irregular binning: We stress again that any additional fixed-order distribution alongside these histograms may be added as in a usual Matrix run, but that they will not be filled for the resummation and the matching of the calculation.