NNLO QCD+NLO EW with MATRIX+OpenLoops: precise predictions for vector-boson pair production

We present the first combination of NNLO QCD and NLO EW corrections for vector-boson pair production at the LHC. We consider all final states with two, three and four charged leptons, including resonant and non-resonant diagrams, spin correlations and off-shell effects. Detailed predictions are discussed for three representative channels corresponding to $W^+W^-$, $W^{\pm}Z$ and $ZZ$ production. Both QCD and EW corrections are very significant, and the details of their combination can play a crucial role to achieve the level of precision demanded by experimental analyses. In this context we point out nontrivial issues that arise at large transverse momenta, where the EW corrections are strongly enhanced by Sudakov logarithms and the QCD corrections can feature so-called giant K-factors. Our calculations have been carried out in the MATRIX+OpenLoops framework and can be extended to the production of an arbitrary colour singlet in hadronic collisions, provided that the required two-loop QCD amplitudes are available. Combined NNLO QCD and NLO EW predictions for the full set of massive diboson processes will be made publicly available in the next release of MATRIX and will be instrumental in advancing precision diboson studies and new-physics searches at the LHC and future hadron colliders.


Introduction
The production of a pair of massive vector bosons, W + W − , W ± Z and ZZ, plays a crucial role in various areas of the LHC physics programme. Experimental studies for this class of processes lead to precise tests of the Standard Model (SM) at energies that range from the electroweak (EW) scale up to the TeV regime. In particular, due to the high sensitivity to anomalous trilinear couplings, differential measurements at high transverse momenta probe the gauge symmetry structure of EW interactions and provide an excellent chance to detect indirect effects of physics beyond the Standard Model (BSM). Diboson final states are widely studied also in the context of direct BSM searches and in Higgs-boson measurements, where they are non-trivial irreducible backgrounds, as, for example, in the H → W + W − and H → ZZ analyses [1][2][3][4][5][6][7][8]. The increasing level of precision of experimental measurements calls for continuous improvements in the theoretical description of diboson production at the LHC. In particular, SM measurements of vector-boson pair production [9][10][11][12][13][14][15][16][17][18] are approaching percent-level precision, which can be matched on the theoretical side only by highest-order computations in both QCD and EW perturbation theory.
Leptonically decaying vector-boson pairs yield clean experimental signatures with charged leptons and missing energy. Theoretical predictions for on-shell W + W − [19], W ± Z [20] and ZZ [21,22] production, and their full leptonic final states [23][24][25][26], are available up to next-to-next-toleading order (NNLO) in QCD in the Matrix framework [27]. The effect of the resummation of the large logarithmic terms at small transverse momenta up to next-to-next-to-leading logarithmic accuracy (NNLL) has been included in ref. [28] in the same framework, while the matching of the NNLO calculation with parton showers has been presented in ref. [29]. For charge-neutral final states, loop-induced contributions from gluon fusion play an important role, and they are known at O(α 3 S ) [30][31][32].
In order to reach the level of precision required by present and future experimental analyses, higher-order QCD predictions need to be supplemented by the effect of EW corrections. In general, the dominant EW corrections arise from QED radiation effects in distributions of the final-state leptons, and from large Sudakov logarithms in the high-energy tails of kinematic distributions [33,34]. Calculations of diboson production at NLO EW have first been carried out for stable vector bosons [35][36][37][38], and subsequently extended to the case of off-shell vector bosons decaying into any final state with four [39,40] and three charged leptons [41], i.e. and ν with = or = , as well as final states with two charged leptons of different flavour plus missing energy [42], i.e. ν ν with = . Results at NLO QCD+EW accuracy for any final state with two charged leptons plus missing energy, i.e. ν ν and ν ν with = or = , have been presented in ref. [43]. For the case of on-shell W + W − production, EW two-loop Sudakov corrections have been studied up to NNLL accuracy [44].
In this paper, we present the first combination of NNLO QCD and NLO EW corrections for massive vector-boson pair production processes in their leptonic decay modes at the LHC. We consider all the W + W − , W ± Z and ZZ production processes with two, three and four finalstate charged leptons. Spin correlations, interferences and off-shell effects are included throughout. While QCD corrections have no effect on leptonic vector-boson decays, we note that the NLO EW corrections affect both the production of the virtual vector-boson pairs and their subsequent decays. Moreover, the off-shell pp → 4 fermion processes at hand involve also irreducible NLO EW effects that cannot be separately attributed to vector-boson pair production or decay.
The combination of QCD and EW corrections is discussed in detail. In particular, by comparing different combination prescriptions we gain insights into theoretical uncertainties associated with higher-order effects beyond NNLO QCD+NLO EW. In this context we devote particular attention to the tails of various kinematic distributions that feature giant QCD K-factors [37,45] together with large Sudakov-type EW corrections, resulting in sizeable higher-order uncertainties. As we will show, such uncertainties can be tamed by requiring that the two vector bosons are comparably hard in the high-p T regime. This can be achieved either via direct selection cuts on the vector bosons or through a jet veto.
The calculations presented in this paper have been carried out using Matrix+OpenLoops. Our implementation of the NNLO QCD and NLO EW corrections is completely general, and can be applied to the production of an arbitrary colour singlet in hadronic collisions, provided that the required two-loop QCD amplitudes are available. This paves the way to extend the various process libraries that are available in the current public version of Matrix [27] to NNLO QCD+NLO EW accuracy.
The paper is organised as follows. In section 2 we outline the calculation: We describe the employed methods and tools (section 2.1), the leptonic channels contributing to vector-boson pair processes (section 2.2) and the diagrams contributing to their QCD (section 2.3) and EW (section 2.4) corrections. In section 2.5 we discuss the phenomenologically important issue of giant K-factors, and in section 2.6 we present different procedures to combine QCD and EW corrections. Our numerical results for vector-boson pair production processes are presented in section 3, in the inclusive case (section 3.3) and with a veto against hard QCD radiation applied (section 3.4). Finally, in section 4 we draw our conclusions.
2 Vector-boson pair production in NNLO QCD and NLO EW

