Resummed predictions for jet-resolution scales in multijet production in $e^+e^-$ annihilation

We present for the first time resummed predictions at NLO + NLL' accuracy for the Durham jet-resolution scales $y_{n,n+1}$ in multijet production in $e^+e^-$ collisions. Results are obtained using an implementation of the well known CAESAR formalism within the SHERPA framework. For the 4-, 5- and 6-jet resolutions we discuss in particular the impact of subleading colour contributions and compare to matrix-element plus parton-shower predictions from SHERPA and VINCIA.


Introduction
Jet-production processes provide direct means to investigate the dynamics of the strong interaction. Especially multijet production poses a severe multi-scale problem complicating theoretical analyses. Besides integrated jet rates, differential distributions of jet-resolution scales give insight into the mechanism of jet production for a given jet-clustering algorithm. While for well-separated energetic jets a fixed-order QCD estimation might be adequate, in particular for small jet-resolutions this is bound to fail and all-orders resummation techniques need to be employed to provide reliable theoretical predictions. These account for soft and collinear QCD radiation that dominates the emission spectrum.
While jet rates are of phenomenological importance both at lepton and hadron colliders, we here focus on the resummation of jet-resolution scales in multijet production in e + e − annihilation. In particular, we consider the Durham jet algorithm [1], where the jet-resolution parameter is determined by Here, E i,j denote the energies of (pseudo) particles i and j, θ ij the angle between them, and Q 2 the squared centre-of-mass energy. This measure is used to successively cluster particles into jet objects. In what follows, we are interested in the resolution scales where an (n + 1)-jet final state is clustered in an n-jet final state, i.e., where the emission of an additional jet off the n-jet configuration gets resolved. More precisely, we consider the differential 3 → 4, 4 → 5 and 5 → 6 jet resolutions at next-to-leading logarithmic accuracy, i.e., resumming leading logarithms (LL) of type α k s L k+1 and next-to-leading logarithms (NLL) α k s L k with L ≡ − ln y n,n+1 appearing in the exponent of the observable's cummulant distribution. For many observables in e + e − annihilation, there exist results at NLL, NNLL or even N 3 LL order matched to exact NLO or NNLO QCD matrix elements, see for instance [2,3,4,5,6,7,8,9,10,11]. However, these are limited to two-or three-jet final states. Here, we are in particular exploring high-multiplicity processes, i.e., jet emission from 4-and 5-jet matrix elements that feature Born processes with non-trivial colour configurations. For the three-jet resolution y 23 , NLL resummed predictions were calculated in [12]; NNLL+NNLO results were presented in [13]. This observable in particular has been used in extractions of the strong coupling constant α s from LEP data, see for instance [14,15,16]. Recently, in Ref. [17], the extraction of α s from a simultaneous fit of the two-and three-jet rates has been presented. However, in this study the three-jet rate has been considered at fixed-order NNLO QCD without resummation of logarithmic enhanced terms. Up to now, for Durham three-and higher jet rates only the next-to-double-leading logarithms of order α k s L 2k−1 have been calculated [1]. In [18] this has been extended to the generalised class of k T -type algorithms defined in [19]. Besides for the extractions of α s , integrated and differential jet rates prove to be very useful for theoretical validation of parton-shower algorithms [18,20]. Corresponding experimental measurement data from the LEP experiments, see for example [21,22], form an integral ingredient to event generator tuning efforts, as they provide handles to constrain the parameters of phenomenological hadronisation models, see Ref. [23] for a review.
To achieve NLL resummation for jet resolutions, we rely on an independent implementation of the CAESAR formalism [24,25] in the SHERPA [26,27] event-generator framework, presented in [28]. A brief introduction to the CAESAR approach and details on the implementation of automated NLL resummation in the SHERPA framework are discussed in Sec. 2. This includes a brief presentation of the approach used to match our resummed predictions to exact QCD matrix elements. In Sec. 3, we give specific algorithmic details on the actual resummation of jet-resolution scales and discuss the validation of the soft function and the multipleemission contribution in particular.
In Sec. 4, we present our resummed predictions for jet production off 3−, 4− and 5−jet final states in e + e − annihilation. To this end we evaluate the jet resolutions y 34 , y 45 and y 56 at NLL accuracy, matched to exact 4−, 5− and 6−jet NLO QCD matrix elements, respectively. We give an account of subleading-colour contributions, by comparing our full-colour results to the strict and an improved large-N C approximation, the latter corresponding to what is typically used in QCD parton showers. Finally, we compare the NLO + NLL results to matrix-element-plus-parton-shower simulations from SHERPA and VINCIA, where we also address the impact of hadronisation corrections on the observables. We present our conclusions and an outlook in Sec. 5.

