Antenna Showers with One-Loop Matrix Elements

We consider the probability for a colour-singlet qqbar pair to emit a gluon, in strongly and smoothly ordered antenna showers. We expand to second order in alphaS and compare to the second-order QCD matrix elements for Z ->3 jets, neglecting terms suppressed by 1/NC^2. We give a prescription that corrects the shower to the matrix-element result at this order for both soft and hard emissions, thereby explicitly reducing its dependence on evolution- and renormalization-scale choices. We confirm that the choice of pT for both of these scales absorbs all logarithms through order alphaS^2, and contrast this with various alternatives. We include these corrections in the VINCIA shower generator and study the impact on LEP event-shape and fragmentation observables. An uncertainty estimate is provided for each event, in the form of a vector of alternative weights.


Introduction
Modern QCD descriptions of hard-scattering events at particle colliders can be roughly divided into two broad categories. In the first, fixed-order descriptions, matrix elements are computed for all allowed initial states with a given final state, F , plus a limited number of additional partons. The leading-order (LO) description has the minimal number (often zero) of additional partons. For improved accuracy, one includes matrix elements with one extra parton beyond leading order and one loop correction (next-to-leading order) and so forth. The squared matrix elements are numerically integrated over the allowed phase space, after accounting for any divergences. Given the accuracy, i.e. order of the description, the possible number of additional final-state particles is in essence predetermined, and can take one (LO), or two (NLO) etc, values. In the second, parton-shower descriptions, one also starts from matrix elements for the desired hard process, F , but additional radiation is now generated stochastically via a shower algorithm, which is essentially Markovian. This is a unitary process, with probability one, and therefore does not change the probability of the underlying hard process to occur. The number of final-state partons is now not predetermined, and can take an infinity of different values. The two approaches have complementary strengths and weaknesses (for a review, see e.g. [1]). When hard extra emissions (e.g., hard jets) are important to model well, one looks to descriptions of the first category. However, the calculation is then unpredictive for near-collinear and soft radiation (e.g., jet substructure and soft wide-angle jets). The obverse holds for the second category.
Even from this very cursory summary, it is clear that methods to unite the twocombining strengths and eliminating weaknesses -are very important. Two longstanding and very successful approaches for combining one-loop matrix elements with parton showers are mc@nlo [2][3][4][5][6] and powheg [7][8][9]. An important restriction of both of these is that only the spectra of the basic hard-scattering partons are corrected to NLO precision, while those of additional QCD emissions are not. Removing this limitation, fully or partially, has been the focus of much recent effort [10][11][12][13][14][15][16][17][18][19], and is also among the main goals of this paper. While most other approaches employ parton-shower algorithms which are based on 1 → 2 splitting kernels, we develop an approach for matching NLO descriptions to showers based on 2 → 3 splittings [20]. The equivalent of the 1 → 2 splitting kernels are, for us, dipoleantenna functions [20][21][22]. At the practical level, our approach is in the context of the vincia framework [23,24]. Whatever the splitting kernel, the parton-shower approaches rest upon the factorization of both phase space and matrix element when the splitting is either soft or collinear, or both. A technical advantage of our approach is that the (n + 3)-particle phasespace factorizes exactly into a (n + 2)-particle phase space times a 2 → 3 phase space with all momenta on-shell, without need for momentum reshuffling [25]. Phase-space factorization and the antenna-based matrix-element factorization are important to our approach in about equal measure.
The essential bottleneck in such combinations of fixed order and parton shower is how to avoid double counting both real emissions as well as virtual effects. A key aspect of this is how well the NLO emissions are mimicked by parton-shower emissions. Emissions generated by a parton-shower Markov chain in fact produce approximations to tree-level matrix elements up to arbitrary numbers of legs, while the no-emission Sudakov factors generate the equivalent all-orders loop corrections 1 . This all-orders resummation of contributions is ordered in a measure of jet resolution, called the evolution scale, which we denote Q E . It is typically chosen to be a measure of transverse momentum [20] or invariant mass. Its fundamental role is to separate resolved from unresolved emissions, in analogy to a jet-clustering measure. The different evolution variables each have their strengths, depending on the context. As part of our present study, we judge them by how well their fixed-order expansions approximate the NLO matrix elements.
The main purpose of this paper is to define, for e + e − initial states, an antenna-based shower algorithm that incorporates multileg NLO corrections for both soft and hard emissions, and to study the quality of the matching for a variety of evolution variables. We strive for next-to-leading logarithmic (NLL) accuracy, in a way we shall detail below.
The leading-logarithmic (LL) structure of antenna showers was discussed in [27,28], with explicit comparisons of various algorithms to tree-level O(α 2 s ) matrix elements presented in [24,28,29]. A prescription for matching the showers to reproduce tree-level matrix elements exactly (over all of phase space) was developed in [24], with uncertainty variations provided in the form of vectors of alternative weights for each event. In [30] and [31] substantial speedups of the matching algorithm were obtained by dividing phase space into so-called sectors, and by deriving a formalism for using individual helicity amplitudes to correct the shower evolution, respectively. To further probe the structure of antenna showers, at the subleading-logarithmic level, we shall here consider the expansion of exclusive 2 → 3 splitting probabilities to O(α 2 s ), comparing these to one-loop matrix elements [32] and to corresponding second-order antenna functions [22,33].
We shall compare six different types of ordering criteria for the shower evolution: 1) strong ordering in transverse momentum, 2) strong ordering in dipole virtuality, 3 & 4) strong ordering in two variants of emission energy (mostly intended as cross-checks), and 5 & 6) so-called smooth ordering in p ⊥ and in dipole virtuality, as defined by [24]. We also consider several different choices for the renormalization scale µ R used in the tree-level antenna functions and discuss how to systematically absorb contributions proportional to the β-function by this choice, elaborating on earlier arguments [34,35].
Finally, we will present a prescription for how to systematically incorporate the secondorder (one-loop) qq → qgq antenna into the shower evolution, for each of the studied evolution variable choices. This will essentially constitute the NLL accuracy mentioned above. The resulting shower algorithm, whose qq → qgq splitting probability should therefore be correct to O(α 2 s ) over all of phase space, has been implemented in the publicly available vincia plug-in [23] to the pythia 8 event generator [36].
We have organized the paper as follows. In section 2 we discuss introductory aspects of (antenna) shower algorithms, define the various evolution variables, and the implementation of an ordering prescription that rules the shower evolution. In section 3 we present our matching prescription in detail, initially for 3-parton final states in Z-decay, then generally for n partons. In section 4 we discuss details of the Sudakov integrals required in the matching prescription and compare the infrared limits of those integrals to those of the one-loop matrix elements. In section 5 we combine one-loop and tree-level corrections in a single algorithm, perform speed benchmarking, and study the impact on LEP observables, especially in the context of α s (m Z ) extractions. We conclude in section 6 and elaborate on technical aspects in the appendices.

Antenna Showers
In this section, we recap the basic antenna-shower formalism, as used in the vincia shower algorithm. This also serves to introduce the basic notation and conventions that will be used in later sections.

The Formal Basis of Antenna Showers
Antenna showers are based on the factorization of (squared) colour-ordered QCD amplitudes in soft and collinear limits, which can be expressed as follows |M (. . . , p i , p j , . . .)| 2 i||j → g 2 s C P (z) s ij |M (. . . , p i + p j , . . .)| 2 (2.1) |M (. . . , p i , p j , p k , . . .)| 2 j g soft → g 2 s C A g (s ij , s jk , s ijk ) |M (. . . , p i , p k , . . .)| 2 , (2.2) with g 2 s = 4πα s the strong coupling and the subscript g in the second line emphasizing that the soft limit is only relevant for gluons.
In the collinear limit (first line), P (z) are the Altarelli-Parisi splitting kernels [37], z is the energy fraction taken by parton i (with a fraction (1 − z) going to parton j), and C is a colour factor, which we discuss below. This limit forms the basis for traditional parton showers, such as those in the pythia generator [38].
In the soft-gluon limit (second line), the function A has dimension GeV −2 , and is called an antenna function. For unpolarized massless partons 2 , its leading term is the so-called eikonal or dipole factor, where s ik = s ijk − s ij − s jk for massless partons. It was found early on that this factor can be reproduced by a traditional parton shower by imposing the requirement of angular ordering of subsequent emissions [40]. This gave rise to the angular-ordered showers [41,42] in the herwig and herwig++ generators [43,44] as well as the imposition of an angular-ordering constraint [38,45] in the jetset and pythia generators [36,46]. In fixed-order calculations, dipole [47] and antenna [21,22,25] functions are frequently used to define subtraction terms. These functions include additional subleading terms, beyond the eikonal one, which are necessary to correctly describe both soft and collinear limits in all regions of phase space. In the parametrization we shall use, their most general forms, for the branching process IK → ijk, are for gluon-emission and gluon-splitting processes, respectively, with the parent antenna invariant mass, m 2 IK = (p I + p K ) 2 = (p i + p j + p k ) 2 and the scaled invariants,

IK
; y jk = s jk m 2 IK , (2.6) and we use the notation δ ig = 1 if parton i is a gluon and zero otherwise. The functions F Emit and F Split allow for the presence of non-singular terms, which are in principle arbitrary. A logical choice would be F = 0, but this would not be invariant under reparametrizations of the antenna functions across the gluon-collinear singular limits [24]. Since the F functions can anyway be made useful in the context of uncertainty estimates [23,24], we therefore leave them as functions whose forms we are free to choose. In the soft-gluon limit, the eikonal factor is reproduced by the first term in in eq. (2.4). In the collinear q → qg limit, the AP splitting kernel also is reproduced. For collinear g → gg and g → qq branchings, one must sum over the contributions from two neigbouring antennae, which together reproduce the AP splitting kernel. Limits that are both soft and collinear are also correctly reproduced [22].
In the antenna context, the colour factors are 2C F for qq → qgq, C A for gg → ggg 3 , and 2T R for gluon splitting to qq, again using the normalization convention adopted in [24]. However, for qg → qgg there is a genuine subleading ambiguity whether to prefer, say, 2C F , C A , or something interpolating between them [48]. At fixed order, the question of subleading colour can in fact be dealt with quite elegantly, by using C A for all antennae and then including an additional qq antenna with a negative colour factor, −C A /N 2 C , spanned between the two endpoint quarks, for each qg . . .q chain [49]. In the context of an antenna-based shower, however, it is desirable to use only positive-definite antenna functions, and a prescription for absorbing the negative one into the positive ones was given in [24]. In the context of this work, however, we shall largely ignore subleading-colour aspects and, unless explicitly stated otherwise, assign a colour factor C A to the qg → qgg antenna function, thereby overcounting the collinear limit in the quark direction by a factor C A /(2C F ) 1 + 1/N 2 C . The renormalization scale used to evaluate the strong coupling in the antenna function, g 2 s = 4πα s (µ PS ), is typically chosen proportional to p ⊥ (following [34]). As alternatives, we shall also consider µ 2 PS ∝ m 2 D = 2min(s ij , s jk ), and, as an extreme case which connects with fixed-order calculations, the invariant mass of the antenna, µ 2 PS ∝ m 2 IK . A final aspect concerns the phase-space factorization away from the collinear limit. Within the framework of collinear factorization (and hence, in traditional parton showers), the momentum fraction, z, is only uniquely defined in the exactly collinear limit; outside that limit, the choice of z is not unique. In addition, a prescription must be adopted for ensuring overall momentum conservation, leading to the well-known ambiguities concerning recoil strategies (see e.g. [1]). In antenna showers, on the other hand, the antenna function is defined in terms of the unique branching invariants, s ij and s jk , over all of phase space, and the phase space itself has an exact Lorentz-invariant and momentum-conserving factorization, with dΦ Ant = 1 16π 2 m 2 IK ds ij ds jk dφ 2π (2.8) for massless partons (for massive ones, see [39]), with the φ angle parametrizing rotations around the antenna axis, in the CM of the antenna. Note the equality signs; no approximation is involved in this step. The only remaining phase-space ambiguity, outside the singular limits, is present when specifying how the post-branching momenta are related to the pre-branching ones. This is defined by a kinematics map, the antenna equivalent of a recoil strategy, which we here take to be of the class defined by [21,23].

