Analytic integration of soft and collinear radiation in factorised QCD cross sections at NNLO

Within the framework of local analytic sector subtraction, we present the full analytic integration of double-real and real-virtual local infrared counterterms that enter NNLO QCD computations with any number of massless final-state partons. We show that a careful choice of phase-space mappings leads to simple analytic results, including non-singular terms, that can be obtained with conventional integration techniques.


Introduction
Computing QCD cross-sections at next-to-next-to-leading order (NNLO) in the strong coupling is becoming mandatory to provide sufficiently precise fixed-order predictions for many processes of interest at high-energy colliders.This precision goal has led to the development of a host of new techniques in perturbative quantum field theory, ranging from the determination of parton distributions, to jet algorithms and of course to the calculation of high-order scattering amplitudes (for a recent review, see [1]).
One of the problems that need to be efficiently tackled in order to perform multi-parton NNLO QCD calculations is the cancellation of infrared singularities.Indeed, it is well known that, beyond leading order (LO) in QCD, both virtual corrections and real-radiation corrections contribute to any infrared-safe cross section: while these contributions are separately infrared (IR) singular, their sum (after UV renormalisation of virtual corrections) gives finite predictions for physical observables [2,3].This cancellation is well-understood in principle, but the increasing complexity of scattering amplitudes at high orders, and the intricate dependence of many collider observables on experimental cuts and jet algorithms, lead to significant difficulties in the practical implementation of the cancellation.
Subtraction algorithms form a class of proposed solutions to this problem.The basic ingredient of subtraction is the construction of universal infrared counterterms, defined locally in the radiative phase spaces.Such counterterms are required to mimic the behaviour of the radiative squared matrix element in all singular phase-space regions; on the other hand, they must be simple enough to be integrated over unresolved degrees of freedom in d = 4−2 dimensions, in order to analytically cancel the poles in arising from virtual corrections.Given such a set of counterterms, one proceeds by subtracting the counterterms from the radiative squared matrix element, so that the resulting expression can be numerically integrated without encountering singular contributions.One then adds to virtual corrections the integral of the counterterm over the radiative degrees of freedom, thus cancelling all infrared poles, and without having introduced any approximations in the distribution of the chosen infrared-safe observable.
At next-to-leading order (NLO), subtraction is well understood and successfully applied to a vast ensemble of observable multi-parton distributions.The most used subtraction methods at NLO are the Frixione-Kunszt-Signer (FKS) [4] scheme and the Catani-Seymour (CS) [5] algorithm.One order higher in the perturbative expansion (at NNLO), the development of a fully general and efficient subtraction method has been the subject of active research by many groups for several years.The literature is too vast to be comprehensively cited, but the main characteristics and some important applications of the most developed methods can be found in Refs.[6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24].It must also be mentioned that subtraction is not the only possible approach to the problem: an alternative viewpoint is provided by slicing methods, where an infrared cutoff is introduced to isolate the singular regions of the radiative phase space, and approximate expressions for the matrix elements are employed below the cutoff scale.Such methods were successfully used already at NLO [25,26], and have been applied at NNLO to a number of important processes [27][28][29][30][31][32][33][34][35][36].Furthermore, new ideas have been recently proposed [37][38][39], including theoretical developments concerning infrared factorisation [40], and the analysis of the infrared structure of Feynman diagrams [41][42][43], as well as purely numerical methods based on the cancellation of singular contributions at the integrand level, before loop and phase-space integrals are performed [44][45][46][47][48][49].Finally, the first developments for the extension of some of these tools to N 3 LO have been presented [50][51][52][53].This vast activity bears witness to the fact that the problem of subtraction (or more generally the problem of local cancellation of IR divergences) at NNLO is very intricate: available NNLO schemes are often characterised by a remarkable degree of complexity if compared with the NLO ones, and do not always feature desirable aspects such as universality, analytic control, and full locality in phase space.We believe that there is much room for further understanding, especially in view of future extensions to several (possibly massive) partons in the final state, and to higher perturbative orders.
In the present work, building on the results of Ref. [54], we tackle the problem of analytic integration of local subtraction counterterms; in the context of other NNLO subtraction schemes, this issue was addressed in Refs.[55][56][57][58][59][60][61][62].To be more precise, we note that the definition of a set of infrared counterterms has two main ingredients.On the one hand, these local functions in the radiative phase space must, in all unresolved limits, reproduce the factorised soft and collinear kernels which emerge in QCD at leading power in the soft-parton energies and in the collinear-parton transverse momenta.On the other hand, phase space itself must be factorised and parametrised so that the integration over the radiative degrees of freedom can be completely decoupled from the integration over the Born configurations: only when this step has been successfully performed can one claim the universality of the resulting subtraction algorithm.The necessary mappings of phase space have been extensively discussed in Ref. [63]: many choices are possible, and this choice is a crucial ingredient of any subtraction procedure.
Let us consider more carefully the interplay between the choice of infrared counterterms and the choice of phase-space mapping.Any QCD (squared) amplitude with the emission of one or more unresolved partons can be written as a product (to be understood as a matrix product in the colour and helicity spaces) of the (squared) amplitude for the process without the emissions, times a soft or a collinear kernel, containing all dependence on the momenta of the unresolved radiated particle(s).Any definition of subtraction counterterms must have the same factorised structure, and the kernels defined by the counterterms must reproduce the kernels of the QCD factorisation formulae, in all singular regions.Quite naturally, therefore, the first naïve choice is to use in the counterterms the kernels of the QCD factorisation formulae themselves.This is, for example, the case for FKS subtraction [4] at NLO, and for the Colourful subtraction scheme [8] at NNLO.Other well-known choices are the dipoles in the CS subtraction scheme [5] at NLO, and the antennas in the Antenna subtraction method [7] at NNLO.The CS and antenna kernels have expressions that are more involved than the ones of the QCD factorisation formulae, but still reduce to the QCD soft or collinear kernels in all singular limits.When it comes to the choice of phase-space mappings, the FKS and Colourful methods essentially involve the momenta of all outgoing particles of the radiative process, producing rather complicated expressions in the phase-space of the radiated particles, which then need to be integrated in d-dimensions.As a consequence, in the Colourful approach, these expressions in some cases can be integrated just numerically (this is not the case in FKS, because of the simplicity of NLO kernels).An easier solution for the phase-space mappings is the one adopted in the CS and antenna subtractions, where the only momenta involved in the parametrisation are the ones contained in the kernels.This choice overcomes the complexity of the latter (which in these subtraction procedures are more complicated than the QCD soft and collinear kernels), and allows for their analytical integration in the radiative phase space.
In what follows, we pursue a different approach, recently proposed in Ref. [54], which combines a definition of the counterterm kernels as close as possible to the QCD soft and collinear kernels (as is the case for the FKS and Colourful methods), together with phase-space mappings involving only the particles present in the particular kernel being integrated (as is the case for the CS and Antenna subtraction methods).As was shown already in preliminary tests performed in Ref. [54], this approach leads to simpler integrals, that can readily be computed analytically with conventional methods.The goal of this paper is thus to present the analytic integration in d-dimensions of the soft and collinear kernels of QCD factorisation formulae at NLO and NNLO, once a specific choice of phase-space mappings, along the lines of Ref. [54], is adopted.We emphasise that the results we present have a universal aspect: the full integration of NNLO QCD kernels with an exact factorisation of the radiation phase space, such that the on-shellness of the underlying Born configuration is ensured, and momentum conservation is properly enforced.On the other hand, these integrals are essential building blocks for the subtraction procedure of Ref. [54]: indeed, all required integrals for a complete subtraction algorithm for massless final-state partons are either contained in the results presented here, or are significantly simpler than the ones we perform.
The structure of the paper is as follows: in Section 2, for clarity and completeness, we present the exact integration of NLO soft and collinear kernels, which was discussed already in Ref. [54].In Section 3, we turn to the integration of tree-level kernels for double-unresolved radiation, considering explicitly double-soft emission and the case of three partons becoming collinear (which we describe as 'double-collinear' limit).In both cases, we consider un-ordered emissions, where both partons involved become unresolved at the same rate.We emphasise that this is the most intricate configuration in view of integration: hierarchical limits, with one of the two partons becoming unresolved at a higher rate than the other one, lead to a subset of the integrals considered here; similarly, nested soft-collinear limits lead to simpler integrals.In Section 4, we tackle the problem of real-virtual corrections, and integrate the QCD kernels for single real radiation at one-loop.In the process, we display the non-trivial cancellation of all singularities proportional to colour tripoles, which is an essential consistency check, given the absence of such singularities in double-virtual and double-real contributions.Finally, in Section 5, we summarise our results and present perspectives for future work.A number of technical details, including a thorough discussion of the phase-space mappings that we employ, and the treatment of integrals with non-trivial azimuthal dependence, are discussed in the Appendices.

