Matching and event-shape NNDL accuracy in parton showers

To explore the interplay of NLO matching and next-to-leading logarithmic (NLL) parton showers, we consider the simplest case of $\gamma^*$ and Higgs-boson decays to $q\bar q$ and $gg$ respectively. Not only should shower NLL accuracy be retained across observables after matching, but for global event-shape observables and the two-jet rate, matching can augment the shower in such a way that it additionally achieves next-to-next-to-double-logarithmic (NNDL) accuracy, a first step on the route towards general NNLL. As a proof-of-concept exploration of this question, we consider direct application of multiplicative matrix-element corrections, as well as simple implementations of MC@NLO and POWHEG-style matching. We find that the first two straightforwardly bring NNDL accuracy, and that this can also be achieved with POWHEG, although particular care is needed in the handover between POWHEG and the shower. Our study involves both analytic and numerical components and we also touch on some phenomenological considerations.


Introduction
Next-to-leading order (NLO) accurate event generators have become the de facto tool for the simulation of particle collisions at the LHC.This is in large part due to the success of the two most widely-used NLO matching schemes, MC@NLO [1] and POWHEG [2,3], alongside the two respective computer programs MadGraph5 aMC@NLO [4] and the POWHEG-BOX [5,6], as well as the independent implementations within the Herwig [7,8] and SHERPA [9,10] event generators.These and other matching procedures (e.g.[11][12][13]) were all developed at a time when parton showers had only leading logarithmic accuracy, and one of the questions that was relevant was whether they preserved that leading-logarithmic (LL) accuracy.
In recent years, a number of developments towards next-to-leading logarithmic (NLL) accurate parton showers have taken place [14][15][16][17][18][19][20][21][22][23][24][25], and one can now investigate the interplay between fixed-order matching and higher logarithmic accuracy.To do so, and to help highlight some of the essential considerations that arise, we focus here on the simplest possible process, NLO matching for γ * /Z → q q and H → gg (in the heavy-top limit), using the PanScales family of final-state showers, which all fully satisfy the broad NLL criteria set out in Refs.[14,15].
Aside from the obvious requirement that matching should preserve the (NLL) logarithmic accuracy of a given shower, the main question that we ask here is whether it can augment the logarithmic accuracy, at least in some situations.To help understand how this might be the case, it is useful to recall the traditional resummation formula for two-jet event shapes and jet rates in two-body decays [26].Specifically, the probability Σ for an observable O to have a value below some threshold e L (with L < 0) is given by Σ(O < e L ) = 1 + α s 2π C 1 + . . .e α −1 s g 1 (αsL)+g 2 (αsL)+αsg 3 (αsL)+... , |L| ≫ 1 . (1.1) The g 1 function is responsible for LL terms (α n s L n+1 ), g 2 for NLL terms (α n s L n ) and both C 1 and g 3 for NNLL terms (α n s L n−1 ). 1 An alternative way of writing Σ is where the h 1 function is responsible for double-logarithmic (DL) terms (α n s L 2n ), h 2 for NDL terms (α n s L 2n−1 ), h 3 for NNDL terms (α n s L 2n−2 ), and so forth.As is well known in the literature on event shapes, NLL resummation automatically implies NDL accuracy, and further inclusion of the C 1 in Eq. (1.1) is sufficient to achieve NNDL accuracy [26].In analytical resummation, C 1 is typically obtained through matching with a NLO calculation.The natural question in a parton shower context is therefore whether matching NLL parton showers with NLO retains the NLL accuracy and additionally achieves NNDL accuracy for two-jet event shapes and the two-jet rate. 2  In examining this question, we will consider three matching approaches: a straightforward multiplicative matching (similar to the matrix element corrections of Refs.[28,29], generalised in Vincia beyond the first order [30] and related also to the KrkNLO method [13]), which multiplies the event or splitting weight by the ratio of the true matrix element to the effective shower matrix element for the emission of an extra parton; the MC@NLO method, which adds in the difference between the true matrix element and the effective shower matrix element for the emission of an extra parton; and the POWHEG method, which takes full responsibility for the first emission, and then hands the event over to the shower for the remaining emissions.For the first two, retaining NLL and achieving cause the only kinematic region in which they act is the hard region. 3In contrast, because the POWHEG method always generates the hardest emission, which is often in the infrared (large-logarithm) region, the interplay with the shower logarithmic accuracy is more delicate.Note that understanding the matching/shower interplay of hardest-emission generator (HEG) methods, like POWHEG, is important because such methods also underpin the most widespread NNLO matching approaches, MiNNLO [32,33] and Geneva [34] (the other main NNLO matching procedure, UNNLOPS [35], does not fall into this class).
In Ref. [3] it was shown that the POWHEG Sudakov form factor achieves NLL accuracy for the hardest emission, and in most of this paper we shall take this NLL accuracy of the POWHEG Sudakov as a given.It should be clear that if POWHEG is then followed by a LL shower, that NLL accuracy will be squandered for practical observables.Less straightforward, however, is the question of what happens when POWHEG, with its NLL Sudakov, is followed by a NLL shower.
The subtleties that we have found in the POWHEG case are closely associated with discussions in the literature [2,27,36] about how to connect the effective shower starting scale with the scale of the POWHEG emission.In particular, it is standard to veto the shower branching, so as to provide the correct relation between POWHEG and the shower across the whole of the soft and/or collinear phase space, with no under-or double-counting.With the developments of the past few years on parton shower logarithmic accuracy [14][15][16][17][18][19][20][21][22][23][24][25], it is possible to approach these same questions with a range of new techniques, both analytical and numerical, that help analyse that logarithmic accuracy.In particular, concerns raised at the time of Ref. [36] about shower-induced recoil modifying the first (POWHEG) emission are closely connected with the origin of NLL failures discussed in Ref. [14] and, from a NNDL perspective, can be solved once one has NLL shower accuracy: the shower can still modify the first emission, but does so only when the shower emits in the immediate angular and transverse momentum vicinity of the first emission.This ensures that the effect of any modification is NNLL (which also implies that it is beyond NNDL for event shapes).In contrast, the questions of avoiding double-and/or under-counting will relate both with retaining NLL accuracy and the augmentation to NNDL.
In this paper, we will work with an e + e − -specific generalisation of POWHEG, which employs an ordering variable that coincides in the soft limit with the generalised PanScales ordering variable, v ∼ k t θ βps , parametrised by β ps .We will refer to it as POWHEG β . 4his generalisation enables us to avoid the considerable complications associated with the need for truncated showers [2] and the related subtleties of interplay with NLL accuracy.
However it leaves open interesting and important questions about the impact of mismatches in the hard-collinear region.As we shall see, these mismatches can arise not only from differences in kinematic maps, but also from the way in which showers partition the g → gg (and potentially g → q q) splitting functions, which effectively breaks their symmetry.Both aspects are relevant for NLL/NNDL accuracy.
Note that the ingredients that we discuss here, while sufficient for event-shape NNDL accuracy, are not enough to obtain NNDL accuracy for more general observables, e.g. for the recently studied sub-jet multiplicities [37,38].
The paper is structured as follows: In section 2 we briefly recall the relevant features of the multiplicative, MC@NLO and POWHEG methods.In section 3 we examine how kinematic mismatches and splitting function de-symmetrisation can affect NNDL accuracy, with both qualitative explanations and explicit calculations.We then verify numerically in section 4 that PanScales showers that are suitably matched to NLO do reach NNDL accuracy for a large set of global event shapes.In section 5 we briefly investigate the impact of the matching schemes in a phenomenological context.We conclude in section 6. Appendix A collects the results for C 1 coefficients that we use in Eq. (1.1), Appendix B discusses our treatment of spin correlations in the hard region, including fixed-order validation plots, and Appendix C gives further technical details of our matching procedures.

Brief overview of standard matching methods
For an in-depth review of NLO matching procedures, the reader may wish to consult Ref. [39], as well as the original papers.Here we give just a brief overview of the main matching approaches, with notation adapted from Ref. [31] so as to be more explicit about the starting scale for the showering following the matched emission.For simplicity, we will ignore the shower cutoff in all of our expressions.