Constructing a Shower Algorithm
In a shower context, the amplitude and phase-space factorizations above imply that we can interpret the radiation functions (AP splitting kernels or dipole/antenna functions) as the probability for a radiator (parton or dipole/antenna) to undergo a branching, per unit phase- where we use Φ as shorthand to denote a phase-space point. (If there are several partons/dipoles/antennae, the total probability for branching of the event as a whole is obtained as a sum of such terms.) An equally fundamental object in both analytical resummations and in parton showers is the Sudakov form factor, which defines the probability for a radiator not to emit anything, as a function of the shower evolution parameter, Q (i.e., similarly to a jet veto, with Q playing the role of the jet clustering scale; we return to the choice of functional form for the shower evolution scale in section 2.3). In the all-orders shower construction, these factors generate the sum over virtual amplitudes plus unresolved real radiation, and hence their first-order expansions play a crucial role in matching to next-to-leading order matrix elements. We here recap some basic properties. The Sudakov factor, giving the no-emission probability between two values of the shower evolution parameter, Q 1 and Q 2 (with Q 1 > Q 2 ), is defined by where it is understood that the integral boundaries must be imposed either as step functions on the integrand or by a suitable transformation of integration variables, accompanied by Jacobian factors. This description has a very close analogue in the simple process of nuclear decay, in which the probability for a nucleus to undergo a decay, per unit time, is given by the nuclear decay constant, dP (t) dt = c N . (2.11) The probability for a nucleus existing at time t 1 to remain undecayed before time t 2 , is This case is especially simple, since the decay probability per unit time, c N , is constant. By conservation of the total number of nuclei (unitarity), the activity per nucleon at time t, equivalent to the "resummed" decay probability per unit time, is minus the derivative of ∆, In QCD, the emission probability varies over phase space, hence the probability for an antenna not to emit has the more elaborate integral form of eq. (2.10). By unitarity, the resummed branching probability is again minus the derivative of the Sudakov factor, with Q 2 (Φ) the shower evolution scale (typically chosen as a measure of invariant mass or transverse momentum, see section 2.3), evaluated at the phase-space point Φ.
In shower algorithms, branchings are generated with this distribution, starting from a uniformly distributed random number R ∈ [0, 1], by solving the equation, for Q 2 . For an initial distribution of "trial" branching scales, we do not use the full antenna function, eq. (2.4), as the evolution kernel, but only its leading singularity, where p ⊥A is the ariadne definition of transverse momentum [50], which is also the one used in vincia. This reflects the universal 1/p 2 ⊥ behaviour of soft-gluon emissions. In addition to the trial scale, Q, two complementary phase-space variables are also generated (which we usually label ζ and φ [24]), according to the shape of A T over a phase-space contour of constant Q. From these, the model-independent set of trial phase-space variables (s ij , s jk , φ) are determined by inversion of the defining relations Q(s ij , s jk ) and ζ(s ij , s jk ), and the full kinematics (i.e., four-momenta) of the trial branching can then be constructed [23].
To decide whether to accept the trial or not, we note that the function A T differs from the eikonal in eq. (2.3) by the replacement of s ik in the numerator by m 2 IK . By accepting the trial scales generated by A T with the probability the eikonal approximation can be recovered, by virtue of the veto algorithm [1,51,52]. Crucially, any other function that has the eikonal as its soft-collinear limit could equally well be imprinted on the trial distribution by a similar veto. Two particularly relevant choices are the full physical antenna function, eq. (2.4) (which includes additional collinear-singular terms in addition to the eikonal) and the GKS-corrected antenna function (which also incorporates a multiplicative factor that represents tree-level matching in vincia), with A Emit given in eq. (2.4) and R n the n-parton tree-level GKS matching factor [24], to which we return in section 3.1. Note that, for gluon-splitting antenna functions (Xg → Xqq), we use Q = m qq , with a trial function ∝ 1/m 2 qq , and again implement the physical antenna function, eq. (2.5), and LO matching corrections by vetos. We also include the so-called ariadne factor, P Ari , defined by 20) with s N the invariant mass squared of the colour neighbour on the other side of the splitting gluon and s P = m 2 IK the invariant mass squared of the parent (splitting) antenna. This does not modify the singular behavior (as will be elaborated upon below), and was shown to give significantly better agreement with the Z → qqqq matrix element in [30].
Explicit solutions to eq. (2.15) using the trial function defined by eq. (2.16) were presented in [24], for fixed and first-order running couplings. In the context of the present work, twoloop running has been implemented using a simple numerical trick, as follows: given a value of α s (M Z ), we determine the corresponding two-loop value of Λ 2−loop QCD . We then use that Λ value in the one-loop solutions in [24], and correct the resulting distribution by inserting an additional trial accept veto: . (2.21) Due to the faster pace of 2-loop running, α 2−loop s (Q, Λ) < α 1−loop s (Q, Λ), hence this accept probability is guaranteed to be smaller than or equal to unity.
A final point concerns if there are several "competing" radiators (equivalent to several competing nuclei, and/or several competing available decay channels for each nucleus). In this case, the trial with the highest value of Q is selected (corresponding to the one happening at the earliest time, t), and consideration of any other branchings (decays) are postponed temporarily. After a branching, any partons involved in that branching are replaced by the post-branching ones, and any postponed trial branchings involving those partons are deleted. The evolution is then restarted, from the scale Q of the new configuration, until there are no radiators left with trial branching scales larger than a fixed, lower, cutoff, normally identified with the hadronization scale, Q had ∼ 1 GeV.

Evolution and Ordering
In order to solve eq. (2.15) we need to specify the form of eq. (2.10), which takes us from one scale Q 2 1 to a lower scale Q 2 2 . We change variables to parametrize the integral by the ordering variable, Q, and another, complementary (but otherwise arbitrary), phase-space variable which we denote by ζ. The generic evolution integral now reads with |J| denoting the Jacobian of this transformation. For branchings involving gluon emission, we consider three possible choices for the ordering variable: dipole virtuality m D , transverse momentum, and the energy of the emitted parton, E * j (in the CM of the parent antenna), with the following definitions, with the energy fraction x j = 2E * j /m IK . All three options are available as ordering variables in the vincia shower Monte Carlo. They are illustrated in figure 1, where contours of constant value of y E = Q 2 E /m 2 IK are shown for each variable, as a function of y ij and y jk . For completeness, we show both the case of a linear (top row) and quadratic (bottom row) dependence on the branching invariants, for each variable. Since the ordering variable raised to any positive power will result in the same relative ordering of emissions within a given antenna, the distinction between linear and quadratic forms does not affect individual antenna Sudakov factors. It does, however, affect the "competition" between different antennae, and the choice of restart scale for subsequent evolution after a branching has taken place, as will be discussed further below.
In labeling the columns in figure 1, we have also emphasized that mass-ordering, as defined here, corresponds to choosing the smallest of the daughter antenna masses as the "resolution scale" of the branching, whereas p ⊥ and energy correspond to using the geometric and arithmetic means of the daughter invariants, respectively. Naively, each of these could be taken as a plausible measure of the resolution scale of a given phase-space point. We shall see below which ones lead to better agreement with the one-loop matrix elements.
We consider two possible definitions for the complementary phase-space variable ζ, We emphasize that the choice of ζ has no physical consequences, it merely serves to reparametrize the Lorentz-invariant phase space. We may therefore let the choice be governed purely by Contours of constant value of y E = Q 2 E /m 2 IK for evolution variables linear (top) and quadratic (bottom) in the branching invariants, for virtuality-ordering (left), p ⊥ -ordering (middle), and energy-ordering (right). Note that the energy-ordering variables intersect the phase-space boundaries, where the antenna functions are singular, for finite values of the evolution variable. They can therefore only be used as evolution variables together with a separate infrared regulator, such as a cut in invariant mass, not shown here. convenience, and, for each antenna integral, select whichever of the above definitions give the simplest final expressions. The corresponding Jacobian factors, for each of the evolutionvariable choices we shall consider, are listed in tab. 1.
Note that, for the special case of the m 2 D and m 4 D variables, which contain the non-analytic function min(y ij , y jk ), the ζ definitions in eqs. (2.26) and (2.27) apply to the branch with y ij > y jk . For the other branch, y ij and y jk should be interchanged. With this substitution, the Jacobians listed in tab. 1 apply to both branches 4 . For branchings involving gluon splitting, g → qq, we restrict our attention to two possibilities, ordering in p ⊥ , defined as above, and ordering in gluon virtuality, defined as (for gluon splitting) . (2.28) Note that the normalization factors for the ordering variables have in all cases been chosen such that the maximum value of the ordering variable is m 2 IK . Since the phase space for subsequent branchings is limited both by the scale Q E of the previous branching (according to strong ordering) and by the invariant mass of the antenna m j , the effective "restart scale", after a branching in a strongly ordered shower, is given by for each antenna j.
Depending on the choice and value of Q, one or both daughter antennae after a splitting may have a non-trivial restriction on the phase space available for subsequent branching. Conversely, if Q > m j , there is no such restriction. Physically, we distinguish between the case in which the strong-ordering condition implies a non-trivial constraint on the evolution of the produced antennae, eating into the phase-space that would otherwise be accessible, and the case in which the strong-ordering condition does not imply such a constraint.
The regions of qq → qgq phase space in which either zero, one, or both of the daughter antennae (qg and gq respectively) are constrained by the ordering condition are illustrated in figure 2, for each of the choices of evolution variable under consideration. The black shaded areas correspond to regions in which both the qg and gq antennae are restricted, by having Q < m j . The darker gray shaded areas show regions in which only one of the antennae is restricted, while the other will still be allowed to evolve over its full phase space. In the light-gray shaded areas, both of the antennae are allowed to evolve over all of their available phase spaces, equivalent to the ordering condition imposing no constraint on the subsequent evolution. We recall that we are here discussing the upper boundary on the subsequent evolution, hence the infrared 5 (IR) poles are not affected.
To further clarify the meaning of the plots in figure 2, let us discuss panel (e) as an example. The coordinates, (y ij , y jk ), represent the 3-parton state before it evolves to a 4parton state, and each point corresponds to a specific value of the evolution variable at hand, cf. figure 1. Assuming ordering in p ⊥ and using subscript (3) for quantities evaluated in the 3-parton state, the value of the evolution variable for a specific (y ij , y jk ) point is The further evolution of the shower, from a 3-to a 4-parton state, involves a sum over all possible branchings of the qg and gq  antennae. Consider the qg one. Its branchings can again be characterized by two invariants (s 1 , s 2 ), both of which will be smaller than m 2 qg . However, depending on the value of m 2 qg (or, equivalently, y ij ) the p ⊥ of the new configuration, 4p 2 ⊥(4) = 4s 1 s 2 /m 2 qg may actually be larger than 4p 2 ⊥(3) . In a strongly ordered shower, such configurations are not allowed, and would be discarded. Whether this situation can occur or not, for one or both of the qg and gq antennae, as a function of (y ij , y jk ), is what figure 2 reveals, for each type of ordering variable.
The mathematical consequence is that, in the dark-and black-shaded regions, respectively, the upper boundary of one or both of the qg and gq antenna integrals is set by the evolution variable, rather than by phase space. This creates an important difference between the integrals generated by a shower algorithm and those used for IR subtractions in traditional fixed-order applications for which the integrals often run over all of phase space, although some subtraction schemes feature parameters that allow restrictions on the phase space for the subtraction terms [53,54]. In particular, we see that the strong-ordering condition will generate additional logarithms involving s ij /Q 2 E(3) as argument. For a "good" choice of evolution variable, these logarithms should explicitly cancel against ones present in the one-loop 0 1 θ m 2 ant − (s ij +s jk ) 2 s Table 2. Boundaries corresponding to the ordering variables portrayed in figure 1, with m 2 ant corresponding to the active 3 → 4 dipole s qg or s gq , and s = m 2 Z at the Z pole. We have chosen ζ 2 as the energy sharing variable for m D and p ⊥ ordering and ζ 1 for E * ordering, with ζ defined as in eq. (2.26) and eq. (2.27). The energy variable will lead to infinities if the hadronization scale is not imposed as a cut-off. matrix elements, a question we shall return to in detail in section 4.
Several interesting structures can be seen in figure 2. Firstly, the linearized variables imply less severe constraints on the subsequent evolution than the quadratic ones. This is easy to understand given that the linearized variables, Q lin , are related to the quadratic ones, Q qdr , by and hence Q lin > Q qdr , implying a higher absolute restart scale for the linearized ordering variables.
It is also apparent that, for a given choice of linearity, mass-ordering reduces the phasespace for further evolution more than p ⊥ -ordering does, which in turn is more constraining than energy-ordering. In this comparison, however, it becomes important to recall that the traditional ordering variables used, e.g., in vincia, are the linearized mass-ordering and the quadratic p ⊥ and energy-ordering variables 6 . Within that group, p ⊥ -ordering appears to be the most restrictive, followed by energy-ordering, with (traditional, linearized) mass-ordering leading to the most open phase space for the subsequent evolution.
We are now able to fully specify the boundaries of the evolution integrals in eq. (2.22). For each Q E contour (see figure 1), the integration limits in ζ are listed in tab. 2. Combined with a Q E interval and an antenna function, these boundaries account for the integrated tree-level splitting probability when going from one scale Q 2 1 to another Q 2 2 , as expressed by eq. (2.22). The last column in tab. 2 tells when the 3 → 4 ordering condition is active. In figure 2 this corresponds to a region darkening due to the restriction, with its shade determined by the amount of restricted dipoles.
Finally, we note that the dependence on Q in eq. (2.29) causes explicit non-Markovian behavior at the 4-parton level and beyond, since the value of Q then depends explicitly on which branching was the last to occur. A more strictly Markovian variant of this is obtained by letting the min() function act on all possible Q values (corresponding to all possible colourconnected clusterings) of the preceding topology. In that case, a single Q value can be used to characterize an entire n-parton topology, irrespective of which branching was the last to occur. Since the distinction between Markovian and non-Markovian shower restart conditions only enters starting from the 4 → 5 parton evolution step, it will not affect our discussion of the second-order 2 → 3 branching process. For completeness, we note that the strongly ordered showers in vincia are of the ordinary non-Markovian type, while the smoothly ordered ones are Markovian.

