Matching Fully Differential NNLO Calculations and Parton Showers

We present a general method to match fully differential next-to-next-to-leading (NNLO) calculations to parton shower programs. We discuss in detail the perturbative accuracy criteria a complete NNLO+PS matching has to satisfy. Our method is based on consistently improving a given NNLO calculation with the leading-logarithmic (LL) resummation in a chosen jet resolution variable. The resulting NNLO$+$LL calculation is cast in the form of an event generator for physical events that can be directly interfaced with a parton shower routine, and we give an explicit construction of the input"Monte Carlo cross sections"satisfying all required criteria. We also show how other proposed approaches naturally arise as special cases in our method.


Introduction
The past decade has seen substantial improvements in the accuracy of fully exclusive event generators. Matching schemes to simultaneously combine multiple leading-order (LO) matrix elements have been interfaced with parton shower (PS) routines and implemented in many event generators [1][2][3][4][5][6][7][8][9]. It has also become possible to match general next-to-leading-order (NLO) calculations with a parton shower and produce physical event samples that describe sufficiently inclusive distributions at NLO [10][11][12][13][14][15][16]. These NLO+PS event generators are now part of the standard tool set for experimental analyses and have made significant impact on phenomenology. Recently, the merging of NLO calculation of different multiplicities has been addressed by several groups [17][18][19][20][21][22][23][24][25]. Event generators continue to push to higher precision, and the LHC physics program will continue to rely on progress in this area.
The frontier of fixed-order precision is calculations at next-to-next-to-leading order (NNLO) in QCD perturbation theory. Fully differential NNLO calculations exist for several important hadron-collider processes involving W , Z, γ, and Higgs bosons as well as top quarks [26][27][28][29][30][31][32][33][34], and the technology for these calculations is continually being pushed towards more complex topologies [35][36][37]. Although experimental analyses regularly make use of NNLO cross sections and distributions, there are many challenges inherent in directly comparing fixed-order results with data.
An event generator that matches NNLO calculations with a parton shower would be an ideal tool to bridge the gap between pure fixed-order calculations and the needs of experimentalists. It would provide hadron-level events that can be more easily interfaced with an analysis while maintaining NNLO accuracy for the underlying hard process, extending the power and flexibility of an NLO+PS generator to NNLO+PS. An important first step in this direction has been taken in ref. [38], where a MiNLO-improved Powheg simulation for Higgs plus one jet [24] was used to produce an NNLO+PS event sample for Higgs boson production by reweighing the events to the NNLO Higgs rapidity distribution.
In this work we present a general method for combining NNLO calculations with leadinglogarithmic (LL) resummation to produce fully differential cross sections, and for attaching a parton shower routine to produce complete events. We derive the conditions that an NNLO+LL generator must satisfy and provide a construction that satisfies these. We also comment on the approach in ref. [38] and show how it can be derived as a special case of our results.
Theoretically, there are two conceptually very distinct aspects to interfacing a fixed-order calculation with a parton shower event generator. The first aspect is the LL improvement of the fully differential NNLO calculation. This corresponds to matching an LL resummed calculation with an NNLO calculation to obtain a combined NNLO+LL calculation, and doing so at a fully differential level. This aspect is a priori completely independent of any particular parton shower algorithm or program, and can be performed solely at the partonic (or matrix-element) level. Here, the NNLO calculation first needs to be recast in a way that is suitable for fully differential event generation. Beyond leading order, the cross section for a fixed number of partons is infrared divergent and thus ill defined, meaning that to generate physical events with a given number of partons the events must correspond to a physically well-defined and infrared-safe partonic jet cross section. In other words, each four-vector in the event should represent a partonic jet, which includes the contribution of an arbitrary number of unresolved emissions below some jet resolution cutoff. The NNLO calculation written in this way is then matched to a LL resummed calculation to obtain a combined fully differential NNLO+LL calculation.
The second aspect is to attach an exclusive parton shower Monte Carlo to this NNLO+LL calculation. In this step, events with N , N +1, and N +2 partons of the NNLO+LL calculation are handed to a parton shower algorithm, which generates additional emissions. Here, one has to take care of double-counting between the shower emissions and the partonic calculation as well as the compatibility of the LL parton shower evolution with the partonic LL resummation.
The conceptual distinction between these two aspects has already been stressed in refs. [22,39,40]. It becomes particularly important at NNLO. As we will see, the first aspect of obtaining a consistent fully differential NNLO+LL matched calculation is the more challenging one, which is why most of our discussion will focus on it. Once this step has been carried out, the step of attaching a parton shower algorithm is relatively straightforward.
This paper is organized as follows. In section 2, we discuss in detail the general framework for generating physical events beyond leading order. The main outcome of this section will be to identify the "Monte Carlo (MC) cross sections" dσ mc , which are the partonic jet cross sections according to which the different event multiplicities are distributed. In particular, we show how the fixed-order (FO) calculation is cast into this form to make it suitable for event generation. In section 3, we discuss the general procedure and conditions for combining the pure FO and pure LL calculations into a matched FO+LL calculation. As an instructive exercise we review the corresponding MC cross section for the known cases of LO+LL and NLO+LL calculations. In section 4, we then discuss in detail how to construct the MC cross sections for an NNLO+LL calculation. In section 5, we discuss how to interface the NNLO+LL calculation with a parton shower, including the conditions needed to avoid any double counting that might arise. In section 6, we discuss how our method encompasses proposed and existing approaches [22,38,41], and in section 7 we give our conclusions. Consider the cross section for some infrared-safe N -jet measurement M X , which can contain a number of cuts (θ functions) as well as differential measurements (δ functions) of observables, which we collectively refer to as X. At leading order in perturbation theory, the cross section for measuring X is given by where B N (Φ N ) is the tree-level (Born) squared matrix element for N emissions. In case of hadronic collisions we assume that the relevant parton densities (PDF) have already been convolved with the matrix elements and we will therefore avoid writing them out explicitly in our formulae. The measurement function M X (Φ N ) implements the measurement on the N -body phase space point Φ N . In particular, since M X is infrared safe it cuts off any possible IR divergences in B N (Φ N ). To obtain σ(X) from eq. (2.1) one usually performs the phase space integral over Φ N numerically. Due to the large dimensionality of N -body phase space, the typical method of choice is Monte Carlo integration: We generate points Φ N with relative weights such that they are distributed according to B N (Φ N ). 1 For each generated point Φ N , we evaluate M X (Φ N ) and record the result for X into appropriate histograms with the associated weight of the point Φ N . At next-to-leading order in perturbation theory, σ(X) is given by 2) The virtual one-loop contribution V N and the (N +1)-parton real-emission contribution B N +1 are separately IR divergent. A convenient way to handle these divergences is the standard subtraction method, where one writes 2 Here, V C N denotes the virtual contribution including the appropriate integrated subtraction terms to render it IR finite. The C m N +1 are the corresponding real-emission subtraction terms. Written in this way, the Φ N and Φ N +1 integrals are separately IR finite and can each be performed numerically by Monte Carlo integration.
The Φ N integral in eq. (2.3) can be performed as before at LO, except that the Φ N points are now distributed according to B N + V C N . The Φ N +1 integral is more involved now due to the presence of the subtraction terms. Their precise form is not important for our discussion. What is relevant is that generically several subtraction terms are needed to remove all possible IR singularities in B N +1 , and that in each subtraction term the measurement must be performed on a (in principle) different projected N -body phase space pointΦ m N (Φ N +1 ). 1 To be precise, if ΦN points are generated according to a probability distribution P (ΦN ), each point gets As a result, each generated point Φ N +1 contributes multiple times to each histogram with multiple weights distributed according to B N +1 and C m N +1 , which are separately IR divergent. As we approach any IR-singular region, the different X values obtained for the real emission term and the relevant subtraction terms approach each other and eventually fall into the same histogram bin, where the IR-divergent contributions of real emission and subtractions cancel each other.

