Matching N3LO QCD calculations to parton showers

The search for new interactions and particles in high-energy collider physics relies on precise background predictions. This has led to many advances in combining precise fixed-order cross-section calculations with detailed event generator simulations. In recent years, fixed-order QCD calculations of inclusive cross sections at N3LO precision have emerged, followed by an impressive progress at producing differential results. Once differential results become publicly available, it would be prudent to embed these into event generators to allow the community to leverage these advances. This note offers some concrete thoughts on ME+PS matching at third order in QCD. As a method for testing these thoughts, a toy calculation of $e^+e^-\rightarrow u \bar u$ at $\mathcal{O}(\alpha_s^3)$ is constructed, and combined with an event generator through unitary matching. The toy implementation may serve also as blueprint for high-precision QCD predictions at future lepton colliders. As a byproduct of the N3LO matching formula, a new NNLO+PS formula for processes with"additional"jets is obtained.


I. GENERAL INTRODUCTION
High-energy collider physics tries to provide insights into a consistent quantum field theory of nature by accumulating immense amounts of experimental data, and confronting it with precision calculations -the expectation being that enough data on as many scattering processes as possible might expose flaws in our current understanding. Unfortunately, since calculations in the Standard Model of particle physics quickly become tedious, hints of new interactions are often explained by phenomena that were omitted in the background calculation. This has resulted a large and vibrant sub-community producing precise and detailed calculations of background processes -especially for experiments at the Large Hadron Collider. The typical model for calculations used by experimentalists is a combination of higher-order QCD (or electroweak) fixed-order calculations (to obtain the best available prediction of the highestenergy scattering process), combined with an event generator handling the evolution of the high-energy particles into composite hadrons, energetic jets and other remnants of the colliding beams.
The first general methods to improve the precision and/or accuracy of such a computational model emerged at the turn of the century, when the matching of NLO QCD calculations with event generators [1,2] and the merging of multiple leading-order multi-jet calculations [3][4][5] were formulated. This kick-started steady progress of incorporating higher QCD orders and/or a more accurate treatment of the other SM interactions [6,7]. Apart from a few impressive examples [8], this progress has slowed in recent years, and is being complemented with improvements of all-order parton showers. Parton showers form the backbone of matching fixed-order calculations to event generators, by e.g. furnishing differential subtractions to make unweighted (or unweightable) fixed-order event generation by Monte-Carlo methods possible. Fully-differential matching procedures are presently only available in NLO QCD, since the (fully differential) singularity-structure of QCD is not yet captured by any available parton shower. This has however not prevented several successful combinations of NNLO predictions and event generators [7,[9][10][11]. Phase-space slicing inspired unitarized merging methods offer a convenient stepping stone towards high accuracy.
This note aims to present ideas towards matching N3LO QCD calculations to event generators. The overall philosophy of the approach is straight-forward and based on enforcing the desired target precision through all-order subtractions inspired by the unitarity of parton showers. This allows to combine several calculations executed with minimal phasespace cuts (to avoid the most singular regions) into a precise matched calculation, similar in spirit to e.g. q ⊥ [24] or Njettiness-sliced [25][26][27][28] fixed-order results. Section II will be used to set the scene for the third-order 1 method from a pedagogic angle, while sec. III explains the construction of a toy fixed-order calculation. In the absence of differential N3LO calculations that produce finite-weight events, it is convenient to work with a toy example to validate the third-order matching in a very controlled environment. Details about third-order matching are set out in sec. IV, and the suggestions are studied with the help of the toy calculation in sec. V. The note also contains detailed appendices with background information on parton showers and the implementation of matching factors.

II. INTRODUCTION TO MATCHING UP TO NNLO PRECISION
Integrating fixed-order calculations into event generators almost always requires removing the overlap between the parton shower approximation embedded in the event generator and the fixed-order result. This overlap can be removed by subtraction or by diving phase-space into disjoint fixed-order and shower regions. Both cases yield the benefit of allowing for event generation, i.e. the production and storage of finite-weight phase-space points allowing an analysis independent of the computational details. For the purposes of this note, it is convenient to take the perspective of a (conventional) parton-shower calculation. This entails that • n-parton phase-space points Φ n can be obtained from parton-shower evolution; the parton shower sums logarithmic enhancements in perturbation theory • the "separation" of partons is determined by the parton-shower ordering variable t calculated from their fourmomenta; a sequence of ordering variables (as is necessary to assign parton-shower factors) is defined through a parton-shower history from the lowest-to the highest-multiplicity state; if necessary to determine an ordering sequence, one amongst all possible histories is chosen probabilistically (cf. Appendix A 2); the entries in a sequence of ordering variables are called t(Φ 0 ), t(Φ 1 ) . . . t(Φ n ) or abbreviated t 0 , t 1 . . . t n • unitarity of the evolution holds, i.e. the action of the parton shower does not change inclusive cross sections; due to unitarity, parton-shower resummation assigns finite (or vanishing) cross sections even to phase-space points with collinear or soft partons; parton-shower all-order factors can be used to regularize fixed-order cross sections evaluated at phase-space points containing collinear or soft partons.
The aim of matching a fixed-order calculation to a parton shower is to combine the strengths of either approximation, while ensuring that the resulting combination is fixed-order accurate and does not impair the accuracy of the partonshower resummation. The term "accuracy of the parton shower" is not immediately obvious. The heavily constraining definition applied in this note is given in Appendix A 3. Before diving into specifics of matching, fig. 1 introduces the features of the notation used throughout this note.
f n (i) [M] (Φn) function of interest; typically, the differential cross section σ, or the Sudakov factor ∆ parton multiplicity; the number of additional partons is defined relative to the lowestmultiplicity process; e.g. e + e − → uū entails n = 0, e + e − → uūg means n = 1. . . perturbative order; the order is defined relative to the lowest order necessary to define the function f ; combinations may be used; e.g. (0 + 1) signifies the Born contributions and first-order corrections; a value ∞ signifies an all-order factor additional tag to further specify the function; e.g.
[inc] for "inclusive" or [tc] for "exclusive, with veto on real emissions above tc"; tags will be defined in the text.
phase space point Φ at which the function is evaluated.
FIG. 1: Reference for the sub-and superscript notation used in this note. Some examples of this notation that will appear repeatedly are given in Table II. Using this notation, the action F of the parton-shower on an ensemble of particles Φ n with distribution dσ (0) n (Φn) either leads to no change in any observable O (i.e. O will be still evaluated at the phase-space point Φ n , i.e. O = O(Φ n ) ≡ O n ), or to the decay of one or more particles: where P (t, z, φ) is the sum of naive parton-shower branching kernels 2 , and ∆ n (t + , t − ) is the Sudakov factor encapsulating the no-branching rate. The evolution of incoming partons typically introduces ratios of parton-distribution functions [29]. Such ratios are suppressed for the sake for readability in this note. The second equality 2 is due to the exponential form of the Sudakov factor, The notation w n+1 (Φ n+1 )∆ n (t + , t) is a symbolic representation of the actual parton-shower result, which includes a sum over all possible "paths" to arrive at the real-emission phase space point. Details of parton-shower factors are given in Appendix A 7. The main text will use the shorter, symbolic notation, to avoid over-crowding. Equation 3 shows that by construction, no interval [t + , t − ] will add to the overall cross section -simply because the integrated radiation pattern between any two scales is subtracted from the semi-exclusive lowest-multiplicity contribution, such that integrating over the last line will lead to the original distribution dσ (0) n (Φn). Viewed slightly differently, this "parton-shower unitarity" can also be interpreted as calculating a radiation pattern in the collinear approximation, regularizing the radiation pattern through all-order Sudakov factors (i.e. the third line), and subtracting the regularized radiation pattern from the next-lower multiplicity (i.e. the second line). This reasoning can be exploited to obtain a more precise calculation, by following the steps 1. Choose a precise target lowest-multiplicity prediction 2. Correct the radiation pattern to the desired accuracy, regularize with Sudakov factors and remove undesired higher orders introduced by this reweighting 3. Subtract the integrated form of the improved radiation pattern from the lowest-multiplicity contribution This "unitarized merging" approach [6] can be performed for an arbitrary t − value. The details of this procedure can become rather intricate -but are manageable, as elaborated on in Appendix A. The unitarized merging paradigm rests on the assumption that the parton shower does not change the inclusive cross section. Certain types of threshold enhancements are known to violate this condition [30], and are thus typically not included. However, an appropriate redefinition of the "target prediction" may circumvent this concern. As a first example of unitarized matching, it is possible to obtain an NLO-correct calculation from the previous result, by performing the operations 3 is the inclusive NLO cross section differential (only) in the variables Φ n . The cut Q n+1 > Q c is necessary to regularize a tree-level calculation of the real-emission configurations. The physical interpretation of the "inclusive cross section" should not be over-stretched, since it only provides a suitable model for observables that cannot resolve any effect of an additional real-emission parton, irrespective of its hardness. Only very few measurements in a realistic environment have this trait. However, inclusive cross sections serve, together with differential real-emission cross sections, as building blocks for more realistic predictions. The steps in Table I for example suggest that an NLO matched rate is given by From this form, it becomes clear that the most natural functional definition of the cut is Q = (parton shower evolution variable), making for the most natural value Q c = t − . In this case, the division into an O n -component and a F is -by virtue of eq. 2 -completely equivalent to the parton shower result. However, eq. 4 is still appropriate if Q n+1 = t n+1 and the parton shower has (effectively) vanishing support for Q n+1 < Q c . This means that the calculation is sub-divided into an n-parton contribution and an n+ ≥ 1 part. The latter are only included fully differentially above Q c , so that Q c should be chosen as small as possible. The presence of many high-energy n-parton contributions might lead to spuriously large hadronization effects, since these parton ensembles are directly hadronized, without dressing the states with further soft or collinear Pqq + Pqq + Pgg + Pgq , with the splitting kernels Pij defined in [31]. Unless expressly necessary, the abbreviated symbol ∆n(tn, tc) will be used synonymously.
n (tn, tc) : the mth term in the expansion of the parton-shower Sudakov factor pertaining to no emission off the particle ensemble Φn between the scales tn and tc.
: the combination of all parton-shower all-order factors except Sudakov factors that the shower would have applied to the ensemble Φn+1. For e + e − → jets annihilation, w where µ is a fixed value of the renormalization scale. If necessary in the interest of brevity, the notation wn+1(Φn+1) will be used synonymously.
: the mth term in the expansion of the expansion of the combination above.
: the inclusive NLO cross section to produce n additional partons, differential only in the variables Φn. The real-emission correction is integrated over the single-parton phase space.
The Bcross sections in the Powheg method [2] fall under this category.
Φn where I is the set of all singular limits of the real corrections. Since various treatments of singular limits are possible, this may lead to slightly different definitions of the inclusive cross section. The matching method developed in this note should be flexible enough to handle differing definitions.

Q(Φn) > Qc
: constraint on phase space points Φn: Φn should lead to Q-values above Qc (cf. Appendix A). The notation Qn ≡ Q(Φn) will also be used, for brevity.  radiation. Similar concerns may be raised in any event generator prediction that employs a parton-shower cut-off, since the no-emission contribution in eq. 1 also consists of high-energy few-parton states. In that case, though, no-emission events are typically rare, such that their hadronization only produces marginal effects. In eq. 4, two contributions that are numerically sampled with many events cancel to produce a small effective no-emission event count. This procedure may lead to a higher sensitivity to statistical outliers in the hadronization procedure, i.e. concerns about hadronization are of technical, not theoretical nature. Nevertheless, for jet observables that guarantee that states with no emissions or emissions with Q n+1 < Q c are treated inclusively, it is permissible to replace O n with the action of a parton shower F (∞,Qn+1<Qc) n (Φ n , t + , t − ) that omits emissions with Q n+1 > Q c . This might improve the transition to the non-perturbative regime. This note will continue to use the notation "O n ", leaving the possibility of such a replacement implicit. Further comments on this "sliced" or "binned" approach are given in Appendix B.
The extension to a unitarized NNLO matching scheme is possible once an NLO-accurate rate beyond lowest multiplicity is available. This can be combined with an inclusive NNLO cross section dσ (0+1+2) [INC] n (Φn) that is differential (only) in the lowest-multiplicity variables Φ n . It is useful to introduce the short notation 4 for the integrations in the all-order subtraction terms. Assuming that O i is a measurement of all degrees of freedom of the phase-space point Φ i (so that the second integral in eq. 4 can be omitted for lighter notation), the UN 2 LOPS matching formula then reads n+2 (Φ n+2 , t n+2 , t − ) . (5) Upon integration over both one-and two-parton states, the terms in (and multiplying) the boxes · , · and · cancel pairwise, leaving only the desired NNLO inclusive cross section. Integrating only over two-parton states, the terms in · cancel, so that for observables that require (n + 1)-parton ensembles, NLO precision is guaranteed. Note that for these (n + 1)-parton observables, the impact of (parton-shower) resummation at higher orders remains if t + t n+1 , as desired. The accuracy of the resummation procedure is not jeopardized by the matching, as the all-order factors are only shifted by finite, multiplicative, fixed-order factors (e.g. dσ (Φ n+1 )). Also, no assumption on the logarithmic accuracy of the parton shower needs to be made.
To avoid over-complicating the discussion at this point, eq. 5 assumes that all two-parton states admit an interpretation as ordered sequence of parton-shower transitions at scales t n+1 and t n+2 < t n+1 . This assumption heavily depends on the details of the parton-shower, as well as the hard-scattering process, and does not necessarily hold in all regions of phase space. The bulk of the cross section is captured by configurations with parton-shower interpretation, but additional "non-shower" configurations (called unordered in [6,33], and exceptional in [9]) need to be considered for an accurate prediction of sub-dominant phase-space regions. This complication (and related complications when matching three-parton states) are neglected in the main text, as well as the for the closure test in sec. V. A detailed discussion based on parton-shower ideas is given in Appendix B 1. In a full-fledged matching implementation, the treatment of non-shower configurations will also depend on the details of the necessary fixed-order (input) calculations.
The use of inclusive cross sections allows for particularly simple pairwise cancellation between configurations differing always by one parton. However, cross sections like dσ (0+1) [INC] n (Φn) or dσ (0+1+2) [INC] n (Φn) are often difficult to obtain due to the required integrations over the real-emission phase space. While inclusive cross sections are not available, it is easy to adjust the unitarized prescription to rely on "jet-vetoed" ("exclusive") cross sections instead 5 . If the veto scale is sufficiently small, these can be obtained by expanding resummed analytic calculations. A sufficiently small veto scale further allows to combine these jet-vetoed cross sections with higher-multiplicity cross sections to form precise results that differ from the exact fixed-order results only by very small power corrections. This strategy is e.g. employed by q ⊥ -or Njettiness-sliced fixed-order calculations. When using only jet-vetoed cross sections, and introducing another short-hand, the unitarized NNLO matching prescription becomes The main difference to the previous result is that the improved radiation patterns are only partially subtracted from the lower-multiplicity configurations, so that the full inclusive cross section is recovered upon integration. The couplingindependent kinematic factors 1 i+j i ≈ 1 ensure that the method to generate complementary i-parton contributions from i + j-parton contributions does not introduce Φ i -dependent biases that are not present in an inclusive fixedorder calculation. These factors need to be included to reconcile the requirement that the parton-shower accuracy is preserved (see Appendix A 3 for the definition of "shower accuracy"), and the requirement that the result is unbiased when expanded to fixed order. It is sufficient to think of 1 i+j i = 1 on first reading 6 . The cancellation mechanism between matched contributions in eq. 6 is now slightly more involved than in eq. 5, since the two-jet contributions · now also cancel between one-and zero-jet observables. Note that eq. 6 deviates slightly from the matching formula presented in [9], since the original proposition builds on an MC@NLO-matched 6 The reason for non-unity 1 i+j i factors is explained in Appendix A 4. While inclusive predictions using eq. 6 are prone to methodological bias (since exact complements need to be generated), exclusive predictions using eq. 5 can exhibit related issues, since the first term in the expansion of the unitarity subtraction should remove integrated hard jet configurations -which again has should avoid unwanted bias. Any unitary matching or merging method is prone to these issues. The similarity of using unitarity subtractions to the subtractions in the recent Projection-to-Born method [34] suggests that, if a selection between several Born-level states were required in that method, similar considerations may also become relevant. An assessment of 1 i+j i factors required for the present note is given in Appendix A 5.
one-jet radiation pattern, and handles the running QCD coupling differently. The matching formula 6 is closely related to the UNLOPS-PC prescription in [35]. Either method has benefits and down-sides, as discussed Appendix B.
Generalizations of the unitarized matching procedure to higher precision will also follow the logic suggested in Table  I. Section IV employs this reasoning to introduce a third-order prescription, reusing eq. 6 in the process.

III. CONSTRUCTING A TOY FIXED-ORDER IMPLEMENTATION FOR CLOSURE TESTING
To extend the precision of an event generator calculation to higher order, precise fixed-order calculations are required. This serves both as testing ground for new developments as well as physics deliverable. This note will only be concerned with the testing aspect of the method, since no differential N3LO calculations are available at present. A very controlled testing environment is in itself quite useful when developing a concrete matching implementation. To this end, a toy third-order calculation is constructed by combining and rescaling readily available tree-level results. This allows maximal control, and enables tests that would not be possible in a theoretically more rigorous calculation.
Gluonic corrections to e + e − → uū form the simplest laboratory for third-order matching ideas, due to the uncolored, non-composite initial state. To be applicable to hadron collider predictions, the matching prescription presented in section IV would require a detailed -and currently absent -understanding of the interplay between the factorization of parton distribution functions (PDFs) and renormalization within parton shower resummation 7 . The toy third-order calculation for e + e − → uū is constructed by 1. Producing tree-level event samples dσ For n > 0, very minimal regularization cuts S(Φ n ) > S c are applied on the projection of the (sum of) gluon four-momenta onto the four-momenta of the other partons, see Appendix A 5.
2. Constructing a sequence of toy fixed-order calculations. A toy exclusive NLO calculation for four-parton states is constructed from where the "TOY INC" contribution is a short-hand for term ∝ αs 2π [· · · ], dΦ R are the degrees of freedom for one additional parton, and y ij = 2p i p j /M 2 Z . The method to produce dΦ R integrals and the application of Q-constraints are explained in Appendix A 5. The term in · serves as proxy of logarithmic contributions of loop integrals, the · term as real contribution below the jet veto scale, and the rescaling αs 2π [· · · ] mimics the deformation of the spectra at NLO, after cancellation of logarithms between real and virtual corrections. The · term in eq. 8 only depends on five-parton states away from the phase space boundaries (since Q(Φ 2 ) > Q c ∧ Q(Φ 3 ) > Q c ) -the cancellation of real and virtual corrections close to the boundary is explicit, as is the miscancellation away from the boundary due to the veto condition. This also means that a dependence on the regularization criterion S(Φ n ) is effectively absent. Overall, this toy prediction has features typical for fixed-order jet-vetoed cross sections, e.g. a strongly negative contribution.
A toy NNLO calculation for three-parton states is constructed from 7 Also, current N3LO calculations for hadron collider processes employ NNLO PDF fits (as N3LO PDF fits are not yet available), leading to further PDF-related ambiguities, see e.g. [13]. Lepton colliders provide a more solid environment for developing N3LO+PS methods.
Again, the term in · approximates loop integrals, the · term implements jet-vetoed real-contributions, and a rescaling that mimics finite NNLO corrections is included. Again, the miscancellation due to the veto condition is again explicit. Finally, the toy N3LO calculation for two-parton states is assembled from dσ (0) 0 (Φ 0 ) and The argument for including the terms in colored boxes is identical to the previous cross sections. The peculiar functional dependences multiplying the coefficients a 0 , b 0 or c 0 would never arise from QCD corrections, but are chosen to allow for very clean closure testing.
3. Choosing the coefficients a X n , b X n and c X 0 in the toy calculation to produce exaggerated higher-order effects. The reproduction of these effects in the full third-order matched calculation will serve as affirmation of the correctness of the matching implementation. Explicit values are given in Table V. The resulting toy calculations exhibit the typical features of jet-vetoed cross sections, and thus serve as a realistic laboratory to test third-order matching ideas.

IV. THIRD-ORDER MATCHING
Having constructed a toy third-order calculation, one step towards concrete third-order event generators has been taken. This calculation should now be combined with parton shower evolution in a manner that guarantees that neither the accuracy of the fixed-order inputs nor of the parton shower are impaired. A unitarized matching ansatz offers a convenient route, since it does not put stringent (and not yet obtainable) restrictions on the parton shower implementation. The only relevant constraint is that on-shell intermediate states should exist at each step in the parton shower evolution, and that the parton-shower rate of any (pre-generated) phase-space point can be calculated numerically. The outlook in sec. VI will theorize on other potential third-order matching methods. A third-order matching method should fulfill the criteria 1. 3rd-order precision for inclusive zero-jet observables 2. 2nd-order precision for one-jet observables, combined with resummation when the jet becomes unresolved 3. 1st-order precision for two-jet observables, combined with resummed effects when either of the two jets becomes unresolved individually 4. 0th-order precision for three-jet observables, combined with resummed effects when any of the three jets become unresolved individually 5. parton-shower resummation of any observable sensitive to unresolved partons should not be impaired The starting point for a simple third-order matching is the availability of an NNLO+PS-matched +1-parton calculation F (Φn, t + , t − ) performed using eq. 6, and the availability of an N3LO jet-vetoed +0-jet cross n (Φ 0 ) for brevity). The recent emergence of N3LO differential cross sections in q ⊥ subtraction and Projection-to-Born methods suggests that this is not a particularly far-fetched requirement. The construction of the simple third-order matching method proceeds in the same spirit as the UN 2 LOPS matching prescription. Concretely, the method is constructed in the following steps: A. Regularize a one-jet UN 2 LOPS calculation with ∆ n (t + , t n+1 ) factors so that the hardest jet can become unresolved.
B. Remove unwanted NNLO terms from the regularized the one-jet spectrum.
D. Include and N3LO jet-vetoed zero-jet cross section.
The result is a valid third-order matching method, and will be referred to as Tomte (for third-order matched transition events). It should be noted that the prescription does not depend on the toy fixed-order calculation introduced before -the Tomte method would yield an N3LO+PS prediction given appropriate inputs. Furthermore, the scheme does not depend on the details of the parton shower, meaning that improved showers with reduced uncertainty will directly improve the Tomte method and decrease its uncertainty in regions dominated by partons becoming unresolved.

