Drell-Yan Production at NNLL'+NNLO Matched to Parton Showers

We present results for Drell-Yan production from the GENEVA Monte-Carlo framework. We combine the fully-differential NNLO calculation with higher-order resummation in the 0-jettiness resolution variable. The resulting parton-level events are further combined with parton showering and hadronization provided by PYTHIA8. The 0-jettiness resummation is carried out to NNLL', which consistently incorporates all singular virtual and real NNLO corrections. It thus provides a natural perturbative connection between the NNLO calculation and the parton shower regime, including a systematic assessment of perturbative uncertainties. In this way, inclusive observables are correct to NNLO, up to small power corrections in the resolution cutoff. Furthermore, the perturbative accuracy of 0-jet-like resummation variables is significantly improved beyond the parton shower approximation. We provide comparisons with LHC measurements of Drell-Yan production at 7 TeV from ATLAS, CMS, and LHCb. As already observed in $e^+e^-$ collisions, for resummation-sensitive observables, the agreement with data is noticeably improved by using a lower value of $\alpha_s(M_Z) = 0.1135$.


I. INTRODUCTION
There are several different ways to obtain theoretical predictions for collider processes involving strongly interacting particles, namely fixed-order (FO) perturbation theory, resummed perturbation theory, and predictions using parton shower (PS) algorithms.
In FO perturbation theory, the perturbative expansion is carried out to a given order in the strong coupling constant α s . The leading order (LO) corresponds to the lowest order in α s required to obtain the desired process. The next-toleading order (NLO) then includes all terms of Oðα s Þ relative to the LO result, and so on.
When the process or observables in consideration involve large ratios of physical scales, FO perturbation theory can become unreliable or even break down, since for each order in α s one can encounter up to two powers of logarithms of large scale ratios. In this case, the most precise perturbative predictions are obtained by resumming the logarithmic terms to all orders in α s . At leading logarithmic (LL) order one includes the most dominant terms of order α n s ln 2n . The next-to-leading logarithmic (NLL) order also includes the next-largest subleading logarithmic terms, and so on. Like at FO, the resummed perturbative expansion happens in a systematic fashion, which is described in more detail later on.
PS algorithms formally work in the strongly ordered limit, where each subsequent emission happens at a smaller resolution scale than the previous one. They start from a given LO prediction and effectively multiply it by splitting probabilities for each additional emission. This allows in principle any observable to be predicted at LL (and including some NLL effects) using a probabilistic Markov-Chain algorithm.
Each type of prediction has its advantages as it provides the most accurate description in the limit where other predictions are not applicable. FO predictions are necessary for a precise description of additional hard (i.e. energetic wide-angle) emissions. These are not correctly described even at LO by either resummation or parton showers, since both of these approaches rely on a kinematic expansion in a small resolution variable. For observables sensitive to many soft and collinear emissions generating large logarithms, FO predictions are not suitable and resummation and parton showers are necessary. Resummed perturbation theory allows one to systematically carry out the resummation beyond the LL or strongly ordered limit and is therefore more accurate than parton showers. On the other hand, the higher-order resummation applies at the level of observables that are still sufficiently inclusive, i.e., that marginalize over many emissions and are only sensitive to a small number of physical scales. To obtain a prediction of the final state that is fully exclusive in all emissions, parton showers are required. In particular, parton shower predictions allow one to further attach a hadronization model to generate fully exclusive hadron-level events. These are an essential requirement for experiments to be able to simulate their data and study detector effects, and allow for the most direct comparison with experimental data.
Clearly, to obtain the best possible predictions, it is desirable to combine the different theoretical descriptions in such a way that one benefits from the advantages of each. In particular, this allows for reliable predictions also in the transition regions in between different parametric limits, which are often important in practice. The (often difficult to answer) question of which parametric regime is the most appropriate in a given phase-space region then becomes much less relevant for practical purposes.
Combining resummed and fixed-order perturbation theory is a standard procedure in higher-order resummed calculations, where it is well understood how to systematically match the pure resummed result to the fixed-order result in the appropriate fixed-order limit.
In the above methods, the primary goal is to improve the perturbative accuracy of parton shower Monte Carlo (MC) programs for inclusive FO observables, while the description of the resummation region is left to the parton shower. Formally, this combination amounts to matching the FO calculation to the LL resummation in the parton shower's resolution variable. The perturbative accuracy criteria that should be satisfied at NNLO þ PS have been discussed in detail in Ref. [33]. GENEVA goes beyond this by also incorporating the higher-order resummation in a suitable resolution variable that separates the FO and parton shower regimes. The perturbative matching essentially happens now from FO to resummation to parton shower. By carrying out the resummation to the appropriate order consistent with the FO accuracy, it effectively mediates between the FO calculation and the parton shower and bridges the perturbative gap between them. This provides several key benefits: (i) A systematic combination of NNLO calculations with parton showers, which extends to arbitrary processes. (ii) An improved perturbative accuracy of (sufficiently inclusive) resummation observables. (iii) An improved description of the matching/transition between resummation and fixed-order regimes. (iv) A systematic assessment of perturbative uncertainties in the resummation and matching. An implementation of our approach for e þ e − collisions was presented in Ref. [24], where NLO calculations for e þ e − → 2 jets and e þ e − → 3 jets, and the NNLL 0 resummation for thrust (2-jettiness) as resolution variable are combined with the parton showering and hadronization provided by PYTHIA8 [34,35]. For all the tested two-jet-like resummation-sensitive observables, the predictions from GENEVA closely match the corresponding exact NNLL ð0Þ resummed calculations, and after hadronization are in excellent agreement with LEP data. Although not emphasized there, the results in Ref. [24] already include all singular NNLO corrections, and are formally accurate to NNLO up to power corrections in the resolution cutoff.
In this paper, we present results for Drell-Yan production pp → γ=Z → l þ l − . Our implementation combines the general construction of Ref. [33] with the inclusion of higher-order resummation in the jet resolution variable similar to Ref. [24]. We combine the fully differential NNLO calculation for Drell Yan with the NNLL 0 resummation for 0-jettiness T 0 (also know as beam thrust [36]) with PYTHIA8. The FO accuracy of the generated events is NNLO for pp → l þ l − , NLO for pp → l þ l − j, and LO for pp → l þ l − jj (which we refer to as NNLO 0 , NLO 1 , and LO 2 , respectively) while the T 0 distribution is exactly reproduced at NNLL 0 . We show that by including the higher-order resummation for T 0 as the resolution variable, the perturbative accuracy for other zero-jet-like observables in the resummation region is significantly improved, and our predictions agree well with dedicated higher-order resummed calculations.
Comparing our results to Drell-Yan measurements from ATLAS, CMS, and LHCb, obtained during the 7 TeV run at the LHC, we obtain overall good agreement with the data for rapidity distributions, inclusive and exclusive jet cross sections, as well as the transversemomentum distribution of the hardest jet. For the transverse momentum of the vector boson and the related ϕ Ã distribution between the leptons, our default settings with α s ðM Z Þ ¼ 0.118 predict a somewhat harder spectrum. We observe a much better agreement with the data when using a lower value of the strong coupling constant α s ðM Z Þ ¼ 0.1135, as was obtained from fits of N 3 LL 0 resummed calculations of event shapes to LEP data [37][38][39][40]. The same was also observed for resummation observables in our e þ e − results in Ref. [24]. This paper is organized as follows. In Sec. II we discuss the theoretical framework of GENEVA. In Sec. II A we provide an overview of the general construction and the master formulas for the partonic jet cross sections that serve as the main perturbative input, and in Sec. II B we discuss specific choices and implementation aspects that are of relevance in more detail. In Sec. III, we describe how the partonic events with up to two partons are interfaced with a parton shower in a way that avoids double counting and such that the parton shower does not destroy the perturbative accuracy of the partonic calculation. In Sec. IV, we compare GENEVA with dedicated higher-order perturbative calculations, and in Sec. V we compare GENEVA's results with measurements from the LHC. We conclude in Sec. VI.

