Matching NNLO to parton shower using N 3 LL colour-singlet transverse momentum resummation in G ENEVA

,


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 Run3 and the future upgrade to the High-Luminosity LHC, the amount of data that will be collected will increase significantly.This groundbreaking machine will thus be able to detect extremely rare phenomena and further improve the precision of measurements.Theoretical predictions of the Standard Model (SM) 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 (PS) in fully exclusive Monte Carlo (MC) event generators.The advantage of MC event generators is the ability to produce hadronlevel events that can be directly interfaced to detector a simone.alioli@unimib.itb alessandro.broggio@unimib.itc a.gavardi@campus.unimib.itd stefan.kallweit@unimib.ite matthew.lim@unimib.itf riccardo.nagar@unimib.itg davide.napoletano@unimib.ith cwbauer@lbl.govi luca.rottoli@physik.uzh.chsimulations.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-of-theart corrections available.For this reason, in recent years there has been an increasing interest in obtaining theoretical predictions where fixed-order predictions at next-to-next-to-leading order (NNLO) accuracy are consistently matched with parton shower 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 [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 t t 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 parton shower.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 colour-singlet 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 nextto-next-to-next-to-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 colour 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 behaviour 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 impact-parameter 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 renormalisationgroup 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 ⊥ are 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 Drell-Yan (DY) production at the LHC (pp → Z/γ * → ℓ + ℓ − ) as a case study, but this approach can immediately be applied also to the charged current case and to other colour-singlet processes.Differ-ential distributions of electroweak gauge bosons play a paramount role in the precision programme at the LHC.These observables are measured at the level of their leptonic decays and are typically characterised by particularly small experimental uncertainties, which can reach the few permille 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.On the theory side, fixed-order predictions for neutral DY production have been known at NNLO accuracy for quite some time at the inclusive [42][43][44][45][46][47][48][49][50] and at the differential level [51][52][53][54][55][56].Due to the outstanding precision of the experimental data the theoretical calculations are currently being pushed at next-to-next-to-next-to-leading order (N 3 LO), and the inclusive cross section for lepton-pair production through a virtual photon has been recently calculated at this accuracy in Ref. [57].Electroweak (EW) corrections for this process are known [58][59][60][61][62][63][64][65][66].The computation of mixed QCD-EW corrections is an active area of research [67][68][69][70][71][72], and they were recently computed for the production of an on-shell Z boson and its decay to massless charged leptons has become available [73].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 hadronise 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 centre-of-mass energy of 13 TeV.The manuscript is organised as follows.In Sec. 2 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. 3 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. 4 and finally draw our conclusions in Sec. 5.

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, focussing 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-byevent 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 colour-singlet production, one has Φ 0 events: Φ 2 events: where by using the notation r N > r cut N we indicate that the events are differential in the N -jet resolution variable, whilst 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 with the jet bins commonly used in experimental analyses, which instead define jets according to a particular jet algorithm. 1efining 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 0-jet and the 1-jet regimes and the associated resummation of the 0-jet resolution variable.The expression for the MC cross sections in the exclusive 0-jet bin and in the inclusive 1-jet bin are and respectively.The label 'res' stands for 'resummed' and the nonsingular contributions labelled 'nons' must only contain terms which are nonsingular in the r 0 → 0 limit.When 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 ′ 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 ′ in SCET.In this work we extend the Geneva framework to use the transverse momentum of the colour-singlet q ⊥ as 0-jet resolution parameter, which we resum at N 3 LL accuracy, see paragraph 2 B. It is worth emphasising 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 power suppressed 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.
We will discuss the value of q cut ⊥ necessary to make the power-suppressed terms sufficiently small that they may be neglected in Sect. 3 B.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 normalisation condition dΦ 1 dΦ 0 dr 0 P(Φ 1 ) = 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. ( 3) and ( 4) are and − dσ res dΦ 0 dr 0 P(Φ 1 ) 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 parametrise 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 parametrise the 2-parton contribution B 2 in the same terms.This means that the mapping used in the projection dΦ 2 /dΦ r0 1 should preserve r 0 , i.e.
to guarantee the pointwise cancellation of the singular r 0 contributions in Eq. ( 9).The Θ r0 (Φ 1 ) in Eq. ( 10) limits the integration to the singular points for which it is possible to construct such a projection.The non-projectable 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 case of the 1 → 0 case we use the FKS mapping [75], which implies that the only non-projectable events are those which would result in an invalid flavour configuration that is not present at the LO level (e.g.qg → ℓ + ℓ − q projected to gg → ℓ + ℓ − ).By denoting the non-projectable region with the symbol Θ FKS map , we supplement the formulae 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 1 dΦ 1 (r 0 > r cut 0 ; r cut 1 ) = dσ res dΦ 0 dr 0 P(Φ 1 ) − dσ res dΦ 0 dr 0 P(Φ 1 ) dσ res dΦ 0 dr 0 P(Φ 1 ) − dσ res dΦ 0 dr 0 P(Φ 1 ) In these formulae U 1 (Φ 1 , r cut 1 ) is the Sudakov form factor that resums the dependence on r cut 1 at NLL, U ′ 1 (Φ 1 , r 1 ) denotes its derivative with respect to r cut 1 , and U (1)′ 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 behaviour of B 2 .The term V C 1 is defined from the relation The normalised splitting function P(Φ 2 ) satisfies the unitarity condition (cfr.the 0-/1-jet case Eq. ( 5)) Analogously to the case of the normalised 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 parametrised 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 Eq. ( 16) and Eq.(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 formulae above should be supplemented by including non-projectable 2-jet events with r where we use the symbol Θ r0 map (Φ 2 ) to denote the region of phase space which cannot be projected on the physical Φ 1 phase space via the Φ r0 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 1-jet resolution variable.We shall use the notation Geneva q ⊥ to label our predictions, while we will use the notation Geneva T0 when referring to the original Geneva implementation using beam-thrust.In the next section we discuss the resummation of the transverse momentum distribution in colour-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 fulfils recursive infrared collinear (rIRC) safety [102].The resummation is formulated directly in momentum space by exploiting factorisation properties of squared QCD amplitudes.The resummation is then evaluated numerically via MC methods.The RadISH formulae 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 colour 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 i-th 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 colour 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 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 hardvirtual corrections H (n) as well as the coefficient functions ci , which have been evaluated to second order for q qinitiated processes in Refs.[88,89,91].The factors P (0) denote the regularised splitting functions.The interested reader is referred to section 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[{R ′ , k i }] to denote an ensemble describing the emission of n identical independent blocks.The average of a function G(Φ 0 , {k i }) over the measure dZ 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 right-hand side of Eq. ( 24) .Eq. ( 23) has been obtained by expanding all the ingredients around k ⊥,1 since , see Ref. [24].Thanks to rIRC safety, blocks with k ⊥,i ≪ k ⊥,1 are fully cancelled by the term exp{−R ′ (k ⊥,1 ) ln(1/δ)} 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 with α S = α S (µ R ) and µ R ∼ M being the renormalisation scale, where M is the hard scale of the process.
Here L = ln(Q/k ⊥,1 ), where Q ∼ M 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 formulae are implemented in the code RadISH, which evaluates them using MC methods.We refer the reader to section 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 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 by L. 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 .

IMPLEMENTATION DETAILS
In this section we discuss the interface of the RadISH code with Geneva and the implementation of the Drell-Yan process in Geneva using q ⊥ as 0-jet resolution variable.We validate our framework by comparing our results with NNLO predictions obtained with Matrix [107] and discuss the interface with parton showers.

A. Interfacing GENEVA with RADISH
The resummation formulae discussed in Sect. 2 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 Sect.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 an 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 runtime parallelisation.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 parametrise the Born phase space using two variables2 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 discretisation 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. 3However, since the expansion is computed semi-analytically in RadISH (for additional details we refer the reader to Sect.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 fixedorder 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 artefacts 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: 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 ′ ); 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)) FIG.2: Nonsingular cumulant at order α 2 S as a function of q cut ⊥ .
and q max ⊥ is the maximum value allowed by the kinematics for each Φ 0 point.Finally, | α 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 ).However, Matrix includes an estimate of this 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 → ℓ + ℓ − production at 13 TeV in an inclusive setup by applying a cut only on the invariant mass of the produced colour 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 [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. 3 E.

C. NLO1 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. 2 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. 2 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 artefacts in a few NLO 1 differential distributions when used also in the T 1 resummation.In particular, one observes 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, 4 More precisely, the mapping preserves a recursive definition of T 0 dubbed T FR 0 , see Ref. [3] for additional details.
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 momentum pµ 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 = pz , as described in Appendix A. 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 characterised by a particular value of q ⊥ , is always retained throughout the calculation.As shown in [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 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 to fiducial cuts are included only up to order O(α 2 S ). 5he 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 phasespace mapping Φ r N (Φ N +1 ).For the following validation, the only process-defining cut is applied to the invariant mass of the resulting colour 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].We set the central renormalisation and factorisation scales to m T = m 2 ℓℓ + q 2 ⊥ , while the central resummation scale is set to Q = m ℓℓ /2.The uncertainty band is constructed as the envelope of a canonical 7-point variation of the renormalisation (µ R ) and factorisation (µ F ) scales and of two additional variations of Q by a factor of two 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 ℓℓ ⊥ 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. 3 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 regularise 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 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 ℓℓ and the rapidity y ℓℓ of the dilepton system, the rapidity of the hardest lepton y ℓ,hardest and the transverse momentum of the hardest lepton p T,ℓ,hardest .We use the same settings as specified in the previous section, but we now use a canonical 7-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,ℓ,hardest we expect possible differences because the distribution is NNLO accurate only below the Jacobian peak at m ℓℓ /2.The presence of resummation effects in the Geneva results improve the description of this observable above the peak with respect to the pure fixed-order result, smearing out the unphysical behaviour 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,ℓ,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 parton shower.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 summarise 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 at preserving as far as possible the NNLL ′ 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 ′ [3].However, due to the vectorial nature of q ⊥ , one can not 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 coloursinglet 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 parton shower 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 analysing the interface to the Pythia8 k ⊥ -ordered parton shower, 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 event 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 re-showered.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. 2 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 remains sizeable and must be handled with a 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 logarithimic order. 6In the T 1 ≪ q ⊥ limit it is however possible to approximate the correct behaviour 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 ), see Ref. [11].The new 1-jet exclusive cross section becomes Ratio to showered FIG.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 hadronisation and MPI (below).The peak (left), transition (centre) and tail (right) regions are shown.
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 permille 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 truncated-vetoed 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.245parton shower and we use the general-purpose 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 colour-singlet kinematics less compared to the standard recoil scheme [7,124].This scheme does not change the transverse momentum of the colourless 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 parton showers [125][126][127][128][129][130] and might induce spurious Ratio to partonic FIG.8: Comparison between parton level and showered results (left) and between parton level and results after shower, hadronisation and MPI (right) for the transverse momentum distribution of the dilepton pair using Geneva q ⊥ and Geneva T0 .
effects at NLL.Nonetheless, since in this work we assume that the accuracy of our predictions is dictated by the parton shower 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 parton shower.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 an 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.In the CMS Monash-Star 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 hadronisation and multi particle interactions (MPI) effects.The difference on 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 hadronised 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 hadronised 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.On the left panel, we observe that for the rapidity distribution the NNLO accuracy is mantained at a very precise level after the showering and the hadronisation stages, as should be expected by the inclusive nature of this distribution.On the right panel, the same excellent agreement is also seen in the case of the p T,ℓ,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 hadronisation and MPI (lower row).We present three separate 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 hadronisation (and MPI) stage is also added, some differences arise, localised in the small q ⊥ region.There we observe that the peak becomes more pronounced and the spectrum is more suppressed at small q ⊥ , while at large q ⊥ the showered and hadronised 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, hadronisation 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 Fig. 8 and Fig. 9 for the q ⊥ and T 0 spectra, respectively.There, the N 3 LL+NLO 1 (NNLL ′ +NLO 1 ) results at the parton level for each 0-jet resolution variable are compared to the results after showering and after hadronisation obtained using q ⊥ and T 0 in Geneva.
For the transverse momentum distribution we observe differences at the 10 − 20% level between the results obtained with Geneva q ⊥ and Geneva T0 after the shower up to values of q ⊥ close to 30 GeV, though the uncertainty bands always overlap.The Geneva T0 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 T0 spectrum is somewhat harder than that of Geneva q ⊥ .A similar pattern is visible after hadronisation 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 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 ′ +NLO 1 result.The hadronisation and especially the addition of MPI effects significantly modify the T 0 distribution, and brings 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 T0 result after hadronisation 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 T0 result displays instead differences similar to those seen for the q ⊥ distribution, both before and after hadronisation.

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. 3 D regarding the choice of PDF and central scales, and we use the shower settings as described in Sec. 3 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-flavour opposite-sign (SFOS) leptons in the final state.We then apply the following cuts to our events: p ℓ ⊥ > 27 GeV, |η ℓ | < 2.47, m ℓℓ ∈ [66,116] GeV.(37) In Fig. 11 we compare our predictions for the normalised 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 hadronisation 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 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 normalisation 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 normalised ϕ * distribution, defying a naïve 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 localised in the first few bins where hadronisation effects are more significant.

