On the factorization of overlapping singularities at NNLO

Real and virtual corrections in NNLO QCD require multi-dimensional integrals with overlapping singularities. We first review ideas and methods which have been proposed for performing such computations. We then present a new method for the factorization of overlapping singularities based on non-linear integral transformations. We apply this method for the evaluation of all integral topologies which appear in double real radiation corrections in cross-section calculations for the production of a heavy system at hadron colliders. Finally, we demonstrate with typical examples that two-loop virtual corrections are amenable to the same method.


Introduction
Modern accelerator experiments require precise perturbative calculations for the event rate of a variety of physical processes. Jets, electroweak gauge bosons and heavy quarks are being produced copiously at the Tevatron and the LHC. The precision of the measurements of physical masses, coupling parameters and the structure of colliding hadrons depends significantly on theoretical uncertainties which are better controlled at higher orders in perturbation theory. The exclusion of hypotheses for novel particles and interactions is more significant when candidate signal processes are predicted accurately. With the arrival of new discoveries, the nature of physics laws will be deciphered more confidently with the aid of solid quantitative theory predictions.
Our abilities to simulate complicated physical processes beyond the leading order (LO) have been improved dramatically in the last few years. At next-to-leading-order (NLO), previously inaccessible calculations with up to five particles in the final state are now possible [1]. Basic collider processes with fewer particles have also been computed at nextto-next-to-leading (NNLO) order in QCD [2][3][4][5][6][7][8][9][10][11][12][13]. For hadron colliders, the experimental frontier in particle physics, only cross-sections for 2 → 1 processes have been computed at NNLO. Such computations must be extended to 2 → 2 processes which are relevant to the experimental program. These include top and bottom quark pair production, inclusive jet production, electroweak diboson production, electroweak gauge boson and Higgs production in association with jets, single top production and beyond the Standard Model signals. It is unclear whether existing methods are suited to this task, and refinements of traditional methods or the development of new ones are required in order to face the increased complexity of such calculations.
A fundamental technical difficulty in NNLO calculations is the appearance of multidimensional integrals over the momenta of up to two additional real or virtual particles with respect to the Born process. These integrals are separately infrared divergent and only their sum is finite. Higher order computations are performed almost exclusively within dimensional regularization, where real and virtual corrections are expanded in a dimensional regulator = 2 − D 2 , where D is the number of dimensions. Laurent expansions in are intricate in the presence of overlapping singularities. In this paper we present a new method for the calculation of the Laurent series in of multidimensional integrals which typically appear in NNLO computations and generic higher order computations.
Existing methods which tackle or bypass the problem of overlapping singularities are based on differential equations [14][15][16], Mellin-Barnes representations [17][18][19][20] and sector decomposition [21][22][23]. The first two approaches can be applied to the calculation of virtual or inclusive real radiation corrections. A subtraction method can reduce the problem of fully differential real radiation calculations at NNLO to inclusive phase-space calculations [10][11][12] enabling the method of differential equations and Mellin-Barnes for fully differential calculations. Sector decomposition can be applied universally, for virtual inclusive phase-space integrations and fully differential integrations of real radiation matrix-elements.
Sector decomposition has been the first successful method for performing fully differential NNLO calculations for hadron collider processes [7,8]. This is largely attributed to the conceptual simplicity of the method and its algorithmic nature which permits a full automatization. The algorithm eliminates overlapping singularities by slicing the integration volume such that variables which contribute to an overlapping singular limit are ordered. In this way, the singularity is always factorized and it appears as a single singular limit of only the largest variable. While this algorithm leads to numerically stable evaluations, it requires a large number of integrals due to the slicing of the integration volume. This hinders the application of the method to processes with more complicated matrix-elements.
We present here an alternative method for the factorization of overlapping singularities. We have observed that a factorization is possible by means of simple rescaling of the integration variables and non-linear transformations which preserve the geometry of the integration boundaries. Our method leads to a rather small number of numerically stable integrations.
We apply our technique to all singular integral topologies which appear in the evaluation of NNLO double real radiation corrections to production processes of a massive system at hadron colliders. We present suitable phase-space parameterizations, analyze the singularity structure of the matrix-elements, and demonstrate how to obtain their expansion in with simple changes of integration variables. We then demonstrate how our technique can be applied to very complicated and maximally singular two-loop master integrals. In massless two-loop boxes overlapping singularities are very hard to treat with non-linear transformations, and we have not been able to find suitable ones which factorize them completely. On the other hand, a hybrid approach of non-linear transformations combined with sector decomposition is straightforward and more efficient than applying only sector decomposition.
In Section 2, we review existing methods for the Laurent expansion in the dimensional regulator of integrals in higher perturbative orders. In Section 3 we present our method and we demonstrate it on typical examples of integrals with overlapping singularities in Section 4. In Section 5 we discuss phase-space parameterizations and the singularity structure of double real radiation at NNLO. In Section 6 we present the numerical evaluation of integrals from all topologies which appear in double real radiation corrections at NNLO. In Section 7 we apply our method to maximally singular two-loop integrals, the crossedtriangle and the crossed-box. Finally, we present our conclusions in Section 8.