Smooth Ordering
In addition to traditional (strongly ordered) showers, we shall also consider so-called smooth ordering [24]: applying the ordering criterion as a smooth dampening factor instead of as a step function. This is not as radical as it may seem at first. Applying a jet algorithm to any set of events will in general result in some small fraction of unordered clustering sequences. This is true even if the events were generated by a strongly ordered shower algorithm, if the jet clustering measure is not strictly identical to the shower ordering variable. An example of this, for strong ordering in p ⊥ and in m D , clustered with the k T algorithm, can be found in [55].
In smooth ordering, the only quantity which must still be strictly ordered are the antenna invariant masses, a condition which follows from the nested antenna phase spaces and momentum conservation. Within each individual antenna, and between competing ones, the measure of evolution time is still provided by the ordering variable, which we therefore typically refer to as the "evolution variable" in this context (rather than the "ordering variable"), in order to prevent confusion with the strong-ordering case. The evolution variable can in principle still be chosen to be any of the possibilities given above, though we shall typically use 2p ⊥ for gluon emission and m qq for gluon splitting.
In terms of an arbitrary evolution variable, Q, the smooth-ordering factor is [55] where Q is the evolution scale associated with the current branching, andQ measures the scale of the parton configuration before branching. A comparison to the strong-ordering step function is given in figure 3, on a log-log scale.
In the strongly-ordered region of phase-space, Q Q , we may rewrite the P imp factor as Applying this to the 2 → 3 antenna function whose leading singularity, eq. (2.16), is proportional to 1/Q 2 , we see that the overall correction arising from the Q 2 /Q 2 and higher terms is finite and of order 1/Q 2 ; a power correction. The LL singular behaviour is therefore not affected. However, at the multiple-emission level, the 1/Q 2 terms do modify the subleading logarithmic structure, starting from O(α 2 s ), as we shall return to below. In the unordered region of phase-space, Q >Q, we rewrite the P imp factor as which thus effectively modifies the leading singularity of the LL 2 → 3 function from 1/Q 2 to 1/Q 4 , removing it from the LL counting. The only effective terms ∝ 1/Q 2 arise from finite terms in the radiation functions, if any such are present, multiplied by the P imp factor. Only a matching to the full tree-level 2 → 4 functions would enable a precise control over these terms. Up to any given fixed order, this can effectively be achieved by matching to tree-level matrix elements, as will be discussed in section 3.1. Matching beyond the fixed-order level is beyond the scope of this paper. We thus restrict ourselves to the observation that, at the LL level, smooth ordering is equivalent to strong ordering, with differences only appearing at the subleading level. The effective 2 → 4 probability in the shower arises from a sum over two different 2 → 3 ⊗ 2 → 3 histories, each of which has the tree-level form thus we may also perceive the combined effect of the modification as the addition of a mass term in the denominator of the propagator factor of the previous splitting. In the strongly ordered region, this correction is negligible, whereas in the unordered region, there is a large suppression from the necessity of the propagator in the previous topology having to be very off-shell, which is reflected by the P imp factor. Using the expansion for the unordered region, Figure 4. Illustration of scales and Sudakov factors involved in an unordered sequence of two 2 → 3 branchings, representing the smoothly ordered shower's approximation to a hard 2 → 4 process. eq. (2.33), we also see that the effective 2 → 4 radiation function, obtained from iterated 2 → 3 splittings, is modified as follows, in the unordered region. That is, the intermediate low scaleQ, is removed from the effective 2 → 4 function, by the application of the P imp factor. Finally, to illustrate what happens to the Sudakov factors, we illustrate the path through phase space taken by an unordered shower history in figure 4. An antenna starts showering at a scale equal to its invariant mass, √ s, and a first 2 → 3 branching occurs at the evolution scalê Q. This produces the overall Sudakov factor ∆ 2→3 ( √ s,Q). A daughter antenna, produced by the branching, then starts showering at a scale equal to its own invariant mass, labeled √ s 1 . However, for all scales larger thanQ, the P imp factor suppresses the evolution in this new dipole so that no leading logs are generated. To leading approximation, the effective Sudakov factor for the combined 2 → 4 splitting is therefore given by in the unordered region. Thus, we see that a dependence on the intermediate scaleQ still remains in the effective Sudakov factor generated by the smooth-ordering procedure. Sincê Q < Q in the unordered region, the effective Sudakov suppression of these points might be "too strong". The smooth ordering therefore allows for phase space occupation in regions corresponding to dead zones in a strongly ordered shower, but it does suggest that a correction to the Sudakov factor may be desirable, in the unordered region. A study of Z → 4 jets at one loop would be required to shed further light on this question. Having presented introductory aspects of (antenna) showers, we now turn to a detailed discussion of how we match them to higher fixed-order calculations.

Tree-Level Matching
The strategy for matching to tree-level matrix elements used in vincia was defined by GKS in [24] and is tightly related to the veto algorithm outlined above. The philosophy is to view the shower produced by the smoothly ordered antenna functions as generating an all-orders approximation to QCD, which can be systematically improved, order by order, by including one more factor in the accept probability, called the matrix-element correction. For a given trial branching, the full trial accept probability, up to the highest matched number of partons, is then obtained as a product of the ordinary trial-accept probability in the shower, multiplied by this extra correction factor.
Since the shower is already correct in the soft and collinear limits, the matrix-element correction factor will tend to unity in those limits, but it can deviate on either side of unity outside those limits. As long as the combined accept probability is still smaller than unity, a probabilistic accept/reject step can still be applied, without having to worry about reweighting the events (which would be required if the total accept probability should exceed unity). It is also important to define the factor only in terms of physical cross sections (here represented by LO matrix elements), which guarantees positivity. (Again, if it were allowed to become negative, one would have to introduce negative-weight events, but this is avoided in the GKS strategy as defined in [24]).
As we have seen, the shower furnishes an all-orders approximation to QCD. The aim is, for each resolved parton/jet multiplicity, to match the answer provided by the shower to an, ideally, all-orders exact expression, by applying a multiplicative correction factor, schematically [24,38] Matched = Approximate Exact Approximate . (3.1) At tree level, we in fact know only the first term in the expansion of the numerator, and we therefore expand the shower approximation to the same level. For n partons, assuming the approximation has already (recursively) been matched to the preceding order, where the subscript "T" indicates trial quantities (cf. section 2.2), we have suppressed the dependence on phase-space points, Φ, and the subscript j in the (n−1)-parton matrix element indicates the configuration obtained by performing the inverse shower step that removes parton j from the n-parton state. The factors in eq. (3.3) are easy to calculate if a tree-level matrix-element (ME) generator is available to provide the |M | 2 factors. The total ME-corrected accept probability is thus simply eq. (2.19), (3.4) As mentioned above, this factor should be positive and smaller than unity, in order to avoid having to reweight any events. In practice, we have found the trial function defined in eq. (2.16) to guarantee this for all processes we have so far considered, mainly consisting of Z → n and H → n partons. As shown in [24], it is also possible to absorb subleading-colour corrections into this matching factor in a positive-definite way, but since subleading colour goes beyond the scope of our study we do not reproduce the arguments here. The fact that these factors change the distribution of the final set of generated events to the correct matrix-element answer can be explained by following the steps of the algorithm and summing over all possible branching histories. We start from Born-level matrix-element events, and generate trial shower branchings, for a trial approximation to the (B + 1)-parton matrix element of: with i running over all possible single-parton clusterings that correspond to allowed shower branchings. Applying the LO accept probability, eq. (3.4), changes this to That is, summed over shower histories, numerators and denominators are designed to cancel exactly, leaving only the LO matrix element for B+1 partons, as desired. Due to the full phasespace coverage and explicitly Markovian nature of the smoothly ordered shower algorithm, this procedure is straightforward to iterate for Born + 2, 3, etc partons 7 .
To provide a connection with antenna subtraction, which will be useful when we extend to NLO matching below, we can rewrite the ratio in eq. (3.1) by a trivial rearrangement, That is not the case for ordinary strongly ordered frameworks, due to the presence of dead zones in phase space and to the generally non-Markovian shower restart conditions. For such algorithms, addition of events in the dead zones [56], with CKKW-like Sudakov-factor prescriptions for multi-leg matching [57,58], would presumably be necessary.
The numerator in this equation is very similar to a standard antenna-subtracted matrix element, though we emphasize that our antennae are of course modified by the presence of the P imp and P Ari factors. Let us finally re-emphasize that since we apply the correction factor to the antenna functions themselves, thereby modifying the probability for a branching to occur, the probability for a branching not to occur is also modified. These corrections will therefore also be present in the Sudakov factors generated by the corrected shower evolution. The fact that the correction factor, R n , is unity in all LL singular limits (since the shower functions are guaranteed to match the matrix-element singularities there) implies that the LL properties of the Sudakov factors are not affected by this modification. However, the tower of subleading logarithms is changed, for better or worse. While it is known that finite terms do not exponentiate our corrections here also include a more subtle aspect, namely a resummation of the subleading logarithms present in the higher-order matrix elements. At this level, however, we cannot be sure that this procedure produces the correct subleading logarithms of a formally higherorder resummation. Therefore, we view it at present merely as an interesting, and hopefully beneficial, side-effect of unitarity-based matching. The examination of formally subleading terms carried out in this paper is a first step towards a more rigorous study of these aspects.

One-Loop Matching at the Born Level
For the Born level, at NLO, the GKS matching strategy [24] reduces to the powheg one [7][8][9]. We nonetheless begin by recapitulating the steps used to derive the one-loop correction to the Born-level matrix element, in our notation. We then continue to higher multiplicities.
As our basis for one-loop matching we take the tree-level strategy described in section 3.1. Since the corrections are applied as modifications to the branching probabilities, without adding, subtracting, or reweighting events, the total inclusive rate after tree-level matching to any number of partons, is still just the leading-order, Born-level one. By the same token, after one-loop matching, at the integrated level, the total NLO correction to the inclusive rate must therefore just be the NLO "K-factor", For processes like Z decay, where the NLO correction has no dependence on the Born-level kinematics, this is trivial to implement as an overall reweighting factor on the Born-level events, where we have introduced the notation V for the NLO correction term, anticipating a similar notation for the multileg case below. Note that one could equally trivially normalize to NNLO or to data, as desired for the application at hand (we note though that such a normalization choice does not, by itself, ensure NNLO precision for any quantity besides the total inclusive rate).
However, when the amount of final state particles exceeds two, the NLO correction depends on the Born-level kinematics, therefore it is worth illustrating the general procedure for deriving a fully differential K-factor, for each phase-space point. This also serves as a useful warm-up exercise for the multi-leg case below.
At NLO, we may distinguish between inclusive and exclusive rates for the first time. Either can in principle be used to derive matching equations between showers and fixedorder calculations, but the exclusive one is best suited for deriving expressions at the fully differential level. We recall that the exclusive n-jet cross section is defined as the cross section for observing n and only n jets, while the inclusive n-jet cross section counts the number of events with n or more jets. One therefore has the trivial relation with Q the resolution scale of whatever (IR safe) algorithm is used to define the jets.

Inclusive Born
The total inclusive rate produced by the tree-level matched shower is just the Born-level matrix element, where the subscript indicates the parton multiplicity (2 for Z → qq decay) and the superscript indicates the loop order beyond the Born level (0 indicates the Born loop order). Because cancellation of real and virtual corrections is exact in both the unmatched shower as well as in the tree-level matching scheme described above, there are no further corrections to consider for the inclusive rate. In other words, the total integrated cross section produced by the shower is obtained merely by integrating eq. (3.11) over all of the Born-level phase space. We now seek a correction term, V 2 , such that gives the correct inclusive NLO rate. From eq. (3.9), we know that the correction term for Z decay is A systematic way of deriving this result, which can be applied to arbitrary processes, is provided by considering the cross section at the exclusive level.

Exclusive Born
The shower expression for the exclusive Z → qq rate (defined at the hadronization cut-off, which is the lowest meaningful resolution scale in the perturbative shower) is where we have expanded the Sudakov factor ∆ to first order. Due to the presence of the hadronization scale, this expression is IR finite and can be defined in 4 dimensions. We remark here on the validity of this expansion in α s for the exclusive cross section. For the purpose of constructing the matching factor to order α s the expansion is a parametric one. In the ratio of the exact and approximate exclusive cross section, since the singularities match to the shower accuracy, divergences or large logaritms (depending on whether one choose zero or finite resolution scale) cancel and the resulting factor has a well-behaved expansion in α s .
The colour factor for qq → qgq is and we assume that the antenna function, A, is either the one derived from Z decay [20] or has been matched to it, using LO matching. That is, (3.16) We first consider the limit Q had → 0, in which case the expression becomes which can only be defined in the presence of an IR regularization scheme. We shall here use dimensional regularization, working in d = 4 − 2 dimensions. Below, we rederive the matching equations in 4 dimensions, for Q had = 0, and show that the same final matching factors are obtained in both cases. At NLO, the exclusive Z → qq rate at "infinite" perturbative resolution is where we have written the right-hand side in a form similar to eq. (3.17), in d dimensions.
Because the resolution scale has been taken to zero, there are no unresolved 3-parton configurations to include. The virtual matrix element is with the function I qq used to classify the divergences [22,33,59]. Note that we have modified the definition of I to make it explicitly dimensionless, see appendix A. On the shower side, the integral of the Z → qgq antenna in eq. (3.17) is [22] s 0 and, not surprisingly, the difference comes out to be exactly α s /π × |M 0 2 | 2 . Writing this correction as a multiplicative K-factor, we obtain eq. (3.9).
As a cross-check, we now repeat the derivation in 4 dimensions, reinstating the hadronization scale. The fixed-order side is then where the integral that has been added corresponds to unresolved 3-parton configurations, with A again given by eq. (3.16). Though eq. (3.14) is now defined entirely in 4 dimensions, we still need dimensional regularization to regulate the two last terms in the fixed-order expression. In principle, the integral in the last term could be carried out explicitly, but it is simpler to rewrite it as where the first term is just the full antenna integral, eq. (3.20), and the second term is identical to the one appearing in eq. (3.14), with which it cancels completely, cf. the definition of the tree-level matching, eq. (3.16). The final correction term is therefore again exactly equal to α s /π × |M 0 2 | 2 . Note that the scale and scheme dependence of the α s /π correction is not specified since its ambiguity is formally of order α 2 s . For definiteness we take the renormalization scale for this correction to be proportional to the invariant mass of the system, , with k inc µ thus representing the free parameter that governs the choice of renormalization scale for the total inclusive rate for Z → hadrons. We shall consider both one-loop and two-loop running options. The value of α s (m Z ) will be determined from LEP data in section 5.