Multiplicative matching
We start with the simplest kind of matching, multiplicative matching.We write the differential cross section dσ in such a way as to make explicit the structure associated with the first (i.e.hardest) emission, schematically Here Φ is the full Born plus one parton phase space, Φ B is the underlying Born phase space, and we define the radiation phase space, Φ rad , through dΦ = dΦ B dΦ rad .In the shower, for a given Φ B and Φ rad , there is an associated value of a shower ordering variable v ps Φ .The Born squared matrix element is given by B 0 (Φ B ), the true Born plus one-parton matrix element is R(Φ).The parton-shower approximation to it, R ps (Φ), is expected to coincide with R(Φ) in the soft and/or collinear limits.The parton-shower Sudakov form factor S ps (v, Φ B ) is given by the NLO normalisation factor can be written 5 where V (Φ B ) are the virtual corrections and the factor I ps (v ps Φ , Φ) represents the iterations of the parton shower branching from a value of the shower ordering variable v ps Φ onwards.Finally, we use the notation ⊗ inside the square brackets to indicate that the first emission is accepted with probability R(Φ)/R ps (Φ), and if it is not accepted then the attempt to create the first emission continues to lower v ps Φ .This also effectively replaces R ps → R in the integrand of the Sudakov form factor.The main practical difficulty to be aware of in multiplicative matching is the requirement R(Φ) ≤ R ps (Φ), in order for the first emission acceptance probability to be bounded below 1. 6 From the point of view of logarithmic accuracy, an important feature of Eq. (2.1) is that when v ps Φ is much smaller than the shower starting scale, the matching brings no modifications, because R ps = R in that region.

MC@NLO matching
Next we consider the MC@NLO procedure, which can be written as where Bps (Φ B ) = B 0 (Φ B ) + V (Φ B ) + R ps (Φ)dΦ rad . (2.5) The interpretation of Eq. (2.4) is that one generates normal shower events (with a suitable NLO normalisation, Bps ) and supplements them with a set of hard events, with weights distributed according to [R(Φ) − R ps (Φ)].This last term vanishes in the infrared (and so has a finite integral) because R ps coincides with R in the soft and/or collinear regions. 7This is important for logarithmic accuracy, because as for multiplicative matching, it ensures that there is no modification to the showering in the infrared region.Note that there is freedom in MC@NLO as to the choice of shower starting scale in the additional hard events, and here we will use v max , the largest accessible value.Modifying this, e.g.taking v ps Φ instead of v max , only affects the α 2 s contribution in the hard region.Assuming that this 5 Given the simplicity of two-body decays (we always take the decaying object to be unpolarised), NLO accuracy will not require an explicit calculation of the B (or, later, Bps) function of this section, but can instead be obtained simply by imposing the correct overall normalisation of the event sample, e.g. for γ * /Z → q q a factor of 1 + 3 4 CF αs π relative to the Born cross section. 6This condition may not always be satisfied, but as long as there are no holes in the shower phase space generation, one can generally find simple numerical workaround solutions, e.g. as discussed in Appendix C.1.2for the PanGlobal shower. 7As long as the shower has the correct full-colour structure of divergences, which it does in our segment and NODS colours schemes [16] for the simple cases of γ * and Higgs decays.Note also that below the shower cutoff, we take [R(Φ) − Rps(Φ)] to be zero.
-5 -multiplies a double-logarithmic Sudakov associated with further emissions, it is expected to affect results at α n+2 s L 2n , which is beyond our target NNDL accuracy.Our actual implementation of the MC@NLO procedure will differ slightly from Eq. (2.4) in terms of how we normalise different contributions.The details are given in Appendix C.2, and the differences correspond to NNLO corrections that are anyway beyond the control of NLO matching procedures.

POWHEG matching
A simple version of the POWHEG method can be written where the Sudakov form factor S heg is defined in analogy with Eq. (2.2) and R heg ≡ R. In this simple version of POWHEG, note the use of v heg Φ in I ps (v heg Φ , Φ), i.e. the ordering variable is deduced from the phase space map used with POWHEG (e.g. the FKS [40] map in the POWHEG-BOX), and adopted directly as a shower starting scale for the remaining shower emissions.Throughout this work, when the POWHEG method is used in conjunction with a NLL shower, we will assume that the HEG Sudakov factor S heg (v heg Φ , Φ B ) is also evaluated with NLL accuracy, i.e. using two-loop running of α s and the CMW scheme [41].
In POWHEG usage with standard transverse-momentum ordered showers, v heg Φ in POWHEG and v ps Φ in the shower will coincide for a given phase space point when the emission is simultaneously soft and collinear.This is essential in order to achieve leading logarithmic accuracy.This is not, however, naturally the case when using β ps > 0 showers in conjunction with the POWHEG-BOX, because of the latter's choice of transverse momentum as an ordering variable (which corresponds to β ps = 0).Accordingly, we will adopt a suitably generalised phase space map and ordering variable that satisfies this property also when used with β ps > 0 showers, cf.Appendix C.3.We will call this POWHEG β .It will be implicit throughout this work that the β value in POWHEG β is always taken equal to β ps .However, even with POWHEG β , the two values of v may still differ, for example, when the emission is collinear and hard.This creates a mismatch in the infrared between the phase space covered by the POWHEG β hardest emission generation and the phase space subsequently covered by the shower.As we shall see in detail in section 3.1, the mismatch has the potential to create NNDL issues, because the logarithmic phase space region associated with the mismatch is of order 1, the radiation probability in that region comes with one factor of α s and it can then multiply some part of a double logarithmic Sudakov factor, giving an overall α n+1 s L 2n .The aim of Ref. [36] was to eliminate the mismatch.Given the POWHEG first-emission event, the recommended approach within Pythia 8 starts the shower at the maximum allowed scale, but then for each shower emission i, based on the kinematics of i, the code deduces the equivalent value of the ordering variable in the POWHEG map, v heg i , and Φ .We write this procedure as where the additional showering condition is indicated after the vertical bar (|) in the I ps shower iteration factor.This should ensure the absence of holes or double-counting in the infrared phase space (at least when emissions are all well separated).With an important proviso concerning the handling of gluon splitting, discussed below in section 3.2, we expect this to be sufficient to restore NNDL accuracy for event shapes, as long as the underlying shower is NLL accurate.
3 Matching and event-shape NNDL accuracy In this section we will explore how O (1) phase-space mismatches in the collinear and/or soft region affect NNDL accuracy for event shapes, in particular in unvetoed POWHEG style matching, cf.Eq. (2.6).We will assume that both the shower and the hardest emission generator (HEG) have an ordering variable that coincides, in the soft-collinear region, with the PanScales family of ordering variables, parametrised by β ps , and written in terms of the emission transverse momentum k t and η (which for a back-to-back dipole is essentially a rapidity), as well as s ĩȷ = 2 p i • p j , s ĩ = 2 p i •Q, with Q the total event momentum and the tildes used to represent pre-branching momenta.Recall that β ps = 0 corresponds to transverse-momentum ordering, and that PanGlobal showers are NLL accurate for 0 ≤ β ps < 1, while PanLocal showers are NLL accurate for 0 < β ps < 1.In Eq. (3.1), ρ is a quantity that is equal to 1 for a back-to-back dipole in the event centre-of-mass frame.
As well as parameterising the shower evolution variable in terms of β ps , it is useful to parameterise event-shape observables in terms of a β obs (cf.Ref. [42], where it was called b).Specifically, when the event has just a single soft-collinear emission, the recursively infrared and collinear safe observable O is given by where k t is again the transverse momentum of the emission and η is its rapidity with respect to the emitter.For example, for the thrust [43,44], β obs = 1, and for the Cambridge √ y 23 3-jet resolution scale [45], β obs = 0.As discussed in the introduction, the quantity we will calculate is the probability Σ(O < e L ) for the observable to be below some threshold (where L is negative).Our discussion will be in two parts.In section 3.1 we will consider the purely kinematic differences between the HEG and shower maps, which brings a logarithmic analysis to the arguments of Ref. [36].In section 3.2 we will examine additional important considerations Figure 1: Lund plane [46] representation of the phase space for soft and/or collinear emission, and illustration of the interplay between a hardest emission generator (HEG), a parton shower and a constraint on an observable (black line).Dimensionful quantities (k t , v) are taken normalised to the centre-of-mass energy Q.
that arise for gluon splitting, connected with the widespread procedure in parton showers of partitioning the symmetric g → gg and g → q q splitting functions into two asymmetric versions, one for each dipole (in the g → gg case, each containing just a single soft divergence).
3.1 Kinematic mismatch between HEG and shower maps

