Four-dimensional formulation of the sector-improved residue subtraction scheme

Four years ago, one of us introduced a novel subtraction scheme for the evaluation of double-real radiation contributions to cross sections at next-to-next-to-leading order (NNLO) in QCD. This approach, named SecToR Improved Phase sPacE for Real radiation (STRIPPER), has already found several non-trivial applications. In particular, it has allowed for the determination of NNLO corrections to hadronic top-quark pair production, fully differential top-quark decays, inclusive semileptonic charmless b-quark decays, associated Higgs boson and jet production in gluon fusion, muon decay spin asymmetry, and t-channel single-top production. Common to these calculations was the use of conventional dimensional regularization (CDR). In this publication, we present a complete formulation of the subtraction scheme for arbitrary processes with any number of colored partons in the final state, and up to two partons in the initial state. Furthermore, we modify the integrated subtraction terms of the double-real radiation to enable the introduction of the 't Hooft-Veltman version of dimensional regularization (HV), in which resolved states are four-dimensional. We demonstrate the correctness of our approach on the example of top-quark pair production in the gluon fusion channel.


Introduction
Perturbative calculations beyond the next-to-leading order of Quantum Chromodynamics are notorious for their complexity. Already the next-to-next-to-leading order presents tremendous obstacles. One of them is the evaluation of infrared and ultraviolet divergent two-loop virtual amplitudes. This is an unsolved issue in the general case despite the existence of a few analytic results for low multiplicity processes. Another obstacle is the evaluation of double-real and mixed real-virtual contributions containing infrared singular phase space integrals. In principle, these contributions can be obtained by means of Monte Carlo techniques, once suitable subtractions have been introduced to generate numerically integrable functions. The form of the subtractions defines a subtraction scheme. Unfortunately, it turns out that subtraction schemes are extremely complex. At present, there are several on-going multi-year efforts at the construction of general solutions. Antenna subtraction [2], and q T subtraction [3], are amongst the most advanced initiatives, which have already found several non-trivial applications. Another scheme under construction has been introduced in Ref. [4]. Here, we will be concerned with the subtraction scheme STRIPPER introduced by one of us in Ref. [1]. Inspired in some aspects by ideas of Frixione, Kunszt and Signer [5], and in other aspects by ideas of Binoth and Heinrich [6], the scheme has proven its worth in many applications: it has allowed for the determination of NNLO corrections to hadronic top-quark pair production [7,8,9,10,11], fully differential top-quark decays [12], inclusive semileptonic charmless b-quark decays [13], associated Higgs boson and jet production in gluon fusion [14], muon decay spin asymmetry [15], and t-channel single-top production [16]. The listed advanced STRIPPER applications performed independently of the inventor have been preceded by the much simpler case of QED corrections to Z-boson decay into a pair of massless leptons [17]. The purpose of the present work is to complete the construction of the scheme in order to allow for the evaluation of cross sections for arbitrary processes.
The original idea of Ref. [1] was to concentrate on the numerical calculation of the coefficients of the Laurent expansion in (dimensional regularization parameter, with space-time dimension d = 4 − 2 ) of the double-real cross section contribution. The latter requires a phase space integral over the momenta not only of the partons present in the Born approximation, but also of two additional massless partons. It seemed obvious that other cross section contributions are much easier to obtain, since their kinematics is, in the worst case, the same as that of next-to-leading order real-radiation contributions. Furthermore, the concept was to refrain from (almost) any analytic integration. By inspection of other efforts, it was clear that it is the insistence on analytic integration that makes the subtraction schemes difficult to develop. Finally, in order to avoid as many complications as possible, the construction was performed uniformly in d dimensions. This corresponds to conventional dimensional regularization, where both momenta and spin degrees-of-freedom of external particles are d dimensional. This implies that even Born matrix elements will have a non-trivial expansion in . While the other basic ideas of the scheme stood the test of time, the use of CDR is now an important drawback. Indeed, software implementations of tree-level matrix elements only provide them at = 0, i.e. in four dimensions. This makes it necessary to recalculate the matrix elements for each project from scratch. Furthermore, the need to parameterize an increasing number of dimensions depending on the multiplicity of the process seems not only a major annoyance, but also a source of inefficiency. In this publication, we will solve the problem by introducing a number of corrections to the integrated subtraction terms in the double-real radiation contribution. As a result, we will formulate the scheme in 't Hooft-Veltman regularization. In particular, we will only need four-dimensional external momenta and polarizations in the evaluation of actual matrix elements. There will still be a trace of higher dimensionality in the integration over unresolved momenta. However, we will only have to consider six dimensional unresolved momenta in the worst case. On the example of top-quark pair production, we will also demonstrate that the introduced improvements of STRIPPER fulfill their purpose. We stress that we only present the algorithm to obtain the expressions needed for the implementation of the subtraction scheme. This algorithm requires the knowledge of soft and collinear limits of QCD amplitudes and provides the expressions for the subtraction and integrated subtraction terms by simple substitutions. Due to the number and size of the resulting final formulae, we do not reproduce them here. Instead, we plan to provide a software package in the future. Nevertheless, existing calculations can be converted to 't Hooft-Veltman regularization with moderate effort following this publication.
The paper is organized as follows. In the next section, we present an outline of the subtraction scheme together with a complete list of cross section contributions to evaluate in a general next-to-next-to-leading order calculation. Subsequently, we proceed with the construction of the scheme by considering its major elements: phase space decomposition (Section 3), phase space parameterization (Section 4), and derivation of the subtraction and integrated subtraction terms (Section 5). This part is common to the CDR and HV formulations. Afterwards, we discuss the modifications necessary to introduce HV regularization: average over the azimuthal angles (Section 6), separation of finite contributions (Section 7), and finally the regularization itself (Section 8). This part is followed by the example of top-quark pair production in the gluon fusion channel (Section 9). The main text is closed with conclusions. Appendices form an important part of the publication. They are self-contained, and the information provided therein is sufficient to construct the complete subtraction scheme in the most general case. Our notation is summarized in Appendix A. Spherical angles in d dimensions are discussed in Appendix B. Divergences of virtual amplitudes are considered in Appendix C, whereas limits of amplitudes are contained in Appendix D and Appendix E. Appendix F contains the well-known Altarelli-Parisi splitting kernels necessary for initial-state collinear renormalization of partonic cross sections.

