1 Introduction

In many experimental measurements a veto on hard jets is imposed to suppress backgrounds. Such a veto is particularly useful to suppress top-quark backgrounds to processes involving \(W\) bosons, since the \(W\) bosons from the decay of the top quarks come in association with \(b\)-jets, which are rejected by the jet veto. For example, a jet veto is crucial to measure Higgs production with subsequent decay \(H\rightarrow W^+W^-\). It is imposed by rejecting events which involve jets with transverse momentum above a scale \(p_T^{{\mathrm{veto}}}\), which is typically chosen to be \(p_T^{{\mathrm{veto}}}\approx \) 20–30 GeV. Since the veto scale is much lower than the invariant mass \(Q\) of the electroweak final state, perturbative corrections to the cross section are enhanced by Sudakov logarithms of the ratio \(p_T^{{\mathrm{veto}}}/Q\). There has been a lot of theoretical progress over the past two years concerning the resummation of jet-veto logarithms in Higgs-boson production. Using the CAESAR formalism [1], these logarithms were first computed at next-to-leading logarithmic (NLL) order in [2], and this treatment was later extended to next-to-next-to-leading logarithmic (NNLL) accuracy [3]. In between these papers, an all-order factorization formula derived in soft-collinear effective theory (SCET) [4–7] was proposed [8], and a resummed result which includes almost all of the ingredients required for N\(^3\)LL accuracy was presented [9]. A third group of authors performed an independent analysis in SCET [10] and also combined the results for different jet multiplicities [11–13].

The jet veto is not only necessary in \(H\rightarrow W^+ W^-\) but also in the measurement of the diboson cross section itself. The fact that LHC measurements [14–17] yield values of the \(W^+ W^-\) cross section that are higher than theoretical predictions has triggered discussions as to whether this excess could be due to New Physics [18–20]. To be sure whether there indeed is an excess, it is important to have reliable theoretical predictions not only for the total cross section, for which the next-to-next-to-leading-order (NNLO) result has been obtained recently [21], but also for the cross section in the presence of experimental cuts, most importantly in the presence of a jet veto.Footnote 1 Several recent papers have addressed this issue and have come to somewhat different conclusions. In [23], the Sudakov logarithms associated with the jet veto were resummed at NNLL accuracy. It was claimed that resummation effects increase the cross section and bring the Standard-Model prediction in agreement with the experimental measurements. On the other hand, based on a study of transverse-momentum resummation, the authors of [24] concluded that resummation effects are small for the relevant values of \(p_T^{{\mathrm{veto}}}\). Most recently, the effect of using a matched parton shower to predict the fiducial cross section, as is done in the experimental analyses, was analyzed in [25]. These authors concluded that resummation effects are small and that a fixed-order computation of the fiducial rate would lead to theoretical predictions in agreement with the measurements, but that the matched parton shower overestimates the Sudakov suppression of the rate and leads to systematically lower theoretical predictions when extrapolating back to the total rate.

In the present paper, we present an automated method to perform resummations for arbitrary vector-boson-production processes involving jet vetoes. Instead of computing resummed cross sections analytically, on a case-by-case basis, we obtain them in an automated way using the MadGraph5_aMC@NLO framework [26]. Our method yields results which are accurate at NNLL and are matched to NLO fixed-order results. Such an automated procedure is obviously much more efficient and less error prone than computing the ingredients by hand or extracting them from the literature. Most importantly, our approach allows us to also include the decay of the vector bosons, along with cuts on the leptons in the final state.

We have implemented two different methods to perform the resummation. The first one is based on reweighting tree-level events generated by MadGraph. It yields jet-veto cross sections accurate at NNLL order. The event weight includes universal resummation factors as well as the process-specific one-loop virtual corrections, which are computed using Madgraph5_aMC@NLO. In the second method, we modify the NLO fixed-order computation in such a way that the end result is accurate at both NNLL and NLO. In this second method not only the hard function, which encodes the virtual corrections, but also the beam functions, which describe the emissions at small transverse momentum, are computed by Madgraph5_aMC@NLO.

Our paper is organized as follows. We start in Sect. 2 by reviewing the resummation formula for cross sections in the presence of a jet veto. We also discuss non-perturbative corrections and point out that they could be sizable, similar in magnitude as the recently calculated NNLO corrections. We then explain in Sect. 3 how the automated resummation can be implemented in the Madgraph5_aMC@NLO framework. In Sect. 4 we use our method to compute cross sections for different boson-production processes and discuss in detail the scale and scheme choices and the resulting theoretical uncertainties. We compare our resummed predictions to fixed-order results for the cross sections, to the results obtained from a matched parton shower, and to the NNLL results of [23]. We also match our resummed result to fixed-order NLO predictions. The relevant matching corrections turn out to be very small, which indicates that the bulk of the NLO result is already captured by the factorization formula evaluated with NNLL accuracy. This remains true after imposing cuts on the leptonic final state in the decays of the electroweak bosons. We compare predictions for the final states \(Z\), \(W^+W^-\), and \(W^+W^-W^+\) and consider ratios of cross sections, which have small uncertainties if they are properly defined. We then discuss the implications of our results on the value of the \(W^+W^-\) cross section and conclude in Sect. 5.

2 Factorization theorem for jet-veto cross sections

We focus on electroweak-boson-production processes with a veto on jets with transverse momentum above a cut \(p_T^{{\mathrm{veto}}}\). The large logarithms which arise in the presence of the jet veto have the form \(\alpha _s^n \ln ^m(p_T^{{\mathrm{veto}}}/Q)\) with \(m\le 2n\), where \(Q\) denotes the invariant mass of the boson system. Our goal is the resummation of these logarithms to all orders in perturbation theory and at leading power in the small ratio \(p_T^{{\mathrm{veto}}}/Q\). For concreteness, we will discuss the resummation for \(W^+ W^-\) pair production in the following, but the formalism applies to any number of massive vector bosons and Higgs bosons or other massive color-singlet particles in the final state. The resummation is based on a factorization theorem which arises in the limit \(p_T^{{\mathrm{veto}}}/Q \rightarrow 0\) [8]. Its schematic form is shown in Fig. 1. The main ingredients of the theorem are hard functions \(\mathcal{H}_{ij}\), which encode the virtual QCD corrections to the partonic hard-scattering processes \(i+j\rightarrow W^+ W^-\), and two beam functions \(\bar{B}_i\) and \(\bar{B}_j\), which describe the low-\(p_T\) emissions collinear to the two beams.

Fig. 1
figure 1

Structure and kinematics of the factorization theorem for the \(W^+ W^-\)-production cross section in the presence of a jet veto

Full size image

Before writing out the factorization theorem in more detail, let us specify the kinematics of the process at low \(p_T^{{\mathrm{veto}}}\). The momenta of the incoming protons are \(p_1\) and \(p_2\). The partons emerging from the parton distribution functions (PDFs) carry momenta \(z_1 p_1\) and \(z_2 p_2\). After possible emissions (described by the beam functions \(\bar{B}_i\)), the momenta \(\xi _1 p_1\) and \(\xi _2 p_2\) are left to produce the boson pair through a hard interaction \(\mathcal {H}_{ij}\). In the limit of small transverse momenta we can neglect recoil effects, so that the partons are still collinear to the proton momentum after the emissions. We define

$$\begin{aligned} \hat{s}= & {} (q_1+q_2)^2=(\xi _1 p_1+\xi _2 p_2)^2=Q^2 , \nonumber \\ \hat{t}= & {} (\xi _1 p_1-q_1)^2 , \quad \hat{u} = (\xi _1 p_1-q_2)^2 , \end{aligned}$$
(1)

with \(\hat{s}+\hat{t}+\hat{u}=2M_W^2\). Note that our definition of the variable \(\hat{s} \) differs from the standard choice \((z_1 p_1+ z_2 p_2)^2\). The quantity \(\hat{s}\) we define is the one relevant for the boson-production process, i.e. the one that enters the hard function. In the small transverse-momentum limit of the emissions, we obtain

$$\begin{aligned} \xi _1= & {} \frac{\bar{n}\cdot q}{\bar{n}\cdot p_1}=\frac{Q}{\sqrt{s}}\,e^{-y}\Rightarrow \xi _1p_1=\left( \bar{n}\cdot q\right) \frac{n}{2}\nonumber \\ \xi _2= & {} \frac{n\cdot q}{n\cdot p_2}=\frac{Q}{\sqrt{s}}\,e^{y}\Rightarrow \xi _2p_2=\left( n\cdot q\right) \frac{\bar{n}}{2}, \end{aligned}$$
(2)