Semi-automated resummation within the SHERPA framework
In Ref. [28], an implementation of the CAESAR formalism in the form of a plugin to the SHERPA eventgenerator framework has been presented, which we employ for resummation of Durham jet-resolutions. A particular focus has there been put on validating the colour decomposition of hard-scattering matrix elements and the colour-insertion operators for multi-parton processes. The resummation plugin has since been used for NLL resummation of the thrust event-shape variable in e + e − , ep and pp collisions [28], with results matched to tree-level real-emission matrix elements using a point-wise local subtraction technique. Recently, the implementation has also been used to resum the soft-drop thrust shape in e + e − annihilation [29].
In this section, we briefly review the CAESAR formalism and comment on specific implementation aspects. Details on the application to resummation of jet-resolution scales will then be provided in Sec. 3.

The CAESAR formalism in a nutshell
The CAESAR formalism has originally been presented for final-state resummation in [24]; its extension to hadronic initial-states has been worked out in [30,31]. For an extensive review we refer to [25]. The method provides all necessary ingredients to perform NLL resummed calculations of recursively infrared and collinear safe global observables in a largely automated manner. The method is based on the observation that for the wide class of global event-shape observables, the resummed cumulative distribution for the observable value V ≤ v can, to NLL accuracy, be expressed in the simple form where the phase-space integral extends over the Born configurations B for each partonic channel δ. Dependence of the various contributions on those will implicitly be assumed in the following, although labels will be dropped. The jet-function H implements kinematic cuts on the Born events and ensures that only sufficiently hard configurations yield non-zero contributions. In the exponent, the collinear radiators for all hard legs l are summed. The function F accounts for the effect of multiple emissions, whereas the soft function S captures colour correlations. Its logarithmic dependence is defined via where the parameter λ is given by λ = α s (µ 2 R )β 0 L. In [25] the functions R l have been evaluated for the general class of observables V , where the impact of one additional, arbitrarily soft, emission with momentum k off leg l can be parameterised as i.e., in terms of the transverse momentum k T,l , rapidity η l , and azimuth φ l of the emission, measured with respect to leg l. Here, µ Q denotes the, in principle arbitrary, resummation "starting"-scale, to be distinguished from the centre-of-mass energy Q. Independence of the observable from the unphysical scale µ Q implies d l ∝ µ a Q . Explicit results for the radiators can be found in the appendix of [25].

Aspects of automated NLL resummation
The CAESAR formalisms is ideally suited for the automation of resummation of appropriate observables. With the observable parametrisation in terms of a, b l , d l and g l , the radiator functions are known. The colour decomposition of the Born process in a suitable basis and the corresponding soft-gluon correlators, needed for constructing the soft-function S, are independent of the actual observable to be resummed and can thus be pre-computed. The required Born-process partial amplitudes can be obtained from an automated matrix-element generator such as COMIX [32] for, in principle, arbitrary processes. However, the colourspace dimensionality quickly grows with the number of external partons [28], making calculations for highmultiplicity processes memory-intense. This motivates the construction of optimal, i.e., minimal-dimensional and orthogonal, bases, see Ref. [33,34,35]. An observable-dependent component that needs special attention is the multiple-emission function F, for that one often has to resort to numerical evaluations, which can, however, be pre-computed and tabulated. Here, we briefly introduce S and F in general terms, before discussing specific aspects and their validation for the resummation of jet-resolution scales in Sec. 3.

The soft function S
The automated generation of the soft-function S as defined in Eq. (3) for arbitrary multiplicity Born processes has been discussed in detail in Ref. [28]. Here, we only want to review some general aspects and briefly discuss the specific colour spaces we face for the case of multijet production in e + e − annihilation. The functional form of S is given by with |B the Born-process matrix element and Γ the soft anomalous-dimension matrix, cf. Eq. (10). By decomposing the Born matrix element over a colour basis 1 , and making a particular choice {|b α } for the colour structures of the Born process, we can define a metric and its inverse c αβ . This allows us to write with the hard matrix H αβ = A α A β in terms of partial amplitudes A γ . In this notation, the soft function may be written as with Here, the first sum runs over all colour dipoles ij, while the second one runs over those consisting only of final-final or inital-initial configurations, corresponding to the exchange of Coulomb gluons [36]. Explicit results for H, c, and Γ are known for up to four hard (coloured) legs; the form of S and the structure of Γ, however, holds more generally [37]. We obtain the hard matrix H from the matrix-element generator COMIX [32] that is part of SHERPA. All calculations of colour correlators relevant here contain at least one quark-antiquark pair and are performed in the trace basis. An all-orders trace basis is obtained by connecting each quark with an antiquark in all possible ways and subsequently attaching all gluons in all possible ways. Colour correlators are then evaluated by explicitly employing SU(N C ) identities. We note that this slightly differs from the approach in the COLORFULL [38] package, which uses specialised replacements in the trace basis instead.  One problem arising in the context of trace bases is that these are generally overcomplete for n qq +n g > N C , such that the metric defined in Eq. (7) does not possess a unique inverse. Although this has generally been solved in Ref. [28], within this study we can circumvent it in an even simpler way. The only critical process is e + e − → qq + 3g (cf. Tab. 1), for which the dimensionality of the basis exceeds the dimensionality of the colour space by one. Thus, we are able to define a basis with lower dimension by combining basis vectors corresponding to the same partial amplitude, i.e., We refer to this as the "reduced trace basis". In Sec. 3.2, we report on the validation of the colour correlators for the specific processes under investigation.