One-Loop Matching for Born + 1 Parton
The approximation to the 3-parton exclusive rate produced by a shower matched to (at least) NLO for the 2-parton inclusive rate and to LO for the 3-parton one, is where M 0 3 is the tree-level Z → qgq matrix element and Q R3 denotes the "restart scale". For strong ordering, Q R3 is equal to Q 3 , while, for smooth ordering, it is given by the nested antenna phase spaces, i.e., by the successive antenna invariant masses. The subscripts on the two Sudakov factors ∆ 2 and ∆ 3 make it explicit that they refer to the event as a whole, see the illustration in figure 5. Again, we have the choice whether we wish to work in 4 dimensions, with a non-zero hadronization scale, Q had , or in d dimensions with the hadronization scale taken to zero. We have maintained the hadronization scale in eq. (3.23), though we shall see below that the dependence on it does indeed cancel in the final result.
The 2-parton Sudakov factor, ∆ 2 , is generated by the (matched) evolution from 2 to 3 partons, with A g/qq again defined by eq. (3.16). Notice that the integral only runs from the starting scale, m 2 Z , to the 3-parton resolution scale, Q 2 3 , hence this integral is IR finite, though it does contain logarithms. In the remainder of this paper, we shall work only with the leading-colour part of the Sudakov and matrix-element expressions, hence from now on we replace 2C F in the above expression by C A , The 3-parton Sudakov factor, ∆ 3 , imposes exclusivity and is given by where the index j runs over the qg and gq antennae, and we use subscripts E and S for gluon emission and gluon splitting, respectively. We have implicitly assumed smooth ordering here, which implies that the upper boundaries on the integrals are given by the respective dipole invariant masses (squared), s j . Note also that we must take into account all modifications that are applied to the LL antenna functions, including P imp , P Ari , and LO matrix-element matching factors. (We do not write out these factors here, to avoid clutter.) I.e., the antenna functions in the above expression must be the ones actually generated by the shower algorithm, including the effect of any modifications imposed by vetos. For strong ordering, there are no P imp factors, and the upper integral boundary is instead min(Q 2 3 , s j ), However, since strong ordering is not able to fill the entire 4-parton phase space [24,29], full NLO matching can only be obtained for the smoothly ordered variant. It is nonetheless interesting to examine both types of shower algorithms, since even in the strongly ordered case, we may compare the Sudakov logarithms arising at O(α 2 s ) to those present in the fixedorder calculation.
On the fixed-order side, the expression for the 3-parton exclusive rate is simply where the last term represents 4-parton configurations in which a single parton is unresolved with respect to the hadronization scale. For Z decay, d-dimensional expressions for the virtual matrix element have been available since long [22,32,33,60]. Details on the calculation and in particular its renormalization, are given in appendix B, in a notation convenient for our purposes.
We now seek a fully differential matching factor, K 3 = 1 + V 3 , such that the expansion of reproduces the exact expression, eq. (3.28), to one-loop order. ("Approximate" here stands for the tree-level matched shower approximation, eq. (3.23).) Trivial algebra yields where we have reinstated d-dimensional forms of the one-loop matrix element and of the divergent 3 → 4 terms. For a shower matched to |M 0 4 | 2 at leading order, the last two terms will cancel, by definition of the matched antenna functions (even for an unmatched shower, the difference could at most be a finite power correction in the hadronization scale, since the matrix element and the shower antenna functions have the same singularities), yielding: Rewriting the remaining integrals in terms of a set of standardized antenna subtraction terms, writing out the ordering functions for gluon emission and gluon splitting, O E and O S , explicitly, and denoting the ariadne factor for gluon splitting by P A , we arrive at the following master equation for the second-order correction to the 3-jet rate: with the standardized Gehrmann-Gehrmann-de Ridder-Glover (GGG) subtraction terms defined by [22]: whose IR limits and integrated pole structures were examined thoroughly in [22,33,59], though we have here rewritten the IR singularity operators I (1) in explicitly dimensionless forms, see appendix A. (The alphabetical labeling in eqs. (3.33) corresponds to the notation used in [22].) The first line combined with the first term on the second line in eq. (3.32) represent a standard antenna-subtracted one-loop matrix element, normalized to the Born level, with the standardized subtraction terms tabulated in eq. (3.33), and the additional finite term V 2Z originating from the NLO matching at the preceding order; see section 3.2, eq. (3.13).
The subsequent terms express the difference between the simple fixed-order subtraction carried out in the first line and the actual terms that are generated by a matched Markovian antenna shower. Physically, these terms represent the difference between the evolution of a single dipole (the original qq system) and evolution of two dipoles (the post-branching qgq system), with a transition occuring at the branching scale Q 3 . As mentioned above, the O Ej and O Sj factors in the third line represent the ordering criterion imposed in the evolution, either strong or smooth. For smooth ordering, they are with Q Ej the evolution variable used for gluon emissions, while for strong ordering, the factor (1 − O j ) can be removed if the integral boundaries are replaced by [Q 2 3 , s j ] (note: this replacement should only be done in the third line).
The last term in eq. (3.32) is an artifact of the ariadne factor, P Ari , which was introduced in section 2.2 and is applied to gluon-splitting antennae in vincia. Summed over the two "sides" of the splitting gluon, this produces the same collinear singularities as the standard gluon-splitting antenna, but in highly asymmetric configurations in which the splitting gluon is near-collinear to a neighbouring colour line, the ariadne factor produces a strong suppression, which improves the agreement with the tree-level 4-parton matrix element [30], and which then generates an additional logarithm.
Notice that all but the δA terms are defined in terms of standarized antenna functions, and the corresponding integrals can be carried out analytically, once and for all. We give explicit forms for each of these terms, for each choice of evolution variable, in the following section.
The only terms of eq. (3.32) that need to be integrated numerically are thus the δA terms, which express the difference between the standardized antenna functions and those generated by the actual (matched) shower evolution, which may have different finite terms and/or be matched to the LO 4-parton matrix element. Nonetheless, since the previous lines already contain most of the structure, we expect these functions to be relatively well-behaved and numerically sub-leading. Specifically, the δA terms for gluon emission and gluon splitting, respectively, are defined by with A LL the unmatched shower antenna function (as defined in [30,39]) and the second-order LO matching factors, R 4E and R 4S (for Z → qggq and Z → qq q q, respectively), defined as in eq. (3.4), but including only the leading-colour terms in R LC 4E . For strong ordering, similarly to above, the O j factors can be removed by changing the integration boundaries of the δA terms to [0, Q 2 3 ]. Finally, we note that one could in principle equally well have defined eq. (3.32) without the terms on the third line. The δA terms in eqs. (3.36) and (3.37) would then likewise have to be defined without P imp and P Ari factors. However, while this would give a seemingly cleaner relation with standard fixed-order subtraction, the behaviour of the (numerical) integrals over the δA terms would be more difficult, due to over-subtraction in the unordered regions. (Showers without either a strong-ordering condition or a smooth-ordering suppression greatly overestimate the real-radiation matrix elements in the unordered region [24,28,30].) Numerically, it is therefore more convenient to integrate the contributions represented by the third line in eq. (3.32) analytically, leaving only the suppressed terms in eq. (3.37) to be integrated over numerically.
To be specific, the numerical integration over the δA terms is performed by rewriting the δA integrals in dimensionless MC form, as: 38) and similarly for the gluon-splitting terms, with Φ i a number of random (uniformly distributed) antenna phase-space points. The common factor 1/4 arises from combining the prefactor 8π 2 above with the area of the phase-space triangle, 1/2, and the factor 1/(16π 2 ) from the phase-space factorization, dΦ ant . For smooth p ⊥ -ordering with an arbitrary normal- , the ordering factors, O j , reduce to: where we have used y with ijk indices for the scaled invariants in the original qgq topology and x with abc indices for the integration variables in the antenna phase space. Note also that the y values are normalized to the full 3-parton CM energy (squared), while the x values are normalized to their respective dipole CM energies (squared).

The Renormalization Term
A further ingredient to be discussed is the choice of renormalization scale on both the fixed order and parton shower sides of the calculation, as these scales are in general chosen differently in both sides. Hence a translation term arises at second order, which must account for this difference, keeping in mind that, as the scale evolves from one to the other, flavour thresholds are passed. Our aim is to have the flexibility to use fixed order matrix elements renormalized according to their usual scheme, while maintaining the freedom to make a different choice on the shower side. The fixed order calculations for Z-decay to jets to which we match are customarily renormalized in a version of the MS scheme called the Zero-Mass Variable Flavour Number Scheme (ZM-VFNS). In this scheme the bare QCD coupling is renormalized as and n F is the number of light flavours. One thus ignores flavours that are heavier than the scale of the calculation, both in the virtual and in the real calculations. Once all the UV poles are cancelled, one has a running coupling that depends on the number of light flavours for the scale µ R at hand. One then changes the flavour number when a threshold is crossed. For our present case of Z boson decay to jets we take n F = 5 for µ R not too different from the Z-boson mass.
Let us be more specific about the matching of α s across flavour thresholds. At one loop, The value of Λ F depends on the number of active flavours, as follows. When passing flavour thresholds the following one-loop matching conditions are imposed These conditions can be satisfied if Λ F obeys the matching conditions With these conditions one can also express α s values for different flavour numbers into ea- . (3.47) For completeness we briefly review how this n F -dependent UV singularity occurs in the context of the (inclusive) 3-jet rate, in the case where we only consider massless quarks [32,60].
In the virtual contribution, the only one-loop diagram for Z → qqg that is sensitive to the number of flavours is the quark self-energy correction on the external gluon. The self-energy diagram itself, being scaleless, is zero in dimensional regularization. However, renormalization of the coupling amounts to adding a n F counterterm on the exteral gluon line proportional to The real contribution contributes a n F dependent (collinear) 1/ pole as well, from gluonsplitting In the sum over real and virtual contributions the poles cancel, as guaranteed by the KLN theorem, leaving a logarithm of the form On the shower side a related prescription is used, in which the running coupling is evaluated at a shower scale µ PS , such that the scale again depends on the number of flavours.
Depending on the value of µ PS , a corresponding value of n F is chosen, as well as of the QCD scale Λ F . This is often different from that for a fixed order calculation. To give a specific example, matrix elements will typically be renormalized at a scale characteristic of the total CM energy, i.e., µ 2 ME = s an event-independent value, while resummation arguments imply one best chooses a running scale, such as µ PS = p ⊥ , for shower applications [34,35], which can differ per event.
Shifting to a different scale for α s of a given flavour number is quite straightforward. Translating from a shower scale µ PS to a matrix-element scale µ ME amounts to replacing, for an antenna function A further aspect is that shower Monte Carlos normally switch to 4-flavour (3-flavour) running for scales µ < m b (µ < m c ), matching the α s value across the thresholds to obtain a continuous running. For a consistent treatment, such thresholds must be taken into account when translating α s from the shower scale to the matrix-element one. At one-loop order (which is all that is relevant for the NLO expansion), this can be done by inserting an additional term for each flavour threshold in the region [µ PS , µ ME ], with m thres the flavour threshold. Physically, eq. (3.51) expresses running with n F flavours all the way from µ PS to µ ME . The correction term, eq. (3.52), expresses that the number of flavours was effectively lower below each flavour threshold passed on the way. Note that this can also be used to account for a larger number of flavours in the shower calculation, e.g., at scales µ PS > m t , with the sign change of the correction then automatically reflected by the logarithm.
For coherent parton-shower models, the arguments presented in [35] also motivate a change to a "Monte Carlo" scheme for α s , in which Λ QCD is rescaled, for each n F , by the so-called CMW factor ∼ 1.5 (with some mild flavour dependence), relative to its MS value. If the shower model being matched employs this scheme, then a further rescaling of the renormalization-scale argument, µ PS → µ PS /k CMW , should be used in eq. (3.51), with The translation of renormalization scale (and scheme) yields then an additional term to be added to the definition of V 3 in eq. (3.32), plus any additional flavour-threshold correction terms, cf. eq. (3.52). By inserting these terms, which enter at overall order α 2 s ln(µ 2 ME /µ 2 PS ), the two calculations can be compared consistently at one-loop accuracy.
Note that if several different shower paths populate the same fixed-order phase-space point, then each path will in general be associated with a distinct µ PS value. Thus, one V 3µ term arises for each shower path, weighted by the relative contribution of each path to the total. Since for our case there is only one antenna contributing to Z → qgq, this particular complication does not arise here.
We finally alert the reader regarding the use of different flavour number α s 's in the master equation (3.32). In that equation cancellation of 1/ divergences take place, already in the first line of the right hand side. For this cancellation it is important that the subtraction terms, originating from the shower expansion and listed in eq. (3.33), use α (5) s renormalized as in the fixed order calculation. All subsequent terms in the master equation are finite, and constitute differences of unordered and strongly ordered shower based terms, which are also finite, and beyond NLO.