A. General setup
In this subsection we discuss the general construction we use. Our discussion here follows that of Ref. [33], and we refer to that paper for more details.
Drell-Yan production at NNLO receives contributions from partonic processes with up to two final-state partons. In a standard FO calculation, the phase space for these partonic contributions is integrated over in such a way that the cancellation of virtual and real IR divergences happens at the level of observables.
In contrast, an event generator is meant to produce physical events, which means that all IR divergences should cancel on a per-event basis. This implies that an N-parton event should correspond to an IR-finite partonic N-jet cross section, which is fully differential in the corresponding partonic N-jet phase space. That is, each generated event represents a point in an N-jet phase space, rather than an N-parton phase space, and each four momentum in the event provides the energy and direction of a partonic jet.
Therefore, the basis of the GENEVA MC framework is the formulation of the perturbative inputs in terms of "MC cross sections," which are well-defined partonic jet cross sections according to which the events are distributed, and which can be systematically computed to the desired perturbative accuracy in FO and resummed perturbation theory [24,41,42].
The MC cross sections are defined in terms of an N-jet resolution variable T N and formally include the contributions of an arbitrary number of unresolved emissions below a resolution cutoff T N < T cut N . In the present case, we require events with zero, one, and two partons, which are distributed according to The exclusive zero-jet MC cross section is defined by T 0 < T cut 0 , the exclusive one-jet MC cross section by T 0 > T cut 0 and T 1 < T cut 1 , and the inclusive two-jet MC cross section by T 0 > T cut 0 and T 1 > T cut 1 . In this way all of the partonic phase space is covered. Adding the one-jet and two-jet events, we can also define the inclusive one-jet MC cross section as which is defined by T 0 > T cut 0 and does not depend anymore on T cut 1 . In Eq. (2) we have made use of the shorthand notation Since the partonic jets are represented by on-shell quarks and gluons in the partonic Φ N events, the partonic jet definition used here is quite different from an ordinary jet algorithm. It depends on a phase-space map Φ N ðΦ M Þ that projects the M-body phase space of unresolved emissions onto the distributed Φ N points, where N ≤ M. This phasespace map must be IR safe, such that the MC cross sections are IR finite. For example, at fixed NNLO the exclusive zero-jet cross section is given by where B N contains the N-parton tree-level contributions, V N the N-parton one-loop contributions, and W 0 the twoloop contribution. There is considerable freedom in the precise definitions of T N and the Φ N ðΦ M Þ map, and we choose particular definitions based on their theoretical properties and simplicity.
Using the events in Eq. (1), the cross section for any observable X is given by where M X ðΦ N Þ is the measurement function that computes the observable X for the N-parton final state Φ N . This cross section is not identical to the exact fixed-order result, because for any unresolved emissions the observable is calculated on the projected phase-space point Φ N ðΦ M Þ, rather than the exact Φ M . This is a fundamental limitation inherent to generating IR-finite events beyond LO and for a more detailed discussion we refer to Ref. [33]. The difference vanishes in the limit T cut N → 0. Hence, we want to choose T cut N as small as possible to have a maximally exclusive description. For small T cut N , however, the MC cross sections contain large logarithms of T N and T cut N , requiring resummation to obtain physically meaningful predictions for them. This is precisely what a parton shower would do at LL (also including some NLL effects), in which case T cut N would play the role of the shower cutoff. In GENEVA we are improving over this by resumming these logarithms at higher accuracy.

0=1-jet separation
To discuss the construction of the MC cross sections, we first consider the separation between zero and one jets which is determined by the zero-jet resolution variable T 0 . We write the exclusive zero-jet and inclusive one-jet MC cross sections as [24,33] Here, dσ resum 0 is the spectrum differential in Φ 0 with the dependence on T cut 0 resummed to a given logarithmic accuracy, while the remaining two terms give the matching corrections required to reproduce the desired FO accuracy. Their precise form depends on the fixed-order content of the resummed contribution. The singular matching dσ sing match 0 contains all contributions that do not vanish as T cut 0 → 0, while the nonsingular matching dσ nons 0 only contains contributions that vanish in this limit. The analogous separation is done for dσ MC ≥1 , where dσ resum ≥1 resums the differential dependence on T 0 , the singular matching contains contributions that diverge at least like 1=T 0 for T 0 → 0, and the nonsingular matching contains contributions that contain at most integrable singularities for T 0 → 0.
We carry out the resummation of the T 0 dependence to NNLL 0 . Thus, the resummed contributions are given by Here, dσ NNLL 0 =dT 0 dΦ 0 is the resummed differential T 0 spectrum and dσ NNLL 0 =dΦ 0 ðT cut 0 Þ its cumulative integral; see Sec. II B 2 below. Note that the two-loop virtual corrections, which at FO are proportional to δðT 0 Þ, are fully incorporated at NNLL 0 . This means they are properly spread to nonzero values of T 0 and contribute to the differential spectrum as dictated by the resummation. The T 0 resummation is naturally differential in Φ 0 , so we can directly use the resummed cumulant in dσ resum 0 . To make the T 0 spectrum fully differential in Φ 1 , we have defined a normalized splitting probability PðΦ 1 Þ which satisfies Z dΦ 1 dΦ 0 dT 0 Since T 0 can be considered as part of the radiation phase space dΦ 1 =dΦ 0 , the integration here is effectively over two remaining radiation variables, e.g. an energy splitting and an azimuthal angle. Thus, the integral of the splitting probability over all Φ 1 points restricted to any fixed values of Φ 0 ðΦ 1 Þ and T 0 ðΦ 1 Þ is equal to 1. This will be discussed in more detail in Sec. II B 4. At NNLL 0 , the resummation fully incorporates all singular contributions to Oðα 2 s Þ, implying that the singular matching vanishes, dσ sing match At our desired NNLO accuracy, dσ MC 0 and dσ MC ≥1 must be correct to NNLO 0 and NLO 1 , respectively, which determines the nonsingular matching corrections to be dσ nons The terms in square brackets are the FO expansions to α 2 s of the resummed cumulant and spectrum in Eq. (7). The NNLO 0 result for the cumulant is given in Eq. (4). The NLO 1 result for the fully differential spectrum is given by It depends on the projection where Θ T ðΦ 2 Þ defines the region of Φ 2 that can be projected onto massless Φ 1 via the IR-safe phase-space map Φ T 1 ðΦ 2 Þ. Only this projectable region of Φ 2 is included in dσ nons ≥1 =dΦ 1 , while the remainder will be included in the Φ 2 events below.
The fact that the NNLL 0 resummation reproduces all singular terms of the full FO result to Oðα 2 s Þ implies that the resummed expanded result in dσ nons ≥1 =dΦ 1 acts precisely as a differential NNLO T 0 -subtraction [43] (see also Ref. [44]). The cancellation of the singular terms between the FO and resummed expanded results is not point by point in Φ 1 but nonlocal, which means it only holds "on average" upon integrating over dΦ 1 =dΦ 0 dT 0 . As discussed in detail in Ref. [43], the cancellation can be made more local by considering more-differential resummations (e.g. utilizing the results of Refs. [45][46][47]). For the cancellation to be point by point in T 0 , dσ NLO 1 ≥1 =dΦ 1 has to reproduce the right singular T 0 dependence when projected onto dT 0 dΦ 0 , which means the map Φ T 1 ðΦ 2 Þ has to preserve the value of T 0 , In other words, if this was not the case, we would not obtain the correct nonsingular corrections for the T 0 spectrum, and this would destroy its NNLL 0 accuracy. The underlying reason is that we perform the FO calculation in an event generator form as in Eq. (5), and so to compute the FO T 0 spectrum we always project some Φ 2 points, namely those below T 1 < T cut 1 , onto Φ 1 and only then onto T 0 . More details on this map are given in Sec. II B 6.

1=2-jet separation
Next, we consider the separation of the inclusive one-jet cross section into an exclusive one-jet MC cross section and an inclusive two-jet MC cross section, using the one-jet resolution variable T 1 . We write them as Here, dσ resum 1 and dσ resum

