Towards NNLO+PS Matching with Sector Showers

We outline a new technique for the fully-differential matching of final-state parton showers to NNLO calculations, focussing here on the simplest case of leptonic collisions with two final-state jets. The strategy is facilitated by working in the antenna formalism, making use of NNLO antenna subtraction on the fixed-order side and the sector-antenna framework on the shower side. As long as the combined real-virtual and double-real corrections do not overcompensate the real-emission term in the three-jet region, negative weights can be eliminated from the matching scheme. We describe the implementation of all necessary components in the VINCIA antenna shower in PYTHIA 8.3.


Introduction
To date it is possible to perform most colliderphysics studies with fully-differential NLO+PS matching thanks to two general, well-developed, and widely applied procedures: MC@NLO [1] and POWHEG [2]. By fully differential matching, we understand that the matching is done point by point in both the Born-and real-emission phase spaces, with a parton shower that reflects the correct singular structure of the fixed-order calculation. In this sense, fully-differential matching requires the fixed-order expansion of the shower to develop the same singularities as the fixed-order calculation up to the matched order. At NLO, this is achieved by parton showers that exponentiate terms reducing to the universal DGLAP kernels in any collinear limit and the eikonal factor in soft limits, with the colour dependence in the soft limit requiring special attention [3]. As of today, no fully-differential matching method obeying these criteria is available at NNLO, although significant progress on including higher-order corrections to parton showers has been made [4][5][6][7].
Existing NNLO+PS matching methods either extend existing merging schemes or utilise analytical resum-Email addresses: johnmc@fnal.gov (John M Campbell), shoeche@fnal.gov (Stefan Höche), haitao.li@northwestern.edu (Hai Tao Li), christian.preuss@monash.edu (Christian T Preuss), peter.skands@monash.edu (Peter Skands) mation for the transition between the fixed-order and parton-shower realms. Examples of the first kind are UN 2 LOPS [8], which extends the UNLOPS [9] scheme to the second order, and MiNNLO PS [10,11] as well as other extensions [12][13][14][15][16][17][18] of the MiNLO technique [19]. The UN 2 LOPS scheme has recently been generalised to processes with an additional jet in the context of an UNLOPS-based N 3 LO+PS matching strategy [20]. The MiNLO-based schemes may be seen as a hybrid approach, since they use a combination of analytical and numerical resummation. A noteworthy example of a scheme employing the latter approach is implemented in the GENEVA framework [21,22]. While all of these have enabled impressive phenomenological studies [8,[23][24][25][26][27][28][29][30][31][32] and provide pathways to matching precision calculations to event generators, they do not provide the same level of theoretical control as the fullydifferential matching methods that are available at NLO.
In this letter, we present for the first time a fullydifferential NNLO+PS matching scheme for final-state parton showers, restricting ourselves to the case of two coloured final-state particles. The new method combines NNLO antenna subtraction with the sectorantenna shower in VINCIA [33], suitably extended to include real-virtual and double-real corrections. A key aspect of the new technique is that the parton shower is employed only as an efficient Sudakov-weighted phasespace generator. It does not define the infrared sub-traction terms that are key to MC@NLO type matching strategies.
The letter is structured as follows. We review the matching method at NLO in section 2 before extending it to the NNLO in section 3, retaining a rather general notation. A numerical implementation in the VINCIA sector-antenna shower in PYTHIA 8.3 [34] is described in section 4, featuring a more detailed description of the matching scheme. We conclude in section 5 and provide an outlook on applications beyond the simple cases considered here.

NLO Matching Strategy
Our matching strategy generalises the technique first developed in [35], which nowadays is referred to as the POWHEG scheme [2,36,37]. To start with, it is thus useful to recap the NLO matching strategy, before moving on to the new NNLO technique.
At NLO, the expected value of an infrared-safe observable O defined on a two-particle final state process with a colourless initial state is given by where B and V denote the Born cross section and virtual correction, differential in the two-particle phase space Φ 2 . Similarly, R denotes the real-radiation cross section differential in the three-particle phase space Φ 3 , and S NLO denotes the differential NLO subtraction term in the antenna subtraction method, with its integral over the antenna phase space given by I NLO S . 1 In order to achieve a Born-local cancellation of the subtraction term upon integration over the real-emission phase space, the observable acting on S NLO must be evaluated at the reduced phase-space point Φ 2 (Φ 3 ), where the precise mapping from the three-parton to the two-parton state depends on the subtraction scheme. We can invert this mapping and factorise the phase space into the 2-particle (Born) phase space Φ 2 , and the one-particle radiation phase space Φ +1 , 1 Other well-established NLO subtraction schemes such as FKS [38] or dipole subtraction [39] may equally well be employed here.
By defining a Born-local NLO weight, eq. (1) can be rewritten as The parton-shower evolution, on the other hand, is described by a generating functional, the shower operator, recursively defined for an infrared-safe observable O by is the sum of all leading-order antenna functions 2 competing for the next branching IK → i jk off the n-parton configuration, := t t c j∈{n →n+1} with the sum and shower variables left implicit in our notation. Note that when working in the sector antenna framework [33,[40][41][42][43], Eq. (6) implicitly defines a partitioning of the real-emission term along with the associated subtractions in Eq. (3). This is crucial to avoid double-counting of radiative corrections generated by the parton shower. The associated Sudakov factor is given by Taking only the first shower emission into account, the expected value of the observable O at LO is given by This implies that upon the replacement where we have defined the 2 → 3 LO matrix-element correction factor, the following matching formula is NLO accurate up to terms appearing at order α 2 This can be seen by expanding the result to order α s .

NNLO Matching Strategy
We now turn to the main result of this work, the definition of a fully-differential NNLO matching strategy for processes with two coloured final-state particles. It is applicable to both decays of colour singlets as well as scattering processes as long as all initial-state particles are colourless, for instance as in e + e − → j j.
In the antenna formalism, the expected value for an infrared-safe observable of a process with two coloured final-state particles is given at NNLO by where RR is the differential double-real radiation cross section and RV and VV denote the differential virtual (one-loop) and double-virtual (two-loop) corrections to the real-radiation cross section R and the Born cross section B, respectively. In this context, the term S NLO denotes the differential NLO real antenna subtraction term, S denotes the differential NNLO double-real antenna subtraction term [44][45][46][47], and T denotes the differential NNLO real-virtual antenna subtraction term [44][45][46][47], Their integrated counterparts are given by I NLO S , I T , and I S . In this context, terms labelled with superscript a constitute the double-real/real-virtual subtraction terms with compensating terms labelled with a superscript c that remove spurious single-unresolved singularities. The single-unresolved singularities are captured by the NLO subtraction terms of the +1-jet calculation, labelled with superscript b, The terms labeled with superscript b and c cancel independently in eq. (11). They are constructed such as to make the integrals individually infrared finite and thus amenable to evaluation with Monte-Carlo methods.
As for the NLO case, we define a Born-local weight, which will be used to construct the NNLO matching formula, and which can be used to perform the fixed-order computation in complete analogy to eq. (4). Here, dΦ +2 is the two-particle radiation phase space that enters the factorised n + 2-particle phase space We shall further need to distinguish between an ordered and unordered component of the two-particle radiation phase space, according to the following partition of unity: The ordered part dΦ < +2 corresponds to the region accessible to strongly-ordered shower paths t 0 > t > t , whereas the unordered part dΦ > +2 is inaccessible to strongly-ordered showers because of the larger intermediate scale t 0 > t > t. We will use VINCIA's sector criterion, cf. sec. 3.3 in [33], to distinguish between the two, cf. section 4.2.
In order to be able to match the NNLO calculation with the shower, the shower needs to incorporate virtual corrections to ordinary 2 → 3 branchings as well as new 2 → 4 branchings, accounting for the simultaneous emission of two particles. These new shower terms correspond to the real-virtual and double-real corrections in the NNLO calculation. In addition, we need to incorporate the corresponding parton-shower counterterms. We start by defining the two-particle NLO Sudakov as [4] where we have introduced the 2 → 4 LO matrixelement correction factor, and the 2 → 3 NLO matrix-element correction factor w NLO 2 →3 (Φ +1 ), which we write in terms of a second order correction to the LO 2 → 3 MEC in eq. (9), The coefficientsw are given by matching the O α 2 s terms in the expansion of the truncated shower approximation to the fixed-order result in eq. (11) [4,48]. We find the fixed-order contributioñ and the second-order parton-shower matching term The factor κ is a constant and µ 2 S is the parton-shower renormalisation scale. The two are conventionally chosen such that the logarithmic structure of eq. (21) is reproduced, which leads to µ S = p ⊥ and κ 2 = exp{K/β 0 }, with K the two-loop cusp anomalous dimension [49][50][51][52]. This is known as the CMW scheme [53].
Note that in eq. (18), the integral over A (0) 2 →4 is defined over the range [t, t 0 ], since the "ordered" contribution t < t has been reabsorbed intow FO 2 →3 (Φ +1 ). It should be emphasised that we do not require the NLO three-jet calculation to be provided externally but include the correction directly in the shower evolution. This means that, different to the situation in traditional merging approaches, this correction is exponentiated into a Sudakov factor. Up to the first emission, this agrees with the treatment in [5,7] and implicitly includes the contribution from higher-order matching terms and collinear mass factorization counterterms that are needed to recover the NLO DGLAP splitting functions.
In addition, we need the 3-particle Sudakov, which we describe at LO, with the 3 → 4 LO matrix-element correction factor, Up to the second emission, the shower operator is thus given by and our final NNLO+PS matching formula takes the simple form: (26) When expanding the truncated shower operator S 2 in eq. (26) up to order α 2 s , NNLO accuracy is recovered for the observable O(Φ 2 ), while O(Φ 3 ) and O(Φ 4 ) achieve NLO and LO accuracy, respectively. This is true, because the combination of the iterated 2 → 3 → 4 and the direct 2 → 4 contributions to eq. (25) yields the correct double-real correction RR in eq. (11) by means of the LO MEC factors in eqs. (9), (19) and (24). Moreover the NLO correction eq. (20) recovers the correct real and real-virtual corrections R and RV in eq. (11) by means of eq. (9) and eq. (21).

Numerical Implementation
In this section, we want to present all necessary components of an implementation of our NNLO matching strategy. These are: • a framework to calculate the Born-local NNLO Kfactors in Eq. (15) • a shower filling the strongly-ordered [54] and unordered [4] regions of the single-and doubleemission phase space • tree-level MECs in strongly-ordered [55] and unordered [56] shower paths • NLO MECs in the first emission [48] With the exception of the first point, (processdependent) implementations of these components existed in previous VINCIA versions (not necessarily simultaneously), and have been described in detail in the various references. We have (re-)implemented all components in a semi-automated 3 fashion in the VINCIA antenna shower in PYTHIA 8.3. We access loop matrix elements via a novel MCFM [57][58][59][60] interface presented in [61] and tree-level matrix elements via a new run-time interface [62] to the COMIX matrix element generator [63] in SHERPA [64,65].
Our NNLO matching algorithm can be summarised in the following steps: 1. Generate a phase space point according to the Born cross section B(Φ 2 ). 2. Calculate the Born-local NNLO factor k NNLO (Φ 2 ) and reweight the phase space point by the result. 3. Let the phase-space maximum given by the invariant mass of the two Born partons define the starting scale for the shower, t now = t 0 (Φ 2 ). 4. Starting from the current shower scale, t now , let the 2 → 3 and 2 → 4 showers compete for the highest branching scale. 5. Update the current shower scale to be that of the winning branching, t now = max(t 2 →3 , t 2 →4 ). 6a. If the winning branching is a 2 → 3 branching, calculate the accept probability including the NLO MEC w NLO 2 →3 . • If rejected, continue from step 4.
• If accepted, continue with a LO shower from the resulting three-particle configuration, starting from t now and including the LO MEC w LO 3 →4 when calculating accept probabilities for the 3 → 4 step. When a 3 → 4 branching is accepted (or the shower cutoff scale is reached), continue with step 7. 6b. If the winning branching is a 2 → 4 branching, calculate the accept probability including the LO MEC w LO 2 →4 .
• If rejected, continue from step 4.
• If accepted, continue with step 7.
7. Continue with a standard (possibly uncorrected) shower from the resulting four-particle configuration, starting from t now .
It should be emphasised that the matrix-element correction factors make this algorithm independent of the splitting kernels (i.e. antenna functions in our case) up to the matched order and the shower merely acts as an efficient Sudakov-weighted phase-space generator. Hence, if the algorithm is stopped after step 6, an NNLOmatched result is obtained, which can be showered by any other parton shower, just as is the case for POWHEG NLO matching. Note, that there remains a dependence on the ordering variable, which has to be properly accounted for.

NNLO Kinematics
For both, the unordered shower contributions and the Born-local NNLO weight, new kinematic maps are needed to reflect their direct 2 → 4, i.e. unordered or double-unresolved, nature. We utilise that the n-particle phase space measure may be factorised into the product of a 2 → 3 antenna phase space and the n − 1-particle phase space measure, as well as into the product of a 2 → 4 antenna phase space and the n − 2-particle phase space. This allows us to write the 2 → 4 antenna phase space as the product of two 2 → 3 antenna phase spaces, dΦ +2 (p I + p K ; p i , p j 1 , p j 2 , p k ) = dΦ +1 (p I + p K ;p i ,p j , p k ) × dΦ +1 (p i +p j ; p i , p j 1 , p j 2 ) , (27) corresponding to the kinematic mapping effectively representing a tripole map [66]. In line with the phase space factorisation, the kinematic mapping is then constructed as an iteration of two on-shell 2 → 3 antenna maps given in sec. 2.3 in [33]. We have tested the validity of our kinematic maps by comparing VINCIA's phase-space mappings (doublegluon emission and gluon-emission-plus-splitting) to a flat sampling via RAMBO.

Unordered Shower Contributions
An important part of our proposal is the inclusion of double-unresolved radiation in the shower evolution. To this end, we employ the sector-antenna framework [33] and amend it by direct 2 → 4 branchings as described in [4]. In the sector-shower approach, each branching is restricted to the region in phase space where it minimises the resolution variable, defined for final-state clusterings by ) is a quark-antiquark pair (29) This is achieved by a "sectorisation" of phase space according the partition of unity, which is implemented in the shower evolution as an explicit veto for each trial branching. Since only a single branching kernel contributes per colour-ordered phase space point, sector antenna functions have to incorporate the full singularity structure associated with the respective sector. At LO, this amounts to including both the full single-collinear and single-soft limits in the antenna function. The full set of VINCIA's LO sector antenna functions is collected in [33]. By construction, the default sector shower generates only strongly-ordered sequences 4 , as the sector veto ensures that each emission is the softest (or mostcollinear) in the post-branching configuration. The inclusion of direct 2 → 4 branchings (which look unordered from an iterated 2 → 3 point of view) in the sector shower is facilitated by extending the sector decomposition in eq. (30) by an ordering criterion, where p 2 ⊥ denotes VINCIA's transverse-momentum ordering variable and hatted variables denote the intermediate node in a sequence IL →îĵˆ → i jk . Here, the scales p 2 ⊥ andp 2 ⊥ are uniquely defined by the ordering variable of the sector-shower emission, i.e., that emission which minimises eq. (29). Direct 2 → 4 emissions are thus restricted to the unordered region of the doubleemission phase space, denoted as dΦ > +2 in eq. (18) and defined as For 2 → 4 emissions off quark-antiquark and gluongluon antennae, we use the double-real antenna functions in [44,45,47]. We note that NLO quark-gluon antenna functions appear in the Standard Model at lowest order for three final-state particles and are hence not of interest for our test case of e + e − → j j. We wish to point out, however, that the NLO quark-gluon antenna functions in [46,47] contain spurious singularities which have to be removed before a shower implementation is possible.
As a validation, we show in fig. 1 the ratio of the four-jet to three-jet evolution variable for e + e − → 4 j at √ s = 240 GeV. To focus on the perturbative realm, the shower evolution is constrained to the region between t 0 = s and t c = (5 GeV) 2 . The region > 0 corresponds to the unordered part of phase space to which stronglyordered showers cannot contribute. Due to the use of sector showers, there is a sharp cut-off at the boundary between the ordered and unordered region, as the sector criterion ensures that the last emission is always the softest and therefore, no recoil effects can spoil the strong ordering of the shower. As expected, the inclusion of direct 2 → 4 branchings gives access to the unordered parts of phase space, a crucial element of our matching method.

LO Matrix-Element Corrections
In order for the shower expansion to match the fixedorder calculation, we need (iterated) 2 → 3 tree-level MECs and (direct) 2 → 4 tree-level MECs. Both take a particularly simple form in the sector-antenna framework, as will be shown below.
At leading-colour, tree-level MECs to the ordered sector shower can be constructed as [55,67] denote squared leading-colour colour-ordered amplitudes with the index i denoting the respective permutation σ i (the number of permutations depends on the process). The sector veto Θ sct j/IK ensures that only the most singular term contributes in the denominators, rendering the fraction exceptionally simple.
Direct 2 → 4 branchings can be corrected in an analogous way, replacing the sum over 2 → 3 antenna functions with a sum of 2 → 4 ones, The full-colour matrix element can be recovered on average by multiplication with a full-colour to leadingcolour-summed matrix-element weight, For gluon splittings, multiple histories contribute even in the sector shower, because all permutations of quark lines have to be taken into account. To ensure that the MEC factors remain finite for final states with multiple quark pairs, an additional quark-projection factor has to be included. Since we only deal with a maximum of two quark pairs, it is given by for 2 → 3 branchings and for 2 → 4 branchings.

NLO Matrix-Element Corrections
Making the antenna subtraction terms explicit, the fixed-order correction to the NLO matrix-element correction eq. (20) reads with the differential NLO antenna subtraction terms S NLO (Φ 2 , Φ +1 ), S NLO (Φ 2 , Φ +1 , Φ +1 ) and their integrated counterparts I NLO S (Φ 2 ), I NLO S (Φ 2 , Φ +1 ) cf. eqs. (11) and (14). Based on the argument of the last subsection, we construct the full-colour NLO matrix-element correction as The integration over the radiation phase spaces denoted Φ +1 in eq. (38) is done numerically, utilising antenna kinematics to map 3-parton configurations to 4parton configurations (similarly for 2-parton configurations). This phase-space generation approach will be described in detail in the next subsection in the context of the NNLO Born weight. Note that the radiation phase space Φ +1 in eq. (38) is generated by the shower.

NNLO Born Weight
The Born-local NNLO weight can be calculated numerically using a "forward-branching" phase-space generation approach [36,37,68,69,71], which has previously been applied to unweighted NLO event generation, using Catani-Seymour dipole subtraction [73]. The application to NNLO corrections to e + e − → 2 j using antenna subtraction has been outlined in [74].
Given a Born phase space point, the real-radiation phase space is generated by uniformly sampling the shower variables (t, ζ, φ) for each antenna, which represent integration channels in this context. As for the shower evolution, every phase space point is restricted to the sector in which the emission(s) correspond to the most-singular clusterings. The momenta of the Born+1 j point are constructed according to the same kinematic map as the shower uses, summarised in sec. 2.3 in [54]. Since antenna functions are azimuthally averaged, they do not cancel spin-correlations in collinear gluon branchings locally. To obtain a pointwise pole cancellation, the subtracted real correction R − S can be evaluated on two correlated phase space points, t, ζ, φ , t, ζ, φ + π/2 which cancels the collinear spin correlation exactly, as it is proportional to cos 2φ . To obtain double-real radiation phase space points for the subtracted double-real correction RR − S, this procedure can be iterated, yielding four angular-correlated phase space points which cancel spin correlations in double single-collinear and triple-collinear limits. Due to the bijective nature of the sector-antenna framework, each 3-or 4-particle phasespace point obtained in this way can be mapped back uniquely to its 2-particle origin, making the NNLO weight exactly Born-local. For e + e − → 2 j this procedure is identical to the one in [74].
We have implemented the NNLO antenna subtraction terms for processes with two massless final-state jets, cf. e.g. [44], in VINCIA in a semi-automated fashion. As a validation, we illustrate the convergence of the double-real radiation subtraction term eq. (12) in the triple-collinear and double-soft limits for the process e + e − → qggq in fig. 2. Phase space points are sampled according to the kinematic map in section 4.1 and we do not make use of the azimuthal averaging alluded to above.
It should be noted that a numerical calculation of the Born-local NNLO weight is not necessary for coloursinglet decays, as the inclusive K-factors are well known from analytical calculations, cf. e.g. [44,75] for Z → qq (with massless quarks), [76][77][78][79] for H → bb (with massless bs), and [45,80] for H → gg (in the Higgs effective theory).

Conclusions and Outlook
We have presented a technique to match final-state parton showers fully-differentially to next-to-next-toleading order calculations in processes with two finalstate jets. To our knowledge, this is the first method of its kind.
We have outlined a full-fledged numerical implementation in the VINCIA antenna shower in the PYTHIA 8.3 event generator. Phenomenological studies employing our strategy will be presented in separate works. We want to close by noting that, while we here focused on the simplest case of two massless final-state jets, the use of the NNLO antenna subtraction formalism facilitates its adaption to more complicated processes such as e + e − → tt or e + e − → 3 j. Considering the latter, spurious singularities in the quark-gluon NNLO antenna subtraction terms need to be removed before exponentiation in the shower. For future work, an extension of our method to processes with coloured initial states can be envisioned, given the applicability of the NNLO antenna subtraction to hadronic collisions. e + e ! qggq @ p s = 240 GeV Vincia+Comix Figure 2: Test of the convergence of the double-real subtraction term S(Φ 2 , Φ +2 , O) in eq. (12) in e + e − → qggq. Top row: progression of weight distributions from x = 10 −2 to x = 10 −4 in the triple-collinear limit (s 134 /s 1234 < x) and double-soft limit (s 134 s 234 /s 2 1234 < x). Bottom row: trajectories x · s 134 , x · s 234 , x → 0 approaching the two triple collinear limits. Phase space points are not azimuthally averaged.