Leading-Colour One-Loop Correction for Z → 3 Jets
Combining the results above, in particular eqs. (3.32), (3.33), and (3.54), we obtain the complete expression for the leading-colour 8 one-loop correction for Z → 3 Jets, where: 8 We use the usual MC definition of leading colour and include terms ∝ CA and ∝ nF but neglect ones ∝ 1/CA.
• the first line contains the full (leading-colour) one-loop matrix element, the V 2Z correction from one-loop matching at the preceding order, and the V 3µ term from the choice of shower renormalization scale; • the second line contains the standardized subtraction term arising from the qg → qgg and gq → ggq antennae; • the third line contains the standardized subtraction term arising from the qg → qq q and gq →q q q antennae; • the fourth to last lines contain the terms arising from the difference between the (matched) shower evolution and the standardized subtraction terms, including the consequences of ordering choices and modification factors such as those arising from the Ariadne factor and from matching to the LO matrix elements.
We denote the singular subtracted 1-loop matrix element by SVirtual In section 4, we compute the analytical integrals corresponding to each of the showergenerated terms, for different choices of evolution variable, ordering criterion, and antenna functions.
With the one-loop matrix element expressed as in appendix B.2, it is easy to see that the infrared singularity operators in eq. (3.56) cancel, leaving only explicitly finite remainders (which may still contain logarithms of resolved scales). This then constitutes the description of the one-loop matching for Z → 3 jets, having already discussed the case for two jets. In the context of eq. (3.10) we have now corrected the first two terms on the rhs to NLO accuracy.
We round off with a few remarks on the normalizations of the various Z → n-parton rates that are obtained by our procedure, since this is a point on which the various approaches to multileg NLO corrections differ. We make the following observations: 1. The total inclusive Z decay rate: the matrix-element correction scheme derived in this paper maintains a strict unitarity between the real and virtual corrections that are applied beyond Born level. An important consequence is that the total inclusive Z decay rate is not changed by switching on the V 3Z correction 9 . 9 The theoretically most sensible choice would be to normalize the inclusive Z → 2-parton rate to the full NNLO result, but at the level we work at, one could equally well normalize to NLO or to data. In either case, the total normalization of the generated sample is left unchanged by the V3Z correction.
2. The inclusive Z → 3 jets rate: both the virtual (one-loop) correction to the 3-jet rate and the real (tree-level) correction to the 4-jet rate are included here. Hence the inclusive 3-jet rate is NLO correct. Without these corrections, it would only be LO correct. Thus, the 3-jet inclusive rate does change when switching on the V 3Z term 10 .
3. The exclusive Z → 2 jets rate: The strict unitarity imposed by our correction method implies that every 3-jet event begins life as a 2-jet one. Since the the V 3Z term modifies the probability for a 2-jet event to evolve to become a 3-jet one, at the O(α 2 s ) level, the 2-jet exclusive rate receives an equal and opposite correction. This represents an O(α 2 s ) ambiguity on the exclusive 2-jet rate, which is not adressed in our paper (though it could be removed by normalizing to the full NNLO result for Z → 2, cf. the inclusive 2-jet rate above). 4. The exclusive Z → 3 jets rate: for a given 3-parton configuration, the evolution to 4 partons and beyond is not changed by V 3Z (though it is changed by the application of the 4-parton tree-level corrections, which we take to be included throughout this paper). Thus, while switching V 3Z on does change the total amount of 3-jet events (cf. the inclusive 3-jet rate above), it does not directly change the fraction of those events which will develop a fourth or more jets.

One-Loop Correction for Born + 2 Partons
To illustrate how the formalism presented here generalizes to higher multiplicities, we take the case of the NLO correction to Z → 4 partons. For simplicity, however, we continue to restrict our analysis of the correction factor to the leading-colour level. At NLO, the exclusive Z → 4 partons rate at "infinite" perturbative resolution (similarly to above) is Labeling the 4 partons by Z → i, j, k, , there are two possible antenna-shower histories leading to each 4-parton configuration, with j and k the last emitted parton, respectively. Those two contributions both enter in the definition of the tree-level 4-parton matching factor, such that their sum reproduces the full 4-parton matrix element. Note that a separate such factor is applied to Z → qggq and Z → qq q q, and that we have suppressed colour and coupling factors here, for compactness (we ignore the small, non-singular extra interference terms for the special case where all four quarks have the same flavour). The antenna functions, A, are understood to include all such factors, as well as any P imp and P ari factors appropriate 10 In section 5.2 below (comparisons to LEP measurements), this is seen most easily in figure 15, where the "(NLO off)" curves undershoot the "(NLO on)" ones, for observable regions dominated by 3 or more jets.
to the branchings at hand. For a general n-parton matrix element, the denominator contains one term for each possible clustering. Labeling the IK → ijk history by A and the JL → jk one by B, the sum over the two histories yields where it is understood that α is an index, not a power, and the last product factor takes into account the NLO matching at the preceding multiplicities. Expanding the Sudakov factors to first order and using the definition of the NLO correction factor at the preceding multiplicity, eq. (3.31), this becomes 60) which we can rewrite as where we again emphasize that the antenna functions are understood to include all relevant coupling, P imp , and P ari factors. The first term represents the new subtraction that the shower generates at 4 partons, while the second represents part of the NLO correction to the preceding multiplicity. For one of the histories (the one followed by the "current" event), this correction has already been evaluated and can be reused. The contribution from the other history will have to be recomputed, however. In general, there will be one subtraction to perform at the n-parton level, and there will be m ∼ n − n Born − 1 new subtractions that have to be done at the (n − 1)-parton level, in addition to the one that was already done during the evolution of the current event.
Clearly, there is an undesirable scaling behavior built into this, which will make NLO matching at many partons quite computing intensive. An alternative, which eliminates the sum over histories, is that of sector showers, see e.g., [30,61]. Though this is not the main avenue pursued in this paper, we nevertheless give some comments below on how a sectorbased NLO matching scheme could be constructed.

One-Loop Matching for Sector Showers
The matching conditions derived above may also be applied to so-called sector showers [30,61], with a few relatively minor modifications. The expansion of the Sudakov factors generated by the LO matched shower will now contain integrals over ratios of matrix elements (which are the LO matched sector antenna functions), multiplied by sector vetos. The presence of the sector vetos makes analytical phase-space integration more difficult.
However, since the sector approach merely represents a different way of decomposing the same singularities as the global one, we may effectively recycle the integrals carried out for the global case by adding and subtracting the terms produced by a smoothly ordered "standard" shower (i.e., using the GGG functions). The first four lines of eq. (3.32) then remain unchanged. The definition of the δA terms in the last line, however, changes to for the terms arising from the qgq → qggq Sudakov factor, with analogous ones arising from the gluon-splitting contributions. The step function, Θ j , represents the sector veto applied to the sector antenna functions. The sector antenna function, up to the tree-level matched orders, is just (3.63) Since the individual sector and global antenna functions have different singularity structures (they are only guaranteed to have the same singularities at the summed level), the integrals in eq. (3.62) are divergent, and cannot be carried out numerically. In order to obtain numerically convergent integrals, we must divide up the contributions of the global terms onto each sector, and perform a set of correlated integrals in which the singularities explicitly cancel in the divergent limits, where each line now corresponds to the sum of contributions to a single sector, for which the difference between sector and global antennae is finite. The individual integrals are of course still divergent, but they can now be treated numerically by collecting the terms on each line under a single integral sign. Analytically, this is complicated since the two integrals on each line are not associated with the same kinematics map 11 . Numerically, however, we may still ensure a point-by-point cancellation in the singular limits by keeping the two integrals formally separate, but carrying them out simultaneously, in a correlated way, as follows. For each antenna, generate a random uniformly distributed phase-space point,

Sudakov Integrals
In this section, we work out the standardized Sudakov integrals appearing in the second and third line of eq. (3.32), for each choice of evolution variable. We also study the soft and collinear limits of the Sudakov integrals and compare them to those of the one-loop matrix element. This provides an explicit check of whether the first-order expansion of the Sudakov factors generates the correct logarithms present in the fixed-order calculation. Given our choice of the GGG antenna functions as our standard ones, the relevant terms are 1) The general form of the first term, which originates from the 2 → 3 branching step, is where the definitions for the K i and the I i functions are given in appendix C, for each type of antenna function and ordering variable. Their derivation and soft/collinear structure will be discussed more closely below, for each choice of ordering and evolution variable. The form of the 3 → 4 integrals depends on whether we work in the context of strong or smooth ordering.
We shall now consider each of those cases in turn, beginning with strong ordering.

Strong Ordering
For strong ordering, the inverted ordering conditions in eq. (3.32), (1 − O E j /S j ), reduce to step functions expressing integration over the unordered region. The integration surface is thus limited from below by the phase-space contour defined by the evolution scale of the first branching, Q 2 , and from above by the edge defined by the invariant mass of the antenna. The expression generated by the 3 → 4 splitting case for gluon emission is (4.3) where K i and I i are the same as those for the 2 → 3 term above, though they here appear with different arguments. The remaining case is the 3 → 4 gluon splitting defined by with H defined in appendix C and P A j as defined in eq. (2.20). We will discuss the derivation of these terms in more detail in the following three subsections, for strong m D -, p ⊥ , and energy-ordering, respectively.

