Matching NNLO predictions to parton showers using N3LL color-singlet transverse momentum resummation in geneva

: We extend the geneva Monte Carlo framework using the transverse momentum of a color-singlet system as the resolution variable. This allows us to use next-to-next-to-next-to-leading-logarithm (N3LL) resummation via the radish formalism to obtain precise predictions for any color-singlet production process at the fully exclusive level. Thanks to the implementation of two different resolution variables within the geneva framework, we are able to assess the impact of such a choice on differential observables for the first time. As a first application, we present predictions for Drell-Yan lepton-pair production at next-to-next-to-leading order (NNLO) in QCD interfaced to a parton shower simulation that includes additional all-order radiative corrections. We provide fully showered and hadronized events using pythia8, while retaining the NNLO QCD accuracy for observables which are inclusive over the additional radiation. We also show that it is possible to obtain a numerically good agreement between showered geneva predictions and the N3LL resummation for the transverse momentum spectrum by choosing a more local recoil scheme. We compare our final predictions to LHC data at 13 TeV, finding good agreement across several distributions. We extend the GENEVA Monte Carlo framework using the transverse momentum of a color-singlet system as the resolution variable. This allows us to use next-to-next-to-next-to-leading-logarithm (N 3 LL) resummation via the RADISH formalism to obtain precise predictions for any color-singlet production process at the fully exclusive level. Thanks to the implementation of two different resolution variables within the GENEVA framework, we are able to assess the impact of such a choice on differential observables for the first time. As a first application, we present predictions for Drell-Yan lepton-pair production at next-to-next-to-leading order (NNLO) in QCD interfaced to a parton shower simulation that includes additional all-order radiative corrections. We provide fully showered and hadronized events using PYTHIA 8 , while retaining the NNLO QCD accuracy for observables which are inclusive over the additional radiation. We also show that it is possible to obtain a numerically good agreement between showered GENEVA predictions and the N 3 LL resummation for the transverse momentum spectrum by choosing a more local recoil scheme. We compare our final predictions to LHC data at 13 TeV, finding good agreement across several distributions.


I. INTRODUCTION
The rich and vast physics program at the Large Hadron Collider (LHC) has delivered an outstanding amount of data so far. Thanks to these impressive performances, particle physics has entered an era where high precision is ubiquitous. With the forthcoming start of the Run 3 and the future upgrade to the High-Luminosity LHC, the amount of data that will be collected will increase significantly. This ground-breaking machine will thus be able to detect extremely rare phenomena and further improve the precision of measurements.
Theoretical predictions of the Standard Model processes are one of the most fundamental ingredients for the interpretation of collider data. The vast majority of experimental analyses rely on theoretical predictions in the form of perturbative calculations at parton level or in combination with parton showers (PSs) in fully exclusive Monte Carlo (MC) event generators. The advantage of MC event generators is the ability to produce hadron-level events that can be directly interfaced to detector simulations.
In order to completely exploit the precision of present and future experimental data, it becomes therefore mandatory to refine MC event generators including the state-ofthe-art corrections available. For this reason, in recent years there has been increasing interest in obtaining theoretical predictions where fixed-order predictions at next-to-nextto-leading-order (NNLO) accuracy are consistently matched with PS simulations. There are currently four available methods which can reach NNLO þ PS accuracy [1][2][3][4][5][6][7], which have been initially applied to 2 → 1 processes such as Higgs [1,5] and Drell-Yan (DY) [3,4,8] production. The number of applications at hadron colliders has been rapidly growing in the last few years, and the methods have been extended to more complex processes such as Higgsstrahlung [9][10][11], hadronic Higgs boson decays [12,13], W þ W − [14], diphoton [15], Zγ [16], and tt production [17].
In this paper, we focus on the GENEVA framework, which has been developed in Refs. [2,3,18,19]. This method attains NNLO þ PS accuracy by combining fixed-order predictions at NNLO accuracy with the higher-order resummation of an N-jet resolution variable and matching the ensuing predictions with the PS. The N-jet resolution variable partitions the phase space into regions with a different number of resolved emissions in the final state, such that infrared (IR) divergent final states with M partons are translated into finite final states with N-partonic jets (where M ≥ N). This ensures that the IR divergences cancel on an event-by-event basis yielding physical and IR-finite events.
Although the GENEVA method is completely general to the point that several different resolution variables can be used, so far the only practical implementations employed N-jettiness T N [20] to achieve the separation between the N-and (N þ 1)-jet events. Nonetheless, any other resolution variable which can be resummed at high enough accuracy can be used. In this paper, we extend the GENEVA method using the transverse momentum q ⊥ of a colorsinglet system as a 0-jet separation variable. The choice of this variable is principally dictated by the availability of higher-order resummation up to the next-to-next-to-nextto-leading logarithm (N 3 LL) and by the extreme precision at which it is measured by the LHC experiments for different processes. However, the peculiar vectorial nature of this observable deserves further comment, especially in view of its usage as 0-jet resolution. Indeed, at variance with N-jettiness, it does not completely single out the soft and collinear limits when q ⊥ → 0. The reason is that there are two competing mechanisms which can drive q ⊥ to small values: The first is the usual ensemble of soft and/or collinear partonic emissions with k ⊥;i ≃ 0 recoiling against the produced color singlet; the second proceeds instead through a combination of two or more relatively hard emissions balancing each other, such that the resulting vectorial sum of their ⃗ k t;i is zero. It happens that the latter mechanism is the one that dominates at small q ⊥ and is responsible for the behavior of the differential cross section in the small q ⊥ limit: Indeed, in the low q ⊥ region the spectrum vanishes as dσ=dq ⊥ ∼ q ⊥ instead of vanishing exponentially as the Sudakov suppression of the first mechanism would suggest [21].
The presence of these two competing mechanisms seems to prevent the usage of q ⊥ as a 0-jet resolution observable. However, both effects which lead to a vanishing q ⊥ are properly included in the 0-jet cross section by performing the resummation of q ⊥ in the Fourier-conjugate impactparameter space [22]. The problem of resummation in direct space is more delicate, and only recently was it shown that it is possible to directly resum the transverse momentum spectrum in direct space without spoiling the scaling at low q ⊥ [23,24]. This problem was also addressed in Ref. [25] within a soft-collinear effective theory (SCET) approach, where the renormalization-group evolution is solved directly in momentum space. The remaining nonsingular contribution stemming from two or more relatively hard emissions resulting in a system with small q ⊥ is strongly suppressed and can be neglected. As a result, it is possible to use q ⊥ as a proper 0-jet resolution variable. In particular, in this work we consider the RADISH approach of Refs. [23,24], which allows us to resum the transverse momentum spectrum at N 3 LL.
Here we focus on neutral DY production at the LHC (pp → Z=γ Ã → l þ l − ) as a case study, but this approach can immediately be applied also to the charged current case and to other color-singlet processes. Differential distributions of electroweak gauge bosons play a paramount role in the precision program at the LHC. These observables are measured at the level of their leptonic decays and are typically characterized by particularly small experimental uncertainties, which can reach the few per mille level in neutral DY production [26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41]. This further motivates the inclusion of higher-order calculations in MC event generators for this process.
In this work, we improve the fully differential NNLO calculation with the N 3 LL resummation for the transverse momentum, and we subsequently shower and hadronize the events with PYTHIA8 [74], while maintaining NNLO accuracy for the underlying process and a numerically good agreement with the resummed q ⊥ spectrum. We compare our predictions with recent data collected at the LHC by the ATLAS [40] and CMS [41] Collaborations at a center-of-mass energy of 13 TeV.
The manuscript is organized as follows: In Sec. II we introduce the theoretical framework; we first review the GENEVA method by discussing its extension to a different 0-jet resolution variable and recap the resummation of the transverse momentum within the RADISH formalism. In Sec. III, we discuss the details of the implementation and present the validation of our results by performing comparisons against fixed-order predictions at NNLO. We also study the differences with respect to the original GENEVA implementation of Refs. [3,19], which uses the beam thrust as the 0-jet resolution variable. We then compare our predictions with the experimental data in Sec. IV and finally draw our conclusions in Sec. V.