Qualitative discussion of double-counting
The case we shall concentrate on is that where β obs = 0 and β ps > 0. This is represented in Fig. 1.The threshold on the observable is represented by the horizontal black line.We imagine a shower whose ordering variable coincides with the HEG in the soft limit, but in the hard-collinear region, for a given value of the ordering variable, generates a configuration that is somewhat harder than the HEG (by a factor of order one).This is illustrated with the green shower contour bending upwards relative to the red HEG contour.Though schematic, for simple e + e − → q q events, this picture captures the essence of a β ps = 1 2 extension of POWHEG together with the PanLocal β ps = 1 2 shower.It is clear that to have O < e L , the HEG emission must be below the upper blue contour in Fig. 1, i.e. it must satisfy ln v < L. If the HEG emission occurs between the two blue lines (e.g. on the red line), at a rapidity that places it below the black line (i.e. the observable threshold), the event will contribute to Σ(O < e L ) only if the shower does not emit further above the black line.That probability is given by the initial HEG Sudakov multiplied by the part of the remaining shower Sudakov that is above the black line.The fundamental problem is that the shaded region between the HEG (red) and shower (green) curves will be accounted for twice, once in the HEG Sudakov and once in the shower Sudakov, resulting in an over-suppression of Σ(O < e L ).
More generally we expect the following properties to hold: for contour mismatches in the hard-collinear region we anticipate NNDL artefacts for β obs < β ps .We do not expect an NNDL artefact for β obs ≥ β ps .In particular, β obs > β ps would correspond to a version of Fig. 1 where the black observable-constraint contour is steeper than the two blue shower contours, so the mismatch region is always below the observable-constraint contour and so does not affect the observable.An analogous argument can be made concerning potential contour mismatches in the soft large-angle region, for which we anticipate NNDL artefacts when β obs > β ps .For the case of e + e − → 2 jets, in the context of the PanScales showers and our specific POWHEG β formulation, there is no mismatch in this region.However, the issue of soft large-angle mismatches becomes important to address when the Born configuration involves dipoles that are not back-to-back, for example e + e − → 3 jets or pp → Z + jet.Such configurations are beyond the scope of this work. 8elow, we will give explicit calculations in the context of a double logarithmic approximation.Before doing so it is useful, however, to consider why it is that the HEG/shower combination can achieve NLL accuracy.Let us assume that the HEG and shower contours align in the soft and/or collinear region and also that the shower is NLL accurate (as is the HEG Sudakov).A crucial point to keep in mind is that the only difference between running the HEG and then the shower, or instead just the shower, is an O (α s ) relative correction to the overall normalisation, as well as the coefficient of the O (α s ) probability of having an emission in the hard region.Both are NNLL effects.Aside from this, the shower continues after the HEG in exactly the same way as it would had the first emission come from the shower.It is for this reason that the result maintains NLL accuracy.
One comment is that the agreement of contours can be arranged either by a vetoing procedure as in Eq. (2.7), or by adapting the design of the HEG and/or the shower such that their contours match naturally.This latter approach, while requiring more tailoring of the HEG/shower combination, has the advantage that it eliminates complexities associated with the implementation of the vetoing procedure, complexities that are likely to add extra challenges as one extends showers to higher logarithmic accuracy.

Evaluation of the discrepancy for
To make the argument quantitative, we will essentially work in a double-logarithmic (DL) approximation for the event shape distribution, and then specifically evaluate the subleading logarithmic effect of the mismatch.We will need the form of the HEG Sudakov, which we write in pure DL form as where C = C F for a q q event and C A for a gg event, and the HEG uses the same β as the parton-shower β ps .As a first step, it is instructive to consider how we obtain the analogous pure DL form for Σ(O < e L ) in a scenario where there is no mismatch between HEG and shower contours.It is given by where the first term in Eq. (3.4a) is the probability for the HEG emission to be below the lower blue contour in Fig. 1, ln v heg < (1 + β ps )L, i.e. it corresponds to S heg (e (1+βps)L ).
The second term is the integral for the HEG emission to be between the two blue contours, with two additional conditions.A first condition is that the HEG emission should be on the part of the HEG (red) contour that is below the observable constraint (black line) ln k t < L, which corresponds to a requirement where ln v ps ≡ ℓ in Eq. (3.4a).The integral over allowed rapidities gives the first factor in the integrand of Eq. (3.4a).A second condition associated with this integrand is that subsequent shower emissions should also be below the observable constraint (i.e.there should be no emissions in either the grey-shaded triangle, or the pink-shaded triangle).
The combination of HEG Sudakov and effective shower Sudakov for this second condition gives the exponential factors in the second (integral) term of Eq. (3.4a).
As a next step, we consider the impact of the mismatch between the HEG and shower contours, up to NNDL accuracy.The result is almost identical to Eq. (3.4), with just an additional term inside the exponent on the first line, involving a quantity ∆ that represents the effective size of one mismatch region (shaded green in Fig. 1).It arises because of the extra shower phase space that needs to be vetoed when the contours don't match.For e + e − → q q it reads ᾱ∆ = where is the reduced q → qg splitting function and θ ps ik (v, ζ) is the ĩ → ik splitting angle in the parton shower for a given value of the ordering variable v and post-splitting quark momentum fraction ζ; θ heg ik (v, ζ) is the analogue for the HEG.For the HEG/shower combinations that we consider, the ratio of angles in Eq. (3.7) will always be independent of v.
A key feature of Eq. (3.6b) is that the correction is indeed NNDL and starts at order α 2 s L 2 .Furthermore it vanishes for β ps = 0 (i.e.β ps = β obs since we have taken β obs = 0 in -10 -our calculation).This is consistent with our expectation that the discrepancy is present only for β obs < β ps .
To concretely evaluate Eq. (3.7), we need to know the θ ik (v, ζ) functions in the collinear limit.For the PanLocal and PanGlobal showers [15] and for our extension of POWHEG as given in Appendix C.3, we have PanLocal: PanGlobal, POWHEG β : Note that the equations are identical except in the last factor, which differs because in PanLocal, the emitter acquires the transverse recoil, while in PanGlobal that transverse recoil is shared across through a boost of the whole event (which in the collinear limit leaves the ik angle unchanged).In the case where we use PanGlobal or POWHEG as HEG followed by the PanLocal shower, combining Eqs.(3.8) with Eq. (3.7) results in the following expression for ᾱ∆: For H → gg, a similar calculation (with values of w gg = w qg = 0 for the parameters governing the de-symmetrisation of the gluon splitting functions, see section 3.2) gives A final comment is that the NNDL-breaking effects discussed here can also be thought of as a violation of the PanScales conditions that there should not be long-range correlations between emissions at disparate rapidities.Specifically, without vetoing, the presence of a soft-collinear HEG emission is effectively changing the probability of a hard-collinear shower emission.As a result, the HEG/shower combination fails to reproduce the factorisation that is present in physical matrix elements for emissions at disparate angles.This also impacts the analytical resummation structure.To see how, we take Σ from Eq. (3.6b), and extract the highest power L at each order in α s , which gives The result clearly fails to satisfy the exponentiation property of Eq. (1.1), specifically the absence of terms α n s L m in ln Σ with m > n+1.This can be viewed as similar, qualitatively, to the spurious super-leading logarithms seen in Ref. [15] (supplemental material, section 3) for standard dipole showers with observables such as the thrust.Below, when we summarise matched shower results together with their logarithmic accuracy, we will use the notation NLL, to remind the reader that the formal NLL accuracy has been lost.One subtlety, however, is that the difference between Eqs. (3.6b) and (3.4b) is always of relative order α s .This has the consequence that in numerical NLL tests with α s → 0 for fixed α s L, this difference would mimic a NNLL term, i.e.NLL accuracy would appear to be preserved despite the presence of spurious super-leading logarithms.
There are, nevertheless, observables that see a larger relative effect.One example is the invariant mass or transverse momentum of the first SoftDrop splitting when using β SD = 0 [47,48].The special characteristic of this observable is that it is not a standard global event shape, and its resummation does not have double-logarithmic terms, i.e. it starts from g 2 in Eq. (1.1).In the fixed-coupling approximation that we have effectively used in this section, the SD cross section has the following single-logarithmic structure, where c is a constant that depends on SoftDrop's z cut parameter, which we take to be small.Using the same strategy as above, one can explore how Eq. (3.12) is modified in HEG/shower combinations with a hard-collinear mismatch.Keeping β ps = 0 for simplicity, one finds where the coefficient ∆ that parameterises the impact of the HEG/shower contour mismatch now depends on z cut .As with Eq. (3.11), this generates spurious α n s L 2n−2 terms.If we examine the derivative of Σ SD (as we will do below in our phenomenology plots), we observe that there is a region, L ∼ 1/ √ α s , where the second term is suppressed relative to the first only by √ α s .Thus in this region, the impact of the HEG/shower mismatch is parametrically larger than the relative O (α s ) correction seen in Eq. (3.6b).

