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

We perform the integration of all iterated singly-unresolved subtraction terms, as defined in ref. [1], over the two-particle factorized phase space. We also sum over the unresolved parton flavours. The final result can be written as a convolution (in colour space) of the Born cross section and an insertion operator. We spell out the insertion operator in terms of 24 basic integrals that are defined explicitly. We compute the coefficients of the Laurent-expansion of these integrals in two different ways, with the method of Mellin-Barnes representations and sector decomposition. Finally, we present the Laurent-expansion of the full insertion operator for the specific examples of electron-positron annihilation into two and three jets.

5 Insertion operator for two-and three-jet production 29 5.1 Two-jet production 29 5.2 Three-jet production 31

Introduction
The most frequently occurring final states in high energy particle collisions contain hadronic jets. Because of their large production cross sections, jet observables can be measured to a high statistical accuracy, and are therefore ideal for precision studies. Examples include: the measurement of the strong coupling, α s , from jet rates and event shapes in electronpositron annihilation; the determination of the gluon parton distribution function (and also α s ) in deep inelastic lepton-hadron scattering into two plus one jets; the measurement of parton distributions in hadron-hadron scattering from single jet inclusive production and vector boson plus jet production. Often the relevant observables are measured with experimental precision of a few per cent or better. Thus, theoretical predictions with the same level of accuracy are necessary to fully exploit the physics potential of these measurements. This usually requires the computation of next-to-next-to-leading order (NNLO) corrections in perturbative QCD. The straightforward calculation of jet cross sections in QCD perturbation theory is however hampered by the presence of infrared singularities in the intermediate stages of the calculation, which must be treated consistently before any numerical computation may be performed. At next-to-leading order (NLO) accuracy, using a subtraction scheme to handle infrared divergences is the approach of choice. Exploiting the fact that the kinematic singularities of QCD matrix elements are universal, one builds process and observable independent counterterms that simultaneously cancel both the kinematic singularities in real-emission phase space integrals and the explicit ǫ-poles in one-loop virtual corrections (here the use of dimensional regularisation in d = 4 − 2ǫ dimensions is implied).
At NNLO accuracy, the calculation of fully differential cross sections is a challenging problem, and various extensions of the subtraction method at NNLO have been proposed, see e.g. refs. [1][2][3][4][5][6][7][8][9][10][11][12]. In very broad terms, when setting up any subtraction algorithm, two quite distinct difficulties must be addressed. First, one must define subtraction terms that properly regularise the real-emission phase space integrals and second, one must combine the integrated form of these counterterms with the virtual contributions, so as to cancel the infrared divergences of the loop matrix elements. In a rigorous mathematical sense, the cancellation of both the kinematic singularities in the real-emission pieces and the explicit ǫ-poles in the virtual pieces must be local. On the one hand, this means that the subtraction terms and the real-emission contributions must tend to the same value in d dimensions, in all kinematic limits where the latter diverge. On the other hand, the cancellation of explicit ǫ-poles between the integrated subtraction terms and the virtual contributions must take place point-wise in phase space. This in particular implies that it is possible to write the integrated counterterms in such a form that they can be explicitly combined with virtual contributions, before phase space integration. From a practical point of view, the full locality of the subtraction scheme is also important to insure good numerical efficiency of the algorithm. Finally, the construction should be universal, i.e. independent of the process and observable being considered. This avoids the need to tediously adapt the algorithm to every specific problem.
However, defining universal subtraction terms that are completely local in the realemission phase space is rather delicate, and there is very little freedom to define these in such a way, that their integration over the unresolved phase spaces becomes convenient. One way to address these difficulties is to use counterterms that are not fully local, but whose complete analytic integration is tractable. For example, the antenna subtraction method [8][9][10][11] builds the subtraction terms from so-called antenna functions. These are simple enough to be integrated analytically, but they do not reproduce azimuthal correlations in gluon splitting, and the cancellation of ǫ poles in the real-virtual contributions is also nonlocal. Because of this, actual numerical computations with the antenna scheme, such as the calculation of total rates [13][14][15][16] and event shapes [17][18][19][20][21] in electron-positron annihilation, require the use of an auxiliary phase space slicing. Another option is to develop dedicated subtraction schemes that are applicable only to some specific processes, such as the production of colourless final states, vector bosons [22,23] or the Higgs boson [24], in hadron collisions. Then, one may even propose to dispense with the subtraction method altogether, and adopt a strategy such as sector decomposition (see [25] and references therein), where the full Laurent expansions (in ǫ) of the real-emission and virtual pieces are computed directly, and their finite pieces combined.
Nevertheless, despite the subtleties, it is possible to define completely local counterterms for real radiation, by first carefully matching the various QCD factorisation formulae for unresolved emission [4,26], and then extending the expressions obtained over the full phase space [1][2][3]. Recall that in the scheme of refs. [1][2][3][4], the NNLO correction to a generic m-jet cross section with no coloured particles in the initial state (work towards an extension to hadron-initiated processes is presented in ref. [27]), is rewritten as a sum of finite integrals are each integrable in four dimensions by construction. In eqs. (1.1)-(1.5) J n (n = m + 2, m + 1 and m) denotes the jet function which defines the infrared safe observable being calculated. All approximate cross sections appearing in eqs. (1.3)-(1.5) above have been defined explicitly in refs. [1][2][3]. To finish the definition of the scheme, we must compute once and for all the one-and two-particle integrals appearing in eqs. (1.4) and (1.5). In previous publications, we have shown that it is possible to adapt and employ well-known techniques of loop integration, such as integration-by-parts identities and solving of differential equations [28,29], the method of Mellin-Barnes (MB) representations with harmonic summation [30][31][32][33][34] and also sector decomposition [25], to compute the integrals that arise, analytically and numerically. Indeed, all one-particle integrals, denoted formally by 1 above, have been evaluated with these methods [35][36][37]. The actual computation of integrated subtraction terms leads to a large number of rather elaborate phase space integrals, however these can to be computed once and for all.
In this paper, we compute the integral of the iterated singly-unresolved subtraction term, 2 dσ RR,A 12 m+2 , over the phase space of the unresolved partons. We find that the integrated approximate cross section can be written as the product of the Born cross section for the production of m partons, times a new insertion operator (in colour space), I (0) 12 . We use the method of MB representations, as developed in this context in ref. [37], to compute the integrals appearing in the various building blocks of the insertion operator. In several cases we find multi-dimensional MB integrals that are very difficult to compute fully analytically, and hence complete analytic expressions cannot be obtained at present. Nevertheless, in these cases direct numerical integration of the appropriate MB representations provides a fast and reliable way to obtain final results. We stress that for phenomenological applications, this is all that is required, hence, we make no severe effort to compute analytic expressions beyond those that are trivial to derive. As a numerical check, all integrals are evaluated using sector decomposition as well. Thus, each integral in this paper is obtained by two independent computations.
Since the paper is quite long and rather technical, readers mainly interested in understanding the general structure of the results or in some applications are advised to first read sections 2, 3.1 and 5. Section 2 sets the notation, and in section 3.1 we present the final expression for the integrated iterated singly-unresolved approximate cross section in eq. (3.13) and the new insertion operator I (0) 12 in eq. (3.14). These two equations are the main results of this paper. In section 5 we discuss some examples, specialising the general formulae to the case of two-and three-jet production processes. The explicit definitions of the integrated counterterms are then presented in sections 3.2 and 4. The technical details of computing the integrated subtraction terms are given in appendices.

Matrix elements
We consider processes with coloured particles (partons) in the final state, while the initialstate particles are colourless (typically electron-positron annihilation into hadrons). Any number of additional non-coloured final-state particles are allowed, too, but they will be suppressed in the notation. Resolved partons in the final state are labeled with letters chosen form the middle of the alphabet, i, j, k, l, . . . , while letters chosen form the end of the alphabet, r, s, t, . . . , denote unresolved final-state partons.
We adopt the colour-and spin-state notation of ref. [38]. In this notation the amplitude for a scattering process involving the final-state momenta {p}, |M({p}) , is an abstract vector in colour and spin space, and its normalisation is fixed such that the squared amplitude summed over colours and spins is This matrix element has the following expansion in the number of loops: where |M (0) denotes the tree-level contribution and the dots stand for higher-loop contributions, which are not used in this paper. Colour interactions at QCD vertices are represented by associating colour charges T i with the emission of a gluon from each parton i. In the colour-state notation, each vector |M is a colour-singlet state, and colour conservation implies where the sum over j extends over all the final state partons of the state vector |M (recall that we are considering processes where the initial state is colourless), and the equation is valid order by order in the loop expansion of eq. (2.2). Using the colour-state notation, we define the two-parton colour-correlated squared tree amplitudes as |M and similarly the four-parton colour-correlated squared tree amplitudes, We will also use the following ⊗ product notation to indicate the insertion of colour charge operators between M (0) | and |M (0) : (2.6) The colour-charge algebra for the product n (T i ) n (T k ) n ≡ T i ·T k is: Here C f i is the quadratic Casimir operator in the representation of particle i and we have C q ≡ C F = T R (N 2 c −1)/N c = (N 2 c −1)/(2N c ) in the fundamental and C g ≡ C A = 2 T R N c = N c in the adjoint representation, i.e. we are using the customary normalisation T R = 1/2.
We also use squared colour-charges with multiple indices, such as 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 using the flavour summation rules: odd/even number of quarks plus any number of gluons gives a quark/gluon, or explicitly for the relevant cases at NNLO: q + g = q , q +q = g , g + g = g , q + g + g = q , q + q +q = q , g + q +q = g , g + g + g = g . (2.8)