II. THEORETICAL FRAMEWORK
In this section, we briefly review the theoretical framework used in this work. We start by providing a short description of the GENEVA method, focusing especially on the separation between 0-and 1-jet events. This is particularly relevant in our case due to the change in the 0-jet resolution variable. We then discuss the resummation of the transverse momentum spectrum in direct space using the RADISH formalism.

A. The GENEVA method
The GENEVA method is based on the definition of physical and IR-finite events at a given perturbative accuracy, with the condition that IR singularities cancel on an event-by-event basis. This is achieved by mapping IR-divergent final states with M partons into IR-finite final states with N jets, with M ≥ N. Events are classified according to the value of N-jet resolution variables r N which partition the phase space into different regions according to the number of resolved emissions. In particular, the GENEVA Monte Carlo cross section dσ MC N receives contributions from both N-parton events and M-parton events where the additional emission(s) are below the resolution cut r cut N used to separate resolved and unresolved emissions. The unphysical dependence on the boundaries of this partitioning procedure is removed by requiring that the resolution parameters are resummed at high enough accuracy. Considering events with up to two resolved emissions necessary to achieve NNLO accuracy for color-singlet production, one has where by using the notation r N >r cut N we indicate that the events are differential in the N-jet resolution variable, while with the notation r cut N , we indicate that the MC cross section contains events in which the resolution variable has been integrated up to the resolution cut. As mentioned in the Introduction, the definition of the partonic jet bins is based on the definition of a suitable phase-space mapping Φ N ðΦ M Þ, where N and M denote the number of jets and partons, respectively. This ensures the finiteness of the dσ MC N cross section; however, they must not be mistaken for the jet bins commonly used in experimental analyses, which instead define jets according to a particular jet algorithm. 1 Defining the events as in Eq. (1), the cross section for a generic observable X is where by M X ðΦ N Þ we indicate the measurement function used to compute the observable X for the N-jet final state Φ N . We stress that the expression σðXÞ is not exactly equivalent to the result of a fixed-order calculation at NNLO accuracy, as for unresolved emissions the observable is calculated on the projected phase space Φ N ðΦ M Þ rather than on Φ M . However, the difference is of nonsingular nature and vanishes in the limit r cut N → 0. The value of the resolution cut r cut N should therefore be chosen as small as possible. In this limit, the result contains large logarithms of the resolution variable r N (and of r cut N ), which one should resum to all orders in perturbation theory to yield meaningful results.
Let us start by discussing the separation between the 0jet and the 1-jet regimes and the associated resummation of the 0-jet resolution variable. The expressions for the MC cross sections in the exclusive 0-jet bin and in the inclusive 1-jet bin are respectively. The label "res" stands for "resummed" and the nonsingular contributions labeled "nons" must only contain terms which are nonsingular in the r 0 → 0 limit. When 1 In principle, nothing prevents the usage of a resolution variable based on a standard jet algorithm to define the GENEVA jet bins, e.g., p T;jet , if not for the lack of the corresponding higherorder logarithmic resummation. working at NNLO accuracy, to ensure that the second contribution on the right-hand side of Eqs. (3) and (4) is genuinely nonsingular, the resummed contribution must include all the terms singular in the resolution variable r 0 at order α 2 s ; i.e., it should be evaluated at least at NNLL 0 accuracy. In the context of hadron collider processes, the GENEVA implementations have so far only employed the beam thrust T 0 as 0-jet resolution variable, whose resummation was performed at NNLL 0 in SCET. In this work, we extend the GENEVA framework to use the transverse momentum of the color-singlet q ⊥ as 0-jet resolution parameter, which we resum at N 3 LL accuracy; see Sec. II B. It is worth emphasizing again that when one chooses q ⊥ as the 0-jet resolution, one is not only singling out configurations without any hard emission, as there exist configurations with two or more partons with transverse momenta balancing each other such that the vectorial sum of their k ⊥ is small. These contributions are included by the resummation in dσ MC 0 dΦ 0 ðq cut ⊥ Þ, but they are only described in the soft and/or collinear approximation, missing powersuppressed nonsingular corrections which are progressively important at increasing values of k ⊥ . However, as long as one keeps q cut ⊥ sufficiently small, the chance of obtaining a small vectorial sum from larger and larger individual transverse momenta is heavily suppressed, and therefore, the contribution to this region stemming from hard jets can be safely neglected. In Sec. III B, we will discuss the value of q cut ⊥ necessary to make the power-suppressed terms sufficiently small that they may be neglected.
In order to make the r 0 spectrum fully differential in the Φ 1 phase space, Eq. (4) contains a splitting function PðΦ 1 Þ fulfilling the normalization condition Z dΦ 1 This function is used to extend the differential dependence of dσ res including the full radiation phase space written as a function of r 0 and two additional variables, which we choose to be the energy fraction z and the azimuthal angle ϕ of the emission. The nonsingular contributions entering in Eqs. where NNLO 0 and NLO 1 indicate the accuracy at which one should compute the fixed-order contributions for each jet bin. The terms in square brackets are the expansion of the resummed result at Oðα 2 S Þ of the cumulant and of the spectrum.
By writing explicitly the expressions for the fixed-order cross sections, one has for the 0-jet bin and for the inclusive 1-jet bin.
In the above equations, B M are the M-parton tree-level contributions, V M are the M-parton one-loop contributions, and finally, W 0 is the two-loop contribution. In Eq. (9),we also introduced the notation to indicate that the integration is performed over a region of Φ 2 while keeping the values of r 0 and of the Φ 1 observables that we use to parametrize the 1-jet phase space fixed. This is required since the expression for dσ MC ≥1 is differential in the 0-jet resolution parameter r 0 , meaning that we should also parametrize the 2-parton contribution B 2 in the same terms. This means that the mapping used in the projection dΦ 2 =dΦ r 0 1 should preserve r 0 , i.e., to guarantee the pointwise cancellation of the singular r 0 contributions in Eq. (9). The Θ r 0 ðΦ 1 Þ in Eq. (10) limits the integration to the singular points for which it is possible to construct such a projection. The nonprojectable region of the Φ 2 phase space, which only includes nonsingular events, must therefore be included in the 2-jet event cross section, as shown later. We finally notice that Eq. (9) should be supplemented by including in dσ MC 1 1-jet events with r 0 < r cut 0 which cannot be mapped into 0-jet events. In the 1 → 0 case, we use the Frixione-Kunszt-Signer (FKS) mapping [75], which implies that the only nonprojectable events are those which would result in an invalid flavor configuration that is not present at the LO level (e.g., qg → l þ l − q projected to gg → l þ l − ). By denoting the nonprojectable region with the symbolΘ FKS map , we supplement the formulas above with Before we move on to discuss the N 3 LL resummation of the transverse momentum, we shall briefly review the separation of the 1-jet cross section into an exclusive 1-jet cross section and an inclusive 2-jet cross section. This separation can be performed in analogy to the 0-/1-jet separation discussed above, with r 1 now being the relevant resolution variable: In these equations, the resummation accuracy can be lowered with respect to the 0-/1-jet separation and is set to NLL. The full expressions for the exclusive 1-jet and the inclusive 2-jet cross section read [3,11] dσ MC In these formulas, U 1 ðΦ 1 ;r cut 1 Þ is the Sudakov form factor that resums the dependence on r cut 1 at NLL, U 0 1 ðΦ 1 ;r 1 Þ denotes its derivative with respect to r cut 1 , and U ð1Þ 1 and U ð1Þ0 1 are their Oðα S Þ expansions. The C 2 term is the singular approximant of the double-real matrix element B 2 , which acts as a standard NLO subtraction reproducing the singular behavior of B 2 . The term V C 1 is defined from the relation The normalized splitting function PðΦ 2 Þ satisfies the unitarity condition [cf. the 0-/1-jet case Eq. (5)] Analogous to the case of the normalized splitting function PðΦ 1 Þ, the function PðΦ 2 Þ extends the differential dependence of dσ res including the radiation phase space for two emissions, where the second emission should now be parametrized in terms of r 1 and of the same two additional radiation variables. Despite having used the same notation, we stress that the mapping used in Eqs. (16) and (18) does not need to correspond to the one used to implement the subtraction in Eq. (15). Indeed the latter does not need to be written as a function of r 1 .
The formulas above should be supplemented by including nonprojectable 2-jet events with r 1 <r cut 1 , namely, where we use the symbolΘ r 0 map ðΦ 2 Þ to denote the region of phase space which cannot be projected on the physical Φ 1 phase space via the Φ r 0 1 ðΦ 2 Þ mapping. In principle, in the r 0 <r cut 0 region also, events with two hard emissions contribute at NNLO. However, we stress again that when the 0-jet resolution variable is q ⊥ and one keeps q cut ⊥ small enough, these nonsingular contributions from events with two hard jets are negligible, and we therefore set This concludes our review of the GENEVA method.
In the following, we will use the transverse momentum as 0-jet resolution variables and 1-jettiness as a 1-jet resolution variable. We shall use the notation GENEVA q ⊥ to label our predictions, while we will use the notation GENEVA T 0 when referring to the original GENEVA implementation using beam thrust. In the next section, we discuss the resummation of the transverse momentum distribution in color-singlet production within the RADISH framework. We refer the reader to Refs. [11,19] for the discussion of the NLL resummation of T 1 within the GENEVA formalism.
The RADISH formalism allows one to resum any transverse observable (i.e., not depending on the rapidity of the radiation) which fulfills recursive infrared collinear (RIRC) safety [102]. The resummation is formulated directly in momentum space by exploiting factorization properties of squared QCD amplitudes. The resummation is then evaluated numerically via MC methods. The RADISH formulas are more conveniently expressed at the level of the cumulative cross section where q ⊥ ¼q ⊥ ðΦ 0 ;k 1 ; …;k n Þ is a function of the Born phase space Φ 0 of the produced color singlet and of the momenta k 1 ; …;k n of n real emissions. For example, in the soft limit, the all-order structure of Σðq cut ⊥ Þ fully differential in the Born phase space can be written as Here, ½dk i and dΦ 0 denote the phase spaces of the ith emission k i and of the Born configuration, respectively, M denotes the matrix element for n real emissions in the soft approximation, and VðΦ 0 Þ is the resummed form factor in the soft limit, which encodes the purely virtual corrections; see [103][104][105]. To obtain the resummation, it is first necessary to establish a well-defined logarithmic counting in the squared amplitude. The counting is established by decomposing the squared amplitude defined in Eq. (22) in n-particle-correlated blocks, which contain the correlated portion of the squared n-emission soft amplitude and its virtual corrections [24,102,106]. In particular, blocks with n particles start contributing one logarithmic order higher than blocks with (n − 1) particles, which allows one to systematically identify all the relevant contributions entering at a given logarithmic order.
Thanks to the RIRC safety of the observables, the divergences appearing at all perturbative orders in the real matrix elements cancel exactly those of virtual origin, which are contained in the VðΦ 0 Þ factor of Eq. (22). The cancellation of the singularities is achieved by introducing a resolution scale q 0 . Radiation softer than q 0 is dubbed unresolved and can be exponentiated to cancel the divergences of virtual origin. Radiation harder than q 0 (resolved) must instead be generated exclusively as it is constrained by the measurement function θðq cut ⊥ −q ⊥ ðΦ 0 ;k 1 ; …;k n ÞÞ in Eq. (22). Recursive IRC safety ensures that neglecting radiation softer than q 0 (unresolved) in the computation of q ⊥ ðΦ 0 ;k 1 ; …;k n Þ only produces terms suppressed by powers of q 0 , thus ensuring that the limit q 0 → 0 can be taken safely.
In the RADISH formalism, the resolution scale q 0 is set to a small fraction δ > 0 of the transverse momentum of the correlated block with the largest k ⊥ , which we henceforth denote as k ⊥;1 . The cumulative cross section at N 3 LL accuracy for the production of a color singlet of mass M fully differential in the Born kinematic variables and including also the effect of collinear radiation can then be written as [24] dσ res dΦ 0 ðq cut ⊥ Þ¼ In the equation above, the first line enters already at NLL, the first set of curly brackets (second to fifth line) starts contributing at NNLL, and the last set of curly brackets (from line six) enters at N 3 LL.
The functions L are luminosity factors evaluated at different orders which involve, besides the parton luminosities, the process-dependent squared Born amplitude and hard-virtual corrections H ðnÞ as well as the coefficient functions C ðnÞ ci , which have been evaluated to second order for qq-initiated processes in Refs. [88,89,91]. The factorsP ð0Þ denote the regularized splitting functions. The interested reader is referred to Sec. 4 of Ref. [24] for the definition of the luminosity factors and their ingredients. We defined ζ si ≡ k ⊥;si =k ⊥;1 and we introduced the notation dZ½fR 0 ;k i g to denote an ensemble describing the emission of n identical independent blocks. The average of a function GðΦ 0 ; fk i gÞ over the measure dZ as (ζ i ≡ k ⊥;i =k ⊥;1 ) is defined as Note that the ln 1=δ divergence appearing in the exponential prefactor of Eq. (24) cancels exactly that in the resolved real radiation encoded in the nested sums of products on the righthand side of Eq. (24).
Equation. (23) has been obtained by expanding all the ingredients around k ⊥;1 since ζ i ¼ k ⊥;i =k ⊥;1 ∼ Oð1Þ; see Ref. [24]. Thanks to RIRC safety, blocks with k ⊥;i ≪ k ⊥;1 are fully canceled by the term expf−R 0 ðk ⊥;1 Þ lnð1=δÞg of Eq. (24). Although such an expansion is not strictly necessary, it makes a numerical implementation much more efficient [24]. Because of this expansion, Eq. (23) contains explicitly the derivatives of the radiator R, which is given by is the resummation scale, whose variation is used to probe the size of missing logarithmic higher-order corrections in Eq. (23). The functions g i are reported in Eqs. (B.8)-(B.11) of Ref. [24]. The above formulas are implemented in the code RADISH, which evaluates them using MC methods. We refer the reader to Sec. 4.3 of Ref. [24] for the technical details.
The resummed result provided by RADISH is valid in the soft/collinear region, i.e., q ⊥ =M ≪ 1, and must be matched with fixed-order predictions at large values of q ⊥ .
Resummation effects should thus vanish in the large q ⊥ region. This is enforced by mapping the limit k ⊥;1 → Q, where the logarithms vanish onto k ⊥;1 → ∞ via modified logarithms where p is a positive real parameter whose value is chosen such that the resummed component decreases faster than the fixed-order spectrum for q ⊥ =M ≳ 1. In the following, we took p ¼ 4. Therefore, the logarithms L in the Sudakov radiator Eq. (26), its derivatives Eq. (25), and the luminosity factors have to be replaced byL. For consistency, one must supplement Eq. (23) with the following Jacobian: This prescription leaves the measurement functions in Eq. (23) unchanged. The final result is modified beyond the nominal accuracy by power corrections in ðQ=k ⊥;1 Þ p .