Additional subtleties for gluon splitting
The purpose of this section is to discuss an issue that can arise even when we have a HEG/shower combination whose kinematic contours (for a fixed value of the evolution variable) are aligned not just in the soft-collinear region, but for any single-emission phasespace point that is soft and/or collinear.The issue is connected with the fact that the g → gg splitting function has two soft divergences, one for ζ → 0 and the other for ζ → 1.This is a consequence of the inherent symmetry ζ ↔ (1 − ζ), which stems from the fact that g → gg corresponds to splitting to two identical particles (hence also the 1/2! factor).However, dipole showers break this symmetry, through the concept of an emitting particle (the "emitter") and a -12 - once the HEG has reached a given v value (v Φ ) without emitting, and an illustration that as the shower continues there may still be phase-space points (such as that marked with a cross) where the Sudakov has only been partially accounted for.The implications are discussed in the text.radiated particle.In particular, in order to help reproduce the correct pattern of largeangle soft radiation, dipole showers de-symmetrise the splitting function so that there is a divergence only when the radiated gluon becomes soft.For example the PanScales showers use 1 2! P asym where the choice of the w gg parameter fixes arbitrariness that arises in partitioning the finite part of the splitting function.It is straightforward to verify that P asym gg The hard matrix element generated by the HEG can be de-symmetrised similarly.The POWHEG-BOX code follows the FKS procedure [40], which introduces so-called S-functions that are used to partition the soft and collinear singularities.The de-symmetrisation discussed above is handled by an additional multiplicative factor h(ζ), cf.Eqs.(2.76)-(2.77) of Ref. [3], with ζ for an ĩ → ik splitting defined as E i /(E i + E k ).One can implement the scheme of Eq. (3.16) by setting The reason that the de-symmetrisation matters is that in many cases the kinematic map is not symmetric under ζ ↔ (1 − ζ).This can be seen in Eqs.(3.8), where the only combination that is symmetric is the PanLocal map for β ps = 0 (this, however, is not NLL -13 -accurate).The asymmetry is illustrated for PanLocal β ps = 1 2 in Fig. 2a, which shows fixedv contours in the physical Lund plane.Specifically for emission of k from an ĩȷ Born event, i.e. ĩȷ → ijk, the plot represents ln k t = ln[min(E i , E k ) sin θ ik ] versus η = − ln tan θ ik /2 in the right-hand Lund plane, and analogously in the left-hand plane with respect to j.A given contour has two parts: the lower branch corresponds to ζ > 1  2 , and contains the soft divergence; the upper branch corresponds to ζ < 1  2 and is free of any soft divergence, so it contributes significantly only in the hard collinear region.
The critical observation is that any specific point X in the hard-collinear region of the Lund plane can be populated by two distinct values of v: first at some v 1 by the lower (ζ > 1  2 ) part of the v 1 contour, and then at a smaller value of v = v 2 < v 1 by the upper (ζ < 1  2 ) part of the v 2 contour.The relative fraction of radiation intensity from the two values of v at a given point in the Lund plane depends on how the splitting function has been de-symmetrised, for example on the value of w gg in Eq. (3.16).For a given shower or HEG, we will refer to f X,1 (w gg ) as the fraction coming from the v 1 contour and f X,2 (w gg ) = 1 − f X,1 (w gg ) as the fraction from the v 2 contour.
Fig. 2b illustrates how this is relevant in the combination of HEG and shower.Suppose that we have a HEG with some w gg = w heg and a shower with some w gg = w ps .We consider a situation where the HEG generated the first emission at some v Φ , and where that emission is deep into the soft-collinear region, so that the branching should have no direct kinematic impact on subsequent hard collinear radiation.Next, let us focus on the point labelled with a cross ("X").That point is associated with two values of the evolution variable, v 1 and v 2 (remember, the map from evolution variable to kinematic contour is identical between HEG and shower across the full soft and/or collinear phase space).The first value, v 1 is already excluded by virtue of the fact that the HEG's emission corresponded to a smaller value of v, i.e. v Φ < v 1 .In effect some fraction f X,1 (w heg ) of the total Sudakov for not emitting at X was the HEG's responsibility, and that fraction depends on the de-symmetrisation parameter w heg of the HEG.The second value of the evolution variable, v 2 , that is associated with kinematic point X is the shower's responsibility, since v 2 < v Φ .If the shower does not emit there, the fraction of the total Sudakov that the shower contributes is f X,2 (w ps ) ≡ 1 − f X,1 (w ps ).If w heg = w ps then those two fractions add up to one.Otherwise they may be larger or smaller than one, which is equivalent to having a partial double-or under-counting of the radiation intensity at X, or more generically in the hard-collinear vicinity of the v Φ contour.
The impact of the double counting will be similar to that of the kinematic-contour mismatch discussed in section 3.1.Indeed all that needs doing to evaluate its NNDL effect is to take Eq.(3.6b) and replace the expression for ᾱ∆ in Eq. (3.7) with ᾱ∆ = where θ ik (v, ζ) is the same function for both the HEG and the parton shower.To understand Eq. (3.18), keep in mind that 0 < ζ < 1 2 is the region corresponding to the upper part of the contour in Fig. 2. If P ps gg (ζ) > P heg gg (ζ) then the combination of HEG and shower is enhancing the Sudakov in that region, much like the double counting of Fig. 1.The impact of the double counting on the event shape is structurally similar to that of the mismatch of contours in section 3.1.It can in fact be thought of as a weighted mismatch of contour, moving the mismatched part of the weight from θ ik (v, 1 − ζ) to θ ik (v, ζ), from which one deduces Eq. (3.18).
Note that there is a similar consideration for g → q q branchings, where the splitting function is sometimes also de-symmetrised.For completeness, using the PanScales form for the de-symmetrised g → q q splitting, we have ) Note that we will usually choose w heg gg = w heg qg ≡ w heg and w ps gg = w ps qg ≡ w ps , which leads to an almost complete cancellation between the C A and n f T R terms for n f = 5 (the cancellation would be exact for n f = 6).Accordingly, when we come to test the above analysis numerically with the full shower below, we will use n f = 0 so as to avoid this cancellation.
As with the kinematic mismatch of section 3.1, the effect that we have just seen corresponds to a violation of the PanScales conditions that there should not be long-range correlations between emissions at disparate rapidities, i.e. the presence of a soft-collinear emission (from the HEG) modifies the probability for subsequent hard-collinear emission (from the shower).This, again, has implications for exponentiation.