≥2
contain the resummation of the T cut 1 and differential T 1 dependence. The dσ match 1 and dσ match ≥2 contain the matching corrections required to achieve the desired FO accuracy. For the T 1 resummation we limit ourselves to LL, which is sufficient for matching to the parton shower. At this level, we can write the resummed contributions as dσ resum Here, U 1 ðΦ 1 ; T cut 1 Þ denotes the Sudakov factor that resums the dependence on T cut 1 to at least LL accuracy, and U 0 1 ðΦ 1 ; T 1 Þ denotes its derivative with respect to T cut 1 and resums the differential T 1 dependence. Its details are given in Sec. II B 5. The differential T 1 resummation in dσ ≥2 is evaluated at the projected Φ 1 point Φ T 1 ðΦ 2 Þ. As we see below, this is required so the T 0 resummation is evaluated at the correct T 0 value. The normalized splitting probability PðΦ 2 Þ is defined analogous to PðΦ 1 Þ in Eq. (8) but for Φ 1 → Φ 2 and using the Φ T 1 ðΦ 2 Þ map, with the details given in Sec. II B 4. The T 1 resummation acts on top of the overall dσ C ≥1 =dΦ 1 in Eq. (15), which is the inclusive one-jet cross section in the T 1 -singular limit. Its NLO 1 expansion therefore has to be dσ C ≥1 dΦ 1 (Its full form will include the T 0 resummation and is derived below.) Here, C 2 ðΦ 2 Þ denotes a standard NLO subtraction [48,49] that reproduces the pointwise singular behavior of B 2 ðΦ 2 Þ. We use the Frixione-Kunszt-Signer (FKS) subtractions [49] in our implementation. Note that the subtraction term here uses its own projection Using Eqs. (15)- (16) and requiring that dσ MC 1 and dσ MC ≥2 are correct to NLO 1 and LO 2 , respectively, then determines the matching corrections to be dσ match Here, U ð1Þð0Þ 1 ðΦ 1 ; T cut 1 Þ denotes the Oðα s Þ term in the expansion of the U ð0Þ 1 ðΦ 1 ; T cut 1 Þ Sudakov factor (or its derivative). These terms cancel the leading double-logarithmic terms in T 1 in the FO pieces (the B 2 contributions). As long as we only include the LL T 1 resummation, there will be some subleading single-logarithmic terms present in the matching corrections.
So far, we have discussed the construction of an additive NLO 1 þ LL 1 matching. The important point is now that this should also include the T 0 resummation which is important at small T 0 . The condition for this is that the integral of NLO 1 þ LL 1 result with the Φ T 1 ðΦ 2 Þ map must reproduce the T 0 -resummed result for the inclusive one-jet MC cross section in Eq. (6). That is, where we have used the identity Inserting the expression for dσ MC ≥1 in terms of dσ resum ≥1 and dσ nons ≥1 in Eqs. (7) and (10), we obtain the result for It contains the full differential T 0 resummation via dσ resum ≥1 . The difference to the full inclusive one-jet cross section is that the T 0 -nonsingular terms are now evaluated in the T 1singular approximation, and this mismatch is compensated by the T 1 matching contributions.
The above expressions completely define the fully differential jet cross sections. In the next subsection we provide some additional details on the specific implementation of the various pieces in GENEVA.