The multiple-emission function F
The component F of Eq. (2) captures the effect of multiple emissions on the cumulative distribution Σ.
Recursive infrared collinear safety of the observable guarantees that one can treat emissions below a cutoff as unresolved, such that the cancellation between real and virtual corrections is complete. To eliminate subleading terms, one would usually have to find suitable integral transforms to factorise the contributions of multiple individual emissions. This is for example straightforward for additive observables, where one can insert a simple integral representation of the corresponding Θ-function. In general, this procedure can however yield rather untractable expressions, or there are no such transformations known as is the case for the jet-resolution scales. In those cases, one can resort to numerical evaluations, rescaling momenta to eliminate contributions that vanish in the strict soft limit. In [25] a general form of the multiple-emission function F suitable for numerical evaluation was derived. It was in particular used to resum y 23 to NLL [12] and (with suitable additions) NNLL in [13] accuracy. The general formula for a specific flavour channel δ reads with where κ i (ζ iv ) defines an emission from the Born configuration B with an individual contribution to observable of V (κ i (ζ iv )) = ζ iv according to Eq. (4). Its rapidity with respect to leg l i is given by η i = ξ i ln 1/ζ iv /(a+b l ) such that ξ i denotes the fraction of the maximal rapidity. See also App. A for details on our notation. Note that the normalisation of the ξ integral is trivial if the observable does not scale with the rapidity in the l i collinear limit, b li = 0 ⇒ N li = 1. The channel δ defines the flavours of the Born partons and hence their corresponding Casimirs C l . In general, F depends on these Casimirs and L separately. For most obervables however, it turns out that the dependence is only on the combination R . We hence also use the notation F(R ) where no ambiguities can arise. If only the normalisation d l but not the actual scaling behaviour depend on B, the multiple-emission function can be evaluated for every Casimir permutation of a reference Born configuration B ref and does not have to be calculated on the fly. This is generally desireable for this approach to be useful, as avoiding subleading contributions in F requires the numerical evaluation of the limitsv → 0, → 0 to high accuracy, usually beyond the limits of double-precision. We have implemented Eq. (12) in a Monte Carlo code, independent of the original CAESAR work, and use this to calculate the F function in our framework. We evaluate the limits of the multiple-emission function for a reference Born configuration numerically, making use of multiple-precision arithmetic, and tabulate F for a grid of R values. This grid is interpolated using cubic Hermite splines, making additional use of the known monotonicity of F [39,40]. We have convinced ourselves that our code correctly reproduces F functions for simple additive observables with different parametrisations such as thrust, C-and D-parameter. As a test for non-additive variables, we reproduced known results for two-jet observables like heavy-hemisphere mass and thrust-major as well as the three-jet observable thrust-minor. We further validate the F-functions used here in Sec. 3.2.

Aspects of matching to fixed order
Finally, the matching to fixed-order QCD matrix elements needs to be accomplished in order to obtain reliable predictions outside logarithmically dominated phase-space regions, i.e., for hard, non-collinear emissions. In contrast to the original appproach in [28], we here aim for NLO + NLL accuracy and thus resort to matching the cumulant distribution rather than using a local subtraction-based method. We briefly discuss possible matching prescriptions widely used in the literature and introduce the logR scheme that we employ. To obtain physical predictions over the full range of the observable, the endpoint of the resummed prediction v max needs to be corrected to the actual real-emission kinematic endpoint. To this end we modify all all expressions by subleading contributions (see e.g. [2]). First, for the expansion to vanish at the endpoint v max , we need to modify the resummed result to In the limit v → v max this subtracts the expansion of S to first order in α s /2π, S (1) L, and the part of the expansion of R(v) that is linear in L. In the limit v → 0 this addition vanishes as a power correction, leaving the logarithmic accuracy of the expression unaffected. The remaining terms are forced to vanish at the endpoint v max by modifying all logartithms to such that L(v = v max ) = 0 and L → ln x v /v ∼ ln 1/v when v 1, and subtracting the corresponding change in the leading logarithm. The parameter x v represents the, in principle, arbitrary normalisation of the observable; p is an additional parameter that can be varied to estimate contributions from power corrections. There is a certain ambiguity in how to choose the default value for x v ≈ O(1). We follow the approach in [30] to cancel the d l dependence of the resummation formula. In the simplest form, and sufficient for the analysis presented here, this amounts to ln x v = 1/n l∈δ lnd l , where There are different approaches to write matched expressions at NLO + NLL for the cumulative distribution Σ, i.e., expressions that reproduce the fixed-order result including terms of order α 2 s relative to the respective Born process and reduce to the NLL result in the limit v → 0. At leading order, we might use a simple additive matching. To express more involved matching schemes, we introduce the following notation: Σ δ res,fo,match denotes the cumulative distribution in resummation, at fixed order, or matched between the two, for the (family of) channel(s) δ as defined by a suitable jet algorithm. In the following, labels are omitted in general expressions and we use the shorthand σ = Σ(1). We denote the expansion of any Σ to order α 2 s relative to the n-parton Born process as Practically, at least for Σ (2) fo , we only calculate We first define the multiplicative matching scheme, To order α (n−2) s α 2 s , we obviously recover the NLO cumulative distribution, apart from a missing additive constant σ δ, (2) fo . To the order of our calculations, it is, however, not affecting the normalised differential distributions we are interested in. In the limit v → 0, it reduces to With this scheme, terms of the order α (n−2) s α k s L 2k−2 are reproduced correctly in the expansion. Note that this relies on the fact that the leading-logarithmic terms ∝ α (n−2) s α k s L 2k do not depend on the kinematics or colour structure of the Born event but only on the flavour assignment. We are only interested in the terms where αs 2π C δ 1 multiplies one of those leading logarithms, all other cross terms are of subleading orders that are not computed consistently anyway. Thus we do not need to worry about the dependence of αs 2π C δ 1 on details of the Born event, apart from the channel δ. See also the discussion in [31]. We refer to this achieved accuracy in the expansion together with the presence of all terms of the form α k s L k+1 and α k s L k in ln(Σ res ) as NLL accuracy.
We also define a matching scheme based on the logarithm of the cumulative distributions, for consistency with the existing literature called LogR matching scheme, see e.g. [2,31], This defines our default matching choice in the evaluation of jet-resolution scales. The arguments on the achieved accuracy apply as for the multiplicative matching after expanding. The final distributions are then obtained by adding the different channels where the second sum accounts for channels not corresponding to any structure found in the physical Born process, but permitted by the jet-clustering algorithm. We will later on present normalised differential distributions given by