Cross sections
In this paper we shall need to use only cross sections of producing n partons at tree level with n = m, the Born cross section, and n = m + 2, the so-called doubly-real correction. We have 9) where N includes all QCD-independent factors and dφ n ({p}) is the d-dimensional phase space for n outgoing particles with momenta {p} ≡ {p 1 , . . . , p n } and total incoming momentum Q, dφ n (p 1 , . . . , p n ; The symbol {n} denotes summation over the different subprocesses and S {n} is the Bose symmetry factor for identical particles in the final state. Then the Born cross section and the doubly-real correction are simply The final result will also contain the phase space factor due to the integral over the (d − 3)-dimensional solid angle, which is included in the definition of the running coupling in the MS renormalisation scheme, 1 .
(2.12) 1 In the MS renormalisation scheme as often employed in the literature, the definition of the running coupling includes the factor Sǫ = (4π) ǫ e −ǫγ E . The two definitions lead to the same expressions in a computation at the NLO accuracy. At NNLO these lead to slightly different bookkeeping of the IR and UV poles at intermediate steps of the computation, but the physical cross section of infrared safe observables is the same. Our definition leads to somewhat simpler expressions at the NNLO level.

Momentum mappings and phase space factorisation
The iterated subtraction terms are written in terms of momenta obtained via various combinations of the basic collinear and soft mappings of ref. [1]: where both −→ may label either a collinear or soft mapping. (In general, both R and S are multiple indices.) As the above notation suggests, the final set of m momenta are denoted by tildes, while hats indicate the intermediate set of m + 1 momenta. In kinematic expressions where only the label of a momentum is displayed (we shall discuss several examples below), the tilde and/or hat is inherited by the label, and we write for instance i , ir and irs , where the latter two label single momenta. However, since these mappings affect only the momenta, but not the colour and flavour (apart from the flavour summation rules of eq. (2.8)), we shall omit the hat and tilde from flavour and colour indices.
We also use labels such as (ir) to denote a single momentum that is simply the sum of two momenta, p (ir) ≡ p i + p r .
Importantly, both the collinear and soft momentum mappings lead to an exact factorisation of phase space as follows: where here and below the bar denotes either a hat, when n = m + 1, or a tilde, when n = m. The one-particle factorized phase spaces can be written in the following form. For the collinear mapping we have [dp where the p (ir) appearing on the right hand side is understood to be written in terms of mapped momenta, that is p µ (ir) = (1 − α ir )p µ ir + α ir Q µ . For the soft mapping we find [dp where the momentum K is massive and K 2 = (1 − y rQ )Q 2 . As the notation above indicates, α ir and y rQ will become integration variables, hence their precise definitions are not important and will not be recalled here. (See ref. [1] for details.)

Kinematic variables
Three types of kinematic variables are used to write the iterated subtraction terms. The precise definitions of these were given in ref. [1]. Here we recall only those formulae that are needed for defining every expression precisely in the present paper.
• Two-particle invariants, such as s ir = 2p i ·p r , s kt r = 2 p kt · p r , or s iQ = 2p i ·Q and s ik ⊥ = 2p i ·k ⊥ . (2.18) The two-particle invariants scaled with Q 2 are denoted by y ij = s ij /Q 2 .
• Momentum fractions z i,r and z k , t for the splittings p ir → p i +p r and p kt → p k + p t , with z r,i = 1 − z i,r and z t , k = 1 − z k , t . Momentum fractions for three-particle splittings are denoted by z k,tr = y kQ y kQ + y tQ + y rQ , (2.20) with the expressions for z r,kt and z t,rk obtained by cyclic permutation. In the following, all momentum fractions will be integrated out, and so they will be expressed using the integration variables.
We also use extensively the uncontracted and contracted eikonal factors: As mentioned above, the sum of two momenta is often abbreviated with the two indices in parenthesis, e.g. p µ i + p µ r = p µ (ir) , which is used systematically in other occurrences, such as Finally, we express the integrated iterated subtraction terms as functions of the following (combinations of) invariants: In the centre-of-momentum frame (i.e. the rest frame of Q µ ), we find that is simply the scaled energy of parton i, while Y ij,Q = (1 − cos χ ij )/2, where χ ij is the angle between momenta p µ i and p µ j .
3 Integrating the iterated singly-unresolved approximate cross section

The integrated approximate cross section and insertion operator
We begin by recalling that the doubly-real emission cross section is defined precisely as in eq. (2.9), with n = m + 2. Then the iterated singly-unresolved approximate cross section times the jet function reads [1] dσ RR,A 12 m+2 where the notation ⊙ J m means that the jet function multiplying the different terms in the sum depends on different sets of momenta. Explicitly, the three terms in eq. (3.1) are given by All momentum mappings in eqs. (3.2)-(3.4) lead to the factorisation of the original m + 2 parton phase space into an m parton phase space times two one-parton phase space measures, as discussed originally in ref. [1], and recalled in section 2.3 above. Symbolically, we may write dφ m+2 ({p}) = dφ m ({ p })[dp 1,m ][dp 1,m+1 ] . (3.5) The jet function does not depend on the variables of the factorized one-parton measures, [dp 1,m ][dp 1,m+1 ], so we can compute the integral of eq. (3.1) over the phase space of the two unresolved partons, independently of the jet function J m , that we shall omit in the following. The result of the integration is a long expression of kinematics-dependent functions -each corresponding to a specific iteration of unresolved limits of the squared matrix elements -times colour factors, where the operation ⊗ means insertion of the colour charges between M m , see eq. (2.6). Three types of colour connections appear in eq. (3.6), and the functions on the right hand side -the "non flavour summed" (see below) integrated counterterms -take three different forms accordingly: where e.g. [X 1 ] f i ... represents a function that results in the integration of the counterterm X (0,0) 1 . The quadratic Casimir operators that appear in eqs. (3.7) and (3.8) are factored out to make the integrals [X 3 ] (i,k)(j,l) ) dimensionless in colour-space. As the notation implies, [X .. may also carry flavour dependence. Incidentally, we note that for every integrated counterterm in eq. (3.6), we consider everything inside the square brackets to be simply part of the function's name. In particular, the lower indices inside square brackets loose their meaning. Nevertheless, we choose to keep these in order to exhibit from which counterterm each function derives.
Eq. (3.6) is not yet in the form of an m-parton contribution times a factor. In order to obtain such a form, we must still perform summation over unresolved flavours, rewriting the symmetry factor of an m + 2-parton configuration to the symmetry factor of an m parton configuration. The complete details of this counting are not very difficult, but rather long, and are given in appendix A. As a result of these manipulations, we obtain functionsthe flavour summed integrated counterterms -denoted by X (0) (j,l)... f i ... , which are specific sums of the non flavour summed integrated subtraction terms, symbolically It is important to realise that the flavour indices on the left and right hand sides of the above equation need not match up. Indeed, the non flavour summed functions on the right hand side carry dependence on unresolved flavours, while the flavour summed functions on the left do not, by definition. Similar change in notation was introduced in the dipole subtraction scheme [38], where we find the definitions (see eqs. (7.21) and (7.22)) where the functions V i (ǫ) on the left hand side represent flavour summed integrated counterterms, while the V ij (ǫ) functions on the right hand side are not flavour summed. In the present paper, due to the extra complications of an NNLO subtraction scheme, we believe that it is helpful to make a sharper notational distinction between flavour summed and non flavour summed integrated counterterms.
After summation over unresolved flavours in eq. (3.6), we find that the final result can be written in the form where the insertion operator (in colour space) has three contributions according to the possible colour structures in eqs. (3.7)-(3.9): 14) with f i denoting flavours, and C q = C F ≡ T 2 q , C g = C A ≡ T 2 g as in eq. (2.7). Eqs. (3.13) and (3.14) are the main results of this paper. In the following, we shall define and compute each term appearing in eq. (3.14). First, in terms of flavour summed integrated counterterms discussed above, we get: .

(3.19)
On the right hand sides of these equations the flavour dependent functions are the flavour summed integrated counterterms discussed above. They depend on the kinematics through variables of the type x i and Y ij,Q . The latter dependence derives from integrating an eikonal factor which is always multiplied with a colour-connected squared matrix element. In order to make the results more transparent, we hid the arguments of the functions, but kept the relation to the colour-connected matrix elements, shown as upper indices.

Flavour summed integrated counterterms
Here we list the flavour summed integrated counterterms appearing on the right hand sides of eqs. (3.15)- (3.19), written in terms of the integrated subtraction terms.
3. Soft-collinear-triple collinear-double soft: i.e. it is independent of the flavour f .
6. Flavour-independent soft-collinear-double soft: kt ] (j,l) . Analytic expressions for the expansion coefficients have been obtained to O(ǫ −2 ) accuracy in all cases, and we present these below. However, in the case of the single and double poles as well as the finite terms, we encountered several instances where obtaining complete analytic expressions was not feasible. This being the case, we made no severe effort to derive analytic expressions beyond those presented here. Higher order coefficients in the Laurent expansions will be given numerically in the form of a computer code elsewhere. In the following results we use d ′ 0 = D ′ 0 + d ′ 1 ǫ (see appendix B and especially eq. (B.8)) and the abbreviations Collinear-type terms: 1. Collinear-triple collinear: 2. Collinear-double collinear: (3.43) 3. Collinear-soft-collinear: (3.44) When e.g. i = j, the functions C kt CS different (where i = j, l is understood) but up to this order, the functional dependence on the variables is the same.
4. Collinear-triple collinear-soft-collinear: (3.45) 5. Collinear-double collinear-soft-collinear: (3.46) As implied by the notation, this function depends on both x i and Y ik,Q . However, the dependence on the latter vanishes up to this order.
4. Flavour-dependent soft-collinear-double soft: This function is independent of the kinematics.
6. Flavour-independent soft-collinear-double soft: , we obtain the following explicit expressions for the kinematics dependent functions appearing on the right hand side of eq. (3.14). For C The function C 12,f i f k simply vanishes up to this order in the ǫ-expansion The two-parton colour-correlated soft function, S (0),(j,l) 12 , is if i is distinct from both j and l, while for e.g. i = j we have , reads if only two indices, say i and k are distinct, while for three distinct indices. Furthermore, the above expression is valid in the case of restricted kinematics of three hard partons.

Integrated counterterms
In this section we list the explicit definitions of the functions that appear on the right hand side of an equation among eqs. (3.20)-(3.39).

Integrated collinear-type counterterms
1. Collinear-triple collinear: f k ftfr are the strongly-ordered three-parton splitting kernels averaged over the spin states of the parent parton (see appendix D.1 and especially eq. (D.27)). The subtraction terms contain the spin-dependent kernels, that together with the corresponding kinematic variables can be found in ref. [1]. In appendix D.1 we prove that the integrals of the spin-dependent kernels give the same contribution as those of the spin-averaged ones, therefore, we can use the latter when integrating the subtraction terms.
The f (α 0 , α, d) functions, defined in eq. (B.6), represent simple modifications to the original subtraction scheme of ref. [1]. As discussed in appendix B in detail, these modifications do not destroy the cancellation of singularities, but serve improved numerical control, efficiency and stability, and result in simpler, m independent, integrated counterterms. The rest of the counterterms are modified similarly by appropriate f functions, and below we shall simply include these factors without further comment.
The integrated counterterm is computed in appendix D. 3. In terms of the functions I The various coefficients read: qg,j , and finally 2. Collinear-double collinear: The integrated counterterm is computed in appendix D. 4. In terms of the function I with coefficients given in eq. (4.4).
3. Collinear-soft-collinear: where S j l ( r ) is the eikonal factor defined in eq. (2.21) and P The integrated counterterm is computed in appendix D.5. In terms of the function I with coefficients given in eq. (4.4). The integral I C with full kinematic dependence, as written above, first appears in computing NNLO corrections to four-jet production.
4. Collinear-triple collinear-soft-collinear: (4.10) This integrated counterterm is also computed in appendix D.5. In terms of the function I C of eq. (D.70) we find with coefficients given in eq. (4.4).

5.
Collinear-double collinear-soft-collinear: This integrated counterterm is also computed in appendix D.5. In terms of the function I C of eq. (D.70) we find with coefficients given in eq. (4.4).
The integrated counterterm is computed in appendix D.6. In terms of the functions I with the flavour-dependent constants b (0) f k ft given in eq. (4.3). Of course, only f k + f t = g gives a non-vanishing result (see eq. (3.25)). In obtaining the form (4.16) of the integrated counterterm, we exploited that the integrals of the two expressions in eq. (D.17) are equal.
The integrated counterterm is computed in appendix D. 6. In terms of the functions I
2. Soft-soft-collinear: [dp with coefficients given in eq. (4.4). Note that j and l are always distinct.
3. Soft-triple collinear-soft-collinear: This integrated counterterm is computed in appendix E.1. In terms of the function I (4.24) 2 We note a harmless misprint in the definition of the subtraction term S t C irt S (0,0) rt in eq. (7.38) of ref. [1], where in the last term of the square bracket The integrated counterterm is computed in appendix E.2. In terms of the functions I The integrated counterterm is computed in appendix E.2. In terms of the functions I Note that j and l are always distinct.

(4.28)
This integrated counterterm is computed in appendix E.2. In terms of the function I (10) S of eq. (E.44) we find 7. Soft-double soft: In this paper we do not discuss the case when i, j, k, and l are all distinct, which first appears in computing NNLO corrections to four-jet production. For the specific cases of two and three hard partons in the final state, we compute the integrated counterterms in appendix E.2. In terms of the functions I and

Integrated soft-collinear-type counterterms
1. Soft-collinear-triple collinear: where we used T 2 krt = T 2 kr because the soft parton t can only be gluon, f t = g. For the same reason T 2 kt /T 2 krt = T 2 k /T 2 kr , to be used in eq. (4.39).

Soft-collinear-soft-collinear:
(4.35) 3. Soft-collinear-triple collinear-double soft: The integrals over [dp kt ] (as seen by the simple exchange of indices k ↔ r) and the distinction of these functions is purely formal, and serves bookkeeping purposes only.
4. Soft-collinear-double soft: where the first equation above defines the flavour-dependent soft-collinear-double soft function, while the second equation gives the flavour-independent one. Again, the integrals over [dp kt ] (j,l) (as seen by the simple exchange of indices k ↔ r) and the distinction of these functions is formal, and serves bookkeeping purposes.
5. Soft-collinear-soft-collinear-double soft: with coefficients given in eq. (4.4). This ends the definition of the integrals of the iterated subtraction terms, that can be used to construct the insertion operator I

Insertion operator for two-and three-jet production
In this section, we present illustrative numerical results for the insertion operator I (0) 12 of eq. (3.14), at various specific phase space points, for processes with at most three hard partons in the final state.
While it is true that the most general collinear and/or soft configuration at NNLO accuracy involves four hard partons, if these are all in the final state, which is the subject of this paper, then the two-loop amplitudes needed in the doubly-virtual correction have at least six massless (or four massless and one massive) external legs. Results for such amplitudes are not foreseen in the near future, therefore, we restrict our discussion to computing NNLO corrections to two-and three-jet quantities.
In an explicit computation of a jet cross section at NNLO, we require the expansion coefficients (up to and including the finite part) of the Laurent series in ǫ of the insertion operator. In general, these expansion coefficients are functions of various kinematic variables (and also parameters such as α 0 , y 0 , d 0 and d ′ 0 ) which depend on the particular phase space point. One may either attempt to compute these functions analytically, or numerically. The former is important as a matter of principle only. For practical purposes (from the point of phenomenology) the latter is sufficient. Indeed, the higher order expansion coefficients (starting form O(ǫ −2 )) of the results we will present were obtained numerically.

Two-jet production
Let us consider the process e + e − → 2 jets. The corresponding squared matrix element at tree level is |M i.e. the quark carries label 1 and the antiquark label 2. Both the colour algebra and kinematics are trivial. Colour conservation implies Hence, the insertion operator is a scalar in colour space. On the other hand, momentum conservation requires that the two final state momenta are back-to-back, i.e. in a properly oriented frame we have The insertion operator eq. (3.14) becomes with all arguments being equal to 1. Substituting the Laurent expansions of the kinematic functions, we obtain With this decomposition the Abelian case is obtained by setting C F = 1, x = 0, y = 1. We can compute the two leading terms in the ǫ expansion analytically: The rest of the expansion coefficients are computed numerically. We present the results in table 1. To obtain these numbers, we used the specific values of α 0 = y 0 = 1, Finally, we show the value of the complete insertion operator for the case of QCD with n f = 5 light flavours: In the above equation, the coefficients of 1/ǫ 4 and 1/ǫ 3 were computed by evaluating the appropriate analytic expressions. However, we have also computed these term with the same numerical algorithms that we used for computing the higher order expansion coefficients. It is then instructive to compare this numerical result to the analytic one. We find: Comparing the exact and numerical results, we see first of all that the two match up to the uncertainty of the numerical computation, and second, that the error estimate on the next-to-leading pole is very conservative.

Three-jet production
Next, let us consider the process e + e − → 3 jets. The corresponding squared matrix element at tree level is |M i.e. the quark carries label 1, the antiquark label 2 and the gluon carries label 3. The colour algebra is again trivial. Colour conservation implies Thus, the insertion operator is again a scalar in colour space. On the other hand, the kinematics is no longer trivial, since the relative orientation of the three final state momenta are not fully fixed by momentum conservation.
(Note that the insertion operator is independent of the overall event orientation with respect to the beam.) Since the three-particle phase space in d = 4 dimensions is 5-dimensional, but three of the independent variables just correspond to the three Euler angles needed to specify the overall orientation, we find that out of the six kinematic variables only two are independent. Nevertheless, we choose not to fix the independent ones below, in order to better exhibit the structure of the insertion operator.
The insertion operator eq. (3.14) becomes Substituting the Laurent expansions of the kinematic functions, we obtain . (5.14) We can compute the leading two terms in the ǫ expansion analytically: In order to obtain eq. (5.15), we used ln Y ij = ln y ij − ln x i − ln x j . The rest of the expansion coefficients are computed numerically. For purposes of demonstration, below we present numerical results in three specific phase space points. Since the result is insensitive to overall orientation, we will always choose the event to lie in the x − y plane, with p µ 1 pointing in the positive y direction. In all cases we use the specific values of α 0 = y 0 = 1, Symmetric point. First, we consider the maximally symmetric configuration which leads to the following values for the kinematic invariants (i, j = 1, 2, 3 and i = j): The coefficients of the Laurent expansion of the insertion operator in the symmetric phase space point are shown in table 2. Finally, we show the value of the complete insertion operator in the symmetric phase space point for the case of QCD with n f = 5 light flavours: In this equation the coefficients of 1/ǫ 4 and 1/ǫ 3 were obtained by evaluating the appropriate analytic expressions. As in the case of two-jet production, it is again instructive to compare the exact result with one obtained by numerical computation. We find: As before, the results match within the uncertainty of the numerical calculation, and the error estimate on the next-to-leading pole is again seen to be very conservative. We will reach similar conclusions in other examples as well.
Collinear point. Next, we choose a configuration where (in the rest frame of Q µ ) we have a hierarchy of angles such that i.e. where momenta p µ 2 and p µ 2 are close to being collinear. Specifically, we set The coefficients of the Laurent expansion of the insertion operator in the collinear phase space point are shown in table 3. Finally, we show the value of the complete insertion operator in the collinear phase space point for the case of QCD with n f = 5 light flavours: We present next the comparison of the exact 1/ǫ 4 and 1/ǫ 3 pole coefficients given above  Table 3. Coefficients of the Laurent expansion of the I 12,3j functions appearing in the insertion operator I with ones computed numerically. We find: We note that the two results match up to the uncertainty of the numerical computation, and as in previous examples, the error estimate on the next-to-leading pole is seen to be very conservative.
Soft point. Finally, we consider a configuration where (in the rest frame of Q µ ) we have a hierarchy of energies such that i.e. where momentum p µ 3 is close to being soft. Specifically, we set 12 (p 1 , p 2 , p 3 ; ǫ) for three-jet production in the soft phase space point. The numbers for i = −4, −3 are obtained by evaluating the appropriate analytic expressions. We used the parameters Finally, we show the value of the complete insertion operator in the soft phase space point for the case of QCD with n f = 5 light flavours: We finish by comparing the exact coefficients of 1/ǫ 4 and 1/ǫ 3 that appear above with the values obtained via numerical computation. We find: Our conclusions are identical to those in the symmetric and collinear phase space point: the values match up to the numerical uncertainty and the error estimate on the next-to-leading pole is again shown to be very conservative.
We finish by briefly commenting on the size of numerical uncertainties. The uncertainties relevant for phenomenology are those associated with the complete I (0) 12 insertion operator, in various phase space points. However, the requirements in terms of precision are different for the pole coefficients and the finite part.
On the one hand, the pole coefficients are only relevant for establishing the cancellation of all ǫ-poles between the doubly-virtual cross section and various integrated subtraction terms. As stressed earlier, our subtraction scheme is fully local, hence this cancellation can be checked point by point in phase space for any specific process. From a practical point of view, it clearly suffices to demonstrate pole cancellation in a relatively small number of phase space points, thus the pole coefficients of I (0) 12 have to be computed as precisely as feasible in a small set of points only. Because of this, the runtime of numerical integration is not an issue, and increased precision may be obtained simply by adjusting the parameters of the numerical integration to include more sampling points for each integral.
On the other hand, the precision requirement on the finite part of the insertion operator is essentially set by the relative uncertainty associated with the numerical phase space integration of the doubly-virtual contribution. This is not expected to be below the per mille level, hence from a practical point of view it is pointless to evaluate the finite term of the insertion operator with a precision much greater than this. In all cases discussed above, the relative uncertainty of the finite part of the insertion operator is at the per mille level already.

Conclusions
In this work, we have performed the integration of the iterated singly-unresolved approximate cross section of the NNLO subtraction scheme of refs. [1][2][3]. The final result can be written as the product (in colour space) of the Born cross section times a newly defined insertion operator, I 12 . The insertion operator depends on the colours, flavours and momenta of the final-state partons, and is an elaborate sum of many different terms, each corresponding to the integrated form of a specific iterated singly-unresolved subtraction term of ref. [1].
We have also explicitly evaluated all integrated subtraction terms which are necessary to assemble the insertion operator for processes involving at most three hard partons in the final state. The knowledge of these integrals (i.e. their Laurent expansions in ǫ to O(ǫ) accuracy) is necessary in order to make the subtraction scheme an effective tool, and we have computed them once and for all.
We have achieved this task by deriving Mellin-Barnes integral representations for all integrals under consideration. In principle, it is possible to evaluate all MB integrals via the residuum theorem, and in a subsequent step to obtain fully analytic expressions by performing the summation of nested sums over series of residua. However in practice, we have encountered several cases of higher order expansion coefficients, where the summation cannot be performed analytically with present methods. Therefore, in this paper, we have concentrated on the direct numerical evaluation of the MB integrals in the complex plane. All MB representations for both the numerical and, if available, the analytic expressions have been checked by an independent evaluation of the integrals using sector decomposition as in ref. [35]. We have found that all integrals contributing to the insertion operator are smooth functions of their variables (in the colloquial sense). For practical applications, this means that all integrals (in particular the finite in ǫ contributions) can be given either in terms of interpolating tables or simple fitting functions, which can be computed once and for all. We leave this step for later work. Finally we want to stress again that the tables we have shown here are for demonstration purposes only, and obtaining high resolution interpolating tables needed for the computation of an actual cross section is straightforward. Increasing the accuracy of each entry is feasible as well, the best way of doing this is under investigation.
The integrals discussed in this paper appear when integrating the subtraction terms that regularise the doubly-real NNLO correction to the jet cross section, see ref. [2]. The final step in finishing the definition of the subtraction scheme is the computation of the integrated counterterms corresponding to the the doubly-unresolved approximate cross section (those labeled by A 2 in ref. [1]). In that case, the analytic structure of the integrals is essentially the same as those studied in this paper, though a few are admittedly somewhat more cumbersome. Nevertheless, we are confident that the techniques of the present paper will also be applicable to the computation of these remaining contributions.

A Summation over unresolved flavours
In this appendix, we discuss how to perform the summation over unobserved flavours in eq. (3.6). It turns out that we only need to consider three different cases explicitly, corresponding to sums involving two, three and four partons. All specific results are then easily obtained by appropriate substitutions.

A.1 Generic flavour sums
Consider an m-parton configuration with m f quarks of flavour f , mf antiquarks of flavour f and m g gluons. From this configuration we can obtain an (m + 2)-parton configuration in the following ways.
1. Increasing the number of gluons by two, (A.1) 2. Increasing the number of quarks and antiquarks of flavour f by one each, 3. Increasing the number of quarks and antiquarks by two each. The flavour of the two quarks may or may not be identical. We will refer to these two cases respectively as the 'equal flavour' (e.f.) and 'unequal flavour' (u.f.) configurations, where f and f ′ are understood to be different quark flavours. This case is relevant only for the doubly-collinear-type configuration, i.e. the sum involving four partons.
The ratios of Bose symmetry factors for identical final state particles in the various cases are , and finally rt ] (j,l) . In the latter case, we see a situation where no flavour index is displayed, since both r and t are constrained to be gluons. In what follows, we will discuss the most general case, when both flavour indices are explicit.
Such terms necessarily appear in eq. (3.6) under a double sum: In this configuration, we go from m to (m + 2) partons as in eqs. (A.1) and (A.2). Then, since we are considering iterated singly-unresolved terms, both indices must correspond to unresolved partons, which means that they are either both gluons or a quark-antiquark pair. Now, decomposing the summation over t and k into sums in which the flavour (including specific quark flavour) of each index is fixed, we find where f ′ stands for the explicit summation over specific quark flavours. Performing the summation over t and k simply amounts to counting the number of ways in which we can assign the proper flavours to k and t in the appropriate (m + 2)-parton configuration: where #(f t ) m+2 denotes the number of partons of flavour f t in the (m + 2)-parton configuration, while #(f k ; k = t) m+2 is the number of partons, different form t, of flavour f k in the (m + 2)-parton configuration. Note that t and k are assumed to be distinguishable, which is the generic case. Clearly we have , The qq do not depend on the specific quark flavour (as implied by the notation), and hence the summation f ′ in eq. (A.8) may be performed, yielding the factor of n f .
Defining the flavour summed counterterm as we obtain : where the second line follows, since [X (j,l) f i fr . As before, we will discuss the most general case, when all flavour indices are explicit.
These terms always appear in eq. (3.6) under a triple sum: In this configuration, we again go from m to (m + 2) partons as in eqs. (A.1) and (A.2). Then, we decompose the summation over t, k and r into sums in which the flavour (including specific quark flavour) of each index is fixed. We obtain Next, we use the flavour summation rules to rewrite the summation over the unobserved indices k, t and r in the (m + 2)-parton configurations into a sum over a single index ktr in the m-parton configuration. We have where the notation is the same as in eq. (A.9) and in particular #(f ktr ) m is the number of partons with flavour f ktr in the m-parton configuration. Again, t, k and r are assumed to be distinguishable, which is the generic case. Then we have . .
. q ′q′ q , which implies that the summations f ′ in eq. (A.15) may be performed, yielding the factors of (n f − 1) and n f . Let us define the flavour summed counterterms as follows: Then we find These terms always appear in eq. (3.6) under a four-fold sum: In this configuration, we go from m to (m + 2) partons as in eqs. (A.1)-(A.3), i.e. case 3 must also be considered. Decomposing the summation over t, k, i and r into sums in which the flavour (including specific quark flavour) of each index is fixed, we find Four-index subtraction terms only arise in conjunction with the double collinear limit and are always completely independent of the specific quark flavours. We have used these facts to write eq. (A.22) in the above form. First, since the pairs of indices k, t and i, r will always correspond to true singly-unresolved collinear limits, we have discarded all terms where this is not the case. In effect, we have dropped all terms where both k and t or both ktir ] qgqq and so on. We have used this fact in writing the equation, hence, in (A. 22), f and f ′ are not necessarily distinct flavours. Then using the flavour summation rules, we can rewrite the summation over the unobserved indices k, t, i and r in the (m+2)parton configurations into sums over indices kt and ir in the m-parton configuration. We have where the notation is the same as in eqs. (A.9) and (A.16). We assume that in general k, t, i and r are all distinguishable. Then we find , , #(g) m #(g; ir = kt ) m = (m f ′ + 1)(mf′ + 1) .

(A.24)
In case 3 we must remember that the counting is slightly different for the 'equal flavour' and 'unequal flavour' contributions even when the counterterms are the same. We have , e.f.
, u.f.  We remind the reader that eq. (A.26) was derived by using that the counterterms are independent of specific quark flavours (as the notation implies), and further that the 'equal flavour' and 'unequal flavour' subtraction terms are equal. Then the sums f and f ′ in eq. (A.22) may be performed, and we obtain the factors of n f and n 2 f as shown.
Finally, we define the flavour summed counterterms as  • In eq. (3.6), most counterterms appear with some explicit overall factor, which must be included in the final result. E.g. for the collinear-triple collinear counterterm, this factor is 1/2.
• In certain cases, the ordering of some flavour indices may be meaningless, due to a symmetry of the integrated counterterms. E.g. in the collinear-triple collinear case, the integrated counterterm is symmetric in the first two indices, [C kt C  ktr ] qq ′q′ . Recall that this is already taken into account in eq. (A.28).
With these in mind, we easily find e.g.
for the collinear-triple collinear flavour summed counterterms. The rest of the results in section 3.2 are obtained similarly.

B Modified doubly-real subtraction terms
We outline a simple modification to the NNLO subtraction scheme presented in refs. [1,2]. Parts of these modifications were presented previously: those relevant to the singlyunresolved approximate cross section dσ RR,A 1 m+2 appearing in eq. (1.3), and to the approximate cross sections in eq. (1.4), were presented in ref. [35]. In this appendix we describe the modification of the iterated singly-unresolved approximate cross section dσ RR,A 12 m+2 , which appears in eq. (1.3).
Recall that the iterated singly-unresolved approximate cross section can be written symbolically as dσ where the iterated singly-unresolved approximation A 12 |M which are iterated in various combinations to produce appropriate mappings of m + 2 → m momenta. As discussed in section 2.3, all such mappings lead to an exact factorisation of the m + 2 particle phase space, symbolically written as The exact form of the factorized phase spaces [dp 1,n ] (n = m, m + 1) is given in eqs. (2.16) and (2.17), but their only feature which is relevant presently is that they carry a dependence on the number of partons, n, of the form [dp The subtraction terms, as originally defined in ref. [1] do not depend on the number of hard partons, thus the m-dependence of the factorized phase space measures is carried over to the integrated counterterms, where furthermore this dependence enters in a rather cumbersome way (see e.g. eqs. (A.9) and (A.10) of ref. [3]).
Thus, as in ref. [35], we reshuffle the m-dependence of the integrated counterterms into the subtraction terms themselves, where it appears in a very straightforward and harmless way, through factors of (1 − α) and/or (1 − y) raised to m-dependent powers. For easier reference, we gather the results in table 5, where together with the subtraction terms, we give the momentum mappings used to define the term(s) and the function which multiplies the original counterterm to produce the modified one. The f functions appearing in table 5 are defined as The pattern of modifications is hopefully clear: if the factorized phase space appropriate to a given subtraction term carries m-dependence through factors of (1 − α) and/or (1 − y) respectively, it is multiplied by a factor/factors of f (α 0 , α, d(m, ǫ)) and/or f (y 0 , y, d ′ (m, ǫ)). We emphasise that the form of the exponents d(m, ǫ) and d ′ (m, ǫ), is actually fixed by the prescription in ref. [35] (see eqs. (3.2), (3.12) and (3.13) in particular) and the requirement that the modified subtraction terms should still correctly regularise all kinematic singularities. In fact, we must have where d 0 and d ′ 0 are the same constants which appear in eqs. (3.2), (3.12) and (3.13) of ref. [35], i.e.
where D 0 , D ′ 0 ≥ 2 are integers, while d 1 , d ′ 1 are real. Also, the parameters α 0 and y 0 must have the same values for all subtraction terms, including the singly-unresolved ones of ref. [35].
Finally, we note that the modifications introduced above do not spoil any of the cancellations which take place among the original subtraction terms, hence the modified counterterms are still a correct regulator of all kinematic singularities. This is not particularly hard to check explicitly, and is actually a manifestation of the fact that the various momentum mappings obey several conditions in soft and/or collinear limits. As these were discussed in ref. [1], we do not go into further details here.

C Basic collinear, soft and soft-collinear functions
Certain basic functions appear repeatedly during computations in this paper. They all arise as integrals of various simple factors over factorized collinear or soft phase space measures. Below we define and give explicit integral representations of these functions. Some, notably I, J and K, have been considered previously in ref. [35], albeit in somewhat more general forms. For completeness, we present these here as well, although only in the special cases used in this paper.

Subtraction term
Momentum mapping Function

Iterated soft counterterms
Subtraction term Momentum mapping Function

Subtraction term
Momentum mapping Function Table 5. The modified iterated singly-unresolved subtraction terms are obtained from the original counterterms (first column) by multiplication with an appropriate function (last column). Also shown are the momentum mappings used to define the subtraction terms (middle column). The f (z 0 , z, p) function is defined in eq. (B.6) while d(m, ǫ) and d ′ (m, ǫ) are defined in eq. (B.7).

C.1 Collinear functions
When computing the integral of the azimuthally averaged Altarelli-Parisi splitting functions over the factorized collinear phase space, the following function arises: which, as indicated, is simply the I function of [35], for the special values of parameters κ = δ = 0 and g I = 1 (for their meaning see ref. [35]). The integral is over [dp (ir) 1,m+1 ], i.e. the factorized measure obtained when going from m + 2 to m + 1 partons via the collinear mapping of ref. [1], which can explicitly be written as [dp The general collinear function was first computed analytically in ref. [36]. For convenience, we present the integral representation for the specific case, used in this paper: where the factors correspond to the collinear pole y ir and momentum fraction z r,i , respectively (with x ≡ x ir , α ≡ α ir and v ≡ 1 − v ir ). Among the iterated subtraction terms considered in this paper, we find two other basic integrals over the collinear phase space measure (C.2). One of these is the Lorentz tensor (C.5) with kinematic dependence only on p µ ir and Q µ . The transverse momentum is defined to be orthogonal to both of these, p ir ·k ⊥,i,r = Q·k ⊥,i,r = 0 and contraction with g µν replaces the fraction in the last factor with four. The most general Lorentz structure that obeys these conditions is where the integral in the second line is clearly just The other one is The definition of z i,r [1] implies y iQ = z i,r (y iQ + y rQ ), and from the definition of p ir we have y iQ + y rQ = 2α ir + (1 − α ir )x ir . The momentum fraction can be expressed with the integration variables as in eq. (C.4), therefore, (C.9)

C.2 Soft functions
When integrating the eikonal factor s ik /(s ir s kr ) over the factorized soft phase space, we encounter the following function: which, as shown, is simply the J function of ref. [35], for κ = 0. The integral is over [dp where (here and below) the dots stand for vanishing components, while the notation 'angles' indicates the dependence of p µ r on the d−4 angular variables that can be trivially integrated in all relevant cases. In terms of the scaled energy-like variable y rQ and the angular variables ϑ, ϕ and η, the two-particle phase space dφ 2 (p r , K; Q) is Often the integrand does not depend on all angles and we can integrate out η, and ϕ, (C.14) The J 1 function of eq. (C.10) was first computed analytically in ref. [36]. We recall that it is conveniently evaluated in the frame eq. (C.11), with the orientation fixed by The precise definitions of p µ i and p µ k via the soft mapping [1] imply and expressing all two-particle invariants with integration variables, we find [36] and the function Ω jl denotes the angular integral Presently we need the special case β 1 = β 2 = 1, which we call the 'massless' angular integral. The result of this angular integration is well-known [40] and is proportional to a hypergeometric function. Finally, using some hypergeometric identities and a onedimensional integral representation of the hypergeometric function, we derive the following integral representation for J 1 , to be used in this paper: (C.20) In some cases the eikonal factor involves three momenta as in Then the soft integral cannot be expressed with the soft function J any longer. In eq. (C.22) β ( i k ) is the velocity of the momentum p µ i + p µ k in the centre of mass frame. We evaluate J (1m) in the frame (C.11), with orientation specified by . . . , sin χ, cos χ) .
Using the precise definition of p µ j (j = i, k, l), we find s il + s kl = (1 − y rQ )(s i l + s k l ) , s ir + s kr = s i r + s k r , s lr = s l r .

(C.24)
Then we express all relevant two-particle invariants with integration variables and obtain (C.25) where we used the soft phase space eq. (C.12) with η integrated out (see eq. (C.13)). In eq. 26) and the function Ω jl is defined in eq. (C. 19) above. Now we need the special case β 1 ≡ β < 1 and β 2 = 1, which we call the 'one-mass' angular integral. The evaluation of this integral will be discussed elsewhere, and here we simply indicate that the result is proportional to an Appell function of the first kind. Finally, using a one-dimensional integral representation of the Appell function, we obtain an integral representation for J (1m) , similar to eq. (C.17): (C.27) Setting β = 1, we see that J (1m) (Y, 1; ǫ, y 0 , d ′ 0 ) = J 1 (Y ; ǫ, y 0 , d ′ 0 ), as expected.

C.3 Soft-collinear functions
The following function arises when integrating the collinear limit of the eikonal factor, 2z i,r /(s ir z r,i ), over the factorized soft phase space (recall that z i,r /z r,i = s iQ /s rQ ): ǫ)) . (C.28) As indicated, this is just the K function that was defined and computed in ref. [35], for κ = 0. The integral is over [dp which follow from the precise definition of p µ i via the soft mapping. We find cos ϑ) , which can also be written as where we made the substitution cos ϑ → 1 − 2z. Notice that the soft-collinear function K is independent of the kinematics. In some cases, the collinear limit of the eikonal factor involves three momenta as in Then the soft-collinear integral (C.33) cannot be expressed with the soft-collinear function K any longer. In eq. (C.33), β ( i r ) is the velocity of the momentum p µ i + p µ r in the centre of mass frame. Using the definition of the mapped momenta [1], we have Then we evaluate K (1m) in the frame given in eqs. (C.11) and (C.23) (with the trivial replacement k → r). Expressing all two-particle invariants with integration variables, we find the following integral representation for the 'one-mass' soft-collinear integral: (C.35) where β ≡ β ( i r ) . We make the substitution cos ϑ → 1 − 2z to obtain the final form of the integral representation used in this paper: (C.36) For β = 1, we recover the soft-collinear integral,