A. Constructing the improved radiation pattern
The baseline for producing a precise radiation pattern is eq. 6 after performing the shift n → n + 1. In this case, NNLO accuracy means that no unwanted terms are introduced at O(α 3 s ). Conversely, this means that reweighting most contributions in F (Φ n+1 , t + , t − ) with Sudakov-and running-coupling factors will only result in O(α 4 s ) shifts. There are two terms term that, when reweighted, threaten to introduce undesirable contributions. First first is the is the jet-vetoed 1-parton cross section dσ The potentially problematic term can be isolated by Any reweighting of the first line will only introduce terms of O(α 4 s ) or higher, while the second line requires careful consideration. Further unsafe terms are highlighted through · in eq. 6. An appropriate weighting strategy is outlined in sec. IV B. Performing the replacement in eq. 13, and introducing the placeholders u 1 , u 2 , u 3 and u 4 for the correct weights of dangerous terms, the one-jet contribution for Tomte can be written as where t n+1 was chosen as upper scale in ∆ n+1 to ensure non-overlapping resummation regions, as required to maintain the parton-shower accuracy.

B. Removing undesirable O(α 3 s ) terms from the radiation pattern
To obtain an appropriate radiation pattern, a correct weighting strategy for the terms proportional to has be be established. This weighting should not introduce unwanted O(α 3 s ) terms, while ensuring that the resulting terms are regularized as Q c → 0. Luckily, the O n+1 contribution to the UN 2 LOPS matching formula for n-parton processes (eq. 6) can act as a blueprint, since that calculation can be regarded as an approximation to the O n+1 contribution in Tomte. In the O n+1 contribution of UN 2 LOPS, the tree-level components are multiplied by the subtracted parton-shower factors , to avoid over-counting universal virtual corrections [9]. The same logic can be applied to the parton-shower reweighting for the O(α 2 s ) contributions to eq. 16, i.e.
This determines u 2 . The weighting of the O(α 1 s ) pieces of eq. 16 also need careful consideration. Starting from eq. 6, reweighting with subtracted parton-shower factors is appropriate. However, the subtracted parton-shower terms need to be expanded to second order to avoid spurious O(α 3 s ) terms. This suggests the shifts to ensure an appropriately weighted NLO cross section for exclusive (n + 1)-parton configurations. The last line removes undesirable O(α 2 s ) terms in the expansion of the all-order-reweighted O(α 1 s ) subtractions, and determines u 1 . Finally, the two-jet leading-order contribution (eq. 15) requires care, since it appears with observable dependence O n+2 , and, through unitarization and to complement the jet-vetoed cross section, with O n+1 dependence. This leads to several boundary conditions on suitable weights, as discussed in Appendix B. A suitable strategy is to use as contribution to O n+2 observables, fixing u 4 . The contribution to O n+1 cross sections contains the complement to a jet-vetoed cross section, and the unitarity subtraction of the O n+2 contribution. Thus, This then defines u 3 , and concludes the discussion of terms that cannot be trivially reweighted. After these considerations, the complete (n + 1)-parton spectrum for Tomte is finally given by Although a somewhat more complicated combination than the UN 2 LOPS prescription in eq. 6, the only two new factors that need to be calculated are w (2) n+1 (Φ n+1 ) and ∆ (2) n (t + , t n+1 ). Both factors are simple to generate in the absence of incoming hadrons. Were PDFs to enter the parton-shower evolution, then the generation would become more cumbersome 8 . However, it might be feasible to evaluate the pieces for color-singlet production processes at hadron colliders, since the terms are related to the evolution of zero-parton configurations and dynamical scale choices for the one-jet contributions alone. For example, a hybrid approach similar to the Minlo method [36] could be possible: Employ analytical Sudakov factors and their expansion for the factors related to the evolution of zero-jet configurations (∆ n ), and use numerical parton-shower results everywhere else. It should be noted that eq. 21 is an NNLO+PS matching formula for processes with an "additional" parton. The benefits of eq. 21 over, e.g., UN 2 LOPS matching are that the additional parton is allowed to become soft, or collinear to another parton. Furthermore, eq. 21 offers a clear, resummation-based, strategy for setting the argument of the running coupling for Φ n+1 states.