Resummation of Durham jet-resolution scales
In this section we want to specify the general CAESAR formalism for the case of Durham jet-resolution scales. In particular we discuss the actual observables and introduce our choices for parameters and scales. This is supplemented by the validation of the soft-colour correlators and the F functions we employ in the resummation of y 34 , y 45 and y 56 .

Definition of the observables
The functional form of the jet-resolution scales we consider here has been given in Eq. (1). We restrict ourselves to what is commonly referred to as the E-scheme, i.e., jets are combined by adding their fourmomenta. To completely define the observable y n,n+1 , we first require y n−1,n > y cut to define a hard nparton configuration around which the resummation is performed. From the point of view of an experimental measurement, it would be tempting to choose a rather small value of y cut to get a sizable cross section. However, besides the obvious need to regularise the infrared divergences in the Born event, cf. the jetfunction H in Eq. (2), the formalism we use here is based on the assumption that additional radiation is soft relative to all legs. In particular, configurations with y n,n+1 ≈ y n−1,n need to be described by the achieved fixed-order accuracy, and configurations where the production of one of the Born jets is already logarithmically enhanced need to be avoided. This prompts for larger values of y cut . In practice, for our theoretical studies, we find an observable definition with y cut = 0.02 a good compromise between these considerations, although admittedly smaller cuts would have to be investigated for a realistic experimental analysis, in particular for the highest multiplicities. Starting from this type of events, we are then interested in the next hardest resolution scale. If this (n+1)th jet consists of only one gluon that is soft and/or collinear to one of the legs l = 1 . . . n, one indeed obtains the general form of observables in the CAESAR formalism in Eq.
(2), with the parameters g l = 1, b l = 0 and a = 2, i.e., In this case, all parameters are independent of the leg l relative to which the emission is soft and/or collinear. Note that in Eq. (23), we have explicitly normalised by the squared centre-of-mass energy Q 2 . This corresponds to fixing the coefficients d l in Eq. (4) to d l = µ 2 Q /Q 2 . The explicit form of the radiator R l for an observable with the scaling behaviour given in Eq. (23), is specified in App. A using our conventions.