D Integrating the collinear-type counterterms
In this appendix, we discuss the integration of the collinear-type counterterms of section 4.1.
In particular, we define and give an explicit integral representation of all I (i) C functions (i = 1, . . . , 9).

D.1 Treatment of azimuthal correlations
Collinear subtraction terms contain azimuthal correlations if the factorisation formula in the corresponding collinear limit involves gluon splitting kernels. For the two-parton splitting ij → i + j, these azimuthal correlations involve a transverse momentum k µ ⊥,i,j that we always chose to be orthogonal to the parent momentum p µ ij . This condition is sufficient to prove that the integral of the spin-dependent and spin-averaged splitting kernels over the factorized phase space of the unresolved parton are the same [38]. Therefore, one can always substitute the spin-dependent splitting kernelsP , as done in section 4. In the case of strongly-ordered three-parton splittings, one cannot directly use the same argument. The splitting kernels in the integral of the collinear-triple collinear subtraction, C kt C (0,0) ktr depend on the transverse momentum k µ ⊥,k,t in two ways. One is when the Lorentz index of the transverse momentum coincides with that of the parent gluon as in the explicit k µ ⊥,k,t k ν ⊥,k,t /k 2 ⊥,k,t terms in the gluon splitting kernels: and µ|P s.o. (0) g k gtgr (z k,t , z t,k , k ⊥,k,t , z kt, r , z r , kt , k ⊥, kt, r ; ǫ)|ν = 4C 2 A − g µν z r , kt z kt, r + z kt, r This transverse momentum is not orthogonal to p µ ktr , as defined originally in ref. [1]. Nevertheless, when integrating these subtraction terms, we can still substitute the spin-dependent splitting kernels with the spin-averaged ones, as we now show.
Recall that the strongly-ordered three-parton splitting kernel appears in the collineartriple collinear subtraction term in the form and the bra-ket expression above has the following precise meaning: where d µµ ′ ( p ktr , n 1 ) and d νν ′ ( p ktr , n 2 ) are gluon polarisation tensors with (time-like) gauge vectors, n 1 and n 2 . Hence, the object to integrate is not simply µ|P that is clearly orthogonal to p µ ktr (because of the presence of the polarisation tensors). Thus, by the usual arguments, the azimuthal correlations present in C kt C (0,0) ktr vanish after integration over the phase space of the unresolved parton. However, we must still be careful to compute the average over the polarisations correctly. When k ⊥ · p = 0, we have where the dots stand for terms proportional to p µ or p ν which vanish after contraction with the matrix element, by gauge invariance. Thus we find ( n is a further time-like gauge vector) Eq. (D.7) shows that the advantage of having k ⊥ · p = 0 actually lies in the trivial azimuthal averaging. However, this can also be arranged if k ⊥ · p = 0. For example, choosing gauge vectors n µ 1 and n µ 2 such that n µ 1 + n µ 2 ∝ k µ ⊥ (and of course n 2 1 = n 2 2 = 0), we find that the above average just reduces to −1/[2(1 − ǫ)], which is the usual result. Since k 2 ⊥ < 0, such nonzero n µ 1 and n µ 2 always exist. Indeed, in any Lorentz frame we have the parametrisation k µ ⊥ = K(1, β v), where v 2 = 1 and β 2 = 1. Then setting e.g.
A shorter, though less transparent proof is to observe that if we change this single troublesome term to wherek µ ⊥,k,t = k µ ⊥,k,t − k ⊥,k,t · p ktr p ktr ·Q Q µ , (D.10) then obviouslyk ⊥,k,t · p ktr = 0, and all the usual arguments apply. What is not immediately obvious, is that this modification does not ruin any of the delicate cancellations between the various counterterms in any IR limit, and hence is allowed. This can be established as follows: in the centre of mass frame of Q, the replacement of eq. (D.9) simply amounts to with some ∆. Then, continuing to choose the gauge vectors as in eq. (D.8), we have i.e. only the normalisation of the gauge vectors is changed by the replacement in eq. (D.9). However, this implies that eq. (D.6) is actually unchanged. Indeed, recalling that k µ ⊥ = n µ 1 + n µ 2 , we find (D. 13) and the last expression is clearly seen to be invariant under (independent) rescalings of n 1 and n 2 . Hence, the replacement in eq. (D.9) is completely harmless. (However, notice that the proof requires the specific choice of gauge vectors as in eq. (D.8).) The other occurrence of k ⊥,k,t in the strongly-ordered kernels is in the ratio s 2 also present in the quark splitting kernels. Examining the explicit forms of the stronglyordered splitting kernels, we find that this ratio always appears in the form p µ r p ν r s kt r I µν ( p kt , Q) , (D.14) where the integral I µν is defined in eq. (C.5) and computed in eq. (C.6). Contracting the latter with p µ r p ν r /s kt r , we obtain Observing that s r Q s ktQ = z r , kt z kt, r , (D. 16) we find that when integrating the strongly-ordered splitting kernels over the factorized phase space, the integrals of are equal, so we can substitute the former with the latter, which we implement in the next subsection, where we give the spin-averaged splitting kernels explicitly.