Laurent expansion of Feynman integrals in the dimensional regulator
Loop integrals and phase-space integrals for the calculation of production rates of physical processes in quantum field theory are divergent in four space-time dimensions. In dimensional regularization, d = 4 − 2 , all divergent integrals are computed as a Laurent expansion in the dimensional regulator . This task is tedious due to physical singularities, corresponding to infrared and collinear configurations of real and virtual particles. Singular manifolds in the integration volume are of increasing complexity at higher orders in perturbation theory. We shall consider examples of physical loop and phase-space integrals in later sections of this paper. Here, we shall present illustrative mathematical examples with similar singular behavior as in realistic cases.
The easiest category of singular integrals is when divergences in the integrand occur as poles of a single variable. Consider an integral with f ( x) being a finite function in the integration volume. This integral is divergent for = 0 due to the pole in x 1 = 0. To expand in we use a subtraction technique, isolating the pole contribution, In the first term, we are allowed to perform a Taylor expansion in , given that the integrand is finite in the limit x 1 → 0. In the second term, we can perform the integration in easily. We then have All integral coefficients in the Laurent series of the last expression can be evaluated numerically. In case of many factorized singularities, we can apply readily the same procedure, and obtain a Laurent expansion in with the substitution We note that one may also encounter singularities due to poles of second or higher order, as for example in The subtraction method can be also applied here, writing The integrals in x of the above expression can be computed numerically (as an expansion in ). The extraction of divergences is more complicated for integrals with overlapping singularities. Consider as an example the integral (2.9) For = 0 the integrand becomes divergent when both x 1 , x 2 → 0. Here, it is easy to perform successively both integrations, finding the explicit result (2.10)

The differential equation method
Analogous problems in realistic NNLO calculations are very hard to tackle with direct analytic integrations. A powerful method which has found numerous applications is the method of differential equations. In this approach we find a physical parameter for the integral and formulate a differential equation using integration by parts. In our example, we can write a differential equation with respect to the parameter c, by integrating the total derivative This yields the differential equation, The inhomogeneous term on the right side of the above equation is simpler than I. Specifically, (2.13) and we find 14) The general solution of Eq. 2.12 involves integrals over one variable only, thus bypassing the problem of overlapping singularities. The constant of integration can be determined from knowing the integral at a special value of c or by exploiting a known limiting behavior or scaling. For example, in our case, we could be using that which we can easily derive with a change of variables x 1 ↔ x 2 in Eq. 2.9.