Monte Carlo event generation
The above Monte Carlo phase space integration is how essentially all (N)NLO programs using subtractions operate. Its main feature is that it allows one to obtain the exact result (up to limitations due to numerical precision) for arbitrary IR-safe observables. It can be contrasted with the event generation used in (parton shower) Monte Carlo event generators. In an event generator, the basic goal is to produce physical events that are generated and stored once and that can be repeatedly processed later, e.g., by performing various measurements on them.
Theoretically, performing a measurement M X on the stored events is exactly equivalent to making a theoretical prediction for σ(X). To illustrate this with a trivial example, imagine we want to compute σ LO (X) in eq. (2.1) by generating events. To do so, we take dσ mc We now first generate a number of points Φ N (the actual generation routine can be the same as before), call them "N -parton events", and store them together with their weights. These events are distributed according to the "MC cross section" dσ mc ≥N /dΦ N . In the second step, we run over all stored events, evaluate the measurement M X (Φ N ), and record the result for X into histograms with the associated weight of each event. The result for σ LO (X) obtained in this way is obviously identical to that obtained by performing the Monte Carlo integration of eq. (2.1) as described there. We have merely changed from two operations in a single loop into two separate loops with one operation each. In practice, this separation becomes vital as soon as the additional processing steps performed on the events become very involved (theoretically and/or computing intensive). This is the case when the events are run through a parton shower and hadronization routine, which then also allows one to perform much more detailed measurements, such as propagating them through a complete detector simulation and using them in different experimental analyses. Now, if we try to perform the NLO calculation in eq. (2.3) with the same approach, then for each generated and stored Φ N +1 point with weight proportional to B N +1 we would also have to keep track and store the complete set of associated (correlated) Φ m N events with weights −C m N +1 (Φ N +1 ). In principle, this is possible and would again give the identical result for σ(X) as before (some fixed-order programs can indeed be run in this mode). However, for experimental purposes, e.g. when matching onto parton shower routines, it is impractical to deal with such "effective" events that consist of a number of correlated unphysical events with large and opposite weights. The point is that B N +1 and C m N +1 separately are not physical cross sections. Their individual contributions are IR divergent and the divergences only cancel each other to give a physical result once they are combined into a physical measurement, i.e., a single histogram bin. Therefore, the goal is to generate events that are physical in the sense that the contribution from each event should correspond to an IR-safe cross section, i.e., all IR divergences should cancel on a per-event basis rather than between several unphysical events. 3 Conceptually, this implies that each N -parton event should be considered a "bin entry" in a partonic N -jet measurement which is IR finite and fully differential in the corresponding partonic Njet phase space. In other words, the generated N -parton events really represent points in an N -jet phase space rather than N -parton phase space.
The definition of an N -jet cross section requires the presence of an N -jet resolution variable, which we call T N . It is defined such that in the IR singular region T N → 0. Emissions below T N < T cut N are considered unresolved and T cut N is called the N -jet resolution scale. When generating events with N and N + 1 partons, they are distributed according to the following Monte Carlo (MC) cross sections: The cross section σ(X) measured from these events is given by Physically, dσ mc N /dΦ N (T cut N ) is a fully differential exclusive partonic N -jet cross section. Perturbatively, it is the cross section for the emission of N identified partons plus any number of unresolved emissions below the resolution scale T cut N . (At higher orders this includes the necessary virtual corrections to render it IR finite). Hence, as mentioned already, Φ N really means Φ jet N here, and when specifying the jet resolution variable T N , one also needs to specify how unresolved emissions with T N < T cut N are projected onto the partonic N -jet phase space Φ jet N in which the events are distributed. To avoid cluttering the notation, we suppress the explicit "jet" label in the rest of the paper.
The cross section dσ mc ≥N +1 /dΦ N +1 (T N > T cut N ) in eqs. (2.5) and (2.6) is an inclusive partonic (N +1)-jet cross section. Perturbatively, it is the cross section for the emission of N +1 identified partons above the N -jet resolution scale T cut N . It includes any number of additional 3 Note that the problem is not the use of weighted events to obtain the desired distribution, since as long as the weighted events are statistically independent they can be (partially) unweighted. What is very impractical is to have unphysical events that must be treated as correlated due to their individual weights being IR divergent, since there is no reasonable way to unweight these. One can also have an "intermediate" case, where the final cross section is made up of independent IR-finite parts, some of which still require events with negative weights. This causes much less severe but still important practical complications and so should be avoided if possible.
emissions, which are mapped onto the partonic (N + 1)-jet phase space Φ N +1 ≡ Φ jet N +1 of the N + 1 identified partons (or rather partonic jets). The jet resolution variable T N is part of the full Φ N +1 and we use the argument T N > T cut N to explicitly indicate the fact that dσ mc ≥N +1 only has support for T N above T cut N . This procedure is essentially what every generator of physical events does, either implicitly or explicitly. For example, in a pure parton shower generator, T N corresponds to the shower evolution variable and T cut N is the parton shower cutoff. In this case, dσ mc N /dΦ N (T cut N ) is the no-emission probability, and dσ mc ≥N +1 /dΦ N +1 (T N > T cut N ) is the probability to have at least one emission above T cut N . This is discussed in detail in section 2.3. We now want to cast the FO calculation in eq. (2.2) into a form suitable for event generation by applying the logic in eqs. (2.5) and (2.6) at fixed order. We start by considering the trivial example of an LO calculation. Since at tree level there are no additional emissions, we do not need to specify a resolution variable, the N jets coincide with the N tree-level partons, and measuring the N -jet phase space simply returns the full N -parton information. Thus, at LO the "MC measurement" function defining the MC cross sections is i.e., the partonic phase space Φ N going into the measurement is mapped trivially onto the partonic N -jet phase space Φ N ≡ Φ jet N of the Monte Carlo events. Inserting this into the LO calculation in eq. (2.1), we obtain dσ mc which is the obvious result and corresponds to eq. (2.4). Starting at NLO, the fully differential MC measurement becomes nontrivial. We now need to specify how the measurement function acts on both Φ N and Φ N +1 points. At NLO, the definition of the MC cross sections given below eq. (2.6) corresponds to the fully differential MC measurements For these to be IR safe, T N (Φ N +1 ) can be any IR-safe resolution variable, andΦ N (Φ N +1 ) can be any IR-safe projection from Φ N +1 to Φ N . In particular, T N (Φ N ) = 0, and T N (Φ N +1 ) > T cut N cuts off all IR-singular regions in Φ N +1 . Below the resolution scale T cut N , the additional emission in Φ N +1 remains unresolved and Φ N +1 is projected onto a corresponding Φ N point viaΦ N (Φ N +1 ). Above T cut N , the additional emission is resolved and we measure the full Φ N +1 dependence. Inserting eq. (2.9) into eq. (2.2), we obtain where in the first equation we have abbreviated Using eq. (2.10) as the MC cross sections in eq. (2.5) we can generate physical NLO events.
Of course, to distribute our N -parton events we still have to perform the NLO calculation in dσ mc N /dΦ N (T cut N ) (which may be nontrivial and require subtractions, but which we will assume exists).
We can ask to what extent other measurements M X are reproduced at NLO when using eq. (2.10) together with eq. (2.6), Comparing to eq. (2.2), it is clear that observables are correct to the appropriate fixed order if and only if they are insensitive to the unresolved region of phase space below T cut N where the measurement is evaluated on the projected phase space pointΦ N (Φ N +1 ) rather than the exact Φ N +1 . That is, • N -jet (integrated) observables are correct to NLO N up to power corrections that scale as O(α s T cut N /T eff N ), where T eff N is the typical resolution scale to which the measurement is sensitive to, i.e. up to which it integrates over Φ N +1 . In particular, it should contain the complete unresolved region of Φ N +1 where T N (Φ N +1 ) < T cut N .
• (N + 1)-jet (differential) observables are correct to LO N +1 if they only include contributions in the resolved region of Φ N +1 , i.e., if their M X (Φ N +1 ) completely excludes the unresolved Here, M -jet observables are those that receive their first nonzero contribution from an Mparton final state, and N n LO M refers to the O(α n s ) correction relative to the corresponding tree-level M -parton result.
An example of the effective resolution scale T eff N is in Higgs boson production with a veto on extra jets (requiring p jet T < p cut T ). If the resolution variable T N is chosen to be the transverse momentum of the hardest jet, then T eff N = p cut T . For a different resolution variable, T eff N corresponds to the effective scale in T N to which the cut on p jet T is sensitive to. For example, if T N is chosen to be the p T of the Higgs, then T eff N p cut T . If it is chosen to be beam thrust [42], then T eff The presence of power corrections in T cut N /T eff N clearly highlights the formal limitation fundamental to the event generation method, namely that we inevitably lose the fully differential information below the resolution cutoff. This is the price we have to pay for the event-by-event IR-finiteness. Fortunately, in practice, this is not a problem, since we can always make T cut N small enough such that either power corrections in T cut N are irrelevant or else, if we do probe scales of order T cut N , the FO expansion breaks down and resummed perturbation theory is required to obtain a stable prediction. In this case, the only observables for which we cannot obtain an accurate FO result are those for which we would not want to use the FO calculation in the first place.
One might think that the breakdown of the FO expansion indicates that our events also become unphysical again. However, the important point is that the events (or more precisely the underlying MC cross sections) are still defined in a physical IR-safe way. For very small T cut N we are simply going into an extremely exclusive and thus IR-sensitive region where the FO calculation itself breaks down, irrespectively of how it is performed. This is precisely the region where improving the FO calculation with the parton-shower LL resummation or a higher-order resummation becomes necessary to obtain a meaningful perturbative result. Rewriting the FO calculation in this way forms the basis (and in fact is a necessary precondition) for combining it with a parton shower event generator. As we will see later, after including the LL improvement T cut N will become equivalent to the parton shower cutoff.