C. Unitarizing
After constructing an appropriate radiation pattern, unitarization is necessary to ensure that the target inclusive N3LO zero-jet is retained after matching. Unitarization further embeds the effects of parton-shower resummation at higher orders into zero-jet exclusive predictions. To this end, all observables in eq. 21 are replaced with O n , and the whole contribution is subtracted from the previous result. Using O n in eq. 21 triggers, as desired, the cancellation between the terms highlighted in · , · , · and · . After this, the prototype subtraction for unitarization reads 9 D. Completing the N3LO cross section to obtain the matching formula The last step in construction the Tomte matching method is to complement the matching formula with an N3LO exclusive cross section dσ (0+1+2+3) [EXC] n (Φn), and to ensure that the complementary parts (with real-emission configurations above the veto scale Q c ) are correctly included. This can be achieved by shifting the unitarization subtraction to not remove the necessary terms: where the necessary factors to complement the exclusive cross sections are highlighted in red. It should be stressed that these factors are introduced because the combination with a jet-vetoed N3LO cross section is foreseen. The combination with an inclusive N3LO calculation can be accommodated by simply ignoring the highlighted factors.
n (Φ 0 ) with eq. 23 and eq. 21 allows to construct the Tomte matching formula. As before, pairwise canceling terms will indicated with identical (hyperlinked) boxes. This acts as visual help to allow the reader to confirm that the criteria listed in Table III are indeed fulfilled. The final Tomte matching formula reads This is the main result of this note. As is always the case when matching fixed-order calculations to parton showers, this is not necessarily the only possible N3LO+PS matching scheme, provoking the question if the matching scheme uncertainty will outweigh renormalization scale uncertainties at N3LO. Scale variations in leading-order parton showers can be large [37], so that non-inclusive observables can be burdened with large uncertainties. By design, eq. 24 will immediately apply also when using higher-accuracy showers, since its derivation did not make assumptions on the shower accuracy. Thus, the scale uncertainties will immediately be reduced once higher-order showers become available. Matching scheme variations related to choices of the functional form of renormalization scales have been observed to introduce uncertainties in NLO merging [35]. Similarly, fixed-order NNLO calculations have been shown to exhibit a dependence on the functional form of renormalization scales that may be larger than conventional renormalization scale variations by constant factors [38]. The Tomte method employs very specific functional form of renormalization scales inherited from parton-shower resummation. Other scale setting mechanisms might be possible in N3LO+PS matching, potentially leading to non-negligible scheme dependence. At present, there is no sound way to assess this issue, without first building some intuition about matching at N3LO. Thus, to allow for toy studies, the matching formula 24 has been implemented in the Dire plugin to the Pythia event generator. The fixed-order methods presented in [17,18] appear particularly interesting for a future full-fledged N3LO+PS implementation, as do the fixed-order components of the resummed calculations [21,22]. Since [18] partly relies on an inclusive N3LO calculation, Appendix C documents a rearrangement of eq. 24 using an inclusive N3LO calculation for n-parton states.