The Mellin-Barnes representation method
Mellin-Barnes representations allow a straightforward Laurent expansion of Feynman integrals by using Cauchy's theorem. Such representations are obtained by using the identity, where the contour of integration is a straight line parallel to the imaginary axis, crossing the real axis at a point w 0 such that the real part of the arguments of the gamma functions are positive. Using Eq. 2.17 and integrating x 1 and x 2 for the toy example of Eq. 2.9, we obtain the Mellin-Barnes representation where the representation is valid (all Gamma functions have arguments with positive real part) if we choose, for example, = −0.9 and Re w = −0.2. Notice that we cannot find any value of Rew which renders the integral well defined if we choose = 0. This means that the integral develops a pole in . A Laurent expansion can be achieved with an analytic continuation method, moving the value of from a value that the integral is well defined, = −0.9 in our example, to = 0 and isolating with Cauchy's theorem the poles which arise when the arguments of dependent Gamma functions become zero or negative integers. In our example, we find that Γ(−1 − − w) = Γ(0) develops a pole as = −0.8 (and w = w 0 ). No other pole is encountered by continuing the value of further to = 0. We can then write, I = TaylorExpand(I) =0 + Res w=−1− (I). (2.19)

The subtraction method
The differential equation method and the Mellin-Barnes method bypass the problem of overlapping singularities by integrating out Feynman parameters and phase-space variables and generating equivalent representations where overlapping singularities cannot occur. Both methods rely on the integration volume being well known and free of parameters, other than the space-time dimension. This is the case for loop integrals and inclusive phase-space integrations. An important class of phase-space integrals requires parametric boundaries which are determined according to varied selection criteria for the experimentally measured observables. For such integrals the differential equation and Mellin-Barnes methods are not generally suitable. One approach is to use a subtraction method in order to map the problem of fully differential phase-space integrations onto a problem of fully inclusive phase-space. Consider the toy example, where the function J(x 1 , x 2 ) plays the role of selecting an arbitrary subregion of the integration volume according to, for example, the wishes of the experimentalists. Using subtraction, we can re-write (2.21) The first integral contains only an integrable singularity as → 0 and can be computed numerically. The second integral is a "fully inclusive" integral and may be computable with the Mellin-Barnes or differential equation method.

The sector decomposition method
A different approach is to use sector decomposition. We divide the integration region according to the relative magnitude of the integration variables which are required for the singular limit (in our example x 1 = x 2 = 0), by multiplying the integrand with an appropriate unity, This gives rise to two integration domains (sectors). In the sector with x 1 > x 2 we rescale x 2 = x 2 x 1 , and in the sector with x 2 > x 1 we transform x 1 = x 1 x 2 . We then obtain The singularities in both integrands are now factorized and a Laurent expansion can be easily achieved with a simple subtraction. The method of sector decomposition is suited for all types of loop and phase-space integrals. It is instructive to see how the method is used on a physical example. Let us consider the one-loop box scalar integral, (2.24) The corresponding Feynman parameterization reads, (2.26) To avoid creating poles at the upper limit of the x i integrations we apply first the method of primary-sectors [21]. We write We now rescale and perform the x integration. This yields (2.29) All terms in the sum can be computed in exactly the same fashion. For convenience, although not necessary, we use the special symmetry of this problem, y 1 ↔ y 3 and y 2 ↔ y 4 , and cast the integral as We observe that the integral becomes singular in the following instances We now apply sector decomposition to factorize the entangled singularity structure. We multiply the integrand with In each of the three sectors of the above equation we rescale the smallest variables with respect to the large ones, mapping the boundaries of the sectors to the unit cube. Specifically, Θ(y 2 y 3 < y 1 < y 2 ) : y 1 → y 1 y 2 and y 3 → y 3 y 1 (2.34) Θ(y 2 < y 1 ) : We then obtain a representation of the one-loop box, with a simple, factorized, singularity structure:

38)
(2.39) The resulting integrals I 1 , I 2 , I 3 of sector decomposition can all be expanded in with the subtraction method.