Event generation at NNLO
To implement an NNLO calculation in the form of event generation, we first have to extend eq. (2.5) to include (N + 2)-parton events. To do so, we split dσ mc ≥N +1 into an exclusive dσ mc N +1 and an inclusive dσ mc ≥N +2 using an additional (N + 1)-jet resolution scale T cut N +1 . Events with N , N +1, and N +2 partons are then distributed according to the following MC cross sections: The cross section σ(X) measured from these events is given by (2.14) Here, dσ mc N (T cut N ) is defined as before as an exclusive partonic N -jet cross section, i.e., the IR-finite cross section for N identified partons plus any number of unresolved emissions below Figure 1. Illustration of the N -jet, (N + 1)-jet, and (N + 2)-jet regions in eq. (2.13) for resolution variables that satisfy T N +1 < T N (e.g., the p T of the leading and subleading jet or N -jettiness [44]). The N -jet bin has T N < T cut N and is represented by N -parton events with T N = T N +1 = 0 (shown by the black dot at the origin). The (N + 1)-jet bin has T N > T cut N and T N +1 < T cut N +1 and is represented by (N + 1)-parton events with T N +1 = 0 (shown by the black line on the T N axis). The inclusive (N + 2)-jet bin has T N > T cut N and T N +1 > T cut N +1 and is represented by (N + 2)-parton events.
the resolution scale T cut N . Next, dσ mc N +1 (T N > T cut N ; T cut N +1 ) is an exclusive partonic (N + 1)-jet cross section and also IR finite. It contains N + 1 identified partons plus any number of unresolved emissions below the resolution scale T cut N +1 . The argument T N > T cut N indicates that the cross section only has support above T cut N , which acts as the condition to have one additional resolved parton. Finally, dσ mc ≥N +2 (T N > T cut N , T N +1 > T cut N +1 ) is an inclusive partonic (N + 2)-jet cross section and also IR finite. It contains at least N + 2 identified partons, where two additional partons are required to be above T cut N and T cut N +1 , respectively, as well as any number of additional emissions. Compared to eq. (2.5), where N + 1 was the highest multiplicity and inclusive over additional emissions, now both N and N + 1 are exclusive multiplicities, while the highest multiplicity is N + 2 and again inclusive over additional emissions. In figure 1, we illustrate the regions in T N and T N +1 contributing to each multiplicity.
At fixed NNLO, the cross section σ(X) is given by 15) where W N contains the two-loop virtual corrections for N partons and V N +1 the one-loop virtual corrections for N + 1 partons. In principle, the phase space integrals in eq. (2.15) can again be performed by Monte Carlo integration using subtractions. Since the singularity structure of the real, virtual, and real-virtual contributions is much more complex than at NLO, the required subtractions are far more intricate now. We now want to recast eq. (2.15) in the form of eq. (2.14). At NNLO, the general definition of the MC cross sections given below eq. (2.14) corresponds to the following MC measurement functions: For these measurements to be IR safe, T N and T N +1 can be any IR-safe resolution variables and the variousΦ N (Φ M ) can be any IR-safe phase space projections. These conditions are much more nontrivial at NNLO compared to NLO, since we now need explicit projections from Φ N +2 down to Φ N , and furthermore the condition T N (Φ N +2 ) > T cut N must cut off all double-unresolved IR-singular regions of Φ N +2 . For example, at NLO T N could simply be defined as the p T or virtuality of the one additional emission (which is IR safe at NLO). However, taking T N and T N +1 as the p T or virtuality of each of the two additional emissions is not IR safe at NNLO. Instead, a properly IR-safe NNLO generalization for T N would be to define it as the p T of the additional jet using an explicit jet algorithm with some jet radius R. This corresponds to using a "local" resolution variable. Another choice is to define it as the p T of all additional emissions or N -jettiness [44]. These correspond to "global" resolution variables.
Plugging eq. (2.16) back into eq. (2.15), we obtain the required MC cross sections, where we have defined the generalization of eq. (2.11), Note that the implementation of the constraint T N > T cut N in dσ mc N +1 is nontrivial now. For simplicity, we have not written any subtractions in eq. (2.17), which will be needed in some form when evaluating the cross sections numerically to separate out and cancel the IR divergences in the virtual and real emission contributions. Applying the MC measurement functions in eq. (2.16) to the required subtraction terms is straightforward. The precise form of the subtractions is however not important for our discussion, and one can apply for example the NNLO subtraction techniques in refs. [45][46][47][48].
As at NLO, writing the NNLO calculation in terms of IR-finite MC cross sections as above forms the basis for using it in an exclusive event generator for physical events. Using eq. (2.17) together with eq. (2.14) the cross section for some measurement M X obtained in this way is This has the same inevitable limitations that we already saw in the NLO case. Since N -parton and (N + 1)-parton events correspond to partonic N -jet and (N + 1)-jet cross sections, the measurement is evaluated on the corresponding projected phase space points in the unresolved regions of phase space. Therefore, the cross section σ(X) is correct to the required fixed order (up to power corrections in the resolution scales) for measurements X that are insensitive to the unresolved regions of phase space. This means: • N -jet observables are correct to NNLO N if they integrate over the complete unresolved regions of Φ N +1 and Φ N +2 . [Power corrections are at most of relative O(α s T cut and T eff N are the typical resolution scales up to which the measurement integrates over Φ N +1 and Φ N +2 , and generically T eff • (N + 2)-jet observables are correct to LO N +2 if they only include contributions in the resolved region of Φ N +2 .
As before, M -jet observables receive their tree-level contribution from an M -parton final state, and N n LO M refers to the O(α n s ) correction relative to that. The definition of T eff N can be understood using an example similar to that used when discussing MC cross sections at NLO. These properties are fundamental to the event generation method and are shared by all implementations. In turn, they will also be the necessary conditions on the FO accuracy that should be maintained by the NNLO+LL calculation.
Although T cut N and T cut N +1 are jet resolution scales, they will typically not define jets that are reasonable to measure experimentally. They effectively serve as IR cutoffs below which observables should be inclusive over unresolved emissions (which in fact means they should be smaller than the typical scales probed in the experimental jet measurements). In practice, T cut N and T cut N +1 can again be made sufficiently small such that FO perturbation theory is no longer appropriate to describe observables that probe emissions at or below these scales. As at NLO, at this point we are not losing any relevant fixed-order information and the parton shower or higher-order resummation is required to provide a valid perturbative description.
To conclude this subsection, we stress that so far we have not done any showering, we have simply rewritten the FO calculation in a form suitable to generate physical events. This will be our starting point for obtaining a fully differential NNLO N +LL calculation and defines the partonic jet cross sections that we require as inputs from the FO calculation. We assume these are available to us and we will not discuss the techniques used to compute them. For dσ mc N +1 and dσ mc ≥N +2 these are the same inputs that are required in the corresponding NLO N +1 +LL calculation. The genuine NNLO input required is the cumulant cross section dσ mc N /dΦ N (T cut N ). We assume that it is provided to us by the FO calculation in a form that allows us to obtain a numerical result for any needed Φ N point and T cut N value. This is likely to be a challenging part in the practical implementation, and its availability might restrict the possible choices for the concrete definitions of T N (Φ N +2 ) andΦ N (Φ N +2 ) that can be used.

Event generation at LL
The parton shower produces events whose cross sections include resummed contributions from all orders in perturbation theory. These resummed rates account for the large cancellations between virtual and real emissions in the IR region of phase space. The shower can therefore describe the resummation region of observables more accurately than FO calculations, as well as produce high-multiplicity final states than can be passed through hadronization routines to produce realistic events. In this subsection, we are interested in using the parton shower approximation to obtain a resummed calculation for the Monte-Carlo cross sections at leadinglogarithmic (LL) order. This will serve as the basis for the LL improvement of the FO cross sections to obtain matched FO+LL calculations in sections 3 and 4. Note that here we are not interested in the algorithmic construction of the parton shower. Formulating the LL calculation in a parton-shower-like fashion will facilitate attaching an actual parton shower to the matched FO+LL calculation.
The parton shower directly works as an event generator and is fundamentally based on evolution in a resolution variable T , which characterizes the scale of an emission. Subsequent emissions occur at increasingly smaller values of T , down to a low-scale cutoff T cut ∼ 1 GeV, where the perturbative parton shower description ceases to be valid. Below this cutoff one enters the nonperturbative regime, where hadronization models are used. In the leadinglogarithmic limit, all emissions are strongly ordered, i.e., each emission occurs at a much smaller value of T than the previous one, such that all emissions can be considered independent. Due to this single-emission nature, at LL there is no distinction between global and local resolution variables that are equivalent for a single emission. Hence, we can define the N -jet resolution variable T N as the emission scale T of the N + 1st emission, with the resolution scale T cut N given by the shower cutoff T cut , i.e., To start, we consider an N -jet process (with N partons at the Born level) and are interested in generating events with N and N + 1 partons as in eqs. (2.5) and (2.6). The MC cross sections using the above N -jet resolution variable are then given at LL order as where all ingredients and the notation we have introduced are discussed in detail in the following. To shorten the notation, we will often drop the explicit dependence on Φ N +1 for most objects, as in the last line of eq. (2.21), but one should keep in mind that in general all objects which depend on the emission label m (which is explained below) have Φ N +1 as their argument.
is the N -parton Sudakov factor, which effectively sums the dominant contribution from an arbitrary number of unresolved emission below T cut N at LL, corresponding to the general definition of dσ mc the discussion below eq. (2.6)]. It can be written as 22) where P N (Φ N , T ) is a global N → N +1 splitting function which sums over all possible singleparton emissions from each parton in Φ N at the emission scale T . It arises from projecting the full emission phase space dΦ N +1 /dΦ N , which contains the complete set of splitting variables, onto the resolution variable T : The m labels in eqs. (2.21) and (2.23) run over all the possible (IR-singular) emission channels (q → qg, g → gg, g → qq, etc.), including the information of which parton in Φ N was split and which two partons in Φ N +1 resulted from the splitting. For each emission channel m, T m (Φ N +1 ) determines the relevant emission scale and the splitting function P m N (Φ N +1 ) contains all coupling and kinematic prefactors times the usual Altarelli-Parisi splitting function. For simplicity we keep the upper limit T < T m max on the emission scale T implicit in the definition of P m N . 4 Finally, the projectionΦ m N (Φ N +1 ) can be any IR-safe projection and as before specifies how the partonic Φ N +1 is mapped onto the partonic N -jet phase space point Φ N ≡ Φ jet N in which the N -parton events are distributed. The projection can be different for each m. (As far as the parton shower goes,Φ m N is the inverse of the momentum reshuffling performed when Coming to dσ mc ≥N +1 in eq. (2.21), the differential parton shower rate for the emission with index m is given by its splitting function times the Born contribution, For future use we also define which is the LL approximation of the full real emission contribution B N +1 in the IR-singular limit. The Sudakov factor ∆ N (Φ m N ; T m N ) appearing in dσ mc ≥N +1 in eq. (2.21) is the same as in eq. (2.22) but evaluated at the emission scale T m N . It effectively resums the contributions from arbitrary additional emissions below T m N at LL. The cross section for some measurement M X obtained from the LL MC cross sections in eq. (2.21) is To discuss its perturbative accuracy we define is a typical hard scale in the process. Formally, the resummation corresponds to a reorganization of the perturbative series, which is achieved by expanding in α s while The leading-logarithmic order is O(1) in this counting. For the cumulant cross section integrated up to T cut N , this corresponds to resumming all terms ∼ α n s L 2n cut relative to the Born cross section, while for the cross section differential in T N , this corresponds to resumming all terms ∼ α n s L 2n−1 /T N . For a general measurement this means: • N -jet (integrated) observables are correct to LL resumming all terms ∼ α n s ln 2n (T eff N /Q) where here T eff N is the typical resolution up to which the measurement is integrated. (In particular, for dσ mc The parton shower intrinsically preserves probability, which is a consequence of the fact that it is formulated as a Markov chain process with the probability of each emission given by the exact differential of the integrated probability. Taking the special case where , we precisely reproduce the total leading-order N -jet cross section from eq. (2.26), Here, we used the fact that the differential T N spectrum is the exact derivative of the integrated T cut N cumulant cross section, As a result, the T cut N dependence precisely cancels between the cumulant and the integrated spectrum in eq. (2.29). For a general measurement M X (Φ N +1 ) that cannot be written in terms of the shower projectionΦ m N , the LO cross section is reproduced up to small power corrections ∼ T cut N /Q, which introduce a small residual T cut N dependence.
In the resummation counting of eq. (2.28) the Sudakov factors in eqs. (2.26) and (2.29) are O(1), and in particular 1 − ∆ N (T cut N ) ∼ O(1), despite the fact that its FO expansion would start at α s , which is essential for eq. (2.29) to work out. What happens is that S N +1 ∼ α s L/T N , which upon integration over T N > T cut N becomes α s L 2 cut ∼ 1. In other words, the T N spectrum at small T N is O(1) at LL, even though in fixed order it only starts at α s .

Combining fully differential FO calculations with LL resummation
In this section, we discuss the general conditions to combine the fully differential FO and LL calculations in an event generator. After the general discussion in section 3.1, we will review the LO+LL and NLO+LL cases in the following subsections. The NNLO+LL case is then discussed in detail in section 4.

General discussion
The goal of combining the FO calculation with the LL resummation is to improve the perturbative accuracy in the resummation region, where the FO expansion itself becomes invalid, to attain at least the O(1) accuracy provided by the LL resummation there. At the same time, the perturbative accuracy of the FO calculation must be maintained in the FO region where the resummation is unimportant.
As a necessary precondition, the combined FO+LL calculation must be simultaneously correct to the desired fixed order (LO, NLO, etc.) and resummation order (LL, NLL, etc.). Here, the fixed order is counted as usual by powers of α s , while the resummation order is dictated by the logarithmic counting in eq. (2.28), where L = ln(T N /Q) and L cut = ln(T cut N /Q) [see eq. (2.27)]. Therefore, the MC cross sections of the FO+LL calculation have to satisfy the conditions which require that upon expanding/truncating the MC cross sections to either FO or LL, denoted by [· · · ] FO or [· · · ] LL , the pure FO or LL results appearing on the right-hand sides in eq. (3.1) correctly reproduce the results in section 2. These conditions ensure that the input MC cross sections for each event multiplicity have the desired perturbative accuracy in both the resummation and fixed-order regions. For example, at NLO+LL, where we need events with N and N + 1 partons, the MC cross sections dσ mc N and dσ mc ≥N +1 are correct to NLO N +LL and LO N +1 +LL, respectively. Similarly, for NNLO+LL, where we need events with N , N + 1, and N + 2 partons, the corresponding dσ mc N , dσ mc N +1 , and dσ mc ≥N +2 are correct to NNLO N +LL, NLO N +1 +LL, and LO N +2 +LL, respectively.
We also have to achieve the desired perturbative accuracy at FO and LL for general measurements M X . As discussed in section 2, when generating physical events, σ(X) is predicted at the desired accuracy only up to power corrections in the resolution scale T cut N , which should therefore be as small as possible. At the same time, for integrated N -jet observables the residual dependence on the resolution scale T cut N in the pure FO and LL calculations is at most power suppressed. The important condition is now that the same must also hold for the combined FO+LL calculation. Therefore: • Since T cut N must be taken as small as possible to minimize power corrections, it is imperative that logarithms of T cut N must be counted as in eq. (2.28), for which we adopt the notation O cut , such that α n s L m cut ∼ O cut (α n−m/2 s ).
• For integrated N -jet and (N + 1)-jet observables that in fixed order are predicted at α n s with corrections starting at O(α n+1 s ), any residual logarithmic dependence on the jet resolution scales T cut N and T cut N +1 must be O cut (α ≥n+1 s ), i.e., only give corrections at the level of accuracy (or higher) as expected from higher FO corrections.
To ensure this, the conditions in eq. (3.1) alone are not sufficient. In addition, the MC cross sections for different multiplicities must be consistent with each other and satisfy the relation 6 violations for an N n LO N +LL calculation. (The missing exact dependence on Φ N +1 below T cut N will still introduce the same power corrections in T cut N for general measurements M X as in the pure FO and LL cases.) This condition enforces that after projecting the fully differential Φ N +1 dependence onto {Φ N , T N } the differential T N spectrum is the derivative of the cumulant with respect to T cut N (for any fixed Φ N ). Integrating eq. (3.2) over T N we obtain the equivalent condition for the cumulant being the integral of the T N spectrum. That is, for any T c N (and fixed Φ N ) violations for an N n LO N +LL calculation. In figure 2, we show how the FO and resummed contributions determine the accuracy of the cross sections in different regions of phase space. In table 1, we summarize the perturbative accuracy as well as the size of uncontrolled higher-order corrections from fixed order, resummed, and residual resolution scale dependence for integrated N -jet observables and differential (N + 1)-jet observables for various FO+LL orders.  of the form α n s L 2n−6 (1). Hence, such a T cut N dependence would spoil the desired O(α 2 s ) accuracy of the NNLO+LL calculation. When increasing the FO accuracy, the condition in eq. (3.2) becomes more and more stringent and thus more challenging. As we saw in section 2.3, in the LL calculation the cancellation of the T cut N dependence to all orders is achieved by virtue of the fact that the differential cross section in T N is given by the exact derivative of the cumulant cross section with respect to T cut N . The same is also obviously true for the pure FO calculation. Therefore, a simple and generic method to ensure the cancellation of the resolution scale dependence (up to power corrections) also for the FO+LL calculation is to explicitly construct the spectrum and cumulant by enforcing eqs. (3.2) and (3.3) exactly. There are different choices for doing so, as we will see in section 4, as well as different options for the practical implementation, which we will come back to in section 6.
Note that a priori we do not require the resummation order to match the perturbative accuracy of the fixed order. For example, the NLL terms in an NNLO+LL cross section are allowed to be incorrectly predicted even though in the resummation region they are formally more important than the NNLO terms. These higher-order resummed terms will affect observables in the singular regime at small T eff N but not observables at large T eff N , which are controlled by FO corrections. In section 4 we will explicitly see how the mismatch between the LL resummation and the NNLO calculation enters. A consistent matching of fixed order Table 1. Perturbative accuracy of N -jet (integrated) and (N +1)-jet (differential) observables satisfied at different FO and FO+LL. Here T eff N is the effective scale to which the observables are sensitive. For T eff N ∼ Q, the perturbative accuracy is set by the FO expansion, with corrections from higher FO contributions as well as residual T cut N dependence. (The latter will depend on the details of the matching so we show the minimal required accuracy which has to match the FO level of accuracy, see the discussion of eq. (3.2) for more details.) For T eff N Q, the perturbative accuracy is set by the resummation counting in eq. (2.28). and resummation at the same perturbative accuracy would clearly be a desirable feature. As was shown in ref. [22], by performing the resummation at NNLL, the merging of two NLO calculations with different multiplicities arises as a byproduct. Maintaining the perturbative accuracy with higher-order matrix elements and higher-order resummation is obviously more challenging as more ingredients are required and additional complications arise, e.g., one has to employ a resolution variable that is resummable to the desired order. These issues were thoroughly addressed in ref. [22], and we discuss the connection in section 6.1.

LO+LL
The LL calculation performs the LL resummation in T N and T cut N , as outlined in section 2.3. It naturally contains the full LO N contribution, so it is already LO N +LL correct, but does not include the full contribution from the LO ≥N +1 matrix elements for additional jet multiplicities (beyond the shower approximation). The goal of LO+LL matching is to combine the LO ≥N +1 calculations with the LL resummation, an example of which is the CKKW method [1][2][3]7].
Considering the matching of LO N , LO N +1 , and LL, denoted as LO N,N +1 +LL, the exclusive N -jet and inclusive (N + 1)-jet MC cross sections are Here, the B m N +1 are defined such that B N +1 = m B m N +1 , and whenever an emission m becomes IR singular B m N +1 contains all its divergences. A possible choice would be to take ). For ease of notation, from here on we always group the emission label m on expressions with the notation m {· · · } m to denote that all relevant terms within the curly brackets receive a label m.
The cross sections in eq. (3.4) are correct to LO N and LO N +1 respectively simply because any corrections to B N or B N +1 are of higher fixed order. The Sudakov factors multiplying the Born contributions render the N -jet cumulant correct to LL in T cut N and the (N + 1)-jet spectrum correct to LL in T N .
To discuss the perturbative accuracy of integrated N -jet observables from residual T cut N dependence, we rewrite dσ mc ≥N +1 in eq. (3.4) as The first term on the right-hand side is identical to the pure LL cross section, and when projected onto Φ N and integrated over , which exactly cancels the T cut N dependence in the cumulant dσ mc N (T cut N ) [see eq. (2.30)]. The second term corresponds to the FO matching correction making dσ mc ≥N +1 to be LO N +1 accurate. Its T cut N dependence is determined by the accuracy of B N +1 − S N +1 . If this difference contains subleading singular dependence on T N , which would be terms ∼ α s /T N , then the T cut N dependence in integrated N -jet observables will be of order α n s L 2n−1 cut ∼ O cut (α 1/2 s ). Interestingly, this is not actually sufficient to preserve the 1+O(α s ) accuracy required at LO N (see table 1). In the case that S N +1 does reproduce the full singular structure of B N +1 (which generically will not be the case for parton showers), then the residual T cut N dependence will only appear as O cut (α s T cut N ) power corrections. Improved LO+LL methods that explicitly remove this residual O cut (α 1/2 s ) dependence and restore the LO N accuracy have been discussed in detail in refs. [21,41,49,50]. They essentially enforce the consistency conditions in eq. (3.3).
Finally, we note that at LO N,N +1 +LL another possible valid choice for dσ mc ≥N +1 is to take where compared to eq. (3.5) we have dropped the Sudakov factor in the last line. The T cut N dependence in this case is different numerically but of the same accuracy as for eq. (3.5), depending in the same way on the extent to which S N +1 reproduces the IR singularities of B N +1 .