III. IMPLEMENTATION DETAILS
In this section, we discuss the interface of the RADISH code with GENEVA and the implementation of the DY process in GENEVA using q ⊥ as a 0-jet resolution variable. We validate our framework by comparing our results with NNLO predictions obtained with MATRIX [107] and discuss the interface with PSs.

A. Interfacing GENEVA with RADISH
The resummation formulas discussed in Sec. II B are implemented in the FORTRAN90 code RADISH. For each Born event, the code produces the initial-state radiation and performs the resummation of large logarithmic contributions using MC methods as described in Sec. 4.3 of Ref. [24]. As a result, the cumulative distribution is filled on the fly yielding dσ res =dΦ 0 ðq ⊥ cut Þ. The spectrum can be calculated by differentiating the cumulative histogram numerically.
We note that our need to obtain the spectrum at a given value of q ⊥ as in Eq. (4) poses some technical challenges which need to be addressed in order to interface RADISH with GENEVA. Indeed, the RADISH code has been designed to compute the whole cumulant distribution by generating a MC ensemble of soft and/or collinear emissions. The simplest solution would be to run the MC algorithm with a very large number of events for each Born configuration provided by GENEVA in order to yield sufficiently stable results for the numerical derivative of the cumulant. However, despite the effectiveness of the MC generation, this on-the-fly computation is rather inefficient.
A better approach is to first compute the RADISH cumulant with a very large number of Born points, exploiting run-time parallelization. This is done before starting the main GENEVA runs. From these parallel RADISH runs, we build an interpolation grid, which is then used to provide the cumulant and the spectrum for the GENEVA runs on an event-by-event basis. To this end, we first parametrize the Born phase space using two variables 2 and construct a discrete lattice in these two parameters. We then compute the resummed cumulant in each bin and for each combination of perturbative scales up to the maximal value of q ⊥ kinematically allowed. The resulting grids are then loaded and interpolated on the fly by means of Chebyshev polynomials as implemented in the GSL library [108], yielding a fast evaluation of the cumulant and of its derivative. Since the shape of the resummed result depends only mildly on the Born variables, it is sufficient to rescale the resummation by the Born squared amplitude to obtain a result differential in the Born kinematics. We found that a (10 × 10)-dimensional grid in the Born variables provides a sufficiently fine discretization of the Born phase space.
The other ingredient which is provided by RADISH is the expansion of the resummation to remove the double counting between resummation and fixed order. In this case, the use of interpolation routines is not advisable, since one must ensure the exact cancellation of the large logarithmic terms at low q ⊥ to avoid spurious and potentially large effects on the final results. 3 However, since the expansion is computed semianalytically in RADISH (for additional details, we refer the reader to Sec. 4.2 of [24]), we can avoid the interpolation altogether and compute it on the fly without affecting the speed of the code.
Once the values of the resummation and the expansion are available, the matched results for the spectrum of Eq. (4) can be obtained by means of a standard additive matching. The use of modified logarithms (27) automatically ensures that the effect of the resummation vanishes in the fixed-order region. However, numerical instabilities could arise since, while the expansion is exact, the resummed result is obtained by means of a Chebyshev interpolation. This can induce tiny but visible numerical artifacts in the region where the cumulant flattens and its derivative approaches zero. In order to avoid this unwanted feature, we introduce a function which further suppresses both the resummed and the resummed expanded results smoothly in the fixed-order region, where a complete cancellation between them should be achieved. In this way, undesired numerical instabilities in this region are removed. In particular, we use the following function: 2 Specifically, we use the rapidity of the virtual boson Y ll and a ll ¼ arctanððm 2 ll − m 2 Z Þ=ðm Z Γ Z ÞÞ for DY production. The latter variable is chosen to flatten the m ll distribution; see, e.g., Ref. [9].
with m ¼ðr þ lÞ=2, w ¼ðr − lÞ=2, and where we take l ¼ 2Q and r ¼ 3Q and "erf" is the standard Gaussian error function. We observe that with this choice there is a very tiny discontinuity at q ⊥ ¼ r, which, however, does not give any visible effect since f damp ðr − ϵÞ ≃ 0. The damping function is plotted as a function of q ⊥ =Q in Fig. 1.