Factorization of singularities with non-linear transformations
In this section, we propose a new method for the factorization of overlapping singularities. We consider again the very simple case of an overlapping singularity in the two-dimensional toy example integral We shall perform a rescaling transformation over the entire integration region, This yields, We notice that there is a factorized singularity at x 1 = 0. In this singular point of the x 1 integration, the variable λ ranges up to +∞. However, the λ integration is convergent at = 0, since the integrand scales as 1 λ 2 for very large values of λ. Therefore, we could immediately treat the singularity at x 1 = 0 with the subtraction method We can then evaluate numerically the integrals which are produced after we perform an expansion in .
Alternatively, we could perform a transformation 1 on λ to bring the integration region back to [0, 1]. λ = g(u), Such a transformation maps the integration region to, It is very important to select carefully this transformation. A linear mapping is clearly ineffective, since it undoes the original rescaling of x 2 = λx 1 . However, non-linear mappings, such as are effective. For almost all practical applications in this paper we employ the mapping . (3.11) with δ often chosen equal to 1. Explicitly, the transformation with a Jacobian disentangles the overlapping singularity, transforming the integral of eq. 3.1 as 14) The singularity in the limit x 1 = 0 and = 0 can be subtracted away, and a Laurent series expansion around = 0 is achieved using the expansion In this approach, we have achieved to factorize the overlapping singularity with a simple transformation. In comparison, a factorization with sector decomposition doubles the number of integrals, as we have seen in the previous section. Economizing in the number of integrals is even more significant for physical applications where entanglement of singularities with more variables may take place.
Let us now revisit the one-loop box calculation using the new method instead of sector decomposition. Following the "analytical-transformation" approach, we perform the change of variables on the integral of Eq. 2.30 This yields the integral In this integral, the singularities have nicely factorized in the term (y 2 y 3 ) −1− . In comparison to sector decomposition, we now have to perform one integration rather than three.

Characteristic forms of entangled singularities
In this section, we present some typical examples of integrals with overlapping singularities and the mappings that we use to disentangle them. Our second example is the integral which is a trivial extension of eq. 3.1. We use the mapping where we have also designated the singularity structure of the integral. This mapping leads to where the singularities are factorized in terms of y and z.
Next, let's consider Here we use the simultaneous double mapping which leads to Next, let's consider the integral (4.7) Here we use successively The integral then becomes It is maybe instructive to see how this integral of eq. 4.7 behaves under simple rescaling. Consider z = λ z w and then x = λ x yw. We get (4.10) Only the integral over λ x extends to infinity and in that limit the behavior of the integrand is dλ x /λ 2 x which vanishes at infinity. Next let's consider (4.11) Here we use the successive mappings which brings the integral to the factorized form When overlapping singularities appear together with factorized singularities in the same variable, a slight complication appears. Consider the integral . (4.14) It has a factorized singularity at x = 0 and an overlapping singularity at x = 0 = y. Let us call the singularity at x = 0 active and the one at y = 0 passive. In order to disentangle the singularity we would like to use the same mapping as in the previous examples. It turns out that we can do this, but only as long as we choose to remap the active singularity, i.e.
x(x + y) : The integral then becomes which can be subtracted easily. Note that applying the wrong rescaling y = λx one gets We can see immediately that the dλ integral is logarithmically divergent at the active singularity x → 0 due to the upper limit of the integration region. On the contrary, applying the correct rescaling x = λy one gets which, as y → 0, behaves as dλ/λ 2 which is finite. We see that the simple λ-rescaling works as a guideline, showing when a mapping properly factorizes the singularities of an integral. It will be instrumental in more complicated cases presented below.
Let us now consider the integral 2 . (4.19) Identifying x = 0 = z and y = 0 = z as two independent overlapping singularities, where z is passive and both x and y are active, we know from the previous example that we should not map z from the previous example. So one is left to map x → z or y → z. By the symmetry of the integrand it does not matter which one of these one can choose. Let us choose x and a slightly modified mapping that keeps the expressions simpler: The integrand then becomes We have now managed to "activate" the singularity at z = 0. At the same time the singularities at y = 0 and at x = 0 have remained active as well. However, we now find an overlapping singularity at z = 0 = y and x = 1, where the singularity at x = 1 is passive. Notice that we started with two independent (partially interfering) overlapping singularities, have treated one of them and are now left with only one, which lies at a different point. We shall remap z and y as follows The integrand becomes . (4.23) Note that the remaining singularity of the integrand is integrable. Let's now explore the potential of a slightly different kind of mapping. We have the integral with a, b independent of x but potentially dependent on y i . In the latter case the integral might have overlapping singularities, as a, b → 0 or x, b → 0. We employ and get If N ≥ 2 this mapping factorizes the singularity at b → 0 and, at the same time, exposes the a + b structure of the overlapping a, b → 0 singularity, making it ready for further mappings.
Let's see, as an example, (4.27) We map: to get We can now use the mapping of eq. 4.25 with a = x 3 + x 4 and b = x 4 x 5 to get where F (x i ) is a finite function of x i . Noting that x 4 is an active singularity, we use the mapping of eq. 4.25 again with a = 1 + x 5 and b = x 3 to get (4.31) Let us now see some examples where we employ a hybrid method of one-step of sector decomposition and non-linear transformations to factorize overlapping singularities. A similar singularity structure appears in two-loop massless box integrals.
We consider (4.32) We split this integral in two sectors which has the singularity structure of eq. 4.4 and can be factorized by the mapping of eq. 4.5, and which is of the type of eq. 4.7 and can be factorized with the mapping of eq. 4.8.
We now move to the most complicated example of this section, the integral with A, E finite and B, C, D finite functions of τ 1,2 . We split this integral in two sectors, x 1 , x 4 , and we get which is of the type of eq. 4.1 and can be immediately factorized with the mapping eq. 4.2.
The other sector is This should be further split in x 3 , x 4 to get (4.38) which is of the type of eq. 4.4 and we can use the mapping eq. 4.5 to factorize it. The other sector is (4.39) and requires further splitting. We choose to split in the variables x 4 , τ 1 to get (4.40) which is of the type of eq. 4.11, and (4.41) which requires a final split in x 4 , τ 2 , (4.42) which is finite and (4.43) which is of the type of eq. 4.7. The original integral can be written in terms of its five sectors as I 10 = I 10A + I 10B1 + I 10B2A + I 10B2B1 + I 10B2B2 . Finally, let's consider with A, E, C, D finite and B 1,2 = 2 − τ 1,2 , also finite for all values of τ 1,2 . We begin by splitting the integral in τ 1 ,τ 2 . We get (4.46) and We notice that we can get I 11B from I 11A if we exchange C and D and rename the dummy integration variables τ 1 ↔ τ 2 , so that (4.48) We will use the same decompositions and mappings to factorize I 11B and I 11A (with τ 1 and τ 2 interchanged), so we only describe the latter below. We split I 11A in two sectors with respect to τ 1 and x 4 . We get and then which factorize all singularities. We also have (4.52) We will split with respect to τ 1 ,τ 2 to get (4.53) which can be factorized by and τ 2 > τ 1 : which can be factorized by (4.56) The original integral can, therefore, be factorized in six different integrals:

Double real radiation for final states with massive particles
One of the major challenges at NNLO in QCD has been the computation of the double real emission part of the cross-section. While the computation of the matrix elements with N +2 particles in the final state is not a problem per se, difficulties arise when one integrates over the phase space of the two unresolved particles. The corresponding integrals are infrared divergent in the soft and collinear limits and are dimensionally regulated. The divergences have to be subtracted before the integrals can be numerically evaluated. As long as the singularities are factorized, as they usually are at NLO, it is straightforward to use a Laurent expansion over the singular variables, and evaluate its coefficients numerically. At NNLO, the singularity structure of the integral is more intricate, as line and overlapping singularities appear, and the desired factorization is not straightforward. The method of sector decomposition has already been applied successfully to achieve this factorization for hadron collider [7,8] and decay processes [24][25][26]. A drawback of the method is that it leads to a large number of sectors. The goal of this paper is to replace sector decomposition for double-real radiation integrals with an economical factorization method based on non-linear transformations.

Infrared singularities in double real radiation
We consider double real emission to a generic NNLO 2 → n + 2 process (see Fig 1) with n massive particles and 2 massless partons in the final state. We denote the momenta of the incoming particles by q 1 , q 2 , those of the outgoing massive particles by p 1 ..p n and those of the two unresolved partons by q 3 and q 4 . Infrared singularities in this phase space will occur whenever q 3 and/or q 4 become soft or collinear to q 1 , q 2 or to each other. For the case of double real radiation to the production of a single massive particle (e.g. Higgs, W or Z production) potentially singular propagators can be summarized as Note that s 123 , s 124 are bounded from below. Further soft singularities can be found if there are colored massive particles in the final state, which can radiate off soft gluons. One can then get also the following possibly singular denominators: the soft singularity is factorized in E 3 . Whenever some heavy colored state radiates off two gluons we can also get the denominator t 34i , it can only become singular in the double soft limit when E 3 = 0 = E 4 . However the double soft limit will always be factorized as we will show in the next section. Let us now discuss the denominator structure of the most singular diagrams which one could expect in double real radiation: those where radiation is emitted by initial state particles. We have illustrated the propagator structure of these topologies using some diagrams containing gluons in Fig.2 (diagrams containing just massless quarks correspond to the same topologies).
Diagrams whose propagator structure can be related to the ones in Fig.2 by a simple interchange of q 3 with q 4 or of q 1 with q 2 will also fall into the same topology.
By considering square and interference terms of the topologies C 1 , C 2 and C 3 , we obtain the following list of integrals: 7. Topology C 4 ⊗ C 4 : 10. Topology C 4 ⊗ C 3 : Where dΦ 3 is the differential double emission phase space element for 2 + n final state particles, and N ({s ij }) is in general a finite function of the kinematical invariants. The topology C 4 contains only soft singularities similar to those in C 1 . The topologies C 4 ⊗ C 4 and C 4 ⊗ C 1 are, therefore, easier than C 1 ⊗ C 1 . They can be treated exactly like C 1 ⊗ C 1 and we will not discuss them in what follows.