where \(n^\mu =(1,0,0,1)\) and \(\bar{n}^\mu =(1,0,0,-1)\) are two light-cone vectors in the beam directions, \(y\) denotes the rapidity of \(q=q_1+q_2\) in the laboratory frame, and \(s=(p_1+p_2)^2\). The crucial feature of (2) is that it shows that one can obtain the arguments of the hard function directly from the vector-boson (and proton) kinematics. The same is true for an arbitrary electroweak final state.

At low \(p_T^{{\mathrm{veto}}}\), the differential cross section in the presence of a jet veto has the factorized form [8, 9]

$$\begin{aligned} \frac{d^3\sigma (p_T^{{\mathrm{veto}}})}{dy\,dQ^2\,d\hat{t}}= & {} \sum _{i,j=g,q,\bar{q}} \sigma ^0_{ij}(Q^2,\hat{t},\mu )\,P_{ij}(Q^2,\hat{t},p_T^{{\mathrm{veto}}},\mu )\,\nonumber \\&\times \bar{B}_i(\xi _1,p_T^{{\mathrm{veto}}})\,\bar{B}_j(\xi _2,p_T^{{\mathrm{veto}}}) . \end{aligned}$$
(3)

Here \(i\) and \(j\) are the flavors of the partons which enter the hard-scattering process after initial-state radiation, and \(\sigma ^0_{ij}(Q^2,\hat{t})\) is the Born-level cross section for the production of the electroweak final state. Since the electroweak final state is a color singlet, we either deal with \(q\bar{q}\) or \(gg\). For \(W^+ W^-\) pair production at leading order only the quark channels contribute, but starting from NNLO also the gluon-induced reaction occurs.

The second ingredient in (3) are the beam functions \(\bar{B}_i(\xi ,p_T^{{\mathrm{veto}}})\), which are given by a convolution of a perturbative kernel \(\bar{I}_{q\leftarrow k}(z,p_T^{{\mathrm{veto}}},\mu )\) describing the emissions with the standard PDFs \(\phi _k\):

$$\begin{aligned} \bar{B}_i(\xi ,p_T^{{\mathrm{veto}}}) = \sum _{k=g,q,\bar{q}} \int _\xi ^1\!\frac{dz}{z}\, \bar{I}_{i\leftarrow k}(z,p_T^{{\mathrm{veto}}},\mu )\,\phi _k(\xi /z,\mu ) . \end{aligned}$$
(4)

The bar over these functions indicates that a factor \(e^{h_i(p_T^{{\mathrm{veto}}},\mu )}\) has been extracted from the original definitions of these functions in terms of SCET operators, called \(B_i\) and \(I_{i\leftarrow k}\), such that

$$\begin{aligned} \bar{B}_i(\xi ,p_T^{{\mathrm{veto}}}) = e^{-h_i(p_T^{{\mathrm{veto}}},\mu )} B_i(\xi ,p_T^{{\mathrm{veto}}},\mu ) , \end{aligned}$$
(5)

and analogously for \(\bar{I}_{i\leftarrow k}(z,p_T^{{\mathrm{veto}}},\mu )\). This factor is normalized such that \(h_i(p_T^{{\mathrm{veto}}},p_T^{{\mathrm{veto}}})=1\), and chosen such that the remaining function \(\bar{B}_i(\xi ,p_T^{{\mathrm{veto}}})\) is renormalization-group (RG) invariant. The explicit form of \(h_i(p_T^{{\mathrm{veto}}},\mu )\) as well as the one-loop kernels \(\bar{I}_{i\leftarrow k}(z,p_T^{{\mathrm{veto}}},\mu ) \) are listed in the appendix.

The final ingredient in (3) is the factor \(P_{ij}(Q^2,\hat{t},p_T^{{\mathrm{veto}}},\mu )\), which includes the hard function and the resummation of large logarithms. It has the form

$$\begin{aligned} P_{ij}(Q^2,\hat{t},p_T^{{\mathrm{veto}}},\mu ) = \mathcal {H}_{ij}(Q^2,\hat{t},\mu _h)\, E_{i}(Q^2,p_T^{{\mathrm{veto}}},\mu _h,\mu , R),\nonumber \\ \end{aligned}$$
(6)

where the hard function

$$\begin{aligned} \mathcal {H}_{ij}(Q^2,\hat{t},\mu _h) = 1+\frac{\alpha _s(\mu _h)}{4\pi }\,\mathcal {H}_{ij}^{(1)}(Q^2,\hat{t},\mu _h) + \cdots \end{aligned}$$
(7)

contains higher-order finite virtual corrections to the Born-level cross section. Since these higher-order corrections contain (double) logarithms of \(Q/\mu _h\), the hard matching scale \(\mu _h\) should be chosen of order \(Q\). The evolution of the hard function to a lower scale \(\mu \ll Q\) is controlled by an RG evolution equation. The corresponding evolution function \(U_{i}(Q^2,\mu _h,\mu )\), together with the collinear anomaly [27] and the prefactors extracted from the beam functions, is absorbed into the factor \(E_i\) in (6). The collinear anomaly arises due to light-cone divergences and provides an additional source of large logarithms in processes sensitive to small transverse momenta. The explicit form of the quantity \(E_i\) reads

$$\begin{aligned}&E_{i}(Q^2,p_T^{{\mathrm{veto}}},\mu _h,\mu , R)\nonumber \\&\quad = U_{i}(Q^2,\mu _h,\mu ) \left( \frac{Q}{p_T^{{\mathrm{veto}}}} \right) ^{-2F_{i}(p_T^{{\mathrm{veto}}},\mu , R)}e^{2h_{i}(p_T^{{\mathrm{veto}}},\mu )} . \end{aligned}$$
(8)

The evolution factor at NNLL accuracy is given in the appendix. It differs for quark-initiated (\(i=q\)) and gluon-initiated (\(i=g\)) processes but is independent of the quark flavors. Note that the evolution factor depends on the kinematics of the final state only via the invariant mass \(Q\). The anomaly exponent \(F_{i}(p_T^{{\mathrm{veto}}},\mu , R)\) resums the large anomalous logarithms in the beam and soft functions, which arise from the rapidity difference between the modes which contribute to the individual functions [27–29]. Starting from two-loop order (which is needed for NNLL resummation) this exponent depends on the jet radius \(R\), but it is the same for any \(k_T\)-style sequential jet-clustering algorithm. The explicit form of the two-loop exponent can be found in the appendix. It was calculated in [9] and is related to the function \(\mathcal F\) obtained earlier in [3].

We stress that the factorization theorem holds up to power corrections suppressed by \(p_T^{{\mathrm{veto}}}/Q\), and up to non-perturbative effects suppressed by \(\Lambda _\mathrm{QCD}/p_T^{{\mathrm{veto}}}\). For the weak-boson transverse-momentum spectrum, these corrections depend on \(\mathbf {p}_T^2\) and hence are of second order in \(p_T^{{\mathrm{veto}}}/Q\) and \(\Lambda _\mathrm{QCD}/p_T^{{\mathrm{veto}}}\). The definition of the jet veto, on the other hand, involves an absolute value of the jet transverse momentum, and for this reason there can be first-order power corrections. Non-perturbative corrections to processes involving an anomaly were studied in [30], where it was found that these effects are enhanced by a logarithm of the rapidity difference between the left- and right-collinear emissions and can be viewed as a non-perturbative contribution to the anomaly exponent \(F_i\) in (8). The leading non-perturbative corrections to jet-veto cross sections are therefore expected to scale as

$$\begin{aligned} \sigma _\mathrm{NP}(p_T^{{\mathrm{veto}}})\sim \sigma ^0\times \frac{\Lambda _\mathrm{NP}}{p_T^{{\mathrm{veto}}}}\, \ln \frac{Q}{p_T^{{\mathrm{veto}}}} . \end{aligned}$$
(9)