B. Power-suppressed corrections to the nonsingular cumulant
In Eq. (8), we wrote the expression for the NNLO accurate 0-jet cross section fully differential in the Born phase space, whose implementation would require a local NNLO subtraction scheme. In our implementation, the NNLO 0 accuracy is instead achieved replacing Eq. (8) with the following expression: which only involves a local subtraction at Oðα S Þ and the expansion of the resummation at the same order. The formula assumes an exact cancellation between the fixedorder and the resummed expanded contribution below the value of r cut 0 at order α 2 S . The cancellation is guaranteed for the singular contributions due to the accuracy of the r 0 cumulant (we stress again that in order to achieve NNLO 0 accuracy, the resummation accuracy must be at least NNLL 0 ); however, the formula fails to capture nonsingular contributions at Oðα 2 S Þ. These nonsingular contributions can be expressed as and their integral over the phase space is The nonsingular cumulant vanishes in the limit r cut 0 → 0 since the functions f i ðr cut 0 ; Φ 0 Þ contain at worst logarithmic divergences. As discussed above, our calculation includes the term f 1 ðr cut 0 ; Φ 0 Þ since we implement a NLO 1 FKS local subtraction. On the contrary, f 2 ðr cut 0 ; Φ 0 Þ is not included in Eq. (31). Neglecting the Oðα 2 S Þ power corrections is acceptable as long as we choose r cut 0 to be very small. So far, the GENEVA method has been based on N-jettiness subtraction [109][110][111]. In this work, we take r 0 equal to q ⊥ ; effectively, this corresponds to basing GENEVA on a q T subtraction scheme [112]. The availability of different resolution parameters in GENEVA is beneficial, as the size and the scaling of the power corrections can be different. The size of the missing Oðα 2 S Þ power corrections as a function of q cut ⊥ can be calculated as where δσ NNLO ¼ σ NNLO − σ NLO , Σ N 3 LL asymp is the cumulant of the resummed contribution [see Eq. (21)] and q ⊥ max is the maximum value allowed by the kinematics for each Φ 0 point. Finally, j α n S indicates the expansion up to order α n S . We have calculated δσ NNLO by computing σ NNLO with MATRIX. Note that MATRIX also achieves NNLO accuracy via q T subtraction, and therefore potentially misses power corrections at Oðα 2 S Þ.H o w e v e r ,MATRIX includes an estimate of these power corrections by interpolating the result to q cut ⊥ ¼ 0, and including an estimate of this interpolation procedure in its error.
We show the size of the missing power corrections in Fig. 2, where we consider values of q cut ⊥ down to 1 GeV. We consider pp → l þ l − production at 13 TeV in an inclusive setup by applying a cut only on the invariant mass of the produced color singlet. We observe that the power corrections are below the 2‰ level for log 10 ðq cut ⊥ Þ ≲ 0.5, corresponding to a value of q cut ⊥ ≲ 3 GeV, in accordance with what was observed in Ref. [107]. Motivated by this plot, we choose q cut ⊥ ¼ 1 GeV as our default value. The negligible size of the missing power correction allows us to avoid the need for the reweighting of the Φ 0 events, which was instead the approach followed in the previous application of the GENEVA method. Our choice is further justified by a detailed comparison between GENEVA and an independent NNLO calculation for distributions differential in the Φ 0 variables, as we will show in Sec. III E.

