Jet production in the CoLoRFulNNLO method: event shapes in electron-positron collisions

We present the CoLoRFulNNLO method to compute higher order radiative corrections to jet cross sections in perturbative QCD. We apply our method to the computation of event shape observables in electron-positron collisions at NNLO accuracy and validate our code by comparing our predictions to previous results in the literature. We also calculate for the first time jet cone energy fraction at NNLO.


Introduction
The strong coupling α s is one of the most important parameters of the standard model. A clean environment for determining α s is the study of event shape distributions in electron-positron collisions. Particularly well suited for this task are quantities related to three-jet events, as the leading term in a perturbative description of such observables is already proportional to the strong coupling. Accordingly, three-jet event shapes were measured extensively, especially at LEP [1][2][3][4]. The precision of experimental measurements calls for an equally precise theoretical description of these quantities. Because the strong interactions occur only in the final state, non-perturbative QCD corrections are restricted to hadronization and power corrections. These corrections can be determined either by extracting them from data by comparison to Monte Carlo predictions or by using analytic models. Hence, the precision of the theoretical predictions is mostly limited by the truncation of the perturbative expansion in the strong coupling.
Current state-of-the art includes next-to-next-to-leading order (NNLO) predictions for the three-jet event shapes of thrust, heavy jet mass, total and wide jet broadening, C-parameter and the two-to-three jet transition variable y 23 [5,6], as well as for oblateness and energyenergy correlation [7]. Next-to-leading order (NLO) predictions for the production of up to five jets [8][9][10][11][12] (and up to seven jets in the leading color approximation [13]) are also known. Moreover, logarithmically enhanced contributions to event shapes can be resummed at up to next-to-next-to-leading logarithmic (NNLL) accuracy [14][15][16][17][18] and even at next-to-next-to-nextto-leading logarithmic (N 3 LL) accuracy for some observables [19,20].
In addition to its phenomenological relevance, three-jet production in electron-positron collisions is also an ideal testing ground for developing general tools and techniques for higher-order calculations in QCD. The straightforward evaluation of radiative corrections in QCD is hampered by the presence of infrared singularities in intermediate stages of the calculation which cancel in the final physical results for these observables. Nevertheless, they must be regularized and their cancellation has to be made explicit before any numerical computation can be performed. This turns out to be rather involved for fully differential cross sections at NNLO and constructing a method to regularize infrared divergences has been an ongoing task for many years .
In this paper we present a general subtraction scheme to compute fully differential predictions at NNLO accuracy, called CoLoRFulNNLO (Completely Local subtRactions for Fully differential predictions at NNLO accuracy) [41][42][43][44][45][46][47][48][49][50][51]. The method uses the known universal factorization properties of QCD matrix elements in soft and collinear limits [54][55][56][57][58][59][60][61][62] to construct completely local subtraction terms which regularize infrared singularities associated with unresolved real emission. Virtual contributions are rendered finite by adding back the subtractions after integration and summation over the phase space and quantum numbers (color and flavor) of the unresolved emission. We have worked out the method completely for processes with a colorless initial state and involving any number of colored massless particles in the final state. We validate our method and code by computing NNLO corrections to three-jet event shape variables and comparing our predictions to those available in the literature [5,6]. We also present here for the first time the computation of the jet cone energy fraction (JCEF) at NNLO 1 accuracy. We note that the CoLoRFulNNLO method has already been successfully applied to compute NNLO corrections to differential distributions describing the decay of a Higgs boson into a pair of b-quarks [63], as well as to the computation of oblateness and energy-energy correlation in e + e − → 3 jet production [7].
The paper is structured as follows: after introducing our notation and conventions in section 2, we present the CoLoRFulNNLO method in section 3. In section 4 we describe the application of the general framework to the specific case of three-jet production. In particular, we show that the double virtual contribution is free of singularities. Our predictions for event shape observables follow in section 5. We draw our conclusions and give our outlook in section 6.

Notation and conventions 2.1 Phase space and kinematics
The phase space measure in d = 4 − 2 dimensions for a total incoming momentum Q µ and n massless outgoing particles reads dφ n (Q 2 ) ≡ dφ n p µ 1 , . . . , p µ n ; Throughout the paper, we will use y ik to denote twice the dot-product of two momenta, scaled by the total momentum squared Q 2 . For example, We also introduce the combination Y ik,Q = y ik y iQ y kQ (2.3) for later convenience.

Matrix elements
We use the color and spin space notation of ref. [64] where the renormalized matrix element for a given process with n particles in the final state, |M n , is a vector in color and spin space, normalized such that the squared matrix element summed over colors and spins is given by The renormalized matrix element has the following formal loop expansion where the superscript denotes the number of loops. We will always consider matrix elements computed in conventional dimensional regularization (CDR) with MS subtraction. We introduce the following notation to indicate color-correlated squared matrix elements (obtained by the insertion of color charge operators between M ( 1 ) n | and |M ( 2 ) n ): (2.6) The color charge algebra for the product a (T i ) a (T k ) a ≡ T i ·T k is where C i is the quadratic Casimir operator in the representation of particle i. We use the customary normalization T R = 1/2, and so C A = 2T R N c = N c in the adjoint and in the fundamental representation.