V. NUMERICAL CLOSURE TEST
The toy fixed-order calculation constructed in sec. III allows for detailed tests of the implementation of the matching formula eq. 24. All the conditions for N3LO+PS accuracy listed in Table III can be tested in a controlled environment. These tests are presented in Figure 2 10 . The toy fixed-order calculation makes use of the renormalization scale µ =ŝ = M 2 Z and α s (µ) = 0.118, while all parton-shower factors use the running coupling α s (t) = α s 4 (p radiator · p emission )(p emission · p recoiler ) (p radiator · p emission ) + (p emission · p recoiler ) + (p radiator · p recoiler ) .
The prefactor in the argument of α s is chosen so that α s (t) → α s (µ) in fixed-order dominated phase-space regions.
That the tree-level result for observables requiring three or more additional partons is recovered is illustrated by fig.  2a. This shows that for well-separated jets, the fixed-order result is approached adequately. However, the agreement is largely accidental, since the pure fixed-order region is exceedingly small: the whole distribution receives significant Sudakov suppression (as illustrated by the dashed gray curve in fig. 2a), and the running-coupling effects present in the matched calculation have a significant impact also for well-separated configurations. When approaching phasespace regions with unresolved partons, the all-order factors included in the matched calculation produce a physically meaningful regularize the cross-section. In conclusion, the matched calculation behaves as expected, combining fixedorder accuracy with parton-shower resummation. The NLO accuracy of observables requiring two additional partons is assessed in fig. 2b. Again, the inclusive fixed-order result is approached as desired, although the pure fixed-order region (measured with the absence of Sudakov effects) is small, and running-coupling effects are large. As before, the matched calculation exhibits the desired Sudakov suppression when approaching phase-space regions containing unresolved partons. The two extremes are again consistently matched. The same conclusions can also be drawn for the NNLO accuracy of inclusive observables requiring at least one additional parton (as shown in fig. 2c), for which the fixed-order region is larger, and the agreement in the hard region is sound. Finally, the N3LO accuracy of inclusive observables relying only on the particles present at Born level is checked in fig. 2d. For such observables, the matched calculation should recover the toy N3LO calculation exactly, since no parton-shower factors rescale the inclusive lowestmultiplicity prediction. This is indeed the case, confirming that the implementation of the Tomte N3LO+PS matching method is consistent. Extensions to hadron collisions are conceivable, and would benefit from similar closure tests.

VI. SUMMARY AND OUTLOOK
This note introduces a straight-forward method to match N3LO calculations to parton shower resummation, thus allowing for the construction of N3LO-precise event generators. The construction of the method is based on the simple idea of unitary matching, and is an extension of the UN 2 LOPS method. The final formula eq. 24 appears a bit daunting at first glance, but should be easy to understand with the visual help of pairwise cancelling boxes. The Tomte matching formula had been implemented in the Pythia + Dire generator, assuming that the fixed-order results can be supplied by external calculations. A toy fixed-order calculation was constructed, and a closure test of the matching scheme was performed.
Although the necessity for N3LO+PS accurate event generators is not completely obvious, the development might be an interesting alternative to other avenues of improving event generators. In particular, given the limited application of NLO-merged predictions combining more than three NLO calculations in experimental analyses, it could be argued that there is little need for merging more than two NNLO calculations. If so, then it would seem that an N3LO+PS matching method would be adequate for the foreseeable future. However, if the goal of improved event generators is a decreased uncertainty throughout the spectrum, then improving the fixed-order precision only should not be considered an end on its own, but rather be accompanied by a better understanding of all-order perturbative effects. The derivation of the Tomte matching formula requires no reference to the actual logarithmic accuracy of the parton shower. Thus, eq. 24 is a valid method to combine N3LO calculations with high-precision analytic resummation.
This note offers a first proof of concept for N3LO+PS matching. Future research could apply the Tomte method to processes of interest, using realistic sliced fixed-order calculations. The current proposal could be improved with ease by using an MC@NLO-matched 2-jet component, or tailored to (in spirit similar) Projection-to-Born calculations. Is is also potentially possible to employ a fully differential NNLO+PS calculation as one-jet component. No fully-differential NNLO+PS methods for such processes exist to date, since they would require the parton shower to recover the complete singularity structure of QCD at O(α +2 s ). Nevertheless, such improvements could materialize within the next years. In any case, for simple processes, N3LO+PS matching is feasible with our current understanding of event generators.

VII. ACKNOWLEDGMENTS
This note is supported by funding from the Swedish Research Council, contract numbers 2016-05996 and 2020-04303.
The main text assumed that all factors in the various matching formulae (eqs. 4, 6 and 24) are known and can be generated numerically. The main aim of this appendix is to assess and confirm this assumption for all parton-showerrelated terms. The availability of appropriate fixed-order results is still left assumed.

Fixed-order cross sections
The formulae in the main text rely on differential fixed-order predictions, which are then combined with each other and the parton shower to obtained matched predictions. This section gives symbolic definitions of such calculations. Unless explicitly stated otherwise, n-parton configurations will be assumed to only contain well-separated phase-space regions, i.e.
The inclusive NLO cross section is given by where I n is the analytically integrated analogue of the real-emission subtraction S n+1 . The projection rate w FO n+1 (Φ n+1 , Φ n ) for replacing real-emission kinematics with underlying-Born-kinematics is discussed in Appendix A 4. The exclusive (jet-vetoed) next-to-leading order cross section can be obtained from its inclusive counter-part by restricting the phase-space for real-emission contributions: Note that the regularizing subtractions should be unaffected by the phase-space constraint. At NNLO, the inclusive cross section is given, symbolically, by where it is assumed that the subtractions S RV n+1 , S RR n+2 and the integrated subtractions I VV n add to zero in infrared-safe observables. The exclusive (jet-vetoed) counterpart can be obtained from the inclusive cross section by restricting the phase-space for real-emission contributions: The same relation between inclusive and exclusive cross sections persists at N3LO. In this case, the inclusive cross section is At NNLO, the inclusive cross section is given, symbolically, by where the subtractions S RVV n+1 , S RRV n+2 , S RRR n+3 and the integrated subtractions I VVV n add to zero for infrared-safe observables. The exclusive cross-section then reads

