Transverse-momentum resummation for vector-boson pair production at NNLL+NNLO

We consider the transverse-momentum ($p_T$) distribution of $ZZ$ and $W^+W^-$ boson pairs produced in hadron collisions. At small $p_T$, the logarithmically enhanced contributions due to multiple soft-gluon emission are resummed to all orders in QCD perturbation theory. At intermediate and large values of $p_T$, we consistently combine resummation with the known fixed-order results. We exploit the most advanced perturbative information that is available at present: next-to-next-to-leading logarithmic resummation combined with the next-to-next-to-leading fixed-order calculation. After integration over $p_T$, we recover the known next-to-next-to-leading order result for the inclusive cross section. We present numerical results at the LHC, together with an estimate of the corresponding uncertainties. We also study the rapidity dependence of the $p_T$ spectrum and we consider $p_T$ efficiencies at different orders of resummed and fixed-order perturbation theory.


Introduction
Run 1 of the Large Hadron Collider (LHC) has been a great success for the Standard Model (SM). The collected data are in good agreement with the theoretical predictions so far and led to the discovery [1,2] of a resonance at a mass of 125 GeV, which appears to be fully consistent with the SM Higgs boson. Among the most important reactions at hadron colliders is the production of vector-boson pairs. This class of processes gives access to the vector-boson trilinear couplings which may be modified in a large set of Beyond the Standard Model (BSM) theories. Even small deviations in both the production rate and the shape of distributions could be a signal of new physics. Anomalous couplings related to vector-boson pair production have been constrained first by LEP2, and later by the Tevatron for larger invariant masses. ATLAS and CMS will continue to tighten the bounds on anomalous couplings, especially with increasing sensitivity during Run 2 of the LHC. † On the other hand, vector-boson pair production constitutes an irreducible background to new-physics searches as well as Higgs studies. Particularly important are the off-shell effects below the ZZ and W + W − thresholds, relevant for the Higgs signal region, and the high-mass tail used to extract the width of the Higgs boson [4][5][6]. Furthermore, Higgs boson measurements, in particular in the W + W − channel, strongly rely on the background rejection through specific categories based on the transverse momenta of final-state particles, such as the classification into jet bins or in the Higgs transverse momentum and related variables. An accurate modelling of the respective observables for both signal and backgrounds is crucial for such analyses.
The first precise predictions for ZZ production at hadron colliders in the SM were obtained at the next-to-leading order (NLO) already more than 20 years ago for stable Z bosons [7,8], and the leptonic decays were added in Ref. [9]. The full spin correlations and off-shell effects at the NLO were first included in Refs. [10,11] using the corresponding one-loop helicity amplitudes [12]. An important loop-induced contribution proceeds through gluon fusion; it is enhanced by the gluon densities and was first computed for on-shell Z bosons in Refs. [13,14], while leptonic decays were included later [15][16][17]. All the contributions to ZZ production discussed so far are implemented in the numerical program MCFM [18]. Electroweak (EW) corrections were evaluated in Ref. [19,20], while on-shell ZZ+jet production is known through NLO QCD [21,22]. Recently, the inclusive cross section for the production of ZZ at the next-to-next-to-leading order (NNLO) has been presented [23].
The production of W + W − pairs constitutes the largest cross section among all massive vectorboson pair production modes. On the other hand, its leptonic decay (W + W − → l + l − νν) embodies the most challenging experimental signature since the presence of two neutrinos prohibits a reconstruction of mass peaks. Therefore, a precise understanding of both the signal and the background is required.
The W + W − cross sections for on-shell W bosons at the NLO [24,25] as well as the gluon-fusion component [13,26] have been known for decades. Also in this case, spin and off-shell effects were included in the NLO prediction [10,11] after the relevant one-loop helicity amplitudes had been computed [12]. Leptonic decays for the loop-induced gluon-fusion contribution were considered in Refs. [27,28]. More recently, also interference effects with the Higgs production mode through gluon fusion were determined [29]. Analogously to ZZ production, all the contributions to W + W − production discussed so far are implemented in MCFM [18]. Furthermore, EW corrections have been evaluated [20,30,31]. On-shell W + W − production in association with one jet has been studied through NLO QCD in Refs. [32][33][34]. Detailed Monte Carlo simulations of e + ν e µ −ν µ production in association with up to one jet at NLO have been presented in Ref. [35]. Recently, the first NNLO results for the inclusive W + W − cross section have been obtained [36].
Due to the very recent computation of the qq → V V helicity amplitudes at two-loop order [37,38], the inclusion of the off-shell effects and the leptonic decays in the NNLO cross section is expected in the near future. In the meanwhile, also the calculation of gg → V V helicity amplitudes has been performed at two-loop order [39,40]. This renders the evaluation of NLO QCD corrections to the gluon-fusion channel feasible.
The production of ZZ pairs in hadron collisions has been measured extensively at the Tevatron and the LHC (see Refs. [41][42][43][44][45][46][47] for some recent results). The W + W − cross section has also been measured already at the Tevatron, see e.g. Ref. [48], and at the LHC both at 7 TeV [49, 50] and 8 TeV [46, 51, 52]. The ATLAS collaboration recently reported an excess [51] with respect to the SM prediction, which has drawn a lot of attention on the W + W − process, since the W + W − final state is a typical signature in many BSM scenarios [53]. In the meanwhile, the excess has been alleviated to a significant degree by the recent NNLO computation [36]. The more recent measurement by the CMS collaboration [52] is in good agreement with the NNLO prediction.
The transverse-momentum distributions of ZZ and W + W − pairs are among the most important differential observables for these processes. The p T spectrum has already been measured in the case of ZZ production [47] at the LHC. Transverse-momentum resummation for ZZ and W + W − production has been studied in Refs. [54][55][56][57][58]. In all these calculations the resummed computation is essentially performed up to next-to-leading logarithmic (NLL) accuracy (next-to-next-to-leading logarithmic (NNLL) effects in the Sudakov exponent are considered in Refs. [55,57,58]) and matched to the fixed-order result up to the first order in the QCD coupling α S .
In this paper we consider transverse-momentum resummation for the production of both ZZ and W + W − pairs. We use the formalism of Ref. [59] to perform the first computation of the p T spectrum at full NNLL accuracy matched to the O(α 2 S ) fixed-order result valid at large p T . Although our focus is on the inclusive p T spectrum of the ZZ and W + W − system, our computation is fully differential in the degrees of freedom of the vector bosons and allows us to include their EW decays, once the helicity amplitudes are implemented ‡ . The p T -resummation formalism of Ref. [59] is closely related to the subtraction method of Ref. [63], which was used to compute the NNLO cross section for these processes [23,36]. For W + W − production we employ the four-flavour scheme (4FS) and remove all contributions with final-state bottom quarks from our computation of the W + W − transverse-momentum distribution in order to eliminate the contamination from tt and W t production. The difference to the prediction in the five-flavour scheme (5FS), where such terms have to be consistently subtracted, has been shown to be small for the NNLO inclusive cross section [36]. Furthermore, we neglect the loop-induced gluon-fusion contribution throughout this paper, since, up to NNLO, it contributes only at p T = 0.
We note that in the CMS measurement reported in Ref.
[58] of the p T spectrum has been used to correct the spectrum from the Monte Carlo simulation. Along these lines, the computation reported in this paper will be useful in order to validate the predictions obtained from Monte Carlo simulations for both ZZ and W + W − production, as done in the case of Higgs boson production with the calculation of Ref. [59].
The manuscript is organized as follows. In Section 2 we review the transverse-momentum resummation formalism applied to vector-boson pair production. In Section 3 we report our numerical results, starting with remarks on the choice of the resummation scale in Section 3.1. In Section 3.2 we present our numerical predictions for the inclusive p T spectrum and study the ensuing uncertainties. In Section 3.3 we analyse the behaviour of the spectrum at different rapidities of the vector-boson pair. In Section 3.4 we investigate p T efficiencies at different orders in resummed and fixed-order perturbation theory. In Section 4 we summarize our results.