NLO+LL
The matching of fully differential NLO calculations to parton shower routines has been addressed by several frameworks [10,12,15,22,51,52]. Here we review the general structure of the underlying matched NLO+LL calculation. The MC cross sections underlying the MC@NLO [10] and Powheg [12,13] approaches are given by 7 where  is essentially the inclusive NLO N cross section, but using the real emission given by S N +1 instead of B N +1 . This means that S N +1 must contain the full IR singularities of B N +1 in the limit T N → 0, such that upon integration the virtual IR divergences of V N are canceled in eq. (3.9). We can easily check that eq. (3.7) is correct to NLO and LL, i.e., that it satisfies eq. (3.1). Dropping the NLO corrections, which amounts to taking dσ S ≥N → B N and dropping the dσ B−S N in dσ mc N , we reproduce the LO N,N +1 +LL result in eq. (3.6). Using the fixed O(α s ) expansion of the Sudakov, we see that expanding eq. (3.7) to NLO exactly reproduces eq. (2.10) at NLO N and LO N +1 , where the T N in the NLO calculation is now the same m-dependent resolution variable that is used in the LL calculation.
As written in eq. In MC@NLO, where PS m N +1 denotes the parton shower approximation to B N +1 for channel m as determined by the splitting factors used in an actual parton shower algorithm like Herwig or Pythia, C m N +1 could be used as an NLO subtraction for B m N +1 , and the purpose of G(T N ) is to smoothly join the two. [In principle, G(T N ) ≡ G m N +1 (Φ N +1 ) can depend on m and the full Φ N +1 .] Note that the value of S N +1 for T N < T cut N was not needed in the LL and LO+LL discussions, but is needed here and the expressions we use are specific to the NLO+LL construction. In our formulation of eq. (3.7), the MC@NLO method corresponds to taking G(T N > T cut N ) = 1, since an actual parton shower is used to generate the Sudakov factor and T cut N is identical to the parton shower cutoff. The condition lim T N →0 G(T N ) = 0 is necessary to ensure that all IR divergences cancel in the limit T N → 0, because PS N +1 does not provide a valid NLO subtraction.
Even though there is no explicit T cut N dependence in eq. (3.9), the fact that PS N +1 does not reproduce the full IR singularities of B N +1 causes an implicit logarithmic sensitivity to scales ≤ T cut N in dσ S ≥N . To see this, we rewrite (3.12) The first three terms are IR finite and T cut N independent. The last term is also IR finite since lim T N →0 G(T N ) = 0. However, since G(T N > T cut N ) = 1, the subleading singular dependence in PS N +1 − C N +1 is integrated down to T cut N and only cut off below, which means this last term scales as O cut (α 1/2 s ). 8 Taking into account this implicit T cut s ), while differential (N +1)-jet observables are only accurate to 1+O(α s )+O cut (α 1/2 s ). Formally, this is not sufficient to maintain the perturbative accuracy expected at NLO N and LO N +1 , cf. table 1. In practice, the numerical impact depends on how well the employed parton shower algorithm is able to capture the subleading singular structure of the full real emission contribution. In refs. [10,11], this was shown to be a minor problem.
In Powheg, S N +1 is constructed by dividing the full B N +1 between the IR singular regions for the different emission channels, The conditions imposed on the Θ m N +1 ensure that the full B N +1 is obtained in any singular limit, such that S N +1 reproduces the full IR-singular structure and dσ S ≥N is IR finite. The function F (T N ) is included so the resummation can be turned off by letting can depend on m and the full Φ N +1 .] In this case, since S N +1 contains the full singular structure also above T cut N , there is no implicit T cut N dependence. Strictly speaking, this is true as long as Θ m and F do not introduce a sensitivity to small T N . The full Φ N +1 dependence in dσ mc ≥N +1 in eq. (3.7) is determined by S N +1 (Φ N +1 ) in the resummation term, i.e., by the approximate Φ N +1 dependence in the splitting factor that determines the Sudakov factor. The FO matching correction, dσ B−S ≥N +1 ∼ (B N +1 − S N +1 )(Φ N +1 ), additively corrects the approximate Φ N +1 dependence in S N +1 to the full LO N +1 dependence given by B N +1 . Another possible approach is to also multiply this term by the Sudakov factor, or equivalently, directly use the full B N +1 dependence in the resummed 8 The T cut N dependence becomes explicit if one takes G(TN > T cut N ) = θ(TN > T cut N ), in which case the integral would produce an explicit ln T cut N . For a smooth G this logarithm is smeared out but the integral has the same scaling. spectrum, such that This corresponds to the usual CKKW procedure for LO N,N +1 +LL in eq. (3.4). It is also analogous to the Geneva method in ref. [22], where the Φ N +1 -differential FO calculation is multiplicatively combined with the T N spectrum resummed to higher order. In eq. (3.14), the spectrum is not the exact derivative of the cumulant anymore, resulting in a residual T cut N dependence in the integrated cross section. The effective correction term by which eq. (3.3) is violated and that gets added to the correct NLO N cross section is given by In fixed order this is O(α 2 s ) and beyond NLO N . However, its impact on the perturbative accuracy depends again on the extent to which the IR singularities of B N +1 are correctly reproduced by S N +1 . If S N +1 contains the full IR singularities, so B N +1 − S N +1 is finite for T N → 0, then the leading term in eq. (3.15) scales as T cut . Therefore, in this case the correction can be regarded as a power correction. If S N +1 does not reproduce the full IR singularities, so that B N +1 −S N +1 contains subleading divergences ∼ α s /T N , then the leading term scales as α 2 s ln 3 (T cut N /Q). Hence, in this case the correction is of O cut (α 1/2 s ) and clearly violates the NLO N +LL accuracy, which allows at most O cut (α 2 s ) corrections (see the first column of table 1). Note that the perturbative accuracy of the residual T cut N dependence in either case here is the same as in eq. (3.4) at LO N,N +1 +LL. The reason is that it is determined by the resummation counting and the NLO matching by itself only improves the FO accuracy.