Calculational setup, parameter and scale choices
For 2-jet observables in e + e − annihilation, there is almost no ambiguity in the choice of scales, with the centre-of-mass energy Q = √ s basically being the only physical scale present in the Born process. For the multijet processes we consider here, this is no longer valid; due to phase-space constraints, some of the jets have to be associated with significantly lower scales. Our default choice for the resummation scale for the variable y n,n+1 is µ 2 Q = y n−1,n Q 2 . As d l is the same for all l, i.e., d l ≡ d, we have x v = d = y n−1,n . This means that the logarithms are effectively of the form ln y n−1,n /y n,n+1 . Our considered Born processes feature up to five jet thus constitute a severe multi-scale problem. We address this by choosing the renormalisation scale µ 2 R according to the CKKW prescription [41], i.e., based on the nodal Durham jet-resolutions of the n-parton Born process. For the observable y n,n+1 , this results in which is solved for µ 2 R assuming leading-order running of α s , thus with λ i = α s β 0 ln y i−1,i and Λ 2 QCD = Q 2 exp −1/α s β 0 . The strong coupling α s is evolved at two-loop accuracy with a fixed number of n f = 5 massless flavours.
We fix the endpoint of the resummed distribution to y max = min(y kin , y n−1,n ), where y kin denotes the maximal value of y n,n+1 with equal energies for all legs, E = Q/(n + 1). The second constraint enforces y n,n+1 < y n−1,n . Note, this condition is automatically satisfied for our central scale choices.
To identify the channels δ in the fixed-order calculations, we use the jet algorithm described in [42] for e + e − collisions. In our default setup, we restrict the jet algorithm to produce jets with at most one flavour, and additionally require that the flavour assignment is compatible with a Z decay, i.e., that there is at least one pair of jets identified as quark and antiquark jets with the same flavour. With these requirements, the second sum in Eq. (21) vanishes. As we treat all active quarks as massless, we can ignore the different flavours and only need to identify jets as either quark-or gluon-like. As there is also no dependence on the relative energy ordering of the legs, see in particular the discussion on the F function, we can collect all events with the same number of quarks and gluons into one family of channels δ.
The calculation of Σ res proceeds as described in the previous section. The fixed-order calculation is also done within the SHERPA framework, making use of the Catani-Seymour subtraction scheme [43] as implemented in COMIX [32]. We use OpenLoops [44] for one-loop virtual corrections to 3-, 4-, and 5-parton matrix elements. Virtual corrections to 6-parton matrix elements are generated with Recola [45,46,47]. Note that σ (2) is not needed in our matching formulas, and hence there is no need to compute any purely virtual corrections of order α (n−2) s α 2 s , while still achieving NLO accuracy in the normalised differential distribution.

Validation of soft-colour correlators
We validate all non-trivial colour correlators in Eq. (27) by comparing the eikonal approximation to exact (n + 1) tree-level matrix elements (see also Ref. [28]), in form of the ratio with the squared eikonal current We randomly pick 100 distinct, non-collinear (n + 1)-parton configurations, regularised by y cut = 0.02, and scale the momentum k s of the first gluon in the amplitude by a softness parameter λ s , i.e., k s → λ s k s , λ s → 0. This gluon is assumed to be emitted by the dipole spanned by its direct neighbours s − and s + which absorb its recoil by In the limit of soft, non-collinear kinematics, the ratio in Eq. (26) has to approach unity by the factorisation theorem of QCD. This is indeed observed for all the partonic channels relevant for the resummation of y 34 , y 45 and y 56 in e + e − annihilation. Fig. 1 contains the validation of all non-trivial colour contributions to partonic configurations relevant for 5-and 6-jet production, respectively.

Validation and results for the multiple-emission function
The F function is evaluated numerically for all relevant multiplicities and, where needed, different channels δ.
In the limitv → 0, all emissions become collinear to their hard legs and emissions are only clustered together if they originate from the same parton in the Born event. From this, it is straight-forward to see that the multiple-emission function is independent of details of the Born kinematics and can thus be evaluated for a reference configuration for each flavour channel δ. This amounts to calculate F for every possible number of external gluons. Applying the condition of emissions being only clustered if they were emitted from the same Born parton explicitly removes any dependence on the used kinematics. For reference, the numerical results for F(L) are shown in Fig. 2, for a configuration with µ 2 R = µ 2 Q = y cut Q 2 . This could be converted to F(R ) using Eq. (31). In the same figure, we also show how F(R ) approaches thev → 0 limit for several representative R values. It is not feasible to stick to double precision such that that multiple-precision arithmetics have to be used. However we are observing convergence starting from values ofv ≈ 10 −500 . It appears that, when interpreted as a function of R , the dependence of F on the Born channel is almost insignificant, in particular for y 56 . Note, however, that the two functions are not exactly equal in any case and that this similarity depends on the relatively close numerical values of the Casimirs C F and C A . To match to NLO fixed-order distributions, we also need the expansion of F. The argument of [12] applies to higher multiplicities, and we can generally write Numerically, we checked that these values are correctly reproduced in the corresponding integrals, again validating the limit.

Predictions for Durham jet-resolution scales
In this section, we present our results for the 4−, 5−, and 6−jet Durham resolution scales.