Due to the fact that the correction is of first order and logarithmically enhanced, these effects might not be negligible. For example, assuming \(\Lambda _\mathrm{NP} = 0.5\,\mathrm{GeV}\) and \(p_T^{{\mathrm{veto}}}= 20 \,\mathrm{GeV}\), one ends up with a 6 % effect at \(Q=222\,\mathrm{GeV}\), which is the median \(Q\) value in \(W^+ W^-\) production. Numerically, this is not much smaller than the NNLO correction to the total cross section calculated in [21]. The value of the non-perturbative quantity \(\Lambda _\mathrm{NP}\) is unknown, but it could be obtained from the matrix element

$$\begin{aligned} \mathcal {M}_\mathrm{veto} = \sum \int \limits _{X,\mathrm{reg}} p_T^\mathrm{jet} \left| \langle X | S_n^\dagger (0)\,S_{\bar{n}}(0) |0\rangle \right| ^2 \end{aligned}$$
(10)

of two soft Wilson lines along the beam directions, where \(p_T^\mathrm{jet}\) is the transverse momentum of the leading jet in the final state \(X\). The phase-space integrals in the matrix element \(\mathcal {M}_\mathrm{veto}\) suffer from a rapidity divergence, which needs to be regularized. The parameter \(\Lambda _\mathrm{NP}\) multiplies the rapidity divergence (see [30] for more details). To get an idea of the size of non-perturbative effects, we have computed the hadronization effects on the cross section using Pythia 8 [31] with its default tune. We find that they change the cross section by about 10 % at \(p_T^{{\mathrm{veto}}}=10\,\mathrm{GeV}\) and 3 % at \(p_T^{{\mathrm{veto}}}=20\,\mathrm{GeV}\). Above \(p_T^{{\mathrm{veto}}}>20\,\mathrm{GeV}\), the simple parametrization in (9) with \(\Lambda _\mathrm{NP}=240\,\mathrm{MeV}\) provides a good description of the Pythia hadronization corrections, while a first-order power correction without logarithmic enhancement would underestimate the effects at higher \(p_T^{{\mathrm{veto}}}\) values. However, one should be careful in relying on Pythia hadronization effects in the context of precision calculations. There are other examples, such as the event-shape variable thrust, where Pythia appears to underestimate the size of these effects [32]. In the absence of a non-perturbative evaluation of the soft matrix element (10), the only reliable way to determine the size of the power corrections is to measure jet-veto cross sections at several different low \(p_T^{{\mathrm{veto}}}\) values and for different values of \(Q\) and compare it to the resummed perturbative prediction. It would be interesting to do so, and there should be enough Drell–Yan and \(Z\)-production data to make such a study possible.

3 Automated resummation

We now explain how to automate the resummation by suitably modifying existing fixed-order results. We shall employ two different resummation schemes. In Scheme A, we work with tree-level events obtained from MadGraph5_aMC@NLO [26]. We supply the beam functions from explicit calculations but compute the hard functions automatically and then reweight the events to achieve the resummation. In Scheme B, we use MadGraph5_aMC@NLO in fixed-order mode and compute the NLO cross section with a jet veto. To achieve the resummation we subtract the logarithmically enhanced pieces from the fixed-order cross section and multiply them back in resummed form. In this second scheme, both the hard functions and the beam functions are computed using MadGraph5_aMC@NLO. The second scheme is more convenient for practical computations but limited to NNLL order, while the first scheme allows (in principle) for arbitrary accuracy of the resummation.

3.1 Scheme A: NNLL from reweighting born-level events

The fact that the resummed result (3) has Born-level kinematics in the limit \(p_T^{{\mathrm{veto}}}\rightarrow 0\) makes it possible to achieve the resummation of large logarithms by a simple reweighting procedure. If we use a tree-level event generator such as MadGraph, the resummation can be implemented by rescaling the event weights with the ratio of the resummed to the tree-level cross sections at each kinematic point. Specifically, we need to replace the PDFs \(\phi _i\) used in the leading-order (LO) result with the beam functions \(\bar{B}_i\), and we need to supply the hard matching correction and the resummation factor \(E_i\). For an incoming particle of flavor \(i,j \in \{q, \bar{q}, g\}\), the reweighting factor at NNLL order reads

$$\begin{aligned} d\sigma _{ij}^{\text {NNLL}}(p_T^{{\mathrm{veto}}})= & {} \left( 1+\frac{\alpha _s(\mu _h)}{4\pi }\mathcal {H}_{ij}^{(1)}(Q^2,\hat{t},\mu _h)\right) \,\nonumber \\&\times E_{i}(Q^2,p_T^{{\mathrm{veto}}},\mu _h,\mu , R)\frac{\bar{B}_i(\xi _1,p_T^{{\mathrm{veto}}},\mu )}{\phi _i(\xi _1,\mu _{\text {Mad}})}\nonumber \\&\times \frac{\bar{B}_{j}(\xi _2,p_T^{{\mathrm{veto}}},\mu )}{\phi _j(\xi _2,\mu _{\text {Mad}})} \left( \frac{\alpha _s(\mu )}{\alpha _s(\mu _{\text {Mad}})}\right) ^N d\sigma ^{0}_{ij}(\mu _{\text {Mad}}).\nonumber \\ \end{aligned}$$
(11)

All the kinematic variables are determined by the event kinematics. At leading order \(\xi _1\) and \(\xi _2\) are just the momentum fractions of the incoming particles, \(\xi _i=2\,E_i/\sqrt{s}\). Note that we do not need to adopt the same value of the renormalization scale \(\mu \) as in the Born-level events, which were evaluated at a scale \(\mu _{\text {Mad}}\) inherent to the MadGraph code. However, in cases such as Higgs production, where the Born-level cross section depends on \(\alpha _s\), we have to multiply by the appropriate power \(N\) of the ratio \(\alpha _s(\mu )/\alpha _s(\mu _{\text {Mad}})\), where \(N=2\) for gluon-induced processes. We therefore only run MadGraph once, with a fixed reference scale \(\mu _{\text {Mad}}\). Scale uncertainties can then be estimated by repeating the reweighting with different values of \(\mu \) and \(\mu _h\).

Let us now detail the numerical implementation of the reweighting factor, starting with the beam functions, which are defined in (4) in terms of convolutions of perturbative kernel functions with PDFs. At one-loop order, they are linear in the logarithm of \(p_T^{{\mathrm{veto}}}\), and hence

$$\begin{aligned} \bar{B}_i(\xi ,p_T^{{\mathrm{veto}}},\mu )= & {} \phi _i(\xi ,\mu )+ \frac{\alpha _s(\mu )}{4\pi }\nonumber \\&\times \left( b_i(\xi ,\mu ) + c_i(\xi ,\mu ) \ln \frac{\mu }{p_T^{{\mathrm{veto}}}} \right) . \end{aligned}$$
(12)

To perform the reweighting in an efficient way, we compute and tabulate the convolution integrals for \(b_i(\xi ,\mu )\) and \(c_i(\xi ,\mu )\) for a grid of \(\xi \) and \(\mu \) values. Since the beam functions are independent of the final state, this can be done once and for all. Using the same grid as the underlying PDFs itself, we then use standard PDF interpolation routines to have fast and accurate numerical representations for the beam functions. We have implemented the beam functions and the resummation factor \(E_i(Q^2,p_T^{{\mathrm{veto}}},\mu _h,\mu , R)\) in a small Fortran code, which is called by the event reweighting routine written in Python.

The most complicated component of the reweighting factor by far is the hard function \(\mathcal {H}_{ij}^{(1)}(Q^2,\hat{t},\mu _h)\). This is process dependent and its computation requires a one-loop calculation. Fortunately, the necessary one-loop computations have been automated in the past few years. In particular, the MadGraph5_aMC@NLO framework provides the possibility to evaluate virtual corrections at specific phase-space points [33]. We use this code to evaluate the virtual corrections \(V_{ij}\) for each event. At each phase-space point, the code provides the result in the form of the coefficients \(C_i\) of the double pole, single pole, and finite terms in the expansion in \(\epsilon \), which is written in the form

$$\begin{aligned} V_{ij}= & {} d\sigma _{ij}^0(\mu )\left[ 1 + \frac{\alpha _s(\mu )}{4\pi } \frac{2\,e^{-\epsilon \gamma _E}}{\Gamma (1+\epsilon )} \left( \frac{\mu ^2}{\mu _\mathrm{Mad}^2}\right) ^\epsilon \right. \nonumber \\&\times \left. \left( \frac{C_2}{\epsilon ^2} + \frac{C_1(\mu _\mathrm{Mad})}{\epsilon } + C_0(\mu _\mathrm{Mad}) \right) _{\!ij} \right] . \end{aligned}$$
(13)