Phase-space of double real parton radiation
We would like to point out that different parameterizations of the phase-space can factorize different sets of kinematic invariants. We will derive two such parameterizations which allow for a more convenient numerical evaluation of diverse diagrams, according to their topology.
The phase-space of n massive particles in four dimensions is: where s = (q 1 +q 2 ) 2 . We assume that a 2 → n process exists at leading order in perturbation theory, and a strictly four-dimensional evaluation is therefore sufficient. At NNLO, the double emission phase space is given by including two further massless particles (whose momenta we denote by q 3 and q 4 ) We factorize the double real phase space into a 3-particle phase space times an n-particle phase space as follows and parameterizing s 1..n linearly we get (5.20) The parameter x 5 ∈ [0, 1] then uniquely defines the double soft limit when x 5 = 1. In the following discussion we will use the variable z = s 1..n s (5.21) which in the special case of n = 1 reduces to z = m 2 1 s . Then the variable x 5 is trivially removed and the double soft singularity occurs whenever s = m 2 1 . In the following we will assume that one can parametrize the n-particle phase space dΦ n , and we will focus on the phase-space of the potentially unresolved massless partons dΦ 3 .

Energies and angles parameterization
The three particle phase space element dΦ 3 is Integrating out Q and using that d d qδ (+) (q 2 ) = dEE d−3 dΩ (d−1) /2 we get cos θ 34 )). (5.23) We can solve the delta constraint for the energies in a symmetric way using the following ansatz: We find The double soft limit now appears when z → 1, while the single soft singularities occur as x 1 → 0, 1. After this transformation the phase space volume becomes .
(5.26) Having solved the energy constraint we move on to parametrize the angles. Choosing the z-axis as the direction of q 1 , we directly parameterize the angles which q 3 and q 4 make with the z-axis. Finally we parametrize the angle φ between q 3 and q 4 in the x-y plane leading to the following expressions of the solid angles Suppressing any extra dimensional components our 4-vectors are then fully parametrized as q 3 = E 3 (1, sin θ 3 , 0, cos θ 3 ) and q 4 = E 4 (1, sin θ 4 sin φ, sin θ 4 cos φ, cos θ 4 ). Mapping the remaining angles linearly, i.e. cos θ 3 = 2x 3 − 1, cos θ 4 = 2x 4 − 1 and φ = x 2 π, one obtains The following lists the propagators of massless partons in this parameterization: and The angle between q 3 and q 4 is related tõ We can now apply this parameterization to all integrals of type C 1 ⊗ C 1 ,C 3 ⊗ C 3 and C 1 ⊗ C 3 .