C. NLO 1 calculation and phase-space mapping
In this section, we discuss the implementation of the Φ r 1 ðΦ 2 Þ mapping introduced in Eq. (15) and of that introduced in Eq. (18). As we anticipated in Sec. II A, these phase-space mappings do not necessarily need to coincide, since only the latter needs to be written as a function of the 1-jet resolution variable T 1 .
We start by discussing the mapping used to implement the NLO 1 calculation. Let us first notice that the Φ r 1 ðΦ 2 Þ mapping used in the B 2 term and the Φ C 2 ðΦ 1 Þ used for the C term in Eq. (15) can be different, but provided that both are IR safe, they must be equivalent in the IR singular limit. However, as we discussed in Sec. II A, the Φ r 1 ðΦ 2 Þ mapping must have the additional property that it preserves the 0-jet resolution variable q ⊥ in order to guarantee that the NLO 1 is pointwise consistent with the q ⊥ resummation. Moreover, the mapping must be invertible, meaning that one should be able to reconstruct all the Φ 2 points which can be projected on a given Φ 1 configuration. The additional requirement that the mapping should preserve q ⊥ poses some challenges as discussed in Sec. A6 of Ref. [3]. In addition to preserving the 0-jettiness T 0 which was used as 0-jet resolution variable, 4 the T FR 0 mapping introduced in Ref. [3] also preserves the four-momentum of the vector boson q μ ðΦ 2 Þ¼q μ ðΦ T 1 ðΦ 2 ÞÞ and thus the full transverse momentum ⃗ q ⊥ . Therefore, in our study one could in principle use the same mapping to implement the NLO 1 calculation and the splitting function PðΦ 2 Þ used in the T 1 resummation. It happens, however, that the T FR 0 mapping creates undesired higher-order artifacts in a few NLO 1 differential distributions when used also in the T 1 resummation. In particular, one observes the effects on quantities such as the rapidity separation between the leading jet and the vector boson. For this reason, in all recent applications of the GENEVA method based on the resummation of T 0 , the T FR 0 mapping has been modified, relaxing the condition that it preserves the transverse momentum of the vector boson. Therefore, while one could use the ⃗ q ⊥ -preserving T FR 0 mapping for the subtraction, one needs a new phasespace mapping which preserves ⃗ q ⊥ for the T 1 resummation projection. For the initial-state radiation (ISR), we found that the q μ -preserving mapping introduced in Ref. [113] does not create distortions in any differential distribution, and we have therefore implemented this mapping in the GENEVA code. For final-state radiation (FSR), a mapping preserving the four-momentum q μ of the vector boson is uniquely defined once one decides how to map the (massive) jet with momentum p μ 12 ¼ p μ 1 þ p μ 2 of the 2-jet configuration into a massless jet with momentump μ in the 1-jet configuration. In particular, we choose to conserve the longitudinal component of the jet, i.e., p z 1 þ p z 2 ¼p z ,a s described in Appendix.
We have verified that these mappings can be used both in the NLO 1 calculation [in Eq. (15)] as well as in the T 1 resummation [in Eq. (16)]. We have, however, preferred to maintain the (⃗ q ⊥ -preserving) T FR 0 mapping in the NLO 1 calculation due to slightly better convergence properties, while we resort to the new ISR and FSR mappings in the implementation of the T 1 resummation. In this way, we avoid unnecessarily large higher-order effects in the distributions.

D. Validation of the N 3 LL resummation
In this section, we validate the implementation of the N 3 LL resummation in GENEVA by comparing the matched q ⊥ spectrum with an independent calculation obtained using the MATRIX + RADISH interface [114]. We note that the approach taken in the GENEVA framework is somewhat different from that taken in MATRIX + RADISH. In the GENEVA approach, each event generated is assigned a resummed weight according to the value of the chosen resolution parameter; as a consequence, the information about the physical event, which is characterized by a particular value of q ⊥ , is always retained throughout the calculation. As shown in Ref. [100], by retaining the exact lepton kinematics the GENEVA q ⊥ result effectively resums also the linear power corrections which appear as soon as one imposes fiducial cuts [115] beyond the Oðα 2 S Þ, which is in any case included after the matching to NNLO. This is FIG. 2. Nonsingular cumulant at order α 2 S as a function of q cut ⊥ . 4 More precisely, the mapping preserves a recursive definition of T 0 dubbed T FR 0 ; see Ref. [3] for additional details.
not the case in the MATRIX + RADISH calculation, where the matching is instead performed at the level of the final distribution, and therefore the linear power corrections associated with fiducial cuts are included only up to order Oðα 2 S Þ. 5 The predictions for the q ⊥ spectrum obtained in the two approaches are equivalent in the absence of fiducial cuts on the Born-level variables (see, e.g., the discussion in Sec. 4.3 of [15]). In the presence of cuts, the two approaches are equivalent only in the limit q ⊥ → 0, unless the quantity subject to the cuts is preserved by the phase-space mapping Φ r N ðΦ Nþ1 Þ. For the following validation, the only processdefining cut is applied to the invariant mass of the resulting color singlet, which is preserved by our mappings. Therefore, we expect full agreement above q cut ⊥ between our results and those obtained with MATRIX + RADISH, which we can use to validate our implementation.
To produce our predictions we use the NNLO NNPDF3.1 parton distribution function set [116] with α S ðM Z Þ¼0.118 through the LHAPDF interface [117].W e set the central renormalization and factorization scales to The uncertainty band is constructed as the envelope of a canonical seven-point variation of the renormalization (μ R ) and factorization (μ F ) scales and of two additional variations of Q by a factor of 2 for central μ R and μ F . The comparison is presented in Fig. 3, where we show the resummed and the matched results for the q ⊥ ≡ p ll ⊥ spectrum up to 100 GeV, on a scale which is linear up to 20 GeV and logarithmic for larger values. In both cases, we observe a very good agreement between the two calculations, both for the central prediction and for the theory uncertainty bands. We observe marginal deviations at large q ⊥ in the resummed result, which can be traced back to the additional damping present in the GENEVA implementation as discussed in Sec. III A. As expected, this difference cancels in the matched result. The difference in the first bin of the matched result, on the other hand, originates from the sensitivity to the different IR cutoffs employed in the two calculations necessary to regularize the q ⊥ → 0 limit in the NLO 1 calculation. We observe that after matching, the uncertainty bands are significantly reduced between 10 and 25 GeV and are at the 1%-2% level, which might not reflect the actual size of missing higher-order effects. However, since we are not including further sources of uncertainty such as parton distribution function (PDF) errors, mass effects, etc., one should not regard this as the full theoretical uncertainty. Including all these additional effects in order to provide a fully realistic error estimate is beyond the scope of this work.