Resummed predictions for multijet resolutions
Our central results are the resummed distributions of the Durham jet resolutions y 34 , y 45 , and y 56 , in the definition discussed in Sec. 3.1. In particular, we require the Born events to possess n = 3, 4, 5 jets separated by at least y cut = 0.02 in the calculation for y n,n+1 , respectively. We perform our calculation at the LEP1 centre-of-mass energy √ s = 91.2 GeV. We use the LogR matching scheme as described in Sec. 2.3 to obtain physical distributions over the full range of the observable. Fig. 3 shows our predictions at LO and NLO, matched to NLL by the means of the preceding sections to achieve an accuracy of LO + NLL and NLO + NLL , respectively. To estimate theoretical uncertainties from missing higher-order corrections, we consider independent variations of µ R , µ Q , and p by factors of 0.5 and 2, resulting in the yellow (LO + NLL ) and orange (NLO + NLL ) bands. For higher jet multiplicities, we observe a narrowing of the distributions as well as a growing impact of NLO contributions compared to the LO result, being significant only for y 45 and y 56 . In the peak region, the impact of NLO corrections stays rather small in any case, being of order of a few percent only. Regarding y 34 , the NLO corrections do not have a large impact on the central prediction, we do, however, observe a considerable reduction of the theoretical uncertainties, at least away from the very soft region, where higherlogarithmic corrections become significant. This expected reduction of uncertainties is observed for the higher multiplicities as well. Here, the impact of NLO corrections on the matched distributions become larger when approaching the kinematic endpoint, y n,n+1 ≈ y kin , but consistently stay within the LO uncertainty band.
The remaining statistical uncertainty from the numerical determination of the F-function is propagated through to the final result, indicated by the error bars on the NLO + NLL prediction. It is entirely negligible compared to the overall uncertainty from scale variations.

Impact of fixed-order corrections
In Fig. 4, we compare the fixed-order predictions to the expansion obtained from the resummed results for the three observables y 34 , y 45 , and y 56 . We explicitly check their asymptotic agreement in the limit y n,n+1 → 0. As expected, the difference between fixed order and the expansion of the resummed distribution approaches zero in the soft limit. In the expansion Σ (2) , there are missing terms of order α (n−2) s α 2 s L 2 and α (n−2) s α 2 s L that are present in the fixed-order NLO calculation. Hence, we have a mismatch between fixed-order calculation and expansion which, in the differential distribution, is growing linearly in L, see the orange line in Fig. 4. After including the αs 2π C δ 1 coefficients, which effectively happens in the matching, only missing terms of order α (n−2) s α 2 s L remain. These are subleading with respect to the NLL resummation, although leading to a constant but finite difference in the differential distributions between fixed order and expansion. This is also demonstrated in Fig. 4, where in the purple line we include the matching coefficient as they would effectively appear in the multiplicative matching. As the NLO calcuation is computationally expensive, we do not attempt to carry it out at sufficient precision to allow for an extraction of the actual constant. Note that the matching schemes in Sec. 2.3 are designed such that the matched distributions behave in a physically meaningful way, i.e., are driven to zero by the resummed distribution, cf. Fig. 3, requiring a much less precise calculation in the soft tail of the fixed-order correction. We compared our results presented in the previous section to the ones obtained in the multiplicative matching scheme, but did not find significant differences, except for the region y n,n+1 > y cut . There, however, scale uncertainties also become very large, such that the difference between the two matching schemes is always covered by the NLO scale variation in the LogR matching scheme. We thus do not include an explicit uncertainty related to the matching scheme.

Impact of subleading colour contributions
Our results enable us to quantify the impact of subleading colour corrections on multijet observables. This is of particular interest for example in the context of recent and ongoing developments to systematically include such corrections into parton-shower event generators, see e.g. [48,49,50,51,52,53,54]. To assess the effect of subleading colour contributions, we redo our calculation in the t'Hooft large-N C limit, defined by taking N C → ∞ while keeping α s N C fixed [55]. This approximation, to which we will refer to as leading colour (LC), has a significant impact on the various contributions in Eq. (2). Firstly, all Casimirs corresponding to different legs are simplified to N C = 2C F = C A . Secondly, quark-loop contributions become negligible in both, the anomalous dimensions and the beta function. Finally, all non-planar diagrams vanish, leading to a simplification of colour-insertion operators. As the latter is only relevant to the colour-correlation contribution S, we define an "improved large-N C " approximation (imp. LC), in which we only treat the colour correlators appearing in Γ, and hence S, in the strict N C → ∞ limit, while still including the correct subleading contributions in R l and F. Those have a clear interpretation as sums over contributions from individual legs, such that the proper leg-specific SU(3) Casimir operators can be assigned. As in our main calculation, we match both resulting resummed distributions to the full-colour NLO result in the LogR scheme. The results of this study are presented in Fig. 5 and we collect details on the analytical treatment of subleading colour contributions in App. A.
Durham scale y 34 √ s = 91.2 GeV    Figure 4: Fixed-order predictions for y 34 , y 45 , and y 56 . We show the expansion of Σ res to the relevant orders, and the expansion to second order including the expression approaching the C 1 coefficient in the soft limit. The lower panels show the total difference between fixed order and expansion.
While there is a significant impact on the distributions when taking the t'Hooft large-N C limit, it is negligible for the improved large-N C limit, although in both cases, the results stay entirely within the uncertainty band of the NLO + NLL prediction. Moreover, in the peak region there is almost no difference between the full-colour and the improved large-N C result; the only visible difference being in the soft-tail region, where, however, the effect is below 5% for y 34 and y 45 and only sizeable in the ultra-soft region of y 56 .
An immediate interpretation of these observations in terms of subleading colour corrections in parton showers is not straight forward, as the matching involves an intricate interleave of subleading colour and kinematical corrections. It is, however, reasonable to argue that subleading colour effects are rather small and anyway overwhelmed by the effect of fixed-order corrections. At least for the observables considered here, there is virtually no difference between the calculations at full-colour and in the improved large-N C limit, as long as both are matched to exact, full-colour matrix elements.  Figure 5: Durham splitting scales y n,n+1 at NLO + NLL accuracy, in the N C → ∞ limit, and in the improved LC scheme.