D.2 Explicit forms of the spin-averaged splitting kernels
For the sake of completeness we present the explicit expressions for the spin-averaged splitting kernels used in the integrals of the subtraction terms. The azimuthally averaged Altarelli-Parisi splitting kernels are well known: In our convention the ordering of the labels on the splitting-kernels is usually meaningless, but in eq. (D.20) z refers to the momentum fraction of the second label. In other words P f k ftfr the ordering matters, too. As a result, the same triple-parton splitting function may have different strongly-ordered limits, which can be distinguished by the momentum labels in the kernel, once the ordering of the limits is fixed by the momentum mapping, ktr → kt + r → (k + t) + r in our convention. We always choose z = z t,k and z = z r , kt as independent variables. For quark splitting we have while for gluon splitting we find where the constants b  27) We see that the second term is present only if the three-parton splitting involves a twoparton sub-splitting with parent gluon.

D.3 Collinear-triple collinear counterterm
The collinear-triple collinear counterterm involves two successive collinear mappings, which leads to exact phase space factorisation in the iterated form dφ m+2 ({p}; Q) = dφ m ({ p } ( kt r ,kt) ; Q)[dp ( kt r ) 1;m ( p r , p kt r ; Q)][dp The one-parton factorized phase spaces [dp (kt) 1;m+1 (p k , p kt ; Q)] and [dp ( kt r ) 1;m ( p r , p kt r ; Q)] are given explicitly by eq. (C.2) after appropriate changes in labelling, including the replacement m → m − 1 in the second case. The Altarelli-Parisi splitting functions and the factor z(1 − z) in eq. (D.27) can be expressed as linear combinations of powers of momentum fractions. Consequently, the integral over the factorized phase space [dp (kt) 1,m+1 ] in the integrated collinear-triple collinear counterterm is written in terms of collinear functions I 1 (x kt , ǫ, α 0 , d 0 ; k) of eq. (C.1). In order to compute the subsequent integrals over [dp the variable x kt = y ktQ needs to be expressed in terms of p µ ktr instead of p µ kt (see eq. (C.8) with proper changes in labelling), x kt = y ktQ = α kt r + (1 − α kt r )x ktr v kt r . (D.30) Then using the abbreviations α = α kt r , v = v kt r (the integration variable corresponding to 1 − z r , kt ) and x ≡ x ktr , the integral representations (C.3) and (C.9), we find that the integrated collinear-triple collinear counterterm can be expressed using the following three types of integrals: × I 1 (α + (1 − α)xv, ǫ, α 0 , d 0 ; k) , k, l = −1, 0, 1, 2 , We can use the relations and to reduce the explicit computation of the I C integral to the case l = −1. In terms of the functions I

D.4 Collinear-double collinear counterterm
The collinear-double collinear counterterm also involves two successive collinear mappings, which leads to exact phase space factorisation in an iterated form similar to that in eq. (D.28). The integral over the factorized phase space measure [dp (kt) 1,m+1 ] leads to the same collinear integrals as in eq. (C.1). Then the necessary integrals over [dp (D.37) Again, x kt needs to be expressed in terms of p µ kt instead of p µ kt , x ≡ x kt , y = x ir , and using the v ↔ 1 − v symmetry of the integration measure, we find that the integrated collinear-double collinear counterterm can be expressed as a linear combination of the integrals (D. 39) In terms of the functions I C (x, y; ǫ, α 0 , d 0 ; k, l) we find the result given in eq. (4.7).

D.5 Collinear-soft-collinear-type counterterms
The collinear-soft-collinear-type counterterms in eqns. (4.8), (4.12) and (4.10), involve a collinear mapping followed by a soft mapping of the phase space, which leads to an exact factorisation of the original m + 2-particle phase space in the form where the factorized phase space measures [dp we have to express the invariants of the dependent momenta (with hat) with those of the independent ones (with tilde): with a similar expression for z kt, r /z r , kt .
To write explicit integral representations of I (i) C (i = 5, 6 and 7), we choose the specific Lorentz frame of eq. (C.11) with a different orientation for each function.
Integrated collinear-soft-collinear counterterm. Here and in the following, we will use the partial fraction identity below to disentangle the singularities associated with the factors of 1/s j r and 1/s l r , appearing in the eikonal factor (first term in the braces in eq. (D.41)): This is useful when computing the integral via iterated sector decomposition, while the original form is better suited to derive the Mellin-Barnes representation. Since 'undoing' the partial fractioning is trivial, below we will show the more elaborate form of the integrals, which are directly suited to treatment with sector decomposition.
The convenient frame for integrating the first term in eq. (D.45) is . . , sin χ l , cos χ l ) , p µ kt = E kt (1, . . . , sin φ kt sin χ kt , cos φ kt sin χ kt , cos χ kt ) , while for the second term we choose a frame where j and l are interchanged as compared to eq. (D.46). In terms of the scaled energy-like variable y r Q and the angular variables ϑ, ϕ and η, the two-particle phase space dφ 2 ( p r , K; Q) is given by eq. (C.12). The two-particle invariants y j r , y l r , y kt r have to be expressed in terms of the integration variables, i.e.
. we find the following explicit expression for I (5) C : The second term in the squared brackets in eq. (D.53) corresponds to the second term in the partial fractions of eq. (D.45). The integral representations of the two terms are formally identical, only the kinematic variables Y 2 and Y 3 are interchanged.
In terms of the functions I C we find the result given in eq. (4.9). In writing eq. (D.53), we have tacitly assumed that p µ j , p µ l and p µ kt are all distinct (massless) momenta, whose kinematics is furthermore unconstrained. For processes involving two or three hard final state partons, this is not the case, so the integral I (5) C with full kinematic dependence, as written above, first appears in computing NNLO corrections to processes with at least four hard final state partons.
When there are only three hard partons in the final state, the kinematics of the event is constrained because momentum conservation forces the final state momenta to be coplanar. Thus, in eq. (D.46) we have sin φ kt = 0, and the parametrisation of p j , p l and p kt simplifies accordingly: where we choose cos φ kt = −1, so that we may assume sin χ l and sin χ kt to be nonnegative. As a result of the constrained kinematics, we can first of all perform the cos η integration in eq. (C.12) using eq. (C.13). Second, out of the four kinematic invariants (x kt , Y j l ,Q , Y j kt ,Q , and Y l kt ,Q ) of the general case in eq. (D.53), only two are independent.
(Momentum conservation implies three constraints among the five variables E j , E l , E kt , cos χ l and cos χ kt .) However, it is convenient to leave the formal dependence on all four variables, with the constraints The physical region for Y j l ,Q and Y j kt ,Q is 0 ≤ Y j l ,Q , Y j kt ,Q ≤ 1 and Y j l ,Q +Y j kt ,Q ≥ 1.
The two-particle invariant y kt r becomes independent of φ kt , y kt r = 1 2 y kt Q y r Q (1 + sin χ kt sin ϑ cos ϕ − cos χ kt cos ϑ) , (D.57) therefore, the dependence on Y j l ,Q , Y j kt Q and Y l kt Q enters only through the angles cos χ l = cos χ(Y j l ,Q ) , cos χ kt = cos χ(Y j kt ,Q ) (D.58) for the first term in the partial fraction and the angles cos χ j = cos χ(Y j l ,Q ) , cos χ kt = cos χ(Y l kt ,Q ) (D.59) in the second one (with eq. (D.52) for cos χ(Y )). Integrating out cos η as in eq. (C.13), the integral (D.53) simplifies to in the second one (with eq. (D.52) for cos χ(Y )). Integrating out cos η as in eq. (C.13), the integral (D.53) reduces to (D.70) as in eq. (4.13).
Finally, the integral of the third term in the braces in eq. (D.41) is obtained by setting i = kt in the previous case, which implies Y i kt ,Q → Y i i ,Q = 0, so the integration over ϕ can be evaluated using eq. (C.14) and the integral in eq. (4.10) can be expressed as a linear combination of the integrals (D.71) (k = −1, 0, 1, 2) as in eq. (4.11).

D.6 Integrated collinear-double soft-type counterterms
The collinear-double soft counterterm is defined by an iterated application of a collinear and a soft momentum mapping, which results in an exact factorisation of the original (m + 2)-particle phase space very similar to that in eq. (D.40). The difference is that in the present case the soft measure involves the momentum p µ kt instead of p µ r . The integrand of the collinear integral is the spin-dependent splitting kernel of gluon splitting, with Lorentz structure and b (0) f k ft is defined in eq. (4.3), hence, the integral over the phase space measure [dp (kt) 1,m+1 ] involves both collinear functions I and I µν . Introducing the abbreviation f k ft (z k,t , z t,k , k ⊥,k,t ; ǫ)|ν , Then we can proceed as usual by writing the factorized phase space [dp For the eikonal factor (first term in the braces in eqs. (D.84)-(D.86)) we use the partial fraction identity (D.45) (with the substitution r → kt), and the j ↔ l symmetry to integrate only one term. Thus we are left with the following three integrals: (here Y corresponds to Y j l ,Q ), For the integrated collinear-triple collinear-double soft subtraction in eq. (4.16) the phase space factorisation is the same as for the integrated collinear-double soft subtraction, and we need to evaluate the integrals

E Integrating the soft-type terms
In this appendix, we discuss the integration of the soft-type counterterms of section 4.2.
In particular, we define and give an explicit integral representation of all I (i) S functions (i = 1, . . . , 12).

E.1 Soft-triple collinear-type counterterms
There are three integrated counterterms that involve a soft mapping of the momenta, followed by a collinear one, which leads to an exact factorisation of the original m + 2particle phase space in the form (E.1) The factorized phase space measures [dp Integrated soft-triple collinear counterterm. We begin by noting that the soft functions P (S) f i frft (originally defined in ref. [1]), which appear in eq. (4.18) can be written in the following unified form: ir This general form hides the fact that f t = g, which also implies C f it = C f i and C frt = C fr (used in eq. (4.39)). Furthermore, according to our definition of the momentum fractions, we have for any i, j and k. Therefore, performing the integration over the soft phase space [dp (t) 1,m+1 (p t ; Q)] first, we find that the three terms in eq. (E.2) all lead to known functions. The eikonal term gives (−1) times the soft function J 1 (Y, ǫ, y 0 , d ′ 0 ), while both soft collinear terms give 1 2 times the same soft-collinear function, K 1 (ǫ, y 0 , d ′ 0 ). For these functions we have the integral representations discussed in appendix C. Hence, we can express the integrated soft-triple collinear counterterms as linear combinations of two types of basic integrals, and The soft-collinear function K does not depend on any kinematic variable, therefore, it factorizes completely from the integral over [dp 1,m ]. The remaining integral can be expressed with the collinear function I(x; ǫ, α 0 , d 0 − 1 + ǫ, 0, k, 0, 1) (the parameter d 0 is shifted because this integral corresponds to m partons in the final state), so The soft integral depends on Y i r ,Q , that has to be expressed with the integration variables of [dp where we used the usual abbreviations α = α i r and v = v i r . Clearly, Y i r ,Q (x ir ; α, v), as well as the phase space measure in eq. (C.2) is symmetric in v ↔ 1 − v. Thus, it is sufficient to consider only the following coupled collinear and soft integral, In terms of the integrals I S and I S the integrated soft-triple collinear subtraction term can be written as in eq. (4.19).
Integrated soft-soft-collinear counterterm. To obtain the counterterm in eq. (4.20), we again perform the integration over the factorized soft phase space [dp (t) 1,m+1 (p t ; Q)] first. Now, we must distinguish two cases: (i) j, l and ir are three distinct labels, (ii) j or l coincide with ir. (We always have j = l.) In the first case, the integral over the factorized soft phase space simply leads to the J 1 function of eq. (C.10), and we find that the integrated soft-soft-collinear counterterm can be expressed as a linear combination of the integrals Using the collinear mapping formula p µ j,l = (1 − α i r ) p µ j,l , we find Consequently, the collinear and soft integrals decouple, and we obtain In case (ii) if e.g. j = (ir), the eikonal factor S (ir)l (t) evaluates as in eq. (2.22), and the integral over the soft phase space leads to the one-mass soft function J (1m) , defined in eq. (C.22) and computed in eq. (C.25). Hence we find that the integrated soft-soft-collinear counterterm can be expressed as a linear combination of the integrals .

(E.12)
This time the collinear integral does not decouple because the parameters Y ( i r ) l ,Q and β ( i r ) depend on the collinear integration variables. Using the definition of p µ ir and p µ l , we have p µ (ir) = (1 − α i r ) p µ ir + α i r Q µ and p µ l = (1 − α i r ) p µ l . Consequently, and (E.14) (k = −1, 0, 1, 2). In terms of the integrals I S and I S , the integrated soft-soft-collinear counterterm can be expressed as in eq. (4.21).
Integrated soft-triple collinear-soft-collinear counterterm. When computing the integrals in eq. (4.22), we first perform the integral over the soft phase space. Using the same frame as in eq. (C.23), we find that the integral leads to the one-mass soft-collinear function K (1m) defined in eq. (C.33) and computed in eq. (C.36). Thus the integrated softtriple collinear-soft-collinear counterterm can be expressed as linear combination of the integrals .
The soft-collinear integral K (1m) (β) does not decouple from the collinear one because β = β ( i r ) (y ir Q , α i r ) depends on the collinear integration variable α i r as in eq. (E.13). Thus, (E. 16) In terms of the integrals I S , we find the result as in eq. (4.23).

E.2 The soft-double soft-type counterterms
The remaining soft-type integrated counterterms involve two successive soft momentum mappings, which leads to an exact factorisation of the original m + 2-particle phase space in the form The factorized phase space measures [dp ( r ) 1,m ] and [dp (t) 1,m+1 ] are given in eq. (C.12), with appropriate changes in labelling (including m → m − 1 in the first case). We often need to express the two-particle scaled invariants of the momenta after the first mapping ('hatted momenta') with the final momenta obtained after the second mapping ('tilded momenta'). The relevant formulae, collected here for later reference, are: y i l = (1 − y r Q )y i l , y k r = y k r , y k Q = (1 − y r Q )y k Q + y k r , k = i, j, l .

(E.18)
Integrated soft-triple collinear-double soft counterterms. In computing the integrals in eq. (4.24), we first perform the integration over the soft phase space factor [dp (t) 1,m+1 ]. This integration leads to either the J 1 soft function, or the K 1 soft-collinear function of section C. Then the remaining integral over [dp .

(E.19)
When the result of the first integration is a K 1 function, it decouples from the second integral, which gives a soft-collinear function again. Thus we find that all terms in the integrand of the type 1 s jl integrate to the product of two soft-collinear integrals: On the other hand, when the result of the first integration is a soft function J 1 , the two integrals do not decouple. In order to compute the second integral over the phase space factor [dp ( r ) 1,m ], we choose the usual frame (C.11), with orientation specified by Using eq. (E.18), we compute Y i r ,Q = y i r y i Q y r Q = y i r [(1 − y r Q )y i Q + y i r ]y r Q = 1 − cos ϑ 2 − y r Q (1 + cos ϑ) . (E.23) Then we find that the term in the integrand of the type 24) leads to the integral To write eq. (E.25) in the above form, we made the usual substitution of cos ϑ → 1 − 2z.
Integrated soft-soft-collinear-double soft counterterm. In the case of the integral in eq. (4.26), we begin by integrating the eikonal function 1 2 S jl (t) over the soft phase space [dp (t) 1,m+1 ]. However, whenever j or l is equal to ir (recall that j = l), the eikonal factor evaluates as in eq. (2.22). Therefore, we have to distinguish two cases: (i) j, l and ir are three distinct labels, thus the integration over the first soft phase space leads to a J 1 soft function and we obtain the integral and (ii) j or l coincide with (ir), hence the integration over the first soft phase space leads to a J (1m) one-mass soft function and we find the integral (choosing j = ir for concreteness) (E.27) To proceed, we must express the parameters of the soft functions with independent momenta. In the first case, using eq. (E.18), we can express Y j l ,Q as Y j l ,Q = 4Y j l ,Q 1 − y r Q [2(1 − y r Q ) + 2y j r /y j Q ][2(1 − y r Q ) + 2y l r /y l Q ] .
(E. 28) In the second case the one-mass soft function also depends on the velocity of the momentum p µ i + p µ r , β ( i r ) = 1 − 4y i r y 2 ( i r )Q .
(E.29) Again using eq. (E.18), we find β ( i r ) = 2(1 − y r Q ) + 2y i r /y i Q + 2y r Q /y i Q 2 − 16y i r /y 2 i Q 2(1 − y r Q ) + 2y i r /y i Q + 2y r Q /y i Q , (E.30) Y ( i r ) l ,Q = 4(1 − y r Q )Y i l ,Q + 4y l r /(y i Q y l Q ) 2(1 − y r Q ) + 2y i r /y i Q + 2y r Q /y i Q 2(1 − y r Q ) + 2y l r /y l Q . (E.31) Turning to the integral over [dp ( r ) 1,m ], we use two different orientations of frames in the two cases. In the first case, we fix the orientation such that p µ j = E j (1, . . . , 1) , p µ l = E l (1, . . . , sin χ l , cos χ l ) , (E.32) that implies 2y j r y j Q = y r Q (1 − cos ϑ) , 2y l r y l Q = y r Q (1 − sin χ l sin ϑ cos ϕ − cos χ l cos ϑ) , (E. 33) where cos χ l = cos χ(Y j l ,Q ) (E. 34) with eq. (D.52) for cos χ(Y ). Thus we see that although the integral over [dp ( r ) 1,m ] is of soft-collinear-type, we cannot trivially integrate over ϕ because Y j l ,Q depends on it. Thus where β ( i r ) and Y ( i r ) l ,Q are given in eqs. (E.30) and (E.31), with the ratios of invariants in eq. (E.38) (x corresponds to y i Q , Y to Y i l ,Q and y to y r Q ). Our final result for the integrated counterterm defined in eq. (4.26) is presented in eq. (4.27).
Integrated soft-triple collinear-soft-collinear-double soft counterterm. The integral over the soft phase space factor [dp (t) 1,m+1 ] in eq. (4.28) is a one-mass soft-collinear integral K (1m) of eq. (C.33), which is computed in eq. (C.36), where the velocity of the momentum p µ i + p µ r is found in eq. (E.29). Then the integral over the measure [dp is again of soft-collinear-type, but the two integrals are coupled through β ( i r ) that depends on the integration variables of the second integral as in eq. (E.30). The ratio of the invariants 2y i r /y i Q is 2y i r y i Q = y r Q (1 − cos ϑ) , (E.42) in a frame whose orientation is fixed by setting with x corresponding to y i Q , y to y r Q and β(x, y, z) = (1 − y + yz + y/x) 2 − 4yz/x 1 − y + yz + y/x , (E.45) where, as usual, we set cos ϑ → 1 − 2z. Our final result for the integrated counterterm defined in eq. (4.28) is presented in eq. (4.29).
Integrated soft-double soft counterterms. There are two types of integrated softdouble soft counterterms: an 'abelian' one in eq. (4.30) and a 'non-abelian' one in eq. (4.31).
Let us first consider the 'abelian' case. Performing the integral over [dp (t) 1,m+1 ] first, we obtain a soft function −J 1 (Y j l ,Q , ǫ, y 0 , d ′ 0 ). In order to compute the integral over [dp S ;ik,jl = − 16π 2 S ǫ Q 2ǫ [dp we must first express Y j l ,Q in terms of p µ j and p µ l as in eq. (E.28) (j, l = r), which is seen to depend on the integration variables in [dp ( r ) 1,m ] through the appearance of y j r and y l r in the denominator. Thus, the integration over [dp ( r ) 1,m ], which is of soft-type, is nontrivial. Also, although i = k and j = l, there is no restriction on whether or not i, k is equal to j, l. Thus we must consider the following three cases: (i) all of i, k, j and l are distinct, (ii) only three of the four indices are distinct and, e.g. l = k, and (iii) only two indices are distinct and e.g. j = i and l = k.
Case (i) requires at least four hard partons in the final state. Hence the corresponding integrated counterterm, [S t S (0) rt ] (i,k)(j,l) with all labels distinct, does not enter a computation of two-or three-jet quantities, and we will not consider it in this paper.
In case (ii), we have Y j l ,Q → Y j k ,Q , and this is expressed with the independent momenta as in eq. (E.28), after a l → k replacement. To evaluate the integral over [dp ( r ) 1,m ] in eq. (E.46), we use the partial fraction identity (D.45) for the eikonal factor S i k ( r ) (with the substitutions j → i and l → k). Further, we restrict our attention to the case when there are precisely three hard partons in the final state. As discussed around eq. (D.54), this leads to a constrained kinematics for the three momenta p µ i , p µ k and p µ j , and we take this into account below. It is convenient to introduce two different orientations of the frame eq. (C.11) and integrate the two terms of the partial fraction in these different frames. In the first, we set with cos χ(Y ) given by eq. (D.52). In the second frame we exchange i and k, whose only effect on the integrand is to interchange χ ij and χ jk . Thus we find I (11) S ;ik,jk (Y 1 , Y 2 , Y 3 ;ǫ, y 0 , d ′ 0 ) = −4Y 1 Γ 2 (1 − ǫ) 2πΓ(1 − 2ǫ) Above, Y 1 ≡ Y i k ,Q , Y 2 ≡ Y i j ,Q and Y 3 ≡ Y j k ,Q . The two terms in the square bracket in eq. (E.52) correspond to the two terms of the partial fraction in eq. (D.45). Their integral representations are formally identical, only the kinematic invariants Y 2 and Y 3 are interchanged.
In terms of the integral I S ik,jk the integrated subtraction term defined in eq. (4.30) can be expressed as in eq. (4.32).

(E.62)
In terms of the integrals I S ik,ik and I (12) S ;ik , the integrated subtraction term defined in eq. (4.31) can be expressed as in eq. (4.33).

F Integrating the soft-collinear-type terms
In this appendix, we discuss the integration of the soft-collinear-type counterterms of section 4.3. In particular, we define and give an explicit integral representation of all I (i) CS functions (i = 1, . . . , 3).

F.1 Soft-collinear-triple collinear-type counterterms
The soft-collinear-triple collinear-type counterterms involve a soft momentum mapping followed by a collinear one, which leads to an exact factorisation of the original m+2-particle phase space in the form of eq. (E.1). To evaluate the integrals of eqs. (4.34) and (4.35) over the factorized one-particle phase space measures in eq. (E.1), we first observe that the integral over the soft measure is a K 1 soft-collinear function, and the remaining the integral over the collinear measure can be expressed as linear combination of the integrals [dp

F.2 Soft-collinear-double soft-type counterterms
The soft-collinear-double soft-type counterterms involve two successive soft momentum mappings, which leads to an exact factorisation of the original m + 2-particle phase space in the form of eq. (E.17). To evaluate the integrals of eqs. (4.36), (4.37) and (4.38) over the factorized one-particle phase space measures in eq. (E.17), we first observe that the integral over the soft measure [dp (t) 1,m+1 ] is again a K 1 soft-collinear function, and the remaining integral over the second soft measure contains either an eikonal factor or its collinear limit, The second integral decouples from the first one also for I CS because K 1 is independent of the kinematics. Then, the final integral over [dp CS (Y ; ǫ, y 0 , d ′ 0 ) = J 1 (Y ; ǫ, y 0 , d ′ 0 − 1 + ǫ)K 1 (ǫ, y 0 , d ′ 0 ) , (F.5) where Y corresponds to Y j l ,Q .
In terms of the functions I (i) CS , (i = 1, 2, 3) we find the results given in eq. (4.39).