{'fraction_non_alphanumeric': 0.054241707598083805, 'fraction_numerical': 0.014371475350191013, 'mean_word_length': 4.087613759062163, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Integration of collinear-type doubly unresolved counterterms in NNLO jet cross sections

In the context of a subtraction method for jet cross sections at NNLO accuracy in the strong coupling, we perform the integration over the two-particle factorised phase space of the collinear-type contributions to the doubly unresolved counterterms. We present the final result as a convolution in colour space of the Born cross section and of an insertion operator, which is written in terms of master integrals that we expand in the dimensional regularisation parameter.


Introduction
The high precision of the experimental measurements at the LHC calls for an evaluation of the jet, heavy-quark, vector-boson and Higgs-boson production rates at hadron colliders which is at least as precise. In recent years, we have witnessed fast progress in the evaluation of the production rates mentioned above, in association with many jets, at next-to-leading order (NLO) accuracy in the strong coupling constant α s . That rested upon efficient methods to compute one-loop amplitudes with many legs, and upon subtraction algorithms to evaluate QCD cross sections at NLO at the partonic level [1][2][3][4][5], or at the hadronic level, through interfaces, like MC@NLO [6] or POWHEG [7][8][9] with parton-shower event generators. Such algorithms are based on the universality -i.e., on the independence from a specific scattering process -of the infrared emissions.
No subtraction algorithm to compute cross sections at hadron colliders at NNLO accuracy, though, has been devised yet. In order to do that, one must define subtraction terms that properly regularise the real-emission phase-space integrals and then one must combine the integrated form of those counterterms with the virtual contributions, so as to cancel the infrared divergences of the loop amplitudes, in such a way that the cancellation of both the kinematic singularities in the real-emission pieces and the explicit ǫ-poles in the virtual pieces be local. This implies that the subtraction terms and the real-emission contributions must tend to the same value in all kinematic limits where the latter diverge, and that the cancellation of explicit ǫ-poles between the integrated subtraction terms and the virtual contributions must take place point-wise in phase space. In particular, that implies that it is possible to write the integrated counterterms in such a way that they can be explicitly combined with virtual contributions, before phase-space integration. Practically, the locality of the subtraction scheme is also important to ensure good numerical efficiency of the algorithm. In broad outline, we remind how this occurs at NLO. After fixing the leading-order m-jet cross section as the integral of the fully differential Born cross section dσ B m of m final-state patrons over the available m-parton phase space defined by the jet function J m , can be rewritten as a sum of finite integrals, where each of the two terms on the right-hand side, is integrable in four dimensions by construction [1][2][3][4][5]. Likewise, we can write the NNLO correction to the cross section as a sum of three contributions, the doubly unresolved, the one-loop singly unresolved, and the two-loop doubly virtual terms, and rearrange it as follows, where each of the three terms on the right-hand side, is integrable in four dimensions by construction [20,22,23]. The counterterms which contribute to dσ NNLO m+2 and to dσ NNLO m+1 were introduced in Refs. [22] and [23], respectively. The integral of the real-virtual counterterms (the last two terms of Eq. (1.10)) was performed in Refs. [25][26][27]. The integral of the iterated singly unresolved counterterm (the third term of Eq. (1.10)) was performed in Ref. [31]. In this paper, we compute the integral of the collinear-type contributions to the doubly unresolved counterterm (the second term of Eq. (1.10)). The soft-type contributions to the same integral are presented in a companion paper [51].

Notation
In this paper, we use the notation introduced in Ref. [31], which we recall in this Section.

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 from the middle of the alphabet, i, j, k, l, . . . , while letters r, s denote unresolved final-state partons.
We adopt the colour-and spin-state notation of Ref. [2]. In that notation the amplitude |M for a scattering process 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 (2.1) In this paper we only need |M in the Born approximation, denoted by |M (0) . Using the colour-state notation, we define the two-parton colour-correlated squared tree amplitudes as and similarly the four-parton colour-correlated squared tree amplitudes, where we also introduced the ⊗ product notation to indicate the insertion of colour charge operators between M (0) | and |M (0) . 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 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.5)