Combining NNLO calculations with LL resummation
As we saw in section 2.2, at NNLO we need events representing N , (N + 1), and (N + 2) partonic jets, defined through the N -jet and (N + 1)-jet resolution variables T N and T N +1 . The same is therefore also the case at NNLO+LL. Hence, we need to construct expressions for the corresponding fully differential MC cross sections [see eqs. (2.13) and (2.14)] As discussed in section 3.1, at NNLO+LL we require that N -jet observables are correct to NNLO N +LL, (N + 1)-jet observables to NLO N +1 +LL, and (N + 2)-jet observables to LO N +2 +LL, provided that any observable built from these cross sections is sufficiently inclusive over the unresolved regions of phase space. Since the FO calculation is supplemented with the LL resummation of the jet resolution variables T N and T N +1 , the perturbative accuracy of the prediction in the IR-singular regime is improved relative to the pure FO calculation, which breaks down in this region. The required perturbative accuracy at NNLO+LL in the FO and resummation regions is summarized in table 1.
To construct the NNLO+LL MC cross sections, it will be convenient to proceed in two steps. In section 4.1, we first consider the separation between the exclusive N -jet and inclusive (N + 1)-jet cross sections using T N and construct the corresponding exclusive dσ mc N (T cut N ) and an inclusive dσ mc ≥N +1 (T N > T cut N ). In section 4.2, we then consider the further separation of dσ mc ≥N +1 (T N > T cut N ) into the final exclusive dσ mc N +1 (T N > T cut N ; T cut N +1 ) and inclusive dσ mc ≥N +2 (T N > T cut N , T N +1 > T cut N +1 ) using T N +1 . To make the notation as transparent as possible, we drop the emission labels m throughout this section. They can be inserted straightforwardly into all formulae giving the different contributions to the cross sections.

The Exclusive N -jet and Inclusive (N + 1)-jet Cross Sections
As we have already seen at LO and NLO, it is convenient to divide the full FO exclusive N -jet cross section, dσ FO N (T cut N ), into a singular and a nonsingular contribution 9 , At NNLO, dσ FO N (T cut N ) is given in eq. (2.17). Its singular approximation is given by where C N +1 , V C N +1 , and C N +2 reproduce the exact IR singularities of B N +1 , V N +1 , and B N +2 , respectively, i.e., they correspond to a valid set of NNLO subtractions, such that eq. (4.3) is IR finite. The full logarithmic T cut N dependence arises from integrating B N +1 , V N +1 , and B N +2 , over the IR-singular region, which is fully reproduced by the C N +1 , V C N +1 , and C N +2 contributions in eq. (4.3). Therefore, dσ C N (T cut N ) contains all logarithms in T cut N , while the remainder dσ B−C N (T cut N ) in eq. (4.2) is a power correction in T cut N . To identify the relevant terms, we rewrite the N -jet MC cross section in terms of a resummed contribution and FO matching corrections. As we have seen at NLO+LL in section 3.3, the LL resummed contribution can be obtained by multiplying an inclusive cross section by the LL Sudakov factor for T cut N . The resulting expression in general differs from 9 To be precise, singular terms in the cumulant contain logarithms of T cut N or constants, while nonsingular terms vanish as T cut N → 0. In the spectrum, singular terms contain plus distributions or delta functions of TN , while nonsingular terms contain no singular distributions and at most integrable singularities. the correct FO result by both singular and nonsingular terms in T cut N , which are accounted for by adding corresponding FO singular and nonsingular matching corrections. This gives The first term is the resummed contribution, where dσ C ≥N is the singular approximation of the inclusive FO N -jet cross section, obtained by dropping the θ(T N < T cut N ) in eq. (4.3). It is by construction T cut N independent, so all dependence on T cut N in the resummed term resides in the Sudakov factor ∆ N (Φ N ; T cut N ), which sums the LL series in T cut N . The remaining two terms are FO matching corrections to ensure the correct FO expansion of eq. (4.4).
The last term in eq. (4.4), labeled B − C, is the FO nonsingular term from eq. (4.2). It contains the difference between the full FO contribution and its singular limit, As discussed above, it contains no logarithmic dependence on T cut N . The second term in eq. (4.4), labeled C − S, is the singular FO matching correction. It contains the difference between the singular approximation containing the full logarithmic T cut N dependence and that obtained by expanding the resummed term in fixed order, i.e., Hence, it supplies the FO singular terms in T cut N that are not contained in the resummed contribution. In the second line we show the NLO result for illustration. As already discussed in section 3.3, since the splitting function S N +1 generically only reproduces the leading singularities in C N +1 , dσ C−S N (T cut N ) can in general contain logarithmic dependence as large as α s L cut at NLO and α 2 s L 3 cut at NNLO, which contribute at O cut (α 1/2 s ) with the counting of eq. (2.28).
A potential problem with implementing eq. (4.4) is the presence of explicit logarithms in dσ C−S N (T cut N ), which become large as T cut N is reduced, and in particular dσ C−S N (T cut N ) diverges for T cut N → 0. While by construction this divergence cancels in physical observables, it could give rise to events with large or even negative weights. To circumvent this and regulate the logarithmic divergence, we can alternatively choose to multiply the singular matching terms with the Sudakov factor and write where the FO singular matching corrections are now given by Note that while multiplying with the Sudakov factor helps to suppress the FO T cut N logarithms in d σ C−S N (T cut N ), this choice does not amount to an actual resummation of these logarithms. A downside of this choice is that it introduces a more complicated T cut N dependence at all orders that must be canceled in inclusive N -jet observables. Since d σ where the various ingredients are discussed in detail in sections 4.1.1 and 4.1.2. For case 1, the FO singular and nonsingular matching terms are pure FO corrections and to obtain them it is sufficient to enforce that dσ mc ≥N +1 expands to the correct NLO cross section. For case 2, the singular matching correction is more complicated, and its T N dependence is obtained by taking the derivative of d σ C−S N (T cut N ) ∆ N (T cut N ) in eq. (4.7) with respect to T cut N . This ensures that the singular matching corrections in the spectrum correctly integrate up to cancel the corresponding T cut N dependence in the cumulant. 10 Before we give the detailed expressions for all ingredients required to construct eqs. (4.4), (4.7), (4.9), and (4.10), it is instructive to see how the NLO+LL case arises from this notation. At NLO, we have The corresponding results for the differential spectrum are . For MC@NLO, the splitting function is given in eq. (3.11). It depends on a function G(T N ), which for the sake of illustration we can choose as G(T N ) = θ(T N > T cut N ) (even though this is not the choice made in the MC@NLO implementation). In this case, the expression for dσ S ≥N given in eq. (3.12) is equivalent to dσ S ≥N = dσ C ≥N + d σ C−S N , which corresponds to case 2 in eq. (4.7) for the cumulant. However, the corresponding spectrum in eq. (3.7) is not that of case 2 in eq. (4.10). This is the origin of the residual T cut N dependence in MC@NLO discussed below eq. (3.12).
It should be clear from the discussion so far that the expressions in eqs. (4.4) and (4.9) for case 1 or alternatively eqs. (4.7) and (4.10) for case 2 provide a completely general result for the FO+LL matching valid to any fixed order. The explicit NNLO+LL expressions are given in detail below in section 4.1.1 for case 1 and section 4.1.2 for case 2. Besides the 10 Notice that there might be points in ΦN+1 for which BN (ΦN ) = 0 due to either kinematical or PDF effects. To avoid that the ratio SN+1(ΦN+1)/BN (ΦN ) goes to infinity, one has to define SN+1 such that it vanishes for these points. This implies that the contributions from these phase space regions are contained in dσ C−S or d σ C−S . choice one has between the two cases, different implementations can be obtained by making different choices for the C N +1 , V C N +1 , and C N +2 contributions that are used to approximate the singular behavior of the full theory, as well as for the splitting function S N +1 that is used to define the Sudakov factor. This amounts to shifting nonsingular corrections or subleading logarithms between the resummed contribution and the FO matching corrections.