Tree-level infrared kernels with one real emission
In this section we recall methods and results for the integration of the tree-level factorisation kernels with a single unresolved real emission, as performed in Ref. [54], and we introduce notations that we will use in the rest of the paper.We consider a generic process with a colour-singlet initial state, producing n massless coloured particles in the final state at lowest order.We will therefore be interested in scattering amplitudes involving up to n, n + 1 and n + 2 final-state coloured particles at LO, NLO and NNLO, respectively.We will denote the sets of momenta of coloured particles by {k}, where the number of particles involved will be clear from the context.Furthermore, we will adopt the notation {k} / i for the set obtained from {k} by omitting the i-th particle, and {k} [ij] for the set obtained from {k} by removing particles i and j, and introducing in their stead a single particle with momentum k i + k j .We note from the outset that, if the set {k} involves n + 1 onshell momenta k i satisfying k 2 i = 0 and i k µ i = q µ , then the set {k} / i does not satisfy the same momentum sum outside the strict soft limit k i = 0, while in the set {k} [ij] the momentum k i + k j is off-shell outside the strict collinear limit k µ i ∝ k µ j .A crucial concern in what follows will be, therefore, to choose a parametrisation of the radiative phase space factorising a lowest-order parton configuration with n on-shell partons and enforcing momentum conservation.
We expand perturbatively the amplitude for the emission of n partons as and we will use the notation B ({k}) for the Born-level squared matrix element, B({k} At NLO, we will also need the colour-connected Born squared matrix elements, n , where we use the standard notation [5,64] for the colour-insertion operators T i , responsible for the radiation of a gluon from Born-level parton i, and the spin-connected Born squared matrix elements, B µν , obtained by stripping the spin polarisation vector of a selected parton from the Born amplitude and from its complex conjugate.In this language, the virtual correction at NLO is given by V ({k}) = 2 Re(A With these definitions, we can write the well-known factorised expressions for R({k}) in the limits where one particle becomes unresolved, as follows.Defining the Mandelstam invariants of the process as we can introduce a soft-limit operator S i , extracting the leading power of R({k}) as s im → 0, uniformly for all m = i, i.e. taking all ratios of the form s il /s im to be of order one; similarly, the collinear-limit operator C ij extracts the leading power of R({k}) as s ij → 0, with all ratios s im /s jm , for m = i, j, taken to be independent of m in the limit.Under these limits, R({k}) factorises as where the normalisation factor N 1 is given by with µ the renormalisation scale and γ E the Euler-Mascheroni constant.In order to express the infrared kernels in a compact and flavour-symmetric way, we introduce flavour Kronecker delta functions: if f i is the flavour of parton i, we define for example δ fig as δ fig = 1 if parton i is a gluon, and δ fig = 0 otherwise; in similar vein, we define δ f {q,q} ≡ δ f q + δ f q , and δ {fifj }{q q} ≡ δ fiq δ fj q + δ fi q δ fiq .The soft limit is then expressed in terms of the eikonal kernel lm , which is given by In order to characterise precisely the collinear limit for partons i and j, on the other hand, we select a massless reference vector k r , which is conveniently chosen among the momenta {k} of the outgoing particles; we then introduce ratios of Mandelstam invariants, that can be interpreted as longitudinal momentum fractions along the collinear direction, as and a transverse-momentum vector We can now write the Altarelli-Parisi kernels P ij , for collinear emissions in a generic flavour configuration, in the form represents the flavour contribution with k radiated collinear gluons (k = 0, 1, 2), and can be written explicitly as The azimuthal tensor kernel Q µν ij , on the other hand, is (2.9) The task is now to introduce a parametrisation of the (n + 1)-particle phase space in terms of n on-shell massless momenta, carrying the same total momentum as the original set of n + 1 partons, and factorising the integration over the degrees of freedom of the unresolved parton.A broad set of solutions to this problem, inspired from [5], is described in Appendix A.1, and we apply it below, with the goal of simplifying as much as possible the subsequent integration.

Phase-space mappings and integration for the soft kernel
For the eikonal kernel lm , we perform the mapping described in Appendix A.1, choosing the momenta {k a , k b , k c } differently for each term in the sum in Eq. (2.2), as (2.10) Promoting the set {k} / i (which preserves momentum conservation just in the soft limit) to the momentum-conserving set { k} (ilm) of Appendix A.1, we define the mapped soft limit of R({k}) as which manifestly satisfies the condition S i S i R = S i R, necessary to ensure a local cancellation.Eq. (2.11) can be exactly integrated in d = 4 − 2 dimensions over the radiative phase space.One writes where the soft integral depends on the kinematics of particles i, l, m only through the radiative soft function J s , with argument s(ilm) lm .Substituting the expression for the Mandelstam invariants given in Eq. (A.4), J s can be trivially calculated to all orders in , with the result 2.2 Phase-space mappings and integration for the collinear kernels For the collinear kernels, we choose the momenta {k a , k b , k c } of the mapping of Appendix A.1 in the most natural way as We promote the set {k} [ij] (where the momentum k i + k j is on-shell only in the collinear limit) to the set of on-shell momenta { k} (ijr) of Appendix A.1, and we define the mapped collinear limit of R({k}) as which can easily be shown to satisfy the locality condition Proceeding with the integration, we first notice that the azimuthal kernel Q µν ij integrates to zero [5], because of its tensor structure, taking into account that kij • k(ijr) j = 0.The remaining terms, involving the P ij kernels, can again be integrated exactly in the radiation phase space.We write where the collinear integral depends on the kinematics of particles i, j, r only through the radiative collinear functions J with argument s(ijr) jr .Using again the expression for the Mandelstam invariants given in Eq. (A.4) one finds the following results.The radiation of a collinear quark-antiquark pair gives the radiation of a collinear gluon from a quark or an antiquark, on the other hand, yields finally, the radiation of two collinear gluons from a gluon yields which completes the required NLO calculations.To be precise, in order to build a complete NLO subtraction procedure one also needs to introduce and integrate the soft-collinear kernel, extracted from the combined limits S i C ij R = C ij S i R, after introducing an appropriate mapping.This presents no further difficulties, as discussed in detail in Ref. [54].
3 Tree-level infrared kernels with two real emissions In this section we consider the integration of tree-level infrared kernels with two real emissions.We first rewrite the factorisation formulae derived in Ref. [65,66] for the emission of two soft particles and three collinear particles.Indicating with n+2 | 2 the tree-level squared matrix element for the emission of two extra partons, the general structure of the double-soft limit S ij , where both particles i and j become uniformly soft, can be written as e =i,j,c,d On the other hand, in the collinear limit C ijk , where particles i, j and k become uniformly collinear, we have the general structure In Eqs.(3.1-3.2),N 1 is given by Eq. ( 2.3), while the momentum sets {k} / i/ j and {k} [ijk] are obtained from {k} by removing k i , k j , and by combining In Eq. (3.1), furthermore, we have introduced the doubly-colour-connected Born squared matrix element n , which is multiplied times the product of two eikonal factors, defined in Eq. (2.4).In the latter expression of Eq. (3.1), we have rearranged all sums in such a way that each term features only unequal colour indices.The (singly-)colour-connected squared amplitude B cd , on the other hand, multiplies the pure NNLO soft kernel, which can be written as with the explicit expressions [66] s ic s jd + s id s jc (s ic +s jc )(s id +s jd ) .
In the collinear factorisation formula, Eq. (3.2), the collinear kernels can be organised by flavour structure as where q is a quark of flavour equal to or different from that of q; similarly, the azimuthal tensor kernel can be written as In Eqs.(3.5-3.6)we introduced δ {{faf b }fc}{q q} = δ faq δ f b q δ fc q + δ fa q δ f b q δ fcq , and, as before, the superscripts (kg) refer to the number of final-state gluons featuring in the various kernels.
The expressions for P ijk , P ijk , and ijk can be extracted from Ref. [66], and can be written as where we defined and k r , as before, is a massless reference vector, which can be chosen among the momenta of the outgoing particles.From Ref. [66] we can also obtain the expressions for the azimuthal tensor kernels , which we report in our notation for completeness.They are where all terms are proportional to azimuthal tensors of the form and, in analogy with Eq. (2.6), we defined a transverse-momentum vector where the momentum is the parent momentum of the three collinear particles.We stress that the symmetry of under exchange of i, j, k, and of all other kernels under exchange of i, j, guarantees that the kernels P ijk and Q µν ijk defined in Eqs.(3.5-3.6) are totally symmetric under permutations of i, j, and k.

Double-soft kernel
In order to integrate the double-soft kernel in Eq. (3.1) we introduce different phase-space mappings according to the number of different momenta involved in the various contributions to the kernel.
For the terms containing B cd and B cdcd , where only the four particles i, j, c, d are present, we use the mapping described in Appendix A.3.1, with the identifications For the terms with B cded involving the five particles i, j, c, d, e, we use the mapping given in Appendix A.3.2, with Finally, for the terms proportional to B cdef , we use the mapping of Appendix A.3.3, with The mapped double-soft limit of RR({k}) is then + 4 e =i,j,c,d and its integral in the n + 2 phase-space is given by + 4 e =i,j,c,d where the radiative integrals of products of eikonal kernels are defined by while the radiative integral of the pure double-soft kernel is The kinematic dependence of these integrals is described by the five radiative double-soft functions J s⊗s , J s⊗s , J s⊗s , J (qq) ss and J (gg) ss .The integrals defining J s⊗s and J s⊗s are factorised, so their calculation is trivial, and can be performed to all orders in , analogously to the case with one emission.We find and The integrals defining J (2) s⊗s , J ss , J (q q) ss are not factorised, and are thus more involved.They have been performed following the procedure described in Section 3.2, with the results (3.25)

Double-collinear kernel
In order to integrate the double-collinear kernel, we perform the phase-space mappings described in Appendix A.3.1, with the choices The mapped double-collinear limit of RR({k}) is then As was the case for the single-collinear limit at NLO, the integrals of the azimuthal tensor kernel Q µν ijk vanish because of its Lorentz structure: rad,2 q µν a = 0, for a = i, j, k , which relies on the fact that ka • k(ijkr) k = 0 for a = i, j, k.The remaining terms, featuring the P ijk kernels, can be integrated in the (n + 2)-particle phase-space, and the result can be written as where the radiative integral kr admits a flavour decomposition following from Eq. (3.5), and has a kinematic dependence described by the five radiative double-collinear functions cc , J cc , J with argument s(ijkr) kr .Here we have introduced symmetrised flavour delta functions, according to where again q is a quark of flavour equal to or different from that of q.The integration of J is the computation of the highest complexity among those presented in this paper.It can however be performed analytically following the procedure described in Section 3.2.The results are This completes the integration of the factorised kernels for tree-level double-unresolved radiation.
Once again, in order to build a complete subtraction procedure at NNLO, one needs to consider both strongly-ordered and composite limits, mixing soft and collinear configurations.For such limits, it is important to find a consistent set of phase-space mappings, which need to be mutually consistent when the relevant limits are taken, in order to guarantee a local cancellation of singularities: a procedure to do so is described in Ref. [54].When it comes to the phase-space integration, however, all the composite and strongly-ordered limits are either contained in the results we just stated, or lead to significantly simpler integrals.We have thus provided all the key ingredients necessary for the integration of local counterterms (for massless final-state partons) at NNLO.

Details of the integration procedure
In this section we describe in detail how the integration of the kernels J s⊗s , J ss , J ss , J cc , J has been performed.We note that the procedure we follow, while certainly non-trivial, does not require the deployment of advanced techniques such as integration by parts or the use of differential equations (see, for example, [67][68][69][70][71][72][73]): in this sense our method, at NNLO, allows for a complete analytic integration of all subtraction counterterms, by means of relatively simple tools.
The integration procedure is simplified by a careful analysis of the symmetries of the relevant integrals under exchanges of particle labels.When integrating in the two-body radiative phase space dΦ (abcd) rad,2 , the freedom in choosing k a , k b , k c , k d does not stem from the symmetries of the kernel itself, but from those of the four-body phase space.In particular, following Ref.[57], we note that the four-body phase space for momenta k a , k b , k c , k d is symmetric under the permutation of the four momenta, as well as under the following permutations of Mandelstam invariants: These symmetries are reflected in our parametrisations of phase space, in particular when moving from the set , y, z, φ, y , z , w }, and this is crucial to simplify the analytic integration.
In order to exploit these symmetries for the integration of the soft and collinear kernels, after assigning the momenta k a , k b , k c , k d according to the discussion of Section 3.1, we apply the following transformations: • in the terms containing 1/(s ad + s bd )/(s ad + s cd ), all permutations of the invariants s ab ↔ s cd , s ac ↔ s bd , s ad ↔ s bc are performed; • in the terms containing 1/(s ad + s cd ) (but not 1/(s ad + s bd )), the permutation k b ↔ k c is performed; • in the terms containing 1/(s bd + s cd ) (but not 1/(s ad + s bd )), the permutation k a ↔ k c is performed; • in the terms containing 1/(s ad s bd ), the partial fractioning is performed, and in the first term the permutation k a ↔ k b is applied.
• in the terms containing 1/s ad (but not 1/s bd ), the permutation k a ↔ k b is performed.
After these transformations, the denominators of all integrands feature only the following combinations of invariants:  We now detail the integration procedure, focusing on one variable at a time.In Section 3.2.1 we analyse the trivial integration over y, and the first non-trivial structure arising from the w integration.Then, the subsequent integrations over z and z are detailed in Section 3.2.2,including a discussion on how we linearise the argument of the resulting hypergeometric functions.Finally, Section 3.2.3concerns the -expansion of intermediate results, and the last integration step.

Integration on y and on the azimuthal variable w
Since in all denominators in the list (3.35) the dependence on y is factorised, the integration in the y variable is always of the form which clearly gives B(n + 2 − 2 , m + 2 − 2 ).We now switch to the integration over the azimuthal variable w .According to Eq. (3.36), the only denominator containing the azimuthal variable w is s bd , while the presence of the w in the numerator uniquely stems from linear combinations of s ad and s bd , see Eq. (A. 19).Thus, terms without s bd in the denominator are of the form Terms containing the ratio s ad /s bd can be simplified according to therefore, no dependence on w in the numerator is left in the presence of the denominator s bd .The only non-trivial integration involving the azimuthal variable w is then with A = y z (1 − z) and B = z(1 − z ).Note that, as already discussed at the beginning of this section (see Eq. (3.37)), the y dependence is trivially factorised.Therefore, from now on, we understand the y dependence to be integrated out.
The integral I w is of the type described in Appendix B.1, with a = 1 and b = 1 + .From Eq. (B.10) we get then

Integration of the variables z and z
After integrating over y and w , one is left with three integrations over the variables z, z and y .We now analyse the z and the z integrations.While all numerators have a polynomial dependence on z, z , the denominators manifest a richer structure.In particular, • the invariants s ab , s ac , s bc , s cd , s ac + s bc feature a trivial dependence on z and z, as they are just products of powers of z , (1 − z ), z and (1 − z); • the structure s ad + s bd does not depend on z , while it depends on z through the factor y + z − y z; analogously, s ab + s bc depends only on z , through the factor 1 − z + z y ; • when the denominator is s bd , the z, z dependence is confined to the arguments and the prefactors of the hypergeometric functions of in Eq. (3.41), as well as in the accompanying Θ functions; the latter are to be understood as constraints on the integration region for either z or z .
The soft and collinear kernels feature products of the invariant structures described above.Among them, a non-trivial dependence on z and z arises from the following building blocks: In contributions proportional to the first structure in Eq. (3.42), the z integration gives Beta functions, while the z integration takes the form where we used 2 F 1 (a, b, c, x) = (1−x) −a 2 F 1 (a, c−b, c, −x/(1−x)).Note that m, n stand for generic powers of z, arising from the numerators.Similarly, in terms that embed the second structure in Eq. (3.42), the z integration is trivial (Beta functions), while the z integration takes the form (3.44) In the third (fourth) structure of Eq. (3.42) the whole z (z) dependence is contained in I w , and this variable is integrated first.Finally, in the fifth structure of Eq. (3.42), where no denominator depends on z nor on z , the order of integration of z and z is irrelevant.Accounting for generic numerators, whose dependence upon z and z is polynomial, we can cast all integrals to be performed as combinations of the following building blocks 1 : where n is an integer such that n ≥ −1.
Because of the symmetries of I w upon z ↔ 1 − z , the results for w z and J (n) w z can be inferred from those for I (n) w z and J (n) w z , respectively.We then proceed with the computation of the latter two integrals, which are of the type described in Eq. (B.18) of Appendix B.2 with b = 1 + .Specifically

.48)
We now show the result for specific values of n, and in particular we distinguish between n = −1 and n ≥ 0. For n = −1, Eq. (3.48) reads In some cases it is necessary to apply the partial fractioning where in the second integral of Eq. (3.49), we have inverted the argument of the hypergeometric function by means of Eq. (B.20).For n ≥ 0 the hypergeometric functions are of the class 2 F 1 (1, c + n, c, x), with c = 1 − , and can therefore be written as a finite sum in the form The integrals of Eq. (3.48) for n ≥ 0 are then given by Notice that the two integrals coincide for n = 0, evaluating to After the first z (z ) integration has been performed, all non-trivial dependence on the remaining z (z) variable is encoded in one of the following structures: where the integers n, p, q, m are such that n, p, q ≥ −1, while m = 0, 1.For later convenience, we recursively use the following partial fractioning until the condition p + q ≥ m is satisfied.
Using the symmetry under the exchange z ↔ 1 − z , the integrals I (n,p,q,m) w zz , J (n,p,q,m) w zz can equivalently be written in terms of To proceed with the computation, we choose the representation of I (n,p,q,m) w zz , J (n,p,q,m) w zz in terms of I (n) w z and J (n) w z , according to Eq. (3.53).Thanks to the results in Eq. (3.51), the case n ≥ 0 is trivial, and yields For n = −1, we exploit the integral representation of the hypergeometric functions in Eq. (3.49), introducing the auxiliary integration variable t, and write The second expression makes sense only if p ≥ 0, which is the case in all soft and collinear kernels.For m = 0, the z integration gives For m = 1, before performing the remaining z integration, we employ the partial fractioning with the results Summarising, for n ≥ 0 we still have to perform the last integration over the y variable.Conversely, for the case n = −1, we are left with two integrations, one over the physical variable y , the other over the auxiliary variable t stemming from the integral representation of hypergeometric functions.Notice that, so far, all our results are exact in : only while performing these last steps we resort to an expansion2 in powers of .