Histories
Parton-shower histories are crucial for calculating the matching terms. This note also employs parton-shower histories for the construction of cuts and integrated contributions in the toy fixed-order calculation described in sec. III. Some background on histories may thus be helpful. Showers generate multi-parton configurations through a sequence of parton branchings. A key realization of showercentered matching -and, particularly, merging -methods is to invert this picture when treating fixed-order calculations that contain multiple partons as input to the event generator. This leads to the concept of parton-shower histories, which are then utilized for many of the necessary tasks. A parton shower history is the set of all possible sequences of branchings that may have lead to the multi-parton state. Figure 3 shows the shower history for the production of a qgq final state: Two distinct sequences -starting from different states Φ Thus, the parton-shower prediction for an observable O 1 is an admixture of the rate of two paths. Assuming that Φ 1 can be parametrized by the dimensionless variables (τ, ζ, ϕ), and that the phase-space points fulfill τ < τ + , then the superscript " [1]" : factors pertaining to transitions where parton "1" emits ( · in fig. 3) superscript " [2]" : factors pertaining to transitions where parton "2" emits ( · in fig. 3) with the symbols defined in Table IV. Note that a spurious α s (µ)-factor was inserted for later convenience. Equation A9 highlights that parton-shower resummation produces a particular admixture of functional forms of the argument of α s , and of Sudakov factors with differing integration regions ([t(τ + ), t [1] (τ )] vs. [t(τ + ), t [2] (τ )]). The same admixture of all-order factors should be employed when calculating parton-shower factors for fixed-order matching. The partonshower accuracy of the prediction may otherwise be in danger.
• Calculate the factors highlighted by · and · , and include the factors for both paths according to these proportions. For example, this can be achieved by probabilistically picking one path according to these factors 11 .
• Calculate the ∆ 0 (t(τ + ), for the chosen path (i), and include this as rescaling of the cross section for the phase-space point Φ 1 .
This reasoning extends beyond this one-emission example [4,6]. Relevant histories for two-gluon and three-gluon states are shown in fig. 4 and fig. 5, respectively. The weight of a path in a history for more than one additional parton is given by the product of the weights for all the individual transitions in the path, potentially including shifted rates due to matching or merging fixed-order matrix elements, and including matrix-element correction factors and ordering constraints [33]. If the history H(Φ n ) is the collection of all paths p 0 , p 1 , . . . , p m from any Φ 0 to Φ n , then the mixing weight for path p a becomes where the variables z i−1 . All parton-shower factors necessary for matching n-parton configurations are calculated by constructing the full parton-shower history of Φ n and employing the method above. For example, the calculation of the · terms in eq. 24 require sequences , while the calculation of the · terms assumes knowledge of ( The sequences for products of all-order factors are thus constructed in analogy to the CKKW-L merging scheme [4], as admixture of contributions from all paths in the history. Subtracted factors (e.g. · in eq. 24) also rely on the construction of the parton-shower history, since the all-order factors being subtracted are determined from the history. More details on parton-shower factors can be found in Appendix A 7.

Definition of parton-shower accuracy
The logarithmic accuracy of parton showers is notoriously difficult to define, since parton showers are tools to produce event samples, i.e. are not observable-specific. Thus, defining the accuracy by discussing individual observables or small sets of observables can be misleading. Furthermore, the parton shower is, ultimately, a numerical algorithm. This algorithm should remain self-consistent even in matching procedures. This note thus employs an operational definition of the term "parton-shower accuracy" that is of algorithmic nature, and stricter than typical log-counting arguments: The all-order factors of the parton-shower should be reproduced exactly, such that no measurement could distinguish the parton-shower and the matched prediction if all fixed-order cross section were calculated using the same approximations employed to derive parton-shower splitting kernels. When using exact fixed-order calculations, three levels of defining a criterion for "maintaining the parton-shower accuracy" may be considered:

Strict parton-shower accuracy criterion
All parton-shower all-order factors are reproduced identically, particularly including their admixture (cf. eqs. A11 and A11), and their (approximated fixed-order) prefactors. No other sources of higher-order contributions exist. This would also entail that sub-leading or finite fixed-order contributions should not multiply all-order factors, and could only enter as corrections at a fixed coupling power. This strict condition may be considered undesirable, since it e.g. ignores arguments about the treatment of hard virtual corrections [40]. Consequently, no matching or merging method in the literature fulfills this overly strict condition. This strict criterion is not used in this note.

Balanced parton-shower accuracy criterion
All parton-shower all-order factors are reproduced identically, particularly including their admixture (cf. eqs. A11 and A11). The fixed-order prefactors of all-order factors may differ from strict parton-shower result by sub-leading contributions. Thus, complete fixed-order calculations (including finite and power-suppressed contributions) may multiply parton-shower all-order factors. This is the norm for matching methods. This note will apply this criterion when assessing changes to the all-order behavior due to matching.

Weak parton-shower accuracy criterion
Parton-shower all-order factors are reproduced in the strongly ordered limit, where a single parton-shower path saturates the rate to produce the configuration. This condition allows subtleties about the admixture of all-order factors, and is thus closely related an accuracy definition by enumerating logarithms of the shower evolution variable. In this approximation, it may be possible to replace numerical parton-shower factors by analytic expressions, as e.g. advocated in [3]. This weak criterion is not invoked in this note.
This note will employ the balanced parton-shower accuracy criterion. It follows that the "symbolic parton-shower factor notation" of the main text is defined by Terms multiplied by an expansion of shower factors should be calculated by multiplying each term in the sum by its own expansion. In order to claim compliance with parton-shower accuracy, the value of such "parton-shower factors" for a fixed Φ i+1 should be numerically identical, irrespective of Φ i+1 having been sampled by the shower, or if Φ i+1 had been pretabulated (by a fixed-order calculation) prior to showering, and the "shower factor" had been calculated post facto.

Generation of real-emission integrals and bias correction factors
Understanding the mixed weighting of the "radiation pattern" |M 1 (Φ 1 )| 2 in eq. A11 also leads to a definite method how to generate the (all-order) subtractions necessary for unitary matching. Since the shower produces an admixture of factors for fixed emission states Φ 1 , and the mixture is applied when matching an externally generated radiation pattern, the mixing should also be employed when projecting Φ 1 onto lower-multiplicity phase space points Φ [i] 0 . This ensures eq. 4 (as an extension of eq. 2) are accurately reproduced. The integration over the spectrum that is necessary to produce the (and , ) integrals in subtractions for unitarization, and to produce the dΦ R integrations for the toy fixed-order calculation are thus produced numerically. To give a concrete example, the (n + 1)-parton dependence of the third-to-last line of eq. 6 (highlighted by · ) can be obtained by a) tabulating two-gluon phase space points and constructing the history in fig. 4, b) calculating the probability w q n+2 of each path q, and picking a path according to its probability, c) multiplying the all-order factors for the chosen path p, and performing the replacement Φ n+2 → Φ

[p]
n+1 , and subtracting the result. This method will lead to the desired result Note that the requirement to recover the correct admixture relies on the weights w p n+2 , which depend on the details of both the (n + 1)-and the (n + 2)-parton configurations. That the probability to assign an (n + 1)-parton states depends on the (n+1)-parton configuration is necessary from the shower perspective. However, this is not appropriate for generating the complements to exclusive cross sections, since it can introduce artificial deformations of inclusive spectra. Take · in eq. 6 as an example. Since the sum of all w p n+2 is unity, any method to assign (n + 1)-parton states is equivalent if the degrees of freedom of the(n + 1)th parton are integrated over. This is however not the case for differential distributions. In fact, using eq. A11 to assign underlying states will favor (n + 1)-parton states with small inter-parton separation, due to collinear enhancements. Thus, the transverse momentum spectrum of the (n + 1)th parton will be deformed such that low transverse momenta become more likely than high values. Similar deformations were observed in [6], and erroneously attributed to mismatched phase-space mappings, casting doubt on that implementation. Furthermore, in the integration of real corrections to form B-cross sections in the Powheg method, similar artifacts may naively arise when using Catani-Seymour subtraction, and if the Born matrix element values vary significantly across phase space, and would have to be regularized [2].
The introduction of correction factors 1 i+j i allows to overcome such issues. These factors should ensure that a) inclusive fixed-order i-parton cross sections do not accumulate undue biases, and that b) a potential admixture of all-order factors for transitions from 0 to i additional partons is identical for all contributions to the inclusive cross section. For example, the factor 1 n+3 n+2 in the · terms in eq. 24 arranges that the integration matches the real-emission integral in the fixed-order prefactors of · terms. At the same time, it establishes that the method to produce all-order factors in · and · is equivalent.
It will be assumed that all parton-shower factors multiplying a specific fixed-order contribution will be calculated simultaneously. This calculation will proceed by constructing the history of phase-space points entering the fixedorder term. To obtain the correct admixture of shower factors, a path is chosen probabilistically according to eq. A11. If the path is also used to replace phase-space points with underlying configurations (e.g. for the sake of creating a subtraction), then an undesirable bias (relative to the fixed-order inclusive calculation) has been introduced. For Tomte matching, the bias correction factors 1 n+1 n , 1 n+2 n , 1 n+3 n 1 n+2 n+1 1 n+3 n+1 and 1 n+3 n+2 have to be defined. The most straight-forward bias-correction factors is 1 n+1 n , since its coefficients in eq. 24 does not include partonshower factors. Introducing the symbol w FO n+m for the rate at which a real-emission phase-space point Φ n+m would have been replaced by an underlying phase-space point Φ [p] n in an inclusive fixed-order calculation, this bias correction is This is obtained by reading off from the expected w FO weighting of the complementary low-separation real-emission n ), and dividing by the rate at which the shower method would have suggested the replacement Φ n+m → Φ which are again obtained by dividing the expected w FO weighting by the rate at which the shower method would have replaced the state. All other bias-correction factors (1 n+2 n+1 , 1 n+3 n+2 , 1 n+3 n+1 ) should also guarantee the correctly weighted mixture of partonshower factors for the underlying states. This is luckily already guaranteed by analogous definitions (A18) 12 For example, both paths Fig. 3 may lead to the leftmost underlying state in Fig. 4.
This can be seen by combining the rate with which the parton-shower would have "mixed in" the chosen path and bias-correction factor, for example leading to , by virtue of eq. A11, thus leading to the desired result. The same conclusion applies to 1 n+3 n+2 and 1 n+3 n+1 . While the factors 1 i+j i guarantee the correct fixed-order behavior, they unfortunately complicate the assessment of changes to the parton-shower accuracy for exclusive (i > n)-parton observables. To determine their impact, it is useful to note that terms related to 1 i+j i typically enter exclusive predictions in the form i+j (Φ i+j ) removes hard real-emission terms, while the remaining higher orders resum the effect of the jet veto. Thus, this term does not impair the accuracy. The (1 i+j i − 1) factor does not contain any large logarithms. It measures the difference between the fixed-order bias and the admixture of shower paths, and would vanish for an ideal fixed-order calculation that employs the method of Appendix A 2 to produce the real-emission integrals in the inclusive cross sections, or if there is a single dominant path from a phase-space point Φ [p] i to Φ i+j . Thus, log-counting suggests that the term can be omitted when discussing the all-order accuracy of exclusive observables. The term however changes the "accuracy of the parton shower" in the strictest sense, i.e. when requiring that the admixture of all-order factors is completely equivalent to the shower. The · term basically suggests that in a matched calculation, the definition of a fixed-order exclusive cross section should ideally be determined by the parton-shower method. Thus, the change in all-order accuracy is at the same level as attaching parton-shower resummation to finite parts of fixed-order cross sectionsand may thus be deemed acceptable. A detailed study of the impact of different Φ i -selection procedures could be an interesting addition to a future publication.
The functional form of the w FO i+j factor may differ between fixed-order methods that produce inclusive cross sections. The toy calculation described in Appendix A 5 will employ a non-trivial w FO i+j to allow the validation of the correctness of the 1 i+j i implementations.