Cross sections
In this paper we shall need to use only tree-level n-parton production cross sections, with n = m, the Born cross section, and n = m + 2, the so-called doubly real correction. We have 6) 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, The symbol {n} denotes summation over 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.9)

Momentum mappings and phase-space factorisation
The subtraction terms are written in terms of the momenta obtained by the mappings of Ref. [22], where X RR −→ may label a triple collinear, double collinear or soft collinear mapping 2 . In particular, the soft collinear mapping is obtained by the iterated application of a basic collinear mapping, followed by a basic soft mapping, (2.11) while the former two cases do not have such an iterated form. 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 is inherited by the label, and we write for instance i , ir and irs , where the last 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.5)), we shall omit the tilde from flavour and colour indices. 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 . In a computation at NLO accuracy, the two definitions lead to the same expressions. At NNLO they lead to slightly different bookkeeping of the IR and UV poles at intermediate steps of the computation, our definition leading to somewhat simpler expressions. Of course the physical cross section does not depend on these details. 2 Integration of the subtraction terms obtained with the double soft mapping are presented in the companion paper [51].
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 . The momentum mappings are chosen as to lead to an exact factorisation of the phase space. For the triple collinear phase-space mapping, we have, (2.14)

Kinematic variables
The following types of kinematic variables are used to write the doubly unresolved subtraction terms.
• Two-particle invariants, such as (2.15) Two-particle invariants scaled with Q 2 are denoted by y ij = s ij /Q 2 .
• Momentum fractions z i,r for the splittingsp ir → p i + p r and z i,rs for the splittings p irs → p i + p r + p s , z i,r = y iQ y iQ + y rQ and z i,rs = y iQ y iQ + y rQ + y sQ , (2.16) with z i,r + z r,i = 1 and z i,rs + z r,si + z s,ir = 1 (z r,si and z s,ir are computed by obvious permutations).
• We also use the eikonal factors, S jk (s) = 2s jk s js s ks . (2.17) As mentioned above, the sum of two momenta is often abbreviated with the two indices in parentheses e.g., p µ i + p µ r = p µ (ir) , which is also used in other occurrences, such as Finally, we express the integrated subtraction terms as functions of the following (combinations of) invariants: In the rest frame of Q µ , y i Q = 2Ẽ i /Q is simply the scaled energy of partonĩ, while Y i j ,Q is related to χ ij , the angle between momentap µ i andp µ j , by cos 3 Integrating the doubly unresolved approximate cross section

The integrated approximate cross section and insertion operator
The doubly unresolved approximate cross section times the jet function is defined as, where the ⊙ notation reminds us that the set of momenta entering J m is different from term to term, All terms on the right-hand side of Eq. (3.2) are defined in Ref. [22].
Using that all momentum mappings lead to the exact factorisation of phase space, dφ m+2 ({p}; Q) = dφ m ({p}; Q)[dp 2;m ], and that the jet function J m does not depend on the variables of the factorised two-parton measure, [dp 2;m ], we can compute the integral of Eq. (3.2) over the two-parton factorised phase space independently of J m , In order to make the integrated subtraction terms dimensionless in color space, we have factored out quadratic Casimir operators in Eq. (3.3).

Eq. (3.3)
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 counting is performed as in the Appendix of Ref. [31]. As a result of the summation over the unresolved flavours, we obtain functions -the flavour-summed integrated counterterms -denoted by X (0) (j,l)...   After summing over unobserved flavours, the integrated doubly unresolved approximate cross section can be written as, where the insertion operator (in colour space) has three contributions according to the possible colour structures, . In terms of flavour-summed integrated counterterms discussed above, the kinematic functions of the insertion operator can be written as, .