Dipole Virtuality
We begin with dipole virtuality as evolution variable, which is perhaps the simplest case. We start by repeating the integrals of eq. (3.32) with the one-particle phase space defined as in eq. (2.8  . (4.7) The integration is over a triangular surface. The lower integration boundary cuts off the evolution variable at the value of the 3-parton evolution scale. The other boundary is determined by the total energy of the dipole before branching which here is √ s. We use the integral To discuss the 3 → 4 Sudakov terms, let us for definiteness assume that we are in a 3-parton phase-space point with s qg > s gq . The opposite case is symmetric. Again we only include the T = 2 terms explicitly here, with the details of the full antenna integrals relegated to appendix C. (4.10) The integration is again over a triangular surface. The total energy of the dipole before branching is now s qg . The integral in eq. (4.10) corresponding to the sum over antenna integrals only contains one d 0 3 integral because the other has equal upper and lower integration boundaries. Note that this integral actually vanishes for s qg ≤ Q 2 3 , which amounts to the dipole-virtuality ordering allowing the 3 → 4 branchings to populate their full respective phase spaces (i.e. no correction term is necessary).
Focusing on the case s qg > 2s gq for which the second integral is nonvanishing (which now amounts to the ordering condition imposing a nontrivial restriction on the 3 → 4 phase space, see figure 2.3), we obtain, including the 2 → 3 term with lower-transcendentality terms again available in appendix C. For the mirror case s gq > 2s qg essentially symmetric expressions are obtained, while for the intermediate cases in which the two invariants are within a factor 2 of each other, the second integral in eq. eq. (4.11) simply vanishes. The full double-logarithmic term from the expanded Sudakov terms in eq. (4.5), for strong ordering in dipole virtuality, is then where the Θ function ensures that the second term is only active if s max = max(s qg , s gq ) > 2 min(s qg , s gq ) = Q 2 3 , (4.14) so that the expression applies over all of phase space.
We shall now consider the infrared limits of this result, and compare them to those of the one-loop matrix element. For this comparison we keep only terms that involve logarithms of the invariants. The soft limit corresponds to vanishing Q 2 3 (ξ min → 0). The first line of eq. (4.13) represents the contribution of the 2 → 3 expanded Sudakov. To find the contribution in the soft limit, we choose to approach the limit along the diagonal of the phase space triangle. Parametrizing this by s qg /s = s gq /s → y we find for this term The contributions of the 3 → 4 Sudakovs in the soft limit are examined in two separate cases corresponding to the two regions in figure 2.3. In the first case given by s max < 2s min , corresponding to the light grey area in the figure, the step function in eq. (4.13) yields zero. In the second case given by s max > 2s min , corresponding to the dark grey area in the figure, the step function is equal to one. The double logs and dilogarithms now yield a finite contribution that does not diverge in the soft limit. We can understand this by parametrizing the soft limit by λ s qg = λs s gq = pλs s 1 = λκs s 2 = λµs p > 2 , This implies that the integration variable scales with the integration limits and is independent of the soft limit. We can also expect this behaviour from examing figure 2.3. The shape of the different regions does not change for different values of Q 2 3 , in contrast with the case of transverse momentum, as we will see below.
After the poles cancel in eq. (3.55), the pole-subtracted version of the one-loop matrix element, SVirtual , defined in eq. (3.56), contains all the relevant terms on the matrix-element side. The transcendentality-2 terms of SVirtual are given by − R(y 1 , y 2 ) = Li 2 (y 1 ) + Li 2 (y 2 ) − π 2 6 − ln y 1 ln y 2 + ln y 1 ln(1 − y 1 ) + ln y 2 ln(1 − y 2 ) . (4.17) Including the transcendentality-1 terms (see appendix B), taking the soft limit by sending s qg /s = s gq /s = y → 0, and keeping only logarithmic terms, the pole-subtracted matrix element (ME) reduces to ME: The single log proportional to C A originates from the renormalization term and the single log of the closed quark loops (proportional to n F ) arises due to the definition of the infrared singularity operator, defined in the appendix in eq. (A.3).
Taking the same limit of the Sudakov integrals for dipole virtuality eq. (4.5), but omitting for the time being the renormalization term, V 3µ , we find for the parton shower (PS),

−PS:
s qg /s = s gq /s = y → 0 α s C A 2π ln 2 y + 3 2 ln(y) . We see that the soft limit almost cancels against eq. (4.18). For an NLL-accurate shower, however, all divergent terms should match precisely, leaving at most a finite remainder in the final matching correction, eq. (3.55). In the expressions above, this holds for the ln 2 (y) term but not for the single logarithms (different coefficient). Interestingly, the remainder is proportional to the QCD β function, specifically ME − PS → − α s 2π It can therefore be absorbed in the choice of renormalization scale by solving for µ PS in V 3µ , which yields: µ 2 PS ∝ y s . This tells us that, in the soft limit, the specific choice of a renormalization scale that is linear in the branching invariants will absorb all logarithms up to and including α 2 s ln(y). Interestingly, this reasoning would rule out µ 2 R ∝ p 2 ⊥ , since our p ⊥ -definition is quadratic in the invariants, p 2 ⊥ = s ij s jk /s. A better choice of renormalization scale would appear to be µ R ∝ m D , specifically Taken at face value, this seems to contradict the standard literature [34] on p ⊥ as the optimal universal renormalization-scale choice. However, as we shall see below in figure 6, there is in fact no choice of renormalization scale that absorbs all logarithms for this particular evolution variable; the choice µ R ∝ m D merely manages to reabsorb the additional logarithms that are generated by the ordering condition as y → 0, but leftover logs in other parts of phase space will remain uncanceled, ruining the NLL precision. In that sense, choosing µ R ∝ p ⊥ would simply leave a different set of uncanceled logs, nonvanishing as y → 0. Before we show the results over all of phase space however, we first investigate a complementary interesting limit, the hard-collinear one, which is characterized by one of the invariants becoming maximal while the other vanishes. In this limit, the pole-subtracted one-loop matrix element, SVirtual , becomes ME: There are no log-squared terms in this limit, and both of the single-log terms are half as large here as they were in the soft limit. The Sudakov integrals for m D -ordering yield one divergent term, − 1 6 C A ln(y), in the hardcollinear region, modulo a factor α s /(2π). The Sudakov integral for gluon splitting in the neighbouring antenna, represented by the first term on the second-to-last line of eq. (3.55) is specified in the last line of eq. (4.5). The step function is only non-zero for the first term in the hard-collinear limit s qg → s, s gq → 0 and produces a term 1 6 P A j n F ln(y). The numerator of the corresponding Ariadne factor contains the invariant of the neighboring dipole s gq which vanishes in this limit and causes the dipole splitting contribution to reduce to zero. The n F -dependent contribution is instead shifted to the last term of eq. (3.55), which has the same limit but without the Ariadne pre-factor. The hard-collinear limit of the shower terms, including only terms involving logarithms of the invariants and not including the V 3µ term, is therefore Again, the combination (ME−PS) relevant for computing the correction factor is proportional to the QCD β function, and in fact has exactly the same form as eq. (4.20). The conclusion is therefore that, also in this limit, all logarithms through α 2 s ln(y) can be absorbed by choosing a renormalization scale which is linear in the vanishing invariant. The particular choice which is linear in both the soft and collinear limits is µ PS ∝ m D . To illustrate this, we show the full NLO Z → 3 jets correction factors, (1 + V 3Z ), for m D -ordering with a few different choices of renormalization scale and scheme, in figure 6. Note that the axes are logarithmic, in ln(y ij ) = ln(s ij /s), to make the infrared limits clearly visible.
Without the V 3µ term, the correction factor looks as depicted in the top left-hand plot in figure 6. The increasing contours towards the axes indicate uncanceled logarithms in the correction factor. The middle pane shows the correction factor derived for µ PS = p ⊥ . As discussed above, there is an uncanceled logarithm in the soft limit (lower left-hand corner of the plot), since p ⊥ is quadratic in the vanishing invariants there. However, in the hardcollinear limits (upper left-hand and lower right-hand corners), p ⊥ is linear in the vanishing invariant, and hence the contours remain bounded there. In the right-hand pane, we show the choice µ PS = m D , which can be seen to lead to bounded correction factors (below ∼ 1.3) in all three phase-space corners. Nonetheless, there is still an uncanceled divergence between the soft and hard collinear limits. We shall see in the section on p ⊥ -ordering below that the cure for this is basically to choose a better evolution variable.
In the bottom row of figure 6, we show a few variations on µ PS = m D , specifically we include the CMW rescaling of Λ QCD , as defined by eq. (3.53), and show how a variation of the renormalization scale by a factor of 2 affects the correction factor. In the left-hand pane, we show µ PS = 1 2 m D and in the right-hand one µ PS = m D . Of these, the choice µ PS = 1 2 m D , with CMW rescaling, leads to the smallest correction factors (best LO behaviour), and this could therefore be taken as a useful default for m D -ordering, e.g. for uncertainty estimates.

Transverse Momentum
For a shower ordered in p ⊥ , the antenna phase-space integrals in eq. (3.32) are performed over contours such as those depicted for p T squared in figure 1. The curved contours motivate a coordinate transformation from (s 1 , s 2 ) to a basis defined as the dimensionless evolution variable y = Q 2 s = 4s 1 s 2 s 2 , complemented by an energy-sharing variable, which we define as z = s 1 s . Note that the coordinate transformation depends explicitly on the total invariant mass s of the 2 → 3 dipole. For the 3 → 4 integrations, the invariant mass s is replaced by the invariant mass of the appropriate dipole (either s qg or s gq ). The integration boundaries in z are determined by the intersection of the invariant mass of the dipole with the evolution parameter Q 2 . The choice of y and its integration boundaries make the effect of strong ordering especially clear since we see integration from Q 2 to the total invariant mass of the dipole (the unordered region). As before, the integration over the gluon-splitting antenna e 0 3 makes use of a different phase space integration, in m qq , and only uses the evolution parameter as a cut-off in the singularity of the corresponding dipole.
The contributing terms are: For n = 1 we set m 2 IK = s, for n = 2 m 2 IK = s qg and for n = 3 m 2 IK = s gq . The Ariadne factor P A j is defined in eq. (2.20). The max condition on the outer integration boundary of A 2 and A 3 reflect that the correction term disappears if the generated Q 2 3 is larger than the invariant mass of the dipole. As for m D -ordering, we here work out the most divergent behavior explicitly, by focussing on the double log terms arising from the eikonal term in the antenna, and relegate the full form of the antenna integrals to appendix C. The double poles give rise to terms which lead to the following generic transcendentality-2 integrals, The double logarithm in the shower expansion is generated by a combination of the 2 → 3 and 3 → 4 Sudakov integrals, with the respective pieces adding up to We see that a partial cancellation arises between the first two terms (which come from the 2 → 3 Sudakov expansion) and the last three (which come from the 3 → 4 expansion). What remains is a log squared in both invariants ln (s qg /s) ln (s gq /s). At the single-log level, the 3 → 4 terms give a numerically larger coefficient than the 2 → 3 ones, leading to a single log remainder. The gluon-splitting term also reduces to a single log. The overall result in the soft limit is then − PS: Comparing with the result of the virtual correction in the soft limit, eq. (4.18), we see that the shower generates the double log terms correctly, and, similarly to the case of m D -ordering, there is a single-log remainder which is proportional to the QCD β function. However, for p ⊥ -ordering the constant of proportionality is 1, rather than 1 2 , a difference which translates to the optimal renormalization-scale choice being quadratic in the invariants in this case, rather than linear. Before commenting further on this difference, let us first consider the complementary, hard-collinear, limit.
In the hard-collinear limit, we find the same as for m D -ordering, −PS: . Integration over the s gq dipole, however, is associated with an Ariadne factor carrying s qg in the numerator and therefore reduces to zero. As before, we can write the remainder as half the QCD β function, which implies that a renormalization scale linear in the vanishing invariants can absorb the logarithm.
To summarize, for p ⊥ -ordering we find that the optimal renormalization-scale choice must be quadratic in the vanishing invariants in the soft limit and linear in the hard-collinear limit. Those conditions are fulfilled by p ⊥ itself, thus here µ PS = p ⊥ which leaves the correction factor free of logarithms. Indeed, we see that the combination of evolution and renormalization in p ⊥ leads to a rather flat correction factor over all of phase space, showing that this combination is indeed "better" than m D -ordering.
In the bottom row of plots in figure 7, we include the CMW factor and show the correction factors for µ PS = p ⊥ (left) and µ PS = 2p ⊥ (right). In particular on the left-hand pane, the NLO correction factor is essentially unity in the soft limit, while the corrections in the hardcollinear regions remain less than ∼ 20%. This gives some additional weight to the arguments for p ⊥ -ordered showers with p ⊥ as renormalization scale being the best default choice for strongly ordered dipole-antenna showers. It also provides some rationale why one typically finds a rather large value of α s (m Z ) ∼ 0.13 (with CMW rescaling, or α s (m Z ) ∼ 0.14 without it) when tuning such models to LEP event shapes; there is still a genuine order 20% NLO correction in the hard resolved region (upper right-hand corner). We return to this in more detail in the context of full LO + NLO matching in section 5.

Energy
To put the differences between m D and p ⊥ in context, we now briefly examine the case of energy ordering, which is known to produce the wrong DGLAP evolution in the collinear limit [27,28,62], and hence we should find larger (possibly divergent) NLO corrections.
Slicing phase space with the energy variable Q 2 3 = s ijk (y ij + y jk ) 2 , see figure 1, requires the use of an explicit infrared cut-off because the contours otherwise allow for the invariants to hit singular regions for every value of the contour. We will here use a cut-off in transverse momentum (a cut-off in dipole virtuality is also possible). The cutoff motivates us to switch to a different choice of integration variables. Therefore integration is transformed from (s 1 , s 2 ) to the dimensionless evolution parameters y 2 E = Q 2 s = (s 1 +s 2 ) 2 m 2 IK and completed with the energy   sharing variable ζ = s 2 m 2

IK
. The interesting integrals arising from expanding the Sudakov form factor then are The inner integral has been rescaled to make it independent of the outer integral with ζ = y E ζ . To establish the cut-off, we use the relation 4 s 1 s 2 s 2 = 4p 2 ⊥ /s, which we demand to be larger than the cut-off ∆. In terms of our variables we then have the condition (4.34) The upper and lower limits on ζ are then Focussing on the eikonal integral the result for this integral is (4.37) In the soft limit y 2 E = 4y 2 → 0 this reduces to Thus we see that there are explicit non-cancelling double-logarithmic terms that involve the hadronization cutoff, ∆. Depending on the ratio between the dipole mass and the cutoff, these would lead to asymptotically divergent correction factors. One might wonder whether using a linearized form of energy ordering would make a difference, see figure 1. Rather than go through the derivations again, we merely show the full NLO corrections in figure 8, for both linear (top row) and squared (bottom row) energy ordering, for an (arbitrary) dimensionless cutoff value of ∆ = 10 −7 . From left to right in both rows, we show the three renormalization-scale choices, µ PS = p ⊥ (left), µ PS = m D (middle), and µ PS = Q E (right), with the latter equal to linear energy in the top row and squared energy in the bottom one. Interestingly, the contours in the linear case are increasing towards the soft region, while they ultimately decrease in the squared case. It is clear, however, that no intelligent choice of renormalization scale can absorb the infrared divergences. Moreover, any other choice of ∆ would have led to different contours, due to the explicit ln(∆) terms, hence even if a "least bad" choice was found, it would not be universal.
As mentioned above, the main point of showing these comparisons is to place the comparison between m D and p ⊥ in the previous subsections in perspective. Thus, while we saw that p ⊥ was generating a better-behaved correction factor than m D , the one for m D is still far better behaved than is the case for energy ordering. From this perspective, we still believe it could make sense, e.g., to use m D -ordering, with the NLO correction factor included, as a conservative uncertainty variation for a central prediction based on p ⊥ -ordering.

Smooth Ordering
We will now discuss the same Sudakov integrals as in the previous subsections but for the case of smooth ordering (section 2.4). This is especially interesting given that smooth ordering is the way vincia is able to fill all of phase space without significant under-or overcounting at the LO level [24]. As discussed in section 2.4, however, this does involve some ambiguity in what Sudakov factors are associated with the unordered points, and the NLO 3-jet correction factors should tell us explicitly whether this ambiguity generates problems at this level.
The Sudakov integrations are actually more straightforward for smooth ordering than was the case for strong ordering, since the P imp factor regulates the integrands on the boundaries. Therefore the integrations always run over the full phase space of the system. The 2 → 3 splitting generates the same terms as in the strong-ordering case, eq. (4.2). Including also the 3 → 4 terms, the expanded Sudakov generates the following antenna integrals, where Q 3 is the evolution scale evaluated on the 3-parton configuration and Q E j (m qq ) is the scale of the 3 → 4 emissions (splittings) being integrated over. The full answer for the 3 → 4 case for gluon emission is (4.40) where K i and L i can be found in appendix C. The full answer for the 3 → 4 case for gluon splitting is where G can be found in the appendix. We will discuss the derivation of these terms in more detail in the following two subsections, for smooth m D -and p ⊥ -ordering, respectively.

Dipole virtuality
Since the 2 → 3 emission terms remain the same as in the case of strong m D -ordering, we only need to rederive the 3 → 4 contributions to V 3Z , which are with Q 2 3 = 2 min(s qg , s gq ) and e 0 3 carrying the singularity in s 1 . We will focus again on deriving the transcendality-2 terms explicitly, with the full expressions given in the appendix. We start by recalling the expression for the strongly-ordered 2 → 3 branching, s .
To this we add the results from the eikonal term sgq sqg ). Taking the limit for the soft region y 2 3 → 2 (since we take the invariants as vanishing simultaneously), we see that the remainder is just a finite term that contains no logarithms of the vanishing invariants, We will receive this contribution twice. Including all divergent logarithmic contributions and disregarding constant terms such as in eq. (4.44) , we find the same as in the strong-ordering case, − PS:  and hence the preferred choice of scale in the soft limit remains one which is linear in the vanishing invariants, such as µ PS ∝ m D .
In the hard collinear limit the Sudakov integrals plus the 'Ariadne Log' reduce to − PS: again the same as in the strongly ordered case, cf. eq. (4.24).
To summarize, we do not expect any qualitatively different limiting behaviour in the smoothly ordered case with respect to the strongly ordered one, though details may of course still vary. To illustrate this, we include the plots in figure 9. In all cases, we use a renormalization scale ∝ m D , but with different prefactors, from left to right: µ PS = m D , µ PS = m D /2, and finally µ PS = m D /2 with CMW rescaling. In particular the latter gives correction factors very close to unity in both the soft and hard collinear limits, while we still see the leftover divergence inbetween those limits that was also present in the case of strong m D -ordering, cf. figure 6. Nonetheless, it is worth noting that for a large region of phase space, say with m ij > 0.1 m (corresponding to ln(y ij ) > −4.6), the corrections are still quite well behaved and relatively small, less than ∼ 20%.