Ultraviolet renormalization
In massless QCD renormalized amplitudes |M n are obtained from the corresponding unrenormalized amplitudes |A n by replacing the bare coupling α B s with the dimensionless renormalized coupling α s ≡ α s (µ) computed in the MS scheme and evaluated at the renormalization scale µ, where and S MS = (4π) exp(− γ E ) corresponds to MS subtraction, with γ E = −Γ (1) the Euler-Mascheroni constant. Although the factor (4π) exp(− γ E ) is often abbreviated as S in the literature, we reserve the latter to denote which emerges in the integration of the angular part of the phase space in d = 4−2 dimensions. If the loop expansion of the unrenormalized amplitude is written as where m is the number of massless final-state partons in the Born process) then using the substitution in eq. (2.8) the relations between the renormalized and the unrenormalized amplitudes are given as follows: and where The role of the factors of (µ 2 /µ 2 0 ) is to change the regularization scale to the renormalization scale so that the renormalized amplitudes in eqs. (2.12)-(2.14) only depend on µ. Furthermore, after the IR poles are canceled in a fixed order computation, we may set = 0, therefore, the factors of (µ 2 /µ 2 0 ) and S MS in C(µ, µ 0 , q; ) do not give any contribution, so we may perform the replacement (2.16)

Jet production in CoLoRFulNNLO
We consider the production of m jets from a colorless initial state as in, e.g., Higgs boson decay or electron-positron annihilation into hadrons. In perturbative QCD the cross section for this process is given by an expansion in powers of the strong coupling α s . At NNLO accuracy we retain the first three terms in this expansion The leading order contribution is simply given by the integral of the fully differential Born cross section dσ B m of m final-state partons over the the available m-parton phase space defined by the observable J, (often called jet function) Here and in the following J m denotes the value of the infrared-safe observable J evaluated on a final state with m partons.