Case 1
Here, we use dσ mc N (T cut N ) as given in eq. (4.4), with its corresponding inclusive dσ mc ≥N +1 (T N > T cut N ) given in eq. (4.9), which we repeat here for completeness: The explicit expressions for all ingredients are given in the following. By construction these are correct to NNLO N and NLO N +1 and include the correct LL resummation for T cut N and T N , respectively. Also, each of the three terms in the cumulant and spectrum separately satisfy the exact consistency relations in eqs. (3.2) and (3.3) without any residual T cut N dependence. The singular inclusive cross section, dσ C ≥N , appearing in the resummed terms is obtained by removing the constraints on T N in eq. (4.3), which gives (4.14) Since dσ C ≥N is explicitly T cut N independent, the resummed terms satisfy eq. (3.2) because [see eq. (2.30)] (4.16) The differential equivalent dσ B−C ≥N +1 (T N > T cut N ) is defined exactly analogously, (4.17) and one can easily see that eqs. (4.16) and (4.17) explicitly satisfy the consistency condition in eq. (3.3). Finally, the singular matching corrections, dσ C−S , are defined as By definition they satisfy eqs. (3.2) and (3.3), because each of the terms on the right-hand sides do so. To obtain their explicit expressions we use the NNLO expansion of the Sudakov factor, which we write as Here, we used S (n) N +1 to denote the α n s contribution to S N +1 , i.e., For convenience, we also define the subtracted one-loop virtual correction, which is the IRfinite NLO term in dσ C ≥N , The differential version is easier to obtain (since it does not explicitly require ∆ (2) N ), and we find The cumulant version is given by The integrals here are explicitly over T N > T cut N , which cuts off all IR singularities that do not cancel between the full FO singular contributions and their LL approximation arising from the Sudakov expansion, which is given by the last lines in eqs. (4.22) and (4.23). Note that C N +2 here fulfills two roles. First, it produces the leading double logarithms α 2 s (L 4 cut + L 3 cut ) (for the cumulant). The α 2 s L 4 cut is always canceled by the square [∆ N produces the correct single logarithm α s L cut at NLO. Second, the (N + 1)-parton virtual IR divergences in V C N +1 are canceled by the T N +1 → 0 limit in the Φ N +2 integral over C N +2 , where the remainder is an α s (α s L 2 cut + α s L cut ) correction. Generically, these are only partially canceled by the corresponding V C N ∆ (1)

Case 2
For this case, we use dσ mc N (T cut N ) as given in eq. (4.7), with its corresponding inclusive dσ mc ≥N +1 (T N > T cut N ) given in eq. (4.10), which we repeat here for completeness: The explicit expressions for all ingredients are given in the following. As for case 1, these are correct to NNLO N and NLO N +1 and include the correct LL resummation for T cut N and T N , respectively. The resummation terms involving dσ C ≥N ∆ N and the nonsingular FO matching terms, dσ B−C , are the same as in case 1 [see eq. (4.14) and eqs. (4.16) and (4.17)] and separately satisfy the consistency relations in eqs. (3.2) and (3.3).
The difference to case 1 is how the singular matching corrections, d σ S−C , are included. For the cumulant, we have where dσ C−S N (T cut N ) is given in eq. (4.22). The corresponding differential result in the spectrum is obtained by requiring eq. (3.3), where dσ C−S ≥N +1 (T N > T cut N ) is given in eq. (4.23). One can easily check that with this result the expression for dσ mc ≥N +1 in case 2 expands to the correct NLO N +1 result.

The Exclusive (N + 1)-jet and Inclusive (N + 2)-jet Cross Sections
The inclusive (N + 1)-jet MC cross section is divided into the exclusive (N + 1)-jet and inclusive (N + 2)-jet MC cross sections using a resolution scale T cut N +1 , dσ mc Note that this is just a special case of the consistency condition in eq. (3.3) applied to T N +1 and taking T c N +1 ≡ T max N +1 . The inclusive dσ mc ≥N +1 already resums the leading logarithms of T N in the (N + 1)-parton phase space. On top of that, we also want to resum the leading logarithms of T cut N +1 and T N +1 appearing in dσ mc N +1 (T cut N +1 ) and dσ mc ≥N +2 (T N +1 ). The LL resummation for T N +1 is obtained using the (N + 1)-parton Sudakov factor, ∆ N +1 , which is defined as where the upper limit on the integration over T N +1 should be chosen of order T N . Note that the (N +1)-parton splitting function S N +2 enters in the Sudakov factor relative to the (N +1)parton Born matrix element B N +1 , which is required to correctly sum the logarithms of T N +1 across the whole range of T N , even for T N ∼ T max N . In terms of the resummation accuracy, achieving (N)LO N +1 +LL implies that the (N + 1)-parton Sudakov factor must multiply the complete B N +1 matrix element to obtain the LL resummation of T N +1 (or T cut N +1 ) in the limit Given these considerations, we again divide the exclusive (N +1)-jet and inclusive (N +2)jet MC cross sections into a resummed contribution and FO matching corrections, This has precisely the structure of the usual NLO N +1 +LL calculation [see eq. (3.7)], but with the dependence on the singular and nonsingular FO matching corrections, dσ C−S and dσ B−C , written out explicitly. Furthermore, dσ C ≥N +1 (T N > T cut N ) is the singular approximation to the full (N + 1)-jet inclusive cross section on which the T N +1 resummation acts. The crucial difference compared to the usual NLO+LL case discussed in section 3.3 is that the NLO N +1 +LL calculation is used down to very small values T N > T cut N , and so dσ C ≥N +1 (T N > T cut N ) now has to include the LL resummation in T N . In terms of the inclusive dσ mc ≥N +1 (T N > T cut N ) [given by either eq. (4.9) or eq. (4.10)] we can write it as where the second term on the right-hand side removes the dependence on B N +2 from dσ mc ≥N +1 , i.e., it removes the last line in dσ B−C ≥N +1 in eq. (4.17). By definition of C N +2 this term has no logarithmic dependence on T N , and therefore does not affect the LL resummation in T N .
Expanding this to fixed NLO N +1 reproduces the (N + 1) version of eq. (4.11), This shows that in the limit of turning off the T N resummation eq. (4.28) reproduces the correct NLO N +1 +LL result as required.
The FO matching corrections are determined by imposing the correct NLO N +1 and LO N +2 expansions of eq. (4.28). The nonsingular matching corrections are given as and (again by definition of C N +2 ) have no logarithmic dependence on T cut N +1 . For the singular matching corrections we then find Here, we can explicitly see the mismatch between the exact definition of T N (Φ N +2 ) required at NNLO N from the shower approximation in the S N +2 term, which inherits theΦ N +1 (Φ N +2 ) dependence from the projection from Φ N +2 to Φ N +1 in the (N + 1)-jet Sudakov factor. Generically, this can introduce a subleading logarithmic dependence on T cut N in dσ C−S (even in the limit S N +2 = C N +2 ), whose coefficient scales as ∼ T cut N . With the above results, we can check that no residual T cut N +1 dependence (beyond power corrections) is introduced in physical observables, because eqs. (3.2) and (3.3) are explicitly satisfied. For the FO matching corrections this is clear from their above expressions. The resummed terms combine correctly to the inclusive dσ C ≥N +1 using the equivalent relation to eq. (2.29) for the (N + 1)-parton Sudakov, (4.33) Using this relation, we can also easily check that eq. (4.26) is satisfied. Upon integration over dΦ N +2 /dΦ N +1 the dσ C−S N +1 and dσ C−S ≥N +2 terms cancel each other, while the dσ B−C N +1 and dσ B−C ≥N +2 terms combine to precisely cancel the second line in eq. (4.29). Hence, we precisely get back dσ mc ≥N +1 (T N > T cut N ), which shows that no residual T cut N dependence is introduced. In the above construction we have the same amount of freedom as in section 4.1 in how to implement the T N +1 resummation and where to put the FO singular corrections. Above we have used the analog of case 1 from section 4.1, where dσ C−S is included at fixed order. Various alternatives are: • One can multiply dσ C−S N +1 by the ∆ N +1 Sudakov, analogous to case 2 in section 4.1. In this case, eq. (4.26) is maintained exactly when the corresponding case 2 version is also used for the differential spectrum.
• One has the freedom in eq. (4.29) and all the results following it to use a different C N +2 than the C N +2 used in section 4.1. This includes whether one uses T N (Φ N +2 ) or T N (Φ N +1 ) to implement the T N > T cut N constraint for the C N +2 contribution. In particular one could use a simpler NLO N +1 subtraction for C N +2 . (In general this can change the logarithmic dependence on T N at the subleading level.) • One can use different choices for S N +2 . In particular, in conjunction with using an alternative C N +2 , one can use a Powheg approach for NLO N +1 +LL such that one can take S N +2 = C N +2 .

Matching the NNLO+LL calculation with a parton shower
In the previous sections we have shown how to consistently combine LO, NLO, and NNLO calculations with LL resummation, and to obtain the MC cross sections dσ mc N , dσ mc N +1 , and dσ mc ≥N +2 . In this section, we discuss how to interface the corresponding N -parton, (N + 1)parton, and (N + 2)-parton events with a parton shower. The resulting NNLO+LL event generator will thus be able to produce events with any parton multiplicity.
The NNLO+LL MC cross sections of section 4 provide resummation in the resolution variables T N and T N +1 , but in general do not explicitly resum large logarithms arising in singular regions of phase space for other observables. In the resummation regime the shape of a generic exclusive observable will therefore only be accurately predicted after the addition of the parton shower, which in general provides LL accuracy. Furthermore, care must be taken when interfacing to the parton shower such that the perturbative accuracy provided by the MC cross sections dσ mc M is maintained. This includes their FO accuracy, the LL accuracy in the evolution variables, and the absence of residual dependence on the resolution scales T cut N and T cut N +1 . Precisely, the matching with the parton shower must satisfy three conditions: 1. Any exclusive observable must be correct to at least LL in the resummation regime. This includes the resolution variables T N and T N +1 , for which the LL accuracy of the MC cross sections must be maintained. Additionally, the LL accuracy requirement extends to observables requiring more than N + 2 jets, for which the parton shower provides the only prediction.
2. The FO accuracy of any observable should be that of the NNLO calculation (see section 2.2), which means: • N -jet observables are correct to NNLO N up to power corrections of relative order O(α s T cut and T eff N are the effective resolution scales to which the observable is sensitive. • (N + 1)-jet observables are correct to NLO N +1 if they only include contributions in the resolved region of Φ N +1 , up to power corrections of relative order is the effective resolution scale to which the observable is sensitive.
• (N + 2)-jet observables are correct to LO N +2 if they only include contributions in the resolved region of Φ N +2 .
Note that no FO accuracy is implied for observables sensitive to the unresolved regions of phase space, T N < T cut N and T N +1 < T cut N +1 , as the parton shower provides the only prediction in these regions (see below).
3. For observables that must be correct to N n LO any residual dependence on the resolution scales T cut N and T cut N +1 must enter at O cut (α ≥n+1 s ).
The conditions above naturally echo those imposed on the MC cross sections in section 3.1. In fact, in cases where the parton shower yields events with ≤ N + 2 partons, the exact phase space constraints implemented by the MC cross section definitions can be used on the shower (see figure 1). In cases with more emissions, one must develop analogous constraints making sure the above conditions remain satisfied.

LL shower constraints
Condition 1 above requires us to maintain the LL accuracy of the event sample and combine it with the parton shower LL resummation for additional emissions. For this purpose, the identical considerations apply to our NNLO+LL calculation as in the case of interfacing a merged LO N,N +1,N +2 +LL calculation with a parton shower [1][2][3][4][5][6][7][8][9]. The reason is that as far as the LL structure is concerned, the only relevance of the higher FO accuracy in our case is that it imposes a tighter constraint in condition 3 above. However, since the parton shower is formulated such that the probability of an emission is the exact differential of the no-emission probability [i.e. of the Sudakov factor, see eq. (2.30)], condition 3 will be satisfied as long as any additional constraints imposed on the parton shower do not spoil this relation.
The simultaneous LL resummation of T N and T N +1 in the NNLO+LL calculation can be achieved by choosing both variables to be equivalent (at the single-emission/LL level) to the same local shower evolution variable T [see eq. (2.20)], in which case we can assume that they are ordered as T N +1 < T N .

Equivalent resummation and shower evolution variables
The simplest case is when the evolution variable of the parton shower is equivalent to T (i.e. it has the same LL structure). The event sample with N , N + 1, and N + 2 partons can then be viewed as the result of the first two steps in the normal parton shower evolution in T , and attaching the parton shower simply corresponds to continuing this evolution down to the shower cutoff, where the relevant starting scale, T res , is given by the scale of the last emission or the resolution scale, namely • T res ≡ T cut N for the N -parton events • T res ≡ T cut N +1 for the (N + 1)-parton events Note also that in principle one can choose T cut N = T cut N +1 to be equal (or very close) to the actual shower cutoff T cut , such that no (or very few) additional emissions need to be generated for the N -jet and (N + 1)-jet samples.

Different resummation and shower evolution variables
If the local evolution variable T of the parton shower differs in its LL structure from the variable T used to implement the LL resummation in the partonic FO+LL calculation, one has to utilize a veto procedure on the shower to achieve condition 1. In principle, two approaches may be used here, using either a vetoed shower algorithm or a global veto procedure. Additionally, one has to specify the starting scale of the shower evolution.
The use of a vetoed parton shower was discussed in detail in refs. [1,12], for the case where T is the p T of an emission and using an angular-ordered parton shower where T is the emission angle. The same veto procedure can be applied here. The vetoed shower works by evolving in T and in each emission step only emissions satisfying the constraint T < T res are allowed, where T res is given as above. If an emission at some T violates this constraint, it is vetoed and the evolution continues from T . This vetoed shower exponentiates the T < T res constraint, which effectively transforms the shower evolution variable from T into T .
In the global veto procedure one lets the evolution proceed undisturbed. After the showering is done, the showered event is accepted if the condition T < T res is satisfied for all emissions. If this is not the case, the showering is repeated from the start on the same partonic event, and this is done until an acceptable showered event is generated. This second approach is certainly less efficient but it has the advantage that one does not need to modify the parton shower algorithm at all.
In either vetoing approach one has to choose appropriate starting scales for the T evolution. First, one determines the maximal starting scale T max , which should be either the value T max (Φ N ) that one would normally choose when starting the shower directly from B N (Φ N ), or the maximum value of T kinematically allowed for a given T res , whichever is smaller. The simplest approach is then to start the shower for all partons at T max . A somewhat better approach is to choose the starting scale according to the emission history. 11 For partons that had no emissions the shower is started at T max . For the daughter partons of an extra emission step in the (N + 1)-jet and (N + 2)-jet samples, the shower is started from the scale T res of the emission. The possible additional emissions for T max > T > T res are then added by running a truncated shower [12] from T max to T res along the parent parton line of the emission.

FO shower constraints
The constraints on the shower implied by condition 2 are simpler for event samples with higher jet multiplicity, as the desired perturbative accuracy is lower. Therefore, we start by discussing the (N + 2)-jet, working our way down to the N -jet sample. Note that if the shower evolves directly in T and both T cut N and T cut N +1 are set to the shower cutoff, only the (N + 2)-jet sample gets showered, and the additional complications arising for the (N + 1)-jet and N -jet samples become irrelevant.

Showering the (N + 2)-jet event sample
The MC cross section dσ mc ≥N +2 of the NNLO+LL calculation is given in eq. (4.28). Its perturbative accuracy is LO N +2 +LL, which the parton shower can easily maintain by applying constraints analogous to those applied to the highest jet multiplicity in a LO+LL matched event sample. The LO N +2 accuracy of the cross section is automatically guaranteed by the fact that additional emissions from the parton shower are higher order in α s . Therefore, there are no additional FO constraints on the shower. (Strictly speaking, the showered events in 11 The LL resummation in TN and TN+1 is formulated as a consecutive sum over emission channels m when splitting from N to N + 1 partons (in the construction of dσ mc ≥N +1 ) and from N + 1 to N + 2 partons (in the construction of dσ mc ≥N +2 ). Hence, we can naturally associate each contribution in this sum with an emission history for going from the underlying ΦN to the final ΦN+1 or ΦN+2 point. this sample must still satisfy the constraints T N > T cut N and T N +1 > T cut N +1 . If T N +1 < T N , ignoring this gives rise to at most power corrections.)

Showering the (N + 1)-jet event sample
The MC cross section dσ mc N +1 (T N > T cut N ; T cut N +1 ) of the NNLO+LL calculation is given in eq. (4.28). It contains the integrated cross section for T N +1 < T cut N +1 calculated to NLO N +1 +LL. Before adding the parton shower, it is represented by (N + 1)-parton events, which have T N +1 = 0 (see figure 1). By adding emissions, the parton shower distributes the events located at T N +1 = 0 to nonzero T N +1 values. In doing so, it must respect the exclusive (N + 1)-jet definition of the cross section, i.e., the cross section for T N +1 < T cut N +1 after showering has to remain accurate to NLO N +1 +LL. Since the parton shower preserves the total cross section, this means it is only allowed to fill out the region 0 ).] At LL accuracy, this is achieved by vetoing shower emissions with T > T cut N +1 , as discussed in section 5.1. In addition, to satisfy condition 2 it is also necessary that the cross section for T N +1 < T cut N +1 remains correct to NLO N +1 . The veto on single emissions with T > T cut N +1 is sufficient for this purpose as well, so we do not require an additional constraint on the shower. To see this, consider the shower emission with the largest value of T and sum over all other emissions. Strictly speaking we need the emission to satisfy , where Φ rad is the emission phase space andΦ N +2 is the inverse of the phase space projection Φ N +1 (Φ N +2 ) that is used in the NLO N +1 calculation. The single-emission veto in the shower corresponds to imposing the constraint is the phase space map used in the parton shower. In principle, the two constraints are different, since the two phase space maps can be different. However, both maps have to be IR safe and must agree in the IR limit T cut N +1 → 0. Therefore, the difference can be at most a power correction in T cut N +1 . From this discussion it follows that a generic (N + 1)-jet observable receives at most power corrections from showering of O(α s T cut N +1 /T eff N +1 ), where T eff N +1 is the effective scale that the observable is sensitive to. Similarly, since dσ mc ≥N +1 contributes at O(α s ) to generic N -jet observables, they receive at most power corrections of O(α 2 s T cut N +1 /T eff N +1 ). Hence, condition 2 is satisfied. In fact, as long as the T cut N +1 value is kept small, the spectrum for T N +1 < T cut N +1 is correctly described by the shower. The parton shower therefore improves the description of the previously unresolved region T N +1 < T cut N +1 . As a result, the power corrections induced by the shower actually compensate for the power corrections in the partonic calculation arising from the unresolved region below T cut N +1 . Of course, this is only true if the shower cutoff is lower than T cut N +1 .