Details on generating the toy fixed-order calculation
The toy third-order calculation used in this note in constructed from tree-level event samples generated with Mad-graph5 aMc@nlo. It should again be stressed that this is a toy calculation that only serves to assess the implementation of the Tomte formula eq. 24. Since tree-level events are used, the generation of samples containing additional  gluons requires regularization cuts S(Φ n ) > S c . Thus, minimal cuts are applied on the projection of the (sum of) gluon four-momenta onto the four-momenta of the other partons. The notation S(Φ n ) > S c implies very inclusive cuts, with S c = 0.1 GeV 2 . These cuts should be minimal enough to ensure that the majority of the radiation spectra is retained, and allow to assemble reasonable approximations for jet-vetoed fixed-order calculations.
These tree-level samples are then used to construct the toy jet-vetoed cross sections defined in eqs. 8, 10 and 12. Parton-shower histories are used to perform the dΦ R integration, and to apply the jet veto constraints Q(Φ n ) > Q c 13 : 1. A small value of Q c is chosen. The results in sec. V use Q c = 1 GeV 2 .
2. The parton-shower history for a pre-tabulated phase-space point Φ n is constructed (cf. Appendix A 2 and figs.

3, 4 and 5)
3. Only paths that satisfy t n > Q c are retained, where t n is the parton-shower evolution variable assigned to the transition Φ Note that this criterion depends on the path.
4. If no paths exist that fulfill the constraint, the event is discarded.
5. If valid paths exist, and no further constraints should be applied, and the event should not be processed further, then include the phase-space point in the toy calculation.
6. If valid paths exist, and the further cut Q(Φ n−1 ) > Q c should be applied (to produce the · integrals in eqs. 8, 10 and 12), then retain only the paths that satisfy the constraint. If no paths exist that fulfill the constraint, the event is discarded. Otherwise, the dΦ R integral is performed and an appropriate phase-space point Φ n−1 included in the toy calculation.
The dΦ R integrals are obtained with the method described in Appendix A 4, but omitting bias correction factors.
To avoid undue dependence on the parton shower splitting kernels, the requested Φ n−1 points are chosen amongst {Φ [1] n−1 , Φ [2] n−1 , . . . } according to the simpler weights where the variables z [r] and t [r] are calculated from Φ n and knowledge of its production from Φ n−1 , and q sums over all possible replacements. This differs slightly from the prescription for producing unitarity subtractions, since the latter mixes underlying states according to the weight in eq. A11. The prescription A21 is chosen to be less sensitive to the details of lower-multiplicity phase-space points, to allow non-trivial validation, and to increase the similarity to other methods for generating toy calculations such as [41]. 13 Note that the final form of the toy calculations only depend on constraints Q(Φn) > Qc that remove configurations too close to the phase-space boundaries 14 Another strategy would be to ignore this constraint and retain the paths, resulting in all phase-space points passing the slicing cut.
The state would then be employed in matching, and reweighted with parton-shower factors. At this stage, the path selection weight in eq. A11 would favor paths with tn < Qc, so that no-emission factors ∆n(tstart, t(Φn)< Qc) would multiply Φn states. Such contributions would be absent in the parton shower if Qc is identified with the parton-shower cut-off (as is assumed in this note), since the evolution would have terminated before being able to produce an emission at t(Φn), i.e. ∆n(tstart, t(Φn)< Qc) ≡ 0 when applied to Φn states. Thus, retaining the paths will have no impact on the final, matched result for Φn states, nor on their unitarity subtraction. The path may still contribute to the complement for exclusive (n − m)-parton calculations (where 0 < m ≤ n), if t n−m > Qc holds. Inspecting the definitions of exclusive cross sections in section III reveals that the effect of including these contributions to the complements is to cancel contributions to the exclusive cross section derived from paths with t(Φn)< Qc). Thus, the final result is identical to systematically discarding such paths. Only the selection of paths with tn > Qc will lead to a non-vanishing shower weight, i.e. contribute to the final result. The validity of this argument has been verified with the toy implementation to sub-percent level, for several tens of observables.
The · contribution to eq. 8 offers a tangible example for the dΦ R integration. This contribution can be obtained by a) tabulating three-gluon phase space points, b) constructing the history in fig. 5, c) retaining only histories for which the Q-constraints are fulfilled (see discussion above), d) calculating the weight A21 of each valid Φ n−1 . Otherwise, it is challenging to match the Sudakov reweighting between reweighted (toy) fixed-order inputs and Tomte-produced complements. The generation of the latter proceeds simultaneously to the generation of other Φ n -dependent factors, of which some require selecting parton shower paths of ordered emission sequences to Φ n . Thus, to produce matching Sudakov factors for complements and unitarity subtractions, the former should sum over the same paths. To obtain an exact match of the all-order factors generated when reweighting (toy) fixed-order inputs and their complements, the reweighting of the former should consider the value of the clustering scale t [p] when selecting valid shower paths. For the closure test calculation used in this note, not invoking the constraint can lead to ≤ 5% differences when approaching phase-space regions with multiple collinear partons. Retaining an additional clustering scale information is uncomfortable, since it can blurs the line between the fixed-order and the shower programs. However, the documentation of this scale in an inclusive fixed-order calculation should be technically feasible, and allowed within the accuracy of the calculation. Note that auxiliary (veto or shower starting) scale information is routinely deemed acceptable in POWHEG and MC@NLO matching.

Assessing the impact of bias correction factors
The sample generation discussed at the end of the previous section also allows non-trivial studies of the spectrum deformation effects (and the correction via 1 3 2 factors) discussed at the end of Appendix A 2. The current section will focus on such deformations for one instructive example. For other lepton (hadron) collider processes, deformations of similar (larger) size are expected. Figure 6a highlights that, given qqggg states, the method to assign underlying qqgg configurations does indeed lead to a non-negligible bias. The differences are similar in size to the naive expectation for NLO corrections. The deformation of the spectra is somewhat unexpected. Using A21 results smaller y 34 jet separation than the democratic clustering approach, yet the naive expectation is that the probability-based clustering removes the least-separated partons first, potentially yielding harder four-parton states. The bias introduced by eq. A11 is not completely obvious either, since it produces larger y 34 separation, while the naive expectation is that the dependence of eq. A11 on (an approximation of) the four-parton matrix element would tilt the spectrum in favor of less separated partons. The clustering bias is non-trivial. These undesirable biases can be removed by including appropriate corrective weights, (A22) 1 3 2 (FO) = 1/number of possible replacements of a 5-parton state with a 4-parton state as shown in Figure 6b, were different strategies are mapped to a democratic selection baseline. Any corrective weight, to convert mapping strategies into each other, may be obtained once the parton-shower history has been constructed. The price to pay for this flexibility is, as the comparison between Figures 6a and 6b indicates, a slower statistical convergence (given an identical number of events) in the latter due to an additional source of event weights. This is acceptable, since the weights are necessary to reinstate formal correctness. In the development of the Pythia +Dire implementation of Tomte, the removal of biases can explicitly be verified, due to having full control over the toy fixed-order calculations. Overall, it would seem prudent to assess the bias (potentially) introduced by the unitarization strategy also for other unitarized matching and merging approaches.

Parton-shower factors and trial showers
Parton-shower histories are paramount in defining the boundary conditions under which parton-shower and matching factors can be calculated. This appendix will be concerned with the technical implementation of the matching factors appearing in eq. 24. For the sake of brevity, the main text uses a symbolic notation for parton-shower factors. This symbolic notation for products of parton-shower factors f n+m (Φ n+m ) is, as required by the balanced parton-shower accuracy criterion, defined by the mixture where the sums extend over all possible parton shower paths. A detailed example for one contribution in eq. 24 is This mixture can be constructed with the algorithm given in Appendix A 2. From here on, the discussion will assume that a particular path p a of transitions has been chosen, and the path index p a will be omitted. The sequence of reconstructed states will be Φ 0 → Φ 1 → · · · → Φ n , and the evolution scales at which transitions occur are t + , t 1 , . . . , t n (where t + is a parton-shower starting scale). The factors w (∞) i (Φ i ) implement the effect of the running coupling, and are thus straight-forward to calculate, where µ = M 2 Z is used. For second-order running, the expansion of this term reads The second-order expansion always appears in the combination −w For first-order running, the β 1 -dependent term can be omitted. All w i -related factors are simple analytical expressions. The calculation of parton-shower Sudakov factors using analytic results is typically tedious, due to complicated z-integration limits in eq. 3. Luckily, since the sequence of states and evolution scales is known, it is possible to use trial showering to generate Sudakov factors and their expansions. The factors ∆ m (t + , t − ) may be generated by: 1) Initialize ∆ trial = 0, and define the number of trial shower sampling points N .
2) Initialize the parton-shower on the state Φ m , and with a maximal evolution scale t + .
3) Generate a parton-shower sequence. If the first parton-shower emission occurs at an evolution scale t < t − , shift ∆ trial → ∆ trial + 1.
Trial showering allows access to the details of individual parton-shower emissions that produce ∆ m (t + , t − ), and thus also allows to generate terms in the expansion of Sudakov factors. The most prevalent term required for Tomte matching is where P (0) (t, z, φ) is the sum of the first-order expansion of all parton-shower splitting kernels that may lead to a transition from Φ m 15 . Note the evaluation of α s at the reference scale µ. This term can be generated by extracting the average number of emissions from the trial shower [6]: trial + 1 w where 1 w = αs(µ) αs(tem) . The ratio of couplings is necessary if the trial shower is performed with dynamic α s argument, and guarantees that the desired fixed-α s (µ)-result is extracted. 4) Discard the emission, reset the maximal evolution scale to t em and repeat step 3) until t + < t − . 5) Repeat steps 2) -4) for a total of N times. Finally, set ∆ (1) This method is employed heavily in NLO merging methods, which also feature first-order shower expansions. The Tomte matching method further requires the generation of the second-order expansion This term always appears in combination with two other subtractions, The first term is most easily generated jointly with A28, since this avoids separating the splitting kernel into two components. The calculation of the term highlighted by · is a straight-forward extension of generating the average number of emissions -only the replacement 1 w → αs(µ) αs(tem) αs(µ) 2π β0 2 ln (t n+1 /t em ) is necessary to produce the contribution. The term in · is the rate of an emission witht ∈ [t, t n+1 ] given a previous emission at evolution scale t with t ∈ [t + , t n+1 ]. This may also be extracted from generating trial emissions: 1) Define the number sampling points N , initialize DE trial = 0, and initialize an empty set e trial 2) Initialize the parton-shower on the state Φ m , and with a maximal evolution scale t + .
3) Generate a parton-shower emission. If the emission occurs at an evolution scale t em > t n+1 , insert the weight 1 w = αs(µ) αs(tem) at the end of the set, e trial → e trial ∩ 1 w .
4) Discard the emission, reset the maximal evolution scale to t em and repeat step 3) until t + < t n+1 . For later convenience, it is useful to collect the complete form of the O(α 2 s ) subtraction that is required for all-order weighted tree-level n + 1-parton cross section dσ