Expansion in and integration of the y and t variables
After the y, w , z and z integrations have been performed following the steps detailed in the previous sections, the integrations over y and t only involve monomials y , (1 − y ), t, (1 − t), and hypergeometric functions of the types For the first type the constraint n 3 ≥ n 1 + 1 is always achieved, thanks to the condition p + q ≥ m, which comes from the partial fractioning described in Eq. (3.54).We first manipulate these hypergeometric functions by means of the identity Since ), and can be treated recursively using the relation until it reduces to a hypergeometric function of the type We then make use of the relations (properly combined) until all hypergeometric functions are reduced to the following forms Their expansions in is known to all orders and is given by where the Spence functions S n,p (x) are defined by and reduce to standard polylogarithms for p = 1, with S n,1 (x) = Li n+1 (x).At this point, all poles in can be extracted using the standard identities where x can be either y or t.The remaining dependence does not generate any pole and can be safely expanded in Taylor series.Therefore, at this point, the remaining integrals (in t or y ) can be easily performed using standard techniques.Discarding terms that vanish in the → 0 limit, we obtain the final expressions given in Section 3.1.

One-loop infrared kernels with one real emission
To complete the study of NNLO factorisation formulae we are left with the integration of one-loop infrared kernels involving the emission of one soft or two collinear particles at the one-loop level.These kernels are known from the literature [74][75][76][77][78][79][80], and we rewrite them in the most suited form to perform their integration in the radiation phase-space in the context of our method.Indicating with RV ({k}) the renormalised one-loop squared matrix element for the emission of one unresolved parton i, the factorisation formulae for the soft limit S i and for the collinear limit C ij read where the symbols lm , P ij , and Q µν ij were already introduced in Section 2. In addition, here we have introduced the completely antisymmetric tripole-colourcorrelated Born squared matrix element as well as the colour-connected one-loop squared matrix element n , and the spin-connected one-loop squared matrix element V µν , obtained by stripping the spin polarisation vectors of the particle with momentum k i + k j from both the matrix element and its complex conjugate inside V .The one-loop soft kernels are Their collinear counterparts, on the other hand, can be written as While the one-loop kernels are rather intricate, there is only a single further unresolved radiation: the phase space mapping, to which we now turn, is therefore simpler in this case.