Practical implementation of vetoing
For completeness, we give here the specific algorithm that we adopt to achieve the vetoing of parton shower emissions following on from a first HEG step.The combination in which we will use it is with a POWHEG or PanGlobal-like HEG, followed by a PanLocal shower (dipole or antenna).This combination has the property (cf.Eq. (3.8)) that we only have to deal with double-counting, never with holes in soft and/or collinear phase space.The algorithm comes in two parts.
One part keeps track of the indices of the particles that should be considered the descendants of the particles in the Born configuration Φ B .For our cases later of e + e − → q q and H → gg, we label the Born partons a and b.For any branching (HEG or shower) where the emitter is one of the partons currently labelled as Born, e.g.ĩ → ik with ĩ a Born particle, the more energetic of i and k inherits the Born label.Note that for a q → qg splitting we could alternatively have considered the descendent quark to always be the Born particle.We will discuss the relevance of the choice below.
The second part of the algorithm follows the various HEG and shower steps: 1. Allow the HEG to generate the first emission, resulting in a Born+1 configuration Φ, with a value v heg Φ of the HEG ordering variable.The HEG step updates the Born indices as per above.

Start the parton shower with v ps = v heg
Φ .This is allowed because our HEG/shower combinations can lead to double counting, but not holes in the soft and/or collinear phase space.
3. For all subsequent emissions k from a dipole ĩȷ, if the emitter ĩ carries a Born label, check the following veto condition: For antenna showers, where the emitter/spectator distinction is absent, the same check is performed for either end of the dipole, insofar as it has the Born label.
4. If the splitting is accepted, update the Born labels.
A few comments may be helpful.The first concerns the freedom in how we assign the Born label.For NNDL accuracy, the critical element is that when Born-labelled parton ĩ splits as ĩ → ik, then when k is a soft gluon, it should be i that acquires the Born label.When k is a hard gluon, subsequent vetoing only affects triple-collinear configurations (and their associated virtual corrections), i.e. configurations with three partons at commensurate angles and with commensurate energies.Those do not play a role at NNDL accuracy.A second comment concerns the shower starting scale, which we could equally well have taken to be v ps = v max , as written in Eq. (2.7), as long as emissions from non-Born partons are only allowed for v ps < v heg Φ .This should not affect NNDL accuracy, and we have verified that in a phenomenological context the impact is small, at the percent level.Were we to consider HEG/shower combinations that result in IR holes, there would be less freedom in the choice of shower starting scale and it might well be necessary to start with v ps = v max .Finally, it is useful to be aware that the kinematic variables that we adopt for the contour check in step 3b differ from the Lund contours represented in Fig. 2, specifically with the use of E k rather than min(E i , E k ).Insofar as the meaning of ζ > 1  2 is the same in the HEG and the shower, this helps avoid complications with the folding of contours seen with the min(E i , E k ) choice in Fig. 2. However it is still necessary, in the g → gg/q q case, to ensure that the same de-symmetrisation of the splitting functions is used in the HEG and shower steps.

Logarithmic tests
We now turn to logarithmic tests with event shapes.There are quite a few potential combinations of matching scheme and shower.The subset that we explore is listed in Table 1.

mult. MC@NLO POWHEG β HEG PanGlobal
PanLocal β ps = 0.5 (dip.) The main matching and shower combinations for which we will test NNDL accuracy, for both γ * /Z → q q and H → gg.A (v) next to a check mark indicates that we use the vetoing algorithm of section 3.3.We use the NODS colour scheme from Ref. [16].
The combinations that do not have a check mark would also have been straightforward to test, 9 and are left out simply to limit CPU usage and bookkeeping.
We have run the standard NLL tests as in Refs.[15,16] at < 1% target accuracy, to verify that the matched showers continue to reproduce full-colour NLL accuracy for global event shapes (and NDL for multiplicity at ≲ 2% target accuracy).These tests were successful.Given the large number of results, we refrain from showing them.Instead, here we concentrate on the NNDL accuracy tests.
For the NNDL predictions we take the standard NLL formulae as used in earlier work and supplement them with the C 1 coefficients as given in Appendix A, some taken from the literature, some extracted numerically for this work.Slightly adapting the notation of the introduction, we denote by Σ(α s , L) the probability for an event shape observable to have a value O that is less than e L for a given α s (Q) = α s where Q is the γ * → q q centre-of-mass energy, or the Higgs-boson mass for H → gg tests.For a matched shower to qualify as being accurate at the NNDL level, its prediction for the cumulative distribution of a given observable, Σ PS , must clearly satisfy the following criterion, where fixing ξ is equivalent to fixing α s L 2 , as is relevant for isolating different terms in the DL-type expansion of Eq. (1.2).
To perform the tests, we run the showers with up to six values of α s = 0.1 N 2 , with N ∈ {3, 4, 5, 6, 8, 12}.We will show results at ξ = α s L 2 = 1.296 for γ * → q q, and ξ = α s L 2 = 0.791 for H → gg, which correspond to values of the cumulative distributions 0.25 ≤ Σ ≤ 0.6, i.e. in the bulk of the distribution for all observables under scrutiny.In each run, we choose a shower cutoff such that showering continues substantially below the smallest value needed to accurately predict the observable at the given ξ value.The limit in Eq. (4.1) is extracted numerically by extrapolating a fit that is linear or quadratic in powers of √ α s .In γ * → q q events, we use a linear polynomial fit to the three smallest α s values.
We repeat the fit with the four smallest α s values and take the difference in the intercept of both fits as a systematic uncertainty for the extrapolation procedure.In H → gg events, we find that there is still a visible quadratic component for some observables at the α s and ξ values we are considering, and thus fit a quadratic polynomial to all α s values (all but the largest for the systematic uncertainty).We will consider a range of observables with different values of β obs , cf.Eq. (3.2).Given the discussion of section 3, which showed that the presence of potential NNDL artefacts depends in a non-trivial way on the choice of both β obs and β ps , it is important to ensure that for each type of matching, we have shower/event-shape combinations with β obs < β ps , β obs = β ps and β obs > β ps .The observables that we take are • The Cambridge √ y 23 resolution parameter [45], total and wide jet broadenings B T and B W [26], which all have β obs = 0 (i.e. the same LL and NDL structures), but differ at NLL and NNDL.
• Three sets of parameterised observables for β obs = 0, 1 2 , 1: the fractional moments FC 1−β obs of the energy-energy correlations [42], as well as the sum and maximum, i k T,i /Qe −β obs |η i | and max i k T,i /Qe −β obs |η i | respectively, among primary Lund declusterings i [14,52].The first two observables are equivalent to each other at NLL accuracy but differ at NNDL, while the latter differs also in its NLL terms.
We start by showing results for Eq.(4.1) for showers without matching, in order to gauge the size of the NNDL discrepancy.This is shown in Fig. 3. Points are coloured in green if the central value is consistent with zero within 2σ (and in red otherwise), where the 1σ-uncertainty band is given by the statistical uncertainty and a systematic fit uncertainty added linearly. 10With the exception of the PanGlobal shower in H → gg, all showers without matching are clearly inconsistent with the NNDL result and the discrepancy can be significant, notably for the PanLocal showers where it is of order 2−3.The PanGlobal H → gg results are marked in amber, because the agreement is fortuitous: while the shower's effective 3-jet matrix-element is different from the exact result, we found that a coincidental cancellation leads to a seemingly correct result at n f = 5. 11 For each shower, the discrepancy is independent of the observable, because the test is carried out for small values of the observable (since α s L 2 is fixed with α s → 0), while the mismatch between the effective shower matrix element and the true matrix element is limited to the hard region. 10Note that in contrast to earlier PanScales work where statistical and systematic uncertainties were added in quadrature, this is a looser criterion.This choice reflects the significantly larger number of tests being performed here, and the correspondingly higher chance that at least one "correct" shower-matching scheme combination is mislabelled as having failed.Additionally, as in earlier work, in situations where the tests yield a result that is expected to be consistent with a given logarithmic accuracy but differs by more than 2σ, we extend the runs (either with further statistics or additional αs values) so as to clarify whether there is a genuine failure or not. 11We have also performed runs with n f = 0, 9 and verified that there is a non-zero NNDL discrepancy, which coincides with expectations from a simple semi-analytic calculation.2 ) and PanGlobal (β ps = 0 and β ps = 1 2 ) showers, without 3-jet matching, for γ * → q q (top) and H → gg (bottom).In these and subsequent plots, points marked green (red) show agreement (disagreement) with the NNDL results at the 2σ-level.Amber points manifest coincidental agreement for n f = 5 as explained in the text.The plots of this figure have no green points.