5) Shift
Although the phase-space generation of modern parton showers is based on exact phase-space factorization formulae, and thus should allow full coverage, there are various restrictions that prevent parton showers from reaching all phase space regions. Possible restrictions are due to the parton-shower starting scale t + , and due to the requirement that parton-shower emissions are ordered by decreasing values of the evolution variable. The former restriction is absent for e + e − → jets, since the starting scale is typically chosen to beŝ = (p e + + p e − ) 2 . The ordering constraint is ameliorated if the shower includes the emission of multi-parton clusters, since the partons within one such cluster need not adhere to specific ordering constraints. This is e.g. the case for the two-parton emission clusters of the NLO parton shower presented in [42,43]. Configurations that cannot be reached by the parton shower will be called non-shower configurations.  Configurations that cannot be reached by showering require special considerations when matching to the parton shower. For N3LO+PS matching, configurations with one, two or three additional partons need to be considered. Possible one-parton states are listed in Table VI. In states with evolution variable above the shower starting scale (t(Φ n+1 ) > t + ), the additional parton cannot become unresolved, so that no Sudakov resummation is necessary. Thus, the b 1 contributions enter as fixed-order corrections. However, a sensible scale choice for the argument of α s is mandatory. The factorization scale value also needs to be considered carefully, since it doubles as the starting scale of a subsequent shower. A detailed scale setting mechanism for generic multi-parton states that is independent of the shower accuracy and is suitable for matched calculations can be found in [33].
Beyond these consideration -which apply irrespective of the desired accuracy -it is important to determine how non-shower states contribute the fixed-order prediction. If e.g. non-shower configurations are included in inclusive fixed-order cross sections, then a double-counting of the contributions when matching several multiplicities can be avoided by simply including non-shower configurations in the unitarization procedure. The b 1 contributions in Table  VI would e.g. by projected onto a Φ n phase-space point, and included as fixed-order subtraction. Alternatively, if non-shower states are not (or should not) be included in inclusive fixed-order predictions, then the contributions can be included without explicit unitarization.
States with two additional partons allow for many more evolution scale orderings, see Table VII. Non-shower states with one and two partons need to be considered, as well as the emission of a "soft" two-parton cluster. The former again demand suitable choices of renormalization and factorization scales. The effective scales introduced in [33] can serve this purpose. The evolution variable of double-emission contributions in an NLO parton shower offers a theoretically appealing scale definition for configurations with t(Φ n+2 ) > t(Φ n+1 ), since it may be applied irrespective of the value of t + . In Table VII, it was assumed that showering off hard Φ n+1 corrections would be initiated at t(Φ n+1 ), and that any scale hierarchies above t + are not large enough to warrant resummation. The latter assumption could easily be relaxed.
The inclusion of non-shower configurations with two hard jets (d 2 , f 2 , g 2 ) in inclusive predictions follows similar arguments as outlined for one-parton states above. If the contributions are included, in integrated form, in inclusive cross sections, then the contributions should be included in the unitarization procedure. For the configurations f 2 and g 2 , a fixed-order subtraction would have to be generated by projected the configurations onto Φ n+1 phase-space points, and then subtracting. Similarly, an all-order subtraction of b 2 should be generated with Φ n+1 dependence. Note that as an all-order reweighted contribution, b 2 should enter on equal footing to the other two-parton tree-level contributions in eq. 24, i.e. the Sudakov factor t(Φ n+1 ) and t(Φ n+2 ) should be subtracted appropriately. The two-jet  configurations d 2 include Φ n+2 states are fully resolved, while t + t(Φ n+1 ) may produce unresolved (n + 1)-parton states. Guidance on their treatment can be obtained from the double-emission shower of [42,43], which allows t(Φ n+2 ) > t(Φ n+1 ). There, the impact of double-emission clusters on inclusive observables cancels between n-and (n + 2)-parton states, since both emissions become unresolved simultaneously. Similarly, configuration d 2 may be regarded as emission of a hard two-parton cluster. Consequently, the necessary fixed-order subtraction for d 2 should be generated with Φ n dependence. Other treatments of d 2 are conceivable, and might become relevant when NNLO+PS and N3LO+PS matching mature. If inclusive fixed-order corrections do not contain non-shower states, then non-shower two-jet contributions can be included without explicit unitarization. Configurations a 2 and e 2 are reachable by an NLO parton shower and thus do not require special attention.
States with three additional partons allow for yet more evolution scale orderings (Table VIII). Note that table VIII assumes that showers off hard Φ n+1 corrections would commence at t(Φ n+1 ), and showering from hard two-jet configurations was started at t(Φ n+2 ) or t eff (t(Φ n+1 ), t(Φ n+2 )). Furthermore, any scale hierarchies above t + are not considered large enough to warrant resummation.
Multiple ordering combinations reachable by showering (a 3 , d 3 , n 3 , q 3 , t 3 , w 3 ) are possible. Several of these would only be accessible in an NNLO parton shower that includes the emission of correlated triple-parton clusters. Although this is currently beyond reach, the extension of [42,43] suggests that the evolution variable of correlated triple-emission may be ∝ p a (p 1 + p 2 + p 3 )p b (p 1 + p 2 + p 3 )/(p a p b ), where p a and p b are the momenta of (high-energy) radiator and recoiler after the branching, and p 1,2,3 are the emission momenta. Such a definition would constitute an appealing renormalization, factorization, and shower starting scale for non-shower configurations with three hard jets. Alternatively, the method of [33] might be used irrespectively of the shower accuracy.
The treatment of non-shower states with one or two partons was discussed above. Similar configurations appear also at three-parton level (b 3 , c 3 , e 3 , f 3 , s 3 ), and will be weighted by all-order Sudakov factors. The latter will depend on the (effective) scales assigned to the fixed-order configurations. If inclusive fixed-order cross sections include these nonshower configurations in integrated form, then unitarization through all-order subtraction is necessary. This proceeds by subtracting the contributions after projection onto (n + 2)-parton states, with the notable exception of s 3 . That contribution contains potentially unresolved (n+2)-and (n+3)-parton configurations with t(Φ n+3 ) > t(Φ n+2 ). Again taking guidance from NLO parton showering, the impact of these configurations on inclusive observables should cancel between (n + 3)-and (n + 1)-parton configurations, so that the necessary all-order subtraction should be generated with Φ n+1 -dependence.
Finally, configurations containing three hard jets can be treated as fixed-order corrections. If such configurations were included (again in integrated form) in inclusive cross sections, then fixed-order unitarity subtractions are required to retain the correct inclusive results. Depending on the ordering, these subtractions enter with Φ n+2 -dependence (g 3 ,h 3 ,j 3 ,o 3 ,r 3 ), Φ n+1 -dependence (k 3 ,v 3 ), or Φ n -dependence (i 3 , l 3 ,m 3 ,p 3 ,u 3 ,x 3 ). The last treatment is inspired by a hypothetical NNLO shower that would include the emission of internally disordered triple-parton clusters.

Differences in matching at NNLO
Comments on eq. 6 This note used the NNLO+PS matching formula eq. 6 as a starting point to derive N3LO+PS matching. This equation

Ordering
Interpretation Treatment a3) t+ > t(Φn+1) > t(Φn+2) > t(Φn+3) : Three ordered potentially unresolved emissions Sudakov between t+ and t(Φn+1), t(Φn+1) and t(Φn+2), and t(Φn+2) and t(Φn+3) b3) t(Φn+1) > t+ > t(Φn+2) > t(Φn+3) : Hard jet + two ordered potentially unresolved emissions Fixed-order contribution for Φn+1, Sudakov between t(Φn+1) and t(Φn+2) and t(Φn+2) and t(Φn+3) c3) t(Φn+2) > t+ > t(Φn+1) > t(Φn+3) : Two hard jets + potentially unresolved emission Fixed-order contribution for Φn+2, Sudakov between t eff (t(Φn+1), t(Φn+2)) and t(Φn+3) d3) t+ > t(Φn+2) > t(Φn+1) > t(Φn+3) : Potentially unresolved double emission, followed by potentially unresolved emission Sudakov from t+ to t eff (t(Φn+1), t(Φn+2)), and t eff (t(Φn+1), t(Φn+2)) and t(Φn+3)  differs slightly from the UN 2 LOPS prescriptions in the literature. The first change is the more abundant inclusion of running-coupling factors, especially the running-coupling rescaling of virtual and real corrections if additional partons were already present at tree-level. The reason for this choice is encapsulated in eqs. 2 and 4. These define the relation between the emission pattern and the no-emission (Sudakov) factor: the latter is determined by the difference of the integral of the former and unity. Conversely, the integral of the emission pattern is determined by the difference of the Sudakov factor from unity. Theoretical arguments strongly favor a dynamic coupling evaluation in the Sudakov exponent [44]. The same is true if the Sudakov factor is produced from an emission spectrum via explicitly unitarization. Hence, the integrated emission spectrum should employ a dynamic coupling evaluation. This also applies to (N)NLO fixed-order results. Thus, it is arguably more consistent with eq. 2 to include running-coupling factors in the rescaling of virtual and real corrections -though the inclusion of running-coupling effects is beyond the formal accuracy of the method. The second change relative to [9] is that eq. 6 demands and exclusive n + 1-parton NLO cross section, which is complemented by tree-level n + 2-parton distributions, whereas [9] employed an MC@NLO-matched n + 1-parton calculation. The main reason for this differences is technical: The fixed-order cross sections in eq. 6 are not dependent on the parton shower, so that it is straight-forward to produce toy fixed-order calculations for the inputs to eq. 6. The sliced approach of eq. 6 admits a more dynamical scale-setting method for n + 2-parton contributions than would be possible in an MC@NLO calculation, potentially resulting in a more "physical" result especially in the presence of large hardness hierarchies in n + 2-parton states. On the other hand, using an MC@NLO-matched calculation avoids the slicing of n + 1-parton configurations into two disjoint "bins".
Comments on fixed-multiplicity "bins" in the matching formula The matching formulae used in this note (eqs. 4, 6, 24) are combinations of weighted n-parton, n + 1-parton and possibly n + 2-parton and n + 3-parton states. This is indicated by the "observable dependence" O m (m = n . . . n + 3). The highest parton multiplicity sample will be distributed over even higher multiplicities by the action of the parton shower. In all other cases, there is no "smearing" of the description of m-parton states to higher multiplicities m + 1. This "binned" approach is common for matching [6,9]. However, it has the disadvantage that virtual corrections for m-parton states naively do not migrate to real-emission observables, at odds with the arguments of [40]. The implementation in this note accepts this feature. A method to ameliorate the "binned" behavior was discussed in [10].