Showering the N -jet event sample
The MC cross section dσ mc N (T cut N ) of the NNLO+LL calculation is given in eq. (4.4) or eq. (4.7). It contains the integrated cross section for T N < T cut N calculated to NNLO N +LL, which before p cut T p cut T Figure 3. Illustration of the issues in defining an IR-safe phase space separation at NNLO using single-parton variables in case of vector boson production. Limiting each emission to be below p cut T (dashed lines) results in a miscancellation of IR divergences between the tree-level contribution on the left, which would contribute to dσ mc 0 (p cut T ), and the corresponding one-loop contribution on the right, which would contribute to dσ mc ≥1 (p T > p cut T ).
showering is represented by N -parton events with T N = 0. The basic considerations here are similar as for the (N + 1)-jet case. Repeating the discussion in section 5.2.2, the shower must be constrained to not change the cross section for T N < T cut N , but to only fill out the T N spectrum below T cut N . Since the action of the parton shower is entirely within the N -jet cumulant bin, the induced power corrections of O(α s T cut N /T eff N ) are again at the level allowed by condition 2, and will actually improve the predictions of observables, because the unshowered events at T N = 0 are distributed over the previously unresolved region T N < T cut N with an LL-accurate shape. There is a further complication however, that arises starting at NNLO. At NLO+LL, the resolution variable must have two properties: it must realize an IR-safe separation of the phase space at the level of a single emission and it must have an LL resummation. Because LL resummation arises from exponentiating independent emissions, these two properties are essentially one and the same. For example, in an NLO+LL calculation of vector boson production, the resolution variable separating events with 0 jets and 1 jet can be chosen as the transverse momentum of the leading parton, with 0-jet events corresponding to p T < p cut T and 1-jet events corresponding to p T > p cut T . At NNLO+LL, however, the story is different: constraining the shower evolution in terms of independent single-parton variables is no longer sufficient to preserve IR safety in the separation of jet bins. To see how the problem arises, it is instructive to consider again the example of vector boson production with two emissions illustrated in figure 3. Demanding that the transverse momentum of each emitted parton is below p cut T (dashed lines) does not yield an IR-safe definition for the 0-jet cross section. If the two partons are collinear to each other and each satisfies p T > p cut T , this IR-divergent contribution would be included in the 0-jet cross section, while the corresponding IR-divergent virtual diagram on the right would contribute to the 1jet cross section. As already discussed in section 2.2, we must use a resolution variable which is properly IR-safe at NNLO. For example, we can sum over all emissions (T N = p T ), or combine them using an IR-safe jet-clustering procedure (T N = p jet T ). From this discussion, it is clear that the constraint T N < T cut N that the parton shower needs to satisfy, cannot be formulated in terms of individual emissions but must take at least two emissions into account. Generally, it is not sufficient to only consider the two hardest emissions, since they do not necessarily give the hardest jet. Therefore, the NNLO constraint can only be imposed via a global veto after the showering. In case one uses a vetoed shower with a single-emission local veto to enforce the LL constraints as described in section 5.1, the additional NNLO constraint should be enforced separately.

Implementation and relation to existing approaches
In this section, we discuss the relation of our framework to recent related work, and the NNLO+PS implementation given in ref. [38]. This will show that our method is indeed quite general and encompasses these other approaches. It also illustrates that an actual implementation of our results is indeed feasible.

GENEVA
The motivation to build an NNLO+LL event generator is to interface the most precise FO calculations available with a parton shower routine to be able to simulate realistic events with high perturbative accuracy. Whenever higher logarithmic resummation is also available (NLL for several resolution variables, NNLL for certain resolution variables such as N -jettiness, and NNLL for select processes 12 ), it can be implemented to also improve the perturbative accuracy in the resummation region (see figure 2) following the Geneva approach [22].
If NNLL resummation is available, the resummation order matches the fixed NNLO accuracy in the sense that all NNLO singular terms are naturally included in the resummation. Hence, the FO singular matching correction vanishes, because the FO expansion of the NNLL resummed result reproduces the full NNLO singular corrections. The remaining contributions in the N -jet MC cross section can then be associated as follows: That is, the cross section takes the form of a traditional resummed calculation, with the FO nonsingular corrections corresponding to dσ B−C N and the higher-order resummed cumulant replacing the resummation term dσ C ≥N ∆ N (T cut N ). The same relations also apply for the exclusive (N + 1)-jet and inclusive (N + 2)-jet cross sections.
The results in ref. [22] took this approach, using a jet resolution variable for which higherorder logarithmic resummation is available. There, the NNLL resummation for e + e − → jets for small T 2 was used together with the NLO 2 nonsingular terms, combined with the fully differential 3-jet cross section at NLO 3 , and interfaced with a parton shower algorithm. As discussed above, the resummation to NNLL already incorporates the full singular contributions up to NNLO, including the two-loop virtual corrections. Thus, the only missing contributions to make the calculation in ref. [22] correct to full NNLO 2 are the nonsingular corrections at NNLO 2 . Since they scale as a power correction in T cut 2 , one could also take the value of T cut 2 small enough to make their numerical impact small.

NNLO+PS using HJ-MiNLO
Results combining the inclusive NNLO Higgs cross section with a parton shower algorithm were presented recently in ref. [38]. This approach uses the Multi-Scale Improved NLO (MiNLO) calculation for the production of Higgs in association with a jet [53], in which the Powheg HJ calculation [54] is supplemented by an analytic Sudakov resummation factor, which includes logarithmic terms that become large as the transverse momentum of the Higgs boson tends to zero. The Sudakov factor effectively regulates the divergences in the Powheg HJ calculation when the transverse momentum of the Higgs boson, q T , goes to zero. As a result, the HJ-MiNLO sample can be used over the whole phase space even in the limit q T → 0. In practice, it is used down to q T of order Λ QCD ∼ 1 GeV.
It was shown in ref. [24] that by explicitly including NNLL information in the Sudakov factor, the HJ-MiNLO cross section integrates up to the correct inclusive Higgs cross section at NLO 0 . The HJ-MiNLO sample is then reweighted to the differential NNLO 0 Higgs cross section, which is facilitated by the fact that it is only single-differential in the Higgs rapidity. This provides NNLO 0 accurate predictions for 0-jet observables without spoiling the NLO 1 accuracy of 1-jet observables. One feature of this approach is that it does not require a Higgs + 0-jet sample, since the full NNLO 0 information of inclusive Higgs production is explicitly included through the reweighting factor.
While this approach seems at first sight quite different from the discussion in this paper, we will now show that it directly follows as a special case from our results in section 4. Hence, it can be viewed as a specific implementation of the general method developed in this paper. We first write the results of ref. [38] in terms of the MC cross sections dσ mc 0 (T cut 0 ) and dσ mc ≥1 (T 0 > T cut 0 ), corresponding to the exclusive Higgs + 0-jet and inclusive Higgs + 1-jet cross sections. We then show how these expressions follow directly from our general results by making specific choices.
The 0-jet resolution variable used in ref. [38] to separate 0 from 1 or more extra jets is the transverse momentum of the Higgs boson, so We do not need to discuss how to separate the inclusive 1-jet sample into an exclusive 1-jet and an inclusive 2-jet sample. For this purpose, ref. [38] uses the standard Powheg approach, which we have already shown in section 3.3 to be a special case of our approach. As mentioned already, the Higgs + 0-jet cross is not included in ref. [38], since it vanishes in the limit T cut 0 → 0. The inclusive MC cross section for one or more jets is then given by Here, the inclusive 1-jet cross section, dσ HJ-MiNLO