Outline of the subtraction scheme
We are interested in the calculation of next-to-next-to-leading order QCD corrections to scattering cross sections. The subtraction scheme we will describe can also be used for decay processes, since the phase space integrals present there share the same structure. The formulae below would only have to be trivially modified. Furthermore, it is possible to include QED corrections by a modification of the color algebra (see Appendix A). We will not discuss this issue in any more detail. We now need to specify the possible initial and final states. At present, the most relevant case is hadron-hadron scattering. Therefore, we will give expressions with this class of processes in mind. Nevertheless, we will separate the hadron scale physics from the parton scale physics. Indeed, we will only ever manipulate partonic cross sections. This makes it possible to include lepton-hadron and lepton-lepton scattering processes by simply removing factorization contributions. As far as the final state is concerned, we only require that there be at least two particles at the Born level. The reason is that in the special case of a single particle, the partonic cross section becomes a distribution, and we assume that we can work with ordinary functions only. In the following, we will concentrate on external partons carrying a color charge. Any additional particles, which are not charged under the color gauge group play no role whatsoever in the construction of the scheme, and can be trivially included in the phase space integrals.
The hadronic cross section is known to factorize into a convolution of parton distribution functions and the renormalized partonic cross section σ h 1 h 2 (P 1 , P 2 ) = ab 1 0 dx 1 dx 2 f a/h 1 (x 1 where P 1,2 are the momenta of the hadrons h 1,2 , and f a/h (x, µ 2 F ) is the parton distribution function (PDF) of parton a within the hadron h, at the factorization scale µ F . In the following, p 1,2 = x 1,2 P 1,2 will denote the parton momenta. If we are able to numerically evaluate the partonic cross section,σ ab , for an arbitrary initial-state energy, then we can also obtain the hadronic cross section by an additional integration over the parton momentum fractions. This is the reason, why the construction of the subtraction scheme will never make reference to the initial-state hadrons, but only to the initial-state partons. The partonic cross section can be expanded in a series in α s (µ 2 R ), where µ R is the renormalization scale. Up to next-to-next-to-leading order (NNLO) it readŝ σ ab =σ (0) ab +σ (1) ab +σ (2) ab .
Each term of the expansion can be further decomposed according to the multiplicity of the final state. At leading order whereŝ = (p 1 + p 2 ) 2 is the square of the partonic center-of-mass energy, while N ab is the spin and color average factor, defined as the product of the number of spin and color degrees of freedom of the partons a and b. The subscript n points to the number of final states in this contribution. The notation and normalization of phase spaces and matrix elements is specified in Appendix A. F n is the measurement function defining the observable. It is a function of the final state momenta and flavors. Up to the next-to-next-to-leading order, we need three such functions, F n , F n+1 and F n+2 . They must fulfill the requirements of infrared safety, i.e. if the energy of a final state parton vanishes, or if the parton becomes collinear to another one, then F n+m → F n+m−1 for m > 0 and F n → 0. At next-to-leading order there isσ (1) ab =σ R ab +σ V ab +σ C ab , where the splitting functions required are reproduced in Appendix F. Finally, at next-to-next-to-leading order, we haveσ (2) ab =σ RR ab +σ RV ab +σ VV ab +σ C1 ab +σ C2 ab , andσ C1 ab (p 1 , p 2 ) = α s 2π where In general, it is only possible to evaluate the cross section contributions listed above using numerical Monte Carlo integration methods. Unfortunately, phase spaces with n + 1 and n + 2 final states contain infrared singular configurations, which are regulated with dimensional regularization. The latter, if applied naively, spoils the stability of the numerics, since a limit has to be taken once all contributions are combined. This problem is resolved with a subtraction scheme, which extracts the explicit singularities and provides integrable functions, which do not depend on the parameter of dimensional regularization.
The idea behind the construction of the subtraction scheme STRIPPER is to derive Laurent expansions in for each of the cross section contributions independently. The basic algorithm has three stages: phase space decomposition (Section 3): phase spaces with n + 1 and n + 2 final-state particles are decomposed into sectors, in which only certain types of singularities may occur; phase space parameterization (Section 4): in each sector, a special parameterization is introduced using spherical coordinates in d dimensions (Appendix B), in which singularities are only parameterized with 2 variables for n + 1-final-state-particles phase spaces, and with 4 variables for n + 2-final-state-particles phase spaces; generation of subtraction and integrated subtraction terms (Section 5): in each parameterization, subtraction terms in the relevant variables are introduced, which make it possible to obtain an expansion in , the coefficients of which are integrable. The subtraction terms only require the knowledge of the singular limits of QCD amplitudes (Appendix D and Appendix E), and are process independent in the sense that the process dependence is confined to the matrix elements. Furthermore, pointwise convergence of phase space integration is guaranteed.
After application of this algorithm, cross sections may in principle be evaluated numerically. Nevertheless, since dimensional regularization involves infinite dimensional vectors, the effective dimension of the vectors, which actually occur in the calculation increases with multiplicity. In fact, any new vector requires an increase of the effective dimension by one. For two-to-two processes at leading order, one already needs five dimensions at next-to-next-toleading order. Furthermore, matrix elements must be provided as expansions in . In order to simplify the calculation, we introduce the 't Hooft-Veltman version of dimensional regularization, in which resolved particle momenta and spin degrees-of-freedom are four-dimensional. In this case, we also only need four-dimensional matrix elements. The construction proceeds in additional three stages: average over azimuthal angles (Section 6): integrated subtraction terms, which have been derived in relation to a collinear limit, are averaged over the unphysical transverse direction. For most cases, this is equivalent to the use of averaged splitting functions, but there are important exceptions. This step is important in order not to have contractions of four-dimensional matrix elements with d-dimensional transverse vectors; separation of finite contributions (Section 7): the different contributions listed in this section are further decomposed into classes with different kinematics and loop order. The sum of the terms in each class is finite. This requires a modification of the integrated subtraction terms for the double-real radiation. In practice, counterterms are introduced, which are added to one class of contributions and subtracted from another; 't Hooft-Veltman regularization of separately finite contributions (Section 8): the measurement function is modified to contain delta-functions restricting the momenta of resolved particles to be four-dimensional. For most classes of finite contributions, this is already sufficient to fulfill the requirements of 't Hooft-Veltman regularization, and the matrix elements can be evaluated in four dimensions. Nevertheless, one class, the singleunresolved contributions to double-real radiation, requires a further modification of the integrated subtraction terms.
After these steps, the subtraction scheme does not require higher orders of the -expansion of the matrix elements, and all resolved momenta are four-dimensional. The calculation still involves unresolved momenta, which may need up to two additional dimensions, but only occur in soft and splitting functions. In general, two-to-two processes at leading order require five-dimensional unresolved momenta at next-to-next-to-leading order. For higher multiplicity, six-dimensional momenta must be introduced.

Phase space decomposition
The first step in the construction of the subtraction scheme is the decomposition of the phase space into sectors. In each of the sectors, the momenta of selected partons are parameterized with a few relevant variables (see Section 4), which allows for easy extraction of explicit poles in the dimensional regularization parameter , and generation of numerically integrable coefficients of the Laurent series (see Section 5). The decomposition is performed with the help of selector functions. The only essential property required of the selector functions, besides that they form a decomposition of unity, is that they vanish in case of soft and/or collinear limits not covered by the parameterization of a given sector. Cross sections are defined with a measurement function, which specifies an infrared safe observable. This measurement function reduces the number of possible singular configurations. For instance, at next-to-leading order, the measurement function F n+1 only allows a parton to become soft (vanishing energy), two massless partons to become collinear to each other (vanishing relative angle), or a parton to become both soft and collinear to another massless parton. This is achieved by the vanishing of F n+1 for any other singular configuration. F n+1 also enters the definition of cross sections at next-to-next-to-leading order, if they involve n + 1 partons in the final state. The selector functions for n + 1-parton phase space integrals are used to regulate the singularities allowed by F n+1 . Since at most a pair of partons may be involved, the selector functions for this case, S i,k , have two indices. The first index, i, specifies the unresolved parton, which is massless and in the final state, while the second index, k i, specifies the reference parton, which is massless as well, but may be either in the initial or in the final state. S i,k does not vanish if the unresolved and reference partons become collinear to each other, or if the unresolved parton becomes soft. However, it vanishes if any gluon besides i becomes soft, or if any two partons besides i and k become collinear and this kinematical configuration is singular. In consequence, if the reference parton is in the initial state, we only have to consider the following flavor pairs, ( f i , f k ), of the flavor f i of parton i and flavor f k of parton k (g, g), (g, q), (g,q), (q, g), (q, g), (q, q), (q,q) .
If the reference parton, on the other hand, is in the final state, then the number of possible pairs is reduced, because we may use the symmetry between i and k. Nevertheless, if the unresolved-reference-parton pair contains a gluon, the latter has to be the unresolved parton. The list is thus (g, g), (g, q), (g,q), (q,q) .
While the choice of S i,k is by no means unique, we have to provide a valid set. We choose a form very similar to the one described in [18]. To this end, we introduce where α = γ for a gluon and α = 0 for a quark, while β, γ > 0. E i is the energy of parton i, and θ ik is the relative angle between partons i and k. With the help of d i,k , we define Clearly, S i,k form a decomposition of unity We note that the contributions of two sectors, which have the same flavors of the unresolved and reference partons, are equal. This implies that there is no need to calculate the contributions of such sectors independently. In consequence, even for high-multiplicity processes, the number of independent sectors is moderate. For example, the most computationally intensive pure gluon amplitudes in multi-jet production will only require three sectors independently of the multiplicity. The next-to-next-to-leading order measurement function occurring in double-real radiation, F n+2 , allows for more singular cases than F n+1 . In particular, it does not vanish if three partons, or up to two pairs of partons, become collinear. Similarly, it does not vanish if up to two partons become soft. In other words, it allows for two unresolved partons. In order to accommodate the possible triple-collinear (three partons becoming collinear) and doublecollinear (two pairs of partons becoming collinear) limits, we must introduce two types of selector functions: S i j,k for unresolved partons i and j with reference parton k; and S i,k; j,l for unresolved parton i with reference parton k, and unresolved parton j with reference parton l. Of course, all indices must be different in both cases, and the ordering of i and j in S i j,k , as well as the ordering of the two pairs, (i, k) and ( j, l) in S i,k; j,l is irrelevant. In the triple-collinear sector, the flavor sets allowing for singular kinematical configurations are {g, g, g}, {g, g, q}, {g, g,q} {g, q,q}, {q,q, q } , where q can be any quark or anti-quark including q andq, and we have used crossing to define the flavors of all partons as if they were out-going. If k is in the initial state, the possible flavor triples, ( f i , f j , f k ), we have to consider, are obtained by selecting any of the different flavors in each set of the list (15) and crossing it to the initial state. This procedure results in the following list (g, g, g), (g, g, q), (g, g,q), (g, q, g), (g,q, g), (g, q, q), (g,q,q), (q,q, g), (q,q, q ), (q , q, q), (q ,q,q) .
If k is in the final state, on the other hand, we have additional freedom to reorder the partons within the triple. Minding the necessity to regulate final state soft gluons and double-soft quark-anti-quark pairs, we obtain the following list (g, g, g), (g, g, q), (g, g,q), (g, q,q), (q,q, g), (q,q, q ) .
The flavor assignments we have to consider in the double-collinear sector are simply a composition of two next-toleading order cases for the two involved pairs. Let us now introduce where α i, j = γ for gluons, α i = α j = γ for a quark-anti-quark pair, and α i, j = 0 otherwise. We define Note that in the case of S i,k; j,l , if (i, j) is a quark-anti-quark pair, we set α = γ in both d i,k and d j,l , contrary to the next-to-leading order case. The decomposition of unity is The number of sectors to actually evaluate is again reduced by the fact that sectors with identical flavors give equal contributions. Returning to the example of pure gluon amplitudes, we would always only have to calculate double-real radiation contributions of seven sectors: three triple-collinear and four double-collinear.
We must finally point out that the sectors we have introduced are sufficient for any problem with at least three massless partons. In the case of only two massless partons (say heavy-quark pair-production from colorless initial states) in the final state, it is necessary to introduce fictitious reference momenta.

Phase space parameterization
After having decomposed the phase space into sectors using selector functions of Section 3, we now need to introduce appropriate parameterizations for the single-collinear sector (one unresolved momentum), and the triple-and doublecollinear sectors (two unresolved momenta). The phase space integrals always assume the form where n is the number of final state momenta in the Born approximation, while n u is the number of unresolved momenta, and n f r the number of final state reference momenta. Q is the total momentum of the remaining final-state particles. In particular, the phase space dΦ n−n f r (Q) describes the decay process of a state with momentum Q into n − n f r particles with momenta q i Q → q 1 + · · · + q n−n f r .
For future reference, we define the notation for the minimal invariant mass of the final state where m i is mass of the final state particle i in the Born approximation. Furthermore, we introduce the related maximal energy of a massless final state parton whereŝ is the square of the partonic center-of-mass energy. The parameterizations are constructed with several rules in mind. In particular, the phase space is defined in the partonic center-of-mass frame in all sectors but the single-collinear sector for use in the calculation of the factorization contributionσ C1 . The initial state momenta, p 1 and p 2 are always in the center-of-mass system. The unresolved momenta are parameterized by rescaled energy and angle variables, to allow for a straightforward identification of the singular limits of amplitudes. The unresolved phase space integrals contain for a single unresolved momentum, and for two unresolved momenta. In general η are related to angles, and ξ to energies. The factors x a−b , x ∈ {η, ξ} or x ∈ {η 1 , η 2 , ξ 1 , ξ 2 }, regulate singularities at x = 0. This requires a i ≥ 0 and b i > 0. The latter restriction follows from the fact that infrared divergences are regulated by < 0. If the product of a matrix element and the measurement function is integrated with the given phase space parameterization together with the appropriate selector function, there are no other singularities than those located at x = 0, and the singular behavior in the limit is described by the factor x −1−b , i.e. the coefficient of this factor is finite at the limit for any value of . The method of proof of this fact has been discussed in [1], and the reader is referred to that publication for details. The parameterizations make ample use of the rotation invariance of the single-particle measure. They are valid in CDR with d-dimensional momenta of both reference and unresolved momenta. Nevertheless, if reference momenta and additional final state particle momenta are restricted to four dimensions, then one unresolved momentum is at most five dimensional, while a second at most six dimensional. This statement should be interpreted in the sense that the additional dimensions can be trivially integrated out, since nothing depends on them.
A final comment is necessary on the order of integration. The phase spaces are written in such a way that every next integral inherits the parameters of the previous one. In other words, the integration range of a given integral is restricted depending on the particular values of the parameters of all the previous integrals. In order to achieve a pointwise cancellation of singularities during Monte Carlo integration, as discussed in the next section, it is necessary to parameterize dΦ n−n f r (Q) in a specific way. This phase space depends on 3(n − n f r ) − 4 parameters in four dimensions. These parameters must be rescaled to a fixed finite (usually unit) range, as natural for Monte Carlo integration. The rescaling itself depends on Q. During the evaluation, the rescaled parameters are kept fixed for both the unsubtracted matrix element with phase space weight, and for all the subtraction terms with their phase space weights. The same comment also applies to the reference parton energy.
Let us now specify the parameterizations using the notation of Appendix B for vectors and integrals in spherical coordinates in d-dimensions.
4.1. Single-collinear sector parameterization: one reference momentum, one unresolved momentum The resolved and unresolved parton momenta read with the angular parameterizationr For reasons explained in Section 7.4, we will allow for a boosted frame, where the initial state momentum p 1 is rescaled with z. Of course, the symmetric case, where the rescaling is applied to p 2 can be treated by relabeling. In most applications, we will set z = 1. In the case of an initial state reference momentum, the phase space is while for a final state reference momentum where and in the particular case n = 2 r 0 max 0 dr 0 (r 0 ) 1−2 2(2π) 3−2 dΦ 1 (zp 1 Figure 1: Example momentum parameterizations in the case of two unresolved momenta, u 1 and u 2 , and one reference momentum, r, which may be either in the final state (left), or in the initial state (right, with r = p 1 ). The angles θ 1 and θ 2 allow for a straightforward parameterization of the collinear limits with respect to the reference momentum. Q is the total momentum of the remaining final-state particles, while p 1 and p 2 are the initial-state momenta.
The unresolved phase space reads where The behavior of amplitudes in the collinear limit is characterized by the transverse vector 4.2. Triple-collinear sector parameterization: one reference momentum, two unresolved momenta The resolved and unresolved parton momenta read with the angular parameterization (see Fig. 1) 1 (α 1 , α 2 , . . . )n (3−2 ) (θ 1 , φ 1 , ρ 1 , ρ 2 , . . . ) , In the case of an initial state reference momentum, the phase space is while for a final state reference momentum where and in the particular case n = 2 The parameters u 0 1,2 , θ 1,2 , and φ 2 of the unresolved momenta are replaced byξ 1,2 ,η 1,2 , and ζ The unresolved phase space is split according to the ordering of the energies of the unresolved partons Each of the two resulting contributions is further decomposed according to Fig. 2 into five sectors where η 1,2 , ξ 1,2 parameterizeη 1,2 ,ξ 1,2 as in Tab. 1, while µ S i can be found in Tab. 2, with , .
(45) The essential property of these functions is that they do not vanish in any of the singular limits of amplitudes indicated by the vanishing of any of η 1,2 , ξ 1,2 . The behavior of amplitudes in collinear limits is characterized by the transverse vectors Table 1: Original kinematic variables of the triple-collinear sector parameterization,η 1 ,η 2 ,ξ 1 ,ξ 2 , expressed through the sector variables, η 1 , η 2 , ξ 1 , ξ 2 , of the five sectors, S 1 , . . . , S 5 , defined in Fig. 2. The function ξ 2 max is defined in Eq. (45). Figure 3: Example momentum parameterizations in the case of two unresolved momenta, u 1 and u 2 , and two reference momenta, r 1 and r 2 . The left picture shows a case with initial state reference momenta with r 1 = p 2 and r 2 = p 1 , while the right picture shows a mixed case with one final state reference momentum, r 1 , and one initial state reference momentum, r 2 = p 1 . The angles θ 1 and θ 2 allow for a straightforward parameterization of the collinear limits with respect to the reference momenta. Q is the total momentum of the remaining final-state particles, while p 1 and p 2 are the initial-state momenta.
with tanφ ± 2 (θ 1 , ζ) = ± sin θ 1 ∂ + θ 2 φ 2 (θ 1 , ζ) ,φ + 2 ∈ 0, For our particular choice of the dependence of φ 2 on ζ, Eq. (42), there is 4.3. Double-collinear sector parameterization: two reference momenta, two unresolved momenta 4.3.1. General case with n > n f r The resolved and unresolved parton momenta read with the angular parameterization (see Fig. 3) In the case of two initial state reference momenta, the phase space is while for one final state (r 1 ) and one initial state (r 2 ) reference momentum and similarly if the rôles of r 1 and r 2 are reversed. For two final state reference momenta, there is r 0 max is defined in Eq. (40) withr =r 1 . If the second reference momentum is in the initial state and n = 2, then Eq. (41) applies. Furthermore, in the case of two final state reference momenta and in the particular case n = 3 The unresolved phase space is split according to the ordering of the energies of the unresolved partons with where Notice that if α i = β i = 0 for i > 2, thenû 1 ·û 2 is only a function of α 1,2 , β 1,2 , θ 1,2 , φ 1,2 , ρ 1 , and σ 1,2 . This fact will be important for the four-dimensional formulation of the subtraction scheme. The behavior of amplitudes in collinear limits is characterized by the transverse vectors Figure 4: Parameterization of a fully massless four-particle phase space (left), where two of the final states are used as reference vectors, r 1 and r 2 , while the other two, u 1 and u 2 , are considered to be unresolved. p 1 and p 2 are initial state momenta, while r is a reference vector equal to the sum of one reference and one unresolved momentum, r = r 1 + u 1 . If the parameters of one of the unresolved momenta are specified, then the available range of parameters of the other one is split into two disjunct regions (right). There are neither soft nor collinear singularities in region II. withû 1⊥ = lim 4.3.2. Special case with n = n u = n f r = 2 This is a parameterization of a fully massless four-particle phase space. Unlike previous cases, it uses an auxiliary reference vector, r. The configuration is depicted in Fig. 4, and defines the following reference and unresolved momenta where with the angular parameterization We start with a general discussion of the phase space. A parameterization suitable for the derivation of the subtraction and integrated subtraction terms will be provided near the end of the section. We write where Notice that in the collinear limit θ i → 0, alsoθ i → 0. This implies that the integration measure factors do not influence the scaling of the integrand in the limit. The integration region splits into two, region I and region II. We choose the energy and angle of the first parton to provide restrictions on the energy and angle of the second parton, see Fig. 4, thus Region I reads Region II readsŝ Region II does not contain soft or collinear singularities. Indeed, as long as u 0 1 > 0 and cos θ 1 < 1, the endpoints u 0 2 = 0 and cos θ 2 = 1 are not reachable. On the other hand, if either u 0 1 = 0 or cos θ 1 = 1, region II has zero volume. We are now ready to present the final parameterization of the phase space. We have The contribution on the second line represents region II. Since there are no singularities there, we will keep the formulae we have given before without any further modification. The integral can be performed in four dimensions, i.e. = 0. On the other hand, we will write where The factor in the last line of Eq. (81) is regular in all limits. The behavior of amplitudes in collinear limits is characterized by the transverse vectors

Angular integrations beyond four dimensions
The parameterizations we have presented involve d-dimensional angular integrations. In practice, the number of relevant dimensions depends on the number of vectors present in the problem. However, we will later define the scheme in 't Hooft-Veltman regularization, in which the resolved momenta (the final state momenta q i , the reference momenta r i , and possibly up to two of the unresolved momenta), are four-dimensional. Inspecting the explicit parameterizations of the unresolved momenta from the previous subsections, we notice that the scalar products amongst themselves and with four-dimensional vectors involve parameters from at most two additional dimensions. The angular integrations over the parameters, which do not occur in the integrand can be performed explicitly. If the integrand depends on four-dimensional parameters only, then we use If the integrand depends on the parameters of a single unresolved momentum, then Finally, if the integrand depends on the parameters of two unresolved momenta, then the integration over the parameters of the first is done with the previous formula, while the integration over the parameters of the second requires The integrands contain the distribution In the case of two-to-two processes, i.e. processes with n = 2, there is an additional simplification. Indeed, we can then assume that the resolved momenta in the Born approximation are three-dimensional, i.e. they have twodimensional spatial components. This corresponds to scattering on a plane. In such a case, the necessary dimension of the unresolved momenta drops by one. This implies that for n = 2, five-dimensional unresolved momenta are sufficient.

Generation of subtraction and integrated subtraction terms
Using the decomposition of Section 3 and the parameterizations of Section 4, we have succeeded in confining the singular phase space integrations to just two variables in the n + 1 case, and four variables in the n + 2 case. For a single unresolved parton, the phase space measure contains (see Eq. The collinear limit between the unresolved and the reference parton is at η = 0. The relevant approximation of the tree-level matrix element in this limit in terms of a factorization formula is given in Eq. (D.2). This approximation also covers the soft-collinear limit, where η = ξ = 0. The pure soft limit of vanishing unresolved parton energy, ξ = 0, which is only singular for gluons, is described by Eq. (D. 35). The behavior of the matrix elements in the collinear and soft limits is at most as singular as 1 η Tree-level matrix elements with a single unresolved parton are to be found inσ R andσ C1 . We will discuss the construction of the subtraction scheme in the case ofσ R . The other cases follow exactly the same pattern. We havê where we sum over the unresolved, i, and reference, k, partons. f i,k (η, ξ) is independently regular at both η = 0 and ξ = 0. Unless the unresolved parton is a gluon, it even vanishes at ξ = 0. For example, in the case of an initial state reference momentum, using Eqs. (29) and (33) with z = 1, we obtain The phase space of the remaining partons, dΦ n p 1 + p 2 − u , the selector function, S i,k , the measurement function, F n+1 , and the matrix element, M (0) n+1 |M (0) n+1 , depend on η and ξ through the momentum, u, of the unresolved parton. The Laurent expansion ofσ R can be derived using Eq. (92) and the master formula where x is either η or ξ. Since the limits commute, the formula should be used recursively for both variables. Furthermore We will call the contribution of the delta-function, δ(x), the integrated subtraction term, or the pole term in the variable x. For this term, we will say that a pole has been taken in variable x. The end-point subtraction in the plus-distribution will be called the subtraction term. While the application of Eq. (94) indeed leads to an explicit expansion in with numerically integrable coefficients, the non-trivial part is the evaluation of the limits of the matrix elements multiplied with η ξ 2 . These are needed to evaluate f i,k (η, 0), f i,k (0, ξ) and f i,k (0, 0), and take the form As we pointed out at the beginning of this section, these limits can be obtained from factorization formulae, where the process dependent information is contained in matrix elements with n partons in the final state, while the process independent information is to be found in the splitting and soft functions. Therefore, in order to implement the scheme in general, one only needs to determine the limits of the splitting and soft functions. By construction, the resulting integrals will be pointwise convergent, since the amplitude limits were pointwise. These features will be present in the subsequent more complicated cases to be found in this section. The last contribution with n + 1 final-state partons isσ RV . The factorization formula for one-loop matrix elements present inσ RV in the collinear limit can be found in (E.2), whereas the one for the soft limit in (E.13). We notice that the scaling of the matrix elements in these limits is not uniform. In the collinear limit, there are two different terms proportional to 1 η and respectively. The second scaling is an artifact of the virtual integration. In the soft limit, on the other hand, we encounter 1 ξ 2 and The behavior in the soft-collinear limit is a combination of the two, since soft and collinear limits commute. The possible scalings are 1 η While we cannot simply use Eq. (94), we only need a minor modification. Suppose, that there is a function f (x), which behaves as then, we write This is a generalization of Eq. (94) for functions with non-uniform scaling. The formula should be used for η and ξ independently. The order of the variables, in which subtraction and integrated subtraction terms are introduced is irrelevant. We define the terminology of subtraction and integrated subtraction terms by analogy to the previous case. The double-real radiation contribution,σ RR , contains two unresolved partons. Eqs. (44), (63) and (81) show that the phase space always contains (102) The physical limits corresponding to the vanishing of the four variables, η 1 , η 2 , ξ 1 and ξ 2 , in the triple-collinear parameterization can be determined using Eq. (42) and Tab. 1. In the double-collinear parameterization, it is necessary to use Eq. (64) and Eq. (82). The required factorization formulae for the limits of tree-level matrix elements can be found in Appendix D. Multiplying the matrix elements by appropriate powers of η 1 , η 2 , ξ 1 and ξ 2 , we can rewrite any contribution toσ RR in the form The construction of the subtraction scheme amounts to the recursive use of Eq. (94). Notice that some of the nextto-next-to-leading order limits are iterated next-to-leading order limits. For example, if two different pairs of partons become collinear, we would apply formula Eq. (D.2) to each of the pairs independently. Similarly, if two partons become collinear, and another parton becomes soft, we would apply Eq. (D.2) to the collinear pair, and Eq. (D.35) to the soft parton. The last limit, which can be obtained in this way, is when two partons become collinear, and both of them become soft. This is a single-collinear double-soft limit. We would first apply Eq. (D.2) to the collinear pair, and then we would take the soft limit. The approximation to the matrix element is where {a 1 , a 2 } is either {g, g} or {q,q}, and is the soft current. The algorithm described in this section requires a multitude of matrix elements. The complete list including contributions, which do not necessitate subtraction is n .
(106) The one-and two-loop matrix elements can be further decomposed into divergent parts and finite remainders as shown in Appendix C. This, however, does not make the list any longer, and amounts to a replacement of |M (1,2) n,n+1 by |F (1,2) n,n+1 . In the case of gluons, spin correlated matrix elements are needed. Those are indicated by the presence of which just means that the polarization of gluon i is fixed to be λ i in the matrix element and λ i in the conjugated matrix elements. The remaining polarizations are summed over. In case there is the matrix element has a double spin correlation in gluons i and j. The spin correlators are due to the contraction of matrix elements with transverse vectors, which occur in the splitting functions, e.g. in Eqs. (D.4) and (D.5). Indeed, we can decompose any transverse vector as as long as k ⊥ · k = k ⊥ ·k = 0, wherek µ = k µ .

Average over azimuthal angles
The construction of the subtraction scheme presented in the previous sections relied, among others, on spin correlated splitting functions. The latter were necessary to guarantee the pointwise convergence of the numerical integration of the amplitudes and their respective subtraction terms. Nevertheless, explicit divergences of virtual amplitudes do not exhibit spin correlations. This suggests that it should be possible to obtain the integrated subtraction terms from azimuthally averaged splitting functions. There are two loopholes in this argument. First, it might happen that the integrated subtraction terms do not involve spin correlations in the coefficients of the poles in , while still containing them in the finite parts. Second, the splitting functions might not be contracted with matrix elements directly. In this section, we will demonstrate that both of these cases indeed take place. We will also point out when azimuthally averaged splitting functions may be safely used. Let us first consider the single-collinear case described by the parameterization of Section 4.1. In the limit θ = 0, besides the transverse vector u ⊥ , nothing depends on the azimuthal angles φ, ρ 1 , . . . . The spin correlator can be evaluated explicitly with where u ⊥ is given in Eq. (35), andr µ = r µ . This result follows from the vanishing of the µ, ν = 0 components, rotation invariance in the 1 − 2 dimensions parameterized by φ, ρ 1 , . . . , orthogonality to r, and normalization of the trace to unity. Transversality of the amplitudes implies that we can simply replace the spin correlated splitting function by its averaged counterpart. Similar arguments apply in the double-collinear sector, Section 4.3, for collinear limits of both unresolved momenta independently, and in the triple-collinear sector, Section 4.2, for the triple-collinear limit, and for the collinear limit of the unresolved momenta with respect to the reference momentum. The situation is slightly complicated by the fact that the parameterization of the second unresolved momentum depends on the first. In case the latter is collinear to the reference momentum, we have to use rotation invariance to decouple the momenta, and only then perform the azimuthal average. We will now be interested in the possibility to perform an azimuthal average of the splitting function generating the collinear pole related to the two unresolved momenta becoming collinear to each other, but not to the reference momentum. This case occurs in sectors S 4 and S 5 of the triple-collinear parameterization. We need to consider the integration over φ 2 in the triple-collinear sector parameterization with a general parameter ζ, Eq. (49). We thus study the integration measure at θ 2 ≈ θ 1 0. The relevant transverse vector u 3⊥ , Eq. (52), defines the azimuthal angleφ 2 , Eq. (53), The subsequent angular parameters are integrated over with the correct d-dimensional measure. Thus, if the integration measure overφ 2 , which is obtained by a non-linear remapping is correct, we can replace the spin correlated splitting function by its averaged counterpart. Let us proceed by expanding the relevant expressions around We first examine the contribution from the region θ 2 > θ 1 . There is The singularity in the limit is generated by the invariant s 12 which leads to This expression contains a regulated logarithmic singularity in the integration over θ 2 at θ 2 = θ 1 , which will result in a single pole in . Once the second region with θ 1 > θ 2 is added, the pole contribution is proportional to the integral Unfortunately, the azimuthal integration measure overφ 2 contains | tanφ 2 | −2 instead of sin −2 φ 2 . In consequence, replacing the spin correlated splitting function by the averaged one will only be correct at = 0. On the other hand, the integrals can still be performed exactly where u 3⊥ , given in Eq. (52), does not depend on σ 1 , σ 2 , . . . atφ 2 = 0. Moreover,ū µ 1 = u 1 µ , with u µ 1 given in Eq. (37). The result follows from the vanishing of the µ, ν = 0 components, orthogonality to u 1 , rotation invariance in the space parameterized by σ 1 , σ 2 , . . . , normalization of the trace to unity, and finally from the explicit value of the integral when contracted with u µ 3⊥ (φ 2 = 0) u ν 3⊥ (φ 2 = 0). In practice, we can thus replace the spin correlated splitting functions Eqs. (D.4) and (D.5) as followŝ The terms in the square brackets should be compared with the averaged splitting functions Eqs. (D.9) and (D.10). As expected, the functions are different at order . This discussion is relevant first and foremost to the single-collinear pole generated by η 2 in sector 4, and by η 1 in sector 5 of the triple-collinear sector parameterization, as can be verified using Tab. 1. Clearly, there will be no spin correlations in the coefficient of the pole itself, but the matrix element of the finite part will be contracted with u 3⊥ (φ 2 = 0). Note that the subtraction terms to this contribution in the collinear limit of u 1 with respect to the reference momentum r should be derived by using an iteration of splitting functions, instead of the triple-collinear splitting function. The relevant approximation of the azimuthally averaged matrix elements takes the general form where the variables with subscript "12" describe the first limit u 1 ||u 2 , whereas the variables with subscript "r12", the second limit r||u 1 + u 2 . If spin correlations are present in the first limit, then The second contribution, which is affected by the non-trivial averages of the spin correlated splitting functions, is the single-collinear double-soft double pole, generated by the pair of variables {ξ 2 , η 2 } in sector 4, and {ξ 2 , η 1 } in sector 5 of the triple-collinear sector parameterization (see Tab. 1). In both cases, as long as ξ 1 is non-vanishing, the collinear limit of the two unresolved partons will generate spin correlations. The double-pole contribution can be obtained from the matrix element factorization (see Eq. (104)) The splitting functions should be replaced according to Eq. (118). In this case, the difference between the correct replacement and the averaged splitting functions affects the single pole. Thus, using averaged splitting functions would lead to an incomplete cancellation of the divergences in the cross section. There are no other instances, where Eq. (118) must be used, because all other subtraction and integrated subtraction terms contributing to the cases just discussed do not involve spin correlations. Table 3: Single-unresolved contributions to double-real radiation,σ RR SU . Listed are the variables leading to pole (integrated subtraction) terms after application of the algorithm of Section 5. Each pole is accompanied by all possible subtraction terms. Sectors S 1 , . . . , S 5 belong to the triple-collinear sector parameterization of Section 4.2. The variables of the double-collinear sector parameterization are defined in Section 4.3.

Separation of finite contributions
In order to define the subtraction scheme in 't Hooft-Veltman regularization, which will be done in the next section, it is necessary to understand which cross section contributions are separately finite. Clearly, the leading order cross section is finite and can be directly evaluated in four dimensions. At next-to-leading order, the situation is slightly more complicated. From the three contributions,σ R ,σ V ,σ C , we can decompose the first two as followŝ The subscript "F" stands for "finite" here and below, while "U" stands for "unresolved". For real radiation we havê with the notation of Section 2. The subtraction terms are generated using the algorithm of Section 5. Any integrated subtraction terms generated along the way belong toσ R U . Thus,σ R U only involves tree-level matrix elements with n final-state particles. For virtual corrections we havê where the finite remainder |F (1) n and the singular color-space operator Z (1) are defined in Appendix C. The finiteness of the next-to-leading order cross section implies that the following contributions are separately finitê At next-to-next-to-leading order, we begin by decomposing the double-real radiation as followŝ The single-unresolved (SU) contribution,σ RR SU , contains all integrated subtraction terms with n + 1 resolved particles including the necessary subtraction terms. The relevant pole terms are listed in Tab. 3. The double-unresolved (DU) contribution,σ RR DU , contains all integrated subtraction terms with n resolved particles including the necessary subtraction terms. The relevant pole terms are those not listed in Tab. 3. The real-virtual cross section can be decomposed as followsσ RV The limits required to generate the subtraction terms of 2Re The other two contributions,σ RV FR andσ RV DU , contain all the integrated subtraction terms. The latter are distributed such thatσ RV FR involves finite remainders (FR) of one-loop amplitudes only, whereasσ RV DU the left-over Born matrix elements. Both types of matrix elements correspond to amplitudes with n final-state particles. In order to obtain the relevant expressions, it is necessary to use the formulae of Appendix E.1 and Appendix E.2. The double-virtual cross section is decomposed as followŝ Finally, the factorization contributions are decomposed as followŝ The separately finite contributions arê The finiteness ofσ FR can be proven by noticing that the sumσ FR +σ SU +σ DU is finite by the finiteness of the next-to-next-to-leading order cross section, and the analytic structure of the matrix elements inσ FR andσ SU +σ DU is different. Indeed,σ FR may have thresholds due to virtual integrations, which cannot be present inσ SU +σ DU , as the latter only involves tree-level matrix elements. Thus,σ FR must be separately finite. Another way to conduct the proof is to notice thatσ FR is generated from leading order splitting and soft functions (compare with Appendix E.1 and Appendix E.2). It is thus finite by the finiteness of the next-to-leading order cross section, since the unresolved contributions to the latter are derived from exactly the same splitting and soft functions [19]. The finiteness ofσ FR implies, of course, the finiteness ofσ SU +σ DU . Using a suitable measurement function, it is possible to obtain the next-to-leading order cross section for n + 1 well-separated partons from a next-to-next-to-leading order cross section calculation for n well-separated partons. This measurement function would satistify F n = 0. It would set the contributionsσ VV F ,σ FR , andσ DU to zero. On the other hand,σ RR F ,σ RV F andσ SU would correspond toσ R F ,σ V F andσ U as followŝ At this point, there are two remaining contributions, the single-and double-unresolved, which are not finite by themselves. Both involve tree-level matrix elements only. Nevertheless, they have a different physical interpretation. The single-unresolved contribution corresponds to an inclusive phase space integral over n + 1-resolved-particles kinematics. The phase space singularities are regulated with suitable subtraction terms. The double-unresolved contribution, on the other hand, corresponds to a phase space integral over n-resolved-particles kinematics. Due to Table 4: Separately finite parts of a next-to-next-to-leading order cross section divided into: finite (F), unresolved (U), finite-remainder (FR), single-unresolved (SU), and double-unresolved (DU) contributions. Precise definitions can be found in the text.
appropriate separation cuts inherent in the measurement function, there are no further phase space singularities. For the purpose of a four-dimensional formulation of the subtraction scheme, we must render both contributions separately finite. Since their sum is finite, we can concentrate on one of them. Once it is rendered finite by specifically designed counterterms, it is sufficient to subtract the same counterterm from the other. Clearly, the single-unresolved contribution is substantially less complex than the double-unresolved. We will thus proceed with its analysis. Once our task is finished, we will have the separately finite contributions summarized in Tab. 4. The single-unresolved contribution consists of three parts: double-real radiation,σ RR SU , real-virtual radiation,σ RV SU , and factorization,σ C1 SU . In a first step, we notice that if, instead of a next-to-next-to-leading order measurement function, we would use a next-to-leading order one, which would require n + 1 well separated particles, the doubleunresolved contribution would vanish. Thus, the single-unresolved contribution would be finite. Once we return to the original measurement function, divergences do not cancel in the sumσ RR SU +σ RV SU +σ C1 SU anymore. This lack of cancellation is, therefore, due to the subtraction terms. Let us now concentrate on the real-virtual and doublereal corrections causing final state divergences. Neglecting irrelevant integration variables, the essential part of any contribution toσ RV SU in any phase space sector has the following form in the parameterization of Section 4.1 (see Section 5) The function f (η, ξ) is a product of a selector function, the matrix element of the Z (1) operator, and the measurement function. The terms with η and/or ξ vanishing are the subtraction terms. The essential point is that they are integrated over the full range of variation of the resolved particle kinematics given by η and ξ, i.e. over the unit interval. Contributions toσ RR SU can be classified according to the number of poles taken. Double-pole contributions always correspond to the disappearance of the second unresolved parton with momentum u 2 , as can be checked by direct inspection of the parameterizations of Section 4. They take the form where g(η 1 , ξ 1 ) is a product of a selector function, splitting or soft function, matrix element, and the measurement function. η 1 and ξ 1 can be identified directly with η and ξ from (140) due to the way they enter the matrix element. g(η 1 , ξ 1 ) gives a contribution to the divergences of − f (η, ξ). In the same way, any of the subtraction terms corresponding to vanishing η 1 and/or ξ 1 constitutes a contribution to the subtraction terms of the divergences of − f . If we could write all contributions toσ RR SU in the same way, not only the divergences of the integrals of f (η, ξ) and g(η 1 , ξ 1 ) would cancel, but the same would also be true of their subtraction terms. Unfortunately, the situation is complicated by the left-over integration present in the single-pole contributions to double-real radiation. We have to deal with integrals of the form 1 0 dy y 1+a dx 1 where {y, x 1 , x 2 } ⊂ {η 1 , η 2 , ξ 1 , ξ 2 }, and the pole has been taken in the variable {η 1 , η 2 , ξ 1 , ξ 2 } \ {y, x 1 , x 2 }. It can be checked, that in all contributions, there is always one variable, which can be directly identified with either η or ξ from (140). We denote this variable by y. The second kinematic variable, denoted by x (e.g. if y corresponds to η, then x corresponds to ξ), is a function of x 1 and Using the same arguments as in the double-pole case, one can convince one-self that there is nothing to correct, if x = x 1 or x = x 2 . However, if the transformation between these variables is non-linear, the divergences of the subtraction terms built of g may not correspond to the divergences of the subtraction terms built of f . Alternatively, it may turn out that they are integrated over a different range, i.e. not over the unit interval. In each of these cases, we must introduce counterterms to compensate for the difference. It is important that when x → 0, the divergences of the subtraction terms must constitute correct contributions to the divergences of the subtraction terms of − f . The reason is that the unsubtracted integral gives a correct contribution to the divergences of the integral of − f (η, ξ), and the subtraction terms are defined as its limits at vanishing y and/or x. This allows to derive the counterterms by inspection of the limit x → 0. We now discuss different cases separately. Initial state divergences will be treated last.

Case I
We assume that x is related to x 1 , x 2 by where c is some function. Notice that it is then the difference in (142) that gives a contribution to the unsubtracted divergences, since neither of the two terms vanishes in the presence of a next-to-leading order measurement function. Nevertheless, each of them has a different relation between x 1 , x 2 and x. Indeed, for g(y, x 1 , 0) there is x = x 1 . In consequence, subtraction terms with x 2 = 0 will not require counterterms. Let us then consider terms with a non-trivial dependence on x 2 and define h(x 2 ) = − g(y, 0, x 2 ) − g(0, 0, x 2 ) .
Omitting the irrelevant integration over y, the subtraction terms take the form 1 0 dx 1 where If the second integral has any dependence on x, it must be corrected. Such a dependence may be induced by a nontrivial dependence on x of x 2 max (x) or of the integrand. As discussed above, the correct behavior is defined by the limit The counterterm to be added to the single-unresolved contributions is, therefore The integral over x can be performed explicitly by a change of the order of integration, with the result 1 0 dx 2 Thus finally, the counterterm takes the form The same counterterm will have to be subtracted from the double-unresolved contributions in order to render them finite. Below, we give specific values for three different cases. They depend on the scaling of x 2 , i.e. on the exponent b 2 . We provide the latter for CDR and HV regularizations. The latter is given in anticipation of the discussion of the next section. Notice that double-unresolved contributions will always use the CDR values, while it will only be the single-unresolved contributions that will be affected by the regularization scheme change.
7.1.1. Triple-collinear parameterization sector S 4 , collinear pole in η 2 The variable assignments are y = η 1 , In this case, the subtraction term (147) takes the form We notice that the integrand of the x 2 integration depends on x in the second term, while the integration range depends on x in both terms. This leads to the necessity of compensation. The counterterm (150) is Returning to the original variables, the counterterm is parameterized by with the scaling b CDR 7.1.2. Triple-collinear parameterization sector S 5 , collinear pole in η 1 This case is identical to the previous upon the interchange η 1 ↔ η 2 . Thus 7.1.3. Triple-collinear parameterization sector S 5 , soft pole in ξ 2 The variable assignments are y = ξ 1 , The subtraction term (147) takes the form We notice that, while the integrand of the second integral is independent of x, the integration range is not. In order to compensate for this fact, we introduce the counterterm Returning to the original variables, we obtain with the scaling b CDR 7.2. Case II: triple-collinear parameterization sector S 2 , soft pole in ξ 2 The variable assignments are y = ξ 1 , The reason for separate treatment is that the resolved parton momentum is parameterized in a symmetric way by the sector variables It is, therefore, only g(y, x 1 , x 2 ) that contributes to the unsubtracted divergences. The symmetry requires to consider counterterms with a dependence on either variable, x 1 and x 2 . We treat explicitly the x 2 case with a generic function h(x 2 ). The relevant contribution to the subtraction term can be rewritten as We notice that the range of the resolved parameter x is restricted to [0, 1/2] although it should be [0, 1]. At the same time the integration over x 2 has an x dependent range. We first determine the correct behavior at x → 0. This can be achieved by writing the right-hand side of Eq. (168) as The first two terms in the square bracket have a uniform scaling in x as required of subtraction contributions matching the divergences of the real-virtual single-unresolved cross section. The last term must be removed by a counterterm. It is also necessary to extend the integration range of x. In consequence, the complete counterterm has the form The first term in the square bracket can be integrated over both variables with the result The second and third terms in (170) can be combined together with the result Using the above result and symmetry of the expressions with respect to the interchange x 1 ↔ x 2 , we obtain the complete counterterm for sector S 2 expressed through the original variables with the scaling b CDR 7.3. Case III: double-collinear parameterization and triple-collinear parameterization sector S 3 , collinear pole in η 1 The variable assignments are y = η 2 , x 1 = ξ 2 , x 2 = ξ 1 .
This case is different from those treated until now, because at η 1 = 0 the unresolved parton with momentum u 1 is collinear to the reference parton, which is assumed to be in the final state. In the notation of Section 4, the reference momentum is r 1 in the double-collinear parameterization, and r in the triple-collinear parameterization. Here, we will denote it by r irrespective of the case. As far as the double-collinear parameterization is concerned, we will treat the general case of subsection 4.3.1 throughout, and only return to the special case of subsection 4.3.2 at the end. Since both partons, the unresolved and reference, are moving in the same direction, the reference momentum for comparison with the real-virtual corrections is not r, but rather Let us introduce rescaled variables for the energy of the reference parton: ξ r for the original one, and ξ r for the composite one. Taking into account the different variation ranges, we write where ξ r max denotes the maximum of r 0 /E max , and similarly for ξ r max . By definition, ξ r , ξ r ∈ [0, 1]. The variation range of r 0 is the same as in the respective real-virtual contributions. Since the comparison betweenσ RV SU andσ RR SU involves the reference momentum, we define a new function h through g(y, x 1 , x 2 ) = 1 0 dξ r h(ξ r , y, x 1 , x 2 ) .
Therefore, we only need to consider the subtraction terms related to the functions The result for the correction to the second one must be the limit at y = 0 of the result for the correction to the first one. Thus, we concentrate on h(ξ r , y, 0, x 2 ). Due to the special kinematics We can now specify the relationship between the variables of the double-real contribution, and those of the real-virtual contribution where The variablez is integrated over in order to obtain a contribution, which cancels the poles of the matrix element of Z (1) present inσ RV SU for fixed kinematics specified by x, y and ξ r . The soft limit in the integration over the phase space of the unresolved parton with momentum u 1 corresponds toz = 0 (this is the reason for the bar in the notation, since usually the soft limit is at z = 1). The function x 2 max depends on y through the scalar productû 1 ·û 2 . In the triple-collinear case, this dependence is simple, sinceû 1 ·û 2 = 2y. On the other hand, the dependence in the double-collinear case is only indirect and involves the angles of the other reference momentum. For convenience, we define It turns out that The integral of the subtraction term takes the form 1 0 dξ r dy y 1+a dx 1 where we have introduced a shorthand notation,h, for the integrand. To derive the behavior of the subtraction term at small x, we must work with the appropriate order of the integration variables. For now, we shall neglect y, which is irrelevant to this problem. In view of the parameterization of the contributions toσ RV SU , the integration order must be dx dξ r dz .
Using Eq. (186), we obtain The correct subtraction term matching the poles of the respective subtraction term of the real-virtual contribution, is obtained from this expression by taking the limit x → 0 of the square bracket, and extending the integration range over x. It reads where we have made use of the fact that a potential singularity at ξ r = 0 is regulated by the selector function. The counterterm we are seeking is the difference between Eq. (190) and Eq. (189). The contribution containingh(ξ r , y,z) is easily obtained by returning to the original order of integration variables, and noticing that it is only necessary to extend the integration range over This expression is not integrable at = 0, but we can extract the explicit pole in by a subtraction 1 0 dξ r dx 2 1 0 dξ r h(ξ r , y, 0, 0) The second contribution to the counterterm is obtained from Eq. (190), by taking only the terms involvingh(ξ r , y, 0, 0) As expected, the divergence cancels in the sum of Eqs. (192) and (193). The counterterm is thus finite and integrable. Returning to the original variables, the counterterm to be added to the single-unresolved double-real radiation contribution, and subtracted from the double-unresolved contribution is Finally, we note that the counterterm for the special case of the double-collinear parameterization of subsection 4.3.2 is the same as above with the replacement

Case IV: initial state divergences
The corrections needed in the case of initial state divergences have their origin in the same pole as Case III. We will thus consider the collinear pole in η 1 in the double-collinear parameterization, and in sector S 3 of the triple-collinear parameterization. However, the reference momentum will now be in the initial state. Furthermore, we must compare the contributions toσ RR SU with those toσ C1 SU , and not those toσ RV SU . Instead of Eq. (140), we have where z is the variable used in the convolution with the splitting functions in Eq. (8), and the latter have been included in f . In order to facilitate the comparison, we assume that the kinematics in the calculation ofσ C1 SU is the same as in σ RR SU , which means that z is defined by the emission of a fictitious parton from the initial state, such that the parton momentum entering the matrix element is zr, r being the same reference momentum as for the unresolved parton with momentum u 1 inσ RR SU (one of the initial state momenta by assumption). Of course,σ C1 is boost invariant, and we can evaluate it in any system we like. The question whether its partial contributionσ C1 SU also has this property will be discussed at the end of this subsection. The relevant contribution toσ RR SU has the form of Eq. (142). We rewrite it through the original variables of the contribution . (198) ξ 1 and ξ 2 are related to z and ξ through Just as in Case III (compare to Eq. (179)), we note that We thus only need to consider g(η 2 , ξ 1 , 0) , g(0, ξ 1 , 0) .
We will work with g(η 2 , ξ 1 , 0), as the correction to g(0, ξ 1 , 0) can be derived as the limit to that of g(η 2 , ξ 1 , 0). Neglecting the irrelevant integration over η 2 , the integral of the subtraction term is The necessary counterterm amounts to the extension of the integration range over ξ. It is thus As it stands, this correction is not integrable in ξ 1 if g(η 2 , 0, 0) 0. The reason is that ξ 2 max /ξ 2 max = ξ 1 + O(ξ 2 1 ). The singularity at ξ 1 = 0 corresponds to z = 1. There is a subtraction at z = 1 in the convolution with the splitting functions in the factorization contribution,σ C1 SU . In the double-real contributions we are considering, an analoguous subtraction is not necessary, since the integration range of ξ vanishes at z = 1, as seen in Eq. (202). The counterterm we have derived reintroduces the singularity at the endpoint. This singularity must be subtracted in order to match the factorization contributions. This amounts to replacing We thus finally arrive at the complete counterterm We have derived the counterterm for this case by comparing the single-unresolved double-real radiation contribution to the single-unresolved factorization contribution in a boosted frame. Unfortunately, the separation into single-and double-unresolved contributions is not Lorentz invariant, because it is based on taking poles in energy and angle variables. In the next section, we will manipulate differently both types of contributions. In consequence, the construction will only be correct, if the factorization contributions,σ C1 SU andσ C1 DU , will be evaluated in the boosted frame as assumed here. This is the reason for the necessity of the unusual parameterization of Section 4.1.

't Hooft-Veltman regularization of separately finite contributions
The final stage of our construction is the introduction of the 't Hooft-Veltman regularization. The latter differs from the conventional dimensional regularization in the description of the resolved partons. In HV, their momenta and polarizations are four-dimensional, while in CDR they are d-dimensional. The differences in the number of spin degrees-of-freedom only affect the gluons. In HV, tree-level matrix elements do not have any expansion in . Indeed, the dependence is due to spin sums over squared matrix elements. If these sums are only restricted to fourdimensional degrees-of-freedom, there can be no dependence on . Due to the virtual integrations, one-and two-loop matrix elements do have a non-trivial dependence on . Nevertheless, our goal will be to only work with the first term of the -expansion of the finite remainders. These can also be viewed as four-dimensional loop corrections. We note that the two features of HV regularization, four-dimensional momenta and four-dimensional polarizations, are quite different. Since our subtraction scheme makes extensive use of phase space parameterizations, the restriction to four-dimensional momenta will be achieved by modifying the phase spaces. On the other hand, the restriction to fourdimensional polarizations amounts to the removal of higher order terms of the -expansion of the matrix elements. We now consider the different contributions described in Section 7 in increasing level of complexity.
The simplest contributions do not involve any singularities in . These arê We can simply set = 0 in the phase space integrals and in the matrix elements. We thus directly obtain HV regularized contributions, as if the calculation were done completely in four-dimensions without ever making reference to dimensional regularization. The second group of contributions only has n resolved partons. These arê We first note that since these contributions are separately finite, we can drop the higher order terms in the -expansion of the matrix elements. Indeed, the cancellation of divergences is due to the unresolved phase space integrals of the soft and splitting functions, and the form of the divergences of the virtual amplitudes, which are contained in the Z (1,2) operators. The cancellation is not related to any particular functional dependence of the amplitudes on the kinematics. Thus, it occurs separately for any order of the expansion in of the matrix elements. However, while the leading order will give a finite contribution, the -suppressed subleading orders will give vanishing contributions in the limit → 0. In consequence, we can remove them from the very beginning. We stress that this mechanism only works after azimuthal averages have been performed. Otherwise, spin correlators would mix different orders of the -expansion, and at least the finite parts of the results would not be correct. We point out that the cancellation of divergences happens even at the level of independent color correlated amplitudes. We do not see, however, much use of this fact for our purposes, besides testing of course.
At this point, we have four-dimensional polarizations. We now note that the contributions under consideration are finite for any measurement function, as long as it is infrared safe. F n is the only measurement function present, and we can use it to restrict the resolved momenta to be four-dimensional with the replacement where q i are the momenta of the final-state resolved partons, and the delta-functions remove their -dimensional components. We note that it is only possible to explicitly restrict the momenta of all but one of the final-state particles. This is due to momentum conservation. On the other hand, the same momentum conservation will always provide a restriction on the last final-state momentum. Thus, the result is independent of the subset of momenta chosen in (209). In the simplest case, where there are either no reference momenta or they are in the initial state, the replacement amounts to the following change of the resolved parton phase spaces In other words, the calculation can be performed with a four-dimensional phase space for the resolved partons. The situation is more complicated, if reference momenta are in the final state. The rule (210) is still valid for all the final-state particles not explicitely parameterized in the given sector. There will appear, however, delta-functions of the form δ (−2 ) (r + u) = (r 0 + u 0 ) 2 δ (−2 ) (r) for the single-collinear sector δ (−2 ) (r + u 1 + u 2 ) = (r 0 + u 0 1 + u 0 2 ) 2 δ (−2 ) (r) for the triple-collinear sector, δ (−2 ) (r 1 + u 1 ) δ (−2 ) (r 2 + u 2 ) = (r 0 1 + u 0 1 )(r 0 2 + u 0 2 ) 2 δ (−2 ) (r 1 ) δ (−2 ) (r 2 ) for the double-collinear sector. (211) These delta-functions will restrict the angular integrations for the reference momenta to be four-dimensional. The energy factors are, nevertheless, non-trivial. In the case of collinear poles, the integration over the reference momentum energy is not simply four-dimensional, since the actual resolved momentum is a sum of the reference and unresolved momenta. It remains to argue that the restriction (209) does not affect the finite part of the result. This is a true assertion, since the result is finite independently of the measurement function. The missing -dimensional integrations are finite, but evaluated over a volume of the phase space, which vanishes at = 0. Thus the missing contribution vanishes as well and can be neglected from the start. We have thus shown that the HV regularization for this group of contributions amounts to evaluating matrix elements and resolved-parton phase spaces in four dimensions. The integration over the unresolved partons is nevertheless performed in d dimensions. In particular, if there is one unresolved parton and the contribution is due to a soft limit, the parton kinematics will be integrated non-trivially over the fifth dimension. Similarly, if there are two unresolved partons and the contribution is due to a double-soft limit, one parton will be effectively integrated in five dimensions, while the other in six. This is, however, the upper bound on the number of required dimensions independently of the multiplicity.
The last contribution to consider is the single-unresolved cross section It differs from the previous case in one aspect. It involves both n + 1-and n-parton matrix elements. Let us first assume that the calculation is performed with a next-to-leading order measurement function, i.e. F n = 0. In this case, we can just repeat the previous discussion with n → n + 1, and simply evaluate n + 1-parton matrix elements and n + 1-parton phase spaces in four-dimensions. However, once we turn to the next-to-next-to-leading order measurement function, i.e. F n 0, the subtraction terms will not match the singularities of the amplitudes anymore. The reason is that the limits originally used in the construction of the subtraction terms in Section 5 apply to d-dimensional amplitudes. Forσ RV SU andσ C1 SU , this problem is easily resolved. In the subtraction terms, we take the splitting functions at = 0. The soft functions do not depend on and do not require modifications. The calculation ofσ RV SU andσ C1 SU can now be performed with a four-dimensional phase space and four-dimensional matrix elements.
The double-real single-unresolved contribution,σ RR SU , requires more care. The pole contributions must be evaluated with complete splitting functions, i.e. at 0. However, when subtraction terms to these pole terms are derived, we take into account that the limits are iterated. This means, for example, that instead of the triple-collinear splitting function, Eq. (D.20), we can use a product of two splitting functions. The first one will generate the pole, and it must be taken azimuthally averaged and at 0. The second one corresponding to the subtraction term must be taken spin-correlated and at = 0. An example of this procedure is the pole at η 2 due to a quark and a gluon in the triple-collinear parameterization sector S 4 , to which we would generate a collinear subtraction term with a gluon as reference parton. The corresponding limit would be a special case of Eq.
Since next-to-leading order soft functions do not depend on , there is nothing special to do in their case. This concerns the cases of a soft pole and collinear subtraction, and of a collinear pole and soft subtraction. We note, however, that the iterated soft-pole-soft-subtraction limit is easier to obtain form the double-soft limit. The double-soft function does depend on . However, this dependence does not contribute to the iterated limit. After each contribution listed in Tab. 3 has been treated this way, we may evaluate all tree-level matrix elements, both with n + 1-and n-partons in the final state, in four dimensions. The derivation of the subtraction and integrated subtraction terms is performed using the algorithm of Section 5, but taking into account that the phase space restriction (209) now acts on the unresolved partons u 1 and u 2 . Since we use (209) with n → n + 1, the following delta-functions will appear in the triple-collinear parameterization δ (−2 ) (r + u 1 ) δ (−2 ) (u 2 ) for the collinear pole in η 1 in sector S 3 , δ (−2 ) (r + u 2 ) δ (−2 ) (u 1 ) for the collinear pole in η 2 in sector S 1 , δ (−2 ) (r) δ (−2 ) (u 1 + u 2 ) for the collinear pole in η 1 in sector S 5 , and in η 2 in sector S 4 , δ (−2 ) (r) δ (−2 ) (u 1 ) for the soft pole in ξ 2 in sectors S 1 , S 2 , S 4 and S 5 . (214) In the double-collinear parameterization, on the other hand, there will be δ (−2 ) (r 1 + u 1 ) δ (−2 ) (r 2 ) δ (−2 ) (u 2 ) for the collinear pole in η 1 , δ (−2 ) (r 2 + u 2 ) δ (−2 ) (r 1 ) δ (−2 ) (u 1 ) for the collinear pole in η 2 , Delta-functions involving a reference momentum are only present, if it is in the final state. Notice that in all collinear pole cases, there are no d-dimensional integrations left. However, there is a d-dimensional integration over the unrestricted direction of u 2 in the soft-pole case. In practice, most integrations can be performed analytically, because nothing depends on the respective angles. There will, however, remain a single integration beyond the dimensions of the resolved momenta. Thus, in the general case, there will be a five-dimensional integration. The delta-functions, which do not involve the reference momentum, influence the scaling in the singular variables. This, in turn, has consequences for the counterterms derived in Section 7. There, we gave the necessary scaling exponents for both CDR and HV cases, the latter following from the application of the above listed delta-functions. Whenever the exponents concerned energy variables (as for collinear poles of subsections 7.1.1, 7.1.2, 7.3 and 7.4), ξ 1,2 , the difference between the exponent in CDR and that in HV was 2, simply because the delta-functions provide an additional factor of ξ 2 1,2 . On the other hand, the scaling of the angular variables in the soft-pole cases of subsections 7.1.3 and 7.2 was affected by the change of the angular measure. Indeed, the delta-function δ (−2 ) (u 1 ) introduces a factor of which changes the exponent of the relevant angular variables by 1 (η 2 in sector S 5 , and both η 1 and η 2 is sector S 2 ). Notice that the difference in scaling only influences the finite parts of cross sections. In consequence, if one would proceed along the algorithm of this section, but without separating the single-and double-unresolved contributions, and without applying the counterterms of Section 7, the result would be finite, but incorrect. This remark completes the discussion of the 't Hooft-Veltman regularization of the single-unresolved contributions, and by the same of the complete cross section.
Finally, we remind that the factorization contributionsσ C1 SU andσ C1 DU must be calculated in a boosted frame as discussed in subsection 7.4. One could wonder, if this means that the complete calculation is frame dependent. This is not the case, since after all the modifications, the result for the O( 0 ) cross section is exactly the same as in CDR. In CDR, however, each cross section contribution listed in Section 2 is separately Lorentz invariant. 9. Example: g g → tt + ng, n = 0, 1, 2 In this section, we present a comparison between results obtained using conventional dimensional regularization and 't Hooft-Veltman regularization for an example cross section at next-to-next-to-leading order. The aim is to demonstrate that the modifications described in Sections 7 and 8 lead to correct results in a realistic calculation. We select inclusive top-quark pair production in the gluon fusion channel with up to two gluons in the final state. This choice is motivated by the fact that collinear limits involving gluons have the most involved structure of spin correlations. We are thus able to verify the azimuthal average corrections we have provided in Section 6. We neglect all contributions involving finite remainders of one-and two-loop amplitudes. As explained before, their independence of the regularization scheme can be proven explicitly by noticing that they enter the calculation in exactly the same manner as Born amplitudes in a next-to-leading order calculation. Finally, we do not include the contribution of the subtracted six-point Born amplitude, which is trivially independent of the regularization, as it is finite by definition.
The relevant partonic cross section is rendered dimensionless and independent of the value of the strong coupling with the normalizationσ ( whereσ (2) is the total cross section contribution at O(α s 4 ). Furthermore, we set Our results are obtained at the point Tables 5 and 6 contain results for partial double-unresolved contributions calculated in CDR and HV regularizations respectively. In each case, the last row gives the total double-unresolved contribution. Whenever errors are quoted, they are due to Monte Carlo integration. Contributions involving two-parton kinematics have been computed with a deterministic integration method, which implies that their error is negligible, and, therefore, not specified in the tables. By construction, double-unresolved contributions should be finite in both regularizations. Indeed, coefficients of the poles in are consistent with 0 within one standard deviation for all but the leading singularity in CDR, where consistency at a level below two standard deviations is observed. The calculation has been performed with optimization based on the behavior of the finite part, which is one reason for the slightly lower quality of the leading pole contribution. Notice, nevertheless, that there is no difference between CDR and HV at the leading pole of each contribution as far as the actual integrand is concerned. Therefore, divergence cancellation in the HV case within one sigma is sufficient to claim divergence cancellation in both regularizations. On the other hand, analytic cancellation of the coefficient of the 1/ 4 pole has already been shown in [7]. Due to severe cancellations between the different partial contributions, it was necessary to use large Monte Carlo samples. For instance, the quoted precision of the doublereal contributions required nearly 10 11 points. In all cases, the convergence in the HV regularization was noticeably better. As far as the finite parts of the results are concerned, we aimed at about 1% precision for the total contribution. Nevertheless, the agreement between evaluations in CDR and HV regularizations is at the level of one permille.    Table 6: Double-unresolved (DU) contributions to the partonic cross section gg → tt + X, with X consisting of up to two gluons, evaluated in 't Hooft-Veltman regularization (HV). The error estimates quoted in parentheses are due to Monte Carlo integration. The definition of partial contributions is given in the text.
The partial results quoted in the tables have been obtained by integrating tree-level amplitudes only. They are: σ VV DU : Double-virtual contributions obtained by integrating the two-loop and one-loop squared amplitudes for gg → tt without their finite remainders.
σ RV DU : Real-virtual contributions obtained from the integrated subtraction terms of the one-loop amplitude for gg → tt + g, without the contribution of the finite remainder of the one-loop amplitude for gg → tt.
σ RR DU : Double-real contributions obtained from the double-unresolved integrated subtraction terms of the Born amplitude for gg → tt + gg, including corrections described in Section 8, which make the total double-unresolved contribution finite.
σ C1 DU : Factorization contributions obtained from the convolution of the leading order splitting function with the cross section contribution of the integrated subtraction terms of the Born amplitude for gg → tt + g. σ C2 DU : Factorization contributions obtained from the convolution of the leading order splitting function with the cross section contribution of the one-loop amplitude for gg → tt without its finite remainder, and the convolution of the next-to-leading order splitting function as well as two leading-order splitting functions with the Born cross section for gg → tt.
The results in CDR have been obtained without azimuthal averaging, i.e. with splitting functions containing full spin correlations. The results in HV, on the other hand, have been obtained with azimuthal averaging, i.e. with averaged splitting functions whenever possible. Details can be found in Section 6. The counterterms of Section 7 have been applied to averaged splitting functions. We finally note that in the HV regularization, the finite parts receive nonvanishing contributions from integrated subtraction terms only. Indeed, there is no contribution fromσ VV DU andσ C2 DU . The latter will, however, contribute if µ R µ F .    Results for the single-unresolved partial contributions are displayed in the Tabs. 7 and 8 for CDR and HV regularizations respectively. We again observe finiteness of the total contributions given in the last row of each table. The coefficients of the poles are consistent with 0 at the one sigma level in all cases but the leading singularity evaluated in HV regularization. There, consistency is only observed at a level better than two standard deviations. In this respect, the same comments apply as in the discussion of the double-unresolved contribution. The single-unresolved contributions in both regularizations have a precision of better than two permille relative error. The agreement between them is at the same level, and does only slightly exceed one standard deviation.
The partial cross section contributions in the single-unresolved case are: σ RR SU : Double-real contributions obtained from the single-unresolved integrated subtraction terms of the Born amplitude for gg → tt + gg, including corrections described in Section 8, which make the total single-unresolved contribution finite. The splitting functions used in the derivation of the integrated subtraction terms are given by the azimuthally averaged expression Eqs (D.9). The correct result is obtained after addingσ A SU . σ RV SU : Real-virtual contributions obtained by integrating the one-loop amplitude for gg → tt + g together with its subtraction terms, after removal of all finite remainders.
σ C1 SU : Factorization contributions obtained from the convolution of the leading order splitting function with the cross section contribution of the Born amplitude for gg → tt + g together with its subtraction terms.
σ A SU : Difference between the single-unresolved contributions obtained with spin-correlated and azimuthally-averaged splitting functions in integrated subtraction terms as explained in Section 6.
We notice thatσ RV SU andσ C1 SU only contain poles in HV regularization. However,σ C1 SU would develop a finite part if µ R µ F .
Finally, we can compare the single-and double-unresolved contributions to the full partonic NNLO cross section, which was calculated in Ref. [11]. At the chosen value of β, there is The single-unresolved part contributes 20% to the full cross section, whereas the double-unresolved part only contributes 2% in this specific case.

Concluding remarks
We have presented a complete construction of the sector-improved residue subtraction scheme in four dimensions.
It is now possible to evaluate next-to-next-to-leading order cross sections using ordinary tree-level matrix elements without higher order terms of the -expansion. This is crucial, since it allows to use the myriad of publicly available software designed for efficient tree-level calculations. Of course, it is also necessary to have access to one-and twoloop amplitudes. Fortunately, at least the former are also available from open access packages. The problem is thus currently reduced to the two-loop virtual corrections. Recent progress shows that a breakthrough may be possible on the scale of the next few years. At present, our construction is an algorithm, which leads directly from various soft and splitting functions collected in appendices of this paper, to process independent subtraction and integrated subtraction terms necessary for a Monte Carlo implementation of the phase space integration. In the future, we intend to provide an optimized code that will handle the burden of bookkeeping and will contain built-in tree-level Standard Model amplitudes, similarly to the next-to-leading software packages of Refs. [20] and [21]. Virtual amplitudes will still have to be provided by the user.
Of course, the work on the scheme does not stop here. There are several possible improvements. For instance, we imagine that it would be quite advantageous to allow for random polarization in the integration of the most computationally intensive n + 2 tree-level, and n + 1 one-loop amplitudes. This can definitely be achieved by polarized splitting functions. Another issue to consider is the introduction of cutoffs on the subtraction phase space. In any case, more experience has to be accumulated in order to decide, which modifications to include first. We leave this to future work.

Phase spaces
Tr

Appendix B. Spherical coordinates in d dimensions
Let d d r be the Euclidean integration measure in R d . We can decompose it into a radial and an angular part with the help of a δ-function insertion, if we rescale the r vector as r = rn We have thus defined a rotationally invariant measure, dΩ, on the unit (d − 1)-sphere, S d−1 1 . Notice that we will, from now on, include the dimensionality in the notation of the versorsn. Let us introduce a recursive parameterization in terms of anglesn The volume of the unit (d − 1)-sphere is (B.6) We will also need the following result which implies the correct reduction of the dimensionality of space We further introduce a representation of the angular versor parameterization through rotations of a basis vector. To this end, we definen and a d × d rotation matrix transforming the coordinates i and j where we have introduced the shorthand notation Due to the commutation properties of the rotation matrices, there is (B.14)

Appendix C. Infrared divergences of virtual amplitudes
We consider renormalized on-shell virtual amplitudes including wave function renormalization factors, which are non-trivial for external massive quarks. The strong coupling is assumed to be renormalized in the MS scheme with decoupling of massive quarks. Infrared divergences can be factorized from virtual amplitudes as follows where the infrared (IR) renormalization constant Z is an operator in color space, and depends on the momenta {p i } = where the anomalous dimension operator Γ is given by [24,25,26,27,28,29] The triple color correlations given in the third and fourth lines of Eq. (C.6) cannot contribute to the divergences of spin and color summed amplitudes at next-to-next-to-leading order, as long as the Born amplitudes do not contain complex couplings or masses [30]. The explicit solution of the RGE Eq. (C.5) can be found in Ref. [31], and reads up to order α 2 where the leading beta-function coefficient is with n l the number of massless quark flavors. The expression contains the anomalous dimension Γ and its derivative Γ is given in terms of the anomalous dimensions γ cusp , γ q , γ Q , γ g , and two functions F 1 and f 2 . We give explicit formulae for the coefficients of the expansion in α s which we have taken literally from Refs. [25,31]. The massless cusp anomalous dimension is T F n l . (C. 12) In the massive case the cusp anomalous dimension can be written as For massless quarks (anti-quarks) we also have whereas the massive quark (anti quark) anomalous dimension is The anomalous dimension for gluons reads Finally, we take the functions F 1 and f 2 from Ref. [28] The matrix element factorizes as  where the averaged splitting functions are We are also interested in the case without summation over the polarization of the final state gluon. Let us assume that the latter has momentum p 1 . The splitting functions depend on the polarization vector ε µ 1 , which, for our purposes, may be assumed to be real. The polarized splitting functions read We recover the unpolarized splitting functions, if we sum over the gluon polarizations spin ε µ 1 ε ν 1 = −g µν + p µ n ν + p ν n µ p · n , The initial state collinear limit can be recovered from the given formulae with minor replacements. Both collinear configurations are depicted schematically in Fig. (D.5). The essential difference is in the direction of the momenta. All splitting functions require the following replacement P a 1 a 2 −→ − 2s a +2s a 1P a 1 a 2 , (D. 16) where s a and s a 1 are the spins of partons a and a 1 respectively. The splitting variable z can be obtained in the collinear limit from the energies of the involved partons. The crossing amounts to the replacement Let us now turn to the triple-collinear limit. Consider the set of three vectors where as before p 2 = n 2 = p · k ⊥i = n · k ⊥i = 0. This configuration fulfills no other constraints, but rather the limits are expressed through derived variables The factorization formula is obtained in the limit k ⊥i → 0 and reads with s 123 = (p 1 + p 2 + p 3 ) 2 and x = x 1 + x 2 + x 3 . In the following we will drop the superscript, (0), inP (0) a 1 a 2 a 3 in order not to clutter the notation beyond the necessary. The flavor of the parton a is obtained by flavor conservation. The complete set of splitting functions is taken from Ref. [32] (see also [33,34]). In the case of spin conservation, we only give the averaged splitting functions P a 1 a 2 a 3 P ss a 1 a 2 a 3 = δ ss P a 1 a 2 a 3 .

Appendix F. Splitting functions
For the collinear factorization contribution, we need the splitting functions up to O(α s ) [43] P q i q j (x, α s ) = δ i j P (0) qq (x) + The leading order contributions are Beyond leading order one writes the splitting function P q i q j in terms of a flavor singlet (S) and non-singlet (V) contribution P q i q j (x, α s ) = δ i j P V qq (x, α s ) + P S qq (x, α s ) , (F.9) P q iq j (x, α s ) = δ i j P V qq (x, α s ) + P S qq (x, α s ) . (F.10) The next-to-leading order contribution to the splitting functions are P V(1) qq (x) = C 2 F − 2 ln x ln(1 − x) + 3 2 ln x p qq (x) − 3 2 + 7 2 x ln x − 1 2 (1 + x) ln 2 x − 5(1 − x)