(3.7)
On the right-hand side of Eq. (3.7), the flavour-summed functions depend on the kinematics through variables of the type x i and Y ij,Q . The latter dependence stems from integrating an eikonal factor which is always multiplied by 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
On the right-hand side of Eq. (3.7), the flavour-summed functions are the following combinations of the integrated subtraction terms, 1. Triple collinear: (3.8) The factor (n f − 1), which appears after summing over the unresolved flavours, can be traded for n f by introducing the integrated subtraction term [C This definition matches the decomposition of the same-flavour triple splitting function computed in Ref. [52]. Thus, we obtain 2. Triple collinear -soft collinear: i.e., it is independent of the flavour f .

Double collinear -double soft:
i.e., it is independent of the flavours f 1 and f 2 .
In this paper, we compute the collinear-type counterterms, which do not involve S (0) rs , i.e., those in Eqs. (3.10)-(3.14). The integration of the subtraction terms which are obtained by the double soft mapping is presented in a companion paper [51].

Integrated counterterms
In this section we define the integrated counterterms and compute them in terms of master integrals.

Integrated triple collinear counterterm
The integral of the triple collinear counterterm is Here P (0) f i frfs (z i,rs , z r,is , z s,ir , s ir , s is , s rs ; ǫ) are the complete 3 spin-averaged three-parton splitting functions, recalled in Appendix B, where we explain why some of these are different from those presented in Ref. [52]. For consistency of the complete subtraction scheme one must use the forms given in Appendix B. The function f (α 0 , α irs , d(m, ǫ)) is defined and its role is explained in Appendix A. The other subtraction terms discussed in this paper will also contain such harmless modifications as compared to the original ones in Ref. [22].
Following the decomposition of the triple collinear splitting functions into abelian and non-abelian pieces [52], we also decompose the integrated functions [C Table 1. Non-zero coefficients of the I They are computed in Appendix C. The result can be written as where the constants a f i frfs denote colour-factor ratios as follows: The non-zero coefficients are listed in Tables 1-5, while the integrals I (j,k,l,m) 2C ,n (x, ǫ; α 0 , d 0 ) (n = 1, . . . , 5) are defined and computed in Appendix C. Table 2. Non-zero coefficients of the I

Integrated double collinear counterterm
The integral of the double collinear counterterm is The integrated counterterm is computed in Appendix D. The result is where the constants a f i fr denote the colour factor ratios and the coefficients c f i fr are those already defined in Ref. [27] stripped of their colour factors, as presented in Table 6. Table 3. Non-zero coefficients of the I
The integrated counterterm is computed in Appendix E. The result is if e.g., j = (ir), where a f i fr and the coefficients are given in Eq. (4.7) and Table 6 respectively.

Results
After evaluating the integrals as explained in the Appendices, we obtain the kinematic functions defined in Eq. (3.7) as Laurent expansions in ǫ. We compute those expansions analytically up to O(ǫ −2 ), while the remaining coefficients up to O(ǫ 0 ) are calculated numerically via sector decomposition, see [53] and references therein. In obtaining our results, , as well as α 0 and y 0 are left symbolic throughout the analytic computation and sector decomposition.
To present analytic results for the non-soft counterterms up to O(ǫ −2 ), we use the following abbreviations, where γ q (n f ) is actually independent of the number of light flavours, but introducing the flavour dependence formally makes possible a flavour-independent notation. Furthermore, the kinematic functions are dimensionless in colour space, hence these n f -dependent γ f (n f ) constants are related to the constans γ f often used in the literature [54] by γ f → C f γ f (n f ).
The remaining coefficients in the Laurent expansion are computed numerically. By way of illustration, we present results for the flavour-summed counterterms in Tables 7-19 for two kinematic points: one relevant for 2-jet production and the other corresponding to the fully symmetric configuration of final state momenta in 3-jet production. In terms of the invariants introduced in Eq. (2.19), these two phase-space points correspond to 2-jet : 3-jet symmetric : In the numerical computations, we chose the following values for the phase-space cut parameters: Tables 7-19, a displayed error estimate of ±0.000 implies that the numerical uncertainty for that particular value is smaller than 10 −3 .

Conclusions
We have computed the integrals over the two-particle factorised phase space of the collineartype contributions to the doubly unresolved counterterm of the NNLO subtraction formalism, defined in Refs. [22,23]. We presented those integrals in terms of parametric rep-colour resentations that are suitable for evaluation with sector decomposition. After evaluating them, we checked the numerical results with the publicly available code SecDec [55], always finding agreement within the numerical uncertainty of the integrations. The soft-type contributions are presented in a companion paper [51]. By these two papers, we complete the integration of the subtraction terms over the unresolved phase spaces, and the computation of the finite cross section in Eq. (1.9) becomes feasible for electron-positron annihilation into two and three jets. Although the formalism is complete, for a higher number of jets some more work is required, because some integrals were evaluated specifically for three-jet kinematics.

A Modified doubly real subtraction terms
We outline a simple modification to the NNLO subtraction scheme presented in Refs. [22,23]. Parts of these modifications were introduced previously: those relevant to the singly unresolved approximate cross section dσ RR,A 1 m+2 in Eq. (1.8), and to the approximate cross sections in Eq. (1.9), were introduced in Ref. [25], while those relevant to the iterated doubly unresolved approximate cross section dσ RR,A 12 m+2 appearing in Eq. (1.8) were presented in Ref. [31].
Recall that the doubly unresolved approximate cross section can be written symbolically as in Eq. (3.1) dσ where the doubly unresolved approximation is a sum of terms (see Eq. (3.2)). The precise definition of these terms involves the momentum mappings discussed in Section 2.3. All such mappings lead to an exact factorisation of the m+2-particle phase space, symbolically written as The only feature of the factorised phase spaces [dp 2;m ] that is relevant presently is that they carry a dependence on the number of partons, m, of the form [dp and finally [dp The subtraction terms, as originally defined in Ref. [22] do not depend on the number of hard partons, thus the m-dependence of the factorised phase space measures is carried over to the integrated counterterms, where this dependence enters in a rather cumbersome way (see e.g., Eqs. (A.9) and (A.10) of Ref. [4]). Thus, as in Refs. [25,31], we reshuffle the m-dependence of the integrated counterterms into the subtraction terms themselves, where it appears in a very straightforward and harmless way 4 , through factors of (1 − α) and/or (1 − y) raised to m-dependent powers. We gather the results in Table 20, where together with the subtraction terms, we give the momentum mappings used to define the term(s) and the expression which multiplies the original counterterm to produce the modified one. The function f in Table 20 is defined as It is important to note that the functions d(m, ǫ), d ′ (m, ǫ) and constants α 0 , y 0 in Table 20 Doubly unresolved counterterms  Table 20. The modified doubly 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).
are the same as those in all other modified subtraction terms, discussed in Refs. [25,31].
The form of the exponents d(m, ǫ) and d ′ (m, ǫ) is actually fixed by the prescription adopted in Ref. [25] (see in particular Eqs. (3.2), (3.12) and (3.13) therein) 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. [25] i.e., where D 0 , D ′ 0 ≥ 2 are integers, while d 1 , d ′ 1 are real.

B Spin-averaged splitting kernels
In this Appendix, we list the spin-averaged splitting kernels. Although some of these already appeared elsewhere, for the sake of completeness, we give all the splitting functions of our computations.

B.1 Two-parton kernels
The azimuthally averaged two-parton splitting kernels are well known, In our convention the ordering of the labels on the splitting kernels is usually meaningless, but in Eq. (B.3) z means the momentum fraction of the second label, i.e., P gq (z; ǫ) = P (0) qg (1 − z; ǫ). The other two cases are symmetric with respect to z ↔ 1 − z.

B.2 Three-parton kernels
The spin average of the splitting kernels was computed in Ref. [52], however the forms presented there for gluon splittings are not suitable for us, as it was explained in Ref. [22]: in the gluon splitting kernelsP g i qrqs andP g i grgs , the terms that depend on the transverse momenta must always be written in the form k µ ⊥,j k ν ⊥,k /k ⊥,j ·k ⊥,k (k can be equal to j). Otherwise the collinear behaviour of the counterterm cannot be matched with that of the singly collinear counterterm in the singly unresolved region of phase space. The correct azimuth dependence can be achieved by substitutions according to the following replacements, where {k, l} = {i, r, s} \ {j} and j can be i, r or s. With these forms azimuthal averaging amounts to the simple substitutions, where . . . denote spin averaging, and j = k is also allowed. For quark splitting into unequal and equal quark flavours we have, and respectively (hence Eq. (3.9)), where For splitting into a quark and a gluon pair we find, We call attention to the normalisation (in colour space) of the abelian and non-abelian parts of the splitting function in Eq. (B.10), which is the same as in Ref. [20]. Notice that a factor of C F , respectively C A , is made explicit in the definition of the abelian, respectively non-abelian, part as compared to the complete splitting function. However, we prefer to define [C (nab) f i frfs . In particular, we must remember to include the factors of C F and C A explicitly with P (ab) f i frfs and P (nab) f i frfs , i.e., we must set P (0) f i frfs to obtain the correct normalisation.
For gluon splitting into a gluon and a quark pair we have (B.14) and 1 s 2 Note again that the normalisation (in colour space) of the abelian and non-abelian parts of the splitting function in Eq. (B.13) matches that of Ref. [20], i.e., a factor of C F and C A is made explicit in the definition of the abelian and non-abelian piece respectively. Hence, the comments below Eq. (B.12) apply in this case as well.
Finally for splitting into three gluons we have 2s irs z i,rs z s,ir − s ir z s,ir − s rs z i,rs + s is z r,is 2 + (5 permutations) .

C.1 Master integrals
The triple collinear momentum mapping leads to the exact factorisation of phase space as given by Eq. (2.12), where the two-particle factorised phase space can be written as [dp with p µ (irs) = (1 − α)p µ irs + αQ µ . When writing Eq. (4.1), we have used the fact that the spin correlations present at the level of the factorisation formula vanish upon azimuthal integration via the usual arguments, since k µ ⊥,j,k (j, k = i, r and s) as defined in Ref. [22] is orthogonal top µ irs . Therefore, the integrals of the spin-dependent and spin-averaged splitting functions are equal.
The spin-averaged triple collinear functions P (0) f i frfs depend on six variables: three twoparticle invariants: s ir , s is , s rs and three momentum fractions: z i,rs , z r,is and z s,ir , and are given in Appendix B. When counting the number of independent kinematical structures in Eqs. (B.8)-(B.16), we make two observations. Firstly, of the six variables, only four are independent, since s ir + s is + s rs = s irs and z i,rs + z r,is + z s,ir = 1. Secondly, from Eq. (C.1) it is clear that the factorised phase-space measure [dp [dp (irs) 2;m (p r , p s ,p irs ; Q)]P (p σ(i) , p σ(r) , p σ(s) ) , (C.2) where σ ∈ S 3 , with P an arbitrary function of momenta (a term in the splitting functions). This further reduces the number of independent structures.
Making use of the constraints among variables and the S 3 permutation symmetry, we find 46 structures to integrate 5 , that can be grouped into five classes, where we introduced the scaled two-particle invariants t kl = s kl /s irs (k, l = i, r, or s). The exponents take values as given in Tables 1-5. With these five classes, we can give new forms of the spin-averaged splitting functions that lead to the same integrated triple collinear subtraction terms, we map the region of integration onto the unit hypercube such that physical singularities are on the borders. Solving the constraints in a particular way, we find a parametrisation over the unit hypercube where all physical singularities are on the border except z s = 0, which introduces a line singularity (see below). We find [dp In terms of the variables t is , τ rs , v r and w s we have and (C. 13) We see that z s = 0 corresponds to α = 0, w s = 0, 1 and τ rs = v r , hence the line singularity. Since the only integral involving 1/z s is I 2C ,4 , this is the only place where the line singularity has to be resolved.
Using the parametrisation of Eq. (C.11) and Eqs. (C.12) and (C.13), we find the following explicit parametric integral representations of the basic integrals. (C.14) The expressions for I 2C ,n , n = 1, 2 and 3 in Eqs. (C.14)-(C.16) are suitable for evaluation with general purpose sector decomposition codes as they stand. In I 2C ,5 , Eq. (C.18), the terms in the braces require some care. Although this factor is just 1/(1 + z s ), which is finite and hence in principle it can be simply carried through sector decomposition as it is, nevertheless depending on the precise internal implementation, it may lead to undefined arithmetic expressions in the (symbolic) computation. We find that no such trouble arises either with our private implementation or with SecDec, if we gather the two terms in the braces over a common denominator.
I 2C ,4 requires special attention, as could be anticipated by the presence of the factor 1/z s . There are two problems: firstly, the presence of the square roots in the vanishing denominator in the braces on the last line of Eq. (C.17) interferes with the treatment of overlapping singularities. (This could be solved by simply making the expression squarefree.) Secondly, as we have indicated, there is a line singularity inside the integration region at τ rs = v r . We address both issues by deriving an alternative representation for I 2C ,4 , which is free of both square roots and line singularities. The price to pay is that the new representation is quite cumbersome, The derivation of this alternate form relies on rewriting the integration over v r and w s , as an angular integral, , (C. 21) where in a suitable frame, (the . . . denotes vanishing components) and q µ = (1, ..'angles'.., sin ϑ sin ϕ, sin ϑ cos ϕ, cos ϑ) .
(C.24) (Here ..'angles'.. represents those angular variables that are trivial to integrate, since the integrand does not depend on them.) Clearly both p 1 and p 2 are massive and time-like, p 2 1 , p 2 2 > 0, and p 1 · p 2 > 0. The angular integral can be written in terms of a Mellin-Barnes representation according to Eq. (60) of Ref. [56]. Finally, we can perform the Mellin-Barnes integrations at the expense of reintroducing two real integrations over the unit interval. We obtain a new two-dimensional real integral representation, free of square roots and singularities inside the integration region.

D.1 Master integrals
The double collinear momentum mapping leads to the exact factorisation of phase space as given by Eq. (2.13), where the two-particle factorised phase space can be written as, [dp (ir;js) 2;m (p r , p s ,p ir ,p js ; where p µ (ir) = (1 − α − β)p µ ir + αQ µ and p µ (js) = (1 − α − β)p µ js + βQ µ . When writing Eq. (4.5), we have used that spin correlations generally present at the level of factorisation formulae cancel after azimuthal integration, since k µ ⊥,i,r and k µ ⊥,j,s as defined in Ref. [22] are orthogonal top µ ir andp µ js respecitvely. Examining the actual form of the spin-averaged Altarelli-Parisi functions, and using the symmetry of the phase-space measure under the momentum exchanges p i ↔ p r , or p j ↔ p s , and under exchange of the pairs (p i , p r ) ↔ (p j , p s ), we see that we must compute the following integrals, Using the parametrisation of Eq. (D.7) and the expressions in Eqs. (D.5) and (D.6), we find the following parametric integral representation for the master integral, Eq. (D.8) is directly suitable for treatment with general purpose sector decomposition codes.

E.1 Master integrals
The consecutive collinear and soft mappings lead to the exact factorisation of phase space as given by Eq. (2.14), where the one-particle factorised phase spaces can be written in the following form. For the collinear mapping we have [dp For the soft mapping we find [dp where y ≡ yŝ Q and the momentum K is massive with K 2 = (1 − y)Q 2 . As the notation above indicates, α and y are integration variables. To integrate Eqs. (4.8)-(4.13) over the factorised phase spaces of Eqs. (E.1) and (E.2), we can pass to the azimuthally averaged Altarelli-Parisi splitting functions, since the azimuthal correlations generally present in Eqs. (4.8)-(4.13) vanish after integration by the usual argument (the transverse momentumk µ ⊥,i,r in the splitting kernels is defined to be orthogonal to the parent momentump µ ir ). The Altarelli-Parisi functions can be expressed as linear combinations of powers of momentum fractions, so we need to compute integrals of the form, We perform the integration over [dp The successive singly soft mapping of the momenta rescales those two-particle invariants that do not involve the soft parton, while those involving the soft parton s become yĵŝ = yjŝ. Finally, we have for j = s, while yŝ Q = y is simply an integration variable. Therefore, the integrands are expressed as follows. For the eikonal factor, we find if j and k are distinct from (ir), while if e.g., j coincides with (ir). Finally, we have Hence, we define five soft collinear master integrals, [dp (ir) 1;m+1 (p r ,p ir ; Q)] [dp [dp (ir) 1;m+1 (p r ,p ir ; Q)] [dp Turning to [dp (ŝ) 1;m (p s , K; Q)] in Eq. (E.2), we use the (scaled) energy and the angular variables ofp µ s in some specific Lorentz frame to write the two-particle phase space dφ 2 (p s , K; Q) explicitly. In particular, we work in the rest frame of Q µ , where Q µ = √ s(1, . . .) andp µ s is parametrised as, p µ s =Ê s (1, ..'angles'.., sin ϑ sin ϕ sin η, sin ϑ sin ϕ cos η, sin ϑ cos ϕ, cos ϑ) . (E.21) The specific orientation of this frame is chosen below according to what is most convenient for computing each integral. However, independent of orientation, in terms of the scaled energy-like variable, and the angular variables ϑ, ϕ and η, the two-particle phase space dφ 2 (p s , K; Q) reads We are now ready to display the master integrals.

(E.24)
Then, expressing the various invariants in terms of the integration variables in the chosen frame of Eq. (E.24), we find x ir y(1 − sin φ ir sin χ ir sin ϑ sin ϕ cos η − cos φ ir sin χ ir sin ϑ cos ϕ − cos χ ir cos ϑ) , where cos χ k = cos χ(Y j k ,Q ) , cos χ ir = cos χ(Y j ir ,Q ) The factor ykŝ in the denominator of the eikonal factor vanishes at cos ϕ = 1 and cos ϑ = cos χ k . Hence, the integrand has a line singularity that we remove by performing partial fractioning as in Refs. [31,57], Then we find where yjŝ, ykŝ and y irŝ are understood to be functions of the integration variables as given in Eqs. (E.25)-(E.27). We also used y ir Q = x ir , and performed the integration over εŝ with the δ-function in Eq. (E.23). The choice of frame in Eq. (E.24) is convenient for integrating the first term in the partial fraction, while the second term is more straightforward to integrate in a frame where k and j are exchanged. Upon performing this exchange, we find that the functional form of the master integral is unchanged, and we must simply interchange Y j ir ,Q with Y k ir ,Q to obtain one term from the other.

(E.39)
Note that Y k ir ,Q is also easily expressed with Y j k ,Q and Y j ir ,Q , Y k ir ,Q = Y j k ,Q + Y j ir ,Q − 2Y j k ,Q Y j ir ,Q + 2 Y j k ,Q (1 − Y j k ,Q )Y j ir ,Q (1 − Y j ir ,Q ) .
(E.40) The physical region in Y j k ,Q and Y j ir ,Q is simply given by 0 < Y j k ,Q , Y j ir ,Q < 1 and Y j k ,Q + Y j ir ,Q > 1.
In the case of two-jet production, this integral does not appear at all because we do not have three independent hard final state momenta.
As for I (E.57) The parametric forms for the integrals I