The scale \(\mu _\mathrm{Mad}\) can be chosen when running the MadLoop code. In the \(q\bar{q}\) channel, the double-pole coefficient is \(C_2=-C_F\gamma ^\mathrm{cusp}_0/2=-2C_F\), while the other two coefficients depend on the choice of \(\mu _\mathrm{Mad}\). For \(\mu _\mathrm{Mad}=Q\), the coefficient of the single-pole term is \(C_1(Q)=\gamma ^q_0=-3C_F\), and the finite part in the expansion of the above expression in \(\epsilon \) directly yields the hard function

$$\begin{aligned} \mathcal {H}_{q\bar{q}}^{(1)}(Q^2,\hat{t},\mu )= & {} 2 C_0(Q)+ C_F \left( \frac{\pi ^2}{3}-2 \ln ^2\frac{Q^2}{\mu ^2}+ 6 \ln \frac{Q^2}{\mu ^2} \right) .\nonumber \\ \end{aligned}$$
(14)

For \(Z\)-boson production one has \(C_0(Q)=-32/3+4\pi ^2/3\). For other choices \(\mu _\mathrm{Mad}\ne Q\) this result gets modified to

$$\begin{aligned} \mathcal {H}_{q\bar{q}}^{(1)}(Q^2,\hat{t},\mu )= & {} 2 C_0(\mu _\mathrm{Mad})+ C_F \left[ \frac{\pi ^2}{3}+2\ln ^2\frac{\mu _\mathrm{Mad}^2}{\mu ^2} \right. \nonumber \\&+\left. \ln \frac{\mu _\mathrm{Mad}^2}{\mu ^2} \left( 6-4\ln \frac{Q^2}{\mu ^2}\right) \right] . \end{aligned}$$
(15)

In practice, we first compute the hard function at some value of the reference scale \(\mu _\mathrm{Mad}\) for each event and write the result in the event record. The result at a different scale can then be obtained using the above relation. The reweighting script uses the result for the hard function and combines it with the beam functions and the resummation factor.

To obtain the best possible prediction, we match our result to the NLO fixed-order result for the cross-section. This matching allows us to also include terms which are power suppressed as \(p_T^{{\mathrm{veto}}}\rightarrow 0\). The simplest way to achieve the matching is to subtract from the resummed result its expansion to NLO and to then add back the full NLO result

$$\begin{aligned} \frac{d\sigma ^\mathrm{NNLL+NLO}}{dp_T^{{\mathrm{veto}}}}= & {} \frac{d\sigma ^\mathrm{NNLL}}{dp_T^{{\mathrm{veto}}}}- \left. \frac{d\sigma ^\mathrm{NNLL}}{dp_T^{{\mathrm{veto}}}}\right| _{\text {expanded to NLO}}\nonumber \\&+ \frac{d\sigma ^\mathrm{NLO}}{dp_T^{{\mathrm{veto}}}} . \end{aligned}$$
(16)

Our final NNLL \(+\) NLO result resums higher-order terms that are logarithmically enhanced, but also includes the full NLO result. To obtain the expansion of the resummed result, we simply do the reweighting with the fixed-order expansion of the reweighting factor in (11). The NLO result can be obtained from running MadGraph5_aMC@NLO in fixed-order mode. The difference between the full NLO result and the expansion of the resummed result is called the matching correction. By definition, this correction vanishes as \(p_T^{{\mathrm{veto}}}\rightarrow 0\) and is expected to scale as \(p_T^{{\mathrm{veto}}}/Q\). As we will discuss in Sect. 4.2, it is numerically very small for the values of \(p_T^{{\mathrm{veto}}}\) which are experimentally relevant.

3.2 Scheme B: NNLL \(+\) NLO with automated computation of the beam functions and matching corrections

In the reweighting scheme discussed above, we use MadGraph5_aMC@NLO to compute the hard functions but supply the beam functions from an explicit calculation. One can go even further and also compute the beam functions and the matching corrections automatically and in a single step. This is done by first factoring out the hard corrections and then performing a NLO run in the presence of the jet veto. An advantage of this second approach is that the beam functions are computed on the fly and it is therefore easy to use different PDF sets without any need to recompute the beam functions. A slight disadvantage is that one has to run MadGraph5_aMC@NLO in NLO mode. One can thus no longer work with events and will have to perform a new run when changing the cuts. However, if the matching is included in Scheme A described above, then a NLO run is needed also in this case. Note also that Scheme B only works at NNLL accuracy, while Scheme A allows for arbitrary precision if the necessary reweighting factor is supplied.

In order not to contaminate the matching corrections with the large logarithms contained in the hard function, we factor out the prefactor \(P_{ij}\) in (6) and define a reduced cross section \(\tilde{\sigma }_{ij}\) by

$$\begin{aligned} d\sigma _{ij}(p_T^{{\mathrm{veto}}}) = P_{ij}(Q^2,\hat{t},p_T^{{\mathrm{veto}}})\,d\tilde{\sigma }_{ij}(p_T^{{\mathrm{veto}}}) . \end{aligned}$$
(17)

The reduced cross section has the form

$$\begin{aligned} d\tilde{\sigma }_{ij}(p_T^{{\mathrm{veto}}}) = d\sigma _{ij}^0(Q^2,\hat{t},\mu )\, \bar{B}_i(\xi _1,p_T^{{\mathrm{veto}}})\,\bar{B}_j(\xi _2,p_T^{{\mathrm{veto}}}) + \Delta \tilde{\sigma },\nonumber \\ \end{aligned}$$
(18)

where \(\Delta \tilde{\sigma }=\mathcal{O}(p_T^{{\mathrm{veto}}}/Q)\) contains the power corrections and is given by the matching correction (16) divided by the prefactor. The function \(P_{ij}\) receives one-loop corrections from the hard function and the evolution factor \(E_{i}\) so that we can write

$$\begin{aligned} d\tilde{\sigma }_{ij}(p_T^{{\mathrm{veto}}})= & {} d\sigma _{ij}^{\text {NLO}}(p_T^{{\mathrm{veto}}},\mu )- \frac{\alpha _s(\mu )}{4\pi } \left( \mathcal {H}_{ij}^{(1)}(Q^2,\hat{t},\mu ) \right. \nonumber \\&+\left. E_{i}^{(1)}(Q^2,p_T^{{\mathrm{veto}}},\mu ) \right) d\sigma _{ij}^0(\mu ) . \end{aligned}$$
(19)

Provided we choose \(\mu \sim p_T^{{\mathrm{veto}}}\) in the reduced cross section \(\tilde{\sigma }\), all large logarithms are resummed in the RG-invariant prefactor \(P_{ij}\). Multiplying back the prefactor then yields the full NNLL \(+\) NLO cross section in the form

$$\begin{aligned}&d\sigma _{ij}^{\text {NNLL+NLO}}(p_T^{{\mathrm{veto}}})\nonumber \\&\quad = P_{ij}(Q^2,\hat{t},p_T^{{\mathrm{veto}}})\times d\tilde{\sigma }_{ij}(p_T^{{\mathrm{veto}}})\nonumber \\&\quad = \left( 1+\frac{\alpha _s(\mu _h)}{4\pi }\mathcal {H}_{ij}^{(1)}(Q^2,\hat{t},\mu _h) \right) E_{i}(Q^2,p_T^{{\mathrm{veto}}},\mu _h,\mu , R) \nonumber \\&\quad \quad \times \bigg [ d\sigma _{ij}^{\text {NLO}}(p_T^{{\mathrm{veto}}},\mu ) - \frac{\alpha _s(\mu )}{4\pi }\left( \mathcal {H}_{ij}^{(1)}(Q^2,\hat{t},\mu )\right. \nonumber \\&\qquad \quad \left. +\, E_{i}^{(1)}(Q^2,p_T^{{\mathrm{veto}}},\mu )\right) d\sigma _{ij}^{0}(\mu ) \bigg ]. \end{aligned}$$
(20)