The NLO correction
The NLO correction is a sum of the real radiation and one-loop virtual terms, both divergent in four dimensions. These two contributions can be made finite simultaneously by subtracting and adding back a suitably defined approximate cross section dσ Several prescriptions are available for the explicit construction of the approximate cross section [42,47,[64][65][66]. Specifically in the CoLoRFulNNLO framework, it is written as where the approximate matrix element for processes with m + 1 partons in the final state is given by [42,43], On the right-hand side of eq. (3.6), C (0,0) ir and S (0,0) r denote counterterms which regularize the p i || p r collinear limit and the p µ r → 0 soft limit in arbitrary dimensions. The role of the C ir S (0,0) r soft-collinear counterterm is to make sure that no double subtraction takes place in the overlapping soft-collinear phase space region. These counterterms were all defined explicitly in refs. [42,43]. In our convention the indices of C ri . As the sums in eq. (3.6) over i and r are likewise not ordered, the factor of 1 2 is needed so that we count each collinear limit precisely once. Finally, the meaning of the superscript ( 1 , 2 ) is the following: The corresponding counterterm involves the product of an 1 -loop unresolved kernel (an Altarelli-Parisi splitting function or a soft eikonal current) with an 2 -loop squared matrix element (in color or spin space). Specifically, (0, 0) means that in these subtraction terms a tree level collinear or soft function acts on a tree level reduced matrix element. Such superscripts will appear also for other counterterms throughout the paper. Importantly, the approximate matrix element in eq. (3.6) takes into account all color and spin correlations in infrared limits, and hence it is a completely local and fully differential regulator of the real emission matrix element over the m + 1-particle phase space. The complete locality of the subtraction is a necessary condition for the regularized real contribution, to be well-defined in four dimensions. As pointed out long ago [64], when the subtraction terms are not fully local, for instance because spin correlations in gluon decay are neglected, the evaluation of the difference m+1 dσ R m+1 J m+1 − dσ . Rather, the definition must be supplemented by the precise specification of an integration procedure which must be shown to give the correct numerical values for all integrals that are finite away from d = 4, but whose four-dimensional value is ill-defined. As in CoLoRFulNNLO the subtractions are completely local, eq. (3.7) is well-defined in four dimensions as it is and may be computed with whatever numerical procedure is most convenient. These remarks apply also to the regularized double real and real-virtual cross sections in eqs. (3.14) and (3.15) which enter the NNLO correction.
Turning to the virtual contribution, the Kinoshita-Lee-Nauenberg (KLN) theorem ensures that the integral of the approximate cross section precisely cancels the divergences of the virtual piece for infrared-safe observables, so adding back what we have subtracted from the real correction, the virtual contribution becomes finite as well. We have performed the integration of the various subtraction terms analytically in ref. [42] and here we only quote the result, which can be written as, where the ⊗ product is defined in eq. (2.6) and the insertion operator is in general given by [42] 1 The variables Q µ , y iQ and Y ik,Q were defined in section 2.1. The kinematic functions C 1,i (y iQ ; ) and S (0),(i,k) 1 (Y ik,Q ; ) have been computed as Laurent expansions in in ref. [42]. They are needed up to finite terms in a computation at NLO accuracy and up to O( 2 ) in a computation at NNLO accuracy. We present these kinematic functions explicitly up to O( ) in appendix A, which is sufficient for checking the cancellation of the -poles at NNLO analytically. We note that there is no one-to-one correspondence between the unintegrated subtraction terms in eq. (3.6) and the kinematic functions that appear in eq. (3.9). The latter are obtained from the former by integrating over the unresolved momentum as well as summing over all unobserved quantum numbers (color and flavor), and organizing the result in color and flavor space. Loosely speaking, the integrated form of C . We are, however, free to assign the integrated form of C ir S (0) r to either of the integrated counterterms. This final organization was performed differently in ref. [42] and in this paper. In ref. [42] we grouped the integrated form of C ir S , while here we find it more convenient to group it into C It is straightforward to check that the poles of this expression coincide with those of the I({p}; ) operator of ref. [64], hence 1 dσ R,A 1 m+1 as given in eq. (3.8) correctly cancels all -poles of the virtual cross section dσ V m . Thus the regularized virtual contribution, is finite and integrable in four dimensions.

The NNLO correction
The NNLO correction to the cross section is a sum of three contributions, the tree level double real radiation, the one-loop plus a single radiation and the two-loop double virtual terms, which are all divergent in four dimensions. In the CoLoRFulNNLO method, we render these terms finite by the rearrangement where, (3.16) The right-hand sides of eqs. (3.14) and (3.15) are integrable in four dimensions by construction [41,43,44], while the integrability in four dimensions of eq. (3.16) is ensured by the KLN theorem.
Equation (3.14) includes the double real (RR) contribution that is singular whenever one or two partons become unresolved. In order to regularize the two-parton singularities, we subtract an approximate cross section, ir;s and S (0,0) rs in eq. (3.18) are subtraction terms which regularize the p i || p r || p s triple collinear, the p i || p r , p j || p s double collinear, the p i || p r , p µ s → 0 one collinear, one soft (collinear+soft) and the p µ r → 0, p µ s → 0 double soft limits. The rest of the counterterms appearing in eq. (3.18) account for the double or triple overlap of limits, hence multiple subtractions are avoided in overlapping double unresolved regions. The role of each specific counterterm is suggested by the notation. For instance, C irs CS (0,0) ir;s accounts for the triple collinear limit of the collinear+soft counterterm, with the rest of the counterterms having similar interpretations. All functions appearing in eq. (3.18) were defined explicitly in ref. [43]. The factors of 1 6 , 1 8 , etc., in eq. (3.18) appear so that each limit is counted precisely once, since the collinear indices of counterterms and the sums over them are not ordered in our convention.
After subtracting the double unresolved approximate cross section, the difference is still singular in the single unresolved regions of phase space. To regularize it, we also subtract where A 1 has been defined in eq. (3.6). To avoid double subtraction in overlapping single and double unresolved regions of phase space, we must also consider dσ RR,A 12 m+2 where the iterated single unresolved approximate matrix element reads with the three terms above given by [43], (3.25) The notation in eqs.  [43]. Hence the overlap of single and double unresolved subtractions is properly taken into account. All of the counterterms appearing in eqs. (3.23)-(3.25) were defined in ref. [43] explicitly. As before, the collinear indices and sums over them in eqs. (3.22)-(3.25) are not ordered, so factors of 1 2 appear at various instances. The combination of terms appearing in eq. (3.14) was shown to be integrable in all kinematic limits in ref. [43]. Thus, the regularized double real contribution to the m-jet cross section is finite and can be computed numerically in four dimensions for any infrared-safe observable.
Turning to eq. (3.15), it describes the emission at one loop of one additional parton, the real-virtual (RV) contribution. In addition to explicit -poles coming from the one-loop matrix element, the RV contribution has kinematical singularities when the additional parton becomes unresolved. The explicit poles are cancelled by the integral of the single unresolved subtraction term in the double real emission contribution to the full NNLO cross section, which is simply given by eqs. (3.8) and (3.9) after the obvious replacement of m → m + 1 As shown above in eq. (3.10) the combination, is finite in . Nevertheless, eq. (3.27) is still singular in the single unresolved regions of phase space and requires regularization. We achieve this by subtracting two suitably defined approximate cross sections, dσ which matches the kinematic singularity structure of dσ RV m+1 . The general definition of the real-virtual counterterm is [44] The basic organization of this subtraction in terms of unresolved limits is identical to the tree level single unresolved counterterm in eq. (3.6). However in eq. (3.29) we have terms with tree level collinear or soft functions multiplying (in color or spin space) one-loop matrix elements (those with the (0, 1) superscript), as well as terms with one-loop collinear or soft functions multiplying tree level matrix elements (denoted with the (1, 0) superscript). This reflects the structure of infrared factorization of one-loop QCD matrix elements [59][60][61][62]. The functions appearing in eq. (3.29) are defined explicitly in ref. [44].
Then we consider the counterterm, which regularizes the kinematic singularities of 1 dσ . This counterterm is given by [44] The structure of this subtraction in terms of unresolved limits is again the same as the tree level single unresolved counterterm in eq. (3.6). However we have two types of terms for each limit, labeled by the different superscripts. The reason is the following. This counterterm is constructed from factorization formulae describing the behavior of the product of a QCD squared matrix element times the I (0) 1 insertion operator of eq. (3.9) in the collinear and soft limits. (The existence of a universal collinear factorization formula for the product |M 1 is not guaranteed by the factorization properties of QCD matrix elements. The requirement that such a formula exists puts highly non-trivial constraints on the form of I (0) 1 , i.e., on the specific definition of the single unresolved approximate cross section. See section 4.1.1 of ref. [44] for a discussion of this point.) These factorization formulae were computed in ref. [44] and turn out to be sums of two pieces. Both pieces involve the product of a tree level collinear or soft function times a tree level matrix element. One piece is further multiplied by the I (0) 1 insertion operator appropriate to the reduced matrix element, while the other is multiplied with a well-defined scalar (in color space) remainder function R. The superscripts on the various terms in eq. (3.31) are meant to reflect this structure. The combination of terms appearing in eq. (3.15) is both free of -poles and integrable in all kinematically singular limits [44]. Hence, the regularized real-virtual contribution to the m-jet cross section is finite and can be computed numerically in four dimensions for any infrared-safe observable.
Finally, the two-loop double virtual (VV) contribution to the NNLO corrections appears in eq. (3.16). The VV contribution has explicit infrared poles that cancel against the poles of the four integrated counterterms, which are shown in eq. (3.16). The integral of the real-virtual counterterms (the last two terms of eq. (3.16)) was computed in refs. [45,46,48] and can be written as and (3.35) Again, there is no one-to-one correspondence between the unintegrated double unresolved subtraction terms in eqs. (3.29) and (3.31) and the kinematic functions that appear in eqs. (3.34) and (3.35). The latter are obtained from the former after integration over unresolved momenta and summation over unobserved colors and flavors. This remark applies also to the rest of the insertion operators to be discussed below.
The integral of the iterated single unresolved counterterm (the third term of eq. (3.16)) was evaluated in ref. [49] yielding the result The insertion operator has five contributions according to the possible color structures, Finally, the integration of the collinear-type contributions to the double unresolved counterterm (the second term of eq. (3.16)) was performed in ref. [50]. The soft-type contributions to the same integral were presented in ref. [51]. We find where the structure of the insertion operator I  The resulting expressions are rather lengthy and involve, in addition to logarithms, dilogarithms and trilogarithms of rational arguments in the variables y iQ and Y jk,Q . For the finite parts, we computed analytically all terms that diverge logarithmically on the boundaries of the phase space (i.e., when y iQ → 0 and/or Y jk,Q → 0), while the remaining regular contributions were computed numerically. We stress that our method is generic and we can construct counterterms for processes with an arbitrary number m of jets in the final state. The only missing ingredients are the corresponding two-loop matrix elements, and current only the two-loop matrix elements for two and three-jet production are available. Since the poles of all integrated counterterms are known analytically, we can demonstrate explicitly that the regularized double virtual contribution to the m-jet cross section is finite and free of -poles. For m = 2 this was done in ref. [63], while the m = 3 case will be discussed in the next section.

Electron-positron annihilation into three jets
We consider e + e − → 3 jet production through the exchange of a photon or a Z boson of momentum Q in the s channel. Through NNLO in QCD, this production cross section receives contributions from the following partonic subprocesses: where we show the four-momenta of the particles in parentheses. The tree level matrix elements for the production of five jets were first obtained in refs. [67][68][69], while the one-loop corrections to four jet production have been computed in refs. [59,[70][71][72]. The two-loop matrix elements for γ * /Z * → qqg are also available both in squared matrix element form [73] and as helicity amplitudes [74]. In the CoLoRFulNNLO framework the subtraction terms correctly account for all spin and color correlations in the various infrared limits. Hence, we also need the threeparton and four-parton matrix elements including color and/or spin correlations. When there are only three partons in the final state, the color correlations factorize completely (see eq. (4.6) below), so computing the color-correlated three-parton matrix elements is trivial at any loop order. This is no longer the case for the four-parton matrix elements. In our computation, we only need the four-parton color-correlated matrix elements at tree-level and these are given in ref. [75] 2 . The required spin-correlated matrix elements on the other hand are rather easy to implement starting from helicity amplitudes.
The sum of the three-, four-and five-parton contributions is finite for any infrared-safe observable, but the four-and five-parton contributions to three-jet observables contain infrared singularities associated with unresolved real emission, which must be subtracted and cancelled against the infrared singularities coming from loop integrals in the three-and four-parton final states. We accomplish this cancellation with the CoLoRFulNNLO method as outlined in the previous section.

e + e − → 3 jet production at NLO
It is instructive to spell out the computation of the NLO correction in some detail. The fourparton real emission contribution to the differential cross section for three-jet production is (4.1) The integral over the phase space is divergent in four dimensions because of the singularities in the regions where one parton is collinear and/or soft. In order to regularize those singularities, we subtract where the approximate matrix elements are defined in eq. (3.6). The counterterms are explicitly defined in refs. [42,43] in a form that is immediately suitable for inclusion in a general purpose computer code. By generating sequences of phase space points tending to each infrared limit, we have checked that the sum of subtractions correctly reproduces the real emission differential cross section point-by-point in any single unresolved region of phase space. As a consequence the difference dσ NLO is integrable in four dimensions and the regularized real contribution can be computed using whatever numerical procedure is most convenient.
Turning to the three-parton virtual contribution, we have qqg . The insertion operator I (0) 1 is given in eq. (3.9). For three-jet production, as there are only three partons in the final state, the color connections that appear in the generic case in eq. (3.9) factorize completely, Thus, (4.7) Using eq. (3.10) (or the expressions in appendix A), it is straightforward to check that and that the combination is free of -poles. Thus eq. (4.9) is finite in four dimensions and the regularized virtual contribution can be computed using standard numerical techniques for any infrared-safe observable.

e + e − → 3 jet production at NNLO
Turning to the NNLO correction, we consider first the double real emission contribution to the differential cross section for three-jet production, (4.10) The integral over the phase space is divergent in four dimensions because of infrared singularities in regions of phase space where one or two partons are collinear and/or soft. In order to regularize those singularities, we subtract approximate cross sections dσ as explained in section 3.2. The counterterms are defined in ref. [43] explicitly, in a form directly suited for implementation into a general purpose computer code. We have checked in all kinematic limits that the difference is integrable in four dimensions by generating sequences of phase space points tending to each infrared limit. Hence the double real emission differential cross section is regularized point-bypoint in phase space. The complete locality of the subtractions then ensures that the integral of eq. (4.11) is well-defined and finite in four dimensions for any infrared-safe observable and can be computed with any suitable numerical technique.
The real-virtual contribution to the differential cross section is . The calculation in ref. [42] for general m assures us that the combination is finite in . (Of course, this can be checked explicitly using eq. (3.10) or the expressions in appendix A as well.) However, eq. (4.13) is still singular in the single unresolved regions of phase space. We regularize these singularities by subtracting the approximate cross sections dσ is both free of -poles and integrable in all kinematically singular limits in four dimensions (as usual, by generating sequences of phase space points tending to all infrared limits). Thus, since the subtractions are fully local, the regularized real-virtual contribution to the 3-jet differential cross section is well-defined and finite and can be computed numerically in four dimensions for any infrared-safe observable.
Finally, the double virtual contribution to the differential cross section reads and contains explicit -poles coming from the two-loop matrix element and the square of the one-loop matrix element. The structure of these poles was presented explicitly in ref. [73] which we reproduce here using our conventions for the notation: where the universal constants are and the three-parton insertion operator is qqg (s 12 , s 13 , s 23 , µ 2 ; ) = α s 4π with γ q = 3 2 C F and γ g = β 0 2 . (4.21) The signs of the imaginary parts of the (−s ik ) − factors are fixed by the usual s ik +iε prescription on the Feynman-propagators, Hermitian conjugation flips the sign of the imaginary parts. We note that the poles of this operator are closely related to those of the I    form a remarkably simple expression: It is easy to convince oneself that only the universal pole parts of the I (0) 1 operator (given in eq. (3.10) for general m) enter the computation of the poles of J 2 . Furthermore, looking at the explicit definition of I (0) 1 in eq. (3.9), we see that the J 2 operator in eq. (4.27) can be written by simply counting the radiating partons in the event (two quarks and one gluon in our example). This additive nature of J 2 , which is also valid for two-jet production, hints that in general is free of -poles, although to perform the algebra for the 1/ 2 and 1/ poles still requires some effort. Hence eq. (4.29) is finite in four dimensions and we can compute the regularized double virtual differential cross section for any infrared-safe observable numerically.

Event shapes old and new
The CoLoRFulNNLO method provides a robust subtraction scheme for computing NNLO corrections to processes with a colorless initial state (for the moment) and any number of final state jets, provided all necessary matrix elements are known. We have implemented the method in a general purpose, automated parton-level Monte Carlo code which can be used to compute any infrared-safe observable at NNLO accuracy in e + e − → 3 jets. To demonstrate the validity of our code, we compute NNLO corrections to six standard event shape variables (thrust, heavy jet mass, total jet broadening, wide jet broadening, C-parameter and the two-to-three jet transition variable y 23 in the Durham algorithm) and compare our predictions to those available in the literature [5,6]. We also present here for the first time the computation of jet cone energy fraction (JCEF) at NNLO accuracy. Predictions from CoLoRFulNNLO at this order in perturbation theory for oblateness and energy-energy correlation (EEC) were presented in ref. [7].

Definition of event shapes
Thrust [76,77] is defined as where the three-vectors p i denote the three-momenta of the partons and n defines the direction of the thrust axis, n T , by maximizing the sum on the right-hand side. For massless particles thrust is normalized by the center-of-mass energy, i | p i | = Q. In general 1/2 ≤ T ≤ 1, with T = 1/2 for spherically symmetric events, and T → 1 in the case of two back-to-back jets (the dijet limit). For three-particle events, we have 2/3 ≤ T ≤ 1.
Heavy jet mass [78][79][80] is defined by dividing the event into two hemispheres, H L , H R , by a plane orthogonal to an axis which can be chosen to be the thrust axis n T . Then the hemisphere invariant mass is where E vis is the total visible energy measured in the event, which is equal to the center-of-mass energy in perturbation theory with massless partons, E vis = Q. The heavy jet mass is In the dijet limit, we find ρ → 0. For three-particle events we have 0 ≤ ρ ≤ 1/3. At leading order in perturbation theory the distributions of heavy jet mass ρ and τ ≡ 1 − T are identical.
Jet broadening [81,82], like heavy jet mass, is also defined through the two hemispheres H L , H R . First, hemisphere broadening is given by The total and wide jet broadening are then defined as In the dijet limit, both B T and B W vanish, while for spherically symmetric events B T = 2B W = π/8. For three-parton events we have The C-parameter [83,84] is defined through the eigenvalues λ 1 , λ 2 , λ 3 , of the infrared-safe momentum tensor, where i runs over all final state particles. As Θ is a symmetric non-negative tensor with unit trace, the eigenvalues λ i are real and non-negative, with i λ i = 1. Therefore, 0 ≤ λ i ≤ 1, with i = 1, 2, 3. The value of the C-parameter is then defined as C par = 3 (λ 1 λ 2 + λ 2 λ 3 + λ 3 λ 1 ) . (5.7) In the dijet limit the C-parameter vanishes, while for spherical events C par = 1, so 0 ≤ C par ≤ 1.
For events with three-partons in the final state we have 0 ≤ C par ≤ 3/4.
Jet transition variables specify how an event changes from a n-jet to a (n + 1)-jet configuration. For example, given a jet resolution parameter y cut , the two-to-three jet transition variable y 23 [85][86][87][88] is defined as the value of y cut for which an event changes from a two-jet to a three-jet configuration, within some jet algorithm. Here we focus on the Durham algorithm [88], which clusters particles into jets by computing the variable, for each pair (i, j) of particles. The pair with the lowest value of y ij is replaced by a pseudoparticle whose four-momentum is computed in the E recombination scheme, i.e., it is simply the sum of the four-momenta of particles i and j. This procedure is iterated until all pairs have y ij > y cut and the remaining pseudo-particles are the jets.
Finally, jet-cone energy fraction [89] is defined as the energy deposited within a conical shell of the opening angle χ between a particle and the thrust axis n T , whose direction is defined to point from the heavy jet mass hemisphere to the light jet mass hemisphere, In principle 0 o ≤ χ ≤ 180 o , but hard gluon emissions typically contribute only to the region 90 o ≤ χ ≤ 180 o , which is plotted in the data [90].

Event shapes revisited
In this section we present the predictions of the CoLoRFulNNLO method for the event shapes considered also in refs. [5,6]. To begin, we write the perturbative expansion of the differential distribution of an event shape observable O at the default renormalization scale (not to be confused with the regularization scale of section 2.3) µ 0 = Q 2 (the total center-of-mass energy) as where α s = α s (µ 0 ) and σ 0 is the leading-order perturbative prediction for the total cross section of the process e + e − → hadrons. The LO and NLO perturbative coefficients A(O) and B(O) for thrust, heavy jet mass, total and wide jet broadening, C-parameter and the jet transition variable y 23 in the Durham algorithm were computed a long time ago [91], while predictions for the NNLO coefficients C(O) were presented in [5,6] 3 . However, experiments measure the distributions normalized to the total hadronic cross section, σ, thus physical predictions should be normalized to that. At the default renormalization scale µ 0 , distributions normalized to the total hadronic cross section can be obtained from the expansion in eq. (5.10) above by multiplying with the inverse of where [92][93][94] The renormalization scale dependence of a three-jet event shape distribution normalized to the total hadronic cross section can be computed as (5.14) with ξ R ≡ µ/µ 0 . Using three-loop running, the scale dependence of the strong coupling is given by with t = ln(µ 2 /Λ 2 QCD ). The first two coefficients in the expansion of the β function, In order to compare to published predictions, we use α s (m Z ) = 0.118, corresponding to Λ QCD = 208 MeV.
We present physical predictions at the first three orders in perturbation theory for the distributions of the six event shapes in figures 1-3. In the upper panels we show our fixed-order predictions as well as those of the publicly available code EERAD3 [96] 4 , together with the measured data by the ALEPH collaboration. We present our predictions at LO and NLO accuracy as smooth curves and as histogram at NNLO to represent the numbers of the numerical integration as precisely as possible. We observe a very good numerical convergence of our method at NNLO. The bands in the upper panel correspond to the variation of the renormalization scale in the range ξ R ∈ [0.5, 2]. In order to make the scale dependence at NNLO accuracy more visible, we show the relative scale uncertainty on the middle and bottom panels of each figure. It is remarkable that the relative scale dependence is below 5 % for most of the distributions in the ranges that are most relevant for measuring the strong coupling. Nevertheless, there is still a sizable difference between the NNLO predictions and the data for most of the distributions, which we attribute to parton shower (or resummation) and to hadronization effects.
The dependence on the renormalization scale increases significantly beyond kinematical regions of three-parton contributions, for instance for τ > 1/3 or C par > 3/4 and thus those are not shown on the ratio plots. At these three-parton kinematical limits large logarithms appear inside the physical region which have to be resummed similarly to the large logarithms that appear for small values of the event shapes (at the boundary of the physical region) [97]. The lower panels show the ratio of the (updated but unpublished -see text) predictions of [6] (SW) and EERAD3 [96] (GGGH) to CoLoRFulNNLO (this work). The bands on the lower panels show the relative scale uncertainty of our NNLO results.  figure 1 for the C-parameter (C par ) and the two-to-three jet transition parameter (y 23 ) in the Durham clustering algorithm.
As mentioned above, predictions for these six event shapes were presented in refs. [5,6]. In order to quantify the level of agreement across the available perturbative predictions, we also show the ratio of the (updated but unpublished -see below) results of ref. [6] (denoted by SW) and those of EERAD3 (denoted by GGGH) normalized to ours in the middle and bottom panels of each figure. Since the published predictions of [6] are known to be affected by an issue with the phase space generation in the code used to compute those results [98], we have made comparisons to updated but unpublished results which where provided to us by S. Weinzierl 5 . The general conclusion one may draw is that our predictions are in agreement with the updated predictions of SW except for very small and large (beyond the kinematic limits at LO) values of the event shapes, up to the estimated statistical uncertainties. A qualitatively similar statement can be made about the comparison to the GGGH predictions, although the deviations from our results at small values of event shapes are in general more pronounced than for the SW predictions. This is especially apparent for the C-parameter distribution below C par = 0.1. Examining the plots in figure 4, we see that the agreement is generally quite good between the predictions of SW and CoLoRFulNNLO and reasonably good between GGGH and CoLoRFulNNLO. However the precise comparison to GGGH predictions is hampered by the somewhat large integration uncertainties and bin-to-bin fluctuations of those results. Also, significant deviations among the three predictions are visible for small and large values of the event shapes. For example for τ = 1 − T the differences between the CoLoRFulNNLO results and the other two computations grow up to a factor of two for τ > 1/3. However, in this region, the contribution from three-particle final states vanishes and the thrust distribution is determined by a four-jet final state. Thus, C(τ ) is determined by the NLO corrections to four-jet production, which have been known for a long time [8,9] and can also be computed with modern automated tools such as MadGraph5 aMC@NLO [99]. We have checked that our predictions are in complete agreement with those of MadGraph5 aMC@NLO. The same is true for the tails of the other distributions beyond their respective kinematic limits. For small values of the event shapes we have checked that our predictions agree with the resummed predictions obtained from SCET [17,20,100] expanded to O(α 3 s ). coefficients of the thrust, heavy jet mass, total and wide jet broadening, C-parameter and two-to-three jet transition variable y 23 distributions. Lower panels show the ratio of the (updated but unpublished -see text) predictions of [6] (SW) and EERAD3 [96] (GGGH) to CoLoRFulNNLO (this work). The shaded bands on the lower panels represent the relative statistical uncertainties of our predictions due to Monte Carlo integrations.

Jet cone energy fraction
The jet cone energy fraction defined in eq. (5.9) is a particularly simple and excellent observable for the determination of the strong coupling. The smallness of hadronization corrections, detector corrections as well as perturbative corrections allows a specially wide fit range to be used for the extraction of α s [90]. JCEF was computed at NLO accuracy for the first time in ref. [89]. Here we present the first result for the JCEF distribution at NNLO accuracy in perturbative QCD for collider energy of Q 2 = 91.2 GeV. In figure 5 we show physical predictions for JCEF, as well as the measured data by the DELPHI collaboration. As previously, the uncertainties due to the variation of the renormalization scale in the range [0.5, 2] times our default scale choice (the total center-of-mass energy) are shown as bands on the upper panel. We indicate the relative scale uncertainty at NNLO on the bottom panel. To better appreciate the impact of the NNLO corrections, we show in figure 6 the distribution of the NNLO coefficient C(χ) directly. Also for these distributions, we observe a good numerical convergence of our code.

Conclusions and outlook
In this paper we presented the CoLoRFulNNLO framework to compute higher order radiative corrections to jet cross sections in perturbative QCD. CoLoRFulNNLO is a completely local and fully differential subtraction scheme based on the known infrared factorization properties of QCD matrix elements in soft and collinear limits. Since the subtraction terms explicitly take all color and spin correlations into account, the regularized real emission terms (both double real and real-virtual) are well-defined and can be computed in four dimensions with whatever numerical procedure is deemed most convenient. We have shown analytically that explicit infrared -poles coming from loop amplitudes cancel against the integrated forms of subtraction terms both in the real-virtual contribution (for any number of jets) and double virtual contribution (with up to three jets in the final state).
We have also reported on the computation of NNLO corrections to three-jet event shape observables in electron-positron collisions using CoLoRFulNNLO. We observe a very good numerical convergence of our method, which we attribute at least in part to the complete locality of the subtraction terms.
We compared both our physical predictions as well as the NNLO contribution only with similar predictions published earlier (in ref. [96] denoted by GGGH and in ref. [6] denoted by SW) for thrust, heavy jet mass, total and wide jet broadening, the C-parameter and the twoto-three jet transition variable y 23 in the Durham jet clustering algorithm. We find agreement with the updated (unpublished) predictions of SW within the statistical uncertainty of the numerical integrations except for very small and large values of the event shapes, beyond the kinematic limits at LO. The measured data in these regions are limited by statistics and so the phenomenological relevance of the differences is negligible. The same is true for the comparison to the physical predictions of GGGH, however the deviations from our results for small values of event shapes are generally more pronounced than for the predictions of SW. This is especially apparent for small values of the C-parameter, below C par = 0.1. When comparing our physical predictions at the NNLO accuracy to experimental data, we still find large differences. An important source of discrepancy is the neglected large logarithmic contributions which require all-order resummation. Work towards matching the fixed-order predictions to resummed ones is in progress.
Finally, we have shown for the first time perturbative predictions for jet cone energy fraction at NNLO, thereby providing a new observable from which the value of the strong coupling can be extracted at this accuracy. This observable has the remarkable property that the NNLO corrections are very small and the fairly good agreement between data and predictions already at NLO become only marginally better with the inclusion of the NNLO corrections. This stability of the perturbative predictions makes JCEF a good candidate for the extraction of the strong coupling.
We emphasize that our framework is not restricted to three-jet production, but it can be easily applied to study differential distributions for four-or more jet production once the necessary two-loop matrix elements become available. CoLoRFulNNLO is completely worked out for processes with colorless initial states at present. The inclusion of initial state radiation is a conceptually straightforward although substantial task and is work in progress.
A The I  − double virtual -poles at NNLO. However, in the numerical integrations over the three-parton phase space one also needs the finite part, which is rather lengthy and in fact, we have only computed its asymptotic expansion for small kinematic invariants analytically. Here we record this expansion and comment on the remaining regular part, which we compute numerically.
For the Born process e + e − → qqg the J 2 operator can be written in the following explicit form: which then also defines the finite part Fin J 2 ( ) unambiguously. We decompose the finite part into an asymptotic piece which collects all logarithmic contributions that become singular on the borders of the three-parton phase space and a piece which is regular over the whole phase space, i.e. finite on the borders: Fin J 2 ( ) = Fin J asy 2 ( ) + Fin J reg 2 ( ) .

(B.2)
The asymptotic part can be written as We stress that this form as well as the explicit expressions for the asymptotic functions presented below are known to be appropriate only for e + e − → 3 jet production. The analytic expressions We do not have analytic expressions for the regular part. However, computing this piece numerically on a grid over the three-parton phase space, we find that it is in fact flat across the whole phase space (within the uncertainty of the numerical integrations). Hence it can be described by a single number whose numerical value we find to be for N c = 3, T R = 1/2 and n f = 5.