E. Comparison with NNLO predictions and validation
In this section, we proceed with the validation of the GENEVA implementation by comparing our results with those obtained with an independent NNLO calculation using MATRIX. In Fig. 4, we compare four distributions already defined at the Born level, which allows us to assess the agreement of the NNLO corrections. In particular, we compare the invariant mass m ll and the rapidity y ll of the dilepton system, the rapidity of the hardest lepton y l;hardest , and the transverse momentum of the hardest lepton p T;l;hardest . We use the same settings as specified in the previous section, but we now use a canonical seven-point scale uncertainty prescription by keeping the resummation scale fixed to its central value Q. For the first three observables we find an excellent agreement between GENEVA and MATRIX, with differences at the few percent level, which is within the size of the statistical fluctuations. In the case of p T;l;hardest , we expect possible differences because the distribution is NNLO accurate only below the Jacobian peak at m ll =2. The presence of resummation effects in the GENEVA results improves the description of this observable above the peak with respect to the pure fixed-order result, smearing out the unphysical behavior of the distribution around the Sudakov shoulder [118]. Indeed, we observe that the two predictions are in good agreement up to ∼40 GeV. Above this value, the two predictions start to differ, and the two uncertainty bands barely overlap between 50 and 70 GeV. At larger values of p T;l;hardest , the two predictions approach each other again.
In conclusion, the comparison between the GENEVA and the MATRIX predictions provides a robust check of our implementation. In particular, the agreement observed at the differential level fully justifies the value of the q cut ⊥ of 1 GeV which we shall use to produce our results.

F. Interface with the PYTHIA8 parton shower
In order to make a calculation fully differential at higher multiplicities the GENEVA, partonic predictions must be matched to a PS. The shower adds extra radiation to the exclusive 0-and 1-jet cross section and extends the inclusive 2-jet cross section by including higher jet multiplicities. For the sake of definiteness in this study, as in the previous ones, we focus on the PYTHIA8 shower. The GENEVA interface to PYTHIA8 for hadronic processes is discussed in detail in Ref. [3]; here we briefly summarize the most relevant features, and we discuss the modifications to the interface needed to accommodate the change in the 0-jet resolution variable. In the previous GENEVA implementations based on T 0 resummation, the shower constraints aim to preserve as far as possible the NNLL 0 accuracy of the T 0 distribution. In particular, it was shown that for the majority of the events the first emission of the shower produces a distortion of the T 0 distribution at the Oðα 3 S =T 0 Þ level, i.e., beyond NNLL 0 [3]. However, due to the vectorial nature of q ⊥ , one cannot readily apply the same argument in this case. In general, one expects that any emission of the shower could alter the transverse momentum distribution of the color-singlet system, and ultimately, the logarithmic accuracy for the transverse momentum spectrum after the shower is therefore dictated by the shower accuracy. Hence, in this work we refrain from making any claims about the formal accuracy of the predictions for the q ⊥ spectrum after the showering. However, we will show below that it possible, with a suitable choice of the shower recoil scheme, to obtain an excellent numerical agreement between the analytic N 3 LL resummation and the GENEVA showered results. Lastly, one can get an excellent description of the data at small q ⊥ by tuning the PYTHIA8 nonperturbative parameters. However, in any calculation obtained by matching higher-order calculations with PS, one has to carefully evaluate which parameters are truly encoding nonperturbative effects and should therefore be tuned.
In order to discuss the details of the matching, let us start by analyzing the interface to the PYTHIA8 k ⊥ -ordered PS, starting with the 0-jet event case. For this set of events, the shower should simply restore the emissions which were integrated over in the construction of the 0-jet cross section. In our implementation, we set the starting scale to these events to the natural scale q cut ⊥ ; to avoid double counting with events above the cut, we require that after the shower the transverse momentum of the boson does not exceed q cut ⊥ , though we allow for a small spillover if the showered event has q ⊥ > 1.05 × q cut ⊥ to avoid an hard cutoff. Events which do not fulfill this constraint are reshowered. In practice, this spillover has a negligible effect, since 0-jet events account for Oð1%Þ or less of the total cross section, and are therefore a very small fraction of the total; moreover, since the starting scale for the showering of these events is ∼1 GeV; the majority of them automatically satisfy the constraint.
The showering of 1-and 2-jet events is more delicate. As discussed in Sec. II A, the separation between 1-and 2-jet events is achieved by means of a Sudakov form factor U 1 ðΦ 1 ; T 1 Þ, which suppresses the 1-jet cross section at small values of T 1 . The mapping used preserves the 0-jet resolution parameter r 0 ,i . e . ,T 0 in the original GENEVA implementation and q ⊥ in this work. Even with T cut 1 ¼ 1 GeV, σ MC 1 remains sizeable and must be handled with certain care. In the original GENEVA implementation, it was necessary to reduce its size to avoid further distortions of the T 0 distribution after the shower. Here we prefer to follow the same approach, which also helps to reduce the size of the nonsingular contributions in Eq. (19). Properly, to fully account for configurations where q ⊥ ∼ T 1 ≪ Q, one should extend the resummation framework used here and include a joint resummation for the two jet resolutions. Unfortunately, this is not yet available at the required logarithmic order. 6 In the T 1 ≪ q ⊥ limit, it is however possible to approximate the correct behavior and suppress the size of the 1-jet cross section σ MC 1 by multiplying it by an additional LL Sudakov form factor U 1 ðT cut 1 ; Λ 1 Þ;s e eR e f . [11].T h en e w1 -j e t exclusive cross section becomes while the 2-jet inclusive reads The parameter Λ 1 is set to be much smaller than T cut 1 and constitutes the ultimate 1-jet resolution cutoff (in the present work, we set this cutoff to 0.0001 GeV). The 1-jet cross section is therefore extremely suppressed and accounts for a negligible fraction of the cross section, at the few per mille level.
As a result, the inclusive 2-jet cross section accounts for almost all events. The starting scale of the shower for this set of events is chosen to be sufficiently large to completely fill the phase space available to the k ⊥ -ordered shower (in particular, in the present work we choose the largest relative k ⊥ with respect to the direction of the sister parton for the second emission). The events are then showered, vetoing events which have a T 2 value larger than the initial T 1 value once showering is complete. Events which do not pass this veto are reshowered again until they do. This in practice corresponds to the implementation of a truncatedvetoed shower [123] and thus preserves the LL accuracy of the shower for observables other than q ⊥ .
Having discussed the matching conditions of the shower, we now move to discuss the results. We match the NNLO calculation to the PYTHIA8.245 PS, and we use the generalpurpose CMS MonashStar Tune (Tune:pp = 18). As previously discussed, it would be desirable to maintain an agreement between the precise parton level predictions and the showered result. This is possible by choosing a more local recoil scheme in PYTHIA through the option SpaceShower:dipoleRecoil = on, which affects the color-singlet kinematics less compared to the standard recoil scheme [7,124].
This scheme does not change the transverse momentum of the colorless system when there is an emission involving an initial-final dipole, and it was shown to improve the agreement between parton level and showered predictions for the transverse momentum of a diphoton pair in our previous work [15]. We stress that the choice of the recoil scheme has important implications on the accuracy of PSs [125][126][127][128][129][130] and might induce spurious effects at NLL. Nonetheless, since in this work we assume that the accuracy of our predictions is dictated by the PS, we leave a study of the formal accuracy to future work.
Nonperturbative effects are another possible source of modification of the q ⊥ distribution induced by the matching with a PS. In general, the values of the nonperturbative parameters in PYTHIA8 have been tuned starting from predictions which were lacking the higher-order effects that are instead included in this work. Therefore, it is not clear whether their values reflect real nonperturbative effects or simply the lack of perturbative ingredients. When matching to a NNLO computation, it would be advisable to investigate this and eventually perform a retuning. While the determination of an optimal tune is beyond the scope of this work, we observed that the description of the spectrum in the small q ⊥ region strongly depends on the value of the intrinsic (or primordial) k ⊥ transverse momentum of the incoming partons determined by the PYTHIA option BeamRemnants:primordialKThard.I nt h eC M S MonashStar tune, the value for this parameter is set to 1.8 GeV; nonetheless, we observed that this value is responsible for a shift in the spectrum outside scale variation bands up to values of q ⊥ ∼ 10 GeV, i.e., in a region where one would naively expect nonperturbative effects to play a minor role. Therefore, we preferred to lower this value to ∼0.6 GeV, such that the uncertainty bands before and after the shower overlap. In the following, we shall use this value as our default.
The effect of changing the recoil scheme is shown in Fig. 5 for the rapidity distribution of the dilepton pair and for the transverse momentum of the vector boson. In the plots, we compare the parton level predictions with the results after the shower with different recoil schemes, in the absence of hadronization and multiparticle interactions (MPI) effects. The difference in the rapidity distribution is negligible almost everywhere: The dipole recoil scheme is indistinguishable from the parton level results, while there are tiny discrepancies for the default recoil scheme, but these appear only at very large rapidities. On the other hand, the effect on the transverse momentum distribution is somewhat more pronounced, although the two schemes are in good agreement within the scale uncertainty bands, with differences at the few percent level. In particular, when the more local recoil scheme is chosen, the agreement between the showered result and the N 3 LL þ NLO 1 result at the parton level is excellent across the whole spectrum.
We now compare the partonic, showered, and hadronized result for a selected set of distributions. In order to keep the analysis simple, we do not include QED effects in the showered results. In the following, we always include MPI effects in our hadronized results. We start by presenting this comparison in Fig. 6 for the rapidity distribution of the lepton pair and for the transverse momentum of the hardest lepton. In the left panel, we observe that for the rapidity distribution, the NNLO accuracy is maintained at a very precise level after the showering and the hadronization stages, as should be expected by the inclusive nature of this distribution. In the right panel, the same excellent agreement is also seen in the case of the p T;l;hardest distribution case, with differences well inside statistical fluctuations.
Next, in Fig. 7 we focus on the transverse momentum distribution of the lepton pair, comparing the partonic result to the result after showering (upper row) and after hadronization and MPI (lower row). We present three separate FIG. 5. Comparison between the parton level and the showered results with and without the PYTHIA8 DIPOLERECOIL option for the rapidity distribution (left) and for the transverse momentum distribution (right) of the lepton pair. plots, in the peak, transition, as well as in the tail region of the distribution. An extremely good agreement between the showered and the partonic results is observed in all three regions, both for the central predictions and for the scale uncertainty bands. When the hadronization (and MPI) stage is also added, some differences arise, localized in the small q ⊥ region. There we observe that the peak becomes more pronounced, and the spectrum is more suppressed at small  7. Comparison of the q ⊥ spectra between the partonic N 3 LL þ NLO 1 and the showered results, after interfacing to PYTHIA8, before the inclusion of nonperturbative effects (above), and after hadronization and MPI (below). The peak (left), transition (center), and tail (right) regions are shown. q ⊥ , while at large q ⊥ the showered and hadronized results are again in agreement, as one expects from factorization.