Note that the matching procedure differs from the other scheme. In (16) above, we performed a purely additive matching, while in (20) the resummation factor \(E_i\) appears as an overall factor. This multiplicative matching generates higher-order logarithmic terms also for the power-suppressed contributions of order \(p_T^{{\mathrm{veto}}}/Q\) and higher. These additional terms are not controlled by the factorization theorem (3), which holds only at leading power, but one can hope that at least some of the logarithmic terms at subleading power are universal and will be captured by this treatment. For the case of Higgs production, the multiplicative matching scheme is preferred, since the perturbative corrections to the hard function are very large. In (20) they are extracted as a overall factor. For the \(q\bar{q}\)-initiated processes we study in this paper, the two schemes give almost indistinguishable results, as we will see in Sect. 4.2 below.

To implement (20) in MadGraph5_aMC@NLO we have directly modified its Fortran code by including the logarithmically enhanced terms. The expanded logarithmically enhanced terms, i.e. the second term on the right-hand side of (19), is similar to the compensating Sudakov factor introduced in the FxFx merging prescription, see (2.46) of [34], and it is therefore implemented at the same place in the code. In MadGraph5_aMC@NLO each real-emission phase-space configuration has corresponding Born kinematics defined by the FKS mapping [35]. Therefore we can always compute the prefactor \(P_{ij}\) using Born kinematics, and it can multiply the complete reduced cross section, including the real-emission contributions. In order to improve the run time, the time-consuming one-loop matrix elements are computed only once for each phase-space configuration, cached in memory, and used also for the (expanded) hard function. However, compared to normal running of MadGraph5_aMC@NLO, we cannot reduce the number of calls to the virtual corrections by using suitable approximations of it, as described in Sect. 2.4.3 of [26], because the reduced cross section is multiplied by them, resulting in positive feedback loops in setting up the approximations. When running MadGraph5_aMC@NLO in fNLO mode, setting the parameter ickkw in the run_card.dat to -1 turns on in the inclusion of the logarithmically enhanced terms and sets the hard and soft scales to \(Q\) and \(p_T^{{\mathrm{veto}}}\) (given by the ptj parameter in the run_card.dat), respectively. Hard and soft scale variations, as well as PDF uncertainties, can be computed at minimal CPU costs by reweighting [36]. This addition to the MadGraph5_aMC@NLO will become public with the next release of the code.

4 Phenomenological results

We now proceed to give numerical results for different electroweak-boson-production cross sections. Before presenting our final results, we discuss a variety of issues such as the proper choice of matching and factorization scales, the size of the matching corrections and the difference between the two resummation schemes discussed in the previous section. We then present results for the \(W^+W^-\) cross section as well as the cross section including the decay of the \(W\) bosons with cuts on the final-state leptons. Since the published measurements [14, 15] were taken at \(\sqrt{s} = 7\,\mathrm{TeV}\), we will present our results for this center-of-mass energy. For the electroweak parameters we use MadGraph5 default values, in particular \(\alpha _\mathrm{em}=1/132.5\), \(G_F=1.166\times 10^{-5}\,\mathrm{GeV}^{-2}\), \(M_W=80.42\,\mathrm{GeV}\), and \(M_Z=91.19\,\mathrm{GeV}\).

In all of our results below, we work with the MSTW2008 NNLO PDF set and its associated value \(\alpha _s(M_Z)=0.1171\) [37]. The choice of a NNLO PDF set seems appropriate, because we believe that the resummation captures the most important part of the NNLO corrections. In order to illustrate the size of the higher-order terms captured by resummation, we will also evaluate the NLO corrections using the NNLO PDF set, which increases the NLO prediction at \(p_T^{{\mathrm{veto}}}= 20\, \mathrm{GeV}\) by about 2 % in the case of \(W^+W^-\) production. We will define our jets using the anti-\(k_T\) algorithm with a jet radius of \(R=0.4\). The only quantity sensitive to the jet radius at NNLL \(+\) NLO accuracy is the anomaly exponent \(F_{ij}(p_T^{{\mathrm{veto}}},\mu ,R)\), and it is the same for all \(k_T\)-style clustering algorithms. As the default scheme for our plots we use Scheme A, since it is easier to disentangle and discuss the individual ingredients of the calculation (NLL versus NNLL resummation, matching to fixed-order perturbation theory) in this scheme. However, we find that both schemes give almost indistinguishable numerical results at NNLL \(+\) NLO level.

Fig. 2
figure 2

Resummed cross sections for \(Z\)-boson production (top) and \(W^+ W^-\) pair production (bottom) obtained at NLL (red) and NNLL (blue) order. The bands are obtained by varying the hard matching scale \(\mu _h\) and the factorization scale \(\mu \) by factors of 2 about their default values \(|\mu _h|=Q\) and \(\mu =p_T^{{\mathrm{veto}}}\). The gray bands show the fixed-order NLO results with scale variation \(\mu _r=\mu _f\in [p_T^{{\mathrm{veto}}}/2,\,2Q]\) for comparison. The panels on the left refer to the standard choice \(\mu _h^2>0\), while those on the right show results obtained using \(\mu _h^2<0\)

Full size image

4.1 Resummed results and choice of the hard scale

In Fig. 2 we show the results for the resummed \(Z\)-boson and \(W^+ W^-\)-pair-production cross sections at \(\sqrt{s}=7\) TeV, obtained with \(n_f=5\) light quark flavors and jet radius parameter \(R=0.4\). Here and below the two scales \(\mu \) and \(\mu _h\) are varied independently by factors of 2 about their default values \(\mu =p_T^{{\mathrm{veto}}}\) and \(\mu _h=Q\), where \(Q\) is the invariant mass of the electroweak final state, i.e. \(Q^2=M_Z^2\) for \(Z\)-boson production and \(Q^2=(q_{1}+q_{2})^2\) for the \(W^+ W^-\) final state (defined on an event-by-event basis). The resulting uncertainties are then added quadratically. In addition to the standard scale choice \(\mu _h^2\approx Q^2\) we consider using an imaginary value for the hard matching scale, such that \(\mu _h^2\approx -Q^2\). The corresponding results are shown on the right-hand side of Fig. 2. For comparison, we also show the NLO fixed-order results, which will be discussed in more detail in the next section. In all cases, we observe that going from NLL to NNLL accuracy improves the stability of the predictions significantly. Also, the NNLL bands are closer to the fixed-order NLO results than the NLL bands.

The use of an imaginary value of the hard matching scale \(\mu _h\) has been advocated in the context of Higgs production, because it maps the relevant hard function onto the space-like gluon form factor [38, 39]. This Euclidean quantity shows a much better perturbative behavior than the time-like form factor, which suffers from large numerical corrections \(\sim \) \((\alpha _s\pi ^2)^n\) due to imaginary parts from Sudakov double logarithms, which arise in time-like kinematics. The same arguments apply to the case of \(Z\) production. In [23], the choice \(\mu _h^2<0\) was applied to \(W^+ W^-\) pair production, and it was argued that this leads to a significant enhancement of the cross section, bringing the theoretical prediction in agreement with LHC measurements. Indeed, one can observe from Fig. 2 that the resummed results for the cross sections obtained with \(\mu _h^2<0\) are significantly larger than those obtained with the standard choice \(\mu _h^2>0\). For \(W^+ W^-\) production with \(p_T^{{\mathrm{veto}}}=25\) GeV, the increase in the central value of the NNLL \(+\) NLO cross section is about 4.8 % (which is of the same order as the recently calculated NNLO corrections [21]).Footnote 2 We stress, however, that in the case of multi-particle final states such as \(W^+W^-\) the hard function depends on several kinematic scales (\(\hat{s}\) and \(\hat{t}\) in the present case), some of which are time-like and some of which are space-like. Unfortunately, it is impossible to adopt a suitable choice of the hard matching scale, which would map the hard function onto a Euclidean quantity, such that all \((\alpha _s\pi ^2)^n\) terms can be resummed by means of RG evolution equations. It is therefore not clear whether the convergence of the perturbation series can be improved by using the choice \(\mu _h^2<0\). This problem was discussed in detail in the context of Higgs plus jet production in [40]. Even though the convergence in the right panels of Fig. 2 looks somewhat better than in the plots shown on the left, we have decided to adopt the conventional prescription \(\mu _h>0\) for the hard matching scale. Perhaps a more conservative way to assess the scale uncertainty would be to allow for arbitrary complex scale choices \(Q/2<|\mu _h|<2Q\) and then give the resulting uncertainty, as was recently proposed in [41].