Choice of the jet resolution variables
We choose N-jettiness [50] as our N-jet resolution variable. It is defined as where the sum over k runs over all final-state particles, excluding the vector boson and all its decay products. We use a canonical geometric measure, whereq i ¼ n i ¼ ð1;ñ i Þ are lightlike reference vectors along the jet and beam directions. N-jettiness is a global event shape that is explicitly designed as an N-jet resolution variable, as it measures the degree to which the final state is N-jet-like for a given N and automatically clusters the final state into N-jet and beam regions in an IR-safe way and without any dependence on an additional jet clustering algorithm. Furthermore, it is theoretically well understood and its all-orders singular structure and resummation are known.
For Drell-Yan, we need 0-jettiness T 0 and 1-jettiness T 1 , and we always evaluate Eq. (22) in the longitudinally boosted frame where the vector boson has zero rapidity. The reference vectors for the beam directions are always along the beam directions, soñ a ¼ẑ andñ b ¼ −ẑ. For T 1 , we also need to define a jet directionñ J . For our calculation we only need to define it for a Φ 2 event, in which case we can chooseñ J by minimizing T 1 . We define a clustering metric and find the smallest of dðp 1 Þ, dðp 2 Þ, and dðp 1 ; p 2 Þ, where p 1 and p 2 are the two parton momenta in the Φ 2 event. If dðp i Þ is smallest, we can think of p i being clustered with one of the beams, and thenñ J ¼p j =jp j j with j ≠ i. If dðp 1 ; p 2 Þ is the smallest, we can think of p 1 and p 2 being clustered together, and thenñ J ¼ ðp 1 þp 2 Þ=jp 1 þp 2 j. It is not hard to show that this minimizes T 1 and its final value is determined by the minimum metric as 2. The T 0 spectrum at NNLL 0 from SCET The all-orders parton-level factorization theorem for 0jettiness (aka beam thrust) is given by [36,51] Here, dσ B ij =dΦ 0 is the Born cross section for the ij → Z=γ Ã → l þ l − hard scattering. The hard function H ij ðQÞ contains the corresponding Born and virtual squared matrix elements, and the sum runs over all possible qq pairs ij ¼ fuū;ūu; dd;dd; …g. The B i ðt; xÞ are inclusive (anti) quark beam functions [36], with the momentum fractions x a;b given in terms of the total rapidity Y and invariant mass Q ≡ m l þ l − of the vector-boson final state by The beam functions are computed perturbatively in terms of standard parton distribution functions (PDFs) f j , sche- Finally, SðkÞ is the quark hemisphere soft function for beam thrust.
Equation (25) is derived using soft collinear effective theory (SCET) [52][53][54][55][56]. Its key feature is that each of H, B, S only depends on a single characteristic scale. Therefore, there are no large logarithms in their perturbative expansion when each is evaluated at its own characteristic scale, which is given by In Eq. (25), all ingredients have to be evaluated at a common scale μ, which is arbitrary and whose dependence exactly cancels between the different functions. The renormalization group evolution in the effective theory is then used to evolve each function from its own scale to μ. This gives the resummed T 0 spectrum used in Eq. (7) as where we have now written the convolutions between the different functions in schematic form. The renormalization group evolution factors U X ðμ X ; μÞ resum the large logarithmic terms, and the functions evaluated at their own scale provide the boundary conditions for the evolution. At NNLL 0 accuracy, we need the boundary conditions at two-loop order, and the evolution to three (two) loops in the cusp (noncusp) anomalous dimensions. All required expressions are known in the literature [43,[57][58][59][60][61][62][63][64][65] and we do not reproduce them here.
In practice, the canonical scales in Eq. (27) are appropriate in the resummation region, where the singular corrections dominate. In Fig. 1, we compare the singular, nonsingular, and full results for the T 0 distribution at NNLO 0 . We can see that up to T 0 ≲ 50 GeV the singular dominate and the nonsingular are suppressed by an order of magnitude or more. On the other hand, in the FO region for larger values of T 0 ∼ Q, the resummation must be turned off since here the singular terms being resummed become meaningless and there are large cancellations between the singular and nonsingular terms, which must be preserved to reproduce the correct FO result. This is clearly visible in Fig. 1 for T 0 ≳ 80 GeV, where the magnitude of the singular and nonsingular are larger than the full result (the singular turn negative at the dip at T 0 ≃ 70 GeV). For this reason the resummation must be turned off in this region. This is done by taking all scales to be equal to the A smooth transition between the canonical and FO limits is achieved by using profile scales [37,66], where μ B and μ S are smooth functions of T 0 , interpolating between the exact canonical scaling in Eq. (27) and the FO scale. The uncertainties from the resummation and the transition between resummation and FO can be estimated by choosing different sets of profile scales, which provide a sensible variation around the central scale choice. In other words, for each event we perform the calculation of its MC cross section with a different set of profile scales, producing different weights for the event, and the variation in these define our event-byevent uncertainties as described next.
Profile scales and their variations are used in many SCET calculations and have been shown to provide reliable perturbative uncertainty estimates for resummed and matched predictions for observables in a variety of contexts; see e.g. Refs. [37,[67][68][69][70][71]. We define the central profile scales using the same common profile function f run ðxÞ introduced in Ref. [68], namely This implies strict canonical scaling below x 1 and the resummation being turned off above x 3 . From Fig. 1 we choose for our central scale the parameters For the FO variations we vary μ FO up and down to 2Q and Q=2, taking the maximal absolute deviation in the results from the central value as the FO uncertainty. For our resummation uncertainties, we use the profile variations designed for T -like resummation variables from Ref. [70], which give us four independent up and down variations for μ B and μ S that probe deviations from the canonical scaling but never violate the parametric scaling μ 2 B ∼ μ S μ H and turn off beyond x 3 . In addition, we include two more profiles where we vary all x i transition points by AE0.05, giving us a total of six resummation profile variations. We then take the maximal absolute deviation in the result from the central value among all six profiles as the resummation uncertainty. The total perturbative uncertainty is obtained by adding the FO and resummation uncertainties in quadrature. Note that the profile scales depend on the underlying Φ 0 point through the value of Q, but we choose them to be independent of Y.
The same resummation formalism can be used to obtain an expression for the resummed cumulant of the T 0 spectrum in Eq. (7), Here, we first integrate the resummed expression for the T 0 distribution in Eq. (28), and then choose the scales using the same profile scales but as a function of T cut 0 . In particular, in the resummation region, the canonical values are Since the operations of integrating the factorization theorem and choosing the scales do not commute with one another, the resulting expression for the cumulant is not exactly equal to the integral of the T 0 spectrum, where the scales would be chosen before performing the integration. The resulting differences are formally always of higher order, beyond NNLL 0 in our case (see e.g. Ref. [72] for a detailed discussion). However, they can be numerically relevant, especially if one aims at preserving the NNLO 0 total cross section as accurately as possible.
In principle, one could enforce the equality between the cumulant and the integral of the spectrum in Eq. (32) exactly, by simply defining the T 0 spectrum as the derivative of the T cut 0 cumulant. However, this can give rise to unreliable uncertainty estimates in the spectrum as well as a poor description of the region of T 0 where the resummation is turned off and FO predictions are valid [37]. In GENEVA we circumvent these problems by adding to the T 0 spectrum the quantity which is explicitly beyond NNLL 0 order. The profile scales μ h are chosen to turn off the resummation somewhat earlier than the profile scales for the rest of the resummed calculation, in order to preserve the accurate spectrum prediction in the tail of the distribution. Furthermore, μ h and its derivative are smooth functions of T 0 , which is needed because the difference of resummed terms in Eq. (34) is proportional to dμ h =dT 0 [72]. The function κðT 0 Þ is a function that smoothly interpolates from a constant value for T 0 ≪ Q to zero at the point where μ h ðT 0 Þ ¼ Q, and further reduces the effects of the higher-order correction terms in the tail of the distribution. By adjusting the constant value at small T 0 within Oð1Þ factors, the overall size of these terms can be adjusted (in practice this constant value is approximately 2). These higher-order terms may be used to capture the correct FO scale uncertainties in the inclusive cross section upon integrating over the full spectrum. We use the same profile for each of the scale variation described above, but scaled by the value of μ FO =Q (so 1 for the central scale and resummation scale variations, 2 for the up FO scale variation, and 1=2 for the down FO scale variation). By adjusting the constant value of κðT 0 → 0Þ separately for each μ FO =Q value, we can correctly obtain the total integrated FO cross section for each FO scale variation.
Technically, each of the resummation scale variations should yield strictly the same total cross section, such that the resummation uncertainties vanish for inclusive observables. However, in practice, the profile scale variations are primarily meant and designed to capture the resummation uncertainties in the peak and transition regions. Implementing a constraint that they all reproduce exactly the same numerical integral of the resummed spectrum as the central profile is very challenging, especially since the resummation uncertainties in the peak region can be sizeable, while the inclusive cross section uncertainties are very small. [Alternatively, we could separately fine-tune the κðT 0 Þ function for each of the different resummation profiles.] Instead, what we do is simply drop the resummation uncertainties whenever we consider fully inclusive observables.

NNLO nonsingular corrections
As discussed in Sec. II A 1, the fact that we are resumming the T 0 spectrum and cumulant to NNLL 0 accuracy implies that there are no singular terms in the matching to the NNLO calculation. Therefore, GENEVA directly computes the fully differential inclusive zero-jet cross section at NNLO, up to power corrections in T cut 0 . The power-suppressed missing terms are precisely the zerojet nonsingular contributions in Eq. (10), which we can write as The functions f k ðT cut 0 ; Φ 0 Þ can still have divergences as T cut 0 → 0, but these divergences are at most logarithmic. Thus, these corrections vanish in the limit T cut 0 → 0, and for small enough T cut 0 they can be neglected, allowing us to obtain NNLO accuracy without doing the full NNLO calculation [which would determine the func- We show the size of the total correction terms Σ nons ðT cut 0 Þ in Fig. 2 for the NLO and pure NNLO contribution. At our default value of T cut 0 ¼ 1 GeV these terms are very small but also not completely negligible. In the current implementation we include the full NLO term, i.e., the f 1 ðΦ 0 ; T cut 0 Þ piece in dσ nons 0 ðT cut 0 Þ by performing the corresponding zero-jet NLO calculation on the fly. For the NNLO term, we neglect the Φ 0 dependence and only include its total integral. Technically, this is done by obtaining the pure NNLO contribution to Σ nons ðT cut 0 Þ from the known dependence of the total cross section on T cut 0 (i.e. analogous to how we obtain Fig. 2). This correction is then included through a simple rescaling of dσ MC 0 =dΦ 0 ðT cut 0 Þ.

The normalized splitting probabilities
In Eq. (8) we introduced a normalized splitting probability PðΦ 1 Þ, which is used to make the resummed T 0 spectrum fully differential in Φ 1 in Eq. (7). A similar splitting probability PðΦ 2 Þ is used in Eq. (15) to make the T 1 resummed term fully differential in Φ 2 . As mentioned before, the splitting probabilities are normalized, and provide the dependence on the two extra emission variables required besides T N to define a full splitting from Φ N → Φ Nþ1 . The two extra variables are chosen as usual as an energy ratio z and an azimuthal angle ϕ. In the collinear and soft limit, the z variable denotes the ratio z ¼ E M =E S for a splitting M → DS, where a mother particle M splits into a daughter particle D and a sister particle S. For initial state radiation (ISR) splittings, the daughter particle is the final-state particle, while for final state radiation (FSR) splittings it is always the gluon for q → qg splittings, the quark for g → qq ones, and the softer gluon for g → gg splittings.
We first need to decide which two particles in Φ Nþ1 are interpreted as arising from the splitting and which are then labeled with D and S. Here we have to choose the two particles which set the value of T N , which are the two particles that are closest according to the N-jettiness clustering metric, as described in Eq. (23).
The normalized splitting probability is then given by where the splitting function p sp ðz; ϕÞ depends on the chosen splitting sp which as explained above depends on the phase-space point sp ≡ spðΦ Nþ1 Þ.
With this definition, the splitting probability is indeed normalized: with AP sp ðz; ϕÞ denoting the unregularized Altarelli-Parisi splitting function, and fðx; μÞ the standard parton distributions function evaluated at the momentum fraction x and the factorization scale μ. The variables x M and x S denote the momentum fractions of the mother and sister particles in the ISR splitting M → DS.

T 1 resummation
The Sudakov factor U 1 ðT max 1 ; T 1 Þ in Eq. (15) resums the T 1 dependence to LL accuracy. We use the expression that is obtained from the T 1 factorization theorem [50] analogous to the one given for T 0 in Eq. (25), where all ingredients are calculated at tree level, and the running is performed at LL only. The Oðα s Þ expansion of this Sudakov factor has the simple expression The above minimal LL resummation is sufficient for our current purposes, where we focus on zero-jet resummation variables and inclusive observables and where we primarily need it to have a proper matching to the parton shower. A more precise description of one-jet resummation variables will require an improved T 1 resummation, which we leave for future work. Extending the T 1 resummation beyond LL accuracy is of course possible. One can also consider including the resummation through a POWHEG-like Sudakov factor, where the ratio of the real matrix element B 2 over the born B 1 would be exponentiated. We have verified that our underlying NLO 1 þ LL 1 matched calculation reasonably agrees with a POWHEG implementation of V þ 1 jet [73] in regions where the zero-jet resummation is not relevant. 6. NLO 1 calculation and phase-space map The last point addressed in this section concerns the phase-space map Φ T 1 ðΦ 2 Þ that is used to project a partonic phase-space point Φ 2 onto the Φ 1 point in the NLO 1 calculation; see Eqs. (11) and (17). It is also used in the T 1 resummation in Eq. (15). Just requiring having some IRsafe map as usual for fixed-order calculations [like for example the FKS map Φ C 1 ðΦ 2 Þ used in the C 2 subtraction term in Eq. (17)] is not enough anymore.
As discussed around Eq. (13), we have the additional constraint that the Φ T 1 map must preserve the value of T 0 in order for the NLO 1 calculation to be consistent with the T 0 resummation. This property was also important for the 1=2jet separation to preserve the T 0 resummation.
Technically, the only real-emission integral that we have to perform as part of the NLO 1 calculation is the integral in Eq. (17). We perform this integral numerically by Monte Carlo integration over the radiation phase space for a given Φ 1 event. The B 2 and C 2 terms in Eq. (17) are allowed to use different maps, since all IR-safe maps agree in the IR-singular limit. However, we have to be able to invert both maps, i.e., we have to reconstruct all Φ 2 points that would project to the given Φ 1 point. We then parametrize the radiation phase space for both maps in terms of a common set of radiation variables, which are then Monte Carlo sampled. While the usual FKS maps are designed to be invertible, the invertibility turns out to be a very nontrivial requirement on the Φ T 1 ðΦ 2 Þ map. Formally, we only need to preserve T 0 for the singular terms, i.e., we can use any map which preserves a variant of T 0 with the same singular structure to the order we want to carry out the resummation, while any nonsingular difference is captured by the nonsingular corrections. Therefore, to enable us to invert the map, we hold a recursive definition of T 0 fixed, which is effectively defined by the action of the map, The map itself is constructed as follows: We first cluster (by simply adding/subtracting four momenta) the two partons in Φ 2 that set the T 1 value, i.e., that are closest according to the N-jettiness metric. If the two final-state partons are clustered (FSR clustering), then this yields a Φ 1 point with an off-shell final-state momentum. If one of the final-state partons in Φ 2 is closest to a beam direction (ISR clustering), then this yields a Φ 1 point with one of the incoming momenta off shell (and off axis). Next, we compute T 0 on this off-shell Φ 1 point. Finally, we construct an on-shell Φ 1 point that has the same T 0 value and the same q μ of the vector boson. With q μ and T 0 given, the resulting Φ 1 point is fully determined by four-momentum conservation and on-shell conditions. Hence, by construction, the map also preserves the momentum of the vector boson q μ ðΦ 2 Þ ¼ q μ ðΦ T 1 ðΦ 2 ÞÞ. It is easy to see that this map is IR safe, and it is sufficiently simple that it can be inverted. For all ISR clusterings as well as FSR clusterings with both partons in the same hemisphere (defined by the vector-boson rapidity) this definition of T FR 0 is identical to the exact value of T 0 ðΦ 2 Þ. The only case where it differs is for FSR clusterings where the two final-state partons are in opposite hemispheres. In the T 0 -singular limit, this can only happen for soft emissions, and hence will in principle affect the constant terms in the two-loop soft function. (Collinear emissions in the beam functions always correspond to ISR clusterings.) However, this can only happen in a small region of phase space where the soft emitter was already accidentally close in rapidity to the vector boson. We therefore expect this to be a small effect, which we have verified numerically. We can therefore safely ignore it for our resummation and simply take it into account as part of the nonsingular matching corrections.
Finally, we need to comment on the Θ T ðΦ 2 Þ appearing in Eq. (12). It is easy to see that an on-shell massless Φ 1 point must satisfy T 0 ≤ jq T j. However, there can be Φ 2 points that can have a larger T 0 value, and hence such Φ 2 points cannot be projected onto a massless Φ 1 point. The Θ T ðΦ 2 Þ is defined to be 1 for points that can be projected and 0 otherwise. Equivalently, these points also cannot be reached by the inverse map from splitting a Φ 1 point, and so are never part of the integral in Eq. (17). Since the map is IR safe, the nonprojectable Φ 2 region is nonsingular and is included in Eq. (18) by the 1 − Θ T ðΦ 2 Þ constraint.

III. INTERFACING WITH A PARTON SHOWER
In the previous section we have given all required formulas for the jet cross sections dσ MC 0 , dσ MC 1 , and dσ MC ≥2 . As discussed in detail there, these jet cross sections include the contributions of higher multiplicity phase-space points, as long as the jet resolution variable for these points in phase space is below the T cut k value. In Table I we summarize how the phase-space points of different multiplicities contribute to the given jet cross section. Only the first three columns were necessary for the partonic calculation of the previous section. The last column shows instead the constraints to be imposed on events with higher multiplicities. We detail how these affect the development of the parton shower below.
The purpose of the parton shower is to make the calculation fully differential in the higher multiplicities, which can be viewed as filling the exclusive zero-and 1-jet bins with radiation, as well as adding extra jets to the inclusive two-jet multiplicity. A parton shower acts on a partonic event by adding extra radiation in a recursive and unitary fashion. If the evolution variable of the parton shower were chosen to be N-jettiness, 1 a single emission would be given by where T max N ∼ OðT N−1 Þ is the maximum value of T N that can be reached given the value T N−1 that characterizes the hardness of the configuration before the emission and Λ is the shower cutoff. This means that the parton shower either keeps the N-body event unchanged with the no-emission probability determined by the Sudakov factor U N ðT max N ; ΛÞ, or it adds an extra emission according to the probability which we define here in terms of the derivative of the Sudakov form factor and the normalized splitting probability given in Sec. II B 4. From Eq. (42) and using the analog of Eq. (20) it is clear that the parton shower is unitary. Defining the recursive shower can now be written as Note that the definition of the parton shower requires a phase-space map. This is clear from Eq. (44), which depends on dΦ Nþ1 =dΦ N . The above equation assumes that the evolution variable of the parton shower is N-jettiness. PYTHIA8, which we interface to for showering, as well as any other parton showers currently existing, uses a different evolution variable, e.g. the transverse momentum or the angular separation of the emission. However, one can imagine taking the output of any of these showers and reclustering the partons according to the N-jettiness metric defined in Eq. (23). The resulting history of splittings is equal at LL with the recursive expression given in Eq. (45). In the remaining discussion, we work with Eqs. (42) and (45), and it is understood that this is valid only after the aforementioned reclustering. Note that this construction implies a recursive definition of N-jettiness, which in general is not identical to the standard definition given in Sec. II B 1. It does agree, however, with the T FR 0 definition discussed in Sec. II B 6, which is used in GENEVA.
The parton shower is not allowed to affect results of the jet cross sections at the accuracy they were calculated in the previous section. Starting from this condition, let us first consider the constraints we would have to impose on a parton shower that is strongly ordered in N-jettiness, as discussed above. Strong ordering in T N means that the T N value at each emission is much smaller than the T max N upper limit determined by the T N−1 value of the previous emission. In this case, we would need to ensure that at each stage in the shower each partonic phase-space point satisfies the constraints given in Table I. For example, when showering a zero-jet event defined by Φ 0 , the phase-space point Φ 1 reached after the first emission must have T 0 ðΦ 1 Þ < T cut 0 and Φ 1 must be projectable onto Φ 0 . Each subsequent emission only needs to satisfy T 0 ðΦ N Þ < T cut 0 . When showering a one-jet event defined by Φ 1 , the Φ 2 point obtained after the first emission would need to have T 1 ðΦ 2 Þ < T cut 1 , and Φ 2 would have to be projectable onto Φ 1 with the T 0 -preserving map discussed in Sec. II B 6. This implies that the Φ 2 point obtained after one emission has to satisfy T 0 ðΦ 2 Þ ¼ T 0 ðΦ 1 Þ as well as q μ ðΦ 2 Þ ¼ q μ ðΦ 1 Þ, q μ being the vector-boson four momentum. Each subsequent emission then only needs to satisfy T 1 ðΦ N Þ < T cut 1 . For a shower ordered in a different evolution variable, very similar restrictions hold; however the special conditions regarding the phase-space maps are not required on the first emission, but rather on the emission with the largest value of the jet resolution scale, which can happen at a much later stage in the shower. Failing to properly account for this could affect the accuracy at which observables are predicted and ultimately spoil the leading-logarithmic terms and the color-coherence effects that are built into the shower Sudakov factors.
A solution to this problem is to ensure that the first emission has indeed the largest value of the jet resolution scale. This can be achieved by performing the first emission in GENEVA, with a probability obtained from a simple LL Sudakov factor and using a splitting that preserves the map that was so carefully constructed. One would then only need to ensure that subsequent emissions from the parton shower have a lower jet resolution scale than the first, 1 Here and in the following we treat T N as the jet resolution scale of a single emission. Since T N is a global variable that depends on all emissions, this warrants some explanation. What we intend here proceeds from the condition that for massless particles T N ðΦ Nþ1 Þ is always the total plus component of the closest pair of particles, using e.g. the metric defined in Eq. (23), relative to the direction of the sum of their momenta. Moreover, T N → 0 when Φ Nþ1 → Φ N . Near the singular limit, where the shower is important and a good approximation of the underlying physics, one can thus assume that the direction of the axes entering the definition of T N is aligned to the direction of the N hard partons. In this limit, a singular emission does not change these directions. Therefore, the value of T N does in fact represent the hardness of the emission and can be used as the evolution variable. which can be achieved by properly setting the starting scale for the showering and applying standard veto techniques.
The general idea of this approach is to split an N-jet differential cross section into two pieces using a Sudakov factor U N and splitting probability P to define a branching and no-branching probability as in a unitary parton shower. The first term is now an exclusive N-jet cross section with a jet resolution at Λ N , which can be much smaller than the T cut N value used in GENEVA. This process can also be iterated further. Applying this to our Φ 0 and Φ 1 events, and adding all emissions up to two partons, we find TheT max 1 in the theta-function in the last line above is either T cut 1 or T max 1 , depending on whether the derivative is on the first or the second term of Eq. (49). Once the first emissions have been performed as described above, the only remaining constraint on the parton shower is what is given in the last column of Table I, namely the constraint on the jet resolution variable. In practice, we only apply this to the Φ 1 events, since as discussed below, the showering of the Φ 0 events is not affecting the accuracy of the T 0 distribution.
In the remainder of this section we give an argument of why the parton shower acting on the jet events as described above does not affect the NNLL 0 þ NNLO accuracy. The Φ 0 events include all phase-space points for which T 0 < T cut 0 , and by definition GENEVA is only predicting the normalization, and not the distribution in T 0 . Therefore, as long as we constrain the parton shower to only shower these events up to T cut 0 , the shower will exactly provide the missing T 0 shape below T cut 0 which is integrated over by GENEVA. On the other hand, due to the unitarity constraint, the shower will not change the normalization of the events below T cut 0 , which is thus still the NNLL 0 þ NNLO 0 one predicted by GENEVA. By performing the first emission of our Φ 1 events as described above, and taking Λ 1 as small as possible (we use a value of 0.1 GeV), the remaining one-jet cross section is only about 0.1% of the total cross section. The theoretical problem of how the shower affects the accuracy of this small contribution to the total cross section, which requires imposing all the constraints detailed above, can then be ignored for practical purposes.
What is therefore left to show is that the shower does not affect the T 0 distribution when showering Φ 2 events. The T 0 distribution of the original two-body phase space point is given by while the expression after the first emission done by the shower is given by Clearly, thanks to the infrared safety of T 0 , in the soft and collinear limit of Φ 3 one has If the phase-space map of the parton shower were to preserve the value of T 0 , the relation Eq. (53) would remain exact even away from the collinear and soft limits. The first emission done by the shower would then not change the T 0 spectrum If the phase-space map of the shower does not preserve T 0 , instead, a difference between these two values exists. We can assume this difference to be proportional to T 2 ðΦ 3 Þ, which measures the distance from the soft and collinear limits Here aðΦ 3 Þ is well behaved in the singular limit of Φ 3 , because all the singular behavior is incorporated in T 2 ðΦ 3 Þ → 0. For the sake of simplicity and to avoid rendering the notation too heavy we drop the ðΦ 3 Þ dependence of a and treat it as constant in the remainder of this section. The reader can easily convince themselves that the argument given below will not be affected by this approximation. This allows us to write Thanks to the normalization of the splitting probabilities, one can easily perform the integrals over the two radiation variables other than T 2 to obtain Next, the delta function is Taylor expanded to obtain where we have defined the average T 2 value In the last lines of Eq. (58), we have written hT 2 iðT 0 Þ to make the dependence of the average value of 2-jettiness on T 0 explicit. In fact, one can derive this dependence using a simple LL Sudakov factor for the parton shower. Let us first consider a single emission governed by with C being the appropriate color factor. We can now easily show that Iterating this over two emissions, we find where we used the scaling T max N ∼ T N−1 . We can also numerically extract the dependence of hT 2 i on T 0 from an actual shower in the following way: Starting from the showered expression dσ S 3 =dΦ 3 defined as in Eq. (44), we have This gives Taking the Φ 2 events from GENEVA, running them through PYTHIA8, and calculating this ratio for various values of T 0 , we obtain a very good fit to a straight line with Taking advantage of the linear scaling of hT 2 i with T 0 , and dropping the explicit dependence in the notation from here on, we can now rewrite Eq. (58) as This leads to where Since the dominant contribution to dσ=d ln T 0 is given by dσ d ln T 0 ∼ −α s ln T 0 e −α s ln 2 T 0 ; ð69Þ we find This gives our final result for the change of the T 0 spectrum after the first emission done by the parton shower Since the parton shower is strongly ordered in T N , in the sense discussed at the beginning of this section, all the subsequent emissions will not alter this result to the order we are working. This result can be compared to the relative uncertainty in our NNLL 0 calculation of the T 0 spectrum, which can be estimated considering the dominant term beyond the NNLL 0 order. This term scales as α 3 s =T 0 . Comparing this to the first dominant one that we include in the NNLL 0 prediction, we find the relative uncertainty to be Upon inspection of Eqs. (71)-(72), we can see that the shift in the spectrum due to PYTHIA's showering is of the same size as the missing higher-order terms in the perturbative calculation. Therefore, the showering only affects our prediction beyond the claimed accuracy.

IV. COMPARISON WITH DEDICATED PERTURBATIVE CALCULATIONS
In this section, we present the results of GENEVA and the comparisons with dedicated FO and resummed predictions. All calculations are performed for pp collisions with 7 TeV center-of-mass energy. We use the CT10NNLO [74] set from LHAPDF6 [75] as our default PDFs, together with its default value of α s ðM Z Þ ¼ 0.118. For the tree-level matrix elements, including color and spin-correlations, we interface to NJET [76,77] through the BLHA interface [78,79]. The one-loop matrix elements are instead taken from the POWHEG-BOX [7,73] (with original routines from MCFM [80]). We allow for the full interference effects between Z=γ Ã and restrict the invariant mass for the dilepton pair to the range 60 < m l þ l− < 120 GeV during the generation of events. The other parameters relevant for our calculation are The Z couplings are given by where l=q denotes the given left or right component of a lepton or a quark and e ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 4πα em ðM Z Þ p . For the parton shower, we use PYTHIA8.205 [34,35], with a slightly modified version of Tune 4C, namely with α s ¼ 0.118 for both the ISR and FSR shower, p ref T 0 ¼ 4 GeV, and the hard primordial k T ¼ 2.5 GeV. Unless otherwise stated, the GENEVA predictions include showering and hadronization effects, but presently do not include multiple parton interactions.
We first present the T 0 spectrum, and show that we exactly reproduce the input resummed spectrum, and that the parton shower does not modify this distribution as required.
To validate the FO accuracy of GENEVA, we next consider the rapidity distribution of the Z boson, as well as that of the negatively charged lepton from the Z decay. These distributions are known to NNLO accuracy. As discussed in Sec. II, GENEVA reproduces all Oðα 2 s Þ contributions up to small power-suppressed effects, and these observables allow us to test the size of these effects. In addition, it allows us to validate that the procedure described in Sec. II B 2 to turn off the resummation and numerically reproduce the total cross section from the integrated resummed spectrum works.
Finally, to study the accuracy of our resummation, we show the results for various observables that are sensitive to the higher logarithms included in GENEVA. The first is the zero-jet cross section, which only contains events with jets below a given p T cut and is known analytically to NNLL. The corresponding observable in Higgs production is an important ingredient in Higgs analyses. The other two resummation-sensitive observables we consider are the transverse-momentum distribution of the Z boson, as well as the ϕ Ã distribution [81]. For all these observables we compare our results to the analytic resummed results at NNLL matched to NNLO 0 or NLO 1 .

A. The T 0 spectrum
The zero-jet resolution variable chosen in GENEVA is 0jettiness, such that partonic configurations with T 0 < T cut 0 are part of the exclusive zero-jet MC cross section, while partonic configurations with T 0 > T cut 0 contribute to the inclusive one-jet cross section. Since the dependence on the zero-jet resolution variable at NNLL 0 accuracy is used as input in GENEVA, it must reproduce the T 0 spectrum at NNLL 0 accuracy both for the central value and for the event-by-event scale uncertainties.
In Fig. 3 we show the T 0 spectrum in the peak, transition, and tail regions. All the curves are cut off at the T cut 0 ¼ 1 GeV value used in the simulation. We show the analytic expressions with a red band, where the uncertainties are obtained using the same profile scale variations discussed in Sec. II B 2. The green histogram with error bars shows the results of GENEVA before attaching the parton shower. One can see that GENEVA is reproducing the analytic results exactly, as it should, for both the central value and the scale variations. Finally, the blue histogram shows the results of GENEVA after the parton FIG. 3 (color online). Validation of the T 0 spectrum in GENEVA. The analytic NNLL 0 þ NNLO T 0 resummation is compared to the partonic and showered results of GENEVA, but before the addition of nonperturbative effects.
FIG. 4 (color online). Comparison of the showered and hadronized T 0 spectra in GENEVA. shower. As required, the parton shower changes the T 0 spectrum only within the perturbative uncertainties. In fact, over most of the T 0 range the spectrum is essentially unchanged.
The nonperturbative effects are shown in Fig. 4, where we compare the predictions for the T 0 spectrum after the shower (again in blue) and after the addition of proton remnants, intrinsic k T smearing, and hadronization by PYTHIA8 in orange. These effects change the T 0 distribution significantly in the peak region, while they become power corrections in the transition and tail regions. This behavior is precisely as dictated by factorization, from which one expects that these effects should behave for T 0 analogous to thrust in e þ e − . Comparing to our e þ e − results [24], this is precisely what we observe. A benefit of the GENEVA approach is that it allows us to directly combine the higher-order analytic resummation with these nonperturbative corrections provided by the hadronization model in PYTHIA8.

B. Partonic NNLO 0 observables
We now show that GENEVA reproduces fully inclusive observables at NNLO accuracy by comparing to dedicated NNLO calculations. In this section, we only consider the profile scale variations that reproduce the FO scale variations, as described in Sec. II B 2.
In Figs. 5-6 we show the result for the rapidity distribution of the vector boson and the negatively charged lepton from its decay, respectively. The orange band shows the NNLO result from DYNNLO [82]. We show the results of GENEVA as a black histogram, with the error bars representing FO uncertainties as described above. In the lower part of each plot, we show the ratio to the DYNNLO central value.
The central value of GENEVA agrees very well with the fixed NNLO prediction. The perturbative uncertainties from GENEVA are also in reasonable agreement with those from DYNNLO. The few fluctuations observed in the plot are of a statistical nature, as evidenced by the fact that they grow larger toward larger rapidities, where the statistics is poorer. The rapidity distribution of the vector boson has also been validated against the independent NNLO calculation provided by VRAP [83].
In Fig. 7 we show the result of the transverse-momentum distribution of the negatively charged lepton. For p Tl < m l þ l − =2 this observable is a true NNLO distribution and GENEVA agrees very well with the NNLO prediction. The region above and below m l þ l − =2 is very sensitive to Sudakov shoulder logarithms [84]. It is known that the FO calculations perform very poorly in this region and fail to provide an accurate description of the physics. On the other hand, GENEVA will have some of these logarithms resummed and should therefore provide a more accurate prediction. Finally, the region above p Tl > m l þ l − =2 can only be populated if the transverse momentum of the vector boson is nonzero, making this region only NLO accurate. Close to p Tl ≳ m l þ l − =2 GENEVA lies above the FO prediction and converges back to the FO result for large values of p Tl . This is likely still the effect from the Sudakov shoulder at m l þ l − =2. Since the cross section above is much lower, any spillover from below m l þ l − =2 caused by the resummation can have a large relative effect.

C. Zero-jet observables with resummation
Since we resum our zero-jet resolution variable T 0 to NNLL 0 , it is interesting to study how accurately GENEVA is able to predict other observables sensitive to the 0=1-jet separation. Since the distribution for other zero-jet resolution variables is only obtained indirectly, via the T 0 spectrum made fully differential by using splitting functions and a parton shower that are only correct to LL, we cannot formally claim the same NNLL 0 accuracy for these observables. However, since the overall distribution of events in the 0=1-jet region of phase space is clearly improved, we expect some of this accuracy to carry over to other observables as well. Hence, we expect to get predictions that are numerically much closer to NNLL resummed results also for other zero-jet observables. The same behavior was already observed for e þ e − [24].
In Fig. 8 we show the transverse-momentum distribution of the vector boson compared to its analytic NNLL resummation from DYQT [85,86]. The predictions of DYQT have been manually switched to agree with the FO results in the tail, according to the procedure outlined in Ref. [87]. We see that GENEVA agrees reasonably well within the perturbative uncertainties with the NNLL resummed result. 2 The NLL result from DYQT has a significantly different shape, and the GENEVA prediction is certainly in much better agreement with the NNLL predictions than the NLL ones.
A similar prediction for the transverse-momentum distribution, but fully differential on the momenta of vector boson decay products has been presented in Ref. [88]. This allows the direct comparison including the acceptance cuts used by experimental analyses. Using the same lepton cuts as in that study, we show our comparison in Fig. 9. Again, we observe a fairly good agreement with the analytic NNLL prediction matched to NLO 1 . Another variable, quite similar to the transverse momentum of the vector boson, is the ϕ Ã between the two leptons, with the precise definition of ϕ Ã given in [81].  Although GENEVA's perturbative uncertainties appear smaller than DYQT's at very low q T , this should not be misinterpreted as GENEVA being more accurate here. This is mostly an artifact of smearing out the uncertainties from a range of T 0 values, which tends to reduce their numerical impact. In addition, we have not included here any uncertainties associated with the interface to the parton shower and the showering itself. Such uncertainties should eventually also be included as an additional part of the (resummation) uncertainties. As this will require a dedicated study, we leave this for future work.
to the NNLL þ NLO 1 calculation of Ref. [88] is shown in Fig. 10, and we again observe good agreement.
Finally, we show the result for the exclusive zero-jet cross section as a function of p cut T in Fig. 11, where the zero-jet sample is defined as all events containing no jets with p T > p cut T . The jets are reconstructed with the anti-k T algorithm [89] as implemented in FASTJET [90,91], within a radius R ¼ 0.4. We find good agreement between GENEVA and the dedicated NNLL þ NNLO calculation given by JETVHETO [92] within the perturbative uncertainties. For this plot, we use the FO scale uncertainties discussed in Sec. II B 2, such that the uncertainties at large p cut T are estimated correctly and thus precisely reproduce those of JETVHETO. At small p cut T they are now underestimated and here the resummation uncertainties should be added. The better agreement with the NNLL þ NNLO prediction compared to the lower order NLL þ NLO one, especially in the large p cut T region, is of course driven by the correct inclusion of the NNLO corrections in GENEVA.

V. COMPARISON TO LHC DATA
In this section, we compare the results from GENEVA with measurements using data collected during the 7 TeV LHC run I. We use GENEVAþPYTHIA8 to produce fully showered and hadronized events in HEPMC format, and then use the RIVET routines provided by the experimental collaborations to produce the histograms that can be compared with the experimental measurements. By using RIVET [93], we ensure that event selection cuts and object definitions agree with what was used in the experimental analyses. We use the same run that was used for the comparison to the perturbative calculations, such that all parameters are identical to those given in Sec. IV.
In addition, we also include in our data comparisons GENEVA results with a lower value of α s ðM Z Þ ¼ 0.1135, since as we will see, the value of α s ðM Z Þ can have a noticeable effect on resummation-sensitive observables. In the following, the results with α s ðM Z Þ ¼ 0.118 will be shown by the orange bands, while those with α s ðM Z Þ ¼ 0.1135 are shown by a green band. For the latter results, we only change the α s ðM Z Þ value, holding all other parameters including the PDF set fixed. 3 This value is motivated by the fact that all α s ðM Z Þ determinations from LEP data that include higher-order resummation together with some form of analytic treatment of nonperturbative hadronization effects find considerably lower values of α s ðM Z Þ than 0.118 [37][38][39][40][94][95][96]; see also Ref. [97] for determinations from DIS. Consistent with these findings, in Ref. [24] it was observed that the higher-order resummation using α s ðM Z Þ ¼ 0.1135 coupled with PYTHIA8's hadronization model gave excellent agreement with data for e þ e − → 2-jet resummation observables, namely thrust, C-parameter, heavy-jet mass, and jet broadening. On the other hand, the predictions with α s ðM Z Þ ¼ 0.118 gave consistently spectra that were somewhat too hard, i.e., shifted from lower to higher values compared to data. We observe precisely the same qualitative trend here.
The rapidity distribution of the vector boson was measured by ATLAS [98], CMS [99], and LHCb [100]. Since we already saw that GENEVA agrees very well with the NNLO predictions, these comparisons are mostly driven by the PDFs we use. We show the comparison to the rapidity spectrum from ATLAS and CMS in Figs. 12-13. The normalized results from CMS are reproduced very well for the entire Y range. The ATLAS results are not normalized, and while GENEVA agrees very well at low Y values, it slightly overshoots the data (by about two standard deviations) for larger Y. A similar discrepancy was observed in Ref. [98] for the MSTW08, HERAPDF1.5, and ABKM09 PDF sets. The comparison with LHCb is shown in Fig. 14. The agreement at moderate values of Y is good here too. For forward events with Y > 3.25 our predictions show a trend to be lower than the data. This discrepancy was FIG. 12 (color online). Comparison of GENEVA with the rapidity of the vector boson from the ATLAS study in [98]. The GENEVA results with default values of α s ðM Z Þ and with α s ðM Z Þ ¼ 0.1135 are shown in orange and green, respectively, while the ATLAS points are shown in black.
FIG . 13 (color online). Comparison of GENEVA with the rapidity of the vector boson from the CMS study in [99]. The GENEVA results with default values of α s ðM Z Þ and with α s ðM Z Þ ¼ 0.1135 are shown in orange and green, respectively, while the CMS points are shown in black. 3 We hasten to add that using a different α s ðM Z Þ than the default value provided by the PDF set is not problematic as it may seem at first. In the resummation-dominated region, this only leads to a small inconsistency inside the beam functions, where the terms that are supposed to cancel the PDFs μ F dependence will now use a different α s value leading to a small noncancellation in the coefficients of some single-logarithmic terms. It is safe to assume that this will not be numerically relevant compared to the dominant α s dependence coming from the double-logarithmic beam function terms, as well as the hard and soft function, and all the resummation kernels. All of these pieces are completely unrelated to the PDFs, and are in fact the same as for the corresponding T 2 resummation in e þ e − collisions. already noted in the LHCb paper [100] when comparing with the NNLO calculation provided by FEWZ [101]. As we would expect, the lower α s ðM Z Þ value makes almost no difference for the inclusive predictions.
A measurement of the production cross section of a Z boson in association with jets was presented by ATLAS in Ref. [102], for both inclusive and exclusive jet cross sections. We compare our results with the measurements having up to two jets in Fig. 15. We choose to limit ourselves to up to two jets because any additional jet would only be provided at LL accuracy by PYTHIA8. Our predictions for the zero-jet cross section agree well with data, not only for the inclusive cross section, but also for the exclusive one, where resummation plays a role. The uncertainties for our predictions of the inclusive Z þ 1jet cross section are larger, since they are only predicted to NLO 1 , though they can still benefit from the resummation of the 0=1-jet boundary. The separation into exclusive Z þ 1 jet and inclusive Z þ 2 jet is only at LL accuracy, and the FO accuracy of the inclusive Z þ 2-jet cross section is only LO 2 , with correspondingly larger perturbative uncertainties. GENEVA agrees well with the data, somewhat better for lower α s ðM Z Þ.
In Ref. [102], ATLAS also presented a measurement of the transverse-momentum distribution of the hardest jet, and the comparison to these results is shown in Fig. 16. The predictions from GENEVA are in good agreement with the data. Below p jet T ≲ m Z the agreement is noticeably better for lower α s ðM Z Þ. Above p jet T ≳ 300 GeV the predictions tend to be higher than the data, but still consistent within the larger uncertainties. This could be due to the fact that we use a renormalization scale of μ ¼ m l þ l − , while at such large transverse momenta a better choice might be μ ¼ p jet T .
FIG. 14 (color online). Comparison of GENEVA with the rapidity of the vector boson from the LHCb study in [100]. Next, we compare our predictions to the ϕ Ã between the leptons, and the transverse momentum of the vector boson. Both of these observables are zero-jet resummation variables, since at low values they are dominated by events without any hard emissions. Both of these distributions are quite sensitive to the choice of parameters of the parton shower and the nonperturbative model used in PYTHIA8. Without higher-order perturbative corrections included, the MC tune will partially adjust the available parameters to mimic missing higher-order perturbative effects. This implies that when using PYTHIA8 in conjunction with GENEVA, which includes much more higher-order perturbative information, a retuning of the parameters becomes necessary. We stress that no attempt at a systematic retuning of PYTHIA8 has been done for this work, and we expect that a dedicated tune of PYTHIA8 together with GENEVA will improve the data agreement.
The ϕ Ã distribution has been measured by ATLAS in Ref. [103] and LHCb in Ref. [100]. The measurement from ATLAS is a normalized spectrum, while LHCb quotes the unnormalized distribution. The comparisons with these two measurements are shown in Figs. 17-18. GENEVA agrees well with the results from LHCb, though the measurement has relatively large uncertainties. Comparing with the much more precise results from ATLAS, GENEVA predicts a wider distribution, such that our predictions are below the data in the peak region for ϕ Ã ≲ 0.1, and above the data for larger values of ϕ Ã .
Discrepancies with this same trend and at a similar level were observed for other MC predictions in Ref. [103]. The agreement with data is considerably better when choosing a lower value of α s ðM Z Þ ¼ 0.1135.
As a final comparison we consider the transverse momentum (q T ) distribution of the Z boson. The q T FIG. 17 (color online). Comparison of GENEVA with the ϕ Ã distribution from the ATLAS study in [103]. The GENEVA results with default values of α s ðM Z Þ and with α s ðM Z Þ ¼ 0.1135 are shown in orange and green, respectively, while the ATLAS points are shown in black.
FIG. 18 (color online). Comparison of GENEVA with the ϕ Ã distribution from the LHCb study in [100]. The GENEVA results with default values of α s ðM Z Þ and with α s ðM Z Þ ¼ 0.1135 are shown in orange and green, respectively, while the LHCb points are shown in black.
FIG . 19 (color online). Comparison of GENEVA with the transverse momentum of the vector boson from the CMS study in [99]. The GENEVA results with default values of α s ðM Z Þ and with α s ðM Z Þ ¼ 0.1135 are shown in orange and green, respectively, while the CMS points are shown in black. distribution was measured by CMS in Ref. [99] and by ATLAS in Ref. [104], where the latter also measured the q T spectrum in different bins of Y. The comparison of our results with these measurements is shown in Figs. 19-20. The situation is similar to the ϕ Ã distribution discussed above. For α s ðM Z Þ ¼ 0.118, GENEVA is below the data in the peak region q T ≲ 10 GeV, while it is above the data in the transition and FO regions for q T ≳ 10 GeV. This is of course not unexpected, since the observables q T and ϕ Ã are highly correlated. The same level of disagreement was already observed in Ref. [104], where the data were compared with the NNLL þ NLO 1 results of Ref. [92]. As for the ϕ Ã distribution, the agreement between GENEVA and the data is noticeably improved for a lower α s ðM Z Þ ¼ 0.1135.
FIG. 20 (color online). Comparison of GENEVA with the transverse momentum of the vector boson both inclusive (top-left) and in different bins of vector-boson rapidity Y from the ATLAS study in [104]. The GENEVA results with default values of α s ðM Z Þ and with α s ðM Z Þ ¼ 0.1135 are shown in orange and green, respectively, while the ATLAS points are shown in black.

VI. CONCLUSIONS
We have presented a combination of the fully differential NNLO calculation for Drell-Yan production pp → γ=Z → l þ l − combined with the NNLL 0 resummation of 0jettiness and interfaced with the parton shower provided by PYTHIA8 within the GENEVA Monte Carlo framework. The starting point of GENEVA is the formulation of Monte Carlo cross sections, which are IR-safe partonic jet cross sections, defined using a jet resolution variable T N , which is chosen to be N-jettiness in our implementation. We include zero-jet, one-jet, and two-jet cross sections at NNLO 0 , NLO 1 , and LO 2 . Furthermore, the dependence on the zero-jet resolution variable T 0 is resummed to NNLL 0 accuracy, while that of the one-jet resolution variable T 1 is presently resummed to LL.
Interfacing these partonic results with a parton shower such as PYTHIA8 requires a careful treatment. First, one needs to deal with the fact that the parton shower evolution variable is different from N-jettiness. Second, care has to be taken that the parton shower's lower resummation accuracy and different phase-space map does not destroy the higher logarithmic accuracy of the calculated Monte Carlo cross sections. We discussed in detail how to solve these issues by performing the first two shower emissions by hand using the jet resolution parameter as the ordering variable. We then showed that any subsequent showering by PYTHIA8 does not affect the formal perturbative accuracy included in GENEVA. We have validated the FO perturbative accuracy by comparing our results for inclusive zero-jet observables to DYNNLO. We agree with this dedicated NNLO calculation within the small perturbative (and statistical) uncertainties. We also studied how GENEVA's improved perturbative accuracy in the resummation region from resumming 0-jettiness to NNLL 0 translates to other zerojet resummation variables, such as the transverse momentum of the vector boson q T , the ϕ Ã between the two leptons, and the exclusive zero-jet cross section as a function of the jet p cut T . Since GENEVA partially relies on the parton shower for these observables, they formally do not have NNLL 0 accuracy. Nevertheless, we find that GENEVA reproduces dedicated NNLL resummations for these three observables rather well, and in particular the GENEVA predictions are much closer to the exact NNLL þ NLO 1 results than to the NLL þ LO 1 results. This is a clear indication that the gain in resummation accuracy for T 0 , when implemented into a fully exclusive prediction also translates into more accurate predictions for other zero-jet resolution variables.
Finally, we presented a comparison of GENEVA with measurements by ATLAS, CMS, and LHCb using the 7 TeV LHC data. Since GENEVA agrees with the NNLO results for the rapidity distribution of the vector boson, the agreement between GENEVA and the LHC data is very similar to that observed when comparing to other NNLO calculations. Here, we note that experimental measurements of the lepton p Tl spectrum would be very valuable.
We find good agreement between GENEVA and the ATLAS measurement of exclusive jet cross sections and of the transverse-momentum distribution of the hardest jet. For the transverse momentum of the vector boson or the ϕ Ã between the leptons our predictions when using the default value of α s ðM Z Þ ¼ 0.118 are lower in the peak region and higher in the transition region behind the peak compared to the measurements. The same trend has been observed before by other MC predictions and also dedicated resummed calculations for these observables. When using a lower value α s ðM Z Þ ¼ 0.1135, we observe a noticeable improvement in agreement with data for essentially all resummation-sensitive observables. A similar effect was observed previously for a variety of e þ e − 2-jet resummation observables. We believe this deserves further attention.
We encourage the experiments to also measure T 0 -like jet-based variables (see Ref. [70]) and perform differential measurements of resummation observables at much higher m l þ l − values. Since the resummation type and regime for these cases are very different, they would add valuable complementary information for theory comparisons. In general, continued precise measurements of resummation observables at the LHC are essential for our understanding of higher-order QCD effects in exclusive and differential observables, which ultimately will also provide important inputs for the theoretical interpretation of Higgs measurements and new-physics searches.