One-loop soft kernel and cancellation of colour tripoles
As done for the tree-level infrared kernels with one real emission, for the soft kernels we perform the mapping described in Appendix A.1, choosing the momenta {k a , k b , k c } as the momenta {k i , k l , k m } present in each term of the eikonal kernel, according to Promoting the set {k} / i to the momentum conserving set { k} (ilm) of Appendix A.1, we define the mapped soft limit of RV ({k}) as The integral of this function in the (n + 1)-particle phase-space can be written as where J ilm s = δ fig J s , defined and computed in Eqs.(2.13-2.14),must here be expanded up to order O( 2 ).One gets Also the integral Jilm s , defined below, can be easily computed after substituting the expression for the Mandelstam invariants in our parametrisation, Eq. (A.4).The result is whose kinematic dependence is described by the one-loop radiative soft function Js given by We now discuss the last and most interesting contribution to Eq. (4.9): the soft integral proportional to the triple-colour-correlated Born amplitude.It is defined by mp .As can be guessed from the more intricate kinematic dependence, this part of the soft one-loop kernel requires more refined techniques to be analytically integrated: the reason is its peculiar kinematic structure, involving two eikonal kernels linking four particles, which leads to a non-trivial azimuthal dependence.With the phase-space mapping {k a , k b , k c } → {k i , k l , k m }, we can use the results of Appendix A.2 to parametrise the Mandelstam invariants present in lmp in the form since the first and second terms vanish separately upon summation over the indices m and l, respectively.For the remaining structure we get where individual terms vanish thanks to the same symmetry arguments used in Eq. (4.19).This completes the proof that colour tripoles do not contribute to infrared counterterms at NNLO, except for subtraction-scheme-dependent finite contributions.In our approach, these are given by l =i,m =i,l p =i,l,m For the one-loop collinear kernel we choose the momenta {k a , k b , k c } of the phase-space mapping in Appendix A.1, as was done for the tree-level collinear kernel with one real emission.Thus we pick We promote the set {k} [ij] to the set of on-shell momenta { k} (ijr) of Appendix A.1, and get Once again, the integration of the collinear kernels is simplified by the fact that terms proportional to Q µν ij and Q µν ij integrate to zero, because of their Lorentz structure.For the remaining pieces, containing P ij and P ij , integration in the (n + 1)-particle phase-space leads to The integral J ijr c is defined in Eq. (2.18), where it is expressed in terms of the radiative collinear functions J (kg) c (with k = 0, 1, 2).These must now be computed to O( 2 ), and yield where the kinematic dependence can be described by the one-loop radiative collinear functions where n, m can take only the integer values −1, 0, 1.For these values, the integral can be expressed in terms of a generalised hypergeometric function of type 3 F 2 , evaluated at unit argument.More precisely, This, in turn, can be expanded in powers of , using for example the package HypExp [88,89].The integration over the remaining radiative phase-space variables is then straightforward.The final results for the three contributions to Jijr   This completes the list of all the integrals associated with factorised soft and collinear kernels at NNLO.These integrals form the basis for the construction of all integrated infrared counterterms for single-and double-unresolved real radiation at NNLO.