Line singularities in the energy and angles parameterization
One can use a non-linear transformation to get rid of the overlapping structure inx 34 [27]. A convenient way to derive such a mapping is remappingx 34 fromx − 34 =x 34 (φ = 0) tõ x + 34 =x 34 (φ = 1) usingx It is then apparent thatx 34 will vanish wheneverx − 34 orx + 34 will, for any value of x 2 . And that the overlapping line singularity is then re-casted into just a line singularity. To aid numerical stability we perform the mapping x 2 → (1 − cos(x 2 π))/2, such thatx 34 becomesx .
This is in fact identical to the mapping in [28]. The phase space volume then becomes To factorize the line singularity in s 34 (at x 3 = x 4 ) we are forced to split the integration region in two, separating x 3 < x 4 from x 4 < x 3 .

Hierarchical parameterization
Since in the energy and angles parameterization the invariants s 34 , s 134 , s 234 had line and overlapping singularities, it is worth having a second parameterization which factorizes these, but may not factorize the others. Our second parameterization closely resembles the features of the rapidity parameterization published in [7], however it is somewhat simpler. In this parameterization the three particle phase space element dΦ 3 is 39) is first factorized into a product of two 2-particle phase spaces We can parameterize dΦ 2 ( √ s, √ s 34 , √ s 1..n ) in terms of s 134 , yielding where φ is the angle between (p 3 ) ⊥ andQ ⊥ . We fulfil the constraint (p 3 ) ⊥ sin φ ≥ 0 to find the limits of integration for s 13 and then for s 23 . Parameterizing s 134 , s 34 , s 13 and s 23 linearly we arrive at The invariants in this parameterization are and We see that the only invariants which are not factorized are s 13 and s 14 . The variable s 13 contains overlapping singularities at x 3 = 0 = x 2 and x 3 = 1 = x 2 as well as an overlapping line singularity at x 4 = 0, x 1 = 1, x 3 = x 2 , while s 14 contains overlapping singularities at x 3 = 0, x 2 = 1 and x 3 = 1, x 2 = 1 as well as an overlapping line singularity at x 4 = 1, and we are left with just overlapping singularities, which can be treated as explained in the following section. This trick was first discovered by Frank Petriello [29] and it has been used in the implementation of the program FEHiP described in [7], it has been also been used in the evaluation of doublereal counterterms in [30].

Numerical evaluation of double-real radiation phase-space integrals
In this section, we present a numerical evaluation of all types of scalar phase-space integrals which appear in NNLO calculations. To evaluate our integrals numerically we choose the point (s = 1, z = 0.1). We will use the notation where the x i ∈ [0, 1] are parameters of integration.
1. Topology C 1 ⊗ C 1 : . (6.13) This contains the following substructure with A, B, C finite. This becomes singular when x 3 = 0 = x 2 where x 3 is active. We factorize this singularity by applying with A, B, C finite. This becomes singular whenx 3 = 0 =x 2 withx 3 being active. We disentangle this singularity by applyinḡ We then obtain with the numerator structure as in [7]. The singularities factorize in energies and angles. We immediately obtain 7. We will now consider interferences of C 4 with C 2 and C 3 . One can evaluate these interferences in the energies and angles parameterization. In the following we will use t 13 ∼ E 3 ∼ (s 13 + s 23 ) and t 24 ∼ E 4 ∼ (s 14 + s 24 ).
(a) Topology C 2 ⊗ C 4 : The integral I 24 = dΦ 3 s 34 s 134 (s 13 + s 23 )(s 14 + s 24 ) (6.35) has the following singularity structure in the energy and angle parameterization after the mapping (see 5.3.1) is applied. We first split the integration region into two sectors which we define as x 3 < x 4 (sector 1) and x 4 < x 3 (sector 2). After this sector decomposition we are still left with overlapping singularities at x 3 = 0 =x 1 in sector 1 and at x 4 = 0 = x 1 in sector 2. These can be disentangled usinḡ in sector 1 and in sector 2. We then obtain has the following singularity structure 1 It contains no line singularity but several overlapping ones located at x 3 = 0 = x 4 , x 4 = 0 = x 1 and at x 3 = 0,x 1 = 0. To separate the two singularities we first partial fraction the soft singularities by multiplying by 1 = x 1 +x 1 . We then treat the two terms with different nonlinear transformations. For the first term we apply the mapping since x 3 is the only active singularity, it is clear that we had to remap it. The second term is more difficult, since both x 3 and x 1 are now active. We apply the following sequence of mappings: First let and then .