C-parameter
Next, we examine results of the matching with the multiplicative scheme.Fig. 4 displays the results of NNDL tests for the PanLocal dipole (β ps = 1  2 ) and the PanGlobal (β ps = 0, 1  2 ) showers matched multiplicatively.These are all in agreement with the NNDL result for all observables.
In Fig. 5, we show results of the NNDL tests where the PanScales showers are matched with the MC@NLO scheme.Here as well, the matched showers correctly reproduce all observables resummed to NNDL accuracy.
or antenna) shower is used as the main shower, we veto emissions according to the algorithm presented in section 3.3.For all combinations presented in Fig. 6, we also align the choice of de-symmetrisation in the gluon splitting functions, w heg = w ps = 0, following the analysis of section 3.2.Results are all in agreement with NNDL.Finally, we showcase the NNDL discrepancy arising from a failure to take the considerations of section 3 into account when matching showers with the POWHEG scheme.In Fig. 7, we display results of the NNDL tests for two combinations of HEG/shower which require vetoing due to kinematic mismatch (PanGlobal + PanLocal antenna β ps = 1 2 , and POWHEG β + PanLocal β ps = 1 2 , see Table 1), but where we deliberately do not apply the veto algorithm of section 3.3.As anticipated, disabling the veto has a visible effect for observables with β obs < β ps .The expected discrepancy for β obs = 0, which can be obtained by inserting Eq. (3.9c) (for γ * → q q) or Eq.(3.10) (for H → gg) into Eq.(3.6b), is plotted as a blue dotted line.That value is in agreement with the observed discrepancy from the showers.
Similarly, we can investigate whether the numerical results confirm the expected NNDL * qq, s L 2 = 1.296 (MC@NLO matching) H gg, s L 2 = 0.791 (MC@NLO matching) 2 ), is used as a HEG, followed by the same shower for subsequent emissions.In order to see the discrepancy of section 3.2, we choose different values of the de-symmetrisation parameter w gg , see Eq. (3.16), for the first emission (w heg ) and for the rest of the showering (w ps ).As can be seen from Eqs. (3.20a) and (3.20b), the discrepancy in both cases is proportional to C A − n f T R .In order to avoid the large cancellation of this effect with n f = 5, in Fig. 8 we run with n f = 0, and we set w heg = 1 and w ps = 0.While we could have extracted the C 1 coefficients numerically for n f = 0 as well, in Fig. 8 we only show observables for which the analytic form of C 1 is available.Recall that we only expect a discrepancy for β obs < β ps .Though the NNDL discrepancy associated with mismatched w heg ̸ = w ps is numerically smaller than for the case of the kinematic mismatch (note the different scale on the x-axis of Fig. 8), we find that the results agree with our analytic predictions for them, both when there are discrepancies and when there are none.

Phenomenological considerations
In this last section, we briefly explore the interplay of matching and logarithmic accuracy with physical choices for the coupling and values of observables, as opposed to the asymptotic values used in the preceding sections.The intent at this stage is not to be exhaustive, nor to compare to data (for that we would still like to have finite quark mass effects and an interface to hadronisation), but rather to get some insight into how logarithmic-accuracy improvements affect practical distributions.
We will show parton-level results with a preliminary estimate of uncertainties, specifically taking an envelope of two sources of uncertainty: (1) renormalisation scale variation and (2) uncertainties associated with residual lack of control of shower matrix elements beyond the matched emission.
The scale uncertainties are calculated according to Ref. [20]'s adaptation of the prescription by Mrenna and Skands [53].Specifically, for showers that are NLL accurate, we * qq, s L 2 = 1.296 (HEG matching, no veto)  exp.

H gg,
s L 2 = 0.791 (HEG matching, w HEG 2 ), using the same shower as both a HEG and for the subsequent parton showering steps, but with different choices of the de-symmetrisation parameter in the gluon splitting function, w heg = 1, and w ps = 0 respectively.The expected value of the NNDL discrepancy for β obs = 0 from Eqs. (3.20a)-(3.20b) is shown in dotted blue.take the emission intensity to be proportional to -23 -Here z is the fraction of the emitter momentum carried away by the radiation, 12 while b 0 and K are the usual β-function and CMW [41] coefficients, and µ central r is the emission transverse momentum (κ ⊥ ) as defined in the shower.The factor of 1 − z ensures that NLO scale compensation is present for soft-collinear emissions, but turned off for hard emissions. 13Scale variation will be probed by taking x r = { 1 2 , 1, 2}.We will use this scale variation also in the matching, e.g.so that HEG-style matching has the correct scale compensation in the infrared.Throughout, we use the NODS colour scheme [16] (i.e.full-colour NLL for global event shapes), a two-loop coupling, with α s (m 2 Z ) = 0.118, 5 light flavours and an infrared cutoff implemented such that α s (µ 2 r ) is set to zero for µ r < x r × 0.5 GeV.
We will also show results with our PanScales implementation [15] of the Pythia 8 shower [54] (which we call PSPythia 8).Since it is a LL shower, we will not include the scale compensating terms in Eq. (5.1) for shower emissions (nor for the HEG emission), however we do include a two-loop running and the CMW constant term K.
To estimate the uncertainty associated with lack of control of shower matrix-elements in the hard region, we modify the default shower splitting probability according to where x hard = 1 reproduces the default splitting probability and we take as variations As with the scale variation, there is some arbitrariness in these choices and their detailed implementation, whose full investigation we leave to future work together with that of other potential sources of uncertainty.For unmatched shower results, we apply Eq. (5.2) to all splittings, while for matched shower results we apply it to all but the first emission. 14The results are shown without spin correlations, but we have verified for the multiplicative matching procedure that they have a numerically negligible impact on the event-shape type observables shown here (and no impact up to NLL/NNDL).Note that fixed-order tests of spin correlations with matching are given in Appendix B.
Fig. 9 shows parton-level results for the thrust and Cambridge ln y 23 at √ s = m Z .It also features a SoftDrop (SD) ln k t /Q distribution, with z cut = 0.25, β SD = 0.0 [47,48] and k t for a ĩ → ik splitting defined as min(E i , E k ) sin θ ik .The choice of z cut = 0.25 is larger than commonly used, and helps concentrate on the hard-collinear region.The SoftDrop procedure is applied to each of the two jets as obtained from a 2-jet Cambridge clustering.This observable is shown for √ s = 2 TeV insofar as it is intended to be illustrative of an LHC jet substructure observable.
12 Specifically, in POWHEG β we take z equal to ξ in Eq. (C.13).For the PanGlobal and PanLocal showers, we take it equal to a k and b k in Eq. ( 4) of Ref. [15], respectively, when generating the g(η) and g(−η) terms. 13One might argue that scale compensation should be completely turned off earlier, e.g. for z > 1 2 .We leave exploration of different possible schemes to future work. 14In the MC@NLO procedure, x hard variation should ideally be done at all stages of the shower, with a compensating x hard dependence in the (R − Rps) term of Eq. (2.4).We have not yet implemented this, and accordingly we do not perform the x hard variation for the MC@NLO runs.
-24 -  The top row of Fig. 9 shows results for our implementation of the Pythia 8 shower.Recall that since this shower is LL rather than NLL we do not include the scale compensation terms of Eq. (5.1) when varying the renormalisation scale (neither in the shower -25 -itself, nor in the POWHEG β stage). 15The remaining rows show the PanGlobal shower with β ps = 0 and the PanLocal (dipole) shower β ps = 1 2 .Without matching, there are large uncertainties in the 3-jet region, mostly dominated by the x hard variations (dotted lines). 16 Once matching is turned on, it is the scale uncertainties that dominate the uncertainty bands over the full range shown for the event shapes.Note that the scale compensation that is used in the NLL shower brings a visible reduction in uncertainty as compared to the case of LL showers.
Observe also that POWHEG β matching with the PanLocal shower without the kinematic veto, shown in the lower row (green band), induces a noticeable shift in the kinematic distributions.Curiously, this is true not just for the ln y 23 distribution (where there is a NNDL effect), but also for the thrust (where there is no NNDL effect, but still a kinematic double-counting).The clearest effect is seen for the SoftDrop observable at intermediate values of ln k t /Q.This is not unexpected: recall from the discussion at the end of section 3.1 that for moderately large values of the logarithm, L ∼ 1/ √ α s , the impact of the HEG/shower contour mismatch is expected to be of relative order √ α s rather than at most α s for standard global event shapes.In Fig. 9 it was possible to show only a limited number of matching/shower combinations.To help visualise characteristics of a wider range shower/matching combinations, we select a specific bin in each of two distributions, 0.14 < 1 − T < 0.15, and the SoftDrop −3.1 < ln k t /Q < 3.0, chosen so as to be in a transition region between the large-logarithm and hard 3-jet regimes.We then examine the ratio of a range of showers and matching schemes to the multiplicatively-matched PanGlobal β ps = 0 result.This is shown in Figs. 10 and 11, with one row for each shower/matching combination.
Features to note beyond those observed for Fig. 9 are: (1) The matched NLL showers have uncertainties that are broadly similar across matching and shower combinations (somewhat larger for PanLocal β ps = 1  2 ), and consistent with each other to within uncertainties.Note that the residual 10% variation between matched showers could have a significant impact on tuning, and ideally one should include 3-jet NLO matching, e.g. as done in Ref. [55].(2) The characteristics seen in Fig. 9 for POWHEG β matching with Pan-Local (dipole) β ps = 1 2 without the kinematic veto appear to be replicated also when using the PanGlobal shower as a HEG (HEG PG ) in combination with the PanLocal (antenna) shower without a veto, suggesting that they are genuinely an impact of the lack of veto rather than a coincidence.
Overall the results in this section highlight the importance of both matching and NLL   accuracy and of bringing them together consistently, lending support to the practical value of pursuing the programme to improve shower logarithmic accuracy.