Comparison to parton-shower predictions
Our resummed NLO + NLL predictions for the Durham jet resolutions provide highly non-trivial benchmarks for corresponding predictions from QCD Monte-Carlo event generators. They allow to gauge the results from parton-shower simulations and the methods used to combine these with exact higher-order QCD matrix elements and might guide the way to further improving showering algorithms. On the other hand, Monte-Carlo simulations provide means to include non-perturbative corrections from the parton-tohadron transition, that ultimately need to be taken into account for a realistic comparison with experimental data. As a first step in this direction we compare our resummed results against predictions from two distinct parton-shower implementations: the dipole parton shower CSSHOWER [56] as implemented in the SHERPA event generator and the VINCIA 1 antenna-shower plug-in [57] to the PYTHIA 8 event generator [58]. In Fig. 6 we present our results obtained using matrix-element plus parton-shower simulations at parton level for y 34 , y 45 , and y 56 in comparison to the NLO + NLL predictions. Furthermore, in the right panels of Fig. 6 we present hadron-level predictions, compared to the respective parton-level results.
To account for the kinematics and colour correlations of the respective Born processes, both showers are corrected to exact LO and NLO calculations, using two different strategies. One the one hand, in SHERPA we consider an MEPS@LO [59] calculation with exact tree-level matrix elements for e + e − → 2, 3, 4, 5 final-state partons. We also compare to a calculation using the MEPS@NLO merging strategy [60], including one-loop QCD corrections for e + e − → 2, 3, 4 partons via the MC@NLO method [61,62] as implemented in SHERPA and e + e − → 5 partons at tree level. We again make use of the matrix-element generator COMIX, and use OpenLoops to compute the virtual corrections. The merging parameter is set to y cut = Q cut /E CMS 2 = 10 −2 for the both the MEPS@LO and the MEPS@NLO calculation. The VINCIA antenna shower, on the other hand, is matched to e + e − → 2, 3 matrix elements at one-loop level as presented in [63] and e + e − → 4, 5, 6 tree-level matrix elements obtained from the MADGRAPH 4 matrix-element generator [64] via matrix-element corrections in the unitary GKS formalism [65]. Matrix-element corrections for 5− and 6− parton processes are smoothly regularised at a matching scale of Q match /E CMS = 0.05. For comparability, we restrict the phase space to strongly ordered branchings only. In both showers, the strong coupling is evolved at two-loop order in the CMW scheme, assuming an MS value of α s (m 2 Z ) = 0.118. We observe that the MEPS@LO and VINCIA samples agree well with each other and are close to the analytic calculation in the peak region as well as in the hard region, apart from the immediate neighbourhood of the endpoint. The MEPS@NLO sample is generally not improving the agreement between analytic and parton-shower prediction. In fact for y 34 and y 45 it yields somewhat larger deviations. However, we do not determine an explicit uncertainty estimate for the Monte Carlo predictions here, but their size can be considered similar to those of the analytic calculation. Accordingly, all the presented predictions are indeed in very good agreement. Comparing to our results in the previous section, the effects of subleading colour contributions are qualitatively different and smaller than the differences we observe between our resummed results and the parton-shower predictions. Apparently, for the particular observable definition we consider here, missing subleading colour contributions are not driving these differences but rather ambiguities related to recoil schemes, phase-space constraints and the treatment of subleading contributions in the running of α s , see for example Refs. [66,67,68].
As our observable definition differs from the usual experimental definition of jet-resolution scales by the hard cut we impose on the Born event, we need to gauge the influence of the parton-to-hadron transition on this observable. To this end we compare generator predictions on parton and hadron level. For SHERPA, we invoke SHERPA's default hadronisation model, based on cluster fragmentation [69] and furthermore, hadronise SHERPA's parton-level events with the Lund string fragmentation [70] as implemented in PYTHIA 6.4 [71]. Parton-level predictions of the VINCIA antenna shower are hadronised using the Lund string model in PYTHIA 8.2.
We obtain sizeable hadronisation corrections in all multiplicities. Their impact is very similar for the cluster-and string-fragmentation model applied to SHERPA's parton-level results. A qualitatively similar effect can be seen in the string hadronisation of VINCIA's parton-level results. In all cases, the soft tail of the distributions is significantly suppressed, leading to a narrowing of the distribution and hence a more pronounced peak. Reassuringly, the hard tail is only mildly affected by either of the fragmentation models. Quantitatively, however, there are differences of up to 10 − 20% remaining between the hadron-level predictions from SHERPA and VINCIA in the central peak region. We note that this is largely compatible with the deviations seen in the comparison at parton level. As the two predictions from SHERPA with different hadronisation models agree relatively well, it can be expected that generic non-perturbative uncertainties are rather moderate. Moreover, taking into account that we do not include uncertainties regarding the variation of hadronisation parameters, all hadron-level predictions agree reasonably well with each other.