For Higgs production, the resummation of jet-veto logarithms was performed to higher accuracy by including the two-loop hard and beam functions as well as the RG evolution factor at approximate N\(^3\)LL order [9]. The only missing ingredients for full N\(^3\)LL \(+\) NNLO accuracy are the three-loop anomaly exponent and the four-loop cusp anomalous dimension, whose effects have been estimated and included in the error budget. It was observed in this reference that the two-loop beam functions decrease the cross section, and we expect a similar effect in the present case. In the future, it should be possible to reach the same level of accuracy also for \(W^+ W^-\) production and related processes. The corresponding two-loop hard functions can be extracted from the two-loop virtual corrections, which have recently been obtained in [21, 42]. The product of beam functions integrated over rapidity could be extracted numerically from NNLO fixed-order codes for \(Z\)-boson production such as [43, 44], following the procedure employed in [9]. This is sufficient to obtain the inclusive \(W^+ W^-\) cross section, while a two-loop computation of the beam functions would be required for more exclusive cross-section predictions. Once (approximate) N\(^3\)LL+NNLO predictions for the \(W^+ W^-\) cross sections are available, the above-mentioned ambiguities, related to the choice of the hard matching scale, will be reduced significantly.

Fig. 3
figure 3

Left NLO predictions for the \(W^+ W^-\) production cross section obtained with a conservative estimate of scale uncertainties (gray), and with scale variations about high (green) and low (magenta) default values; see text for further information. Right kinematic distribution in the variable \(Q\) of the leading-order cross section

Full size image

4.2 Fixed-order results and matching

In order to obtained the best possible predictions, we need to match our resummed results for the cross sections with fixed-order expressions at NLO. The scale dependence of the NLO expression for the for the \(W^+ W^-\)-production cross section at \(\sqrt{s}=7\) TeV is shown in the left panel in Fig. 3. We set the factorization and renormalization scales equal (\(\mu =\mu _r=\mu _f\)) and vary them from \(\mu =p_T^{{\mathrm{veto}}}/2\) up to \(\mu =2Q\). This is a much larger scale variation than is usually considered, but this wider range seems appropriate since the problem at hand involves physics at both scales. For comparison, we also show the bands one would obtain from a variation of \(\mu \) by a factor of 2 around either a high default value \(\mu =Q\) or a low default value \(\mu =p_T^{{\mathrm{veto}}}\). Our broad scale variation is obviously more conservative, since it covers both options. Nevertheless, fixed-order computations usually adopt the high scale \(\mu =Q\) as the default value, and from Fig. 2 it appears that such a choice indeed leads to smaller higher-order corrections. A similar behavior is found for all cases studied in this paper. The invariant-mass distribution of the \(W\)-boson pair is shown in the right panel of Fig. 3. Defining the average hard scale \(\tilde{Q}\) by the median value of this distribution, one obtains \(\tilde{Q}=222\,\mathrm{GeV}\). This value will be useful in our phenomenological discussion below.

As discussed earlier and shown in (16), in Scheme A this matching is purely additive, i.e.

$$\begin{aligned} \sigma _{\text {NNLL+NLO}}= & {} \sigma _{\text {NNLL}}(\mu ,\mu _h)+ \Big (\sigma _{\text {NLO}}(\mu _m)\nonumber \\&- \sigma _{\text {NNLL}}(\mu _m)\big |_{\text {expanded to NLO}} \Big ) . \end{aligned}$$
(21)

The expansion of the resummed result is obtained by performing the reweighting with the reweighting factor expanded to NLO. If the resummation is performed with NNLL accuracy (or higher), the matching correction inside the parentheses is power suppressed in \(p_T^{{\mathrm{veto}}}/Q\). Note that we are free to use a different scale \(\mu _m\) for the matching correction than for the resummed result, since the power corrections in \(p_T^{{\mathrm{veto}}}/Q\) must be separately scale invariant. To obtain our uncertainty bands, the scales \(\mu \), \(\mu _h\), and the matching scale \(\mu _m\) are all varied independently. We then add the resulting uncertainties quadratically. We choose the number of flavors for the resummed results as \(n_f=5\), but since MadGraph5_aMC@NLO cannot produce five-flavor NLO results for \(W^+W^-\) due to the presence of top-quark resonant contributions in the NLO corrections, we calculate the matching corrections with \(n_f=4\) light flavors.

Fig. 4
figure 4

Resummed and matched predictions for the \(W^+ W^-\) production cross section (obtained by varying the matching scale about the default value \(\mu _m=p_T^{{\mathrm{veto}}}\) and \(\mu _m=Q\)) compared with the fixed-order result at NLO. The panels below the plots indicate the relative size of the power-suppressed matching corrections at NNLL order

Full size image

While the appropriate scale choice is clear for the case of the beam functions which describe emissions near the scale \(p_T^{{\mathrm{veto}}}\), the correct choice of \(\mu _m\) is not immediately obvious, because the matching corrections receive contributions associated with both the low and the high scale. The result for the cross section obtained with a high and a low matching scale is shown in Fig. 4, along with the corresponding relative size of the NNLL matching corrections. The matching corrections are well behaved in both cases. They are very small at the low \(p_T^{{\mathrm{veto}}}\) values shown in Fig. 4 and are therefore difficult to extract numerically. At larger values of \(p_T^{{\mathrm{veto}}}\) they grow linearly up to 3 % at \(p_T^{{\mathrm{veto}}}=80\) GeV. At NNLL order, the matching corrections are small enough that they could be safely ignored for values up to \(p_T^{{\mathrm{veto}}}=35\) GeV. At NLL order, on the other hand, not all leading-power NLO contributions are included in the resummed result, and therefore the predictions depend strongly on the matching scale \(\mu _m\). Figure 2 shows that the NNLL results lie rather close to the NLO results at the high scale \(\mu =Q\). Since, as we have pointed out above, the fixed-order perturbative expansion appears to work better with a high scale choice, we adopt \(\mu _m=Q\) as our default matching scale for all later predictions.

Fig. 5
figure 5

Left comparison of the resummed and matched NNLL \(+\) NLO predictions for the \(W^+ W^-\) cross section obtained in Scheme A (additive matching) with Scheme B (multiplicative matching). Right comparison of the NNLL \(+\) NLO predictions with the NLO result matched to Pythia using aMC@NLO

Full size image

In Scheme B, we do not have the freedom to choose the matching scale separately, since the matching corrections are not separated out, see (20). Numerically, we find that the results of Scheme A and Scheme B are almost indistinguishable, as can be seen in the left panel of Fig. 5. In the right panel of the same figure we show a comparison between our NNLL \(+\) NLO prediction for the \(W^+ W^-\) cross section and the result obtained after combining the NLO prediction with a parton shower using the MC@NLO prescription [45]. We observe that the latter prediction is lower than our result, in particular at higher values pf \(p_T^{{\mathrm{veto}}}\). This is astonishing at first sight, since one would expect that showering does not affect the cross section at higher \(p_T^{{\mathrm{veto}}}\) values. However, because the shower is unitary any change of the cross section at low transverse momenta must be accompanied by a compensating change at higher transverse momenta. Looking at the cross section as a function of the \(p_T\) of the leading jet, we find that the showered NLO result is higher than pure NLO result for all \(p_T^j>20\,\mathrm{GeV}\), so that the integral of the cross section for \(p_T>20\,\mathrm{GeV}\) is larger than the fixed-order result. After unitarization, this in turn implies that the jet-veto cross section, which is the integral \(0\le p_T\le 20\,\mathrm{GeV}\), is lower than the fixed-order result. The use of a matched parton shower therefore underestimates the jet-veto cross section. In contrast, we find that our NNLL \(+\) NLO resummed prediction lies closer to the fixed-order result indicated by the gray band. Genuine resummation effects are small as long as the fixed-order result for the cross section is computed with a high value \(\mu \sim Q\) of the renormalization scale.

Fig. 6
figure 6

Resummed and matched predictions for the cross sections for \(Z\), \(W^+ W^-\), and \(W^+ W^- W^\pm \) production, compared with NLO fixed-order predictions. The lower panels show the ratio of the cross section to the default NLO value with scale choice \(\mu =Q\)

Full size image