NNLO QCD+EW predictions with Matrix+OpenLoops
The NNLO QCD and NLO EW computations for vector-boson pair production processes presented in this paper are carried out within the Matrix+OpenLoops framework. The core of the Matrix framework [27] is the Monte Carlo program Munich 1 which includes a fully automated implementation of the dipole-subtraction method at NLO QCD [46,47] and NLO EW [43,[48][49][50][51], combined with an efficient multi-channel integration algorithm. The required tree-level and oneloop amplitudes, including spin and colour correlations for the subtraction of infrared divergences, are calculated with OpenLoops 2 [52]. This new version of the OpenLoops program supports the automated generation of matrix elements for any SM process at NLO QCD+NLO EW. Scattering amplitudes are built with a new numerical recursion that combines the original method of ref. [53] with the on-the-fly reduction technique of ref. [54]. To avoid numerical instabilities in critical phase space regions, this new algorithm is equipped with a sophisticated system that automatically detects the occurrence of spurious singularities in the reduction identities and circumvents them by using algebraic tricks in combination with analytic expansions up to any order in rank-two Gram determinants [52,54]. Residual instabilities are avoided in a systematic way by means of an efficient hybrid-precision system, which carries out certain critical operations in quadruple precision while using double precision for the bulk of the calculation. [52]. These new techniques improve the numerical stability of one-loop amplitudes in a very significant way, especially in the regions where one parton becomes soft or collinear. In OpenLoops 2 scalar integrals from Collier [55] and OneLOop [56] are used. Our NLO EW results for vector-boson pair processes have been validated systematically, starting from an in-detail comparison of all results in ref. [43] between Sherpa and Munich. All relevant EW one-loop amplitudes were validated between OpenLoops and Recola2 [57] at the level of phase space points. Further, in ref. [58] NLO EW amplitudes and cross sections for pp → + − + − and pp → + − ν ν were checked amongst various public tools [51,[59][60][61][62][63].
The required two-loop amplitudes for NNLO QCD calculations need to be implemented on a process-by-process basis. For the calculation of off-shell vector-boson pair processes, Matrix uses the two-loop qq → V V amplitudes of ref. [64] as provided by the VVamp package. The subtraction of IR singularities at NNLO is achieved through a fully general implementation of the q T -subtraction formalism [65]. In this context, Matrix implements an automated extrapolation procedure to calculate integrated cross sections in the limit in which the cutoff parameter in the q T -subtraction procedure goes to zero [27].
The q T -subtraction method and its implementation in Matrix were originally limited to coloursinglet final states, but have been recently extended to provide NNLO QCD predictions for topquark pair production as well [66,67]. Despite the fact that in this paper we focus on massive vector boson pairs, our implementation is completely general, and can directly be used for any other processes with a colourless final state, e.g. pp → V H, V V V , etc., where V refers to either on-shell vector bosons, V = {Z, W ± }, or off-shell leptonic final states, V = { + − , νν, −ν , + ν}. The only restriction is given by the availability of the corresponding two-loop amplitudes.

Overview of leptonic channels in vector-boson pair production
The present version of Matrix implements the full set of pp → 4 fermion processes that involve W + W − , W ± Z or ZZ resonances in all channels with final states containing two, three or four charged leptons. 2 The complete set of processes and corresponding acronyms used in this paper are listed in table 1. The various final states are categorised according to the number of charged leptons and further split into different-flavour (DF) and same-flavour (SF) channels according to contributing shorthand acronym process resonances in this paper Table 1. Complete list of diboson processes that are implemented in Matrix and will be upgraded to NNLO QCD+NLO EW accuracy in the forthcoming code release. The last column indicates the shorthands used for the three representative processes presented in this paper. In this table it is implicitly understood that = .
the flavours of charged leptons. Note that the channels 2l-SF-ZZ and 2l-SF-ZZWW yield identical experimental signatures and should thus be combined for phenomenological applications. 3 We have combined NNLO QCD with NLO EW corrections within the Matrix+OpenLoops framework for all processes listed in table 1. This implementation will be made publicly available with the forthcoming version of Matrix, and in this paper we illustrate the effect of NNLO QCD+ NLO EW corrections on the representative channels 2l-SF-ZZ, 2l-DF-WW and 3l-DF-WZ. For brevity, we will refer to these three channels as ZZ, W W and W Z production, respectively. As pointed out in the introduction, all relevant pp → 4 lepton matrix elements are computed exactly, i.e. without applying any resonance approximation. All Feynman diagrams with double-, singleand non-resonant topologies are consistently included at each perturbative order using the complexmass scheme [68]. Therefore off-shell effects, interferences and spin correlations are fully taken into account throughout.
In figure 1 we show representative LO Feynman diagrams for the selected ZZ, W W and W Z production processes. As illustrated in figure 2, diboson processes with charge-neutral final states, i.e. ZZ and W W production, involve additional photon-induced channels. In Ma-trix+OpenLoops the photon distribution function is treated on the same footing as the QCD parton densities. Thus, photon-induced channels enter at the same perturbative order as the usual qq channels, and both channels are supplemented by NLO EW corrections. This is important for a reliable description of certain phase space regions where photon-induced effects can be significantly enhanced by the opening of quark-photon channels at NLO EW.

Higher-order QCD corrections
For vector-boson pair production processes, higher-order QCD corrections have a sizeable impact. The NLO QCD corrections increase inclusive cross sections by 40-50% for ZZ and W W production and around 70-80% for W Z production [37,[69][70][71][72][73][74][75][76]. The large NLO effect for W Z production originates from an approximate radiation zero appearing in the leading helicity amplitude for W Z  production at LO [77], which is not present at higher orders. Also NNLO QCD corrections have a quite significant impact, at the level of 10% or more, on the various diboson production processes [19-21, 23-26, 78, 79]. Predictions at NLO QCD require the calculation of virtual and real-emission matrix elements, while NNLO QCD corrections involve double-virtual, real-virtual, and double-real contributions. Representative Feynman diagrams are displayed in figure 3 for the case of W + Z production. Similar diagrams contribute also to the other diboson processes. Only for ZZ production diagrams with triple vector-boson couplings are absent. In addition to the contributions illustrated in figure 3, W W and ZZ production involve also a loop-induced gluon-fusion channel that enters at O(α 2 S ), i.e. it is part of the NNLO QCD corrections. The contribution of this gg → V V channel to chargeneutral final states is quite sizeable. It has been computed to one order higher in perturbation theory [30,32,[80][81][82][83][84], which is assumed to be the dominant O(α 3 S ) correction to these processes. In the combination of NNLO QCD and NLO EW corrections presented in this paper, the gg → V V channels are included at O(α 2 S ) as part of the NNLO QCD corrections, i.e. neglecting O(α 3 S ) effects.

Higher-order EW corrections
The impact of NLO EW effects on inclusive cross sections is typically at the few-percent level and thus important in the context of high-precision studies. In kinematic distributions, EW corrections can be more sizeable. In particular, in the tails of distributions that probe high-energy scales Q M W , the EW corrections are enhanced by Sudakov logarithms [33,34] of the form α w log 2 Q 2 /M 2 W , where α w = g 2 w /(4π) denotes the SU(2) coupling strength. The size of EW Sudakov effects depends on Q as well as on the SU(2)×U(1) quantum numbers of the scattering particles. Such logarithmic effects are most pronounced in processes with (multiple) transversely NLO QCD: virtual polarised W and Z bosons, and in the case of ZZ, W Z and W W production they can lead to NLO EW corrections of several tens of percent at the TeV scale. EW corrections have been studied for various vector-boson pair production modes, treating the vector bosons on-shell [35,37,38], and more recently also fully off-shell [39][40][41][42][43]. Including off-shell and non-resonant effects is preferable, and they can play an especially important role in the tails of kinematic distributions [42].
The structure of NLO EW corrections to vector-boson pair production is illustrated in figure 4, where we show representative Feynman diagrams for the virtual and real corrections to W Z production. 4 The virtual corrections enter only through the qq channel and involve one-loop diagrams with various combinations of photons, Z/W ± -bosons, Higgs bosons, light fermions and heavy quarks in the loop. Real-emission contributions consist of a qq channel with an additional final-state photon, a qγ channel with an additional final-state quark and a correspondingqγ channel. The extra photon couples to any external or internal charged fermion or W boson.
In the case of ZZ and W W production, the presence of additional γγ → V V channels at LO gives rise to corresponding virtual and real contributions at NLO EW (not shown in figure 4). The real EW corrections to γγ → V V processes with V V = ZZ, W W involve γγ → V V γ channels as well as qγ → V V q andqγ → V Vq channels. As discussed in ref. [43], the qγ andqγ channels play NLO EW: virtual the twofold role of NLO EW corrections to qq → V V and γγ → V V . In particular, they cannot be uniquely assigned to one or the other LO channel. Thus, if the photon density is treated on the same footing as the other parton densities, both the qq and γγ channels should be supplemented by virtual and real corrections at NLO EW. Triboson production processes of type pp → V V V and pp → V V H, with V = W, Z, contribute at the same perturbative order as the NLO EW corrections to pp → V V . Thus, in principle, triboson production can be regarded as part of the NLO EW corrections to diboson production. However, diboson and triboson production yield different experimental signatures, and are typically handled as separate processes in experimental analyses. For this reason, and in order to avoid double counting between diboson and triboson production, we do not include pp → V V V /V V H in the EW corrections to pp → V V .

Giant K-factors
At large transverse momenta, as shown in section 3.3, the NLO QCD corrections to vector-boson pair production can become as large as O(1) or even O(10). These so-called giant QCD Kfactors [37,45] arise from phase space regions that are kinematically forbidden at LO and are populated through the emission of hard QCD radiation starting at NLO. Giant K-factors appear in inclusive observables that require a vector boson with very large transverse momentum, p T,V1 ∼ Q M W , while leaving the second vector boson unconstrained. At LO, as a result of momentum conservation, the recoil of the first vector boson is absorbed by the second one, i.e. p T,V2 = p T,V1 . Instead, in the NLO QCD radiative process pp → V V j the recoil can also be absorbed by the additional hard jet, while the second vector boson becomes much softer, p T,V2 = O(M W ). Upon integration over the phase space of the second vector boson, this radiative mechanism results in the cross section [37] where σ V j LO describes the hard subprocess pp → V j, while the double Sudakov logarithm on the right-hand side can be understood as the inclusive probability of radiating a soft vector boson.
General real-emission topologies that lead to giant K-factors are depicted in figure 5. They correspond to a hard pp → V j subprocess at the scale Q M W supplemented by soft vector-boson radiation. The corresponding kinematic regions will be referred to as hard-V j regions, and they are characterised by a hard jet with p T,j ∼ Q and a large gap between the leading and subleading vector boson, p T,V2 p T,V1 . Conversely, standard QCD radiation effects correspond to a hard subprocess pp → V V at the scale Q and QCD radiation at scales well below Q. In this case the two vector bosons are comparably hard, and such phase space regions will be classified as hard-V V regions.
Noteworthy, giant K-factors can also arise at NLO EW, where they appear in γq → V V q realemission processes with a hard γq → V q subprocess and soft vector-boson radiation, as well as in crossing-related qq → V V γ processes with a hard qq → V γ subprocess. At NLO EW, in addition to the topologies of figure 5 with gluons replaced by photons, also extra topologies where the soft vector boson is radiated off external photons arise. Here, the giant K-factor mechanism leads to NLO EW effects of order α w log 2 (Q 2 /M 2 W ), and these are dominated by the γq → V V q channel. The appearance of giant K-factors at NLO raises important questions concerning the convergence of the perturbative expansion and the combination of QCD and EW corrections. In this respect, it is important to note that, contrary to QCD logarithmic effects of soft and collinear origin, the large logarithms in eq. (2.1) do not contribute to all orders in α S . In fact, such logarithms do not arise from soft QCD radiation, but from soft vector-boson radiation in combination with the opening of the hard pp → V (V )j channel at NLO QCD. Since this happens only when moving from LO to NLO QCD, higher-order QCD corrections beyond NLO are free from further giant K-factors. 5 Note also that the availability of NNLO QCD corrections makes it possible to verify the stability of the perturbative expansion beyond NLO and to arrive at reliable QCD predictions for observables that feature giant K-factors. 5 Here, we assume that in diboson production at the scale Q M W at least one vector boson with p T, For what concerns the combination of QCD and EW corrections, the presence of giant K-factors raises more serious issues. In particular, the fact that in the relevant high-p T regions the NLO QCD and NLO EW corrections are both strongly enhanced implies sizeable theoretical uncertainties from large unknown mixed QCD-EW NNLO effects. In principle, depending on the observable and the kinematic region, mixed QCD-EW effects can be approximated through a factorised description of QCD and EW corrections (see section 2.6). However, such a factorisation can be justified only in cases where QCD and EW corrections are both dominated by soft corrections with respect to the same hard subprocess. In the case at hand, this condition is not fulfilled since NLO EW effects are driven by logarithmic Sudakov corrections to hard V V production, whereas giant QCD K-factors are dominated by soft EW boson radiation on top of hard V j production. Actually, the leading source of O(α S α) corrections is given by the NLO EW corrections to the enhanced pp → V V j channel, which cannot be captured through a naive factorised combination of the NLO QCD and NLO EW corrections to pp → V V .
When presenting our results in section 3, the problem of giant K-factors in the inclusive phase space will be illustrated. We will show that giant K-factors can be avoided by means of selection cuts that require a similar hardness of the two vector bosons, e.g. by direct requirements on the hardness of the softer vector boson or by imposing a veto against hard QCD radiation. This will restrict the phase space to hard-V V topologies and suppress hard-V j production. Besides reducing the size of mixed QCD-EW higher-order effects and their respective theoretical uncertainties, selecting hard-V V topologies enhances the sensitivity of experimental measurements that aim at extracting new-physics effects in vector-boson pair processes, such as anomalous triple gauge couplings, from the tails of kinematic distributions. On the other hand, a reliable inclusive description of diboson production is indispensable for background simulations in direct searches at the TeV scale. This can be achieved by merging pp → V V and pp → V V j production including NLO QCD and NLO EW corrections as demonstrated in ref. [85]. The extension of this approach to NNLO QCD+EW is beyond the scope of the present paper.

Combination of QCD and EW corrections
When QCD and EW corrections are both large, also NNLO mixed QCD-EW effects of relative O(α S α) and beyond can become important. In order to gain insights into such higher-order effects, we consider a standard additive combination of NNLO QCD and NLO EW corrections and compare it against factorised combination prescriptions. To this end, we express higher-order effects in terms of relative correction factors with respect to the LO differential cross section, which involves O(α 4 ) contributions from the qq and γγ channels. 6 Higher-order QCD contributions can be cast into the form contribution of the loop-induced gg channel, and all other QCD corrections are embodied in the correction factor δ QCD , which includes the O(α S ) and O(α 2 S ) corrections of the qq, qg/qg, gg and qq/qq channels. 7 Similarly, the NLO EW cross section can be written as (2.5) 6 Note that the γγ channel contributes only to ZZ and W W production. The same holds for the gg channel contributing at NNLO QCD. 7 Here and in the following, higher-order contributions (or terms) of O(α n S α 4+m ) are also referred to as corrections (or effects) of O(α n S α m ).
where all O(α) corrections in the qq, γγ and qγ (includingqγ is implicitly understood) channels are incorporated into the factor δ EW . For the combination of QCD and EW corrections we consider three different prescriptions.
NNLO QCD+EW The first prescription amounts to a purely additive combination, NNLO QCD×EW As a possible approximation of the mixed QCD-EW higher-order corrections we consider the factorised combination where the EW correction factor is applied to the entire NNLO QCD cross section except for the loop-induced gg channel, for which the EW corrections δ EW of the qq and γγ channels are not applicable. The prescription (2.7) can also be written in the form Thus, the factorised combination (2.8) generates extra O(α S α) and O(α 2 S α) mixed QCD-EW corrections. Provided that the dominant sources of QCD and EW corrections factorise, such terms can be regarded as a reasonable approximation of mixed QCD-EW effects. For instance, at scattering energies Q M W this assumption is justified when EW effects are dominated by Sudakov logarithms, and the dominant QCD effects arise at scales well below Q, factorising with respect to the underlying hard-V V process. In such cases, the factorised prescription (2.7) should be regarded as a superior prediction as compared to the additive combination (2.6).
NNLO QCD×EW qq As a motivation for an alternative combination, let us highlight the role of individual partonic channels in the factorised formula (2.7). To this end we rewrite the QCD corrections as where δ qq QCD includes the same QCD corrections as δ QCD , but is normalised to the LO cross section in the qq channel. Moreover we split the EW corrections into contributions from the qq and γ-induced channels, Here in the factor δ qq EW we include only O(α) corrections from the qq channel, whereas all other O(α) effects stemming from the γγ and qγ channels 8 are included in the factor δ γγ/qγ EW . Using the notation of eqs. (2.9)-(2.10) we can rewrite the factorised formula (2.7) as where the EW K-factor corresponds to and can be regarded as the weighted average of the corrections in the qq and γγ channels. The representation (2.11) demonstrates that the factorised combination does not induce any O(α S ) effect in the γγ and gg channels. The only nontrivial factorised correction arises from the term δ qq QCD δ EW , where QCD corrections to the qq channel are combined with the average EW corrections in the qq and γγ channels. The latter includes contributions from qγ channels that can give rise to giant EW K-factors, in which case a factorised treatment is not justified (see section 3.3 for a detailed discussion). For this reason we consider the alternative combination formula where the factorisation of EW corrections is restricted to the qq channel, while photon-induced channels and the loop-induced gg contribution are treated in an additive way. In analogy with eq. (2.8), the prescription (2.13) can be rewritten as 9 (2.14) Both multiplicative combinations (2.8) and (2.14) are implemented at the level of individual distributions by computing the relevant differential EW K-factors δ EW and δ qq EW on a bin-by-bin basis.
When QCD corrections are dominated by hard effects that do not factorise with respect to the hard-V V subprocess, like in the case of giant K-factors, the difference between the additive and the modified multiplicative combination can be regarded as a rough indication of the magnitude of potential effects of O(α S α) and beyond. More details on uncertainty estimates of missing mixed QCD-EW corrections will be discussed in section 3. As far as pure QCD uncertainties are concerned, they are estimated through customary variations of the renormalisation and factorisation scales. Uncertainties from missing EW corrections beyond O(α) are not addressed in this paper: the dominant source of O(α 2 ) effects at high energy are two-loop Sudakov logarithms of the form α 2 w log 4 (Q 2 /M 2 W ), which should be included in order to achieve few-percent accuracy at high p T . The expected size of these two-loop EW effects, assuming naive Sudakov exponentiation, is around

Phenomenological results
In this section we present numerical results for the selected diboson processes All cross sections correspond to the contribution from one lepton family , = e or µ, and = .
In the case of W Z production, the QCD and EW corrections are combined at the level of the individual W + Z and W − Z subprocesses, and their cross sections are summed up afterwards.

Setup
In the following we specify the employed input parameters, scale choices, PDFs, and selection cuts.  GeV. This renders real-emission channels with bottom quarks in the final state separately finite, allowing us to remove such channels from our predictions. In this way, the W W cross section can be defined without any contamination from resonant top-quark contributions that would otherwise enter the higher-order QCD corrections, as discussed in detail in refs. [19,24], and also the higher-order EW corrections, as discussed in ref. [43]. Unstable particles are treated in the complex-mass scheme [52,68], where width effects are absorbed into the complex-valued renormalised squared masses The complex-mass scheme ensures a gauge-invariant description of non-resonant and off-shell effects in the entire phase space. The EW couplings are derived from the gauge-boson masses in the G µ -scheme using where G µ is the Fermi constant, and the complex-valued weak mixing angle is defined through The G µ -scheme guarantees an optimal description of pure SU (2) interactions at the EW scale by resumming large logarithms of the light-fermion masses associated with ∆α(M 2 Z ) and m 2 t /M 2 W enhanced EW corrections related to ∆r, which summarises the weak radiative corrections to the muon decay [86]. This scheme is used for all partonic channels including photon-induced ones. In this respect we note that, at variance with on-shell final-state photons, the appropriate coupling for initial-state photons is not α(0), but α at a scale of the order of the factorisation scale [43,87], as is the case in the G µ -scheme. 10 3)-(2.4) and (2.9). 10 In general, any EW scheme available within OpenLoops [52] is accessible through the Matrix+OpenLoops framework, such that also the α(M Z )-scheme, the α(0)-scheme, or mixed α(0)-short distance schemes (relevant for processes with resolved final-state photons) can be employed. The various schemes can be selected through the model.dat input card of Matrix, where the default (and recommended) choice is the Gµ-scheme.
The CKM matrix is assumed to be diagonal. In fact, the negligible mixing of the first two and the third quark generations justifies a block-diagonal approximation with Cabibbo-type mixing. Since all quarks of the first two generations are treated as massless, the unitarity of the CKM matrix ensures the independence of all physical results of the Cabibbo angles for charge-neutral processes. 11 Hence, for ZZ and W W production the choice of a diagonal CKM matrix is fully justified. The effects from a non-trivial CKM matrix for W Z production are also quite small, being below −1% in case of the NLO QCD cross section and negligible for the NLO QCD K-factor [20,25].
Renormalisation and factorisation scales In order to avoid large shape corrections in the tails of kinematic distributions, we use dynamic QCD scales. Specifically, the renormalisation scale µ R and the factorisation scale µ F are set to Here, the transverse energies of the leading and subleading vector bosons, V 1 , V 2 , are defined as where p T,Vi and m Vi are the transverse momentum and the invariant mass of the off-shell + − /νν pair for a Z-boson (or photon in case of + − ) and a + ν / −ν pair for a reconstructed W ±boson. 12 At NLO EW, in order to guarantee infrared safety, the above dynamic scale is determined using the momenta of dressed leptons as defined below in eq. (3.9). Our central scale corresponds to ξ R = ξ F = 1, and theoretical uncertainties due to missing higher-order QCD contributions are estimated through customary 7-point variations The multiplicative combinations of QCD and EW corrections involve EW K-factors of the form δ EW = dσ NLO EW (µ R , µ F )/dσ LO (µ R , µ F ), where scale variations are correlated between the numerator and the denominator. However, the µ R dependence is completely absent, while the µ F dependence cancels almost exactly in the ratio. As a consequence, the multiplicative combinations of QCD and EW corrections results in essentially the same QCD scale uncertainties as for the underlying QCD cross sections.
Parton densities As parton distribution functions (PDFs) for the calculation of hadron-level cross sections we employ, consistent with the chosen n f = 5 or n f = 4 flavour schemes, the corresponding NNPDF31_nnlo_as_0118_luxqed set [88], which is based on the LUXqed methodology [89] for the determination of the photon flux. 13 Selection cuts For the infrared-safe definition of observables at NLO EW, collinear chargedlepton-photon pairs within a cone of radius R rec = 0.1 are systematically recombined into dressed leptons according to 14 This statement is exact for ZZ production at all considered orders, while for W W production a tiny dependence is introduced at NNLO QCD through double-real contributions in the qq/qq/qq channels through subprocesses of the type of ud → W + W − cs. 12 For the different-flavour processes considered in this paper this scale choice is unambiguously defined. For the same-flavour processes, also available in Matrix, corresponding scale choices are implemented via reconstruction of the intermediate vector bosons according to a distance criterion.
13 See [? ] for an alternative treatment of photon-induced processes via direct application of the structure function formalism avoiding parton densities for initial-state photons.
14 In case this condition holds for more than one lepton, the lepton closest to the photon in ∆R γ is dressed.
All cuts, observables and scales are defined in terms of dressed leptons. We apply uniform fiducial cuts for ZZ, W W , and W Z production that mimic the standard selections imposed in experimental measurements. Specifically, we require p T, ± > 25 GeV , |η ± | < 2.5 (3.10) for the transverse momentum and pseudorapidity of each dressed charged lepton. For processes with neutrinos in the final state we additionally require p T,miss > 25 GeV (3.11) for the missing transverse momentum calculated as the vectorial sum of the neutrino momenta. Moreover, the invariant mass of same-flavour + − pairs is restricted to the Z-mass window Reconstructed vector bosons In the following we present differential distributions in the transverse momenta and invariant masses of the vector bosons. Such observables are defined in terms of the reconstructed vector-boson momenta where all charged leptons are potentially dressed. Here, the vector bosons are kept off-shell, and the correctness of the reconstruction is guaranteed by pairing charged leptons and neutrinos of the same generation, selecting the appropriate neutrino and/or or anti-neutrino momenta at truth level. The reconstructed vector bosons are ordered according to their p T , and the leading and subleading boson are labelled V 1 and V 2 , respectively.
Jet veto As discussed in section 2.5, in order to avoid giant K-factors at high p T , we impose selection cuts that single out regions dominated by hard-V V production while suppressing regions dominated by hard-V j production. To this end we apply a jet veto. More precisely, we impose a restriction on the total jet transverse energy, where we include all reconstructed anti-k T jets [90] with R = 0.4, |y| < 4.5, and arbitrary p T . The upper bound for H jet T is defined in terms of the hardness of the diboson system. Specifically, we use the total leptonic transverse energy, The effect of the above jet veto on the relative hardness of the two vector bosons at high p T can be quantified by translating eq. (3.16) into a lower bound for the p T of the softer vector boson. This is easily achieved by combining eq. (3.16) with which leads to This inequality can be further refined by relating H lep T to the transverse momenta of the two vector bosons. To this end we can use where H T,Vi denotes the total transverse energy of the decay products of the V i vector boson. In the following we assume that both vector bosons are nearly on-shell. Moreover, we focus on the region 20) where the products of the decay of the leading boson, V 1 → ab, are nearly collinear. Thus For the decay of the softer boson, V 2 → cd, we can use where the inequality holds for energies E i in any reference frame that is connected to the laboratory frame via a longitudinal boost, while the last identity is based (without loss of generality) on the reference frame where the longitudinal component of p V2 vanishes. In this way we arrive at which confirms that the two bosons are similarly hard. As demonstrated in section 3.4, this bound is violated only by highly suppressed off-shell contributions. Moreover, at very high transverse momenta, the ratio between the p T of the softer and harder vector bosons is typically well above 2/3 and exceeds 0.9 on average.

Fiducial cross sections
Predictions and scale variations for the fiducial cross sections of the diboson processes ( Table 3. The behaviour of QCD and EW corrections in table 3 is consistent with the well-known results in the literature. The NLO EW corrections amount to about −6% for ZZ production, and only −2% and −3% for W W and W Z production, respectively. The NLO QCD corrections range from +36% for ZZ production up to +73% for W Z production. In the latter case, the huge NLO effect is due to the presence of an (approximate) radiation zero at LO [77]. The NNLO QCD corrections are again positive and vary between +11% and +16%. The largest NNLO effects are found in the case of neutral final states, where the contributions from loop-induced gg channels are sizeable.
As discussed in the following, the sizeable impact of QCD corrections has non-negligible implications on their combination with EW corrections. Comparing NNLO QCD and combined predictions, we observe that the multiplicative prescriptions yield, as expected, relative correction factors that are very close to the NLO EW K-factor. Restricting the factorisation to the qq channel tends to increase the net EW effect by at most one percent in the W Z case. In contrast, the additive prescription leads to a significant reduction of the EW effects by up to 2%. This is due to the fact that in the additive combination (2.6) the sizeable higher-order QCD correction dσ LO δ QCD remains free from EW corrections. In this respect, it is reasonable to assume that dσ LO δ QCD should receive a comparable EW K-factor as the LO cross section. For instance, any short-distance higher-order EW effect associated with the SU(2)×U(1) couplings is expected to factorise with respect to the QCD corrections. Thus, using a factorised combination appears to be better motivated than just omitting any mixed QCD-EW effect. Nonetheless, in the absence of more rigorous theoretical arguments, the difference between additive and multiplicative combinations may be taken as a conservative uncertainty estimate at the level of inclusive cross sections. Matrix + OpenLoops baseline cuts

Differential distributions without veto against QCD radiation
In figures 6-10 we present differential distributions that characterise the behaviour of the vector bosons and their decay products in ZZ, W W and W Z production. The inclusive cuts (3.10)-(3.12) are applied throughout without jet veto, and the vector bosons are reconstructed as detailed in eq. (3.13). In order to highlight EW Sudakov effects, transverse-momentum and invariant-mass distributions are plotted up to the multi-TeV range, where event rates are suppressed by several orders of magnitude. In this respect, we note that the relevance of our results goes beyond leptonic diboson signatures at the LHC, where the energy reach at high luminosity is around 1 TeV in transverse momenta and 2 TeV in invariant masses. In fact, in the case of hadronic vector-boson decays, where statistics is significantly higher, the QCD and EW corrections are expected to behave in a rather similar way as for the leptonic channels discussed in this paper. Indeed, the dynamics of the dominant giant QCD K-factors and large EW Sudakov corrections only depends on the diboson production processes pp → V V and pp → V V j, i.e. it is independent of the vector-boson decays. Moreover, since EW Sudakov logarithms are mainly sensitive to the final-state transverse momenta and depend only weakly on the collider energy, the presented results can be regarded as   an indication of the behaviour of EW corrections in the multi-TeV regime at future pp colliders. The distribution in the p T of the softer vector boson is shown in figure 6. This observable is an optimal probe of hard-V V production since requiring large p T,V2 automatically ensures that also p T,V1 is large, thereby excluding giant K-factor effects due to hard-V j production. As far as the QCD corrections are concerned, at low p T,V2 we find similar NLO and NNLO QCD effects as in the case of the fiducial cross section, i.e. the NLO corrections are rather large for W Z production and significantly less pronounced for W W and ZZ production, while NNLO effects are at the 10-20% level for all three processes. At large p T higher-order QCD effects tend to decrease. In particular, the NNLO QCD corrections behave in a very similar way for all three processes and go down to the few-percent level in the TeV range. Also scale uncertainties decrease at large p T , being only a few percent in the high-p T,V2 tail. The EW corrections are negative, and in the tails they grow like log 2 (p 2 T,V2 /M 2 W ), as expected for EW Sudakov logarithms in hard V V production. At 1-2 TeV the (negative) EW corrections to ZZ, W W and W Z production range, respectively, between 60-75%, 50-65% and 35-45%. Electroweak Sudakov effects remain very large also in the combination of QCD and EW corrections, and due to the moderate size of QCD corrections in the tails, their impact depends only weakly on the employed combination prescription. The two multiplicative combinations are very close to each other, and they are both well motivated since the observable at hand is dominated by hard-V V topologies. The difference between additive and multiplicative combinations can be regarded as upper bound for the uncertainty due to missing O(α S α) effects. In the p T range where NLO EW corrections approach the 50% level, the main source of theoretical uncertainty is given by missing EW Sudakov logarithms at relative O(α 2 ). Based on naive Sudakov exponentiation, their size is expected to be around 1 2 δ 2 EW . In figure 7 we present the distribution in the diboson invariant mass, m V V . Again, at small Matrix + OpenLoops baseline cuts   invariant mass we observe similar QCD corrections as for the fiducial cross sections. At very large m V V the NLO QCD corrections tend to decrease. As for the NNLO corrections, in W Z production they increase mildly towards large m V V and reach up to 20% in the TeV regime. In W W production the additional gg channel, whose contribution decreases at large m V V , results in a relatively flat total NNLO correction at around 10%, with residual scale uncertainties at the 5% level. On the other hand, for ZZ production we observe a significant enhancement of the NNLO corrections up to 40% at very high invariant masses in the multi-TeV range. This is due to the effect of acceptance cuts and their interplay with the opening of the qq → V V qq ,qq → V Vqq and qq → V V qq channels at NNLO. The reason is that, at large m V V , in these channels the gauge bosons are mainly emitted in the forward and backward directions. Hence, the rapidity of their leptonic decay products is often too large to pass the acceptance cuts on charged leptons in W W and W Z production. In contrast, in the case of ZZ production the absence of rapidity cuts on the Z → νν decay products leads to significantly larger contributions from these channels. In fact we have checked that the NNLO QCD corrections to ZZ, W W and W Z production behave in a very similar way when applying the same (technical) cuts on neutrinos and charged leptons, and switching off loop-induced gg contributions. The EW corrections to the m V V distribution are negative, and in the tails they grow like double Sudakov logarithms. However, their impact is less pronounced than for the p T,V2 distribution in figure 6. This is due to the fact that diboson production at large m V V is dominated by tand u-channel topologies where the gauge bosons are mainly emitted in the forward/backward regions, and the scales t, u that enter Sudakov logarithms are well below m V V . The largest EW corrections are found in the ZZ channel, where they amount to −15% at 1 TeV. In the combination of QCD and EW corrections the difference between additive and multiplicative prescriptions is similarly large as NNLO QCD scale uncertainties, and depending on the process it can reach up to 10-20% in the multi-TeV region. For W W and W Z production we also find a difference of up to 5% between the two factorised prescriptions. This effect can be attributed to photon-induced γq → W V q channels, where the topologies with t-channel W bosons that couple to the initial-state photons (see e.g. figure 4 l) yield a significant (positive) NLO EW contribution.
The distribution in the transverse momentum of the harder vector boson, presented in figure 8, shows a completely different behaviour of the higher-order effects. 15 At 100 GeV, the NLO QCD corrections are as large as a factor two, and their size grows with p T,V1 reaching five to twenty times the LO cross section at 2 TeV. These giant NLO K-factors are driven by hard-V j subprocesses, where the recoil of the harder vector boson is absorbed by the hard jet, while the second boson becomes soft and generates large logarithmic corrections of the form of eq. (2.1). As discussed in section 2.5, such large logarithms occur only at NLO and do not need to be resummed to all orders. With other words, the perturbative expansion beyond NLO is expected to be well behaved. This is confirmed by the NNLO QCD results in figure 8, where we observe NNLO QCD corrections at the level of 30-50%, which are moderate compared to the giant NLO QCD corrections. At the same time we find a reduction of scale uncertainties down to the 10% level when moving from NLO to NNLO QCD. In practice, at large p T,V1 , the pp → V V cross section at NNLO QCD is largely dominated by pp → V V j at NLO QCD, which provides the first reliable QCD prediction for the tail of the distribution at hand.
The behaviour of the EW corrections in figure 8 depends strongly on the diboson final state. For ZZ production we observe large negative corrections that grow with p T , reaching about −50% with respect to LO at 1 TeV. This indicates that the NLO EW corrections are dominated by Sudakov logarithms, which originate from the virtual corrections to hard-ZZ production. 16 However, the tail of the distribution at hand is dominated by higher-order QCD contributions corresponding to hard-Zj production. As a result, when NLO EW effects are combined with NNLO QCD predictions using the additive prescription, their relative weight goes down by up to a factor of seven. This is due to the fact that the dominant pp → V V j process does not receive any EW corrections in the additive combination. In contrast, in the multiplicative combinations the δ EW correction factor is applied also to pp → V V j, resulting in large EW Sudakov effects.
In the case of W W and W Z production, the behaviour of the EW corrections at large p T,V1 is completely different: instead of negative Sudakov effects, we observe positive NLO EW K-factors that grow with p T,V1 and approach a factor two in the case of W Z production. This is not due to the absence of EW Sudakov logarithms, but to the fact that they are overcompensated by a positive giant EW K-factor stemming from the emission of one additional (anti)quark at NLO EW. More precisely, this large positive contributions originate from photon-induced γq → W V q and γq → W Vq topologies with t-channel W -boson exchange (see figure 4 l) that we already discussed in the context of the m V V distribution. Such topologies contribute only to W W and W Z production and, as a result of the W W γ triple gauge coupling, they are strongly enhanced as compared to γq → V V q topologies that involve t-channel (anti)quark exchange (see e.g. figure 4 k) and contribute also to ZZ production [37]. Due to the interplay between large positive QCD effects, large negative Sudakov EW effects and large positive photon-induced EW effects, the combination of QCD and EW corrections is strongly dependent on the employed prescription. In the additive combination the relative impact of EW effects is minimal, due to both the large cancellation between positive and negative EW effects and the absence of EW corrections to the dominant hard-V j subprocess. This results in a dramatic underestimate of higher-order EW contributions. In contrast, in the factorised combinations the EW K-factor is applied to the full NNLO QCD cross section, resulting in sizeable Matrix + OpenLoops baseline cuts  net EW effects in the tails. In the fully factorised prescription the overall effect is positive due to the dominance of the γ-induced bremsstrahlung contributions to the EW corrections. However, this prescription is not well justified in the case at hand, since the large γ-induced EW effect is due to the opening of the V V q three-body phase space, while in the process pp → V V j, which dominates the (N)NLO QCD corrections, the three-body phase space is already open. This justifies the restriction of the factorisation to the qq channel, which leads to an overall large negative EW correction of Sudakov-type in the tail. Instead, the fully factorised NNLO QCD×EW combination should be discarded.
In general, in the p T,V1 distribution, the modified multiplicative prescription (2.13) yields a Sudakov-type behaviour for all processes. However, in this approach the NLO EW K-factor for hard-V V production is effectively transferred to phase space regions that are dominated by hard-V j production. This can largely overestimate the EW corrections since, based on the general properties of Sudakov logarithms, the EW corrections to hard-V j production are expected to be roughly half as large as for hard-V V production. For this reason, it is reasonable to use as nominal prediction the average of the additive combination (2.6) and the modified multiplicative combination (2.13), while their spread can be interpreted as O(α S α) uncertainty band.
The distribution in the p T of the hardest charged lepton, shown in figure 9, features a qualitatively similar behaviour of the QCD and EW corrections as for p T,V1 , shown in figure 8. This demonstrates that the issues related to giant K-factors are independent of the reconstruction of the vector bosons. This is substantiated by the distribution in the missing transverse momentum, p T,miss , in figure 10. In the ZZ and W Z channels, where p T,miss is identical or strongly correlated with the p T of one of the vector bosons, the QCD and EW corrections feature qualitatively similar giant K-factor effects as for the distribution in p T,V1 . Also in the W W channel we observe such a Matrix + OpenLoops baseline cuts  behaviour. However, in this case the giant K-factor effects are already visible at 100 GeV, where the NLO QCD corrections grow very quickly, reaching O(10) at few hundred GeV. As discussed in detail in ref. [43], this is due to the fact that for on-shell W W production at LO the vanishing transverse momentum of the diboson system implies that p T,miss ≤ M W . Hence, at LO the production of a neutrino pair with p T,νν = p T,miss M W emerges only from strongly suppressed configurations with off-shell vector bosons. In contrast, thanks to the finite p T,W W generated by the radiative process pp → W W j, at NLO QCD the region p T,miss > M W can be populated with on-shell gauge bosons, resulting in a much larger cross section with respect to LO. Analogous giant K-factor effects show up also at NLO EW through the γq → W W q channel, which leads to a positive correction that overcompensates EW Sudakov logarithms and leads to an enhancement of a few tens of percent in the tail. Since such QCD and EW giant K-factors are due to the opening of the p T,W W > 0 phase space, higher-order effects beyond NLO are expected to be free from this kind of enhancements. Indeed this is confirmed by the moderate size of the NNLO QCD corrections over the entire p T,miss range. For what concerns the combination of QCD and EW corrections, similar considerations as for figure 8 apply.
In summary, for observables that are dominated by hard-V V production, such as the distributions in p T,V2 and m V V , the factorised combinations of EW and QCD corrections are well motivated, and O(α S α) uncertainties are expected to be significantly smaller than the differences between the factorised and additive prescriptions. Conversely, in the tails of the distributions that are dominated by hard-V j production, 17 such as for the distributions in p T,V1 , p T, 1 and p T,miss , due to the presence of giant K-factors all combination prescriptions considered in this paper are inadequate. 17 We remind the reader that, as introduced in section 2.5, by hard-V j production we refer to pp → V V j with the second vector boson being soft. Matrix + OpenLoops Matrix + OpenLoops In particular, the fully factorised prediction (2.7) should be discarded, and the band defined by the difference between the additive combination and the modified multiplicative combination can be considered as uncertainty related to O(α S α) corrections and higher, and its average as nominal prediction. In the following we show how giant K-factors and related issues can be avoided by means of appropriate selection cuts, while a more accurate QCD-EW combination in the presence of giant K-factors is deferred to future work.

Results with a veto against hard QCD radiation
As anticipated in section 2.5, giant K-factors can be avoided by restricting the phase space of pp → V V +jets radiative processes to regions dominated by hard-V V production while rejecting regions dominated by hard-V j production in association with soft vector-boson radiation. This can be achieved either via direct selection cuts that require similar hardness for the two vector bosons or through a veto on QCD radiation. In the following we adopt the latter approach by imposing an upper bound on the total jet transverse energy H jet T , as defined in eq. (3.16). The interplay of hard-V V and hard-V j production and the impact of the jet veto are illustrated in figure 11, which shows NLO QCD predictions for the ratio of the p T of the softer over the one of the harder vector boson, r 21 = p T,V2 /p T,V1 . Due to momentum conservation, at LO the two vector bosons are back-to-back, i.e. r 21 = 1, and the phase space region with r 21 < 1 is populated only through QCD (or QED) radiation starting at NLO. If only baseline cuts are applied (light red curve in the upper frames), the cross section is strongly peaked in the r 21 → 1 region, which corresponds to V V j configurations with a hard-V V system accompanied by soft and collinear QCD radiation. The bulk of the cross section originates from r 21 > 0.8 and is largely unaffected by the jet veto (dark red curve in the upper frames). In the hard region p T,V1 > 1 TeV (see lower frames) the shape of the distribution in r 21 is very different. In the absence of a jet veto, the cross section develops a very strong enhancement in the region of small r 21 , which corresponds to hard-V j production in association with a soft vector boson, and is responsible for the giant K-factors discussed in section 3.3. Note that this hard-V j enhancement exceeds the large-r 21 enhancement due to soft QCD radiation in a very broad region, from very small r 21 up to r 21 0.8. For p T,V1 > 1 TeV, we observe that, consistently with the bound (3.24), the jet veto acts as a rather sharp cut-off rejecting hard-V j configurations with r 21 < 2/3. In fact, the region below 2/3 is filled only by strongly suppressed off-shell contributions. As a result of the jet veto, the cross section is dominated by the region of large r 21 . Actually, the average value of r 21 is well above 2/3 and exceeds 0.9, confirming that hard-V V production dominates. At large r 21 the distribution becomes logarithmically divergent, and the singularity is cancelled by the virtual contribution at r 21 = 1. Any tight cut in the vicinity of r 21 = 1 would induce large logarithmic corrections of soft and collinear origin, requiring an all-order resummation. However, the effective r 21 -cut imposed through the jet veto is far enough from the logarithmically enhanced region of the r 21 -distribution such that fixed-order predictions are well behaved.
In figures 12-15 we study the effect of the jet veto on the distributions discussed in section 3.3. The veto has no effect on LO predictions, since it acts only on (N)NLO radiation. For observables that are free from giant K-factors, the qualitative behaviour of QCD and EW corrections remains more or less unchanged, apart from a significant reduction of the (N)NLO cross section. This is illustrated in figure 12 for the distribution in m V V . At the lower end of the plotted m V V range the veto reduces the NNLO QCD cross section by about 50% for all diboson processes, while its impact Matrix + OpenLoops   is less pronounced at large m V V and amounts to a reduction of 20-35% in the multi-TeV region. 18 This demonstrates that the entire distribution at hand is indeed dominated by hard-V V production also in the inclusive case. For this reason, the effect of the jet veto on EW corrections is quite small: in the case of ZZ production it is negligible, while for W W and W Z production it enhances the negative correction in the tail by about 5% due to the suppression of the positive γq → W V q contribution (see the discussion of figure 7). Since NNLO QCD corrections are largest (up to −25%) for m V V 1 TeV, where EW effects are smaller, while NNLO QCD corrections decrease in the tails, their combination turns out to depend only weakly on the employed prescription.
For the distribution in p T,V2 (not shown) the impact of the jet veto turns out to be even milder than in the case of m V V . In contrast, in the p T,V1 distribution, 19 shown in figure 13, the jet veto leads to a drastic suppression of the giant K-factor effects observed in figure 8. At 100 GeV the NLO QCD corrections are reduced to 20-30%, and at large p T they grow much slower than in the inclusive case, reaching only 50% in ZZ and W W production and 100% in W Z production. This corresponds to a reduction by a factor of ten to twenty as compared to the inclusive case. The moderate size of NNLO QCD corrections (10-20%) and the smallness of NNLO scale uncertainties indicate that, in spite of the huge suppression of the hard-V j contributions induced by the jet veto, the perturbative expansion of the hard-V V process remains free from large logarithmic QCD 18 These quantitative statements can directly be inferred from the ratio plots in figure 7 and 12. For example, at low m V V in the case of ZZ production without jet veto (figure 7) the NLO corrections amount to about +30% while the NNLO corrections read +15%, thereby leading to a positive 50% effect. When the jet veto is applied ( figure 12) the NLO and NNLO effects amount to −5% and −15%, respectively. Since the jet veto does not affect the LO result, this amounts to the quoted 50% reduction at low m V V . 19 The quantitative effects observed in the tail of this distribution are documented in more detail in appendix A.
contributions beyond NLO. In the presence of the jet veto also the huge positive EW corrections (see figure 8) are strongly suppressed. Hence, all three diboson processes feature large negative corrections of Sudakov type that are qualitatively and quantitatively similar to those observed in the tail of the p T,V2 distribution without jet veto. This is a further demonstration that, in the presence of the jet veto, the distribution at hand is dominated by hard-V V production with p T,V2 ∼ p T,V1 . As a result of the jet veto, all combinations of QCD and EW corrections yield large negative EW corrections independently of the employed prescription. However, the sizeable NLO QCD effects in the tails lead to non-negligible differences between the additive and multiplicative combinations (from a few percent at 500 GeV to 15-30% at 2 TeV). Nevertheless, given that the full NNLO QCD prediction is dominated by hard-V V contributions, the multiplicative combinations can be regarded as realistic approximations of mixed QCD-EW corrections, and their differences with respect to the additive combination is expected to largely overestimate O(α S α) uncertainties. The jet veto affects the QCD and EW corrections to the distribution in the p T of the hardest charged lepton, shown in figure 14, in a qualitatively very similar way as the distribution in p T,V1 . This holds also for the distribution in the missing transverse momentum, shown in figure 15, where the giant QCD K-factor due to the off-shell suppression of the LO W W cross section (see figure 10) is reduced by up to a factor of twenty. Moreover, the positive EW contributions from the γq → W V q channels disappear, and the EW corrections to all three diboson processes feature the usual Sudakov behaviour as expected for hard-V V production.
In general, the above observations support the application of selection cuts that favour hard-V V topologies while suppressing hard-V j topologies and related giant K-factors at high p T . This can be achieved either via an appropriate jet veto, as illustrated here, or via a hardness requirement on the subleading vector boson. Besides reducing theoretical uncertainties from mixed QCD-EW higher-order effects, the suppression of hard-V j contributions has the advantage of increasing the sensitivity of experimental measurements to genuine hard-V V topologies at high p T , and thus to BSM effects related to anomalous gauge couplings.

Summary and conclusions
In this paper we have presented the first combination of NNLO QCD and NLO EW corrections for W + W − , W ± Z and ZZ production processes in their leptonic decay modes at the LHC. Our calculation describes all the final states with two, three and four charged leptons. Spin correlations, interferences and off-shell effects are included throughout.
We have presented a detailed study of higher-order QCD and EW effects for three representative leptonic channels, one for each of the W + W − , W ± Z and ZZ production modes. The inclusion of QCD and EW corrections is crucial in order to reach percent-level precision for important kinematical distributions, as demanded by the experimental analyses. The combination of QCD and EW corrections can be achieved in different ways. We have compared an additive prescription (NNLO QCD+EW) and two multiplicative prescriptions where the factorisation of the EW corrections is complete (NNLO QCD×EW) or restricted to the qq channel (NNLO QCD×EW qq ). The differences between those prescriptions can be exploited to gain insights into higher-order QCD-EW mixed effects of O(α S α) and related theoretical uncertainties.
The high-energy tails of various distributions turn out to be dominated by the production of a hard vector boson in association with a hard jet and a soft vector boson. In these so-called hard-V j regions, QCD and EW corrections behave in a completely different way as compared to regions dominated by hard diboson production (hard-V V regions). Thus theoretical uncertainties need to be assessed with different prescriptions for these two kinds of regions.
The regions dominated by hard-V j production are characterised by the presence of giant QCD K-factors together with large negative EW Sudakov logarithms and possibly positive γ-induced     corrections. The interplay of these large O(α S ) and O(α) effects results in sizeable QCD-EW uncertainties of O(α S α) that cannot be controlled through a naive factorisation of the QCD and EW corrections. In particular, the EW Sudakov corrections to qq → V V production cannot be factorised since they are not applicable to the hard-V j topologies that dominate the (N)NLO QCD cross section. For these reasons, for observables dominated by hard-V j production, the band between the NNLO QCD+EW and NNLO QCD×EW qq combinations can be regarded as O(α S α) uncertainty, and we suggest to use their average as nominal prediction. Further reduction of these uncertainties requires the inclusion of EW corrections to pp → V V j production. The merging with pp → V V retaining NNLO QCD+EW accuracy is left to future work. Giant K-factors and related O(α S α) uncertainties are encountered only in rather inclusive observables at high p T . Moreover, their appearance can be avoided by means of acceptance cuts that select hard-V V regions while suppressing hard-V j production. This can be achieved through direct requirements on the hardness of the two vector bosons or, as studied in this paper, through a relatively mild jet veto. Such phase space restrictions can also enhance the sensitivity of searches for anomalous gauge couplings at large energies.
In general, irrespectively of the presence of the jet veto, in regions dominated by hard-V V production the perturbative QCD expansion is well behaved, and the large EW Sudakov logarithms in the tails of distributions are expected to factorise. In those cases, the factorised prescriptions can be regarded as superior with respect to the additive one, and the NNLO QCD×EW qq combination can be used as nominal prediction. The difference between the two factorised prescriptions can be regarded as uncertainty and is typically rather small, while the difference between additive and multiplicative combinations can largely overestimate O(α S α) uncertainties. A thorough estimate of O(α S α) uncertainties for the hard-V V regions requires more detailed studies along the lines of ref. [91] and is left to future work.
The NNLO QCD and NLO EW calculations presented in this paper have been carried out within the Matrix+OpenLoops framework, and their implementation paves the way to achieve the same level of accuracy for the hadronic production of an arbitrary colourless final state for which the required two-loop QCD amplitudes are available. Theory predictions at NNLO QCD+NLO EW accuracy for the full set of massive diboson processes in the channels with two, three and four finalstate charged leptons will be supported in the next release of Matrix. Such predictions can be instrumental in boosting the sensitivity of a wide range of precision measurements and new-physics searches at the LHC and future pp colliders.