Transverse momentum
Again we only need to recompute the contributions from the 3 → 4 Sudakov terms, as the 2 → 3 ones are the same as in the strongly ordered case. The 3 → 4 Sudakov integrals are As before we focus on explicitly calculating the transcendentality-2 contribution arising from the eikonal part of the antenna in the first term in the first line of eq. (4.47), where we have transformed y i = s i sqg and y 2 3 = Q 2 3 sqg . In the limit s min /s, s max /s = y → 0 so that y 2 3 → 0, this yields Adding the contributions from the 2 → 3 splitting and transcendentality-1 terms, we find the following result for the soft limit − PS: as in the strongly ordered case. The double logarithm matches with SVirtual and the single logarithm can be absorbed by choosing a renormalization scale that is quadratic in the vanishing invariants, such as µ PS ∝ p ⊥ . In the hard collinear limit, the shower integrals behave as −PS: the same as in all the other cases. This completes the argument that indeed µ PS ∝ p ⊥ is the appropriate choice also for smooth p ⊥ -ordering.
In figure 10, we show the NLO correction factors, (1 + V 3Z ), for smooth p ⊥ -ordering. The top row shows the correction factors without using the CMW rescaling of Λ QCD , and the plots in the bottom row include it. For the left-hand panes, we used a shower renormalization scale of µ PS = p ⊥ , and for the right-hand ones we used µ PS = 2p ⊥ .
We see that all correction factors are essentially well-behaved, with no runaway logs, similarly to the case of strong p ⊥ -ordering. However, for the case of smooth p ⊥ -ordering, it looks as if the CMW rescaling (bottom row) is almost doing "too much" in the soft region. Given that the CMW arguments [35] were derived explicitly by considering the subleading behaviour of strongly ordered (coherent) parton showers, we do not perceive of this as any major drawback. Instead, one should merely be aware of the slight shifts in the NLO corrections that result from applying it or not, recalling that a rescaling of Λ by the CMW factor ∼ 1.5 is within the ordinary factor 2 variation of the renormalization scale that is often taken as a standard for uncertainty estimates.
The shifts caused by CMW rescaling and/or by renormalization-scale prefactors are of  Figure 11. NLO correction factor for smooth p ⊥ -ordering, with GGG antennae: variations of how gluon splittings are interleaved with gluon emissions, see text. We used α s = 0.12, n F = 5, and µ PS = p ⊥ .
course fully taken into account in our implementation in the vincia code, and are thus reabsorbed into the one-loop matching coefficient at the matched order, stabilizing the prediction. Differences at higher orders will result from the fact that the CMW rescaling, if applied, is used for all shower branchings, while the NLO correction derived here is only applied at the Z → 3 stage of the calculation.
Because smooth p ⊥ -ordering is the default in vincia we wish to understand this case as best as we can, and therefore we include some further comparisons with non-default settings in figure 11.
In the left figure of figure 11, we modify the normalization of the evolution variable from the vincia default Q 2 E = 4p 2 ⊥ to the ariadne choice Q 2 E = p 2 ⊥ . Though the normalization factor cancels in the P imp factor for sequential gluon emissions, it is relevant for deciding the relative ordering between gluon emissions and gluon splittings. As this plot shows, however, the modification only produces quite small differences in the NLO correction factor, and with the "wrong" sign. Thus, we retain N ⊥ = 4 as the default in vincia. In middle figure of figure 11, we change the evolution variable for gluon splittings to be the same as that for gluon emissions, i.e., p ⊥ , with similar conclusions as for the previous variation. In the right figure of figure 11, we switch off the Ariadne factor for gluon splittings. We notice that the NLO corrections get slightly larger. There is no change along the diagonal y ij = y jk since the Ariadne factor is unity there, but along the edges of the plots, the NLO corrections become larger, which further motivates the choice of keeping the Ariadne factor switched on by default in vincia.
The overall result is that the infrared limits are generally well-behaved for p ⊥ evolution with µ PS ∝ p ⊥ . Remaining differences amount to small finite shifts of order 10%-20% away from unity. At that level, the effective finite terms of the antenna functions also play a role, hence it is too early to draw definite conclusions just based on the plots presented here. The impact of finite terms will be studied in section 5 in the context of matching to the LO matrix elements for Z → 4 partons, which effectively fixes the finite terms with respect to the pure-shower answers studied here.

Tables of Infrared Limits
The results of the preceding subsections on the infrared limits of the pole-subtracted matrix elements and of the Sudakov integrals generated by the various evolution-scale choices are collected here, in parametric form, for easy reference. The renormalization terms, V 3µ , are not included. Tab. 3 expresses the limits of SVirtual , while tab. 4 contains the Sudakov-integral limits. Table 3. Limits of SVirtual , with L denoting ln(y), where y parametrizes the limit in the soft region taken along the diagonal of the phase space triangle y = s qg /s = s gq /s → 0. The hard collinear limit only takes one invariant s qg /s or s gq /s to the soft limit while the other is set to 1. We have omitted an overall factor of α s /2π. Table 4. Limits of strong and smooth p ⊥ and m D ordering, with naming conventions as defined in tab. 3. Non divergent terms, such as π 2 have been omitted in the calculation of V 3Z , and the renormalization term in V 3Z is set to zero. An overall factor of α s /2π is suppressed.

Results including both LO and NLO corrections
In the preceding section, we focussed on deriving the analytic forms of the shower integrals and comparing their infrared limits to the matrix-element expressions. It is now time to include also the finite terms arising from matching to the 4-parton tree-level matrix element, expressed by the δA terms in eq. (3.55). Our ultimate aim in this section is to include the full leading-colour one-loop corrections through second order in α s (i.e., up to and including Z → 3 partons) and combine these with the full-colour tree-level corrections through third order in α s (i.e., up to and including Z → 5 partons, the default in vincia). However, since we shall perform the δA integrals numerically, adding those terms to the analytic ones derived in the previous section, we first wish to examine the numerical size and precision required on the δA terms themselves.

Finite antenna terms and LO matching corrections
Finite-term variations of the antenna functions (and in particular fixing the finite terms via unitary LO matching corrections, such as is done in vincia [24]) will affect the terms generated by the 3 → 4 Sudakov expansions in the following way. Larger finite terms will cause an increased amount of 3 → 4 branchings, which in turn will decrease the associated Sudakov factor (in the sense of driving it closer to zero). This will feed into the NLO correction factor, which compensates and drives the final answer back towards its NLO-correct value.
(Note that similar variations will not occur for the 2 → 3 branching step, since we treat that as fixed to the LO 3-parton matrix element throughout.) This feedback mechanism is encoded in the δA terms in eq. (3.55). Following the reasoning above, we should expect larger antenna finite terms to increase the NLO correction factor (since, to stabilize the 3-parton exclusive rate, it must compensate for losing more 3-parton phase-space points to 4-parton ones), and vice versa: smaller finite terms should result in a decrease of the NLO correction. At the pure-shower level (i.e., without LO matrix-element corrections to fix the finite terms), this is illustrated by figure 12. For ease of comparison, all plots use the CMW rescaling of Λ QCD , µ PS = p ⊥ , n F = 5, and α s (m Z ) = 0.12. The default antenna functions in vincia 12 are shown in the middle panes, for strong (top row) and smooth (bottom row) ordering, respectively. A variation with smaller finite terms for the 3 → 4 antenna functions is shown to the left, and one with larger finite terms on the right. As expected, the NLO correction factors react by becoming lower for smaller finite terms and higher for larger finite terms, for both strong and smooth ordering.
We emphasise that the plots in figure 12 are shown purely for illustration, to give a feeling for the changes produced by finite-term variations. In the actual matched shower evolution, the constraint imposed by matching to the LO 4-parton matrix elements fixes the finite terms, via the unitary procedure derived in [24], which was briefly recapped in section 3.1. The effective finite terms then depend on the full LO 4-parton matrix elements, and have a more complicated structure than the simple antenna functions we have so far been playing with. We shall therefore not attempt to integrate them analytically, but prefer instead to let vincia compute a numerical MC estimate for us.
Each point in that MC integration will involve computing at least one LO 4-parton matrix element, hence it is crucial to know how many points will be needed to obtain sufficient accuracy. Since everything else is handled analytically, this will be the deciding factor in determining the speed of the NLO-corrected algorithm. We shall perform a speed test below in section 5.4, but first we need to determine the accuracy we need on the δA integral. A first analytic estimate of the size of the δA terms can be obtained by simply computing the ones produced by switching from GGG to the vincia default antennae (summed and averaged over helicities [31]), with the following O(1) finite-term differences: with F Emit and F Split defined in eqs. (2.4) and (2.5). The δA terms produced by these differences are plotted in figure 13, for strong ordering in m D (left) and p ⊥ (center), and for smooth ordering in p ⊥ (right), respectively. As expected, they do come out to be numerically subleading, roughly of order α s /(2π), relative to LO (unity), yielding corrections ranging from a few permille to about a percent of the LO result. Finally, in figure 14, we include the full LO 4-parton matrix elements and plot the distribution of numerically computed δA terms during actual vincia runs, for 100,000 events. The result is now represented by a one-dimensional histogram, with δA on the x-axis and relative rate on the y-axis. On the left-hand pane, the δA distribution with default settings is shown on a linear scale, while the right-hand pane shows the same result on a logarithmic scale, including variations with higher numerical accuracy.
As mentioned above, the integration is done by a uniform Monte Carlo sampling of the δA integrands. We require a numerical precision better than 1% on the estimated size of the term (relative to LO) and, by default, always sample at least 100 MC points for each antenna integral. In the left-hand pane of figure 14, we see that, even with the full 4-parton LO matrix-element corrections included, the size of the δA terms remains below one percent for the vast majority of 3-parton phase-space points.
On the logarithmic scale in the right-hand pane of figure 14, however, it is evident that there is also a tail of quite rare phase-space points which are associated with larger δA corrections. Numerical investigations reveal that this tail is mainly generated by the integrals over the g → qq terms, in particular in phase-space points in which the gluon is collinear to one of the original quarks. This agrees with our expectation that these terms are the ones to which the pure shower gives the "worst" approximation, and hence they are the ones that receive the largest matrix-element corrections. As a test of the numerical stability of the NLO corrections for these points, we increased the minimum number of MC points used for the δA integration from the default 100 (shown with "+" symbols) to 400 ("×" symbols) and 1600 (" * " symbols), cutting the expected statistical MC error in half with each step, at the cost of increased event-generation time. Though we do observe a slight broadening of the distribution between the default and the higher-precision settings, the shifts should be interpreted horizontally and remain well under the required percent-level precision with respect to LO. The default settings are therefore kept at a minimum of 100 MC points, though we note that future investigations, in particular of more complicated partonic topologies, may require developing a better understanding of, and perhaps a better shower approximation to, these integrals, especially the g → qq contribution.
For completeness, we note that the runs used to obtain these distributions were performed using the new default "Nikhef" tune of vincia's NLO-corrected shower, which will be described in more detail in the following subsection. Parameters for the tune are given in appendix D.

LEP Results
Since we have restricted our attention to massless partons in this work, we shall mainly consider the light-flavour-tagged event-shape and fragmentation distributions produced by the L3 experiment at LEP for our validations and tuning, see [63]. We consider three possible vincia settings: • New default (NLO on): uses two-loop running for α s , with CMW rescaling of Λ QCD .
From the comparisons to event-shape variables presented in this section, we settled on a value of α s (M Z ) = 0.122. A few modifications to the string-fragmentation parameters were made, relative to the old default, to compensate for differences in the region close to the hadronization scale. The revised parameters are listed in appendix D, under the "Nikhef" tune.
• New default (NLO off). Identical to the previous bullet point, but with the NLO correction factor switched off. The three main event-shape variables that were used to determine the value of α s (M Z ) are shown in figure 15, with upper panes showing the distributions themselves (data and MC) and lower panes showing the ratios of MC/data, with one-and two-sigma uncertainties on the data shown by darker (green) and lighter (yellow) shaded bands, respectively. The Thrust (left) and C-parameter (middle) distributions both have perturbative expansions that start at O(α s ) and hence they are both explicitly sensitive to the corrections considered in this paper. The expansion of the D parameter (right) begins at O(α 2 s ). It is sensitive to the NLO 3-jet corrections mainly via unitarity, since all 4-jet events begin their lives as 3-jet events in our framework. It also represents an important cross-check on the value extracted from the other two variables.
For a pedagogical description of the variables, see [63]. Pencil-like 2-jet configurations are to the left (near zero) for all three observables. This region is particularly sensitive to nonperturbative hadronization corrections. More spherical events, with several hard perturbative emissions, are towards the right (near 0.5 for Thrust and 1.0 for C and D). The maximal τ = 1 − T for a 3-particle configuration is τ = 1/3 (corresponding to the Mercedes configuration), beyond which only 4-particle (and higher) states can contribute. This causes a noticeable change in slope in the distribution at that point, see the left pane of figure 15. The same thing happens for the C parameter at C = 3/4, in the middle pane of figure 15. The D parameter is sensitive to the smallest of the eigenvalues of the sphericity tensor, and is therefore zero for any purely planar event, causing it to be sensitive only to 4-and higher-particle configurations over its entire range. Finally, as an aid to constraining the Lund fragmentation-function parameters, the L3 study also included two infrared-sensitive observables: the charged-particle multiplicity and momentum distributions, to which we compare in figure 17, with the momentum fraction defined as There is again no noteworthy differences between the old and new default tunes.
Having determined the value of α s (M Z ) and the parameters of the non-perturbative fragmentation function, we extended the validations to include a set of jet-rate and jet-resolution measurements by the ALEPH experiment [64] (now without the benefit of light-flavour tagging), using the standard Durham k T algorithm for e + e − collisions [65], as implemented in the fastjet code [66]. We also compared to default pythia 8 and, for completeness, checked that the relative production fractions of various meson and baryon species were indeed unchanged relative to the old vincia default. Rather than presenting all of this information in the form of many additional plots, tab. 5 instead provides a condensed summary of all the validations we have carried out, via χ 2 values for each of the models with respect to each of the LEP distributions, including a flat 5% "theory uncertainty" on the MC numbers. Already from this simple set of χ 2 values, it is clear that the LO models/tunes are already doing very well 13 . This agreement, however, comes at the price of using a very large ("LO") value for α s , which is not guaranteed to be universally applicable.
The main point of the overview in tab. 5 is that an equally good agreement can be obtained with an α s (m Z ) value that is consistent with other NLO determinations [72], specifically α s (m Z ) = 0.122 , (5.4) once the NLO 3-jet corrections are included. This should carry over to other NLO-corrected processes, and hence the fragmentation parameters we have settled on should be applicable to future NLO-corrected studies with vincia, and can also serve as a starting point for NLO-level matching studies with pythia 8. In the latter context, the 2-loop running in particular could be retained, while the soft fragmentation parameters would presumably have to be somewhat readjusted to absorb differences between vincia and pythia 8 near the hadronization scale 14 .
13 Both vincia and pythia are known to give quite good fits to LEP data [24,31,67,68]. For comparisons including other generators and tunes, see mcplots.cern.ch [69]. 14 The differences in soft fragmentation parameters between existing LO vincia and pythia-8 tunes could be used as an initial guideline for such an effort, see, e.g., appendix D.  Table 5. χ 2 values for: Top: L3 light-flavour event shapes (left) and fragmentation variables [63], and LEP average meson and baryon fractions (right) [70,71]. Bottom: Durham k T n-jet rates, r nj , and jet resolutions, y ij , measured by the ALEPH experiment [64]. For the latter, the χ 2 calculation was restricted to the perturbative region, ln(y) > −8. A flat 5% theory uncertainty was included on the MC numbers. Both default pythia and the vincia (LO) tune use α s (m Z ) = 0.139 while the vincia (NLO) tune uses α s (m Z ) = 0.122.