Conclusions
For the first time, we have here obtained resummed predictions at NLO + NLL accuracy for multijet resolution scales in electron-positron annihilation. We employ the CAESAR formalism in a largely automated manner within the SHERPA event-generator framework. All relevant colour spaces were decomposed over the trace basis to obtain hard-scattering matrix elements and colour correlators that account for the insertion of soft-gluon radiation. Both, the construction of the basis and the calculation of colour insertions, have been performed in an automated way. Multijet matrix elements are obtained from the COMIX matrix-element generator, that is part of SHERPA. For NLO QCD predictions we obtain virtual corrections for e + e − → 3, 4, 5 partons from OpenLoops, while using Recola for virtual corrections to e + e − → 6 partons. For the evaluation of the multiple-emission contribution, represented by the F-function in the CAESAR formalism, we resort to a numerical evaluation using multiple-precision arithmetics.
We have derived predictions for the Durham jet-resolution scales y 34 , y 45 , and y 56 at NLO + NLL accuracy, using a Durham resolution y cut = 0.02 to restrict the respective Born configurations to sufficiently hard kinematics. The inclusion of NLO QCD corrections significantly reduces the theoretical uncertainties, estimated via scale variations.
We studied the impact of subleading colour contributions to our predictions, by repeating our calculations in the LC as well as the improved LC scheme, which we regard as being similar to the colour treatment in parton showers. We observe significant differences to the full-colour prediction only in the (strict) LC limit, while already at NLL the improved LC scheme well approximates the full-colour prediction. At NLO + NLL accuracy, we observe virtually no difference between the improved LC and the full-colour result.
As a benchmark for parton showers and an estimate of non-perturbative effects on the observable, we compared our resummed predictions against two distinctly different parton-shower algorithms, the VINCIA antenna shower plug-in to PYTHIA and the Catani-Seymour dipole shower in SHERPA, matched to LO and NLO. We observe good agreement between the VINCIA and MEPS@LO as well as the MEPS@NLO results. The found effects of subleading colour contributions are qualitatively different and smaller than the differences observed when comparing resummed results and parton-shower predictions. Apparently these will rather originate from ambiguities related to recoil schemes, phase-space constraints and the treatment of subleading contributions in the running of α s in the parton-shower implementations. The effect of nonperturbative corrections was studied by including hadronisation effects for the VINCIA and SHERPA parton showers, employing both, cluster and string fragmentation for the latter. We found good agreement between all hadron-level results. The observed deviations between SHERPA and VINCIA can already be found at the parton level. This confirms the suitability of jet-resolution scales for studies of perturbative QCD dynamics.
A direct comparison of our predictions to existing LEP measurements of jet resolutions is not straightforward, due to the required regularisation of the multijet Born processes. A corresponding reanalysis of LEP data would be desirable. To be able to compare to existing e + e − data, more general hierarchies of the multi-scale problem need to be addressed. By the generality of the approach presented here, our study may be extended to k T jet-resolution scales in hadronic collisions, motivating dedicated measurements at the LHC.

A. Analytic expressions
For completeness, we here collect all analytic expressions of the CAESAR formalism that did not appear in the main text, restricted to the jet-resolution scales studied here and our conventions. In particular, we choose to rescale the arguments of all logs with x v = d(µ Q ), so that no explicit dependence on d appears. We use the short-hand λ = α s (µ 2 R )β 0 L, and fix the normalisation of colour operators to T R = 1/2, resulting in C l = C A = N C for gluons and C l = C F = (N C − 1/N C )/2 for quarks. The radiator of leg l is given by Its derivative with respect to L is single logarithmic, We also use the shorthands R = l R l , R = l R l where this eases the notation. The usual constants in the MS scheme are where we have already employed our normalisation convention. In the t'Hooft limit, N C → ∞ with α s N C = α s,0 = const., all expressions depend only on the finite quantities , α s β 0 → 11 12π α s,0 , When working in the large-N C limit, the F function needs to be re-evaluated in principle. However, it only depends on the ratios of sums of Casimirs for the different legs and on R . Thus only the configurations with mixed quark and gluon content need to be re-computed with C A /C F = 2. When evaluating colour correlators in the large-N C limit, we are interested in Note that we need to normalise the basis vectors in order to obtain a finite result. As expected, the b α | T i T j b β LNC correlators vanish for all non-planar diagrams and give finite contributions otherwise. Practically, we calculate the large-N C colour correlators by comparing the powers of N C in the numerator and denominator of Eq. (34) and keeping only those with the same highest power. This modification of the colour correlators is the only one that is still present in the improved LC scheme.