Conclusions
The simple framework of two-body decays that we have studied here provides a clean and powerful laboratory to explore the interplay of parton-shower logarithmic accuracy and matching.In particular, we have seen that NLO matching can augment the accuracy of NLL showers, so that they additionally attain NNDL accuracy for global event shapes.This was relatively straightforward to achieve with multiplicative and MC@NLO matching methods, because they alter the shower behaviour or add events only in the hard -27 -   region.In contrast, with matching methods such as POWHEG that take responsibility for generating the hardest emission, an extra element is needed, which is to ensure that the hardest-emission generator and shower align in their generation of phase space in the full soft and/or collinear regions.Failing to account for this prevents the HEG/shower combination from attaining NNDL accuracy.Furthermore, it subtly compromises NLL accuracy, generating spurious super-leading logarithms, Eq. (3.11), that resum in such a way, Eq. (3.6b), as to vanish in standard numerical NLL accuracy global event-shape tests (but not necessarily for single logarithmic observables, such as SoftDrop with β SD = 0).In this paper we used the (standard) approach of vetoing shower steps in order to avoid double-counting phase space already generated with the HEG.However, thinking forward to possible approaches to achieving yet higher logarithmic accuracy, it is likely to be advantageous to consider designing HEG tools such that they have the freedom to mimic the lowest order soft/collinear phase-space generation of any given shower.A related and more subtle issue occurs when a given phase-space point can be reached from more than one value of the HEG or shower ordering variable.In our study, this issue arose in the context of de-symmetrisation of gluon splitting functions in the hard-collinear region, cf.section 3.2.However, we expect it to be relevant more generally also in processes with non-trivial soft large-angle structures, for example in the presence of three or more Born partonic legs.The critical observation is that in such situations, the HEG and the shower must have the same relative weights for each of the distinct values of the ordering variable that can lead to that phase-space point.
The numerical tests of section 4 provided extensive validation of our understanding of parton-shower NNDL event-shape accuracy.The results confirmed NNDL accuracy in the situations where we expected to achieve it, and furthermore reproduced the analytic expectations for discrepancies in situations with kinematic or de-symmetrisation mismatches.
In section 5, we took first steps towards the exploration of the phenomenological impact of logarithmically accurate showers, including preliminary uncertainty estimates.Perhaps the main conclusion to be drawn so far is that there is good consistency across different matching schemes and showers, to within the uncertainties, as well as significant reductions in uncertainties relative to LL showers.Furthermore, matching/shower combinations that do not achieve NNDL accuracy appear to give predictions in tension with those from NNDLaccurate showers.Clearly next important phenomenological steps include interfacing with hadronisation, the inclusion of heavy quark masses, comparisons to data and associated exploration of tuning NLL showers. . (B.9)

H → gq q
We denote the amplitude as M (λ g , λ q , λ q).The non-vanishing helicity configurations are M (λ, λ, −λ) = − S λ (p q , p g ) 2 S λ (p q , p q) , (B.10) .11)Note that this corresponds to the H → gq q amplitude as mediated by the heavy-top limit of the effective loop-induced Hgg operator, as opposed to gluon emission following a direct Hqq Yukawa interaction.

B.1 Validation
We validate the implementation of the three-body hard scattering amplitudes at second order by comparing to exact matrix elements of MG5 aMC@NLO [4].We consider configurations where the first matched emission with momentum p k occurs from a dipole with momenta pi and pj at a relatively large, fixed transverse momentum where p i and p j are the post-branching dipole constituents, while varying its rapidity Figure 13: Second-order validation of the implementation of spin correlations in the hard region for e + e − → γ * → q q.Following Eq. (B.14), the ratio a 2 /a 0 is computed for the shower with matching enabled (blue curve).Below, the ratio of the shower a 2 /a 0 to the exact matrix element a 2 /a 0 is shown.
A second collinear branching then occurs with a small, fixed opening angle θ 2 = 10 −5 and collinear momentum fraction z 2 = 0.4.When the 3-jet event is produced via radiation of a gluon, it is the collinear splitting of the radiated gluon that we examine; in H → gg events where the 3-jet event is obtained through a g → q q splitting, it is the collinear splitting of the remaining Born gluon that we examine.The difference between the azimuthal orientations of the branching planes ∆Ψ 12 is then sensitive to spin correlations.An illustration of this observable is shown in Fig. 12.In particular, the cross section has the form dσ d∆Ψ 12 dη 1 ∝ a 0 (η 1 ) + a 2 (η 1 ) cos (2∆Ψ 12 ) .

(B.14)
The values of a 0 and a 2 can be extracted through a Fourier transform.The results are shown in Fig. 13 for e + e − → γ * → q q and in Fig. 14 for H → gg.The ratio (a 2 /a 0 ) is shown for the PanGlobal shower with matching enabled.Furthermore, the ratio between the parton shower (a 2 /a 0 ) and the exact matrix element (a 2 /a 0 ) is shown, confirming that the matched shower reproduces the exact matrix element.

