A subtraction scheme for computing QCD jet cross sections at NNLO: integrating the doubly unresolved subtraction terms

We finish the definition of a subtraction scheme for computing NNLO corrections to QCD jet cross sections. In particular, we perform the integration of the soft-type contributions to the doubly unresolved counterterms via the method of Mellin-Barnes representations. With these final ingredients in place, the definition of the scheme is complete and the computation of fully differential rates for electron-positron annihilation into two and three jets at NNLO accuracy becomes feasible.


Introduction
The high precision of experimental measurements of benchmark processes such as jet, heavy quark, vector boson and Higgs boson production at the LHC calls for an equally precise theoretical evaluation of relevant production rates. In particular, exact perturbative predictions beyond next-to-leading order accuracy are desirable and sometimes essential to reduce the theoretical uncertainty. Accordingly, in recent years considerable effort has beed devoted to computing next-to-next-to-leading order (NNLO) corrections to various production and decay rates. Fully differential cross sections have been evaluated for vector boson [1,2], Higgs boson [3,4], diphoton [5], and Higgs-vector boson [6] production, and for dijet production in the fully gluonic channel at leading colour [7]. Fully differential decay rates for top [8] and bottom [9] quark decay have also been calculated, and the computation of the total cross sections for top-antitop production [10][11][12] and Higgs boson production in association with a jet [13] at NNLO is also underway. Jet rates and event shapes in electron-positron annihilation to two and three jets have also been computed [14][15][16][17][18][19][20][21][22].
The construction of a general subtraction scheme to calculate cross sections at NNLO accuracy has also received significant attention [23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41]. Recall that the the NNLO correction to the cross section for producing m jets is a sum of three contributions, the doubly real, real-virtual and doubly virtual terms, where each contribution on the right-hand side is separately infrared divergent in four spacetime dimensions. In Eq. (1.1) J n is a generic infrared-safe jet measurement function. The essential idea of subtraction is that by subtracting and adding back suitable approximate cross sections, we can reshuffle divergent pieces among the three terms in Eq. (1.1). This reshuffling amounts to writing Eq. (1.1) as follows, (1.5) are each integrable in four dimensions by construction [24,26,27,42]. The various approximate cross sections appearing in Eqs. (1.3)-(1.5) are as follows.
• dσ RR,A 2 m+2 regulates the doubly unresolved limits of the double real emission cross section, dσ RR m+2 .
• dσ RR,A 1 m+2 regulates the singly unresolved limits of the double real emission cross section, dσ RR m+2 .
• dσ RR,A 12 m+2 compensates for the double subtraction due to the overlap of singly and doubly unresolved limits, i.e., it regulates the singly unresolved limits of dσ RR,A 2 m+2 and the doubly unresolved limits of dσ RR,A 1 m+2 .
A 1 regulates the singly unresolved limits of the integrated approximate cross section 1 dσ Importantly, just as at NLO accuracy [31,[42][43][44][45], all approximate cross sections appearing in Eqs. (1.3)-(1.5) can be constructed once and for all, independently of the process or observable being studied. Indeed, the approximate cross sections that contribute to dσ NNLO m+2 and dσ NNLO m+1 in Eqs. (1.3) and (1.4) are precisely defined in Refs. [26] and [27] respectively. We recall that all subtraction terms are constructed using the known universal limit formulae for collinear and soft QCD radiation, and their various iterated forms [24]. Hence, each approximate cross section is a sum of several contributions which may be classified according to the kind of limit as collinear-or soft-type.
In order to make the subtraction scheme useful for practical calculations, one must also compute the integrals of the approximate cross sections over the phase spaces of unresolved partons. The integrated approximate cross section 1 dσ RR,A 1 m+2 that appears in Eq. (1.4) was evaluated in Ref. [42]. The integral of the real-virtual approximate cross sections, denoted formally by 1 in Eq. (1.5), were computed in Refs. [29,30,32], while 2 dσ RR,A 12 m+2 was evaluated in Ref. [36]. Finally, in Ref. [46], we calculated all collinear-type contributions to the integrated doubly unresolved approximate cross section, 2 dσ RR,A 2 m+2 . In this paper, we complete the definition of the subtraction scheme by computing the double soft-type contributions to 2 dσ RR,A 2 m+2 . We use the method of Mellin-Barnes representations -developed in this context in Ref. [32] -to compute the various integrals that appear as building blocks for the complete expression. With the results of the present paper, the computation of the finite m-parton cross section in Eq. (1.5) becomes feasible. Since the finiteness of the m + 2-and m + 1-parton cross sections in Eqs. (1.3) and (1.4) was demonstrated in Refs. [26] and [27] (for m = 3 specifically), we are now in a position to calculate the fully differential rate for two-and three-jet production in electron-positron annihilation at NNLO accuracy within our framework. For four or more jets some more work is required, since in Refs. [32] and [46] (as well as in this paper) a few integrals were computed specifically for three-jet kinematics. Having said that, we stress that at present the two-loop matrix elements relevant for four or more jet production are unknown and these -rather than the few as yet uncomputed integrals in our scheme -are the essential missing ingredients.
The rest of the paper is organised as follows. In Section 2, we set the notation used throughout. The structure of the integrated doubly unresolved cross section is recalled next in Section 3. It is the product (in colour space) of the Born cross section times an appropriate insertion operator, which itself is given in terms of so-called flavour-summed integrated counterterms. In Section 4 we present these counterterms explicitly, recalling their decomposition into non flavour-summed subtraction terms, which in turn are defined in terms of integrals of soft currents or precisely defined limits of soft currents. These integrals are then expressed as linear combinations of a set of basic integrals. Then, in Section 5 we discuss our approach to evaluating these basic integrals. Our final results are presented in Section 6, including analytic expressions for the integrated counterterms up to O(ǫ −2 ) and illustrative results for the complete insertion operator for two-and three-jet production in electron-positron annihilation. Finally, in Section 7 we draw our conclusions.

Matrix elements and cross sections
Our notation was spelled out extensively previously [36,46] and is recalled here only to the extent that we will need in this paper.
The processes we consider are decays of a colourless initial state into any number of massless coloured patrons plus any number of additional non-coloured final state particles, that are however suppressed in the notation. We use letters form the middle of the alphabet, i, j, k, l, . . . , to denote resolved final state patrons, while r and s will label the two soft partons. The matrix element for a process involving n final state coloured partons will be denoted by |M n , using the colour-and spin-state notation of Ref. [44]. The squared matrix element, is understood to be summed over colours and spins, which fixes the normalisation. In this paper, we will only need to use tree level matrix elements, which we denote |M (0) n . Two-and four-parton colour correlated squared tree amplitudes are defined as follows where T i etc., is the usual colour-charge operator associated with the emission of a gluon from parton i. The square of the colour-charge operator, T 2 i , depends only on the flavour of the emitting parton. We emphasise this by introducing the notation i.e., C f i is the quadratic Casimir operator in the representation of parton i. Hence, We also use squared colour-charge operators with multiple indices, e.g., T 2 ir ≡ C f ir and T 2 irs ≡ C f irs . In such cases, the multiple index denotes a single parton with flavour obtained according to usual flavour summation rules: an odd (even) number of quarks and any number of gluons gives a quark (gluon).
Finally, dσ (0) n ({p}) denotes the fully differential cross section for producing n partons at tree level with momenta {p} ≡ {p 1 , . . . , p n }. Its precise definition reads 5) where N collects factors that are independent of QCD, i.e., it is the flux factor times the spin average factor for the incoming particles; the symbol {n} denotes summation over different subprocesses; S {n} is the Bose symmetry factor for identical particles in the final state; and dφ n ({p}; Q) is the phase space associated with final state momenta {p} = {p 1 , . . . , p n } and total momentum Q, In this paper, we will use the cross sections for producing m and m + 2 partons, which we call the Born and doubly real contributions, respectively