Conclusions
In any massless gauge theory, (squared) matrix elements factorise in soft and collinear limits, at leading power in the soft energy and in the small transverse momentum, yielding universal soft and collinear kernels, which multiply the (squared) matrix element for the Born process, without the unresolved particles.Away from the strict limits (or beyond leading power in the resolving variables) this factorisation is not exact: in particular, the factorised Born matrix element does not conserve momentum (near the soft limit), or is not on the mass shell (near collinear limits).In order to integrate the factorisation kernels over the unresolved degrees of freedom in a universal way (i.e.requiring no information on the underlying Born process), one needs to provide a set of phase-space mappings, which must re-express the factorised Born process in terms of an on-shell, momentumconserving set of momenta.This amounts to a specific choice of a set of sub-leading power terms in the factorisation, and such a choice is a necessary ingredient for any infrared subtraction procedure.
In the present paper, we have presented the complete integration of the QCD factorisation kernels at NLO and NNLO, with a set of phase-space mappings selected along the lines suggested in Ref. [54], chosen with the goal of simplifying as much as possible the analytic integration.As a consequence, we have been able to give analytic results for all kernels, including non-singular terms.In particular, all integrals of the double-real counterterms can be written exactly, to all orders in , in terms of hypergeometric functions, with the most intricate cases involving 4 F 3 evaluated at unit argument.We have however chosen to give the expansion of these hypergeometrics in powers of , up to and including O( 0 ) terms, since this is what is required in practical calculations.All our results have been validated against independent numerical integration codes based on sector decomposition [90][91][92].The analytic results of this paper are necessary (and indeed sufficient) ingredients to build all integrated counterterms in the context of the local analytic sector subtraction of Ref. [54].The present work shows that this novel approach allows to use standard techniques to compute an important class of integrals that appear in all NNLO QCD computations, yielding comparatively very simple results.
We believe that achieving the maximum simplicity in the case at hand -NNLO radiation of massless partons in the final state -is important not only for building an efficient and transparent NNLO subtraction algorithm for these processes, but also for future extensions.The method presented here is expected to be generalisable to initial-state QCD radiation without conceptual changes, and the results are sufficiently simple that a generalisation to massive particles at NNLO appears feasible.Furthermore, since the integrations presented in this paper have been performed with conventional techniques, one may reasonably hope that more advanced techniques, such as those involving differential equations for Feynman integrals (see, for example, [67][68][69][70][71][72][73]), might be sufficient to tackle the problem even at the next perturbative order.
and the w integration has the structure of the master integral of Appendix B.1, with A 2 = Cv and B 2 = D(1 − v).Thus one may write The step functions appearing in Eq. (B.13) modify the v integration domain as Since C, D > 0, then 0 < D C+D < 1 and we get In the integration of the two-unresolved tree-level kernels the integration over the azimuthal angle gives rise to the master integral I a,b,β,γ (C, D) with a = 1 (see Section 3.2.1),which deserves a separate analysis.We define The hypergeometric functions with the first argument set to unity satisfy Then we can write the first hypergeometric function using its integral representation, as Once the integer part of the first index is 0, we can then expand in powers of using and then easily perform the remaining integrations.

. 47 )
We see that the integral I 1+ ,− ,n− (1 − z , y z ) is of the special type I b,1−b,γ (C, D) described in Eq. (B.30), while the integral I 1+ ,n− ,− (1 − z , y z ) is of the special type I b,β,1−b (C, D) described in Eq. (B.31).Using the results derived there we have jr .Terms in P ij which are proportional to the simple polynomials N (1g) ij and N (2g) ij (see Eq. (4.4)) can be integrated easily.Less trivial integrals arise from the P ij M ij term in Eq. (4.4), and in particular from structures of the type momenta are left unchanged ( k(acd,bed) n = k n , n = a, b, c, d, e