≥1
, is equivalent to the modifiedB function from HJ-MiNLO, which is obtained from the usualB function in Powheg by multiplying with the Sudakov factor ∆ 0 (T 0 ), and subtracting its first-order expansion to maintain the NLO 1 accuracy, The term in curly brackets contains the full singular T 0 dependence at NLO 1 . The crucial ingredient [24] is the fact that the exponent of the Sudakov factor ∆ 0 (T 0 ) contains the full NNLL set of T 0 logarithms to O(α 2 s ). This causes the spectrum to become the total derivative of the NLO 0 correct 0-jet cumulant, dσ NLO ≥0 ∆ 0 (T cut 0 ), up to nonsingular corrections in T 0 and higher orders in α s . As a result, the spectrum integrates to the correct NLO 0 cross section up to power corrections that vanish as T cut The reweighting factor R(Φ 0 , T cut 0 ) in eq. (6.4) is then given by the ratio and by construction ensures that the Higgs + 1-jet spectrum in eq. (6.4) integrates to the correct NNLO 0 inclusive Higgs cross section. At the same time, because of eq. (6.6), the reweighting factor has the form and therefore does not affect the NLO 1 accuracy of the inclusive 1-jet cross section up to power corrections in T cut 0 . By taking T cut 0 → Λ QCD these become negligible, and the result becomes a valid NNLO+LL implementation.
To derive this result as a special case from our framework, we make the following two choices: 1. Choose all singular terms equal to the exact tree-level and one-loop contributions, 2. Choose the splitting functions as With these two choices, the singular inclusive cross section defined in eq. (4.14) is given by the full NNLO 0 expression, while all FO matching corrections vanish, The choice of the splitting function S 2 (Φ 2 ) is not relevant for this discussion, since its purpose is to determine how to split the inclusive 1-jet cross section into an exclusive 1-jet and an inclusive 2-jet cross section.
Using the results of section 4.1.1 (or section 4.1.2, which are identical in this case), we then find for the exclusive 0-jet and inclusive 1-jet MC cross sections where in the last equation we inserted the explicit expression for S 1 (Φ 1 ) from eq. (6.10). We can now compare this to the HJ-MiNLO result in eq. (6.4). Since the exclusive 0-jet cross section is proportional to the Sudakov factor ∆ 0 (Φ 0 ; T cut 0 ), it vanishes in the limit T cut 0 → 0. Thus, in this limit the entire 0-jet cross section can be obtained by integrating the inclusive 1-jet result over all values of T 0 , precisely analogous to what happens in refs. [24,38]. Since in practice, T cut 0 ∼ Λ QCD ∼ 1 GeV, one could also keep the 0-jet cumulant, which would avoid introducing any additional power corrections in T cut 0 . The term in curly brackets times the Sudakov factor ∆ 0 (Φ 0 ; T 0 ) is equivalent to dσ HJ-MiNLO ≥1 /dΦ 1 in eq. (6.5), except for the additional V C 0 (Φ 0 ) term. By including this term, the prefactor in dσ mc ≥1 becomes simply the inclusive NNLO cross section normalized to the tree-level result, dσ NNLO ≥0 /B 0 (Φ 0 ), without any need to reweight the events.
With the choice C 1 (Φ 1 ) = B 1 (Φ 1 ) from above, V C 0 (Φ 0 ) is the NLO correction to the inclusive cross section [see eq. (4.21)], dσ NLO ≥0 dΦ 0 = B 0 (Φ 0 ) + V C 0 (Φ 0 ) , (6.14) and in particular T 0 independent. Although in principle there is no need to do so, we can rewrite dσ mc ≥1 and pull this term outside into the prefactor, which gives with the rescaling factor The last term in the denominator here is the O(α 3 s ) cross term that arises from pulling V C 0 (Φ 0 ) out into the rescaling factor. It must be kept because it scales as α 3 s (ln T 0 )/T 0 which upon integration over T 0 becomes an α 2 s correction. Equations (6.15) and (6.16) are now the exact equivalent of the expressions in eqs. (6.4), (6.5), and (6.7). By writing the factor in curly brackets in eq. (6.15) as 1 , one can easily check that the denominator in eq. (6.16) is exactly the integral of eq. (6.15) modulo the R(Φ 0 ) prefactor.
As we have seen, with the two choices given above our method gives an expression with an analogous structure as in ref. [38]. In fact, the result in eq. (6.13) that follows immediately from our approach is automatically correct to NNLO 0 without requiring an additional reweighting. Another difference is the precise form of the Sudakov factors, ∆ 0 (Φ 0 ; T 0 ) and ∆ 0 (Φ 0 ; T 0 ). In our approach, ∆ 0 is constructed from the splitting functions S (i) 1 (Φ 1 ), while in ref. [24] ∆ 0 is obtained from the analytic q T NNLL resummation formula. Both expressions have the same logarithmic dependence on T 0 expanded to O(α 2 s ) in the exponent. We also like to point out that in the approach of refs. [24,38] the known NNLL structure of the T 0 = q T spectrum is essential to analytically control all singular logarithms through O(α 2 s ). In this respect, this approach is thus closely related to the Geneva approach [22] discussed in section 6.1.

UNLOPS
In section 4, we have explicitly constructed the required exclusive N -jet and (N + 1)-jet MC cross sections to satisfy all the requirements to obtain a correct NNLO+LL event sample discussed in section 3.1. Alternatively, one could also start from the inclusive FO+LL M -jet cross sections and generate the exclusive MC cross sections numerically, This method has been applied to merge multiple NLO+LL calculations in refs. [23,41,55], where it is referred to as UNLOPS. Using eq. (6.17), the consistency conditions in eqs. (3.2) and (3.3) between different multiplicities is automatically enforced. The inclusive MC cross sections that are used as inputs must be correct at the relevant FO+LL accuracy according to eq. (3.1). For dσ mc ≥N this means it has to be correct to NNLO N , so it is simply given by the inclusive NNLO N cross section, dσ mc The inclusive (N + 1)-jet cross section must be correct to NLO N +1 with the T N dependence resummed to LL, and the inclusive (N + 2)-jet cross section must be correct to LO N +2 with the dependence on both T N and T N +1 resummed to LL, for which our general results in section 4 [see eqs. (4.9) and (4.28)] can be used. The major drawback of subtracting the integrals over the inclusive cross sections in eq. (6.17) numerically is that one has to generate events with negative weights. The advantage is that the expressions for the inclusive cross sections can be simplified substantially by dropping all higher-order dependence inherited from lower multiplicities. For the inclusive (N + 1)-jet cross section one could then use for example which includes the correct LL resummation and expands to the correct NLO N +1 result. One could also have written this result using a singular approximation to the inclusive cross section, and added a FO matching correction, or only have the Born-level result multiply the Sudakov factors, and then add all higher-order terms in the FO matching correction. This last choice corresponds to what is done in refs. [23,41,55]. For the inclusive (N +2)-jet MC cross section one could use the equivalent of the CKKW result, × ∆ N (Φ N ; T N ) ∆ N +1 (Φ N +1 ; T N +1 ) . (6.20)

Conclusions
In this paper we have developed a general method to combine fully differential NNLO calculations with LL resummation in the form of an event generator for physical events that can be directly interfaced with a parton shower. The basic quantities in our construction are Monte Carlo (MC) cross sections representing an exclusive partonic N -jet cross section, calculated to NNLO N +LL, an exclusive partonic (N + 1)-jet cross section, calculated to NLO N +1 +LL, and an inclusive partonic (N + 2)-jet cross section, calculated to LO N +2 +LL. We use N n LL M to refer to the O(α n s ) result relative to an M -parton tree-level result. These MC cross sections are represented in the generator by events with N , N + 1, and N + 2 partons. They are characterized by N -jet and (N + 1)-jet resolution variables T N and T N +1 , with resolution scales T cut N and T cut N +1 defining the separation between them. We stress that these are not jet-merging scales but IR cutoffs equivalent to a parton shower cutoff.
We have formulated the general conditions on the perturbative accuracy that a complete and fully differential NNLO+LL calculation must satisfy. They require that the MC cross sections must have the correct FO expansion (NNLO N for dσ mc N , NLO N +1 for dσ mc N +1 , and LO N +2 for dσ mc ≥N +2 ), as well as include the LL resummation of the resolution variables and scales (T cut N for dσ mc N , T N and T cut N +1 for dσ mc N +1 , T N and T N +1 for dσ mc ≥N +2 ). In addition, the consistent combination of FO and LL requires that all observables that are expected to be correctly predicted at O(α n s ) at fixed order must be independent of the resolution scales T cut N and T cut N +1 up to residual corrections of O cut (α ≥n+1 s ) [using the LL counting in eq. (2.28)] to maintain their expected perturbative accuracy. We have shown that this can be achieved in general by enforcing a derivative relationship between M -jet exclusive and (M + 1)-jet inclusive cross section.
Our main results are given in section 4, where we derive in detail the MC cross sections needed to construct the NNLO+LL event generator. The MC cross sections are explicitly given in terms of the constituent matrix elements used in FO calculations and the parton shower. Our results are general and we make no choices about the techniques used to evaluate the FO contributions in the MC cross sections. The primary and only NNLO ingredients that are required are a singular approximation of the inclusive NNLO N -jet cross section, dσ C ≥N , and the corresponding NNLO subtractions, both of which are naturally part of existing NNLO calculations. All other ingredients are NLO in nature, and therefore obtainable as in existing NLO+LL implementations. We proved that our construction explicitly satisfies all required conditions on the perturbative accuracy of an NNLO+LL event generator.
We have discussed how the partonic NNLO+LL event generator can be interfaced with standard parton showers using existing technologies, as well as the constraints that must be placed on the parton shower routine. This matching must preserve the FO and LL accuracy of the MC partonic jet cross sections, and the parton shower will provide LL accuracy for general N -jet, (N + 1)-jet, and (N + 2)-jet observables, producing events at all parton multiplicities. For the (N +1)-jet and (N +2)-jet samples, which are needed to NLO N +1 +LL and LO N +2 +LL accuracy respectively, the constraints are essentially the same as for the well-known case of NLO+PS matching. For the showering of the exclusive N -jet sample, which is needed at NNLO N +LL accuracy, we showed that the constraints on the parton shower can not be implemented at the level of individual emissions as was possible for the other multiplicities. However, a global veto on the parton shower can still be used in this case. Alternatively, if the shower evolution variable coincides with the T N and T N +1 resummation variables, the resolution scales T cut N and T cut N +1 can be set equal to the parton shower cutoff itself, in which case only the inclusive (N + 2)-jet sample must be showered.
Finally, we have discussed how other methods for matching higher-order perturbative calculations with parton showers fit into our general framework. For the well-known case of NLO+LL matching, the Powheg and MC@NLO approaches naturally follow as special cases. When employing the higher-order resummation at NNLL as in Geneva, the only missing ingredients to achieve full NNLO accuracy are power-suppressed nonsingular contributions. We have also shown explicitly how the recent results for NNLO+PS using HJ-MiNLO arise as a special case from our general results. We also commented how the ideas of UNLOPS fit into our method.
Our results provide a path for combining the precision frontier of fixed-order calculations with the flexibility and versatility of parton shower Monte Carlo programs. There are various steps that should be taken next toward a practical implementation. While the comparison to existing approaches makes it clear that the implementation is feasible, it remains to be seen what the optimal choices are to make the implementation sufficiently generic so that new NNLO calculations can be incorporated with limited effort. Finally, it should be clear from our discussion, that our general setup does not only apply to NNLO calculations, but can be extended to even higher order, should such results become available, though the details remain to be worked out in this case.