Two loop examples
In what follows we show how one can use non-linear mappings to disentangle singularities in two-loop integrals appearing in NNLO virtual amplitudes. We treat only two indicative cases, the massless non-planar triangle with one leg off-shell and the massless non-planar box with all legs on-shell, due to their particularly intricate singularity structure. Integrals involving masses are in general simpler as far as factorization of singularities is concerned.

7.1
The massless non-planar triangle with one leg off-shell.
The two-loop, non-planar triangle with one off-shell leg (see fig. 4) and momenta p 1 , p 2 A Feynman parameterization reads: The first overlapping singularity is at x = 0 or x = 1 and y = 0. We also notice that there is a singularity at y = 1. To avoid infinite looping we must first guarantee, as in sector decomposition, that no singularities occur at the upper limit of integration.
We split x in the two intervals R a = [0, 1/2] and R b = [1/2, 1] and map the integration region back to the unit hypercube. In R a , x → x/2, This gives two identical integrals (the integral in eq.7.2 is invariant under the combined , and we can write Note here that the denominator of eq. 7.3 has the same singularity structure as x + yzx 1 x 2 − yzx(x 1 + x 2 ) (7.4) In particular, there are still singularities at upper corners of the hypercube, when y → 1 and z → 1. This leads us to split the y integration region, y → y/2 and y → 1 − y/2 and the z → z/2 and z → 1 − z/2.
Xtri 11 has a singularity structure equivalent to that of x + yzx 1 x 2 , similar to eq. 4.1 and we use the mapping (directly analogous to eq. 4.2) x → xyzx 1 x 2 1 − x + yzx 1 x 2 . (7.10) Xtri 12 has a singularity structure equivalent to that of x+yx 1 x 2 , also similar to eq. 4.1 and we use the mapping x → xyx 1 x 2 1 − x + yx 1 x 2 .
(7.11) Xtri 21 has a singularity structure equivalent to that of x+zx 1 x 2 , also similar to eq. 4.1 and we use the mapping x → xzx 1 x 2 1 − x + zx 1 x 2 . (7.12) Xtri 22 is a bit more complicated. Its singularity structure is the one of x + x 1 x 2 − x(x 1 + x 2 ) (7.13) It retains singularities at x 1,2 → 1. We therefore split this integral further in x 1 → x 1 /2 and x 1 → 1 − x 1 /2 as well as x2 → x2/2 and x 2 → 1 − x 2 /2. We obtain (7.14) with a singularity structure equivalent to x + x 1 x 2 (i.e. eq. 4.1) for which we will use the mapping (7.16) with a singularity structure equivalent to x 1 + x(x 1 + y + z + x 2 ), similar to eq. 4.7 for which we will use the following sequence of mappings (7.18) with a singularity structure equivalent to x 2 + x(x 2 + y + z + x 1 ),similar to eq. 4.7, for which we will use the following sequence of mappings (7.20) which is finite! We therefore end up with 7 different integrals to be numerically evaluated. This should be contrasted with the 64 number of sectors one arrives using sector decomposition. The numerical convergence of these integrals poses no additional problems and we have checked that the numerical result agrees with the analytic result known in the literature.
2. Xbox a2 : Xbox a where x 2 > x 1,3,4 has the rather intricate singularity structure of the example eq. 4.45. We follow the discussion given there and decompose it to six integrals.