G. Comparison of resolution variables
The implementation of two different 0-jet resolution variables opens up the possibility of comparing the predictions after shower, hadronization, and MPI obtained using GENEVA in the two cases. Moreover, we can also compare the results to the parton level predictions, where the higher-order resummation is guaranteed by construction for each 0-jet resolution variable. We show this comparison in Figs. 8 and 9 for the q ⊥ and T 0 spectra, respectively. There, the N 3 LL þ NLO 1 (NNLL 0 þ NLO 1 ) results at the parton level for each 0-jet resolution variable are compared to the results after showering and after hadronization obtained using q ⊥ and T 0 in GENEVA.F o r the transverse momentum distribution, we observe differences at the 10%-20% level between the results obtained with GENEVA q ⊥ and GENEVA T 0 after the shower up to values of q ⊥ close to 30 GeV, though the uncertainty bands always overlap. The GENEVA T 0 result is more suppressed at small q ⊥ and slightly harder above the peak. The differences are reduced at larger values of q ⊥ , although the GENEVA T 0 spectrum is somewhat harder than that of GENEVA q ⊥ . A similar pattern is visible after hadronization and MPI.
Analogous differences can be observed when comparing the T 0 spectra obtained using the two different implementations. The showered results are in good agreement for values of T 0 larger than 30 GeV, while they start to differ FIG. 8. Comparison between parton level and showered results (left) and between parton level and results after shower, hadronization, and MPI (right) for the transverse momentum distribution of the dilepton pair using GENEVA q ⊥ and GENEVA T 0 . below. The differences are as large as 50% at small values of T 0 , where the GENEVA q ⊥ result is considerably harder than the NNLL 0 þ NLO 1 result. The hadronization and especially the addition of MPI effects significantly modify the T 0 distribution, and bring the two results closer also at small values of T 0 . To accommodate this larger difference, instead of the usual ratio with respect to the partonic result, in the right panel of Fig. 9 we show the ratio with respect to the GENEVA T 0 result after hadronization and MPI.
Finally, we show the same comparison for the ϕ Ã variable [131] in Fig. 10. Although this observable is not fully resummed at N 3 LL accuracy in the GENEVA q ⊥ implementation, we expect to observe good agreement between the GENEVA q ⊥ and the resummed results since the two observables are closely related. Indeed, we observe that after showering, the GENEVA q ⊥ result is close to the N 3 LL þ NLO 1 result obtained with MATRIX + RADISH. The GENEVA T 0 result displays instead differences similar to those seen for the q ⊥ distribution, both before and after hadronization.

IV. RESULTS AND COMPARISON TO LHC DATA
In this section, we compare our predictions against 13 TeV data collected at the LHC by the ATLAS [40] and the CMS [41] experiments. We have generated events using the same settings as detailed in Sec. III D regarding the choice of PDF and central scales, and we use the shower settings as described in Sec. III F. We remind the reader that our predictions thus include hadron decay and MPI effects, but do not include QED shower effects.