Accuracy of the TOMTE matching formula
The formula eq. 24 allows the combination of N3LO calculations with parton showering, and hence with event generation. It fulfills all the criteria set out in Table III. This appendix provides a more fine-grained discussion of the choices leading to the Tomte method, and of its accuracy. Any confirmation of the desired accuracy relies on accurately reproducing fixed-order results for inclusive cross sections, and producing a combination of fixed-order and resummed results for exclusive cross sections.
More details on matching the tree-level two-jet contribution The matching of this contribution relies on its interplay with several other contributions, and thus requires careful consideration of any terms beyond the formal accuracy that are introduced by parton-shower factors. Firstly, the term forms the Born contribution to the two-parton NLO cross section. This suggests that it should contain all-order prefactors matching the prefactors of the real-and virtual corrections, i.e. dynamic arguments of α s for both partons, and Sudakov factors to ensure a suitable description in the presence of scale hierarchies between the partons. On the other hand, the term should complement the the exclusive one-parton NLO cross section dσ (Φ n+1 ) after integration, and should thus contain all-order factors identical to the latter (i.e. dynamic arguments for one power of α s , and only one Sudakov factor). Finally, the unitarization of O n+2 terms should result in an appropriate Sudakov factor to make the O n+1 contribution exclusive. Thus, any all-order weighting of the tree-level two-jet contribution is seriously constrained, leading to eqs. 19 and 20. Some more discussion can also be found in the following points.

Three-additional-parton observables are LO+PS accurate
Only the last term in eq. 24, n+3 (Φ n+3 , t n+3 , t − ) , contributes to O n+3 observables. The unitarity of the parton shower guarantees that for inclusive observables, this reduces to For exclusive n + 3-parton observables (i.e. when vetoing n + 4-parton states with separation larger than Q c ), the action of the parton shower eq. 1 attaches an additional Sudakov factor ∆ n+3 (t n+3 , t c )), and thus reproduces the correct parton-shower resummation of the jet veto. Hence, three-parton states are described at LO+PS accuracy.
Two-additional-parton observables are NLO+PS accurate The description of two-parton states in Tomte is very similar to the UN 2 LOPS scheme. For inclusive n + 2-parton observables, the terms in · in eq. 24 cancel by construction, resulting in a O n+2 contribution After expanding all shower factors, this gives such that the NLO cross section is recovered. Note that the factor 1 n+3 n+2 ensures that the (n + 3)-parton contribution combines with the (n + 2)-parton contributions in the correct manner, to produce an unbiased -integration. In order to assess how eq. B5 influences the all-order description (for one parton becoming unresolved, or both partons becoming unresolved in a strongly ordered manner), it is useful to reorder the terms into a "parton-shower prediction" and a "remnant term" [9], where dσ (1)[REM] n+2 The first term in eq. B7 is the desired parton-shower result. The remnant term in B7 consists of an all-order factor multiplying an O(α 3 s ) correction (eq. B8). This correction contains all non-universal NLO corrections, since all universal parts (as defined by the shower approximation) are removed from the complete NLO corrections. Thus, the "remnant" does not impair the shower accuracy. Instead, it can be considered favorable, since it allows to view the cross section B7 as separated into a "hard" coefficient (in brackets), dressed with identical all-order factors encapsulating the evolution of the state.
For exclusive (n + 2)-parton observables, the terms in · in eq. 24 do no longer cancel since (n + 3)-parton states with separation larger than Q c are vetoed. The (n + 3)-parton contribution can be written as The expansion · term removes hard (n + 3)-parton events at lowest order. At all orders, it combines (by virtue of eq. 1) with the other contributions to form a Sudakov factor reproducing the parton-shower resummation of the jet veto. As argued in Appendix A 4 the 1 n+3 n+2 − 1 term does not impair the accuracy according to the balanced parton shower accuracy criterion, and has thus been dropped in eq. B10. In conclusion, two-parton states are described at NLO+PS accuracy.

One-additional-parton observables are NNLO+PS accurate
The description of n + 1-parton observables O n+1 contains some of the most interesting features of the Tomte method. In inclusive n + 1-parton observables, the construction 24 enforces a cancellation of the terms in · , · , · , as well as in · . After these simplifications, the O n+1 prediction is where the 1 n+2 n+1 and 1 n+3 n+1 correction factors enforce an unbiased -integrations. The inclusive NNLO cross section is thus reproduced correctly. The influence of the matching on the all-order prediction -if the parton becomes unresolved -is best discussed by rewriting eq. B11 in terms of a parton-shower prediction and remnant term, and with the second-order remnant term defined in eq. B14. This again permits interpreting the matched result as a hard production coefficient (in brackets in eq. B13), dressed with the effect of soft-and collinear radiation. In the remnant term B13, all universal corrections are again subtracted, such that only non-universal corrections remain. An identical term appears in the UN 2 LOPS prescription. Thus, this term does not threaten the parton-shower accuracy of the matched result. The second-order remnant is given by These contributions are all interconnected, such that the impact of the second-order remnant on the all-order result is more subtle than for the first-order remnant. However, using eqs. A28 and A31 as well as the identifications the second-order remnant may suggestively be written as dσ (2)[REM] n+1 The NNLO cross section dσ (2)[INC] n+1 (Φ n+2 ) contains all the other terms, which consequently act to remove universal higher orders from the NNLO result. More specifically, in the bracket in the first line, only universal β 0 -dependent terms remain. These, together with the prefactor, act to remove β 0 -dependent terms from the NNLO result. Similarly, the bracket in the second line does not contain universal NLO corrections to soft gluon couplings 16 , so that only Sudakovlike effects remain. Multiplied with the prefactor, this factor acts to remove further universal terms. The remaining terms again remove universal contributions to the NNLO cross section. Thus, in total, dσ (2)[REM] n+1 (Φ n+1 ) contains -from the viewpoint of the parton shower -only non-universal terms, meaning that including dσ (2)[REM] n+1 (Φ n+1 ) as prefactor to all-order terms does not impair the desired parton-shower accuracy.
Finally, for exclusive n+1-parton observables (i.e. sensitive to exactly one additional parton), several all-order terms in eq. 24 do no longer cancel since n + 2-parton states with separation larger than Q c are vetoed. These uncanceled terms mirror, by design, the inclusive n + 2-parton prediction exactly. The latter provide an NLO+PS accurate description when one of the two additional partons becomes much softer (or more collinear) than the other. Thus, the uncanceled all-order terms in exclusive n + 1-parton observables provide a parton-shower accurate resummation of the jet veto.
In conclusion, n + 1-parton observables are described with NNLO+PS accuracy. Incidentally, this also means that omitting the n-parton contribution to Tomte yields an improved NNLO+PS matching for n + 1 parton processes.
Zero-additional-parton observables are N3LO+PS accurate In the prediction of inclusive n-parton observables O n , the cancellation between higher-multiplicity predictions and all-order subtractions take maximal effect. After replacing all observables in eq. 24 with O n (since inclusive n-parton observables are by definition insensitive to additional partons), the matching formula reduces to Thus, the inclusive N3LO cross section is recovered without any higher-order corrections. For exclusive n-parton observables, the all-order subtraction to balance the n + 1-parton prediction observables remains uncanceled, and thus produces a resummation of the jet veto. The inclusive n + 1-parton spectrum recovers the parton-shower resummation when the additional parton becomes soft or collinear. By virtue of the mechanism in eq. 2, this means that the allorder subtractions correctly reproduce the jet veto resummation for exclusive n-parton observables. To summarize, n-parton observables are described with N3LO+PS accuracy.
⊗ ∆ n (t + , t n+1 )∆ n+1 (t n+1 , t n+2 )w This is trivially obtained from eq. 24 by setting 1 n+1 n = 1 n+2 n = 1 n+3 n = 0. Other combinations of inclusive with exclusive calculations may be obtained in a similar manner. This simple change does, however, not guarantee that exclusive fixed-order cross sections are recovered exactly, since the unitarity subtractions might introduce biases. This was the reason why the use of exclusive cross-sections was prioritized in the main text. It would be valuable to assess and address the bias when using inclusive calculations in the future.