CONCLUSION
In this work we have presented an extension of the Geneva method which uses the transverse momentum of the colour singlet as the 0-jet resolution variable.As a first application, we have provided an NNLO+PS description of Drell-Yan pair production, by matching the parton level calculation with the Pythia8 parton shower.The method presented here is fully general and can be readily applied to other colour-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 0-jet resolution variable 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 permille level.We 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 parton shower and hadronisation on our predictions and found that the NNLO description of inclusive observables is preserved after both shower and hadronisation.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 hadronisation on the transverse momentum distribution, finding, as expected, that they are localised below the peak.The construction of an 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 allows us to robustly assess the impact of such choices on differential observables for the first time.We have studied the effect of the shower and hadronisation (including MPI effects) on the higherlogarithmic partonic predictions using both T 0 and q ⊥ as 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 behaviour.Nonetheless, the two predictions are in reasonable agreement, with the size of the differences never exceeding 15% for the q ⊥ and ϕ * distributions.Finally, we compared our predictions, including hadronisation 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 hadronisation 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 Drell-Yan 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 parton showers 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 parton shower 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.

FIG. 1 :
FIG.1: Damping function used to further suppress the resummation at large values of q ⊥ , see text for details.

FIG. 3 :
FIG. 3: Comparison between the resummed N 3 LL (left panel) and matched N 3 LL+NLO 1 (right panel) results obtained with Matrix+RadISH and with Geneva+RadISH.

FIG. 4 :
FIG. 4: Comparison between the NNLO results obtained with Matrix and the results obtained with Geneva+RadISH.Upper panel: invariant mass (left) and rapidity (right) of the lepton pair.Lower panel: rapidity (left) and transverse momentum (right) of the hardest lepton.

FIG. 5 :
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.

FIG. 6 :
FIG.6: Comparison of the partonic, showered and hadronised results for the rapidity of the lepton pair (left) and the transverse momentum of the hardest lepton (right).

FIG. 11 :
FIG.11:Comparison between Geneva predictions and the ATLAS data for the transverse momentum distribution (left) and for the ϕ * observable (right).

FIG. 12 :
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.

FIG. 13 :
FIG. 13: Comparison between Geneva and the CMS data for different observables.Normalised distributions are shown on the right, see text for details.
by the Swiss National Science Foundation (SNF) under contract 200020 188464.CWB is supported by the Director, Office of Science, Office of High Energy Physics of the U.S. Department of Energy under the Contract No. DE-AC02-05CH11231.We acknowledge the CINECA award under the ISCRA initiative and the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DEAC02-05CH11231, for the availability of the high performance computing resources needed for this work.LR would like to thank LBNL and UC Berkeley for the kind hospitality and for financial support during stages of this work.