A. Comparison to ATLAS data
We start by showing the comparison of our predictions with the ATLAS data of Ref. [40]. The Z boson is reconstructed by selecting the two hardest same-flavor opposite-sign leptons in the final state. We then apply the following cuts to our events: In Fig. 11, we compare our predictions for the normalized q ⊥ and ϕ Ã distributions. For q ⊥ , we use a scale which is linear up to 30 GeV and logarithmic for larger values. The former are in very good agreement with the data in the whole q ⊥ range. Below 30 GeV, the central prediction is within a few percent of the data, and only in the first two bins, where hadronization and nonperturbative effects play a prominent role, do the scale uncertainty bands fail to cover the experimental data. Our predictions are also in good agreement with the ϕ Ã measurements, matching them within scale uncertainty bands down to values of ϕ Ã ∼ 0.01; at lower values, the differences reach the 20% level in the first bin, and the perturbative uncertainty does not cover the data. Here, the inclusion of shower and nonperturbative uncertainties as well as the development of a dedicated tuning could help ameliorate the agreement. Finally, in Fig. 12 we can compare our predictions with parton level results at N 3 LL þ NNLO 1 accuracy [97] for the q ⊥ and ϕ Ã distributions obtained by matching RADISH results with fixed-order predictions from NNLOJET [132][133][134] in the ATLAS fiducial region. In the q ⊥ case, we also show the GENEVA parton level predictions at N 3 LL þ NLO 1 accuracy for reference. Our results for the q ⊥ spectrum are in good agreement with the parton level predictions across the whole range. In particular, the results are within a few percent from the N 3 LL þ NNLO 1 result up to 30 GeV; at larger values, where NNLO 1 corrections become dominant, the differences can reach Oð10%Þ. We observe a similar agreement in the ϕ Ã FIG. 10. Comparison between the N 3 LL þ NLO 1 prediction and showered results (left) and between the N 3 LL þ NLO 1 prediction and results after hadronization and MPI (right) for the ϕ Ã distribution using GENEVA q ⊥ and GENEVA T 0 . case, with differences at most at the 10% level around ϕ Ã ∼ 0.01.

B. Comparison to CMS data
Finally, we compare our predictions with the CMS data of Ref. [41]. The fiducial region in this case is defined by the selection cuts jm ll − 91.1876 GeVj < 15 GeV: We compare our predictions both at the absolute level and for normalized distributions. We note that in the analysis of Ref. [41], the latter are defined by dividing by the sum of the weights rather than by the integrated cross section; i.e., it is not normalized by the fiducial cross section. In Fig. 13, we show the comparison for the absolute rapidity jy ll j and for the q ⊥ and ϕ Ã distributions. In the first row, we observe a reasonably good agreement for the rapidity distribution. The theoretical predictions are a few percent higher than the experimental data at the absolute level, and oscillate around the data in the normalized distribution. In the central row, the predictions for q ⊥ are also in good agreement with the CMS data, repeating the same pattern already observed when comparing to ATLAS data. This is especially true at the normalized level.
In the last row, we observe that the theoretical predictions display statistical fluctuations more pronounced than in the ATLAS case. This is due to the very fine binning at small FIG. 12. Comparison between GENEVA results and N 3 LL þ NNLO 1 predictions for the transverse momentum distribution (left) and for the ϕ Ã observable (right) within ATLAS fiducial cuts. values of ϕ Ã . Moreover, the normalized predictions seem to consistently undershoot the data. A possible explanation for this effect could lie in the interplay between the statistical fluctuations at small values of ϕ Ã and the atypical normalization chosen for this particular analysis, which enhances the impact of the low ϕ Ã region across the whole spectrum. This can be further appreciated by noticing that the relative size of the uncertainty bands is somewhat larger in the normalized ϕ Ã distribution, defying a naive expectation. Finally, in Fig. 14 we compare our predictions for the transverse momentum distribution in different rapidity slices with the CMS data. In all cases, we observe a good agreement, with larger differences localized in the first few bins where hadronization effects are more significant.

V. CONCLUSION
In this work, we have presented an extension of the GENEVA method which uses the transverse momentum of the color singlet as the 0-jet resolution variable. As a first application, we have provided an NNLO þ PS description of DY pair production by matching the parton level calculation with the PYTHIA8 PS. The method presented here is fully general and can be readily applied to other color-singlet processes.
In this specific implementation, the transverse momentum resummation is performed at N 3 LL accuracy by interfacing GENEVA with the RADISH code. We have validated our predictions of the transverse momentum resummation by comparing with the N 3 LL þ NNLO 0 results obtained with the MATRIX + RADISH interface. The use of the transverse momentum as the 0-jet resolution variable has allowed us to reduce the impact of missing power corrections in the NNLO calculation. Setting q cut ⊥ ¼ 1 GeV, we have found that the missing power corrections contribute below the per mille level. We have validated our NNLO differential predictions by comparing a few selected distributions against the MATRIX program.
The availability of a fully exclusive event generator at this accuracy allows for the evaluation of any fiducial cross section, also correctly resumming linear power corrections associated with cuts on the lepton kinematics.
We have then studied the impact of PS and hadronization on our predictions and have found that the NNLO description of inclusive observables is preserved after both shower and hadronization. An important outcome of this study is the observation that it is possible to maintain an extremely good agreement between the N 3 LL þ NLO 1 result for the transverse momentum spectrum and the GENEVA results after the shower by choosing a more local recoil scheme in PYTHIA8. We have also quantified the impact of nonperturbative corrections after hadronization on the transverse momentum distribution, finding, as expected, that they are localized below the peak.
The construction of a NNLO þ PS event generator is subject to several assumptions, and the choice of which resolution variables are used is of particular importance. The availability of two different 0-jet resolution variables within the GENEVA framework has allowed us to robustly assess the impact of such choices on differential observables for the first time. We have studied the effect of the FIG. 14. Comparison between GENEVA and the CMS data for the q ⊥ distribution in different rapidity slices.
shower and hadronization (including MPI effects) on the higher-logarithmic partonic predictions using both T 0 and q ⊥ as the 0-jet resolution variable. The major differences are observed below the peak region, where the higher-order resummation of the correct logarithms is necessary to reproduce the correct behavior. Nonetheless, the two predictions are in reasonable agreement, with the size of the differences never exceeding 15% for the q ⊥ and ϕ Ã distributions.
Finally, we have compared our predictions, including hadronization and MPI effects, to LHC data at 13 TeV. A good agreement has been observed both for inclusive and for more exclusive distributions, such as the transverse momentum spectrum of the dilepton system as well as the ϕ Ã variable, across the whole spectrum. The description at very small q ⊥ and ϕ Ã is sensitive to the details of the hadronization procedure, suggesting that dedicated studies are needed in order to determine the optimal value of the nonperturbative parameters in a NNLO þ PS computation. It should be emphasised, however, that at this level of precision, the inclusion of EW corrections and a careful study of the impact of mass effects become relevant. We have refrained from including these effects in this study, leaving them to future work. Another direction in which our predictions could be improved concerns the inclusion of higher-order matrix elements, beyond what is available in a NNLO calculation for DY with no additional jets. This work could be expanded even further by extending the study to the use of alternative 1-jet resolution variables and their resummation, as well as the interplay between the 0-jet and 1-jet resummations. In particular, in view of the recent advancements in the development of PSs beyond LL accuracy, the choice of a suitable 1-jet resolution variable might prove crucial in order to ensure that the matching of NNLO calculations to PS preserves the shower accuracy.
The code used to produce the results presented in this work is available upon request from the authors, and will be made public in a future GENEVA release [135].
where E cm is the hadronic center-of-mass energy. For later convenience, we have introduced light-cone coordinates relative to the beam axis as n μ ≡ n μ a ¼ð1; 0; 0; 1Þ;n μ ≡ n μ b ¼ð1; 0; 0; −1Þ; ðA4Þ such that The momenta of the incoming partons are then p a ¼ x a P a ;p − a ¼ x a E cm ;p þ a ¼ 0; and we define s ¼ðp a þ p b Þ 2 ¼ x a x b E 2 cm : ðA7Þ The phase-space projection is a variable transformation from the N þ 1 phase space onto the underlying Born phase spaceΦ N , such as