C Technical details of the matching implementation
The tests performed in section 4 require us to probe very small values of α s and very large values of the relevant logarithm.In the parton shower, this corresponds to very small values of the cutoff scale, which in turn requires precise control over the numerical accuracy of the parton shower.If we were to make use of the public implementations of POWHEG [5] and MC@NLO [4], we would lack this precision.The results of section 4 are thus obtained  13, but for all possible configurations of H → gg.Note that in the lower plots, it is the gluon of the underlying 3-jet system that undergoes collinear splitting, rather than the hard emission from the 2-jet system (which is a quark or antiquark, and so does not mediate azimuthal correlation).This is also the reason for the lack of symmetry in the matched result between positive and negative η 1 values.
with a simple dedicated implementation of POWHEG and MC@NLO for the processes at hand within the PanScales framework.
In the following, we provide some details concerning the implementation of multiplicative matching in the PanScales showers, and of the implementation of POWHEG and MC@NLO in the PanScales framework.

C.1 Multiplicative matching
The multiplicative matching of the PanScales showers is accomplished by replacing the usual splitting functions, which are only correct in the soft and collinear limits, by R(Φ)/B 0 (Φ B ) for the first emission only.In situations where the shower partitions the total emission weight into multiple dipoles or antennae, the same partitioning is applied to the weight R(Φ)/B 0 (Φ B ). Furthermore, a Jacobian factor must be included to account for the transformation from the radiative phase space to the shower variables.
For completeness, we give here the analytic expressions of the effective one-emission matrix elements for the PanLocal dipole and PanGlobal showers, as well as the Jacobian factor mentioned above.We also briefly discuss a workaround for the issue of undersampling in the PanGlobal shower, which appears in hard configurations.(C.2) The Dalitz variables (x i , x j ) are related to (a k , b k ) by Finally, the Jacobian for the Panlocal dipole shower (for a single dipole end) can be expressed as Here the relation between the Dalitz and the shower variables is given by and the Jacobian for PanGlobal can be shown to be The Jacobian has a divergence on the contour x i + x j = 1, or equivalently at the point ln v = 0.In other words, the inverse Jacobian vanishes on that same contour, and the shower has an emission probability that is exactly equal to zero on that line.We tackle this issue by replacing the shower's sampling of ln v values: with C > 0. This modification only impacts the hard region, ( On the problematic contour, the modified radiation intensity becomes lim and is thus non-zero for C > 0. In practice we take C = 1 2 .

C.2 MC@NLO
The implementation of MC@NLO amounts to the production of a mixture of two sets of events, one distributed according to the Born matrix element, the other according to the real correction term.Rather than implementing a separate phase space generator for the second set, a simpler method is possible when, following section C.1 a multiplicativelymatched version of the shower is already available.In that case, we can use the shower as the phase space generator.Specifically, we employ a first-order expansion of the shower, in which a radiating dipole is selected randomly with equal probability.Then, event weights equal to the difference between the matched and unmatched weights of the radiating dipoles are attached to the resulting event.
Our normalisation convention differs slightly from the standard MC@NLO approach, mainly for reasons of programmatic convenience.Specifically we use where B 0 , B and Bps are all to be understood as being functions of Φ B .It is straightforward to verify that this differs from Eq. (2.4) only starting from order α 2 s relative to the Born cross section, i.e.NNLO.

Figure 2 :
Figure 2: Schematic illustration of the issue associated with gluon asymmetrisation.(a) Contours on the Lund plane, in the PanLocal family of showers, highlighting the fact that a given physical point X in the Lund plane (highlighted with a red cross) can come from two different values of v.The shading of the green curves represents the variation in radiation intensity along the contour.(b) Density plot, at each point in the Lund plane, representingschematically the fraction of the emission intensity at that point that has been excluded once the HEG has reached a given v value (v Φ ) without emitting, and an illustration that as the shower continues there may still be phase-space points (such as that marked with a cross) where the Sudakov has only been partially accounted for.The implications are discussed in the text.

( a )
Compute k t = E k sin θ ik and η = − ln tan θ ik /2 (this is done using exact momenta).(b) Compare these coordinates with the corresponding contour of the HEG at ln v heg Φ , i.e. ln k heg t (v heg Φ , η) = ln v heg Φ + β ps |η| for POWHEG β and for PanGlobal as a HEG.(c) If the emission is above the contour, ln k t > ln k heg t (v heg Φ , η), veto the emission.

Figure 3 :
Figure 3: Results of the NNDL accuracy tests at fixed ξ = α s L 2 for the PanLocal dipole and antenna (with β ps = 12 ) and PanGlobal (β ps = 0 and β ps = 1 2 ) showers, without 3-jet matching, for γ * → q q (top) and H → gg (bottom).In these and subsequent plots, points marked green (red) show agreement (disagreement) with the NNDL results at the 2σ-level.Amber points manifest coincidental agreement for n f = 5 as explained in the text.The plots of this figure have no green points.

Figure 5 :
Figure 5: Results of the NNDL accuracy tests for the PanScales showers (see Fig.3) matched with the MC@NLO scheme.

Figure 6 :
Figure 6: Results of NNDL accuracy tests for the four HEG/shower combinations shown in Table1.

Figure 7 :
Figure7: Results of NNDL accuracy tests for two combinations of HEG and shower whose contours do not match in the hard-collinear region, all with β ps = 1 2 .The showers are matched with the POWHEG scheme, but the vetoing procedure given in section 3.3 is not applied, in order to highlight the NNDL discrepancy expected for observables with β obs < β ps .The expected value from Eqs. (3.9)-(3.10) is shown in dotted blue (and marked by the blue arrow labelled "exp.").

Figure 8 :
Figure 8: Results of NNDL accuracy tests with n f = 0 for the PanGlobal and PanLocal showers (β ps = 1 2 ), using the same shower as both a HEG and for the subsequent parton showering steps, but with different choices of the de-symmetrisation parameter in the gluon splitting function, w heg = 1, and w ps = 0 respectively.The expected value of the NNDL discrepancy for β obs = 0 from Eqs. (3.20a)-(3.20b) is shown in dotted blue.

Figure 9 :
Figure 9: Thrust (left), Cambridge ln y 23 (middle) and SoftDrop ln k t /Q (right) distributions, unmatched (red) and matched (blue).They are obtained with a LL shower (our PanScales implementation of the Pythia 8 shower (PSPythia 8, top row)) and two NLL showers: PanGlobal with β ps = 0 (middle row) and PanLocal β ps = 1 2 (bottom row).The last row also shows the impact of HEG-style matching without the veto discussed in section 3.3.Dotted lines show x hard variation, while dashed lines show x r variations.

Figure 10 :
Figure 10:Ratio of multiple shower/matching scheme combinations to the multiplicatively-matched PanGlobal β ps = 0 shower in a bin 0.14 < 1 − T < 0.15 of the thrust distribution.The error bands represent the scale uncertainty and are colourcoded differently for each "main" shower, while the background colour reflects the accuracy of the matching/shower combination, i.e.LL (red), LL+NLO 2j /NDL ev.shp.or NLL+NLO 2j /NDL ev.shp.(yellow), NLL (blue), and NLL+NLO 2j /NNDL ev.shp.(green).The thin dotted vertical lines indicate the size of the scale uncertainties for the reference results, i.e. multiplicatively-matched PanGlobal β ps = 0. Recall from section 3.1 that the NLL notation indicates a breaking of exponentiation.

Figure 12 :
Figure 12: An illustration of the spin-sensitive observable ∆Ψ 12 as the difference of the azimuthal angles of the branching plane of the matched emission P 1 and the branching plane of a subsequent collinear branching P 2 .

Figure 14 :
Figure14: Same as Figure13, but for all possible configurations of H → gg.Note that in the lower plots, it is the gluon of the underlying 3-jet system that undergoes collinear splitting, rather than the hard emission from the 2-jet system (which is a quark or antiquark, and so does not mediate azimuthal correlation).This is also the reason for the lack of symmetry in the matched result between positive and negative η 1 values.