4.3 Multiple bosons and cross-section ratios

We are now ready to present our final results for a couple of interesting production cross sections involving multiple electroweak gauge bosons. In Fig. 6, we show predictions for the \(Z\)-, \(W^+ W^-\)-, and \(W^+ W^- W^\pm \)-production cross sections at the LHC with \(\sqrt{s}=7\) TeV; it would be straightforward to rerun our code at different values of the center-of-mass energy. In each case, we present our resummed and matched predictions at NLL \(+\) NLO and NNLL \(+\) NLO accuracy and compare them with the fixed-order NLO prediction. Notice that the value of the cross section drops by about a factor \(10^3\) with each additional boson. The triple-boson production cross section is tiny, but it constitutes a background to Higgs production in association with a \(W^\pm \) and subsequent decay \(H\rightarrow W^+ W^-\). The fact that we can obtain predictions for three-boson final states without any additional effort nicely demonstrates the power of our automated resummation scheme.

We find that the scale uncertainties of our NNLL \(+\) NLO predictions for \(W^+ W^-\) and \(W^+ W^- W^\pm \) production are estimated to be of similar size, while we obtain a much smaller uncertainty for the case of \(Z\)-boson production. This small scale variation should perhaps be taken with a grain of salt. At larger \(p_T^{{\mathrm{veto}}}\) values, our resummed cross section becomes similar to the fixed-order result, and its scale variation is similar to the scale variation of the fixed-order cross section obtained by performing a correlated scale variation with \(\mu _r=\mu _f\). An independent variation of \(\mu _r\) and \(\mu _f\), which is standard practice in fixed-order computations, would give an uncertainty that is twice as large. On the other hand, we have checked that the known NNLO corrections for \(Z\)-boson production are indeed compatible with our small uncertainty band. It is also interesting to note that for \(W^+ W^-\) production the scale uncertainties of the fixed-order prediction obtained from correlated and independent variations of \(\mu _r\) and \(\mu _f\) are found to be of similar size.

We also observe that the scale uncertainties of the fixed-order NLO predictions at small \(p_T^{{\mathrm{veto}}}\) values strongly increase with the number of produced bosons. This is not surprising if we consider the relevant scale ratio \(\tilde{Q}/p_T^{{\mathrm{veto}}}\), which governs the size of Sudakov logarithms. Using the median value \(\tilde{Q}\) of the invariant-mass distribution to estimate the hard scale, we find \(\tilde{Q}=M_Z\) for \(Z\) production, \(\tilde{Q}\approx 2.8\,M_W\) for \(W^+ W^-\) production, and \(\tilde{Q}=5.7\,M_W\) for \(W^+ W^- W^\pm \) production. In all cases, the three-momenta at which the bosons are produced scale with the boson mass, but the average scale increases with the number of the produced bosons. Note that after the resummation of Sudakov logarithms has been performed, the width of the uncertainty bands is only weakly dependent on the veto scale.

The relative perturbative uncertainty of our NNLL \(+\) NLO prediction for the \(W^+ W^-\)-production cross section at \(p_T^{{\mathrm{veto}}}=25\) GeV is \(^{+3.9\,\%}_{-3.0\,\%}\). It was advocated in [46] that taking the ratio of the \(W^+ W^-\)- and \(Z\)-boson-production cross sections might be a good way to reduce the uncertainty in the prediction of the jet-veto cross sections. This proposal was adopted in the experimental analysis reported in [14]. We have thus studied this cross-section ratio in some detail. We find that the relative uncertainty in the cross-section ratio is \(^{+5.2\,\%}_{-2.8\,\%}\), which is even slightly larger than the uncertainty in the \(W^+ W^-\)-production cross section itself. This makes it clear that taking the cross-section ratio does not help reducing the perturbative uncertainties, the reason being that the scale uncertainties are much smaller for \(Z\)-boson production than for \(W^+ W^-\) production. Even though the beam functions are the same in both cases, the cross sections involve different hard functions and RG evolution factors, which spoils the cancelation. We will now explain how an improved relation between the two production channels can be obtained, which only suffers from very small theoretical uncertainties. In a first step, it is useful to consider the jet-veto efficiencies defined as \(\sigma (p_T^{{\mathrm{veto}}})/\sigma \) instead of the cross sections \(\sigma (p_T^{{\mathrm{veto}}})\) themselves, because then the virtual corrections encoded in the hard functions largely drop out (even though this cancelation cannot be exact, since the hard corrections do not factor out of the total cross section). The inclusive cross section \(\sigma \) is evaluated at the hard scale \(\mu _h\). We use the NLO (LO) cross section together with the NNLL (NLL) approximation of \(\sigma (p_T^{{\mathrm{veto}}})\). Our resummed predictions for the ratio of the jet-veto efficiencies for \(W^+ W^-\) and \(Z\) production are shown in the left plot in Fig. 7. By construction, the relative uncertainties from varying \(\mu \) in this ratio are the same as in the ratio of the veto cross sections. In order to obtain more accurate predictions one needs to ensure that the RG evolution factors cancel out in the ratio. This can be accomplished by considering the ratio of the \(W^+ W^-\) cross section to the \(Z^*\) production cross section with an off-shell \(Z^*\) boson with invariant mass squared \(q^2=\tilde{Q}^2\), where \(\tilde{Q}\approx 222\,\mathrm{GeV}\) is the median of the invariant-mass distribution for the \(W^+ W^-\) final state shown in Fig. 3. The corresponding ratio of efficiencies is shown in the right plot in Fig. 7. It is close to 1 and exhibits very small scale uncertainties. A different way of relating \(Z\) and \(W^+ W^-\) production cross sections was proposed in [25]. These authors rescale the \(p_T^{{\mathrm{veto}}}\) value used in the \(W^+ W^-\) process by a factor \(M_Z/(2M_W)\) before relating it to the \(Z\)-boson production process. This rescaling is chosen such that the Sudakov logarithms have a similar size in the two cases. While [25] finds a nice agreement for the NLO efficiencies obtained using this rescaling prescription, it is clear that the relation cannot be exact, since QCD is not scale invariant. Furthermore, the agreement becomes worse if one rescales the \(p_T^{{\mathrm{veto}}}\) value with the more appropriate factor \(M_Z/\tilde{Q}\). In the middle plot of Fig. 7, we show the corresponding ratio of efficiencies, which suffers from sizable scale uncertainties.

Fig. 7
figure 7

Resummed predictions for the ratio of the jet-veto efficiencies for \(W^+ W^-\)- and \(Z\)-boson production (left). In the middle plot the \(p_T^{{\mathrm{veto}}}\) value of the \(W^+ W^-\) process is rescaled by a factor \(M_Z/(2M_W)\), as proposed in [25]. The right plot shows the same ratio for \(W^+ W^-\)- and \(Z^*\)-boson production, where the off-shell boson has invariant mass \(\tilde{Q}_{WW}=222\,\mathrm{GeV}\). The bands are obtained by varying the low scale \(\mu \) about its default value \(\mu =p_T^{{\mathrm{veto}}}\), while keeping the hard matching scale \(\mu _h\) fixed

Full size image

4.4 Experimental cuts

An important advantage of our framework is that we can include the decay of electroweak bosons, together with cuts on the leptonic final state. In the experimental measurements of \(W^+W^-\) production, candidate events are selected with two opposite-sign charged leptons, electrons or muons, and missing transverse momentum coming from the neutrinos in \(pp\rightarrow W^+ W^- +X\rightarrow l\,\nu \,l'\nu '+X\). To account for the detector geometry and to suppress the background from Drell–Yan and top production, a number of cuts are applied to the final state in addition to the jet veto. For example, the ATLAS analysis [14] imposes the following cuts in the \(e^+e^-\) channel:

  1. 1.

    lepton \(p_T>20\,\mathrm{GeV}\),

  2. 2.

    leading lepton \(p_T>25\,\mathrm{GeV}\),

  3. 3.

    lepton pseudorapidity \(\eta _e<1.37\) or \(1.52<\eta _e<2.47\),

  4. 4.

    dilepton invariant mass \(m_{e^+e^-}>15\,\mathrm{GeV}\) and also \(|m_{e^+e^-}-m_Z|>15\,\mathrm{GeV}\).