Transverse-momentum resummation for vector-boson pair production
In this Section we recall the main points of the transverse-momentum resummation formalism we use in this paper. For a more detailed discussion the reader is referred to Refs. [59,64,65].
We consider the inclusive hard-scattering process where the collision of the two hadrons h 1 and h 2 with momenta P 1 and P 2 produces the two vector bosons of momenta p 3 and p 4 . In the center-of-mass frame the momentum of the vector-boson pair q = p 3 + p 4 is fully specified by the invariant mass M 2 = (p 3 + p 4 ) 2 , the rapidity y = 1 2 ln q·P 1 q·P 2 , and the transverse-momentum vector p T .
The kinematics of the vector bosons is fully determined by the vector-boson pair momentum q µ = p µ 3 + p µ 4 (with p 2 3 = m 2 V and p 2 4 = m 2 V ), and by the additional and independent variables that specify the angular distribution of the vector bosons with respect to q µ . Throughout this paper we always consider quantities which are inclusive over these angular variables.
According to the QCD factorization theorem, the differential cross section can be written as where f a/h (x, µ 2 F ) (a = q,q, g) are the density functions of parton a in hadron h at the factorization scale µ F ; µ R is the renormalization scale § ; dσ V V a 1 a 2 is the partonic cross section. The rapidityŷ and the center-of-mass energyŝ of the partonic scattering process are related to the corresponding § Throughout the paper we use parton densities in the MS factorization scheme and α S (q 2 ) is the QCD running coupling in the MS renormalization scheme.
hadronic variables y and s byŷ When the transverse momentum p T of the vector-boson pair is of the same order as the invariant mass M , the QCD perturbative expansion is controlled by a single expansion parameter, α S (M ), and fixed-order calculations can be safely applied. In this region, QCD radiative corrections are known to O(α 2 S ) [21,22,[32][33][34]. When p T M the convergence of the perturbative expansion is spoiled by the presence of large logarithmic terms of the form α n S ln m (M 2 /p 2 T ), that need to be resummed to all orders.
The resummation is performed at the level of the partonic cross section, which is decomposed as The first term on the right-hand side of Eq. (4) contains all the logarithmically-enhanced contributions at small p T and has to be evaluated by resumming them to all orders. The second term is instead free of such contributions and can thus be evaluated at fixed order in perturbation theory.
Resummation is based on the factorization of soft and collinear radiation and is viable in the impact-parameter (b) space, where the kinematical constraint of momentum conservation and the factorization of the phase space can be consistently taken into account [66][67][68]. Using the Bessel transformation between the conjugate variables p T and b, the resummed component is expressed as [59,68] dσ where J 0 (x) is the 0-order Bessel function. In the case of fully inclusive p T resummation, the rapidity dependence is integrated out in Eq. (5). In that case it is convenient to consider Mellin moments with respect to the variable z = M 2 /ŝ. However, in order to retain the rapidity dependence in the resummed cross section we take the 'double' (N 1 , N 2 ) Mellin moments with respect to the variables z 1 = e +ŷ M/ √ŝ and z 2 = e −ŷ M/ √ŝ at fixed M , and organize the structure of W V V in the following exponential form [64], where we have defined the logarithmic expansion parameter L as and b 0 = 2e −γ E (γ E = 0.5772... is the Euler number). The scale Q appearing in Eqs. (7,8), named resummation scale in Ref. [59], parametrizes the arbitrariness in the resummation procedure, and has to be chosen of the order of the hard scale M . Variations of Q around its reference value can be exploited to estimate the size of yet uncalculated higher-order logarithmic contributions. Therefore, Q plays a role very similar to µ F and µ R for missing perturbative terms. The function does not depend on the impact parameter b, and therefore includes all the perturbative terms that behave as constants as b → ∞. It can thus be expanded in powers of α S = α S (µ 2 R ): where σ V V ,(0) is the partonic leading-order (LO) cross section. The exponent G (N 1 ,N 2 ) includes the complete dependence on b and, in particular, it contains all the terms that order-by-order in α S are logarithmically divergent as b → ∞. The logarithmic expansion of G N reads where the term L g (1) collects the LL contributions, the function g (N 1 ,N 2 ) controls the NNLL terms and so forth.
The resummation of the large logarithmic terms carried out in Eq. (7), after transforming back to p T space, allows us to obtain a well behaved transverse-momentum spectrum as p T → 0. However, the logarithmic expansion parameter L in Eq. (8) is divergent as b → 0. This implies that the resummation produces higher-order contributions also in the high-p T region, which is conjugated to b → 0 after Fourier transformation. In this region the fixed-order cross section is perfectly viable and any resummation effect is necessarily artificial. To reduce the impact of such contributions, the logarithmic variable L is replaced by [59] L →L,L ≡ ln The variables L andL are equivalent when Qb 1 but they have a very different behaviour as b → 0. When Qb 1,L → 0 and G (N 1 ,N 2 ) → 1. Moreover, since the behaviour of the Sudakov form factor at b = 0 is related to the integral over p T , the replacement in Eq. (11) allows us to enforce a unitarity constraint such that the fixed-order prediction is recovered upon integration over p T .
A well known property of the formalism of Ref. [59] is that the process dependence (as well as the factorization scale and scheme dependence) is fully encoded in the hard function H V V . In other words, the functions g (i) are universal: they depend only on the channel in which the process occurs at Born level (qq annihilation in the case of vector-boson pair production). Their explicit expressions up to i = 3 are given in Ref. [59] in terms of the universal perturbative coefficients A q ,B (1) q,N ,B (2) q,N . In particular, the LL function g (1) depends on the coefficient A (1) q , the NLL function g (2) (N 1 ,N 2 ) also depends on A (2) q andB (1) q [69] and the NNLL function g The hard coefficients H V V depend on the process we want to consider. The first order coefficients H V V , (1) are known since long time [72,73]: they can be obtained from the one-loop scattering amplitudes qq → V V by using a process independent relation. By exploiting the expressions of H (2) for single Higgs [74] and vector-boson [75] production, in Ref. [65] such a relation has been extended to O(α 2 S ): this implies that the two-loop amplitude for qq → V V is the only process-dependent information needed to obtain the coefficient H V V , (2) .
We now turn to the finite component of the transverse-momentum spectrum, i.e. the second term on the right hand side of Eq. (4). Since dσ V V ,(fin.) ab does not contain large logarithmic terms in the small-p T region, it can be evaluated by truncating the perturbative series at a given fixed order. In practice, the finite component is computed starting from the customary fixed-order (f.o.) perturbative truncation of the partonic cross section and subtracting the expansion of the resummed cross section in Eq. (5) at the same perturbative order: At least formally, this matching procedure between resummed and finite contributions guarantees to achieve a uniform theoretical accuracy over the entire range of transverse momenta. At large values of p T , the resummation procedure cannot improve the fixed-order result, and the resummation (and matching) procedure is eventually superseded by the customary fixed-order calculations.
In summary, the inclusion of the functions g (1) , g (1) in the resummed component, together with the evaluation of the finite component at NLO (i.e. O(α S )), allows us to perform the resummation at NLL+NLO accuracy. This is the theoretical accuracy of the calculations of Refs. [55,56]. Including also the functions g (3) (N 1 ,N 2 ) and H V V , (2) , together with the finite component at NNLO (i.e. O(α 2 S )) leads to full NNLL+NNLO accuracy. Using the recently computed two-loop amplitudes for qq → V V , and the process independent relation of Ref. [65], we are now able to present the complete result for the transverse-momentum distribution of the vector-boson pair up to NNLL+NNLO accuracy. We point out that the NNLL+NNLO (NLL+NLO) result includes the full NNLO (NLO) perturbative contribution in the entire p T range. In particular, the NNLO (NLO) result for the double differential cross section dσ/(dM 2 dy) is exactly recovered upon integration over p T of the differential cross section at NNLL+NNLO (NLL+NLO) accuracy.
We conclude this Section by adding a few comments on the way in which our calculation is actually performed. The practical implementation of Eq. (4) is done in the numerical program Matrix ¶ , which is an extension of the numerical program applied in the NNLO calculations of Refs. [23,36,76,77] and based on a combination of the q T -subtraction formalism [63] with the Munich code. Since already performed within the q T -subtraction formalism, the extension of these calculations to compute the resummed cross section is conceptually quite straightforward, and is obtained by replacing the hard-collinear terms in the fixed-order computation by the proper all-order resummation formula of Eq. (5). This procedure is the same that was applied to perform the NLL+NLO calculations for W + W − [55] and ZZ [56] production, and the NNLL+NNLO calculation for Higgs boson production of Ref. [60]. ¶ Matrix is the abbreviation of "Munich Automates qT subtraction and Resummation to Integrate X-sections", by M. Grazzini, S. Kallweit, D. Rathlev, M. Wiesemann. In preparation.
Munich is the abbreviation of "MUlti-chaNnel Integrator at Swiss (CH) precision"-an automated parton level NLO generator by S. Kallweit. In preparation.
To obtain the numerical results presented here, the resummed component of Eq. (5) is evaluated with an extension of the numerical program used for the calculation of Higgs production [60], based on the earlier computations of Refs. [59,64]. The hard-collinear coefficients are obtained by exploiting the implementation of the corresponding virtual amplitudes for the production of on-shell ZZ and W + W − pairs in Ref. [23] and Ref. [36], respectively, and the knowledge of the collinear coefficients relevant to quark-initiated processes [75]. * * The finite component of Eq. (12) is obtained from an NLO calculation of V V +jet, computed with the Munich code, which provides a fully automated implementation of the Catani-Seymour dipole formalism [80,81] as well as an interface to the one-loop generator OpenLoops [82] to obtain all required (spin and color-correlated) tree-level and one-loop amplitudes. For the numerically stable evaluation of tensor integrals we rely on the Collier library [83], which is based on the Denner-Dittmaier reduction techniques [84,85] and the scalar integrals of [86]. To deal with problematic phase-space points, OpenLoops provides a rescue system using the quadrupleprecision implementation of the OPP method in CutTools [87], involving scalar integrals from OneLOop [88].

Results
In this Section we present our results for the resummed transverse-momentum distributions of W + W − and ZZ pairs. We compare our NNLL+NNLO predictions to the results at the NLL+NLO, and discuss the corresponding theoretical uncertainties. Additionally, we also study the rapidity dependence of the p T cross section as well as the p T -veto efficiency.
For the EW couplings we use the so-called G µ scheme, where the input parameters are G F , m W , m Z . In particular we set G F = 1.16639 × 10 −5 GeV −2 , m W = 80.399 GeV, m Z = 91.1876 GeV. We use the NNPDF3.0 sets of parton distribution functions (PDFs) [89] with α S (m Z ) = 0.118. At NLL+NLO and NNLL+NNLO the running of α S is evaluated at two-and three-loop order, respectively. For ZZ production we consider N f = 5 massless quarks/antiquarks. For W + W − production we make use of the 4FS, which allows us to split off all contributions related to bottom-quark final states in order to remove the tt and W t contamination from our computation. This is straightforward in the 4FS, because the W + W − bb process is separately finite.
We consider proton-proton collisions at √ s = 8 TeV. The central values of the factorization and renormalization scales are set to µ F = µ R = µ 0 = 2 m V . The choice of the central resummation scale Q 0 is discussed in the next subsection.

Choice of the central resummation scale
As discussed in Section 2, the resummation scale Q is the scale entering the large logarithmic terms we are resumming (see Eq. (8)), and it plays the role of the scale up to which resummation is effective. In on-shell Higgs [59] and vector-boson [90] production, the scale is typically chosen * * We note that an independent computation of these coefficients in the framework of Soft Collinear Effective Theory has been presented in Refs. [78,79]. The following considerations apply both to ZZ and W + W − production, so we will focus on W + W − production from now on. In Fig. 1 we consider the invariant mass distribution of the W + W − pair at NLO and NNLO. We see that the distribution is strongly peaked in the threshold region, and that it quickly decreases as M W W increases. As a consequence, for most of the W + W − events, M W W 2m W .
We can compare the transverse-momentum distributions obtained with a dynamical resummation scale Q = M W W /2, and a fixed resummation scale Q = m W . In Fig. 2 we show the ratio (blue, solid curve) of the Q = m W result over the Q = M W W /2 result. The bands are obtained by varying the resummation scale around the central value by a factor of two. Considering the ratio of the central curves for p T 250 GeV, the differences between a fixed and a dynamical scale are extremely small and remain at the 1-2% level over the whole range. In this region of transverse momenta the uncertainty bands obtained with the two choices overlap and are similar in size. In fact, since Q = m W leads to slightly larger uncertainties, it appears to be the more conservative choice. Therefore, we can conclude that either choice of the resummation scale is perfectly valid and indeed consistent with each other as expected from the discussion of the invariant mass distribution.
Looking further at the comparison of the high-p T tails in Fig. 2   for large values of the resummation scale the fixed-order cross section (black dotted curve) is not recovered in the tail of the distribution. It is important to recall that transverse-momentum resummation is supposed to improve the perturbative expansion in the low-p T region. At large p T , any large dependence on the resummation scale is necessarily artificial and an unwanted remnant of the matching procedure. This behaviour is precisely what we observe for the dynamical scale choice Q = M W W /2 in Fig. 2. With this choice, in fact, the resummed result loses predictivity, as its uncertainty becomes increasingly large. By contrast, a fixed resummation scale Q = m W , which is always smaller than Q = M W W /2, eventually leads to a more consistent high-p T behaviour of the resummed prediction.
Based on the above results, we make Q 0 = m V our default choice of the resummation scale in what follows.

Inclusive transverse-momentum distribution
We now present our resummed predictions for the inclusive transverse-momentum spectrum of the vector-boson pair and compare them with the corresponding fixed-order results. We concentrate on W + W − production since we observe no saliently different features in the ZZ case. For completeness, we provide the corresponding reference prediction with uncertainties for ZZ below.
Before presenting our resummed predictions, we recall the well known fixed-order results at O(α S ) and O(α 2 S ) [32][33][34]. In Fig. 3  GeV to about 30% at p T ∼ 400 GeV. The NLO (NNLO) uncertainty ranges from about ±15% (±10%) at p T ∼ 50 GeV to about ±20% (±8%) at p T ∼ 400 GeV. We note that the NLO and NNLO bands do not overlap in the region where p T 300 GeV. This implies that, in this region of transverse momenta, the size of the band obtained through scale variations at NLO definitely underestimates the theoretical uncertainty.
We now move on to the resummed results. In Fig. 4 (a) the NLL+NLO spectrum is compared to the fixed-order NLO result and to the finite component of the resummed cross section (see Eq. (4)) in the region between 0 and 80 GeV. As expected, the NLO diverges to +∞ as p T → 0, while the resummation provides a physically well behaved spectrum down to low values of p T , which exhibits a kinematical peak at p T ∼ 4 GeV. The finite component contributes less than 1% in the peak region, where the result is dominated by resummation, and it increases to ∼ 18% at p T = 50 GeV. The lower inset shows the NLL+NLO result normalized to NLO. In Fig. 4 (b) the region between 80 and 400 GeV is displayed. We see that even at large values of p T the NLL+NLO resummed result does not match very well the fixed-order NLO result, with a difference of about 5%.
The analogous results at NNLL+NNLO are shown in Fig. 5. The NNLO has an unphysical (divergent) behaviour as p T → 0, whereas the resummed spectrum is well behaved, with a slightly harder peak with respect to the NLL+NLO.  the peak region, increasing to ∼ 19% at p T = 50 GeV. Comparing the right panels of Fig. 4 and Fig. 5, we see that the quality of the matching at high p T is significantly improved when going from NLL+NLO to NNLL+NNLO, and we find that this behaviour is indeed preserved up to very high transverse momenta. The NNLL+NNLO result thus gives a prediction with uniform accuracy from small to very large transverse momenta and, in fact, provides a sufficiently large region where a hard switching to the fixed-order result is feasible. We point out that, thanks to our unitarity constraint, both at NLL+NLO and at NNLL+NNLO the integral of the resummed spectrum is in excellent agreement with the respective total cross sections; the differences are at the few-permille level.
We now turn to the scale uncertainties of our resummed results. We start our discussion by separately considering factorization and renormalization scale variations. In Fig. 6 we compare the NLL+NLO (red, dashed) and NNLL+NNLO (blue, solid) predictions with their uncertainty bands from µ F and µ R variations (left and right panel, respectively). In both cases, the bands are obtained by varying the factorization (renormalization) scale by a factor of two around its central value, while keeping the other scales at their default values. First of all, we notice that when going from NLL+NLO to NNLL+NNLO the p T spectrum becomes harder. Comparing with the results of Ref. [55], where the NNLL resummation was implemented without O(α 2 S ) matching, we see that the increased hardness of the p T spectrum is a combined effect of both features, i.e. NNLL resummation and NNLO matching at high p T .
We note that neither in the case of the factorization scale, nor in the case of the renormalization scale, the NLL+NLO and NNLL+NNLO bands overlap. Actually, in the case of the factorization scale, there is no reduction in scale dependence when going from NLL+NLO to NNLL+NNLO, and the uncertainty slightly increases with the perturbative order, even if it is always well below 10%, except at very low p T . The renormalization scale dependence instead exhibits the expected reduction when going from NLL+NLO to NNLL+NNLO.
In Fig. 7 we present our resummed predictions with uncertainty bands obtained from simultaneous variations of µ F and µ R (left panel) and the variation of the resummation scale Q (right panel). In the left panel the uncertainty bands are obtained by varying µ F and µ R as in Fig. 3. In the right panel the resummation scale is varied in the range m W /2 ≤ Q ≤ 2m W . As in Fig. 6 we see that the uncertainty bands do not overlap. The uncertainty from µ F and µ R variations is ±10 − 15% at NLL+NLO and is reduced to 8 − 10% at NNLL+NNLO. At NLL+NLO the resummation scale uncertainty is generally about ±15% except in the region of p T ∼ 10 GeV, where it shrinks to smaller values. We find that at NNLL+NNLO the resummation scale uncertainty is reduced roughly by a factor of two in the region of transverse momenta considered in the figure.
In Figs. 8 and 9 we show our reference resummed prediction for W + W − and ZZ, respectively, with an estimate of their full perturbative uncertainty. In order to obtain a combined uncertainty from µ F , µ R and Q variations, we follow Ref. [90] and independently vary µ F , µ R and Q in the ranges m V ≤ {µ F , µ R } ≤ 4m V and m V /2 ≤ Q ≤ 2m V with the constraints 0.5 ≤ µ F /µ R ≤ 2 and 0.5 ≤ Q/µ R ≤ 2. We recall that the constraint on µ F /µ R , which is the same as applied in Fig. 3 and Fig. 7 (left), has the purpose of avoiding large logarithmic contributions from the evolution of parton densities. Analogously, the constraint on Q/µ R avoids large logarithmic contributions in the expansion of the Sudakov form factor.
For W + W − production the perturbative uncertainty at NNLL+NNLO (NLL+NLO) is about ±8% (±12%) at the peak, it decreases to about ±3% (±5%) at p T = 20 GeV, and it increases again to ±10% (±15%) at p T = 200 GeV. In the high-p T region, the difference between the NNLL+NNLO and NLL+NLO predictions is driven by the NNLO effects, which increase the NLO result by about 30%.
For ZZ production the uncertainties have essentially the same pattern in the small-and intermediate-p T region, while at high p T they are larger than for W + W − production, reaching about ±17% at NNLL+NNLO for p T = 200 GeV. We have checked that this effect is entirely driven by the resummation-scale dependence. As previously pointed out, this behaviour is not particularly worrying since, in the large-p T region, the resummed results should be replaced by the corresponding fixed-order prediction. Also in the ZZ case the large enhancement of the NNLL+NNLO distribution in the high-p T tail stems from the fixed-order cross section. (a) (b) Figure 10: (a) Shapes of the W + W − transverse-momentum distribution differential in the rapidity of the W + W − pair at the NNLL+NNLO for |y| < 0.5 (red, solid), 0.5 < |y| < 1 (blue, dashed), 1 < |y| < 2 (black, dotted), 2 < |y| < 3 (magenta, dash-dotted), 3 < |y| (orange, double-dash dotted); and (b) the shape-ratio with respect to the inclusive result.

Rapidity dependence of the transverse-momentum distribution
So far, we only considered p T spectra for on-shell W + W − and ZZ production that are inclusive in the kinematics of the vector-boson pair. Our numerical program, however, allows us to compute arbitrary observables that are differential with respect to the V V phase space.
In the following we study the behaviour of the transverse-momentum spectrum in different rapidity regions of the vector-boson pair. In Fig. 10 we study the shape of the NNLL+NNLO transverse-momentum distribution, i.e. normalized such that its integral yields one, for |y| < 0.5 (red, solid), 0.5 < |y| < 1 (blue, dashed), 1 < |y| < 2 (black, dotted), 2 < |y| < 3 (magenta, dash-dotted) and 3 < |y| (orange, dash-double dotted). The right panel shows the same results normalized to the fully inclusive distribution. We clearly see that the p T shapes become softer as the rapidity increases. In the central region (|y| < 2) the distributions are still quite insensitive to the specific value of the rapidity and only slightly harder than the inclusive spectrum. In the forward rapidity region, on the other hand, the shapes become increasingly softer.
The observed pattern can be understood in the following way: rapidity and transverse momentum are two not completely independent phase-space variables. Indeed, they affect their mutual upper integration bounds. At higher rapidities the kinematically allowed range of transverse momenta is reduced: this squeezes the p T spectrum which consequently becomes softer. This effect has been observed also in previous studies in the case of Higgs boson production [64].  The excess in the W + W − production cross section measured by ATLAS [51] with respect to the SM prediction has drawn a lot of attention to the W + W − process, since the W + W − signature appears in many new physics scenarios [53]. The inclusion of the recently computed NNLO corrections [36] considerably reduces the significance of the excess. However, particular attention must be payed to the modelling of the jet veto [52, 58,93] when extrapolating from the fiducial region to obtain the inclusive cross section. Effects of jet veto resummation have been considered in Refs. [94,95], though still matching to the fixed-order O(α S ) result.

The
In this paper we are dealing with transverse-momentum spectra, and we perform a resummation on a different variable with respect to the jet p T . However, the vector-boson pair p T and the jet p T are clearly related variables (actually, at O(α S ) they indeed coincide). We will therefore study the p T -veto efficiency in W + W − production at different orders in resummed and fixed-order perturbation theory. We define the p T -veto efficiency as In Fig. 11 we show (p veto T ) at the NNLL+NNLO (blue, solid), approximate NNLL+NLO (magenta, dash-double dotted), NLL+NLO (red, dashed), NNLO (black, dotted) and NLO (grey, dash-dotted). The lower inset shows the same curves normalized to our reference prediction at NNLL+NNLO. Our approximate NNLL+NLO is obtained by simply adding the g (3) function in the Sudakov exponent in Eq. (10) at NLL+NLO, and corresponds to the approximation considered in Refs. [55,58].
For reference, the corresponding numerical values of the efficiencies are given in Tab. 1 for p T = 5-40 GeV. The uncertainty bands are obtained by a combined variation of resummation, factorization and renormalization scales as in Fig. 8. The first thing we observe is that the NLO result appears to be well above the others and cannot be really considered a reliable prediction for the efficiency. This is because it is essentially a LO prediction at finite values of p veto T . We also note that in the small-p T region (say below p T ∼ 10 GeV) the fixed-order NLO and NNLO predictions diverge and cannot be trusted. Comparing further the fixed-order results among each other and the resummed results among each other, we observe that higher-order corrections in fixed-order and resummed perturbation theory reduce the p T -veto efficiency.
Both effects can be easily understood in the light of the results presented up to now. As seen in Fig. 3, the inclusion of the NNLO corrections make the p T distribution harder. Furthermore, resummation effects generally harden the spectrum. A qualitatively similar result is obtained when going from NLL+NLO to NNLL+NNLO (see Fig. 8).
It is interesting to compare the approximated NNLL+NLO result with the NNLO and NNLL+NNLO predictions. For values of p veto T ∼ 25 − 30 GeV we see that the approximated result is in between the NNLO one and our best NNLL+NNLO prediction. This means that the effect of NNLL resummation obtained by the inclusion of the g (3) function in the Sudakov exponent in Eq. (10) is quantitatively important. Nonetheless, the efficiency obtained within this approximation is still about 5% higher than the NNLL+NNLO prediction. We also notice that in this region of p veto T , the NNLO and NLL+NLO results differ by less than 1%.
Comparing the NNLL+NNLO and NLL+NLO results, we find that they are compatible within the corresponding uncertainties.
We add few comments on the recent measurement of the W + W − cross section carried out by the CMS collaboration [52]. The result shows good agreement with the NNLO prediction of Ref. [36]. The corresponding analysis, however, is based on a reweighting procedure of the p T spectrum of the W + W − pair. The events generated with POWHEG [96] plus Pythia6 [97] were reweighted by using the calculation of Ref. [58], which corresponds to our NNLL+NLO approximation, and includes neither the second-order hard-collinear coefficient H W W, (2) in Eq. (9), nor the NNLO matching. The results in Fig. 11 show that the NNLL+NNLO p T -veto efficiency is lower than

Summary
In this paper we have studied the transverse-momentum distribution of vector-boson pairs in hadronic collisions. We presented a computation of the p T spectrum in which the logarithmically enhanced contributions at small p T are resummed up to NNLL accuracy and the ensuing result is combined with state-of-the-art O(α 2 S ) (NNLO) predictions valid at large p T . We presented numerical results for W + W − and ZZ production at the LHC together with a study of their perturbative uncertainties.
We found that, up to relatively large transverse momenta, when scale variations are studied around the fixed resummation scale Q = m V we obtain results which are fully consistent with those obtained using the dynamical choice Q = M V V /2. At large transverse momenta the fixed-order result in the tail of the distribution is nicely recovered with the fixed-scale choice. Our new NNLL+NNLO results significantly reduce the theoretical uncertainties obtained through scale variations compared to lower orders in both the peak region and the tail of the distribution.
We have also studied the rapidity dependence of the resummed transverse-momentum distribution. The rapidity dependence at NNLL+NNLO is quite flat in the central region (|y| 2), but signals a substantially softer spectrum in the forward region. Due to phase-space suppression, the effect on the inclusive transverse-momentum distribution is very moderate though.
Finally, we have studied the p T -veto efficiency at different orders in resummed and fixed-order perturbation theory. Both NNLL resummation and the NNLO effects turned out to be important to obtain an accurate prediction for this quantity. We observed that the veto efficiency at NNLL+NNLO is ∼ 5% lower (∼ 3% in absolute terms) with respect to the approximate NNLL+NLO calculation used in the CMS analysis of Ref. [52]. This result suggests that our NNLL+NNLO predictions will be useful to validate the transverse-momentum spectra obtained from Monte Carlo event generators, similarly to what was done for the NNLL+NNLO calculation of Ref. [59] in the case of Higgs boson production.
In this paper we considered the p T spectrum of stable vector-boson pairs. Exploiting the twoloop helicity amplitudes for qq → V V → 4 leptons [37,38] will allow us to extend the calculation to include the leptonic decay of the vector bosons and off-shell effects. The computation of the transverse-momentum spectrum with realistic experimental cuts will then become possible.  [52] CMS Collaboration, Measurement of the W + W − cross section in pp collisions at sqrt(s) = 8 TeV and limits on anomalous gauge couplings, CMS-PAS-SMP-14-016.