Phase-space factorisation and kinematic variables
The double soft-type subtraction terms are all defined precisely via the double soft momentum mapping of Ref. [26], which maps the original set of m + 2 momenta into a set of m tilded momenta {p} m+2 The explicit form of this mapping can be found in Ref. [26]. Here we recall only that it leads to an exact factorisation of phase space in the following form, where the factorised phase-space measure is [dp with K µ a massive momentum such that K 2 = λ 2 rs Q 2 . We will have use for the iterated double soft momentum mapping as well, also defined in Ref. [26]. This mapping again leads to the exact factorisation of phase space, now in the iterated form The factorised phase-space measures read [dp where t = s and n = m + 1 or t =r and n = m, while K µ is a massive momentum with K 2 = (1 − y tQ )Q 2 and y tQ ≡ 2p t · Q/Q 2 . The two-parton phase space dφ 2 (p t , K; Q) can be written as follows in the rest frame of Q, where ε t = 2E t / Q 2 is the scaled energy of parton t in the Q rest frame, while dΩ d−1 (t) is the angular measure in d − 1 dimensions in this frame. In Eq. (2.14), S ǫ is 1 . (2.15) Using the precise definitions of the single and double soft momentum mappings [26], it is easy to show that the two sets of tilded momenta on the right hand sides of Eqs. respectively, are related simply by a three-dimensional rotation. This will be a key element for the integration of counterterms in Section 5.
Finally, we recall the definitions of kinematic variables. First, two-particle invariants, such as s ik , s iQ , sr Q , etc., always denote twice the dot product of two momenta. A double index enclosed in brackets, as in e.g., s i(rs) , indicates a sum of momenta. Hence When these invariants are normalised to the total incoming momentum squared, we use the notation y ik ≡ s ik /Q 2 . The momentum fractions for double and triple parton splitting, z i,r and z i,rs , are defined as follows with z r,i as well z r,is and z s,ir given by obvious cyclic permutations. Clearly z i,r + z r,i = 1 and z i,rs + z r,is + z s,ir = 1. Furthermore, S ik (r) denotes the eikonal factor, Lastly, the integrated counterterms turn out to be functions of the following two types of variables

The integrated doubly unresolved approximate cross section
We recall from Ref. [46] that after summing over the flavours of unobserved partons, the integrated doubly real approximate cross section can be written as Here the insertion operator I where the sums over i, j, k and l run over all external partons. Let us briefly elaborate on the physical content of this formula. As mentioned in the Introduction, the approximate cross section dσ is actually a sum of terms over the various doubly unresolved kinematical limits (triple collinear, double collinear, soft collinear and double soft) and their specific iterated forms (which remove multiple subtractions when limits overlap), ir;js ({p}) Integrating the approximate cross section, we obtain the corresponding sum of integrated counterterms, where we have factored out quadratic Casimir operators to make the non flavour-summed integrated subtraction terms, denoted generically as [ in Eq. (3.4) with the doubly virtual cross section dσ VV m , it must be put in the form of an m-parton contribution times a factor, which involves rewriting the symmetry factor of the m + 2-parton configuration to the symmetry factor of an m-parton configuration, i.e., summing over unobserved flavours. Performing this rewriting, we obtain flavour-summed integrated counterterms, denoted X (0) (j,l)...
, which are specific sums of the non flavour-summed integrated subtraction terms, see Ref. [46] and Section 4 below. Finally, the flavour-summed integrated counterterms may be organised according to their structure in colour and flavour space into the five terms of Eq. (3.2), 2,i and C (0) 2,ij collect terms that come from taking triple and double collinear limits, respectively. They depend on the flavour(s) of the hard mother parton(s) in the collinear splitting(s) and multiply the Born matrix element with no colour correlations. All terms coming from taking soft collinear limits are gathered into CS (0),(j,l) 2,i which depends on the flavour of the hard mother parton in the collinear splitting and multiplies the twoparton colour correlated Born matrix element. Finally, S (0),(j,l) 2 and S (0),(i,k)(j,l) 2 correspond to parts of the double soft limit that multiply the two-and four-parton colour correlated Born matrix element, respectively. They are both independent of the flavours of hard partons. The structure of each term in colour (flavour) space is indicated by upper (lower) indices.
All flavour-summed functions appearing on the right-hand side of Eq. (3.5) were given in terms of non flavour-summed functions in Ref. [46], where we also computed all those functions in Eq. (3.5) which do not involve double soft contributions. In this paper, we compute the remaining flavour-summed double soft functions.

Integrated double soft-type counterterms
In this section we define explicitly all double soft-type integrated counterterms that appear on the right-hand side of Eq. (3.5). In each case, we first present the flavour decomposition of the flavour-summed counterterms. Then, we give the proper definition of the non flavoursummed functions in terms of integrals of soft currents or precisely defined limits of soft currents. Finally, we compute the results in terms of a set of basic integrals.
Double soft: The flavour decomposition of the double soft counterterms reads where n f is the number of light quark flavours. The functions appearing on the right-hand side of Eq. (4.1) are defined as follows The non-abelian part of the double soft gluon current, S ik (r, s) in Eq. (4.3), reads [47] S ik (r, is the form of this function in the strongly-ordered limit (either E r ≪ E s or E s ≪ E r as the expression is symmetric in r and s when summed over i and k) and S ik (rs) is given by .
Finally, the function f (y 0 , y rQ + y sQ − y rs , d ′ (m, ǫ)) appearing in Eqs. (4.2)-(4.4) represents a small but convenient modification of the subtraction terms as compared to their original definitions in Ref. [26]. Its precise role and form are explained in Appendix A, and in the following, we will include this factor without further comment.
In this paper, we do not discuss the case when i, k, j and l in Eq. (4.2) are all distinct, this requiring at least four jets at NNLO. For the specific cases of two and three hard partons, the integrated double soft counterterms are computed in Appendix B and we find The integrals I 2S ,n for n = 1, . . . , 9 are defined in Eqs. (B.6)-(B.14).
Triple collinear -double soft: The flavour decomposition of the triple collineardouble soft counterterm reads We further decompose the triple collinear -double soft gluon counterterm as a sum of abelian and non-abelian pieces, following the decomposition of the triple Altarelli-Parisi splitting kernels in a similar fashion [47]. Then we have the following explicit definitions for the non flavour-summed func- (nab) Double collinear -double soft: The flavour decomposition of the double collineardouble soft counterterm is trivial, × f (y 0 , y rQ + y sQ − y rs , d ′ (m, ǫ)) .

(4.24)
We remind the reader that above j = l, but (ir) may be equal to j or l. The cases when j, l = (ir) and when e.g., j = (ir) lead to different integrals, as discussed in Appendix E.  ǫ)) .

(4.28)
This integral is computed in Appendix F. We find

Simplifying the integrals
Consider a generic double soft-type master integral, where F n ({s jl , z j }; ǫ) is the integrand, which will depend on two-particle invariants s jl and possibly also momentum fractions z j . (We use z j generically to denote both the two-and three-parton momentum fractions defined in Eq. (2.17).) A straightforward treatment of Eq. (5.1) proceeds to first express the integrand in terms of independent (tilded) momenta. Then, choosing some particular Lorentz frame, the integral is written in terms of e.g., the angles and energies of momenta p µ r and p µ s in this frame. The result obtained however turns out to be completely unwieldy. In particular, the extraction of ǫ poles via sector decomposition is not possible, because one typically finds complicated singularities inside the domain of integration. It also turns out that such a parametrisation is not a useful starting point for deriving Mellin-Barnes representations.
A much more manageable form is obtained, however, by observing that the value of the integral in Eq. (5.1) is unchanged if we replace the double soft unresolved phase-space measure [dp (rs) 2;m ] by the iterated single soft phase-space measure, e.g., [dp To understand why this may be the case, note that by using the precise definitions of the single and double soft momentum mappings [26], it is easy to show that the two sets of tilded momenta {p} are related simply by a three-dimensional rotation. This implies that the dot products of momenta in each set are the same for both sets. But the integrated counterterms depend only on these dot products (as opposed to the overall orientation of the tilded momenta), hence we are free to use the more convenient iterated phase-space mapping when computing the integrals, 2 [dp (rs) 2;m ] F n ({s jl , z j }; ǫ)f (y 0 , y rQ + y sQ − y rs , d ′ (m, ǫ)) = 2 [dp (r) 1;m ][dp (s) 1;m+1 ] F n ({s jl , z j }; ǫ)f (y 0 , y rQ + y sQ − y rs , d ′ (m, ǫ)) .

(5.3)
By using the iterated form of the phase-space measure, the integrations over the variables ofp µ r and p µ s may be performed sequentially.

Computing the integrals via Mellin-Barnes representations
Following the discussion of the previous section, we write the generic master integral in Eq. (5.1) as [dp Since the dimension of F n is [F n ] = (Q 2 ) −2 , we have F n ({s jl , z j }) = 1 (Q 2 ) 2 F n ({y jl , z j,kl }) , (5.5) and hence, using Eqs. (2.13), (2.14) and the exact definition of f (y 0 , y rQ +y sQ −y rs , d ′ (m, ǫ)) in Eq. (A.4), we find In the above form, Eq. (5.6) is directly suitable for treatment by Mellin-Barnes techniques.
In particular, the angular integrations over the directions ofp µ r and p µ s may be preformed sequentially using the results of Ref. [48]. Hence, it becomes essentially straightforward to derive a Mellin-Barnes representation for the integral: 2 1. Substitute the specific form of F n .
2. Perform the angular integration over dΩ d−1 (s), using the results of Ref. [48]. The momenta which are to be kept fixed during this integration are the intermediate, hatted momenta of Eq. (2.11). Thus, at this stage, all invariants need to be expressed in terms of these: The result after this angular integration depends on (scaled) dot products of hatted momenta and perhaps Q, i.e., yîk, yîr, yî Q , etc.
3. Repeat the previous step for the angular integrals over dΩ d−1 (r). The independent momenta with regard to thep µ r integration are now the final set of tilded momenta in Eq. (2.11), so we must express all invariants in terms of these: After this integration, the expression depends only on invariants involving tilded momenta, as well as yr Q and y sQ , which are however integration variables themselves.
4. Finally, perform the integrations over yr Q and y sQ to obtain the full Mellin-Barnes representation. If y 0 = 1, the argument of the Θ function in Eq. (5.6) factorises and these last two integrations are automatically in the form of a beta function integral. For 0 < y 0 < 1, additional Mellin-Barnes integrations must be introduced to bring them to this form.
With this procedure, we obtain a multi-dimensional Mellin-Barnes integral representation of Eq. (5.6). It is usually not possible to evaluate these integrals analytically. However, for practical purposes, we are actually interested in the ǫ-expansion of I 2S ,n , rather than the all-orders result. Importantly, this expansion may be performed in an algorithmic way, e.g., with the MB.m [50] or MBresolve.m [51] packages. Then we obtain the expansion coefficients as finite Mellin-Barnes integrals. We are able to evaluate all integrals that contribute to the 1/ǫ 4 and 1/ǫ 3 poles analytically, 3 while the rest of the expansion coefficients can be computed by direct numerical integration of the Mellin-Barnes representation.
Finally, we call attention to the following technical detail. We find that it is best to fix the specific value of d ′ 0 (see Eq. (A.3)) before the ǫ-expansion of the integrals. Therefore, in the following, we set d ′ 0 = 3 − 3ǫ, i.e., D ′ 0 = 3 and d ′ 1 = −3. This is not an essential restriction and recomputing the expansions for different values of D ′ 0 and d ′ 1 is in principle straightforward.

Analytic expressions to O(ǫ −2 )
We have obtained analytic expressions for the integrated counterterms up to O(ǫ −2 ) accuracy. The results below are quoted for d ′ 0 = 3 − 3ǫ, i.e., D ′ 0 = 3 and d ′ 1 = −3 (see Eq. (A.3)), and generic y 0 ∈ (0, 1]. Before presenting the actual formulae, we call attention to the following features. Up to this pole order: • All kinematic dependence enters through logarithms of the variables x i and/or Y ik,Q defined in Eq. (2.19).
• Dependence on the cut parameter y 0 always enters in the same functional form, as follows, We observe that this is simply Σ(y 0 , 2) (recall that throughout we use D ′ 0 = 3), with the Σ(z, N ) function of Refs. [29] and [36], • Additional constants that enter are always recognised to be slight generalisations of γ q /C q and γ g /C g (since our integrated counterterms are dimensionless in colour space), where γ g = 3 2 C F and γ g = 11 6 C A − 2 3 T R n f [52]. In order to emphasise this connection, we define the following two functions of n f [46], The formal n f dependence of γ q (n f ) is introduced in order to make possible a flavourindependent notation in the following.
We find 1. Double soft: note the lack of Y ij,Q dependence to this order, and 2. Triple collinear -double soft: Note that this counterterm does not depend on kinematics.
3. Double collinear -double soft: Note the lack of Y ij,Q dependence to this order.

Soft collinear -double soft:
for j, l = i, and for e.g., j = i.

Triple collinear -soft collinear -double soft:
Substituting Eqs. (6.4)-(6.11) and the corresponding results from Ref. [46] -recalled here in Appendix G for the convenience of the reader -into Eq. (3.5), we obtain explicit expressions for the kinematics-dependent functions entering the insertion operator. We note that since in this paper we set D ′ 0 = 3, for consistency we must also use the results of Ref. [46] with D ′ 0 = 3. The results do not depend on d 0 and α 0 up to this order. 4 Starting with C (0) 2,i , we find However, C 2,ij enters the insertion operator I 2 summed over its flavour indices, thus we are free to symmetrise in i and j. In particular, setting 14) (we use the usual notation of round brackets around indices to denote symmetrisation) we see that up to this order in the ǫ expansion, C 2,(ij) simply vanishes, For CS (0),(j,l) 2,i we find (note the vanishing of Y jl,Q dependence to this order) when j and l are distinct form i, and for e.g., j = i. The two-parton colour-correlated soft function, S and α0 are parameters of the collinear-type counterterms that correspond to the d ′ 0 and y0 used in this paper. In particular, α0 ∈ (0, 1], if smaller than one, restricts subtractions to near collinear regions in phase space. See Ref. [46] for further details. Last, we consider the four-parton colour-connected soft function, S (0)(i,k)(j,l) 2 . We find (note that Y ij,Q dependence is absent to this order) if only two indices are different.

Insertion operator for two-and three-jet production
Using the results of Ref. [46] and the present paper, we can assemble the complete insertion operator I 2 relevant for two-and three-jet production. Let us consider first the process e + e − → 2 jets. The corresponding squared matrix element at tree level is |M (0) 2 (1 q , 2q)| 2 , i.e., the quark carries label 1 and the antiquark label 2. Both the colour algebra and kinematics are trivial. Colour conservation implies and C f 1 = C f 2 = C F , while momentum conservation gives Then the insertion operator, Eq. (3.2), is a scalar in colour space and reads with all arguments being equal to one. Using Eqs. (6.12)-(6.20), we find where we used the notation [53] Notice that up to this order, I 2 is independent of y 0 . The remaining expansion coefficients are computed numerically. We present these in Table 1 for the quantity I  for e + e − → 2 jet production. We by the scalar product of the vector of numbers in the appropriate column of Table 1 with the vector of coefficients forming the first column.
Turning to the process e + e − → 3 jets, the corresponding squared matrix element at tree level is |M (0) 2 (1 q , 2q, 3 g )| 2 , i.e., the quark carries label 1, the antiquark label 2 and the gluon carries label 3. The colour algebra is still trivial and using colour conservation we find and The insertion operator, Eq. (3.2), is again a scalar in colour space and we have where we used ln Y ik,Q = ln y ik −ln x i −ln x k . Note that the dependence on y 0 again cancels up to this order in the expansion. The rest of the coefficients in the ǫ-expansion are computed numerically. By way of illustration, we present results for the fully symmetric configuration of final state momenta in Table 2. In this phase-space point, the various invariants take the following values for e + e − → 3 jet production.
where i, k = 1, 2, 3 and i = k. We have again used d 0 = d ′ 0 = 3 − 3ǫ and α 0 = y 0 = 1 during the calculation. As for the two-jet case, colour factors and the number of light flavours were kept symbolic. The actual values of the expansion coefficients are again scalar products of vectors of numbers in the last three columns of Table 2 with the vector of coefficients which forms the first column. Note that the numbers presented in Table 2 correspond to the expansion coefficients of I

Conclusions
This paper finishes the calculation of the integrated doubly unresolved approximate cross section of the NNLO subtraction formalism of Refs. [26,27], and thus completes the definition of the subtraction scheme.
In particular, here we computed the double soft-type contributions to the integrated doubly unresolved approximate cross section of Ref. [26]. The integrated counterterms were evaluated in terms of a set of basic double soft integrals, which were calculated using Mellin-Barnes techniques. These contributions represented the last missing ingredients needed to evaluate 2 dσ RR,A 2 m+2 , as the collinear pieces of the doubly unresolved subtraction terms were integrated in Ref. [46]. The final result can be written as the product (in colour space) of the Born cross section times the doubly unresolved insertion operator, I 2 . We were able to compute this insertion operator analytically up to O(ǫ −2 ), while the rest of the expansion coefficients were evaluated numerically.
As stressed above, the definition of the NNLO subtraction formalism of Refs. [26,27] is now complete, and the evaluation of the finite doubly virtual cross section of Eq. (1.5) is feasible for electron-positron annihilation into two-and three-jets. (A few integrals were evaluated specifically for three-jet kinematics, so for a higher number of jets, some more work is required.) In particular, all integrated approximate cross sections appearing in Eq. (1.5) are known analytically up to O(ǫ −2 ) and hence, the cancellation of poles in dσ NNLO m may be checked explicitly to this order. Indeed, we have checked that the 1/ǫ 4 and 1/ǫ 3 poles cancel for e + e − → 2, 3 jets, independently of the value of y 0 . The cancellation of the subleading poles must be checked numerically, and the details will be presented elsewhere. As the regularised doubly real and real-virtual cross sections dσ NNLO m+2 and dσ NNLO  4) were previously shown to be finite [26,27] (specifically for m = 3), we are now in the position to compute the fully differential rate for electron-positron annihilation into two and three jets at NNLO accuracy within our framework.

Acknowledgments
I am grateful to Vittorio Del Duca and Zoltán Trócsányi for useful conversations and comments on the manuscript.

A Modified double soft-type subtraction terms
In previous publications, we have presented an easy modification to the NNLO subtraction scheme of Refs. [26,27]. in Eq. (1.3), were spelled out in Ref. [36]. Finally, the modified doubly unresolved approximate cross section was discussed in Ref. [46].
The introduction of this modification serves two purposes. First, it makes the integrated subtraction terms independent of m, the number of hard partons, see Refs. [29,36,46] for a detailed discussion. Second, it allows to constrain the subtractions to near the singular regions in phase space, which improves the efficiency of the numerical implementation.
In this paper, we only need the precise definition of the modified double soft-type subtraction terms. We recall from Ref. [46] that these are obtained from the original subtraction terms in Ref. [26] via multiplication by the simple factor of f (y 0 , y rQ + y sQ − y rs , d ′ (m, ǫ)), where We also reproduce the relevant part of Table 20 of Ref. [46] here, as Table 3. We note that the form of d ′ (m, ǫ) appearing in Table 3 is actually fixed by the prescription in Ref. [29] and the requirement that the modified subtraction terms still correctly regularise all kinematic singularities. We have where, as in Refs. [29,36,46], the parameter d ′ 0 takes the form  Table 3. The modified double soft-type subtraction terms are obtained from the original counterterms (first column) by multiplication with f (y 0 , y rQ + y sQ − y rs , d ′ (m, ε)) (last column). Also shown is the momentum mapping used to define the subtraction terms (middle column).
with D ′ 0 ≥ 2 an integer and d ′ 1 real. Throughout this paper, we use the specific choice D ′ 0 = 3 and d ′ 1 = −3. Finally, we note that in terms of yr Q and y sQ , we have where S ik (r, s) is given in Eq. (4.5).
Let us proceed to identify the independent kinematic structures that we need to integrate. We begin with the abelian double soft gluon contribution, i.e., the first term in the list of Eq. (B.1). Recall that in this term, the indices i, j, k, l are constrained only by the requirements that i = k and j = l, but there is no restriction on whether or not i, k is equal to j, l. Hence, we have the following situations: (i) all of i, k, j and l are distinct, (ii) only three of the four indices are distinct, e.g., l = k and (iii) only two indices are distinct, e.g., j = i and l = k.
For case (i) to occur, there must be at least four hard patrons in the process. Obviously this does not happen in the calculation of two-and three-jet observables, and we will not consider it in this paper. For cases (ii) and (iii), we have simply S ik (r)S jk (s) = 4s ik s jk s ir s kr s js s ks and S ik (r)S ik (s) = 4s 2 ik s ir s kr s is s ks .

(B.2)
Turning to the non-abelian part of the double soft gluon formula, the second term in the list of Eq. (B.1), we make two observations. First, the factorised phase-space measure, [dp (rs) 2;m (p r , p s ; Q)], is clearly symmetric under the exchange of p µ r and p µ s (see Eq. (2.10)).
Second, the double soft current appears in the subtraction terms summed over its indices, i and k in this case. Hence, we are free to exchange i and k in individual terms in S ik (r, s) without changing the value of the total integrated subtraction term, since i and k are merely summation indices. Using the freedom to make r ↔ s and/or i ↔ k replacements whenever convenient, with some algebra 5 we can derive the following form of S ik (r, s), The specific normalisations were chosen for later convenience. Finally we introduce the following notation, [dp [dp s ik s i(rs) s kr s rs × f (y 0 , y rQ + y sQ − y rs , d ′ (m, ǫ)) , (B.9) [dp for j = i.