The cuts applied in the \(\mu ^+\mu ^-\) channel are fairly similar, while those on the mixed final states \(e^\pm \mu ^\mp \) are looser, because they have much smaller Drell–Yan background. In Fig. 8, we show the cross section for the production and decay \(pp\rightarrow W^+ W^- +X\rightarrow e^+ e^-\nu \bar{\nu }+X\) in the presence of these cuts as a function of the jet-veto scale. The experimental analysis in [14] uses the anti-\(k_T\) algorithm with \(R=0.4\) and fixed \(p_T^\text {veto}=25\,\mathrm{GeV}\). Comparing this figure with the lower plots in Fig. 2, we see that the uncertainties of the cross section are similar to the inclusive case and that the matching corrections remain small also in the presence of the cuts.

Fig. 8
figure 8

Resummed and matched predictions for the \(pp\rightarrow W^+ W^- +X\rightarrow e^+ e^-\nu \bar{\nu }+X\) cross section with the cuts on the leptonic final state described in the text

Full size image

The experimental analysis [14] imposes a few additional cuts, in particular a minimum total transverse momentum of the two charged leptons \(p_T^{e^+e^-}\,>30 \, \mathrm{GeV }\) and minimum requirements on the missing transverse momentum \(p_{T, \text {Rel}}^{\nu \bar{\nu }}\,>45\, \mathrm{GeV}\).Footnote 3 The cut \(p_T^{e^+e^-}\,>30 \, \mathrm{GeV }\) is somewhat problematic for the theoretical analysis, especially when it is applied to predict the \(Z\)-boson background to \(W^+ W^-\) production. The difficulty is that we must make sure that the leptonic cuts do not (strongly) affect the hadronic final state. In the case of \(Z\) production the \(p_T^{e^+e^-}\) is equal (and opposite) to the transverse momentum \(p_T^X\) of the hadronic final state. Imposing a lower bound on \(p_T^{e^+e^-}\) is the same as imposing a lower bound on \(p_T^X\). This interferes with the jet-veto cut which at NLO corresponds to an upper cut on \(p_T^X\). The factorization formula in [8] does allow for additional cuts on \(p_T^X\) in the presence of the jet veto, but the relevant beam functions would be more complicated than those needed without such cuts. For the \(W^+ W^-\)-production process, the quantity \(p_T^{e^+e^-}\) is not directly related to \(p_T^X\) because of the presence of the neutrinos, but the corresponding cut still affects the low-\(p_T^X\) region.

4.5 Difficulties associated with photons

Our framework cannot immediately be applied to processes involving photons. The reason is that photons are massless particles and have hadronic substructure. At high energies, a photon thus needs to be treated as a photon jet, or more precisely a photon surrounded by some hadronic radiation. In fact, many photon-isolation requirements necessitate fragmentation functions. This can be avoided using the photon isolation proposed by Frixione [47], but also in this case the photon has a partonic content and a proper description needs to take into account partons emitted collinear to the photon. This implies that our factorization theorem does not apply, since it assumes that all energetic radiation is collinear to the beam. The photon isolation introduces new small scales to the problem (e.g. the hadronic energy around the photon), which give rise to additional large logarithms not associated with the jet veto.

It is nevertheless interesting to see what happens when we apply our resummation scheme to a process involving photons. To this end, we consider \(W^\pm \gamma \) production using the same setup as before (\(\sqrt{s}=7\) TeV, \(R=0.4\), \(n_f=4\)) and imposing the isolation requirement proposed in [47], with associated parameters \(R_0^\gamma =0.4\), \(x_n=1.0\), and \(\epsilon _\gamma =1.0\). The corresponding results are shown in Fig. 9. The \(pp\rightarrow W\gamma \) process suffers from very large NLO corrections (the LO results are similar to the NLL result). The resummed results, on the other hand, are not very different from the LO predictions, so that the matching corrections are huge, indicating that there are indeed other sources of large corrections in this process. Likely these arise due to Sudakov effects associated with photon isolation. However, even the logarithms associated with the jet veto have a more complicated structure once a process involves partons collinear to the photon directions, which becomes possible at NLO. It would be interesting to analyze such photon processes in the context of SCET. In its present implementation our method does not resum all large corrections in these cases.

Fig. 9
figure 9

Theoretical predictions for \(W\gamma \) production obtained from our resummation scheme. The left plot shows the resummed results without matching to NLO, while the right plot shows the results obtained after the matching has been performed. A proper treatment of production processes with high-energy photons in the final state would require a generalization of the factorization formula (3)

Full size image

5 Conclusion

Higher-order logarithmic resummations in collider physics, both in SCET and using traditional methods, are typically done on a case-by-case basis, similar to the way fixed-order calculations were performed a few years ago. In the meantime, several groups have automated NLO computations in a variety of computer codes. This automation saves time, reduces the possibilities for mistakes and offers the flexibility to also study effects beyond the Standard Model. It is desirable to have the same level of automation for higher-order resummations of large logarithmic corrections. In the present paper, we have achieved this goal for electroweak-boson-production cross sections in the presence of a jet veto, at NNLL \(+\) NLO accuracy. This combination is natural because in the Sudakov region, where \(\ln (Q/p_T^{{\mathrm{veto}}})\sim 1/\alpha _s\), NNLL logarithmic terms have the same parametric scaling as NLO corrections in a region where there are no large logarithms. In contrast, taking resummation effects into account using a parton shower gives a lower parametric accuracy, and the unitarization inherent in the shower approach can sometimes be problematic. In the case of the jet-veto cross section for \(W^+ W^-\) production, for example, unitarization leads to cross sections that are systematically lower than the NNLL \(+\) NLO results.

Resummations are relevant in kinematical configurations which are close to the Born-level kinematics and can therefore be obtained by reweighting Born-level cross sections with appropriate factors. The most complicated ingredient for NNLL resummations are the one-loop hard functions, which encode the virtual corrections. Their computation has been automated, and we use the MadGraph5_aMC@NLO framework to obtain the hard function required for our analysis. We have also presented a modified scheme, in which the beam functions accounting for collinear emissions and the matching onto fixed-order results is automated and performed using existing fixed-order codes. This is possible, because the hard function and the resummation of large logarithms are just overall factors in the differential cross section.

We have used our method to perform a detailed analysis of resummation effects for the \(W^+ W^-\)-pair-production cross section, for which experimental measurements found a slight excess compared to theoretical predictions based on NLO computations matched to parton showers. We observe that the NLO result with a high value of the renormalization and factorization scales \(\mu _r \sim \mu _f \sim Q\) is in good agreement with the NNLL \(+\) NLO resummed predictions, while the results obtained with a matched parton shower are systematically lower. This effect, together with the positive NNLO corrections to the total rate which are now known, helps to bring the Standard Model prediction into better agreement with the measurements. It would be important to include the two-loop virtual corrections into the resummation and to also compute and include two-loop beam functions. This improvement, which is beyond the scope of the present work, would lead to very precise predictions, which could be directly compared with the experimental results. This level of accuracy has already been achieved in Higgs production by extracting the beam functions numerically. It was found that the two-loop corrections to the beam functions were sizable, because they are enhanced by logarithms of the jet radius. In Higgs production, the NNLO corrections to the hard function increase the cross section, while the two-loop beam functions lower it. We expect the same behavior for the \(W^+ W^-\) case, and it will be interesting to see the combined effect of these improvements on the final predictions. Also, at NNLO the \(gg\) channel starts to contribute to \(W^+ W^-\) production and could give rise to important corrections. Since this channel has already been implemented into the MadGraph5_aMC@NLO framework, it will be straightforward to perform the corresponding resummation using our method.

It would also be interesting to generalize our methods to processes with jets in the final state. In addition to hard, beam, and soft functions, these processes involve jet functions describing the energetic final-state radiation. Furthermore, the hard function then has a non-trivial color structure. Existing programs which compute virtual corrections for NLO processes currently only supply squared matrix elements summed over colors, but they can be modified to provide the color information needed for SCET-based resummation. This color structure is then contracted with the color structure of the soft function after RG evolution. The soft, beam, and jet functions will in general need a separate calculation. However, since the jet and beam functions are two-point functions and the soft function is given by a single emission from eikonal lines, these computations are much simpler than full-fledged real-emission computations and could be automated as well. We are confident that such automated resummations will become available in the future and provide higher-order logarithmic resummations for a much wider range of observables.