Uncertainties
As in previous versions of vincia, we use the method proposed in [24] to compute a comprehensive set of uncertainty bands, which are provided in the form of a vector of alternative weights for each event. Each set is separately unitary, with average weight one 15 . The difference with respect to previous versions is that each variation now benefits fully from the inclusion of NLO corrections. When setting the parameter Vincia:uncertaintyBands = on, the uncertainty weights are accessible through the method double vincia.weight(int i=0); with i = 0 corresponding to the ordinary event sample, normally with all weights equal to unity, and the following variations, for i =: 1. Default: since the user may have chosen other settings than the default, the default is included as the first variation.
2. alphaS-Hi: all renormalization scales are decreased to µ = µ def /k µ , where µ def = p ⊥ for gluon emission and µ def = m qq for gluon splitting. The default size of the variation (k µ = 2) can be changed by the user, if desired. A second-order compensation for this variation is provided by the renormalization-scale sensitive term V 3µ .
3. alphaS-Lo: all renormalizzation scales are increased to µ = µ def * k µ , with similar comments as for alphaS-Hi above. 15 vincia currently does not attempt to give a separate estimate of the uncertainty on the total inclusive cross section. The uncertainties it computes only pertain to shapes of distributions and the effects of cuts on the total inclusive rate. 11. NLC-Lo: qg emission antennae use 2C F as color factor, with similar comments as above.
We emphasize that these variations are not all independent (for instance the α s and NLO variations are highly correlated) and hence the corresponding uncertainties should not be summed in quadrature. In the vinciaroot plotting tool included with vincia, the uncertainty band is constructed by taking the max and min of the variations. See the vincia HTML manual for more information about the uncertainty bands and [24] for details on their algorithmic construction.
To illustrate these variations, and the effect of the NLO 3-jet corrections upon them, we include the two plots shown in figure 18. We here take the Thrust observable as a representative example. (More such plots can be generated using the vinciaroot interface and the vincia03-root.cc example program included with the vincia code.) Similarly to previous plots in this paper, the top pane shows the normalized distribution, 1/σ dσ/d(1−T ), and the bottom one shows the ratio of theory to data. Now, however, there is also a middle pane, which gives the relative breakdown of the automated uncertainty variations into their respective components (normalized to unity). In each plot, we compare four individual runs of vincia to the automated uncertainty variations, with the latter based on the central run. This provides a useful cross check of whether the variations are indeed well represented by the automated estimates, before (left) and after (right) including the NLO 3-jet corrections.   In both plots, the scale-variation uncertainty dominates over the full range of the observable, highlighting that this is the main component that would need to be improved in order to obtain more precise results. Note, however, that both the central value at large 1 − T and the amount of scale variation, are improved by the introduction of the NLO 3-jet corrections in the right-hand pane. We also note that the distributions obtained from the explicit variation runs are faithfully reproduced by the automated variations, thus validating our confidence in the automated approach.

Speed
Although the CPU time required by matrix-element and shower/hadronization generators is still typically small in comparison to that of, say, full detector simulations, their speed and efficiency are still decisive for all generator-level studies, including tuning and validation, parameter scans, development work, phenomenology studies, comparisons to measurements corrected to the hadron level, and even studies interfaced to fast detector simulations. For this wide range of applications, the high-energy simulation itself constitutes the main part of the calculation. An important benchmark relevant to practical work is for instance whether the calculation can be performed easily on a single machine or not. Higher matched orders are characterized by increasing complexity and decreasing unweighting efficiencies, resulting in an extremely rapid growth in CPU requirements (see e.g. [31]). At NLO, the additional issues of negative weights and/or so-called counter-events can contribute further to the demands on computing power. With this in mind, high efficiencies and fast algorithmic structures were a primary concern in the development of the formalism for leading-order matrix-element corrections in vincia [24,30,31], and this emphasis carries through to the present work. We can make the following remarks.
• The only fixed-order phase-space generator is the Born-level one. All higher-multiplicity phase-space points are generated by (trial) showers off the lower-multiplicity ones. This essentially produces a very fast importance-sampling of phase-space that automatically reproduces the dominant QCD structures.
• Likewise, the only cross-section estimate that needs to be precomputed at initialization is the total inclusive one. Thus, initialization times remain at fractions of a second regardless of the matching order.
• The matrix-element corrected algorithm works just like an ordinary parton shower, with modified (corrected) splitting kernels. In particular, all produced events have the same weights, and no additional unweighting step is required.
• Since the corrections are performed multiplicatively, in the form of (1 + correction), with 1 being the LO answer, there are no negative-weight events and no counter-events.
The only exception would be if the correction becomes larger than the LO answer, and negative. This would correspond to a point with a divergent fixed-order expansion, in which case the use of NLO corrections would be pointless anyway. Moreover, as demonstrated by the plots in the previous sections, our definitions of the corrections are analytically stable (and numerically subleading with respect to LO) over all of phase space, including the soft and collinear regions, for reasonable renormalizationand evolution-variable choices.
• The parameter variations described in section 5.3 can be performed together with the matching corrections to provide a set of uncertainty bands in which each variation benefits from the full corrections up to the matched orders. These are provided in the form of a vector of alternative weights for each event [24], at a cost in CPU time which is only a fraction of that of a comparable number of independent runs.
These attributes, in combination with helicity dependence in the case of the leading-order formalism [31], allow vincia to run comfortably on a single machine even with full-fledged matching and uncertainty variations switched on. The inclusion of NLO corrections will necessarily slow down the calculation. The relative increase in running time relative to pythia 8, is given in tab. 6  of tree-level matching, with and without the NLO 3-jet correction 16 . Without it (but still including the default tree-level corrections which go up to Z → 5 partons), vincia is 5 times slower than pythia. With the NLO 3-jet correction switched on, this increases only slightly, to a factor 7. For a fully showered and hadronized calculation which includes second-order virtual and third-order tree-level corrections, we consider that to still be acceptably fast. Importantly, an event-generation time of a few milliseconds per event implies that serious studies can still be performed on an ordinary laptop computer.

Outlook and Conclusions
In this work, we have investigated the expansion of a Markov-chain QCD shower algorithm to second order in the strong coupling, for e + e − → 3 partons, and made systematic comparisons to matrix-element results obtained at the same order. Using these results, we have subjected the subleading properties of shower algorithms with different evolution/ordering variables and different renormalization-scale choices to a rigorous examination. At the analytical level, we have compared the logarithmic structures at the edge of phase space, and at the numerical level we have illustrated the difference between the expanded shower algorithm and the oneloop matrix element. We find that the choice of p ⊥ -ordering, with a renormalization scale proportional to p ⊥ yields the best agreement with the one-loop matrix element, over all of phase space. This elaborates on, and is consistent with, earlier findings [34,35]. Using the antenna invariant mass, m D , for the evolution variable still gives reasonable results in the hard regions of phase space, but leads to logarithmically divergent corrections for soft emissions, the exact form of which depends on the choice of renormalization variable. In the vincia code, we retain the option of using m D mainly as a way of providing a conservative uncertainty estimate.
With the NLO 3-jet corrections included as multiplicative corrections to the shower branching probabilities, we find that we can obtain good agreement with a large set of LEP event-shape, fragmentation, and jet-rate observables with a value of the strong coupling con-stant of α s (M Z ) = 0.122. This is in strong contrast with earlier (LO) tunes of both pythia and vincia which employed much larger values ∼ 0.14 to obtain agreement with the LEP measurements. The parameters for the NLO tune are collected in appendix D and represent the first dedicated NLO-corrected tune to LEP data. This paper is intended as a first step towards a systematic embedding of one-loop amplitudes within the vincia shower and matching formalism. To arrive at a full-fledged prescription, this will need to be extended to hadron collisions, ideally in a way that allows for convenient automation. A first step towards developing the underlying shower formalism for pp collisions was recently taken [73], and more work is in progress [74].
In addition, further studies should be undertaken of the impact of unordered sequences of radiation that can occur for the smooth-ordering case (it may be necessary to adopt a strategy similar to the truncated showers of the mc@nlo approach), and the mutually related issues of total normalization and how much of the (hard) corrections are exponentiated (similar to the differences between the powheg and mc@nlo formalisms, but here occurring at one additional order, where the relevant total normalization is the NNLO one). Finally, it would be interesting to develop an extension of this formalism that would allow second-ordercorrected antenna functions to be used at every stage in the shower, thereby upgrading the precision of the all-orders resummation, a project that would involve examining the secondorder corrections to branchings of qg and gg mother antennae as well. We look forward to following up on these issues in the near future.

A Infrared singular operators
Here we list the IR singularity operators from [22,33,59] as they are used in section 3. Because a detailed derivation of the calculation of Z → 3 jets can be found in [32] we restrict ourselves to listing the result in form that is convenient for our purpose. Divergences are regulated using dimensional regularization with d = 4 − 2 . Our results, before ultraviolet renormalization, are cross-checked with [32] where one must undo the renormalization in their case. In order to cancel the ultraviolet poles we need to renormalize the coupling according to (see also section 3.4) and S = (4π) exp(− γ E ) contains the factors characterizing the MS scheme. Due to the renormalization, the leading order calculation will generate a term quadratic in α s (µ 2 R ), which directly cancels the ultraviolet poles of the next-to-leading order calculation.

B.2 One-loop Matrix Element
The fixed-order expression relevant to matching in the vincia context is the one-loop matrix element normalized by the tree-level one. We decompose this into leading-and subleadingcolour pieces, as follows: 2 Re M with the LC piece containing the C A part of the gluon loops, the QL one containing the quark loops, proportional to n F T R , and the SL piece containing the subleading gluon-loop corrections, proportional to −1/N C . As usual in MC applications, we usually refer to "Leading Colour" as including both the N C and T R pieces. These are both associated with so-called planar colour flows that are simple to relate to the colour-flow representations used in Monte Carlo event generators, see e.g. [1,75]. The subleading-colour piece is included below for completeness, but has so far been left out of the NLO matching corrections implemented in the vincia code.
The notation of the infrared pole structure of these terms has been written similar to the integrated antenna in [22], with the difference that we have chosen to expand the scale factor µ in the integrated antenna terms in order to obtain explicitly dimensionless logarithms.
The quark has been labelled 1, the anti-quark 2 and the gluon 3. and the infrared singular operators, I (1) , given in appendix A.
With the one-loop matrix element expressed in this form, cancellation of the infrared poles against the integrated antennae (see below) coming from the shower will be particularly simple and will yield an expression purely dependent on the renormalization scale, µ R , and on the kinematic invariants s 12 , s 23 , and s 12 , but not on the scale factor µ.

C Antenna integrals
In this appendix we list the results of antenna integrals over phase space corresponding to the various evolution variables.

C.1.3 Energy ordering
The results for this evolution parameter are

C.2 Strong Ordering Gluon Splitting
The branching of a gluon splitting into a quark antiquark pair can only take place at the 3 → 4 level splitting. The generation of a gluon splitting takes place through an alternative form of phase space generation than the discussed m D , p ⊥ and E n variables. Instead phase space is sampled in a triangular surface comparable to m D ordering, yet in this case using only one cutoff, the Q 2 generated at the 2 → 3 level, to avoid the singular region of the gluon where K 1 = 1 , K 2 = −2 K 3 = 2, K 4 = −δ Ig − δ Kg , K 5 = 1. (C.37) We now turn to specific cases.

C.3.1 Smooth mass ordering
The only term from eq. (C.29) that requires specification is the damping factor = min(s 1 , s 2 ) min(s 1 , s 2 ) + min(s qg , s gq ) . (C.38) The computation of the individual antenna parts will require separating the phase space triangle in two regions (s 1 > s 2 and vice versa) in order to make the damping factor definite. After summing over these two regions we obtain the following values for gluon emission contributions L 1 = 2 ln (2)

C.4 Smooth Ordering Gluon Splitting
Additionally we also need to consider the gluon splitting antenna function for smooth ordering. Similar to the strong ordering case, the separate generation of gluon splitting variables allows for a new choice for evolution variable and thereby a different phase space surface. As in the case of gluon emission we allow for integration over the whole phase space, using the damping factor to limit the accessible area. A general notation is the following aq /qg (s 1 , s 2 ), (C. 50) with the definition for the gluon splitting antenna as in eq. (C.27) and the default choice of evolution variable for gluon splittings being Q 2 Ej = m 2 qq = s 2 .

D NLO Tune Parameters
In tab. 7 below, we list the perturbative and non-perturbative fragmentation parameters for the Nikhef NLO tune of vincia. For reference, we compare them to the current (LO) default Jeppsson 5 tune, which was used for comparisons to LO vincia in this paper.  Table 7. Parameters of the "Nikhef" NLO tune, compared to those of the "Jeppsson 5" LO tune.